jiri/LA_library
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

*** empty log message ***

Browse Source
This commit is contained in:
jiri 2005-02-04 09:58:36 +00:00
parent 0dc73dad47
commit a235d4cb98
2 changed files with 45 additions and 102 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

7
davidson.h
View File

@ -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 <typename T, typename Matrix>
extern void davidson(const Matrix &bigmat, NRVec<T> *eivecs /*input-output*/, NRVec<T> &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<T> &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);

140
nonclass.cc
View File

@ -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<double> &A, NRMat<double> *B, double *det)
static void linear_solve_do(NRMat<double> &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<double> &A, NRMat<double> *B, double *det)
if (r>0 && B) laerror("singular matrix in lapack_gesv");
}
void linear_solve(NRMat<double> &A, NRMat<double> *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<double> &A, NRVec<double> &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<double> &a, NRMat<double> *b, double *det)
static void linear_solve_do(NRSMat<double> &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<double> &a, NRMat<double> *b, double *det)
}
void linear_solve(NRSMat<double> &a, NRMat<double> *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<double> &a, NRVec<double> &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<n; ++i)
{
for(j=i; j<n; ++j)
@ -612,8 +587,7 @@ for(j=i; j<n; ++j)
}
}
// .......... form the transpose of the upper triangle of inv(l)*a
// in the lower triangle of the array a ..........
// form the transpose of the upper triangle of inv(l)*a in the lower triangle of a
for(i=0; i<n; ++i)
{
for(j=i; j<n ; ++j)
@ -623,7 +597,7 @@ for(j=i; j<n ; ++j)
}
}
// .......... pre-multiply by inv(l) and overwrite ..........
//pre-multiply by l^-1
for(j=0; j<n ; ++j)
{
for(i=j;i<n;++i)
@ -641,38 +615,6 @@ for(i=1;i<n;++i) for(j=0; j<i; ++j) a(j,i)=a(i,j);
diagonalize(a,w,1,1,n);
//transform the eigenvectors back
/*
c
c this routine is a translation of the algol procedure rebaka,
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 forms the eigenvectors of a generalized
c symmetric eigensystem by back transforming those of the
c derived symmetric matrix determined by reduc.
c
c on input
c
c n is the order of the matrix system.
c
c b contains information about the similarity transformation
c (cholesky decomposition) used in the reduction by reduc
c in its strict lower triangle.
c
c dl contains further information about the transformation.
c
c a contains the eigenvectors to be back transformed
c in its columns
c
c on output
c
c a contains the transformed eigenvectors
c in its columns
c
c questions and comments should be directed to burton s. garbow,
c mathematics and computer science div, argonne national laboratory
c
*/
for(j=0; j<n; ++j)//eigenvector loop
{
for(int i=n-1; i>=0; --i)//component loop