- meld version: '1.5.3'
- source control software type: 'Subversion'
- source control software version: '1.6.17 (r1128011)'
- the output of 'svn diff --diff-cmd diff matrix_svd.c' is attached
- patch command: 'patch -p0 -R -d /tmp/tmpoUE2pV-meld'
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 <omp.h>
+#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
K<M,N. This involves
+ complex matrix A (so A is M x N). K specifies the inner dimensions which can
be the same as or less than M and N. This involves
discarding the last (MIN(M,N)-K) columns of U and last (MIN(M,N)-K)
rows of VT.
If the matrix 'A' is complex it is stored
- as an array of Complex types (a structure definied in the NAG C header
file) so real and imaginary parts of a matrix element
+ as an array of tntComplex types (a structure definied in the NAG C header
file) so real and imaginary parts of a matrix element
are elements of each structure in the array. Similarly all the outputs
are also
- stored as an array of Complex structures. The S array gives the
singular values
+ stored as an array of tntComplex structures. The S array gives the singular
values
of A which are of number min(M,N). The U and VT arrays are 2D matrices
formed with a column-wise convention, i.e. contiguous strips of the 1D
memory
represent the columns of the 2D matrix. The left singular vectors are
returned
@@ -63,78 +54,69 @@
is returned, not V itself.
Memory for the tnsr_values structure S, U and V have been allocated,
but the
array that forms the U->rearr 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; blk<nBlocks; blk++)
- {
+ for (blk=0; blk<nBlocks; blk++) {
tmpI = blkDims[2*(blk+1)]-blkDims[2*blk];
tmpJ = blkDims[2*(blk+1)+1]-blkDims[2*blk+1];
tmpIJ = MIN(tmpI,tmpJ);
- chk = Aloc_C = malloc(tmpI*tmpJ*sizeof(Complex));
MEMALLOCCHKP
- chk = Uloc_C = malloc(tmpI*tmpI*sizeof(Complex));
MEMALLOCCHKP
- chk = VTloc_C = malloc(tmpJ*tmpJ*sizeof(Complex));
MEMALLOCCHKP
- chk = Sloc = malloc(MIN(tmpI,tmpJ)*sizeof(double));
MEMALLOCCHKP
- chk = cwork = malloc(4*sizeof(Complex)); /* Give WORK one entry
initially. */ MEMALLOCCHKP
- chk = rwork = malloc(10*tmpIJ*sizeof(double));
MEMALLOCCHKP
-
- /* Fill Aloc data */
- scaled = tnt_matrix_BlockLocalMatrixFillComplex(A, Aloc_C, blkDims,
blk);
-
- /* Note that zgesvd now being called to return U as MxM and VT as
NxN */
- jobu = 'A';
- jobvt = 'A';
- lwork = -1; /* Queries LAPACK function for the optimal size for
WORK. */
- zgesvd_(&jobu, &jobvt, &tmpI, &tmpJ, Aloc_C, &tmpI, Sloc,
- Uloc_C, &tmpI, VTloc_C, &tmpJ, cwork, &lwork, rwork, &info, 1,
1);
-
- lwork = (int) cwork[0].re; /* Extract optimal WORK size. */
-
- free(cwork);
-
- chk = cwork = (Complex *) malloc(20*lwork*sizeof(Complex)); /*
Assign WORK its optimal size. */ MEMALLOCCHKP
-
- /* Call the LAPACK function properly with optimal settings. */
- zgesvd_(&jobu, &jobvt, &tmpI, &tmpJ, Aloc_C, &tmpI, Sloc,
- Uloc_C, &tmpI, VTloc_C, &tmpJ, cwork, &lwork, rwork, &info, 1,
1);
- /* Copy block-based values back to the NNZ matrix */
- tnt_matrix_blockwiseBlockToNNZRestoreComplex(A, U, VT, S, Uloc_C,
Sloc, VTloc_C, blk, outMatPos, scaled);
+ chk = Aloc_C = malloc(tmpI*tmpJ*sizeof(tntComplex));
+ MEMALLOCCHKP /* NO_COVERAGE */
- free(cwork);
- free(rwork);
+ chk = Uloc_C = malloc(tmpI*tmpI*sizeof(tntComplex));
+ MEMALLOCCHKP /* NO_COVERAGE */
+
+ chk = VTloc_C = malloc(tmpJ*tmpJ*sizeof(tntComplex));
+ MEMALLOCCHKP /* NO_COVERAGE */
+
+ chk = Sloc = malloc(tmpIJ*sizeof(double));
+ MEMALLOCCHKP /* NO_COVERAGE */
+
+ /* Fill Aloc data */
+ scaled = _tnt_matrix_BlockLocalMatrixFillComplex(A,
Aloc_C, blkDims, blk);
+
+ /* Call the SVD function for full SVD (0 for first
argument) */
+ ret = _tnt_la_zsvd(0, tmpI, tmpJ, Aloc_C, Sloc, Uloc_C,
VTloc_C);
+ if (ret != TNT_SUCCESS) {
+ omp_exit = TNT_ERR_GEN; /* NO_COVERAGE */
+ continue; /* NO_COVERAGE */
+ } /* NO_COVERAGE */
+
+ /* Copy block-based values back to the NNZ matrix */
+ _tnt_matrix_blockwiseBlockToNNZRestoreComplex(A, U, VT, S,
Uloc_C, Sloc, VTloc_C, blk, outMatPos, scaled);
free(Uloc_C);
free(VTloc_C);
free(Aloc_C);
@@ -188,82 +162,52 @@
}
#ifdef TNT_OMP
} /* end of parallel OMP section */
- if (omp_exit==1) return TNT_ERR_MEMALLOC; /* Escape if MEMALLOCCHKP
failed */
-/*Pprintf("Currently there are %d threads in the parallel team\n",
omp_get_num_threads());*/
-#endif
-#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
+ /* Escape if anything in OMP section failed */
+ if (TNT_SUCCESS != omp_exit) {
+ TNT_ERR_CMDS /* NO_COVERAGE */
+ return omp_exit; /* NO_COVERAGE */
+ } /* NO_COVERAGE */
+
/* Now process blocked U and VT to recover correctly ordered triplets */
- tnt_matrix_blockwiseNNZFullToNNZRestoreComplex(U, VT, S, outMatPos,
nBlocks);
+ ret = _tnt_matrix_blockwiseNNZFullToNNZRestoreComplex(U, VT, S,
outMatPos, nBlocks);
+ TNT_RET_CHK /* NO_COVERAGE */
} else { /* If only one block */
- /* Try normalising the block, retaining the scaling factor for later
rescaling */
- scaled = tnt_matrix_matrixNormaliserComplex(A->minDimI, 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; i<tmpIJ; i++)
- S->rearr[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; blk<nBlocks; blk++)
- {
+ for (blk=0; blk<nBlocks; blk++) {
tmpI = blkDims[2*(blk+1)]-blkDims[2*blk];
tmpJ = blkDims[2*(blk+1)+1]-blkDims[2*blk+1];
tmpIJ = MIN(tmpI,tmpJ);
- chk = Aloc = malloc(tmpI*tmpJ*sizeof(double));
MEMALLOCCHKP
- chk = Uloc = malloc(tmpI*tmpI*sizeof(double));
MEMALLOCCHKP
- chk = VTloc = malloc(tmpJ*tmpJ*sizeof(double));
MEMALLOCCHKP
- chk = Sloc = malloc(MIN(tmpI,tmpJ)*sizeof(double));
MEMALLOCCHKP
- chk = work = malloc(4*sizeof(double)); /* Give WORK one entry
initially. */ MEMALLOCCHKP
-
- /* Fill Aloc data */
- scaled = tnt_matrix_BlockLocalMatrixFill(A, Aloc, blkDims, blk);
-
- /* Note that zgesvd now being called to return U as MxM and VT as
NxN */
- jobu = 'A';
- jobvt = 'A';
- lwork = -1; /* Queries LAPACK function for the optimal size for
WORK. */
- dgesvd_(&jobu, &jobvt, &tmpI, &tmpJ, Aloc, &tmpI, Sloc,
- Uloc, &tmpI, VTloc, &tmpJ, work, &lwork, &info, 1, 1);
-
- lwork = (int) work[0]; /* Extract optimal WORK size. */
-
- free(work);
-
- chk = work = (double *) malloc(20*lwork*sizeof(double)); /*
Assign WORK its optimal size. */ MEMALLOCCHKP
-
- /* Call the LAPACK function properly with optimal settings. */
- dgesvd_(&jobu, &jobvt, &tmpI, &tmpJ, Aloc, &tmpI, Sloc,
- Uloc, &tmpI, VTloc, &tmpJ, work, &lwork, &info, 1, 1);
+
+ chk = Aloc = malloc(tmpI*tmpJ*sizeof(double));
+ MEMALLOCCHKP /* NO_COVERAGE */
+
+ chk = Uloc = malloc(tmpI*tmpI*sizeof(double));
+ MEMALLOCCHKP /* NO_COVERAGE */
+
+ chk = VTloc = malloc(tmpJ*tmpJ*sizeof(double));
+ MEMALLOCCHKP /* NO_COVERAGE */
+
+ chk = Sloc = malloc(tmpIJ*sizeof(double));
+ MEMALLOCCHKP /* NO_COVERAGE */
+
+ /* Fill Aloc data */
+ scaled = _tnt_matrix_BlockLocalMatrixFill(A, Aloc,
blkDims, blk);
+
+ /* Call the SVD function for full SVD (0 for first
argument) */
+ ret = _tnt_la_dsvd(0, tmpI, tmpJ, Aloc, Sloc, Uloc, VTloc);
+ if (ret != TNT_SUCCESS) {
+ omp_exit = TNT_ERR_GEN; /* NO_COVERAGE */
+ continue; /* NO_COVERAGE */
+ } /* NO_COVERAGE */
/* Copy block-based values back to the NNZ matrix */
- tnt_matrix_blockwiseBlockToNNZRestore(A, U, VT, S, Uloc, Sloc,
VTloc, blk, outMatPos, scaled);
+ _tnt_matrix_blockwiseBlockToNNZRestore(A, U, VT, S, Uloc,
Sloc, VTloc, blk, outMatPos, scaled);
- free(work);
free(Uloc);
free(VTloc);
free(Aloc);
@@ -351,86 +285,59 @@
}
#ifdef TNT_OMP
} /* end of parallel OMP section */
- if (omp_exit==1) return TNT_ERR_MEMALLOC; /* Escape if MEMALLOCCHKP
failed */
-#endif
-#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
+ /* Escape if anything in OMP section failed */
+ if (TNT_SUCCESS != omp_exit) {
+ TNT_ERR_CMDS /* NO_COVERAGE */
+ return omp_exit; /* NO_COVERAGE */
+ } /* NO_COVERAGE */
/* Now process blocked U and VT to recover correctly ordered triplets */
- tnt_matrix_blockwiseNNZFullToNNZRestore(U, VT, S, outMatPos, nBlocks);
+ ret = _tnt_matrix_blockwiseNNZFullToNNZRestore(U, VT, S,
outMatPos, nBlocks);
+ TNT_RET_CHK /* NO_COVERAGE */
- } else { /* If only one block */
+ } else {
- /* Try normalising the block, retaining the scaling factor for later
rescaling */
- scaled = tnt_matrix_matrixNormaliser(A->minDimI, 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; i<tmpIJ; i++)
- S->rearr[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; j<N; j++) {
for (i = 0; i<K; i++) {
VT->comarr[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; j<N; j++) {
@@ -469,51 +391,25 @@
}
}
chk = VT->rearr = 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
-
_______________________________________________
meld-list mailing list
[email protected]
https://mail.gnome.org/mailman/listinfo/meld-list
