From 450075214600302b3e6453ecce8929f78d016478 Mon Sep 17 00:00:00 2001 From: Jiri Pittner Date: Wed, 10 Dec 2025 17:12:00 +0100 Subject: [PATCH] lanczos: skeleton initial commit --- Makefile.am | 2 +- la.h | 1 + lanczos.h | 155 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 157 insertions(+), 1 deletion(-) create mode 100644 lanczos.h diff --git a/Makefile.am b/Makefile.am index b619ecc..27756f1 100644 --- a/Makefile.am +++ b/Makefile.am @@ -1,5 +1,5 @@ lib_LTLIBRARIES = libla.la -include_HEADERS = LA_version.h tensor.h miscfunc.h simple.h vecmat3.h quaternion.h fortran.h cuda_la.h auxstorage.h davidson.h laerror.h mat.h qsort.h vec.h bisection.h diis.h la.h noncblas.h smat.h numbers.h bitvector.h fourindex.h la_traits.h la_random.h nonclass.h sparsemat.h sparsesmat.h csrmat.h conjgrad.h gmres.h matexp.h permutation.h polynomial.h contfrac.h graph.h reg.h regsurf.h +include_HEADERS = LA_version.h tensor.h miscfunc.h simple.h vecmat3.h quaternion.h fortran.h cuda_la.h auxstorage.h davidson.h lanczos.h laerror.h mat.h qsort.h vec.h bisection.h diis.h la.h noncblas.h smat.h numbers.h bitvector.h fourindex.h la_traits.h la_random.h nonclass.h sparsemat.h sparsesmat.h csrmat.h conjgrad.h gmres.h matexp.h permutation.h polynomial.h contfrac.h graph.h reg.h regsurf.h libla_la_SOURCES = simple.cc tensor.cc miscfunc.cc quaternion.cc vecmat3.cc vec.cc mat.cc smat.cc sparsemat.cc sparsesmat.cc csrmat.cc laerror.cc noncblas.cc numbers.cc bitvector.cc strassen.cc nonclass.cc cuda_la.cc fourindex.cc permutation.cc polynomial.cc contfrac.cc graph.cc la_random.cc reg.cc regsurf.cc nodist_libla_la_SOURCES = version.cc check_PROGRAMS = t test tX test_reg test_regsurf diff --git a/la.h b/la.h index b1663e1..48a2221 100644 --- a/la.h +++ b/la.h @@ -28,6 +28,7 @@ #include "bitvector.h" #include "conjgrad.h" #include "davidson.h" +#include "lanczos.h" #include "diis.h" #include "fourindex.h" #include "gmres.h" diff --git a/lanczos.h b/lanczos.h new file mode 100644 index 0000000..d5d1799 --- /dev/null +++ b/lanczos.h @@ -0,0 +1,155 @@ +/* + LA: linear algebra C++ interface library + Copyright (C) 2025 Jiri Pittner or + + This program is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see . +*/ +#ifndef _lanczos_h +#define _lanczos_h +#include "vec.h" +#include "smat.h" +#include "mat.h" +#include "sparsemat.h" +#include "nonclass.h" +#include "auxstorage.h" + +namespace LA { + +//Lanczos diagonalization of hermitean matrix +//matrix can be any class which has nrows(), ncols(), diagonalof(), issymmetric(), and gemv() available +//does not even have to be explicitly stored - direct CI +//therefore the whole implementation must be a template in a header + + + +template +extern void lanczos(const Matrix &bigmat, NRVec &eivals, NRVec *eivecs, const char *eivecsfile, + int nroots=1, const bool verbose=0, const double eps=1e-6, + const bool incore=1, int maxit=100, const int maxkrylov = 500, + void (*initguess)(NRVec &)=NULL, const typename LA_traits::normtype *target=NULL) +{ +if(!bigmat.issymmetric()) laerror("lanczos only for hermitean matrices"); +bool flag=0; +int n=bigmat.nrows(); +if ( n!= (int)bigmat.ncols()) laerror("non-square matrix in lanczos"); +if(eivals.size() vec1(n),vec2(n); +NRVec::normtype> r(maxkrylov); +NRVec *v0,*v1; +AuxStorage *s0,*s1; + +if(incore) + { + v0 = new NRVec[maxkrylov]; + v1 = new NRVec[maxkrylov]; + } +else + { + s0 = new AuxStorage; + s1 = new AuxStorage; + } + +int i,j; + +if(nroots>=maxkrylov) nroots =maxkrylov-1; +int nroot=0; +int oldnroot; + +//@@@@@@@@@@ +#if 0 + +//default guess based on lowest diagonal element of the matrix +if(initguess) initguess(vec1); +else + { + const T *diagonal = bigmat.diagonalof(vec2,false,true); + typename LA_traits::normtype t=1e100; int i,j; + vec1=0; + for(i=0, j= -1; i::realpart(diagonal[i])::realpart(diagonal[i]); j=i;} + vec1[j]=1; + } + +//initial step +vec1.normalize(); +bigmat.gemv(0,vec2,'n',1,vec1); //avoid bigmat.operator*(vec), since that needs to allocate another n-sized vector +if(incore) v0[0]=vec1; else s0->put(vec1,0); +if(incore) v1[0]=vec2; else s1->put(vec2,0); + + +//iterative Lanczos +int it=0; +for(k=2; k<=maxkrylov;++k) + { + diagonalize3(smallV,r,1,1,krylovsize+1,&smallSwork,1); //for symmetric matrix they have already been sorted to ascending order in lapack + if(target) //resort eigenpairs by distance from the target + { + NRVec::normtype> key(krylovsize+1); + for(int i=0; i<=krylovsize; ++i) key[i] = abs(r[i] - *target); + NRPerm p(krylovsize+1); + key.sort(0,p); + + NRVec::normtype> rp(maxkrylov); + NRMat smallVp(maxkrylov,maxkrylov); + for(int i=0; i<=krylovsize; ++i) + { + rp[i]= r[p[i+1]-1]; + for(int j=0; j<=krylovsize; ++j) smallVp(j,i) = smallV(j,p[i+1]-1); + } + r = rp; + smallV = smallVp; + } + + } +flag=1; +goto finished; + +converged: +AuxStorage::elementtype> *ev; +if(eivecsfile) ev = new AuxStorage::elementtype>(eivecsfile); +if(verbose) {std::cout << "Lanczos converged in "<get(vec2,j); + vec1.axpy(smallV(j,nroot),incore?v0[j]:vec2); + } + vec1.normalize(); + if(eivecs) eivecs[nroot]|=vec1; + if(eivecsfile) + { + ev->put(vec1,nroot); + } + } + } + +if(eivecsfile) delete ev; + +//@@@@ +#endif + +finished: +if(incore) {delete[] v0; delete[] v1;} +else {delete s0; delete s1;} + +if(flag) laerror("no convergence in lanczos"); +} + +}//namespace +#endif