/*--------------------------------------------------------------------------*/
// 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);
// 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);
/*--------------------------------------------------------------------------*/
/*--------------------------------------------------------------------------*/
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);
/*--------------------------------------------------------------------------*/
/*--------------------------------------------------------------------------*/
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;
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
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;
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";
}
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;
- }
- }
}
}
{
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)
{
}
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)
{
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);
}
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)
{
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);
}
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;
}
-
/*
* 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
projects / scilab.git / commitdiff