Index: matrix_svd.c =================================================================== --- matrix_svd.c (revision 853) +++ matrix_svd.c (working copy) @@ -1,7 +1,7 @@ /* -Authors: Sarah Al-Assam, Stephen Clark, Dieter Jaksch (Oxford) and Chris Goodyer (NAG) -Date: September-December 2012 -(c) University of Oxford 2012 + Authors: Sarah Al-Assam, Stephen Clark and Dieter Jaksch (Oxford) and Chris Goodyer (NAG) + Date: $LastChangedDate$ + (c) University of Oxford 2014 */ /*! \file matrix_svd.c \brief This contains functions which perform an SVD on a matrix @@ -12,49 +12,40 @@ The first set of these threshold entries based on a tolerance to avoid the growth of numerical artefacts and reduce the matrix to only non-zero rows and columns. The second set split the matrix into a block structure and slve these blocks independently if that is appropriate. */ -#include "tnt_int.h" +#include "../headers/tnt_int.h" +#include "../headers/dec_matrix.h" +#include "../headers/dec_parallel.h" -#include "dec_matrix.h" -#include "dec_parallel.h" +#include "../headers/dec_matrix_blocks.h" -#include "dec_matrix_blocks.h" +#include "../headers/dec_la.h" -/*! \def MEMALLOCCHK - Macro for testing for a failed malloc when returning an error code */ -#define MEMALLOCCHK if (NULL == chk) return TNT_ERR_MEMALLOC; #ifdef TNT_OMP /*! \def MEMALLOCCHKP Macro for testing for a failed malloc when in an OpenMP parallel section */ -#define MEMALLOCCHKP if (NULL == chk) {omp_exit=1; continue;} -#include +#define MEMALLOCCHKP if (NULL == chk) {omp_exit=TNT_ERR_GEN; strcpy(tnt_err.errinfostr,"Out of memory"); continue;} #else /*! \def MEMALLOCCHKP Macro for testing for a failed malloc when in an OpenMP parallel section, but not compiled in */ -#define MEMALLOCCHKP MEMALLOCCHK -#endif - -#ifdef TIMINGS -static double timer=0.0; -static double timerS=0.0; -static double timerE=0.0; +#define MEMALLOCCHKP TNT_MEM_CHK #endif /*! Definition of external variable used for setting the function used for calculating the truncated SVD. It is initialised to use the 2 norm as a default */ -double (*tnt_matrix_trunc_err)(struct tnsr_values *, int, int) = tnt_matrix_trunc_2norm; +double (*_tnt_matrix_trunc_err)(struct tnsr_values *, int, int) = _tnt_matrix_trunc_2norm; /*! Definition of external variable used for setting tolerance for singular values of the truncated SVD. It is initialised to use -1.0 (i.e. all values kept up to chi) as default */ -double matrix_trunc_tol = -1.0; +double _tnt_matrix_trunc_tol = -1.0; /*! \ingroup matrix This function is a wrapper for the LAPACK SVD routine and also truncates the inner dimension for output matrices if necessary. The arguments M and N specify the dimension of the general - complex matrix A (so A is M x N). K specifies the inner dimensions Krearr or U->comarr (and similarly for S and VT) has not been allocated, so this is done as a first step - The error, a pointer to which is passed as an argument, is calculated using the function pointed to by tnt_matrix_trunc_error(). */ -TNT_int_err tnt_matrix_trunc_svd(int M, /*!< The number of rows in the matrix A */ - int N, /*!< The number of columns in the matrix A */ - int *Kp, /*!< The number of singular values with truncation < MIN(M,N) */ + The error, a pointer to which is passed as an argument, is calculated using the function pointed to by _tnt_matrix_trunc_error(). */ +TNT_int_err _tnt_matrix_svd(int *Kp, /*!< The number of singular values with truncation < MIN(M,N) */ struct tnsr_values *A, /*!< A pointer to the tnsr_values structure which represents A */ - struct tnsr_values *U, /*!< A pointer to the tnsr_values structure which represents U */ - struct tnsr_values *S, /*!< A pointer to the tnsr_values structure which represents S */ - struct tnsr_values *VT, /*!< A pointer to the tnsr_values structure which represents VT */ - double *err) /*!< The truncation error, which is calculated using the function pointed to by tnt_matrix_trunc_error() */ + struct tnsr_values *U, /*!< A pointer to the tnsr_values structure which represents U. No properties assumed to be set. */ + struct tnsr_values *S, /*!< A pointer to the tnsr_values structure which represents S. No properties assumed to be set. */ + struct tnsr_values *VT, /*!< A pointer to the tnsr_values structure which represents VT. No properties assumed to be set. */ + double *err) /*!< The truncation error, which is calculated using the function pointed to by _tnt_matrix_trunc_error() */ { - /* Define the workspace variables - see LAPACK notes below for explanation. */ - char jobu = 'S'; /* Set JOBU flag for the left singular vectors. */ - char jobvt = 'S'; /* Set JOBVT flag for the right singular vectors. */ - Complex *cwork; /* = WORK. */ - double *work; /* = WORK. */ - double *rwork; /* = RWORK. */ - int lwork; /* = LWORK. */ - int info = 0; /* = INFO. */ - int Korig; /* Original inner dimension of SVD */ - int K = *Kp; - int i, j; /* Row and column counter respectively */ - void *chk; /* Pointer for checking result of memory allocations */ - TNT_int_err ifail; /* Error code */ + int M = A->rowsize, N = A->colsize; /* Shorthand for the dimensions of A */ + int Korig = MIN(M,N); /* Original inner dimension of SVD */ + int K = *Kp;/* Shorthand for inner dimension */ + int j; /* Loop counter */ + TNT_ERR_RET_DEFS /* Error code */ int tmpI, tmpJ, tmpIJ; /* Temporary values of dimensions */ int nBlocks, blk, *blkDims, *outMatPos; - #ifdef TNT_OMP - int omp_exit = 0; - #endif + int omp_exit = TNT_SUCCESS; double scaled; /* The factor by which the A matrix is scaled to normalise the maximum entry */ double *Sloc, *Aloc, *Uloc, *VTloc; /* Pointers to temporarily created real blocks */ - Complex *Aloc_C, *Uloc_C, *VTloc_C; /* Pointers to temporarily created complex blocks */ - -#ifdef TIMINGS - struct timespec tp_start, tp_end; - clock_gettime(CLOCK_REALTIME, &tp_start); -#endif + tntComplex *Aloc_C, *Uloc_C, *VTloc_C; /* Pointers to temporarily created complex blocks */ - Korig = MIN(M,N); + /* Assign all the tensor values to have the same complex flag as the original matrix */ + U->comp = S->comp = VT->comp = A->comp; - lwork = -1; /* Queries LAPACK function for the optimal size for WORK. */ - - /* Call LAPACK SVD function with WORK query. */ + /* Set the outer dimensions of U and V */ + U->rowsize = M; + VT->colsize = N; + + /* Set the original inner dimension from untruncated SVD */ + U->colsize = VT->rowsize = Korig; + + /* Set rowsize and colsize of S to indicate that it is not yet a matrix (result passed back from LAPACK function is vector of diagonal values) */ + S->rowsize = S->colsize = -1; + + /* Calculate the total size of the arrays */ + U->sz = U->rowsize * U->colsize; + S->sz = Korig; + VT->sz = VT->rowsize * VT->colsize; if (A->comp) { /* If matrix A is complex */ /* CEG Matrix investigation */ - ifail = tnt_matrix_blockDiagonalisationComplex(M, N, A); - if (TNT_SUCCESS != ifail) - return ifail; - free(A->com_spare); /* Free here as we don't need non-blocked version again */ + ret = _tnt_matrix_blockDiagonalisationComplex(M, N, A); + TNT_RET_CHK /* NO_COVERAGE */ - tmpIJ = MIN(A->minDimI,A->minDimJ); U->minDimI = A->minDimI; U->minDimJ = A->minDimJ; - blkDims = tnt_matrix_FindBlockDimensions(A->minDimI, A->minDimJ, A, &nBlocks); + blkDims = _tnt_matrix_FindBlockDimensionsComplex(A->minDimI, A->minDimJ, A, &nBlocks); - outMatPos = tnt_matrix_FindOutputBlockPositions(blkDims, nBlocks); + outMatPos = _tnt_matrix_FindOutputBlockPositions(blkDims, nBlocks); - if (nBlocks>1) - { + if (nBlocks>1) { /* Allocate memory for new arrays */ /* Note these are enlarged to be square for now. Reduce appropriately after blk loop */ - chk = U->comarr = calloc((A->minDimI)*(A->minDimI),sizeof(Complex)); MEMALLOCCHK - chk = VT->comarr = calloc((A->minDimJ)*(A->minDimJ),sizeof(Complex)); MEMALLOCCHK + chk = U->comarr = calloc((A->minDimI)*(A->minDimI),sizeof(tntComplex)); + TNT_MEM_CHK /* NO_COVERAGE */ + + chk = VT->comarr = calloc((A->minDimJ)*(A->minDimJ),sizeof(tntComplex)); + TNT_MEM_CHK /* NO_COVERAGE */ /* S is always real */ - chk = S->rearr = malloc(MIN(A->minDimI,A->minDimJ)*sizeof(double)); MEMALLOCCHK + chk = S->rearr = malloc(MIN(A->minDimI,A->minDimJ)*sizeof(double)); + TNT_MEM_CHK /* NO_COVERAGE */ + -#ifdef TIMINGS - clock_gettime(CLOCK_REALTIME, &tp_end); - timerS += (double)tp_end.tv_sec-tp_start.tv_sec + 1.0e-9*(tp_end.tv_nsec-tp_start.tv_nsec); - clock_gettime(CLOCK_REALTIME, &tp_start); -#endif #ifdef TNT_OMP #pragma omp parallel default(none) \ shared (A, U, VT, S, blkDims, outMatPos, nBlocks, omp_exit) \ @@ -144,43 +126,35 @@ #pragma omp for schedule(dynamic,1) nowait #endif - for (blk=0; blkminDimI, A->minDimJ, A->comarr); - -#ifdef TIMINGS - clock_gettime(CLOCK_REALTIME, &tp_end); - timerS += (double)tp_end.tv_sec-tp_start.tv_sec + 1.0e-9*(tp_end.tv_nsec-tp_start.tv_nsec); - clock_gettime(CLOCK_REALTIME, &tp_start); -#endif + tmpIJ = MIN(A->minDimI,A->minDimJ); /* Allocate memory for new arrays */ - /* Note these are enlarged to be square for now. Reduce appropriately after blk loop */ - chk = U->comarr = calloc((A->minDimI)*tmpIJ,sizeof(Complex)); MEMALLOCCHK - chk = VT->comarr = calloc(tmpIJ*(A->minDimJ),sizeof(Complex)); MEMALLOCCHK - /* S is always real */ - chk = S->rearr = calloc(tmpIJ,sizeof(double)); MEMALLOCCHK - chk = cwork = malloc(4*sizeof(Complex)); /* Give WORK one entry initially. */ MEMALLOCCHK - chk = rwork = malloc(10*tmpIJ*sizeof(double)); MEMALLOCCHK - - /* Note that for only one block zgesvd now being called to return U as MxMIN(M,N) and VT as MIN(M,N)xN */ - jobu = 'S'; - jobvt = 'S'; - lwork = -1; /* Queries LAPACK function for the optimal size for WORK. */ - - zgesvd_(&jobu, &jobvt, &A->minDimI, &A->minDimJ, A->comarr, &A->minDimI, S->rearr, - U->comarr, &A->minDimI, VT->comarr, &tmpIJ, cwork, &lwork, rwork, &info, 1, 1); - - lwork = (int) cwork[0].re; /* Extract optimal WORK size. */ + chk = U->comarr = calloc((A->minDimI)*tmpIJ,sizeof(tntComplex)); + TNT_MEM_CHK /* NO_COVERAGE */ - free(cwork); + chk = VT->comarr = calloc(tmpIJ*(A->minDimJ),sizeof(tntComplex)); + TNT_MEM_CHK /* NO_COVERAGE */ - chk = cwork = (Complex *) malloc(20*lwork*sizeof(Complex)); /* Assign WORK its optimal size. */ MEMALLOCCHK + /* S is always real */ + chk = S->rearr = calloc(tmpIJ,sizeof(double)); + TNT_MEM_CHK /* NO_COVERAGE */ - /* Call the LAPACK function properly with optimal settings. */ - zgesvd_(&jobu, &jobvt, &A->minDimI, &A->minDimJ, A->comarr, &A->minDimI, S->rearr, - U->comarr, &A->minDimI, VT->comarr, &tmpIJ, cwork, &lwork, rwork, &info, 1, 1); + if ((A->minDimI > 0) && (A->minDimJ > 0)) { - free(cwork); - free(rwork); - - /* Rescale the S values back from their unnormalised ones */ - for (i=0; irearr[i] *= scaled; - -#ifdef TIMINGS - clock_gettime(CLOCK_REALTIME, &tp_end); - timer += (double)tp_end.tv_sec-tp_start.tv_sec + 1.0e-9*(tp_end.tv_nsec-tp_start.tv_nsec); - clock_gettime(CLOCK_REALTIME, &tp_start); -#endif + /* Call the SVD function for reduced SVD (1 for first argument) */ + ret = _tnt_la_zsvd(1, A->minDimI, A->minDimJ, A->comarr, S->rearr, U->comarr, VT->comarr); + TNT_RET_CHK /* NO_COVERAGE */ + } } /* Put entries from blocks back into the correct place in the matrices - i.e. with the zero rows/cols back in place */ - ifail = tnt_matrix_blockwiseRestoreComplex(M, N, A, U, VT, S); - if (TNT_SUCCESS != ifail) - return ifail; + ret = _tnt_matrix_blockwiseRestoreComplex(M, N, A, U, VT, S); + TNT_RET_CHK /* NO_COVERAGE */ - tnt_zeroMatSUVTComplex(M, N, S->rearr, U->comarr, VT->comarr); + zeroMatSUVTComplex(M, N, S->rearr, U->comarr, VT->comarr); /* Free rearrangement arrays */ - ifail = tnt_matrix_blockwiseFree(A); /* Generic for both real and complex */ - if (TNT_SUCCESS != ifail) - return ifail; + ret = _tnt_matrix_blockwiseFree(A); /* Generic for both real and complex */ + TNT_RET_CHK /* NO_COVERAGE */ free(blkDims); free(outMatPos); @@ -272,34 +216,29 @@ } else { /* If Matrix A is Real */ /* CEG Matrix investigation */ - ifail = tnt_matrix_blockDiagonalisation(M, N, A); - if (TNT_SUCCESS != ifail) - return ifail; - free(A->re_spare); /* Free here as we don't need non-blocked version again */ + ret = _tnt_matrix_blockDiagonalisation(M, N, A); + TNT_RET_CHK /* NO_COVERAGE */ - tmpIJ = MIN(A->minDimI,A->minDimJ); U->minDimI = A->minDimI; U->minDimJ = A->minDimJ; - blkDims = tnt_matrix_FindBlockDimensions(A->minDimI, A->minDimJ, A, &nBlocks); + blkDims = _tnt_matrix_FindBlockDimensions(A->minDimI, A->minDimJ, A, &nBlocks); - outMatPos = tnt_matrix_FindOutputBlockPositions(blkDims, nBlocks); + outMatPos = _tnt_matrix_FindOutputBlockPositions(blkDims, nBlocks); - if (nBlocks>1) - { + if (nBlocks>1) { /* Allocate memory for new arrays */ /* Note these are enlarged to be square for now. Reduce appropriately after blk loop */ - chk = U->rearr = calloc((A->minDimI)*(A->minDimI),sizeof(double)); MEMALLOCCHK - chk = VT->rearr = calloc((A->minDimJ)*(A->minDimJ),sizeof(double)); MEMALLOCCHK + chk = U->rearr = calloc((A->minDimI)*(A->minDimI),sizeof(double)); + TNT_MEM_CHK /* NO_COVERAGE */ + + chk = VT->rearr = calloc((A->minDimJ)*(A->minDimJ),sizeof(double)); + TNT_MEM_CHK /* NO_COVERAGE */ /* S is always real */ - chk = S->rearr = malloc(MIN(A->minDimI,A->minDimJ)*sizeof(double)); MEMALLOCCHK + chk = S->rearr = malloc(MIN(A->minDimI,A->minDimJ)*sizeof(double)); + TNT_MEM_CHK /* NO_COVERAGE */ -#ifdef TIMINGS - clock_gettime(CLOCK_REALTIME, &tp_end); - timerS += (double)tp_end.tv_sec-tp_start.tv_sec + 1.0e-9*(tp_end.tv_nsec-tp_start.tv_nsec); - clock_gettime(CLOCK_REALTIME, &tp_start); -#endif #ifdef TNT_OMP #pragma omp parallel default(none) \ shared (A, U, VT, S, blkDims, outMatPos, nBlocks, omp_exit) \ @@ -309,41 +248,36 @@ #pragma omp for schedule(dynamic,1) nowait #endif - for (blk=0; blkminDimI, A->minDimJ, A->rearr); + /* If only one block */ -#ifdef TIMINGS - clock_gettime(CLOCK_REALTIME, &tp_end); - timerS += (double)tp_end.tv_sec-tp_start.tv_sec + 1.0e-9*(tp_end.tv_nsec-tp_start.tv_nsec); - clock_gettime(CLOCK_REALTIME, &tp_start); -#endif + tmpIJ = MIN(A->minDimI,A->minDimJ); /* Allocate memory for new arrays */ - /* Note these are enlarged to be square for now. Reduce appropriately after blk loop */ - chk = U->rearr = calloc((A->minDimI)*tmpIJ,sizeof(double)); MEMALLOCCHK - chk = VT->rearr = calloc(tmpIJ*(A->minDimJ),sizeof(double)); MEMALLOCCHK - /* S is always real */ - chk = S->rearr = calloc(tmpIJ,sizeof(double)); MEMALLOCCHK - chk = work = malloc(4*sizeof(double)); /* Give WORK one entry initially. */ MEMALLOCCHK + chk = U->rearr = calloc((A->minDimI)*tmpIJ,sizeof(double)); + TNT_MEM_CHK /* NO_COVERAGE */ - /* Note that for only one block zgesvd is being called to return U as MxMIN(M,N) and VT as MIN(M,N)xN */ - jobu = 'S'; - jobvt = 'S'; - lwork = -1; /* Queries LAPACK function for the optimal size for WORK. */ - dgesvd_(&jobu, &jobvt, &A->minDimI, &A->minDimJ, A->rearr, &A->minDimI, S->rearr, - U->rearr, &A->minDimI, VT->rearr, &tmpIJ, work, &lwork, &info, 1, 1); + chk = VT->rearr = calloc(tmpIJ*(A->minDimJ),sizeof(double)); + TNT_MEM_CHK /* NO_COVERAGE */ - lwork = (int) work[0]; /* Extract optimal WORK size. */ - - free(work); - - chk = work = (double *) malloc(20*lwork*sizeof(double)); /* Assign WORK its optimal size. */ MEMALLOCCHK - - /* Call the LAPACK function properly with optimal settings. */ - dgesvd_(&jobu, &jobvt, &A->minDimI, &A->minDimJ, A->rearr, &A->minDimI, S->rearr, - U->rearr, &A->minDimI, VT->rearr, &tmpIJ, work, &lwork, &info, 1, 1); - - free(work); - - /* Rescale the S values back from their unnormalised ones */ - for (i=0; irearr[i] *= scaled; - -#ifdef TIMINGS - clock_gettime(CLOCK_REALTIME, &tp_end); - timer += (double)tp_end.tv_sec-tp_start.tv_sec + 1.0e-9*(tp_end.tv_nsec-tp_start.tv_nsec); - clock_gettime(CLOCK_REALTIME, &tp_start); -#endif + /* S is always real */ + chk = S->rearr = calloc(tmpIJ,sizeof(double)); + TNT_MEM_CHK /* NO_COVERAGE */ + + if ((A->minDimI > 0) && (A->minDimJ > 0)) { + /* Call the SVD function for reduced SVD (1 for first argument) */ + ret = _tnt_la_dsvd(1, A->minDimI, A->minDimJ, A->rearr, S->rearr, U->rearr, VT->rearr); + TNT_RET_CHK /* NO_COVERAGE */ + } } /* Put entries from blocks back into the correct place in the matrices */ - ifail = tnt_matrix_blockwiseRestore(M, N, A, U, VT, S); - if (TNT_SUCCESS != ifail) - return ifail; + ret = _tnt_matrix_blockwiseRestore(M, N, A, U, VT, S); + TNT_RET_CHK /* NO_COVERAGE */ - tnt_zeroMatSUVT(M, N, S->rearr, U->rearr, VT->rearr); + zeroMatSUVT(M, N, S->rearr, U->rearr, VT->rearr); /* Free rearrangement arrays */ - ifail = tnt_matrix_blockwiseFree(A); /* Generic for both real and complex */ - if (TNT_SUCCESS != ifail) - return ifail; + ret = _tnt_matrix_blockwiseFree(A); /* Generic for both real and complex */ + TNT_RET_CHK /* NO_COVERAGE */ free(blkDims); free(outMatPos); } - /* Discard the unwanted parts of the matrices, while adding the of their errors squared to the error */ + /* Discard the unwanted parts of the matrices, while calculating the truncation error */ /* Before discarding, check that the smallest value kept is larger than KTOL, and if it isn't keep reducing K until it is */ - for (j = K-1; j > 0; j--) { - if ((S->rearr[j]/S->rearr[0]) > matrix_trunc_tol) { + for (j = K-1; j >= 0; j--) { + if (S->rearr[j] > _tnt_matrix_trunc_tol) { break; } } @@ -441,26 +348,41 @@ if (K==Korig) { *err = 0.0; } else { + ret = _tnt_matrix_truncate(K,U,S,VT,err); + TNT_RET_CHK /* NO_COVERAGE */ + } /* NO_COVERAGE */ if (A->comp) { /* Discarding last columns by simply freeing last portion of memory */ - chk = U->comarr = realloc(U->comarr,M*K*sizeof(Complex)); - MEMALLOCCHK - /* Shuffle elemnts of VT then delete last entries to remove rows. */ + chk = U->comarr = realloc(U->comarr,M*K*sizeof(tntComplex)); + TNT_MEM_CHK /* NO_COVERAGE */ + /* Shuffle elements of VT then delete last entries to remove rows. */ for (j = 0; jcomarr[i + j*K].re = VT->comarr[i + j*Korig].re; VT->comarr[i + j*K].im = VT->comarr[i + j*Korig].im; } } - chk = VT->comarr = realloc(VT->comarr,N*K*sizeof(Complex)); - MEMALLOCCHK + chk = VT->comarr = realloc(VT->comarr,N*K*sizeof(tntComplex)); + TNT_MEM_CHK /* NO_COVERAGE */ + /* S is real vector so remove last singular double values */ + chk = S->rearr = realloc(S->rearr,K*sizeof(double)); + TNT_MEM_CHK /* NO_COVERAGE */ + } + + } else { + + /* If matrix shrunk to zero, simply free arrays */ + if (0==K){ + free(U->rearr); + free(VT->rearr); + free(S->rearr); } else { /* Discarding last columns by simply freeing last portion of memory */ chk = U->rearr = realloc(U->rearr,M*K*sizeof(double)); - MEMALLOCCHK + TNT_MEM_CHK /* NO_COVERAGE */ /* Shuffle elemnts of VT then delete last entries to remove rows. */ for (j = 0; jrearr = realloc(VT->rearr,N*K*sizeof(double)); - MEMALLOCCHK - } + TNT_MEM_CHK /* NO_COVERAGE */ - /* Calculate the truncation error, using the function assigned to the function pointer external variable tnt_matrix_trunc_err() */ - *err = tnt_matrix_trunc_err(S,K,Korig); - - /* S is vector so remove last singular values */ + /* S is real vector so remove last singular double values */ chk = S->rearr = realloc(S->rearr,K*sizeof(double)); - MEMALLOCCHK + TNT_MEM_CHK /* NO_COVERAGE */ + } } - /* Set the size of the arrays */ + /* Re-Set the size of the arrays */ U->sz = M*K; S->sz = K; VT->sz = N*K; - /* Set the new size */ - *Kp = K; - -#ifdef TIMINGS - clock_gettime(CLOCK_REALTIME, &tp_end); - timerE += (double)tp_end.tv_sec-tp_start.tv_sec + 1.0e-9*(tp_end.tv_nsec-tp_start.tv_nsec); -#endif + /* Re-set matrix dimensions */ + U->colsize = K; + VT->rowsize = K; + } return TNT_SUCCESS; } -#ifdef TIMINGS -void tnt_printSVDTimer() -{ - static float cum_tot=0.0, cum_totS=0.0, cum_totE=0.0; - if (0==Parallel.pid) - { - fprintf(costsptr,"\t%8.3f\t%8.3f\t%8.3f", timer, timerS, timerE); - } - cum_tot+=timer; - cum_totS+=timerS; - cum_totE+=timerE; - timer=0.0; - timerS=0.0; - timerE=0.0; - Pprintf("SVD: Timer is %f (Before=%4.2f After=%4.2f)\n", cum_tot, cum_totS, cum_totE); - - return; -} -#endif -