Bug #12138 fixed - eigs(A,B) returned incorrect eigenvectors for dense matrices. 90/10090/6
Guillaume Horel [Tue, 8 Jan 2013 16:10:39 +0000 (17:10 +0100)]
test_run("arnoldi",["eigs","bug_12138"],["no_check_error_output"])

Change-Id: Ic4d02e9f6653ffab25c15b7c10b09d76e3be74e7

scilab/modules/arnoldi/includes/eigs.h
scilab/modules/arnoldi/includes/eigs_dependencies.h
scilab/modules/arnoldi/sci_gateway/c/sci_eigs.c
scilab/modules/arnoldi/src/c/eigs.c
scilab/modules/arnoldi/src/c/eigs_dependencies.c
scilab/modules/arnoldi/tests/nonreg_tests/bug_12138.dia.ref [new file with mode: 0644]
scilab/modules/arnoldi/tests/nonreg_tests/bug_12138.tst [new file with mode: 0644]

index 1586115..469cc7a 100644 (file)
@@ -14,6 +14,7 @@
 #define __EIGS_H__
 #include "doublecomplex.h"
 
+
 /**
  * @TODO add comment
  *
  * @param INFO_EUPD
  * @param eigenvalue
  * @param eigenvector
+ * @param eigenvalueC
+ * @param eigenvectorC
+ * @param RVEC
  * @return <ReturnValue>
  */
 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__ */
 /*--------------------------------------------------------------------------*/
-
index 9be06ff..e770488 100644 (file)
@@ -10,8 +10,8 @@
  *
  */
 /*--------------------------------------------------------------------------*/
-#ifndef __RTIMESRprime_H__
-#define __RTIMESRprime_H__
+#ifndef __PROCESS_DNEUPD_H__
+#define __PROCESS_DNEUPD_H__
 #include "doublecomplex.h"
 #include <string.h>
 #include <math.h>
 /**
  * @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__ */
 /*--------------------------------------------------------------------------*/
-
index e88aa82..2f96351 100644 (file)
@@ -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;
 }
-
index 1e78840..e42e9c5 100644 (file)
 
 /*--------------------------------------------------------------------------*/
 // 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;
 }
-
index 3c97f76..419def3 100644 (file)
 /*
  * 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 (file)
index 0000000..a40050a
--- /dev/null
@@ -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 (file)
index 0000000..c83c251
--- /dev/null
@@ -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);