From 4ced0156e8d66d86690cb6181fe3af3a1fc36d8b Mon Sep 17 00:00:00 2001 From: Guillaume Horel Date: Tue, 8 Jan 2013 17:10:39 +0100 Subject: [PATCH] Bug #12138 fixed - eigs(A,B) returned incorrect eigenvectors for dense matrices. test_run("arnoldi",["eigs","bug_12138"],["no_check_error_output"]) Change-Id: Ic4d02e9f6653ffab25c15b7c10b09d76e3be74e7 --- scilab/modules/arnoldi/includes/eigs.h | 11 +- .../modules/arnoldi/includes/eigs_dependencies.h | 76 +- scilab/modules/arnoldi/sci_gateway/c/sci_eigs.c | 171 ++- scilab/modules/arnoldi/src/c/eigs.c | 1105 ++++++-------------- scilab/modules/arnoldi/src/c/eigs_dependencies.c | 315 ++---- .../arnoldi/tests/nonreg_tests/bug_12138.dia.ref | 81 ++ .../arnoldi/tests/nonreg_tests/bug_12138.tst | 90 ++ 7 files changed, 770 insertions(+), 1079 deletions(-) create mode 100644 scilab/modules/arnoldi/tests/nonreg_tests/bug_12138.dia.ref create mode 100644 scilab/modules/arnoldi/tests/nonreg_tests/bug_12138.tst diff --git a/scilab/modules/arnoldi/includes/eigs.h b/scilab/modules/arnoldi/includes/eigs.h index 1586115..469cc7a 100644 --- a/scilab/modules/arnoldi/includes/eigs.h +++ b/scilab/modules/arnoldi/includes/eigs.h @@ -14,6 +14,7 @@ #define __EIGS_H__ #include "doublecomplex.h" + /** * @TODO add comment * @@ -39,15 +40,17 @@ * @param INFO_EUPD * @param eigenvalue * @param eigenvector + * @param eigenvalueC + * @param eigenvectorC + * @param RVEC * @return */ int eigs(double *AR, doublecomplex *AC, int N, int Acomplex, int Asym, double* B, doublecomplex* BC, int Bcomplex, int matB, int nev, - doublecomplex* SIGMA, char* which, double* maxiter, double* tol, + doublecomplex SIGMA, char* which, double* maxiter, double* tol, double* NCV, double* RESID, doublecomplex* RESIDC, int* INFO, - double* cholB, int INFO_EUPD, doublecomplex* eigenvalue, - doublecomplex* eigenvector); + double* cholB, int INFO_EUPD, double* eigenvalue, double* eigenvector, + doublecomplex* eigenvalueC, doublecomplex* eigenvectorC, int RVEC); #endif /* __EIGS_H__ */ /*--------------------------------------------------------------------------*/ - diff --git a/scilab/modules/arnoldi/includes/eigs_dependencies.h b/scilab/modules/arnoldi/includes/eigs_dependencies.h index 9be06ff..e770488 100644 --- a/scilab/modules/arnoldi/includes/eigs_dependencies.h +++ b/scilab/modules/arnoldi/includes/eigs_dependencies.h @@ -10,8 +10,8 @@ * */ /*--------------------------------------------------------------------------*/ -#ifndef __RTIMESRprime_H__ -#define __RTIMESRprime_H__ +#ifndef __PROCESS_DNEUPD_H__ +#define __PROCESS_DNEUPD_H__ #include "doublecomplex.h" #include #include @@ -24,67 +24,21 @@ /** * @TODO add comment * - * @param result - * @param R - * @param Rprime + * @param DR + * @param DI + * @param Z * @param N + * @param nev + * @param AR + * @param eigenvalue + * @param eigenvector + * @param sigma_imaginary + * @param RVEC */ -void RtimesRprime(double* result, double* R, double* Rprime, int N); -/** - * @TODO add comment - * - * @param result - * @param R - * @param A - * @param Rprime - * @param N - */ -void invR_times_A_times_invRprime(double* result, double* R, double* A, double* Rprime, int N); - -/** - * @TODO add comment - * - * @param result - * @param U - * @param L - * @param E - * @param N - */ -void invU_times_invL_times_E(double* result, double* U, double* L, double* E, int N); - -/** - * @TODO add comment - * - * @param result - * @param R - * @param Rprime - * @param N - */ -void RCtimesRCprime(doublecomplex* result, doublecomplex* R, doublecomplex* Rprime, int N); - -/** - * @TODO add comment - * - * @param result - * @param R - * @param A - * @param Rprime - * @param N - */ -void invRC_times_AC_times_invRCprime(doublecomplex* result, doublecomplex* R, doublecomplex* A, doublecomplex* Rprime, int N); +void process_dneupd_data(double* DR, double* DI, double* Z, int N, int nev, double* AR, + doublecomplex* eigenvalue, doublecomplex* eigenvector, + int sigma_imaginary); -/** - * @TODO add comment - * - * @param result - * @param U - * @param L - * @param E - * @param N - */ -void invUC_times_invLC_times_EC(doublecomplex* result, doublecomplex* U, doublecomplex* L, doublecomplex* E, int N); - -#endif /* __RTIMESRprime_H__ */ +#endif /* __PROCESS_DNEUPD_H__ */ /*--------------------------------------------------------------------------*/ - diff --git a/scilab/modules/arnoldi/sci_gateway/c/sci_eigs.c b/scilab/modules/arnoldi/sci_gateway/c/sci_eigs.c index e88aa82..2f96351 100644 --- a/scilab/modules/arnoldi/sci_gateway/c/sci_eigs.c +++ b/scilab/modules/arnoldi/sci_gateway/c/sci_eigs.c @@ -58,7 +58,7 @@ int sci_eigs(char *fname, unsigned long fname_len) int iColsFour = 0; int iLen = 0; char* pstData = NULL; - doublecomplex* SIGMA = NULL; + doublecomplex SIGMA; int *piAddressVarFive = NULL; double dblMAXITER = 0; @@ -86,11 +86,15 @@ int sci_eigs(char *fname, unsigned long fname_len) int *piAddressVarTen = NULL; int iINFO = 0; - + int RVEC = 0; // Output arguments - doublecomplex* eigenvalue = NULL; - doublecomplex* mat_eigenvalue = NULL; - doublecomplex* eigenvector = NULL; + double* eigenvalue = NULL; + double* eigenvector = NULL; + doublecomplex* eigenvalueC = NULL; + doublecomplex* eigenvectorC = NULL; + + double* mat_eigenvalue = NULL; + doublecomplex* mat_eigenvalueC = NULL; int INFO_EUPD = 0; int error = 0; @@ -173,13 +177,13 @@ int sci_eigs(char *fname, unsigned long fname_len) } sciErr = getVarDimension(pvApiCtx, piAddressVarTwo, &iRowsTwo, &iColsTwo); - if (iRowsTwo * iColsTwo == 1 || iRowsTwo != iColsTwo) + matB = iRowsTwo * iColsTwo; + if (matB && (iRowsTwo != iRowsOne || iColsTwo != iColsOne)) { Scierror(999, _("%s: Wrong dimension for input argument #%d: B must have the same size as A.\n"), "eigs", 2); return 0; } - matB = iRowsTwo * iColsTwo; if (isVarComplex(pvApiCtx, piAddressVarTwo)) { sciErr = getComplexZMatrixOfDouble(pvApiCtx, piAddressVarTwo, &iRowsTwo, &iColsTwo, &Bcplx); @@ -310,9 +314,8 @@ int sci_eigs(char *fname, unsigned long fname_len) return 0; } - SIGMA = (doublecomplex*)malloc(1 * sizeof(doublecomplex)); - SIGMA[0].r = 0; - SIGMA[0].i = 0; + SIGMA.r = 0; + SIGMA.i = 0; } if (iTypeVarFour == sci_matrix) @@ -324,15 +327,14 @@ int sci_eigs(char *fname, unsigned long fname_len) return 0; } - SIGMA = (doublecomplex*)MALLOC(1 * sizeof(doublecomplex)); - if (getScalarComplexDouble(pvApiCtx, piAddressVarFour, &SIGMA[0].r, &SIGMA[0].i)) + if (getScalarComplexDouble(pvApiCtx, piAddressVarFour, &SIGMA.r, &SIGMA.i)) { printError(&sciErr, 0); Scierror(999, _("%s: Can not read input argument #%d.\n"), fname, 4); return 0; } - if (C2F(isanan)(&SIGMA[0].r) || C2F(isanan)(&SIGMA[0].i)) + if (C2F(isanan)(&SIGMA.r) || C2F(isanan)(&SIGMA.i)) { Scierror(999, _("%s: Wrong type for input argument #%d: sigma must be a real.\n"), "eigs", 4); return 0; @@ -349,7 +351,6 @@ int sci_eigs(char *fname, unsigned long fname_len) { printError(&sciErr, 0); Scierror(999, _("%s: Can not read input argument #%d.\n"), fname, 5); - free(SIGMA); return 0; } @@ -494,6 +495,54 @@ int sci_eigs(char *fname, unsigned long fname_len) } } + if ( dblCHOLB ) // check that B is upper triangular with non zero element on the diagonal + { + if (!Bcomplex) + { + for (i = 0; i < N; i++) + { + for (j = 0; j <= i; j++) + { + if (i == j && Breal[i + j * N] == 0) + { + Scierror(999, _("%s: B is not positive definite. Try with sigma='SM' or sigma=scalar.\n"), "eigs"); + return 0; + } + else + { + if ( j < i && Breal[i + j * N] != 0 ) + { + Scierror(999, _("%s: If opts.cholB is true, B should be upper triangular.\n"), "eigs"); + return 0; + } + } + } + } + } + else + { + for (i = 0; i < N; i++) + { + for (j = 0; j <= i; j++) + { + if (i == j && Bcplx[i + i * N].r == 0 && Bcplx[i + i * N].i == 0) + { + Scierror(999, _("%s: B is not positive definite. Try with sigma='SM' or sigma=scalar.\n"), "eigs"); + return 0; + } + else + { + if ( j < i && (Bcplx[i + j * N].r != 0 || Bcplx[i + j * N].i != 0) ) + { + Scierror(999, _("%s: If opts.cholB is true, B should be upper triangular.\n"), "eigs"); + return 0; + } + } + } + } + } + } + /**************************************** * RESID * *****************************************/ @@ -570,20 +619,29 @@ int sci_eigs(char *fname, unsigned long fname_len) } // Initialization output arguments - eigenvalue = (doublecomplex*) malloc((iNEV + 1) * sizeof(doublecomplex)); - memset(eigenvalue, 0, (iNEV + 1) * sizeof(doublecomplex)); - if (Lhs > 1) { - mat_eigenvalue = (doublecomplex*) malloc((iNEV + 1) * (iNEV + 1) * sizeof(doublecomplex)); - memset(mat_eigenvalue, 0, (iNEV + 1) * (iNEV + 1) * sizeof(doublecomplex)); + RVEC = 1; + } - eigenvector = (doublecomplex*) malloc(N * (iNEV + 1) * sizeof(doublecomplex)); - memset(eigenvector, 0, N * (iNEV + 1) * sizeof(doublecomplex)); + if (Acomplex || Bcomplex || !Asym) + { + eigenvalueC = (doublecomplex*)CALLOC((iNEV + 1), sizeof(doublecomplex)); + if (RVEC) + { + eigenvectorC = (doublecomplex*)CALLOC(N * (iNEV + 1), sizeof(doublecomplex)); + } + } + else + { + eigenvalue = (double*)CALLOC(iNEV, sizeof(double)); + /* we should allocate eigenvector only if RVEC is true, but dseupd segfaults + if Z is not allocated even when RVEC is false, contrary to the docs.*/ + eigenvector = (double*)CALLOC(iNEV * N, sizeof(double)); } error = eigs(Areal, Acplx, N, Acomplex, Asym, Breal, Bcplx, Bcomplex, matB, iNEV, SIGMA, pstData, &dblMAXITER, - &dblTOL, dblNCV, RESID, RESIDC, &iINFO, &dblCHOLB, INFO_EUPD, eigenvalue, eigenvector); + &dblTOL, dblNCV, RESID, RESIDC, &iINFO, &dblCHOLB, INFO_EUPD, eigenvalue, eigenvector, eigenvalueC, eigenvectorC, RVEC); switch (error) { @@ -619,8 +677,8 @@ int sci_eigs(char *fname, unsigned long fname_len) return 0; case -3 : - Scierror(999, _("%s: Wrong type for input argument #%d: B must be symmetric or hermitian, definite, semi positive.\n"), "eigs", 2); - PutLhsVar(); + Scierror(999, _("%s: Error with input argument #%d: B is not positive definite. Try with sigma='SM' or sigma=scalar.\n"), "eigs", 2); + ReturnArguments(pvApiCtx); return 0; case -4 : @@ -675,16 +733,30 @@ int sci_eigs(char *fname, unsigned long fname_len) } else { - Scierror(999, _("%s: Error with %s: info = %d \n"), "eigs", "ZNEUPD", INFO_EUPD); + Scierror(999, _("%s: Error with %s: info = %d \n"), "eigs", "ZNEUPD", INFO_EUPD); } PutLhsVar(); - free(mat_eigenvalue); + return 0; + + case -7 : + Scierror(999, _("%s: A - sigma * B is not inversible, try with a different value of sigma.\n"), "eigs"); + PutLhsVar(); return 0; } if (Lhs <= 1) { - sciErr = createComplexZMatrixOfDouble(pvApiCtx, Rhs + 1, iNEV, 1, eigenvalue); + if (eigenvalue) + { + sciErr = createMatrixOfDouble(pvApiCtx, Rhs + 1, iNEV, 1, eigenvalue); + FREE(eigenvalue); + FREE(eigenvector); + } + else if (eigenvalueC) + { + sciErr = createComplexZMatrixOfDouble(pvApiCtx, Rhs + 1, iNEV, 1, eigenvalueC); + FREE(eigenvalueC); + } if (sciErr.iErr) { printError(&sciErr, 0); @@ -697,13 +769,29 @@ int sci_eigs(char *fname, unsigned long fname_len) else { // create a matrix which contains the eigenvalues - for (i = 0; i < iNEV ; i++) + if (eigenvalue) + { + mat_eigenvalue = (double*)CALLOC(iNEV * iNEV, sizeof(double)); + for (i = 0; i < iNEV; i++) + { + mat_eigenvalue[i * iNEV + i] = eigenvalue[i]; + } + sciErr = createMatrixOfDouble(pvApiCtx, Rhs + 1, iNEV, iNEV, mat_eigenvalue); + FREE(eigenvalue); + FREE(mat_eigenvalue); + } + else if (eigenvalueC) { - mat_eigenvalue[i * iNEV + i].r = eigenvalue[i].r; - mat_eigenvalue[i * iNEV + i].i = eigenvalue[i].i; + mat_eigenvalueC = (doublecomplex*)CALLOC(iNEV * iNEV, sizeof(doublecomplex)); + for (i = 0; i < iNEV; i++) + { + mat_eigenvalueC[i * iNEV + i] = eigenvalueC[i]; + } + sciErr = createComplexZMatrixOfDouble(pvApiCtx, Rhs + 1, iNEV, iNEV, mat_eigenvalueC); + FREE(eigenvalueC); + FREE(mat_eigenvalueC); } - sciErr = createComplexZMatrixOfDouble(pvApiCtx, Rhs + 1, iNEV, iNEV, mat_eigenvalue); if (sciErr.iErr) { printError(&sciErr, 0); @@ -711,7 +799,16 @@ int sci_eigs(char *fname, unsigned long fname_len) return 0; } - sciErr = createComplexZMatrixOfDouble(pvApiCtx, Rhs + 2, N, iNEV, eigenvector); + if (eigenvector) + { + sciErr = createMatrixOfDouble(pvApiCtx, Rhs + 2, N, iNEV, eigenvector); + FREE(eigenvector); + } + else if (eigenvectorC) + { + sciErr = createComplexZMatrixOfDouble(pvApiCtx, Rhs + 2, N, iNEV, eigenvectorC); + FREE(eigenvectorC); + } if (sciErr.iErr) { printError(&sciErr, 0); @@ -728,10 +825,6 @@ int sci_eigs(char *fname, unsigned long fname_len) freeAllocatedSingleString(pstData); } - free(SIGMA); - - free(eigenvalue); - if (matB != 0) { if (Acomplex && !Bcomplex) @@ -744,12 +837,6 @@ int sci_eigs(char *fname, unsigned long fname_len) } } - if (Lhs > 1) - { - free(mat_eigenvalue); - free(eigenvector); - } PutLhsVar(); return 0; } - diff --git a/scilab/modules/arnoldi/src/c/eigs.c b/scilab/modules/arnoldi/src/c/eigs.c index 1e78840..e42e9c5 100644 --- a/scilab/modules/arnoldi/src/c/eigs.c +++ b/scilab/modules/arnoldi/src/c/eigs.c @@ -23,21 +23,29 @@ /*--------------------------------------------------------------------------*/ // dgemm performs one of the matrix-matrix operations -extern int C2F(dgemm)(char* transa, char* transb, int* m, int* n, int* k, double* alpha, double* A, int* lda, double* B, int* ldb, double* beta, double* C, int* ldc); - +extern int C2F(dgemm)(char* transa, char* transb, int* m, int* n, int* k, + double* alpha, double* A, int* lda, double* B, int* ldb, + double* beta, double* C, int* ldc); // zgemm performs one of the matrix-matrix operations -extern int C2F(zgemm)(char* transa, char* transb, int* m, int* n, int* k, doublecomplex* alpha, doublecomplex* A, int* lda, - doublecomplex* B, int* ldb, doublecomplex* beta, doublecomplex* C, int* ldc); +extern int C2F(zgemm)(char* transa, char* transb, int* m, int* n, int* k, + doublecomplex* alpha, doublecomplex* A, int* lda, doublecomplex* B, + int* ldb, doublecomplex* beta, doublecomplex* C, int* ldc); + +// dgemv performs matrix-vector operations +extern int C2F(dgemv) (char* trans, int* m, int* n, double* alpha, double* A, int* lda, + double* x, int* incx, double* beta, double* y, int* incy); +extern int C2F(zgemv) (char* trans, int* m, int* n, doublecomplex* alpha, doublecomplex* A, + int* lda, doublecomplex* x, int* incx, doublecomplex* beta, doublecomplex* y, int* incy); // dgetrf computes an LU factorization of a general M by N matrix A (double) using partial pivoting with row interchanges extern int C2F(dgetrf)(int* m, int* n, double* A, int* lda, int* ipiv, int* info); -// zgetrd computes an LU factorization of a general M by N matrix A (complex*16) using partial pivoting with row interchanges +// zgetrf computes an LU factorization of a general M by N matrix A (complex*16) using partial pivoting with row interchanges extern int C2F(zgetrf)(int* m, int* n, doublecomplex* A, int* lda, int* ipiv, int* info); - -// dlaswp performs a series of row interchanges on the matrix A -// one row interchange is initiated for each of rows k1 through k2 of A -extern int C2F(dlaswp)(int* n, double* A, int* lda, int* k1, int* k2, int* ipiv, int* incx); +// dgetrs solves a linear system using the factors computed by dgetrf +extern int C2F(dgetrs) (char* trans, int* n, int* nrhs, double* A, int *lda, int* ipiv, double* B, int* ldb, int* info); +// zgetrs solves a linear system using the factors computed by zgetrf +extern int C2F(zgetrs) (char* trans, int* n, int* nrhs, doublecomplex* AC, int* lda, int* ipiv, doublecomplex* B, int* ldb, int* info); // dpotrf computes the cholesky factorization of a real symmetric positive definite matrix A extern int C2F(dpotrf)(char* uplo, int* n, double* A, int* lda, int* info); @@ -45,11 +53,20 @@ extern int C2F(dpotrf)(char* uplo, int* n, double* A, int* lda, int* info); // zpotrf computes the cholesky factorization of a real hermitian positive definite matrix A extern int C2F(zpotrf)(char* uplo, int* n, doublecomplex* A, int* lda, int* info); -// dgetri computes the inverse of a matrix using the LU factorization computed by dgetrf -extern int C2F(dgetri)(int* n, double* A, int* lda, int* ipiv, double* work, int* lworkl, int* info); - -// zgetri computes the inverse of a matrix using the LU factorization computed by zgetrf -extern int C2F(zgetri)(int* n, doublecomplex* A, int* lda, int* ipiv, doublecomplex* work, int* lworkl, int* info); +// dtrsm solves a triangular linear system +extern int C2F(dtrsm) (char* side, char* uplo, char* trans, char* diag, int* m, int* n, double* alpha, double* A, int* lda, double* B, int* ldb); +// ztrsm solve a triangular linear system +extern int C2F(ztrsm) (char* side, char* uplo, char* trans, char* diag, int* m, int* n, doublecomplex* alpha, doublecomplex* A, int* lda, doublecomplex* B, int* ldb); +// dsyrk does a rank k symmetric update +extern int C2F(dsyrk) (char* uplo, char* trans, int* n, int* k, double* alpha, + double* A, int* lda, double* beta, double* B, int* ldb); +// ztrmm multiply by a triangular matrix +extern int C2F(ztrmm) (char* side, char* uplo, char* trans, char* diag, int* m, int* n, doublecomplex* alphac, + doublecomplex* A, int* lda, doublecomplex* B, int* ldb); +// ztrmv multiply a vector by a triangular matrix +extern int C2F(ztrmv) (char* uplo, char* trans, char* diag, int* n, doublecomplex* A, int* lda, doublecomplex* x, int* incx); +// dtrmv multiply a vector by a triangular matrix +extern int C2F(dtrmv) (char* uplo, char* trans, char* diag, int* n, double* A, int* lda, double* x, int* incx); /*--------------------------------------------------------------------------*/ /*--------------------------------------------------------------------------*/ @@ -65,17 +82,14 @@ extern int C2F(dseupd)(int *rvec, char *howmny, int *select, double *d, int *n, char *which, int *nev , double *tol, double *resid, int *ncv, double *v , int *ldv, int *iparam, int *ipntr, double *workd, double *workl, - int *lworkl, int *info, unsigned long rvec_length, - unsigned long howmany_length, - unsigned long bmat_length, unsigned long which_len); + int *lworkl, int *info); /*--------------------------------------------------------------------------*/ /*--------------------------------------------------------------------------*/ extern int C2F(dnaupd)(int *ido, char *bmat, int *n, char *which, int *nev, double *tol, double *resid, int *ncv, double *v, int *ldv, int *iparam, int *ipntr, double *workd, - double *workl, int *lworkl, int *info, - unsigned long bmat_len, unsigned long which_len); + double *workl, int *lworkl, int *info); /*--------------------------------------------------------------------------*/ /*--------------------------------------------------------------------------*/ @@ -113,13 +127,16 @@ static double beta = 0.; static doublecomplex alphac = {1., 0.}; static doublecomplex betac = {0., 0.}; -int eigs(double *AR, doublecomplex *AC, int N, int Acomplex, int Asym, double* B, - doublecomplex* BC, int Bcomplex, int matB, int nev, doublecomplex* SIGMA, - char* which, double* maxiter, double* tol, double* NCV, double* RESID, doublecomplex* RESIDC, - int* INFO, double* cholB, int INFO_EUPD, doublecomplex* eigenvalue, doublecomplex* eigenvector) +int eigs(double *AR, doublecomplex *AC, int N, int Acomplex, int Asym, + double* B, doublecomplex* BC, int Bcomplex, int matB, int nev, + doublecomplex SIGMA, char* which, double* maxiter, double* tol, + double* NCV, double* RESID, doublecomplex* RESIDC, int* INFO, + double* cholB, int INFO_EUPD, double* eigenvalue, + double* eigenvector, doublecomplex* eigenvalueC, + doublecomplex* eigenvectorC, int RVEC) + { - int index = 0; // GENERAL VARIABLES int i = 0; int j = 0; @@ -127,7 +144,7 @@ int eigs(double *AR, doublecomplex *AC, int N, int Acomplex, int Asym, double* B int l = 0; int INFO_CHOL = 0; int INFO_LU = 0; - int k1 = 1; + int INFO_INV = 0; int iOne = 1; // VARIABLES DSAUPD, DNAUPD, ZNAUPD @@ -136,8 +153,8 @@ int eigs(double *AR, doublecomplex *AC, int N, int Acomplex, int Asym, double* B int LDV = Max(1, N); int ncv = 0; - int* IPARAM = NULL; - int* IPNTR = NULL; + int IPARAM[11]; + int IPNTR[14]; double* V = NULL; doublecomplex* VC = NULL; @@ -153,86 +170,48 @@ int eigs(double *AR, doublecomplex *AC, int N, int Acomplex, int Asym, double* B char* bmat = "I"; // VARIABLES DSEUPD, DNEUPD, ZNEUPD - int RVEC = 0; // compute eigenvalues if RVEC = 1 also compute eigenvalues and eigenvectors char* HOWMNY = "A"; int* SELECT = NULL; - double* D = NULL; double* DI = NULL; double* DR = NULL; - doublecomplex* DC = NULL; + double* Z = NULL; double* WORKEV = NULL; doublecomplex* WORKEVC = NULL; - double* Z = NULL; - doublecomplex* ZC = NULL; - - double SIGMAR = SIGMA[0].r; - double SIGMAI = SIGMA[0].i; + doublecomplex mSIGMA = {.r = -SIGMA.r, .i = -SIGMA.i }; - double* R = (double*)malloc(N * N * sizeof(double)); - double* Rprime = (double*)malloc(N * N * sizeof(double)); + double* R = NULL; + doublecomplex* RC = NULL; double* AMSB = NULL; doublecomplex* AMSBC = NULL; - double* L = NULL; - double* U = NULL; - double* E = NULL; - - doublecomplex* RC = (doublecomplex*)malloc(N * N * sizeof(doublecomplex)); - doublecomplex* RCprime = (doublecomplex*)malloc(N * N * sizeof(doublecomplex)); - - doublecomplex* LC = NULL; - doublecomplex* UC = NULL; - doublecomplex* EC = NULL; - - double* tmp_WORKD = NULL; - doublecomplex* tmp_WORKDC = NULL; - - double* R_Rprime = NULL; - double* invR_A_invRprime = NULL; - double* invU_invL_E = NULL; + int* IPVT = NULL; - doublecomplex* RC_RCprime = NULL; - doublecomplex* invRC_AC_invRCprime = NULL; - doublecomplex* invUC_invLC_EC = NULL; + double* temp = NULL; + doublecomplex* tempC = NULL; - int* IPVT = NULL; + int oldnev = nev; + int N2 = N * N; - IPARAM = (int*)malloc(11 * sizeof(int)); - memset(IPARAM, 0, 11 * sizeof(int)); IPARAM[0] = 1; IPARAM[2] = (int) maxiter[0]; IPARAM[6] = 1; // by default mode = 1 - IPNTR = (int*)malloc(14 * sizeof(int)); - memset(IPNTR, 0, 14 * sizeof(int)); - - tmp_WORKD = (double*)malloc(N * sizeof(double)); - memset(tmp_WORKD, 0, N * sizeof(double)); - tmp_WORKDC = (doublecomplex*)malloc(N * sizeof(doublecomplex)); - memset(tmp_WORKDC, 0, N * sizeof(doublecomplex)); - // END VARIABLES - // Info to compute eigenvalues and eigenvectors - if (eigenvector != NULL) - { - RVEC = 1; - } - // MODE - if (!strcmp(which, "SM") || (SIGMAR != 0 || SIGMAI != 0)) + if (!strcmp(which, "SM") || (SIGMA.r != 0 || SIGMA.i != 0)) { IPARAM[6] = 3; which = "LM"; } // BMAT - if ((matB == 0) || (IPARAM[6] == 1)) // if B = [] or mode = 1 -> bmat = 'I' : standart eigenvalue problem + if ((matB == 0) || (IPARAM[6] == 1)) // if B = [] or mode = 1 -> bmat = 'I' : standard eigenvalue problem { bmat = "I"; } @@ -262,198 +241,50 @@ int eigs(double *AR, doublecomplex *AC, int N, int Acomplex, int Asym, double* B ncv = (int) NCV[0]; if (ncv <= nev || ncv > N) // Error { - free(IPARAM); - free(IPNTR); - free(R); - free(Rprime); - free(RC); - free(RCprime); - free(tmp_WORKD); - free(tmp_WORKDC); return -1; - } } // NEV - if ((!Acomplex && !Bcomplex && Asym == 1 && nev >= N) || (!(!Acomplex && !Bcomplex && Asym == 1) && nev >= N - 1)) + if ((!Acomplex && !Bcomplex && Asym == 1 && nev >= N) || ((Acomplex || Bcomplex || !Asym) && nev >= N - 1)) { - free(IPARAM); - free(IPNTR); - free(R); - free(Rprime); - free(RC); - free(RCprime); - free(tmp_WORKD); - free(tmp_WORKDC); return -2; } - // B must be symmetric (hermitian) positive (semi-) positive if (matB != 0) { - if (cholB[0]) // Comparison between B and upper triangular matrix + if (cholB[0]) // we already have the cholesky decomposition { - if (!Bcomplex) // B is real - { - for (i = 0 ; i < N ; i++) - { - for (j = i + 1 ; j < N ; j++) - { - if (B[j + i * N] != 0) - { - free(IPARAM); - free(IPNTR); - free(R); - free(Rprime); - free(RC); - free(RCprime); - free(tmp_WORKD); - free(tmp_WORKDC); - free(IPVT); - return -3; - } - } - } - - memcpy(Rprime, B, N * N * sizeof(double)); - - // Compute the lower triangular matrix - memset(R, 0, N * N * sizeof(double)); - for (i = 0 ; i < N ; i++) - { - for (j = 0 ; j < N ; j++) - { - R[i * N + j] = B[j * N + i]; - } - } - } - else // if B is complex + R = B; + RC = BC; + } + else + { + if (IPARAM[6] == 1) { - for (i = 0 ; i < N ; i++) + if (!Bcomplex) // B is real { - for (j = i + 1 ; j < N ; j++) + R = (double *)malloc(N * N * sizeof(double)); + memcpy(R, B, N * N * sizeof(double)); + C2F(dpotrf) ("u", &N, R, &N, &INFO_CHOL); // Compute the upper triangular matrix R + if (INFO_CHOL != 0) // Errors { - if (BC[j + i * N].r != 0 || BC[j + i * N].i != 0) - { - free(IPARAM); - free(IPNTR); - free(R); - free(Rprime); - free(RC); - free(RCprime); - free(tmp_WORKD); - free(tmp_WORKDC); - return -3; - } + free(R); + return -3; } } - - memcpy(RCprime, BC, N * N * sizeof(doublecomplex)); - - // Computes the lower triangular matrix BC - memset(RC, 0, N * N * sizeof(doublecomplex)); - for (i = 0 ; i < N ; i++) + else // B is complex { - for (j = 0 ; j < N ; j++) + RC = (doublecomplex *) malloc(N * N * sizeof(doublecomplex)); + memcpy(RC, BC, N * N * sizeof(doublecomplex)); + C2F(zpotrf) ("u", &N, RC, &N, &INFO_CHOL); // Computes the upper triangular matrix + if (INFO_CHOL != 0) { - RC[i * N + j].r = BC[j * N + i].r; - RC[i * N + j].i = (-1) * BC[j * N + i].i; + free(RC); + return -3; } } } - - } - } - - if (!cholB[0] && IPARAM[6] == 1 && matB != 0) - { - if (!Bcomplex) // B is real - { - memcpy(R, B, N * N * sizeof(double)); - memcpy(Rprime, B, N * N * sizeof(double)); - - C2F(dpotrf)("L", &N, R, &N, &INFO_CHOL); // Compute the lower triangular matrix R - if (INFO_CHOL != 0) // Errors - { - free(IPARAM); - free(IPNTR); - free(R); - free(Rprime); - free(RC); - free(RCprime); - free(tmp_WORKD); - free(tmp_WORKDC); - return -3; - } - - C2F(dpotrf)("U", &N, Rprime, &N, &INFO_CHOL); // Compute the upper triangular matrix Rprime - if (INFO_CHOL != 0) - { - free(IPARAM); - free(IPNTR); - free(R); - free(Rprime); - free(RC); - free(RCprime); - free(tmp_WORKD); - free(tmp_WORKDC); - return -3; - } - - for (j = 0 ; j < N ; j++) - { - for (i = 0 ; i < j ; i++) - { - R[i + j * N] = 0; - Rprime[j + i * N] = 0; - } - } - } - else // B is complex - { - memcpy(RC, BC, N * N * sizeof(doublecomplex)); - memcpy(RCprime, BC, N * N * sizeof(doublecomplex)); - - C2F(zpotrf)("L", &N, RC, &N, &INFO_CHOL); // Computes the lower triangular matrix - if (INFO_CHOL != 0) - { - free(IPARAM); - free(IPNTR); - free(R); - free(Rprime); - free(RC); - free(RCprime); - free(tmp_WORKD); - free(tmp_WORKDC); - return -3; - } - - C2F(zpotrf)("U", &N, RCprime, &N, &INFO_CHOL); // Computes the upper triangular matrix - if (INFO_CHOL != 0) - { - free(IPARAM); - free(IPNTR); - free(R); - free(Rprime); - free(RC); - free(RCprime); - free(tmp_WORKD); - free(tmp_WORKDC); - return -3; - } - - for (j = 0 ; j < N ; j++) - { - for (i = 0 ; i < j ; i++) - { - RC[i + j * N].r = 0; - RC[i + j * N].i = 0; - - RCprime[j + i * N].r = 0; - RCprime[j + i * N].i = 0; - } - } } } @@ -464,106 +295,70 @@ int eigs(double *AR, doublecomplex *AC, int N, int Acomplex, int Asym, double* B { AMSB = (double*)malloc(N * N * sizeof(double)); memcpy(AMSB, AR, N * N * sizeof(double)); - - // Compute LU decomposition AMSB = A - sigma*B - if (matB == 0) // if B = [] -> standart eigenvalue problem : A - sigma + if (SIGMA.r != 0) { - for (i = 0 ; i < N ; i++) + // Compute LU decomposition AMSB = A - sigma*B + if (matB == 0) // if B = [] -> standard eigenvalue problem : A - sigma *I { - AMSB[i + i * N] = AMSB[i + i * N] - SIGMAR; - } - } - else // generalized eigenvalue problem - { - if (cholB[0]) - { - if (R_Rprime == NULL) + for (i = 0 ; i < N ; i++) { - R_Rprime = (double*)malloc(N * N * sizeof(double)); - RtimesRprime(R_Rprime, R, Rprime, N); - } - - for (i = 0 ; i < N * N ; i++) - { - AMSB[i] = AR[i] - (SIGMAR * R_Rprime[i]); + AMSB[i + i * N] -= SIGMA.r; } } - else + else // generalized eigenvalue problem { - for (i = 0 ; i < N * N ; i++) + if (cholB[0]) { - AMSB[i] = AR[i] - (SIGMAR * B[i]); + C2F(dsyrk) ("u", "t", &N, &N, &mSIGMA.r, R, &N, &alpha, AMSB, &N); + if (!Asym) //dsyrk does a symmetric update so we need to correct for the antisymmetric part + { + for (i = 0; i < N; i++) + { + for (j = 0; j < i; j++) + { + AMSB[i + j * N] = AMSB[j + i * N] + AR[i + j * N] - AR[j + i * N]; + } + } + } + } + else + { + C2F(daxpy)(&N2, &mSIGMA.r, B, &iOne, AMSB, &iOne); } } } // LU decomposition - IPVT = (int*) malloc(N * sizeof(int)); - memset(IPVT, 0, N * sizeof(int)); + IPVT = (int*) calloc(N, sizeof(int)); C2F(dgetrf)(&N, &N, AMSB, &N, IPVT, &INFO_LU); - - // Computes the lower triangular matrix L - L = (double*)malloc(N * N * sizeof(double)); - memset(L, 0, N * N * sizeof(double)); - - for (i = 0 ; i < N ; i++) + if (INFO_LU > 0) { - for (j = 0 ; j < i ; j++) - { - L[i + j * N] = AMSB[i + j * N]; - } - - L[i + i * N] = 1; + free(IPVT); + free(AMSB); + return -7; } - - // Computes the upper triangular matrix U - U = (double*)malloc(N * N * sizeof(double)); - memset(U, 0, N * N * sizeof(double)); - - for (j = 0 ; j < N ; j++) - { - for (i = 0 ; i <= j ; i++) - { - //if(i <= j) - U[i + j * N] = AMSB[i + j * N]; - } - } - - // Computes the permutation matrix E - E = (double*) malloc(N * N * sizeof(double)); - memset(E, 0, N * N * sizeof(double)); - - for (i = 0 ; i < N ; i++) - { - E[i * N + i] = 1; - } - - C2F(dlaswp)(&N, E, &N, &k1, &N, IPVT, &k1); - - free(AMSB); - free(IPVT); } if (Asym) // DSAUPD { LWORKL = ncv * ncv + 8 * ncv; + WORKL = (double*) calloc(LWORKL, sizeof(double)); - WORKL = (double*)malloc(LWORKL * sizeof(double)); - memset(WORKL, 0, LWORKL * sizeof(double)); } else // DNAUPD { LWORKL = 3 * ncv * (ncv + 2); + WORKL = (double*) calloc(LWORKL, sizeof(double)); - WORKL = (double*)malloc(LWORKL * sizeof(double)); - memset(WORKL, 0, LWORKL * sizeof(double)); } - WORKD = (double*)malloc(3 * N * sizeof(double)); - memset(WORKD, 0, 3 * N * sizeof(double)); + WORKD = (double*) calloc(3 * N, sizeof(double)); + V = (double*) calloc(N * ncv, sizeof(double)); - V = (double*)malloc(N * ncv * sizeof(double)); - memset(V, 0, N * ncv * sizeof(double)); + if (IPARAM[6] == 1 && matB) + { + temp = (double*) malloc(N * sizeof(double)); + } while (IDO != 99) { @@ -573,58 +368,78 @@ int eigs(double *AR, doublecomplex *AC, int N, int Acomplex, int Asym, double* B } else // DNAUPD { - C2F(dnaupd)(&IDO, bmat, &N, which, &nev, tol, RESID, &ncv, V, &LDV, IPARAM, IPNTR, WORKD, WORKL, &LWORKL, &INFO[0], 1L, 2L); + C2F(dnaupd)(&IDO, bmat, &N, which, &nev, tol, RESID, &ncv, V, &LDV, IPARAM, IPNTR, WORKD, WORKL, &LWORKL, &INFO[0]); } - if (INFO[0] < 0) + if (INFO[0] == -1) //non critical error { - free(IPARAM); - free(IPNTR); - free(R); - free(Rprime); - free(RC); - free(RCprime); - free(tmp_WORKD); - free(tmp_WORKDC); - free(WORKD); - free(WORKL); - free(V); - free(U); - free(L); - free(E); - if (R_Rprime != NULL) + sciprint("%s: WARNING: Maximum number of iterations reached. Only %d eigenvalues converged.\n", "eigs", IPARAM[4]); + break; + } + else + { + if (INFO[0] < 0) { - free(R_Rprime); + if (R != B) + { + free(R); + } + free(IPVT); + free(AMSB); + free(WORKD); + free(WORKL); + free(V); + free(temp); + + return -4; } - return -4; } if (IDO == -1 || IDO == 1 || IDO == 2) { if (IPARAM[6] == 1) // mode = 1 { - if (matB == 0) // B = [] -> standart eigenvalue problem + if (IDO == 2) { - // OP = A*x - C2F(dgemm)("n", "n", &N, &iOne, &N, &alpha, AR, &N, WORKD + IPNTR[0] - 1, &N, &beta, WORKD + IPNTR[1] - 1, &N); + memcpy(WORKD + IPNTR[1] - 1, WORKD + IPNTR[0] - 1, N * sizeof(double)); } - else // generalized eigenvalue problem + else //IDO=1 or IDO=-1 { - // OP = inv(Rprime)*A*inv(R)*x - if (invR_A_invRprime == NULL) + if (matB == 0) // B = [] -> standard eigenvalue problem { - invR_A_invRprime = (double*)malloc(N * N * sizeof(double)); - invR_times_A_times_invRprime(invR_A_invRprime, R, AR, Rprime, N); + // OP = A*x + if (Asym) + { + C2F(dsymv) ("u", &N, &alpha, AR, &N, WORKD + IPNTR[0] - 1, &iOne, &beta, WORKD + IPNTR[1] - 1, &iOne); + } + else + { + C2F(dgemv) ("n", &N, &N, &alpha, AR, &N, WORKD + IPNTR[0] - 1, &iOne, &beta, WORKD + IPNTR[1] - 1, &iOne); + } + } + else // generalized eigenvalue problem + { + // OP = inv(Rprime)*A*inv(R)*x + memcpy(WORKD + IPNTR[1] - 1, WORKD + IPNTR[0] - 1, N * sizeof(double)); + C2F(dtrsm) ("l", "u", "n", "n", &N, &iOne, &alpha, R, &N, WORKD + IPNTR[1] - 1, &N); + memcpy(temp, WORKD + IPNTR[1] - 1, N * sizeof(double)); + if (Asym) + { + C2F(dsymv) ("u", &N, &alpha, AR, &N, temp, &iOne, &beta, WORKD + IPNTR[1] - 1, &iOne); + } + else + { + C2F(dgemv) ("n", &N, &N, &alpha, AR, &N, temp, &iOne, &beta, WORKD + IPNTR[1] - 1, &iOne); + } + C2F(dtrsm) ("l", "u", "t", "n", &N, &iOne, &alpha, R, &N, WORKD + IPNTR[1] - 1, &N); } - - C2F(dgemm)("n", "n", &N, &iOne, &N, &alpha, invR_A_invRprime, &N, WORKD + IPNTR[0] - 1, &N, &beta, WORKD + IPNTR[1] - 1, &N); } } else { if (IPARAM[6] == 3) // mode = 3 { - if (matB == 0) // B = [] -> standart eigenvalue problem + if (matB == 0) // B = [] -> standard eigenvalue problem { if (IDO == 2) { @@ -634,199 +449,120 @@ int eigs(double *AR, doublecomplex *AC, int N, int Acomplex, int Asym, double* B else { // workd[ipntr[1]-1:ipntr[1]+N-1] = inv(U)*inv(L)*inv(P)*workd[ipntr[0]-1:ipntr[0]+N-1] - - if (invU_invL_E == NULL) - { - invU_invL_E = (double*)malloc(N * N * sizeof(double)); - invU_times_invL_times_E(invU_invL_E, U, L, E, N); - } - - C2F(dgemm)("n", "n", &N, &iOne, &N, &alpha, invU_invL_E, &N, WORKD + IPNTR[0] - 1, &N, &beta, WORKD + IPNTR[1] - 1, &N); + memcpy(WORKD + IPNTR[1] - 1, WORKD + IPNTR[0] - 1, N * sizeof(double)); + C2F(dgetrs) ("n", &N, &iOne, AMSB, &N, IPVT, WORKD + IPNTR[1] - 1, &N, &INFO_INV); } } else // matB == 1 so B is not empty and bmat = 'G'-> generalized eigenvalue problem { - if (IDO == 2) + if (IDO == 2 || IDO == -1) { - if (cholB[0]) // workd[ipntr[1]-1:ipntr[1]+N-1] = R * Rprime * workd[ipntr[0]-1:ipntr[0]+N-1] + if (cholB[0]) // workd[ipntr[1]-1:ipntr[1]+N-1] = Rprime * R * workd[ipntr[0]-1:ipntr[0]+N-1] { - if (R_Rprime == NULL) - { - R_Rprime = (double*)malloc(N * N * sizeof(double)); - RtimesRprime(R_Rprime, R, Rprime, N); - } - - C2F(dgemm)("n", "n", &N, &iOne, &N, &alpha, R_Rprime, &N, WORKD + IPNTR[0] - 1, &N, &beta, WORKD + IPNTR[1] - 1, &N); + memcpy(WORKD + IPNTR[1] - 1, WORKD + IPNTR[0] - 1, N * sizeof(double)); + C2F(dtrmv) ("u", "n", "n", &N, B, &N, WORKD + IPNTR[1] - 1, &iOne); + C2F(dtrmv) ("u", "t", "n", &N, B, &N, WORKD + IPNTR[1] - 1, &iOne); } else // workd[ipntr[1]-1:ipntr[1]+N-1] = B * workd[ipntr[0]-1:ipntr[0]+N-1] { - C2F(dgemm)("n", "n", &N, &iOne, &N, &alpha, B, &N, WORKD + IPNTR[0] - 1, &N, &beta, WORKD + IPNTR[1] - 1, &N); + C2F(dgemv) ("n", &N, &N, &alpha, B, &N, WORKD + IPNTR[0] - 1, &iOne, &beta, WORKD + IPNTR[1] - 1, &iOne); } } + + if (IDO == -1) + { + // compute workd[ipntr[1]-1:ipntr[1]+N-1] = inv(U)*inv(L)*inv(P)*workd[ipntr[1]-1:ipntr[1]+N-1] + C2F(dgetrs) ("n", &N, &iOne, AMSB, &N, IPVT, WORKD + IPNTR[1] - 1, &N, &INFO_INV); + } else { - if (IDO == -1) - { - if (cholB[0]) // workd[ipntr[1]-1:ipntr[1]+N-1] = R * Rprime * workd[ipntr[0]-1:ipntr[0]+N-1] - { - if (R_Rprime == NULL) - { - R_Rprime = (double*)malloc(N * N * sizeof(double)); - RtimesRprime(R_Rprime, R, Rprime, N); - } - - C2F(dgemm)("n", "n", &N, &iOne, &N, &alpha, R_Rprime, &N, WORKD + IPNTR[0] - 1, &N, &beta, WORKD + IPNTR[1] - 1, &N); - } - else // workd[ipntr[1]-1:ipntr[1]+N-1] = B * workd[ipntr[0]-1:ipntr[0]+N-1] - { - C2F(dgemm)("n", "n", &N, &iOne, &N, &alpha, B, &N, WORKD + IPNTR[0] - 1, &N, &beta, WORKD + IPNTR[1] - 1, &N); - } - // compute workd[ipntr[1]-1:ipntr[1]+N-1] = inv(U)*inv(L)*inv(P)*workd[ipntr[1]-1:ipntr[1]+N-1] - - if (invU_invL_E == NULL) - { - invU_invL_E = (double*)malloc(N * N * sizeof(double)); - invU_times_invL_times_E(invU_invL_E, U, L, E, N); - } - - C2F(dgemm)("n", "n", &N, &iOne, &N, &alpha, invU_invL_E, &N, WORKD + IPNTR[1] - 1, &N, &beta, tmp_WORKD, &N); - memcpy(WORKD + IPNTR[1] - 1, tmp_WORKD, N * sizeof(double)); - } - else + if (IDO == 1) { - if (IDO == 1) - { - // computes workd[ipntr[1]-1:ipntr[1]+N-1] = inv(U)*inv(L)*inv(P)*workd[ipntr[2]-1:ipntr[2]+N-1] - if (invU_invL_E == NULL) - { - invU_invL_E = (double*)malloc(N * N * sizeof(double)); - invU_times_invL_times_E(invU_invL_E, U, L, E, N); - } - - C2F(dgemm)("n", "n", &N, &iOne, &N, &alpha, invU_invL_E, &N, WORKD + IPNTR[2] - 1, &N, &beta, WORKD + IPNTR[1] - 1, &N); - } + // computes workd[ipntr[1]-1:ipntr[1]+N-1] = inv(U)*inv(L)*inv(P)*workd[ipntr[2]-1:ipntr[2]+N-1] + memcpy(WORKD + IPNTR[1] - 1, WORKD + IPNTR[2] - 1, N * sizeof(double)); + C2F(dgetrs) ("n", &N, &iOne, AMSB, &N, IPVT, WORKD + IPNTR[1] - 1, &N, &INFO_INV); } } } } else { - free(IPARAM); - free(IPNTR); - free(R); - free(Rprime); - free(RC); - free(RCprime); - free(tmp_WORKD); - free(tmp_WORKDC); + if (R != B) + { + free(R); + } + free(AMSB); + free(IPVT); free(WORKD); free(WORKL); free(V); - free(U); - free(L); - free(E); + free(temp); return -5; } } } } // END WHILE - - free(L); - free(U); - free(E); - - if (R_Rprime != NULL) - { - free(R_Rprime); - } - - if (invR_A_invRprime != NULL) - { - free(invR_A_invRprime); - } - - if (invU_invL_E != NULL) - { - free(invU_invL_E); - } - - SELECT = (int*)malloc(ncv * sizeof(int)); - memset(SELECT, 0, ncv * sizeof(int)); + free(AMSB); + free(IPVT); + free(temp); + SELECT = (int *)calloc(ncv, sizeof(int)); if (Asym) // DSEUPD { - D = (double*)malloc(nev * sizeof(double)); - memset(D, 0, nev * sizeof(double)); - - Z = (double*)malloc(N * nev * sizeof(double)); - memset(Z, 0, N * nev * sizeof(double)); - - C2F(dseupd)(&RVEC, HOWMNY, SELECT, D, Z, &LDV, &SIGMAR, bmat, &N, which, &nev, tol, RESID, &ncv, V, &LDV, IPARAM, IPNTR, WORKD, WORKL, &LWORKL, &INFO_EUPD, 1L, 1L, 1L, 2L); + C2F(dseupd) (&RVEC, HOWMNY, SELECT, eigenvalue, eigenvector, &LDV, + &SIGMA.r, bmat, &N, which, &nev, tol, RESID, &ncv, V, + &LDV, IPARAM, IPNTR, WORKD, WORKL, &LWORKL, &INFO_EUPD); if (INFO_EUPD != 0) { - free(IPARAM); - free(IPNTR); - free(R); - free(Rprime); - free(RC); - free(RCprime); - free(tmp_WORKD); - free(tmp_WORKDC); + if (R != B) + { + free(R); + } free(WORKD); free(WORKL); free(V); - free(D); - free(Z); free(SELECT); return -6; } else { - for (i = 0 ; i < nev ; i++) - { - eigenvalue[i].r = D[i]; - } - if (RVEC) { - for (i = 0 ; i < N * nev ; i++) + if (matB && IPARAM[6] == 1) { - eigenvector[i].r = Z[i]; + // we need to revert back to the original problem + // since we really solved for (y,\lambda) in R^{-T}Ay=\lambda y + //with y = Rx, so that x = R^{-1}y + C2F(dtrsm) ("l", "u", "n", "n", &N, &nev, &alpha, R, &N, eigenvector, &N); } } } - - free(D); - free(Z); } else // DNEUPD { - DR = (double*)malloc((nev + 1) * sizeof(double)); - memset(DR, 0, (nev + 1) * sizeof(double)); - - DI = (double*)malloc((nev + 1) * sizeof(double)); - memset(DI, 0, (nev + 1) * sizeof(double)); + DR = (double *)calloc((nev + 1), sizeof(double)); + DI = (double *)calloc((nev + 1), sizeof(double)); + WORKEV = (double *)calloc(3 * ncv, sizeof(double)); - Z = (double*) malloc(N * (nev + 1) * sizeof(double)); - memset(Z, 0, N * (nev + 1) * sizeof(double)); + RVEC = RVEC || (IPARAM[6] == 3 && SIGMA.i != 0); - WORKEV = (double*)malloc(3 * ncv * sizeof(double)); - memset(WORKEV, 0, 3 * ncv * sizeof(double)); + if (RVEC) + { + Z = (double *)calloc(N * (nev + 1), sizeof(double)); + } - C2F(dneupd)(&RVEC, HOWMNY, SELECT, DR, DI, Z, &LDV, &SIGMAR, &SIGMAI, WORKEV, bmat, &N, which, &nev, tol, RESID, &ncv, V, &LDV, IPARAM, IPNTR, WORKD, WORKL, &LWORKL, &INFO_EUPD); + C2F(dneupd) (&RVEC, HOWMNY, SELECT, DR, DI, Z, &LDV, &SIGMA.r, + &SIGMA.i, WORKEV, bmat, &N, which, &nev, tol, RESID, + &ncv, V, &LDV, IPARAM, IPNTR, WORKD, WORKL, &LWORKL, &INFO_EUPD); if (INFO_EUPD != 0) { - free(IPARAM); - free(IPNTR); - free(R); - free(Rprime); - free(RC); - free(RCprime); - free(tmp_WORKD); - free(tmp_WORKDC); + if (R != B) + { + free(R); + } free(WORKD); free(WORKL); free(V); @@ -839,203 +575,147 @@ int eigs(double *AR, doublecomplex *AC, int N, int Acomplex, int Asym, double* B } else { - for (i = 0 ; i < nev ; i++) + if (Z && matB && IPARAM[6] == 1) { - eigenvalue[i].r = DR[i]; - eigenvalue[i].i = DI[i]; - } - - if (RVEC) - { - i = 0; - while (i <= (nev - 2)) - { - for (j = 0; j < N; j++) - { - eigenvector[i * N + j].r = Z[i * N + j]; - eigenvector[i * N + j].i = Z[(i + 1) * N + j]; - eigenvector[(i + 1)*N + j].r = Z[i * N + j]; - eigenvector[(i + 1)*N + j].i = -Z[(i + 1) * N + j]; - } - i = i + 2; - } + // we need to revert back to the original problem + // since we really solved for (y,\lambda) in R^{-T}Ay=\lambda y + //with y = Rx, so that x = R^{-1}y + C2F(dtrsm) ("l", "u", "n", "n", &N, &nev, &alpha, R, &N, Z, &N); } + //we use oldnev, because dneupd increases nev by one sometimes. + process_dneupd_data(DR, DI, Z, N, oldnev, AR, eigenvalueC, + eigenvectorC, (IPARAM[6] == 3) && (SIGMA.i != 0)); + free(DR); + free(DI); + free(Z); + free(WORKEV); } - free(DR); - free(DI); - free(Z); - free(WORKEV); } free(V); free(WORKD); free(WORKL); free(SELECT); + if (R != B) + { + free(R); + } } else // A or/and B complex { if (IPARAM[6] == 3) // mode = 3 { AMSBC = (doublecomplex*)malloc(N * N * sizeof(doublecomplex)); - memcpy(AMSBC, AC, N * N * sizeof(doublecomplex)); - if (matB == 0) // standart eigenvalue problem - { - for (i = 0 ; i < N ; i++) - { - AMSBC[i + i * N].r = AMSBC[i + i * N].r - SIGMAR; - AMSBC[i + i * N].i = AMSBC[i + i * N].i - SIGMAI; - } - } - else // generalized eigenvalue problem + + if (SIGMA.r != 0 || SIGMA.i != 0) { - if (cholB[0]) + if (matB == 0) // standard eigenvalue problem { - if (RC_RCprime == NULL) + memcpy(AMSBC, AC, N * N * sizeof(doublecomplex)); + for (i = 0 ; i < N ; i++) { - RC_RCprime = (doublecomplex*)malloc(N * N * sizeof(doublecomplex)); - RCtimesRCprime(RC_RCprime, RC, RCprime, N); - } - - for (i = 0 ; i < N * N ; i++) - { - AMSBC[i].r = AMSBC[i].r - (SIGMAR * RC_RCprime[i].r + SIGMAI * RC_RCprime[i].i); - AMSBC[i].i = AMSBC[i].i - (SIGMAR * RC_RCprime[i].i + SIGMAI * RC_RCprime[i].r); + AMSBC[i + i * N].r -= SIGMA.r; + AMSBC[i + i * N].i -= SIGMA.i; } } - else + else // generalized eigenvalue problem { - for (i = 0 ; i < N * N ; i++) + if (cholB[0]) { - AMSBC[i].r = AMSBC[i].r - (SIGMA[0].r * BC[i].r); - AMSBC[i].i = AMSBC[i].i - (SIGMA[0].i * BC[i].i); + memcpy(AMSBC, BC, N * N * sizeof(doublecomplex)); + C2F(ztrmm)("l", "u", "c", "n", &N, &N, &mSIGMA, BC, &N, AMSBC, &N); + C2F(zaxpy)(&N2, &alphac, AC, &iOne, AMSBC, &iOne); + } + else + { + memcpy(AMSBC, AC, N * N * sizeof(doublecomplex)); + C2F(zaxpy) (&N2, &mSIGMA, BC, &iOne, AMSBC, &iOne); } } } - - // LU decomposition - IPVT = (int*) malloc(N * sizeof(int)); - memset(IPVT, 0, N * sizeof(int)); - C2F(zgetrf)(&N, &N, AMSBC, &N, IPVT, &INFO_LU); - - // Computes the lower triangular matrix L - LC = (doublecomplex*)malloc(N * N * sizeof(doublecomplex)); - memset(LC, 0, N * N * sizeof(doublecomplex)); - for (i = 0 ; i < N ; i++) - { - for (j = 0 ; j < i ; j++) - { - LC[i + j * N].r = AMSBC[i + j * N].r; - LC[i + j * N].i = AMSBC[i + j * N].i; - } - LC[i + i * N].r = 1; - LC[i + i * N].i = 0; - } - - // Computes the upper triangular matrix U - - UC = (doublecomplex*)malloc(N * N * sizeof(doublecomplex)); - memset(UC, 0, N * N * sizeof(doublecomplex)); - for (j = 0 ; j < N ; j++) - { - for (i = 0 ; i <= j ; i++) - { - UC[i + j * N].r = AMSBC[i + j * N].r; - UC[i + j * N].i = AMSBC[i + j * N].i; - } - } - - // Computes the permutation matrix E - EC = (doublecomplex*)malloc(N * N * sizeof(doublecomplex)); - E = (double*)malloc(N * N * sizeof(double)); - memset(E, 0, N * N * sizeof(double)); - - for (i = 0 ; i < N ; i++) + else { - E[i * N + i] = 1; + memcpy(AMSBC, AC, N * N * sizeof(doublecomplex)); } - C2F(dlaswp)(&N, E, &N, &k1, &N, IPVT, &k1); - - memset(EC, 0, N * N * sizeof(doublecomplex)); - for (i = 0 ; i < N * N ; i++) + // LU decomposition + IPVT = (int*) calloc(N, sizeof(int)); + C2F(zgetrf) (&N, &N, AMSBC, &N, IPVT, &INFO_LU); + if (INFO_LU > 0) { - EC[i].r = E[i]; + free(IPVT); + free(AMSBC); + return(-7); } - - free(AMSBC); - free(IPVT); } - LWORKL = 3 * ncv * ncv + 5 * ncv; - VC = (doublecomplex*)malloc(N * ncv * sizeof(doublecomplex)); - memset(VC, 0, N * ncv * sizeof(doublecomplex)); - - WORKLC = (doublecomplex*)malloc(LWORKL * sizeof(doublecomplex)); - memset(WORKLC, 0, LWORKL * sizeof(doublecomplex)); - - WORKDC = (doublecomplex*)malloc(3 * N * sizeof(doublecomplex)); - memset(WORKDC, 0, 3 * N * sizeof(doublecomplex)); - - RWORK = (double*)malloc(ncv * sizeof(double)); - memset(RWORK, 0, ncv * sizeof(double)); + VC = (doublecomplex*) calloc(N * ncv, sizeof(doublecomplex)); + WORKLC = (doublecomplex*) calloc(LWORKL, sizeof(doublecomplex)); + WORKDC = (doublecomplex*) calloc(3 * N, sizeof(doublecomplex)); + RWORK = (double*) calloc(ncv, sizeof(double)); + if (IPARAM[6] == 1 && matB) + { + tempC = (doublecomplex*) malloc(N * sizeof(doublecomplex)); + } while (IDO != 99) { C2F(znaupd)(&IDO, bmat, &N, which, &nev, tol, RESIDC, &ncv, VC, &LDV, IPARAM, IPNTR, WORKDC, WORKLC, &LWORKL, RWORK, &INFO[0]); - if (INFO[0] < 0) + if (INFO[0] == -1) //non critical error { - free(IPARAM); - free(IPNTR); - free(R); - free(Rprime); - free(RC); - free(RCprime); - free(tmp_WORKD); - free(tmp_WORKDC); - free(LC); - free(UC); - free(EC); - free(E); - free(WORKDC); - free(WORKLC); - free(VC); - free(RWORK); - if (RC_RCprime != NULL) + sciprint("%s: WARNING: Maximum number of iterations reached. Only %d eigenvalues converged.\n", "eigs", IPARAM[4]); + break; + } + else + { + if (INFO[0] < 0) { - free(RC_RCprime); + if (RC != BC) + { + free(RC); + } + free(WORKDC); + free(WORKLC); + free(VC); + free(RWORK); + return -4; } - return -4; } if (IDO == -1 || IDO == 1 || IDO == 2) { if (IPARAM[6] == 1) // mode = 1 { - if (matB == 0) // B = I + if (IDO == 2) { - // OP = A*x - C2F(zgemm)("n", "n", &N, &iOne, &N, &alphac, AC, &N, WORKDC + IPNTR[0] - 1, &N, &betac, WORKDC + IPNTR[1] - 1, &N); + memcpy(WORKDC + IPNTR[1] - 1, WORKDC + IPNTR[0] - 1, N * sizeof(doublecomplex)); } else { - // OP = inv(R')*A*inv(R)*x - if (invRC_AC_invRCprime == NULL) + if (matB == 0) // B = I { - invRC_AC_invRCprime = (doublecomplex*)malloc(N * N * sizeof(doublecomplex)); - invRC_times_AC_times_invRCprime(invRC_AC_invRCprime, RC, AC, RCprime, N); + // OP = A*x + C2F(zgemv) ("n", &N, &N, &alphac, AC, &N, WORKDC + IPNTR[0] - 1, &iOne, &betac, WORKDC + IPNTR[1] - 1, &iOne); + } + else + { + // OP = inv(RC')*A*inv(RC)*x + memcpy(WORKDC + IPNTR[1] - 1, WORKDC + IPNTR[0] - 1, N * sizeof(doublecomplex)); + C2F(ztrsm) ("l", "u", "n", "n", &N, &iOne, &alphac, RC, &N, WORKDC + IPNTR[1] - 1, &N); + memcpy(tempC, WORKDC + IPNTR[1] - 1, N * sizeof(doublecomplex)); + C2F(zgemv) ("n", &N, &N, &alphac, AC, &N, tempC, &iOne, &betac, WORKDC + IPNTR[1] - 1, &iOne); + C2F(ztrsm) ("l", "u", "c", "n", &N, &iOne, &alphac, RC, &N, WORKDC + IPNTR[1] - 1, &N); } - - C2F(zgemm)("n", "n", &N, &iOne, &N, &alphac, invRC_AC_invRCprime, &N, WORKDC + IPNTR[0] - 1, &N, &betac, WORKDC + IPNTR[1] - 1, &N); } } else { - if (IPARAM[6] == 3) // si mode = 3 + if (IPARAM[6] == 3) // if mode = 3 { - if (matB == 0) // B = [] -> matB is empty -> standart eigenvalue problem + if (matB == 0) // B = [] -> matB is empty -> standard eigenvalue problem { if (IDO == 2) { @@ -1045,166 +725,80 @@ int eigs(double *AR, doublecomplex *AC, int N, int Acomplex, int Asym, double* B else { // workd[ipntr[1]-1:ipntr[1]+N-1] = inv(U)*inv(L)*inv(P)*workd[ipntr[0]-1:ipntr[0]+N-1] - if (invUC_invLC_EC == NULL) - { - invUC_invLC_EC = (doublecomplex*)malloc(N * N * sizeof(doublecomplex)); - invUC_times_invLC_times_EC(invUC_invLC_EC, UC, LC, EC, N); - } - C2F(zgemm)("n", "n", &N, &iOne, &N, &alphac, invUC_invLC_EC, &N, WORKDC + IPNTR[0] - 1, &N, &betac, WORKDC + IPNTR[1] - 1, &N); + memcpy(WORKDC + IPNTR[1] - 1, WORKDC + IPNTR[0] - 1, N * sizeof(doublecomplex)); + C2F(zgetrs) ("n", &N, &iOne, AMSBC, &N, IPVT, WORKDC + IPNTR[1] - 1, &N, &INFO_INV); } } else // matB == 1 so B is not empty and bmat = 'G'-> generalized eigenvalue problem { - if (IDO == 2) + if (IDO == 2 || IDO == -1) { - if (cholB[0]) // workd[ipntr[1]-1:ipntr[1]+N-1] = RC * RCprime * workd[ipntr[0]-1:ipntr[0]+N-1] + if (cholB[0]) // workd[ipntr[1]-1:ipntr[1]+N-1] = RCprime * RC * workd[ipntr[0]-1:ipntr[0]+N-1] { - if (RC_RCprime == NULL) - { - RC_RCprime = (doublecomplex*)malloc(N * N * sizeof(doublecomplex)); - RCtimesRCprime(RC_RCprime, RC, RCprime, N); - } - - C2F(zgemm)("n", "n", &N, &iOne, &N, &alphac, RC_RCprime, &N, WORKDC + IPNTR[0] - 1, &N, &betac, WORKDC + IPNTR[1] - 1, &N); + memcpy(WORKDC + IPNTR[1] - 1, WORKDC + IPNTR[0] - 1, N * sizeof(doublecomplex)); + C2F(ztrmv) ("u", "n", "n", &N, BC, &N, WORKDC + IPNTR[1] - 1, &iOne); + C2F(ztrmv) ("u", "c", "n", &N, BC, &N, WORKDC + IPNTR[1] - 1, &iOne); } else // workd[ipntr[1]-1:ipntr[1]+N-1] = B *workd[ipntr[0]-1:ipntr[0]+N-1] { - C2F(zgemm)("n", "n", &N, &iOne, &N, &alphac, BC, &N, WORKDC + IPNTR[0] - 1, &N, &betac, WORKDC + IPNTR[1] - 1, &N); + C2F(zgemv) ("n", &N, &N, &alphac, BC, &N, WORKDC + IPNTR[0] - 1, &iOne, &betac, WORKDC + IPNTR[1] - 1, &iOne); } } + if (IDO == -1) + { + // compute workd[ipntr[1]-1:ipntr[1]+N-1] = inv(U)*inv(L)*inv(P)*workd[ipntr[1]-1:ipntr[1]+N-1] + C2F(zgetrs) ("n", &N, &iOne, AMSBC, &N, IPVT, WORKDC + IPNTR[1] - 1, &N, &INFO_INV); + } else { - if (IDO == -1) - { - if (cholB[0]) // workd[ipntr[1]-1:ipntr[1]+N-1] = RC*RCprime*workd[ipntr[0]-1:ipntr[0]+N-1] - { - if (RC_RCprime == NULL) - { - RC_RCprime = (doublecomplex*)malloc(N * N * sizeof(doublecomplex)); - RCtimesRCprime(RC_RCprime, RC, RCprime, N); - } - - C2F(zgemm)("n", "n", &N, &iOne, &N, &alphac, RC_RCprime, &N, WORKDC + IPNTR[0] - 1, &N, &betac, WORKDC + IPNTR[1] - 1, &N); - } - else // workd[ipntr[1]-1:ipntr[1]+N-1] = B * workd[ipntr[0]-1:ipntr[0]+N-1] - { - C2F(zgemm)("n", "n", &N, &iOne, &N, &alphac, BC, &N, WORKDC + IPNTR[0] - 1, &N, &betac, WORKDC + IPNTR[1] - 1, &N); - } - - // compute workd[ipntr[1]-1:ipntr[1]+N-1] = inv(U)*inv(L)*inv(P)*workd[ipntr[1]-1:ipntr[1]+N-1] - if (invUC_invLC_EC == NULL) - { - invUC_invLC_EC = (doublecomplex*)malloc(N * N * sizeof(doublecomplex)); - invUC_times_invLC_times_EC(invUC_invLC_EC, UC, LC, EC, N); - } - - C2F(zgemm)("n", "n", &N, &iOne, &N, &alphac, invUC_invLC_EC, &N, WORKDC + IPNTR[1] - 1, &N, &betac, tmp_WORKDC, &N); - memcpy(WORKDC + IPNTR[1] - 1, tmp_WORKDC, N * sizeof(doublecomplex*)); - } - else + if (IDO == 1) { - if (IDO == 1) - { - // compute workd[ipntr[1]-1:ipntr[1]+N-1] = inv(U)*inv(L)*inv(P)*workd[ipntr[2]-1:ipntr[2]+N-1] - if (invUC_invLC_EC == NULL) - { - invUC_invLC_EC = (doublecomplex*)malloc(N * N * sizeof(doublecomplex)); - invUC_times_invLC_times_EC(invUC_invLC_EC, UC, LC, EC, N); - } - - C2F(zgemm)("n", "n", &N, &iOne, &N, &alphac, invUC_invLC_EC, &N, WORKDC + IPNTR[2] - 1, &N, &betac, WORKDC + IPNTR[1] - 1, &N); - } + /* compute workd[ipntr[1]-1:ipntr[1]+N-1] = inv(U)*inv(L)*inv(P)*workd[ipntr[2]-1:ipntr[2]+N-1] */ + memcpy(WORKDC + IPNTR[1] - 1, WORKDC + IPNTR[2] - 1, N * sizeof(doublecomplex)); + C2F(zgetrs) ("n", &N, &iOne, AMSBC, &N, IPVT, WORKDC + IPNTR[1] - 1, &N, &INFO_INV); } } - } + } //END mode3 } else { - free(IPARAM); - free(IPNTR); - free(R); - free(Rprime); - free(RC); - free(RCprime); - free(tmp_WORKD); - free(tmp_WORKDC); - free(LC); - free(UC); - free(EC); - free(E); + if (RC != BC) + { + free(RC); + } free(WORKDC); free(WORKLC); free(VC); free(RWORK); - if (RC_RCprime != NULL) - { - free(RC_RCprime); - } - - if (invRC_AC_invRCprime != NULL) - { - free(invRC_AC_invRCprime); - } + free(tempC); - if (invUC_invLC_EC != NULL) - { - free(invUC_invLC_EC); - } return -5; } } } } // END WHILE - free(LC); - free(UC); - free(EC); - free(E); - - if (RC_RCprime != NULL) - { - free(RC_RCprime); - } - - if (invRC_AC_invRCprime != NULL) - { - free(invRC_AC_invRCprime); - } + free(tempC); + free(IPVT); + free(AMSBC); - if (invUC_invLC_EC != NULL) - { - free(invUC_invLC_EC); - } - - SELECT = (int*)malloc(ncv * sizeof(int)); - memset(SELECT, 0, ncv * sizeof(int)); - - DC = (doublecomplex*)malloc((nev + 1) * sizeof(doublecomplex)); - memset(DC, 0, (nev + 1) * sizeof(doublecomplex)); - - ZC = (doublecomplex*)malloc(N * nev * sizeof(doublecomplex)); - memset(ZC, 0, N * nev * sizeof(doublecomplex)); + SELECT = (int *)calloc(ncv, sizeof(int)); + WORKEVC = (doublecomplex *) calloc(2 * ncv, sizeof(doublecomplex)); - WORKEVC = (doublecomplex*)malloc(2 * ncv * sizeof(doublecomplex)); - memset(WORKEVC, 0, 2 * ncv * sizeof(doublecomplex)); + C2F(zneupd) (&RVEC, HOWMNY, SELECT, eigenvalueC, eigenvectorC, &LDV, &SIGMA, WORKEVC, bmat, &N, + which, &nev, tol, RESIDC, &ncv, VC, &LDV, IPARAM, IPNTR, WORKDC, + WORKLC, &LWORKL, RWORK, &INFO_EUPD); - C2F(zneupd)(&RVEC, HOWMNY, SELECT, DC, ZC, &LDV, SIGMA, WORKEVC, bmat, &N, which, &nev, tol, RESIDC, &ncv, VC, &LDV, IPARAM, IPNTR, WORKDC, WORKLC, &LWORKL, RWORK, &INFO_EUPD); if (INFO_EUPD != 0) { - free(IPARAM); - free(IPNTR); - free(R); - free(Rprime); - free(RC); - free(RCprime); - free(tmp_WORKD); - free(tmp_WORKDC); + if (RC != BC) + { + free(RC); + } free(WORKDC); free(WORKLC); free(VC); free(SELECT); - free(DC); - free(ZC); free(WORKEVC); free(RWORK); @@ -1212,42 +806,27 @@ int eigs(double *AR, doublecomplex *AC, int N, int Acomplex, int Asym, double* B } else { - if (!RVEC) + if (RVEC) { - for (i = 0; i < nev; i++) + if (matB && IPARAM[6] == 1) { - eigenvalue[i].r = DC[i].r; - eigenvalue[i].i = DC[i].i; + C2F(ztrsm) ("l", "u", "n", "n", &N, &nev, &alphac, RC, &N, eigenvectorC, &N); } } - else // return eigenvalues and eigenvectors - { - memcpy(eigenvalue, DC, nev * sizeof(doublecomplex)); - memcpy(eigenvector, ZC, N * nev * sizeof(doublecomplex)); - } } free(SELECT); - free(DC); - free(ZC); free(WORKEVC); free(VC); - free(tmp_WORKDC); - free(tmp_WORKD); free(WORKDC); free(WORKLC); free(RWORK); + if (RC != BC) + { + free(RC); + } } - free(IPARAM); - free(IPNTR); - - free(R); - free(Rprime); - free(RC); - free(RCprime); - return 0; } - diff --git a/scilab/modules/arnoldi/src/c/eigs_dependencies.c b/scilab/modules/arnoldi/src/c/eigs_dependencies.c index 3c97f76..419def3 100644 --- a/scilab/modules/arnoldi/src/c/eigs_dependencies.c +++ b/scilab/modules/arnoldi/src/c/eigs_dependencies.c @@ -1,224 +1,121 @@ /* * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab * Copyright (C) 2011 - Scilab Enterprises - Adeline CARNIS - * + * * This file must be used under the terms of the CeCILL. * This source file is licensed as described in the file COPYING, which * you should have received as part of this distribution. The terms - * are also available at + * are also available at * http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt * */ #include "eigs_dependencies.h" -// dgemm performs one of the matrix-matrix operations -extern int C2F(dgemm)(char* transa, char* transb, int* m, int* n, int* k, double* alpha, double* A, int* lda, double* B, int* ldb, double* beta, double* C, int* ldc); - -// dgetrf computes an LU factorization of a general M by N matrix A (double) using partial pivoting with row interchanges -extern int C2F(dgetrf)(int* m, int* n, double* A, int* lda, int* ipiv, int* info); - -// dgetri computes the inverse of a matrix using the LU factorization computed by dgetrf -extern int C2F(dgetri)(int* n, double* A, int* lda, int* ipiv, double* work, int* lworkl, int* info); - -// zgemm performs one of the matrix-matrix operations -extern int C2F(zgemm)(char* transa, char* transb, int* m, int* n, int* k, doublecomplex* alpha, doublecomplex* A, int* lda, - doublecomplex* B, int* ldb, doublecomplex* beta, doublecomplex* C, int* ldc); - -// zgetrd computes an LU factorization of a general M by N matrix A (complex*16) using partial pivoting with row interchanges -extern int C2F(zgetrf)(int* m, int* n, doublecomplex* A, int* lda, int* ipiv, int* info); - -// zgetri computes the inverse of a matrix using the LU factorization computed by zgetrf -extern int C2F(zgetri)(int* n, doublecomplex* A, int* lda, int* ipiv, doublecomplex* work, int* lworkl, int* info); - -static double alpha = 1.; -static double beta = 0.; - -static doublecomplex alphac = {1.,0.}; -static doublecomplex betac = {0.,0.}; - -void RtimesRprime(double* result, double* R, double* Rprime, int N) -{ - C2F(dgemm)("n", "n", &N, &N, &N, &alpha, R, &N, Rprime, &N, &beta, result, &N); -} - -void invR_times_A_times_invRprime(double* result, double* R, double* A, double* Rprime, int N) -{ - double* invRxA = NULL; - int* IPVT = NULL; - double* work = NULL; - double* invR = NULL; - double* invRprime = NULL; - int INFO_LU = 0; - - invRxA = (double*)malloc(N * N * sizeof(double)); - invR = (double*)malloc(N * N * sizeof(double)); - invRprime = (double*)malloc(N * N * sizeof(double)); - work = (double*)malloc(N * N * sizeof(double)); - - - IPVT = (int*) malloc(N * sizeof(int)); - memset(IPVT, 0, N * sizeof(int)); - - // inv(R) - memcpy(invR, R, N * N * sizeof(double)); // copy R to invR - C2F(dgetrf)(&N, &N, invR ,&N, IPVT, &INFO_LU); // LU decomposition - - memset(work, 0, N * N * sizeof(double)); - C2F(dgetri)(&N, invR, &N, IPVT, work, &N, &INFO_LU); // Obtain inverse matrix R - - // inv(Rprime) - memset(IPVT, 0, N * sizeof(int)); - memcpy(invRprime, Rprime, N * N * sizeof(double)); - C2F(dgetrf)(&N, &N, invRprime, &N, IPVT, &INFO_LU); // LU decomposition - - memset(work, 0, N * N * sizeof(double)); - C2F(dgetri)(&N, invRprime, &N, IPVT, work, &N,&INFO_LU); // Obtain inverse matrix Rprime - - C2F(dgemm)("n", "n", &N, &N, &N, &alpha, invR, &N, A, &N, &beta, invRxA, &N); - - C2F(dgemm)("n", "n", &N, &N, &N, &alpha, invRxA, &N, invRprime, &N, &beta, result, &N); - - free(invRxA); - free(IPVT); - free(work); - free(invR); - free(invRprime); -} - -void invU_times_invL_times_E(double* result, double* U, double* L, double* E, int N) -{ - double* invUxinvL = NULL; - int* IPVT = NULL; - double* work = NULL; - double* invU = NULL; - double* invL = NULL; - int INFO_LU = 0; - - invUxinvL = (double*)malloc(N * N * sizeof(double)); - invU = (double*)malloc(N * N * sizeof(double)); - invL = (double*)malloc(N * N* sizeof(double)); - work = (double*)malloc(N * N * sizeof(double)); - - IPVT = (int*) malloc(N * sizeof(int)); - memset(IPVT, 0, N * sizeof(int)); - - // inv L -> L obtained with LU decomposition - memcpy(invL, L, N * N * sizeof(double)); - C2F(dgetrf)(&N, &N, invL, &N, IPVT, &INFO_LU); // LU decomposition - - memset(work, 0, N * N * sizeof(double)); - C2F(dgetri)(&N, invL, &N, IPVT, work, &N, &INFO_LU); // inv(L) - - // inv U -> U obtained with LU decomposition - memcpy(invU, U, N * N * sizeof(double)); - memset(IPVT, 0, N*sizeof(int)); - C2F(dgetrf)(&N, &N, invU, &N, IPVT, &INFO_LU); // LU decomposition +extern int C2F(dgemv) (char* trans, int* m, int* n, double* alpha, double* A, + int* lda, double* x, int* incx, double* beta, double* y, int* incy); +extern double C2F(ddot) (int *n, double* x, int* incx, double* y, int* incy); - memset(work, 0, N * N * sizeof(double)); - C2F(dgetri)(&N, invU, &N, IPVT, work, &N, &INFO_LU); // inv(U) - - C2F(dgemm)("n", "n", &N, &N, &N, &alpha, invU, &N, invL, &N, &beta, invUxinvL, &N); - - C2F(dgemm)("n", "n", &N, &N, &N, &alpha, invUxinvL, &N, E, &N, &beta, result, &N); - - free(invUxinvL); - free(IPVT); - free(work); - free(invU); - free(invL); -} - - -void RCtimesRCprime(doublecomplex* result, doublecomplex* R, doublecomplex* Rprime, int N) +void process_dneupd_data(double* DR, double* DI, double* Z, int N, int nev, double* AR, + doublecomplex* eigenvalue, doublecomplex* eigenvector, + int sigma_imaginary) { - C2F(zgemm)("n", "n", &N, &N, &N, &alphac, R, &N, Rprime, &N, &betac, result, &N); -} - -void invRC_times_AC_times_invRCprime(doublecomplex* result, doublecomplex* R, doublecomplex* A, doublecomplex* Rprime, int N) -{ - doublecomplex* invRxA = NULL; - int* IPVT = NULL; - doublecomplex* work = NULL; - doublecomplex* invR = NULL; - doublecomplex* invRprime = NULL; - int INFO_LU = 0; - - invRxA = (doublecomplex*)malloc(N * N * sizeof(doublecomplex)); - invR = (doublecomplex*)malloc(N * N * sizeof(doublecomplex)); - invRprime = (doublecomplex*)malloc(N * N * sizeof(doublecomplex)); - work = (doublecomplex*)malloc(N * N * sizeof(doublecomplex)); - - - IPVT = (int*) malloc(N * sizeof(int)); - memset(IPVT, 0, N * sizeof(int)); - - // inv(R) - memcpy(invR, R, N * N * sizeof(doublecomplex)); // copy R to invR - C2F(zgetrf)(&N, &N, invR ,&N, IPVT, &INFO_LU); // LU decomposition - - memset(work, 0, N * N * sizeof(doublecomplex)); - C2F(zgetri)(&N, invR, &N, IPVT, work, &N, &INFO_LU); // Obtain inverse matrix R - - // inv(Rprime) - memset(IPVT, 0, N * sizeof(int)); - memcpy(invRprime, Rprime, N * N * sizeof(doublecomplex)); - C2F(zgetrf)(&N, &N, invRprime, &N, IPVT, &INFO_LU); // LU decomposition - - memset(work, 0, N * N * sizeof(doublecomplex)); - C2F(zgetri)(&N, invRprime, &N, IPVT, work, &N,&INFO_LU); // Obtain inverse matrix Rprime - - C2F(zgemm)("n", "n", &N, &N, &N, &alphac, invR, &N, A, &N, &betac, invRxA, &N); - - C2F(zgemm)("n", "n", &N, &N, &N, &alphac, invRxA, &N, invRprime, &N, &betac, result, &N); - - free(invRxA); - free(IPVT); - free(work); - free(invR); - free(invRprime); + /* if sigma_imaginary there is an extra step to compute the eigenvalues + as explained in the dneupd user guide */ + + double* temp1 = NULL; + double* temp2 = NULL; + + int i = 0; + int j = 0; + + double alpha = 1; + double beta = 0; + int iOne = 1; + double real_part; + double imag_part; + + if ( sigma_imaginary ) + { + temp1 = (double*) malloc(N * sizeof(double)); + temp2 = (double*) malloc(N * sizeof(double)); + + while (i < nev) + { + if (DI[i] == 0) + { + C2F(dgemv) ("n", &N, &N, &alpha, AR, &N, Z + N * i, &iOne, &beta, temp1, &iOne); + eigenvalue[i] = (doublecomplex) + { + C2F(ddot) (&N, Z + N * i, &iOne, temp1, &iOne), 0 + }; + i = i + 1; + } + else + { + C2F(dgemv) ("n", &N, &N, &alpha, AR, &N, Z + N * i, &iOne, &beta, temp1, &iOne); + C2F(dgemv) ("n", &N, &N, &alpha, AR, &N, Z + N * (i + 1), &iOne, &beta, temp2, &iOne); + real_part = C2F(ddot) (&N, Z + N * i, &iOne, temp1, &iOne) + \ + C2F(ddot) (&N, Z + N * (i + 1), &iOne, temp2, &iOne); + imag_part = C2F(ddot) (&N, Z + N * i, &iOne, temp2, &iOne) - \ + C2F(ddot) (&N, Z + N * (i + 1), &iOne, temp1, &iOne); + eigenvalue[i] = (doublecomplex) + { + real_part, imag_part + }; + eigenvalue[i + 1] = (doublecomplex) + { + real_part, -imag_part + }; + i = i + 2; + } + } + free(temp1); + free(temp2); + } + else + { + for (i = 0; i < nev + 1; i++) + { + eigenvalue[i] = (doublecomplex) + { + DR[i], DI[i] + }; + } + } + + if (eigenvector) + { + i = 0; + + while (i < nev) + { + if (DI[i] != 0) + { + for (j = 0; j < N; j++) + { + eigenvector[i * N + j] = (doublecomplex) + { + Z[i * N + j], Z[(i + 1) * N + j] + }; + eigenvector[(i + 1) * N + j] = (doublecomplex) + { + Z[i * N + j], -Z[(i + 1) * N + j] + }; + } + i = i + 2; + } + else + { + for (j = 0; j < N; j++) + { + eigenvector[i * N + j] = (doublecomplex) + { + Z[i * N + j], 0 + }; + } + i = i + 1; + } + } + } } - - -void invUC_times_invLC_times_EC(doublecomplex* result, doublecomplex* U, doublecomplex* L, doublecomplex* E, int N) -{ - doublecomplex* invUxinvL = NULL; - int* IPVT = NULL; - doublecomplex* work = NULL; - doublecomplex* invU = NULL; - doublecomplex* invL = NULL; - int INFO_LU = 0; - - invUxinvL = (doublecomplex*)malloc(N * N * sizeof(doublecomplex)); - invU = (doublecomplex*)malloc(N * N * sizeof(doublecomplex)); - invL = (doublecomplex*)malloc(N * N* sizeof(doublecomplex)); - work = (doublecomplex*)malloc(N * N * sizeof(doublecomplex)); - - IPVT = (int*) malloc(N * sizeof(int)); - memset(IPVT, 0, N * sizeof(int)); - - // inv L -> L obtained with LU decomposition - memcpy(invL, L, N * N * sizeof(doublecomplex)); - C2F(zgetrf)(&N, &N, invL, &N, IPVT, &INFO_LU); // LU decomposition - - memset(work, 0, N * N * sizeof(doublecomplex)); - C2F(zgetri)(&N, invL, &N, IPVT, work, &N, &INFO_LU); // inv(L) - - // inv U -> U obtained with LU decomposition - memcpy(invU, U, N * N * sizeof(doublecomplex)); - memset(IPVT, 0, N*sizeof(int)); - C2F(zgetrf)(&N, &N, invU, &N, IPVT, &INFO_LU); // LU decomposition - - memset(work, 0, N * N * sizeof(doublecomplex)); - C2F(zgetri)(&N, invU, &N, IPVT, work, &N, &INFO_LU); // inv(U) - - C2F(zgemm)("n", "n", &N, &N, &N, &alphac, invU, &N, invL, &N, &betac, invUxinvL, &N); - - C2F(zgemm)("n", "n", &N, &N, &N, &alphac, invUxinvL, &N, E, &N, &betac, result, &N); - - free(invUxinvL); - free(IPVT); - free(work); - free(invU); - free(invL); -} \ No newline at end of file diff --git a/scilab/modules/arnoldi/tests/nonreg_tests/bug_12138.dia.ref b/scilab/modules/arnoldi/tests/nonreg_tests/bug_12138.dia.ref new file mode 100644 index 0000000..a40050a --- /dev/null +++ b/scilab/modules/arnoldi/tests/nonreg_tests/bug_12138.dia.ref @@ -0,0 +1,81 @@ +// ============================================================================= +// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab +// Copyright (C) 2012 - Scilab Enterprises - Sylvestre Ledru +// +// This file is distributed under the same license as the Scilab package. +// ============================================================================= +// <-- CLI SHELL MODE --> +// <-- Non-regression test for bug 12138 --> +// +// <-- Bugzilla URL --> +// http://bugzilla.scilab.org/show_bug.cgi?id=12138 +// +// <-- Short Description --> +// eigs(A,B) returns incorrect eigenvectors for dense matrices +// ============================================================================= +//non symmetric case +A = rand(10,10); +[d v] = eigs(A); +assert_checkalmostequal(norm(A*v-v*d),0,[], 1D-8); +[d v] = eigs(A,[],8,'SM'); +assert_checkalmostequal(norm(A*v-v*d),0,[], 1D-8); +[d v] = eigs(A,[],8,1); +assert_checkalmostequal(norm(A*v-v*d),0,[], 1D-8); +[d v] = eigs(A,[],8,%i); +assert_checkalmostequal(norm(A*v-v*d),0,[], 1D-8); +//symmetric case +A=rand(10,10); +A = A*A'; +[d v] = eigs(A); +assert_checkalmostequal(norm(A*v-v*d),0,[], 1D-8); +//general eigenvalue problem +B = rand(10,10); +B = B*B'; +A = rand(10,10); +[d v] = eigs(A,B); +assert_checkalmostequal(norm(A*v-B*v*d),0,[],1D-8); +[d v] = eigs(A,B,8,'SM'); +assert_checkalmostequal(norm(A*v-B*v*d),0,[],1D-8); +[d v] = eigs(A,B,8, 1); +assert_checkalmostequal(norm(A*v-B*v*d),0,[],1D-8); +[d v] = eigs(A,B, 8, %i); +assert_checkalmostequal(norm(A*v-B*v*d),0,[],1D-8); +opts.cholB = %t; +R = chol(B); +[d v] = eigs(A, R, 8, 'LM', opts); +assert_checkalmostequal(norm(A*v-B*v*d),0,[],1D-8); +[d v] = eigs(A, R, 8,'SM', opts); +assert_checkalmostequal(norm(A*v-B*v*d),0,[],1D-8); +[d v] = eigs(A, R, 8, 1, opts); +assert_checkalmostequal(norm(A*v-B*v*d),0,[],1D-8); +[d v] = eigs(A, R, 8, %i, opts); +assert_checkalmostequal(norm(A*v-B*v*d),0,[],1D-8); +A=A*A'; +[d v] = eigs(A,B); +assert_checkalmostequal(norm(A*v-B*v*d),0,[],1D-8); +[d v] = eigs(A,B,8,'SM'); +assert_checkalmostequal(norm(A*v-B*v*d),0,[],1D-8); +[d v] = eigs(A,B,8,1); +assert_checkalmostequal(norm(A*v-B*v*d),0,[],1D-8); +[d v] = eigs(A,B,8,%i); +assert_checkalmostequal(norm(A*v-B*v*d),0,[],1D-8); +//complex case +A1 = rand(10,10); +A2 = rand(10,10); +B1 = rand(10,10); +B2 = rand(10,10); +C1 = A1+%i*A2; +[d v] = eigs(C1); +assert_checkalmostequal(norm(C1*v-v*d),0,[],1D-8); +[d v] = eigs(C1,[], 8,'SM'); +assert_checkalmostequal(norm(C1*v-v*d),0,[],1D-8); +[d v] = eigs(C1, [],8, 1+%i); +assert_checkalmostequal(norm(C1*v-v*d),0,[],1D-8); +C2 = B1+%i*B2; +C2=C2*C2'; +[d v] = eigs(C1, C2); +assert_checkalmostequal(norm(C1*v-C2*v*d),0,[],1D-8); +[d v] = eigs(C1, C2, 8, 'SM'); +assert_checkalmostequal(norm(C1*v-C2*v*d),0,[],1D-8); +[d v] = eigs(C1, C2, 8, %i); +assert_checkalmostequal(norm(C1*v-C2*v*d),0,[],1D-8); diff --git a/scilab/modules/arnoldi/tests/nonreg_tests/bug_12138.tst b/scilab/modules/arnoldi/tests/nonreg_tests/bug_12138.tst new file mode 100644 index 0000000..c83c251 --- /dev/null +++ b/scilab/modules/arnoldi/tests/nonreg_tests/bug_12138.tst @@ -0,0 +1,90 @@ +// ============================================================================= +// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab +// Copyright (C) 2012 - Scilab Enterprises - Sylvestre Ledru +// +// This file is distributed under the same license as the Scilab package. +// ============================================================================= + +// <-- CLI SHELL MODE --> + +// <-- Non-regression test for bug 12138 --> +// +// <-- Bugzilla URL --> +// http://bugzilla.scilab.org/show_bug.cgi?id=12138 +// +// <-- Short Description --> +// eigs(A,B) returns incorrect eigenvectors for dense matrices +// ============================================================================= + +//non symmetric case +A = rand(10,10); +[d v] = eigs(A); +assert_checkalmostequal(norm(A*v-v*d),0,[], 1D-8); +[d v] = eigs(A,[],8,'SM'); +assert_checkalmostequal(norm(A*v-v*d),0,[], 1D-8); +[d v] = eigs(A,[],8,1); +assert_checkalmostequal(norm(A*v-v*d),0,[], 1D-8); +[d v] = eigs(A,[],8,%i); +assert_checkalmostequal(norm(A*v-v*d),0,[], 1D-8); + +//symmetric case +A=rand(10,10); +A = A*A'; +[d v] = eigs(A); +assert_checkalmostequal(norm(A*v-v*d),0,[], 1D-8); + +//general eigenvalue problem +B = rand(10,10); +B = B*B'; +A = rand(10,10); +[d v] = eigs(A,B); +assert_checkalmostequal(norm(A*v-B*v*d),0,[],1D-8); +[d v] = eigs(A,B,8,'SM'); +assert_checkalmostequal(norm(A*v-B*v*d),0,[],1D-8); +[d v] = eigs(A,B,8, 1); +assert_checkalmostequal(norm(A*v-B*v*d),0,[],1D-8); +[d v] = eigs(A,B, 8, %i); +assert_checkalmostequal(norm(A*v-B*v*d),0,[],1D-8); + +opts.cholB = %t; +R = chol(B); +[d v] = eigs(A, R, 8, 'LM', opts); +assert_checkalmostequal(norm(A*v-B*v*d),0,[],1D-8); +[d v] = eigs(A, R, 8,'SM', opts); +assert_checkalmostequal(norm(A*v-B*v*d),0,[],1D-8); +[d v] = eigs(A, R, 8, 1, opts); +assert_checkalmostequal(norm(A*v-B*v*d),0,[],1D-8); +[d v] = eigs(A, R, 8, %i, opts); +assert_checkalmostequal(norm(A*v-B*v*d),0,[],1D-8); + +A=A*A'; +[d v] = eigs(A,B); +assert_checkalmostequal(norm(A*v-B*v*d),0,[],1D-8); +[d v] = eigs(A,B,8,'SM'); +assert_checkalmostequal(norm(A*v-B*v*d),0,[],1D-8); +[d v] = eigs(A,B,8,1); +assert_checkalmostequal(norm(A*v-B*v*d),0,[],1D-8); +[d v] = eigs(A,B,8,%i); +assert_checkalmostequal(norm(A*v-B*v*d),0,[],1D-8); + + +//complex case +A1 = rand(10,10); +A2 = rand(10,10); +B1 = rand(10,10); +B2 = rand(10,10); +C1 = A1+%i*A2; +[d v] = eigs(C1); +assert_checkalmostequal(norm(C1*v-v*d),0,[],1D-8); +[d v] = eigs(C1,[], 8,'SM'); +assert_checkalmostequal(norm(C1*v-v*d),0,[],1D-8); +[d v] = eigs(C1, [],8, 1+%i); +assert_checkalmostequal(norm(C1*v-v*d),0,[],1D-8); +C2 = B1+%i*B2; +C2=C2*C2'; +[d v] = eigs(C1, C2); +assert_checkalmostequal(norm(C1*v-C2*v*d),0,[],1D-8); +[d v] = eigs(C1, C2, 8, 'SM'); +assert_checkalmostequal(norm(C1*v-C2*v*d),0,[],1D-8); +[d v] = eigs(C1, C2, 8, %i); +assert_checkalmostequal(norm(C1*v-C2*v*d),0,[],1D-8); -- 1.7.9.5