From a235d4cb98e4ac090b92bd5d93c5ac6a666bbd70 Mon Sep 17 00:00:00 2001 From: jiri Date: Fri, 4 Feb 2005 09:58:36 +0000 Subject: [PATCH] *** empty log message *** --- davidson.h | 7 +-- nonclass.cc | 140 +++++++++++++++------------------------------------- 2 files changed, 45 insertions(+), 102 deletions(-) diff --git a/davidson.h b/davidson.h index 5e46254..5840b3f 100644 --- a/davidson.h +++ b/davidson.h @@ -4,11 +4,12 @@ #include "sparsemat.h" #include "nonclass.h" -//Davidson diagonalization of real symmetric matrix +//Davidson diagonalization of real symmetric matrix (modified Lanczos) //matrix can be any class which has nrows(), ncols(), diagonalof() and NRVec::gemv() available //does not even have to be explicitly stored - direct CI -//n<0 highest eigenvalues, n>0 lowest eigenvalues export template -extern void davidson(const Matrix &bigmat, NRVec *eivecs /*input-output*/, NRVec &eivals, int nroots=1, const double accur=1e-6, int nits=50, const int ndvdmx = 500, const bool incore=1); +extern void davidson(const Matrix &bigmat, NRVec &eivals, void *eivecs, + int nroots=1, const bool verbose=0, const double eps=1e-6, + const bool incore=1, int maxit=100, const int maxkrylov = 500); diff --git a/nonclass.cc b/nonclass.cc index ef87064..a5e5e60 100644 --- a/nonclass.cc +++ b/nonclass.cc @@ -111,6 +111,7 @@ void laerror(const char *s1, const char *s2, const char *s3, const char *s4) if(s4) std::cerr << s4; } std::cerr << endl; +//@@@throw exit(1); } @@ -119,17 +120,14 @@ void laerror(const char *s1, const char *s2, const char *s3, const char *s4) ////////////////////// // A will be overwritten, B will contain the solutions, A is nxn, B is rhs x n -void linear_solve(NRMat &A, NRMat *B, double *det) +static void linear_solve_do(NRMat &A, double *B, const int nrhs, const int ldb, double *det) { int r, *ipiv; if (A.nrows() != A.ncols()) laerror("linear_solve() call for non-square matrix"); - if (B && A.nrows() != B->ncols()) laerror("incompatible matrices in linear_solve()"); A.copyonwrite(); - if (B) B->copyonwrite(); ipiv = new int[A.nrows()]; - r = clapack_dgesv(CblasRowMajor, A.nrows(), B ? B->nrows() : 0, A[0], A.ncols(), - ipiv, B ? B[0] : (double *)0, B ? B->ncols() : A.nrows()); + r = clapack_dgesv(CblasRowMajor, A.nrows(), nrhs, A[0], A.ncols(), ipiv, B , ldb); if (r < 0) { delete[] ipiv; laerror("illegal argument in lapack_gesv"); @@ -144,28 +142,36 @@ void linear_solve(NRMat &A, NRMat *B, double *det) if (r>0 && B) laerror("singular matrix in lapack_gesv"); } +void linear_solve(NRMat &A, NRMat *B, double *det) +{ +if (B && A.nrows() != B->ncols()) laerror("incompatible matrices in linear_solve()"); +if(B) B->copyonwrite(); +linear_solve_do(A,B?(*B)[0]:NULL,B?B->nrows() : 0, B?B->ncols():A.nrows(), det); +} -// Next routines are not available in clapack, fotran ones will b used with an +void linear_solve(NRMat &A, NRVec &B, double *det) +{ +if(A.nrows() != B.size()) laerror("incompatible matrices in linear_solve()"); +B.copyonwrite(); +linear_solve_do(A,&B[0],1,A.nrows(),det); +} + + +// Next routines are not available in clapack, fotran ones will be used with an // additional swap/transpose of outputs when needed extern "C" void FORNAME(dspsv)(const char *UPLO, const int *N, const int *NRHS, double *AP, int *IPIV, double *B, const int *LDB, int *INFO); -void linear_solve(NRSMat &a, NRMat *b, double *det) +static void linear_solve_do(NRSMat &a, double *b, const int nrhs, const int ldb, double *det) { int r, *ipiv; if (det) cerr << "@@@ sign of the determinant not implemented correctly yet\n"; - if (b && a.nrows() != b->ncols()) - laerror("incompatible matrices in symmetric linear_solve()"); a.copyonwrite(); - if (b) b->copyonwrite(); ipiv = new int[a.nrows()]; char U = 'U'; int n = a.nrows(); - int nrhs = 0; - if (b) nrhs = b->nrows(); - int ldb = b ? b->ncols() : a.nrows(); - FORNAME(dspsv)(&U, &n, &nrhs, a, ipiv, b?(*b)[0]:0, &ldb,&r); + FORNAME(dspsv)(&U, &n, &nrhs, a, ipiv, b, &ldb,&r); if (r < 0) { delete[] ipiv; laerror("illegal argument in spsv() call of linear_solve()"); @@ -181,6 +187,23 @@ void linear_solve(NRSMat &a, NRMat *b, double *det) } +void linear_solve(NRSMat &a, NRMat *B, double *det) +{ + if (B && a.nrows() != B->ncols()) + laerror("incompatible matrices in symmetric linear_solve()"); + if (B) B->copyonwrite(); +linear_solve_do(a,B?(*B)[0]:NULL,B?B->nrows() : 0, B?B->ncols():a.nrows(),det); +} + + +void linear_solve(NRSMat &a, NRVec &B, double *det) +{ + if (a.nrows() != B.size()) + laerror("incompatible matrices in symmetric linear_solve()"); + B.copyonwrite(); +linear_solve_do(a,&B[0],1,a.nrows(),det); +} + extern "C" void FORNAME(dsyev)(const char *JOBZ, const char *UPLO, const int *N, double *A, const int *LDA, double *W, double *WORK, const int *LWORK, int *INFO); @@ -547,56 +570,8 @@ if(n==0) n=m; if(n<0 || n>m) laerror("actual dimension in gendiagonalize out of range"); //transform the problem to usual diagonalization -/* -c -c this routine is a translation of the algol procedure reduc1, -c num. math. 11, 99-110(1968) by martin and wilkinson. -c handbook for auto. comp., vol.ii-linear algebra, 303-314(1971). -c -c this routine reduces the generalized symmetric eigenproblem -c ax=(lambda)bx, where b is positive definite, to the standard -c symmetric eigenproblem using the cholesky factorization of b. -c -c on input -c -c nm must be set to the row dimension of two-dimensional -c array parameters as declared in the calling program -c dimension statement. -c -c n is the order of the matrices a and b. if the cholesky -c factor l of b is already available, n should be prefixed -c with a minus sign. -c -c a and b contain the real symmetric input matrices. only the -c full upper triangles of the matrices need be supplied. if -c n is negative, the strict lower triangle of b contains, -c instead, the strict lower triangle of its cholesky factor l. -c -c dl contains, if n is negative, the diagonal elements of l. -c -c on output -c -c a contains in its full lower triangle the full lower triangle -c of the symmetric matrix derived from the reduction to the -c standard form. the strict upper triangle of a is unaltered. -c -c b contains in its strict lower triangle the strict lower -c triangle of its cholesky factor l. the full upper -c triangle of b is unaltered. -c -c dl contains the diagonal elements of l. -c -c ierr is set to -c zero for normal return, -c 7*n+1 if b is not positive definite. -c -c questions and comments should be directed to burton s. garbow, -c mathematics and computer science div, argonne national laboratory -c -c this version dated august 1983. -c -*/ -// .......... form l in the arrays b and dl .......... + +//cholesky decompose in b and dl for(i=0; i=0; --i)//component loop