lanczos: skeleton initial commit

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

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"

lanczos.h

@@ -0,0 +1,155 @@
+/*
+    LA: linear algebra C++ interface library
+    Copyright (C) 2025 Jiri Pittner <jiri.pittner@jh-inst.cas.cz> or <jiri@pittnerovi.com>
+
+    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 <http://www.gnu.org/licenses/>.
+*/
+#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 <typename T, typename Matrix>
+extern void lanczos(const Matrix &bigmat, NRVec<T> &eivals, NRVec<T> *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<T> &)=NULL, const typename LA_traits<T>::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()<nroots) laerror("too small eivals dimension in lanczos");
+
+NRVec<T> vec1(n),vec2(n);
+NRVec<typename LA_traits<T>::normtype> r(maxkrylov);
+NRVec<T> *v0,*v1;
+AuxStorage<T> *s0,*s1;
+
+if(incore)
+	{
+	v0 = new NRVec<T>[maxkrylov];
+	v1 = new NRVec<T>[maxkrylov];
+	}
+else
+	{
+	s0 = new AuxStorage<T>;
+	s1 = new AuxStorage<T>;
+	}
+
+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<T>::normtype t=1e100; int i,j;  
+	vec1=0;
+	for(i=0, j= -1; i<n; ++i) if(LA_traits<T>::realpart(diagonal[i])<t) {t=LA_traits<T>::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<typename LA_traits<T>::normtype> key(krylovsize+1);
+		for(int i=0; i<=krylovsize; ++i) key[i] = abs(r[i] - *target);
+		NRPerm<int> p(krylovsize+1);
+		key.sort(0,p);
+	
+		NRVec<typename LA_traits<T>::normtype> rp(maxkrylov);
+		NRMat<T> 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<typename LA_traits<T>::elementtype> *ev;
+if(eivecsfile) ev = new AuxStorage<typename LA_traits<T>::elementtype>(eivecsfile);
+if(verbose) {std::cout << "Lanczos converged in "<<it<<" iterations.\n"; std::cout.flush();}
+for(nroot=0; nroot<nroots; ++nroot)
+	{
+        eivals[nroot]=r[nroot];
+	if(eivecs)
+		{
+		vec1=0;
+		for(j=0; j<=krylovsize; ++j )
+			{
+            		if(!incore) s0->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

nonclass.cc

@@ -947,10 +947,10 @@ extern "C" void FORNAME(zggev)(const char *JOBVL, const char *JOBVR, const FINT
 //tridiagonal
 extern "C" void FORNAME(dsteqr)(const char *compz, const  FINT *n, double *d, double *d1, double *z, const FINT *ldz,  double *work, FINT *info);
 
-void diagonalize3(NRVec<double> &d, NRVec<double> &d1, NRMat<double> *v, const bool corder)
+void diagonalize3(NRVec<double> &d, NRVec<double> &d1, NRMat<double> *v, const bool corder, int n0)
 {
-	FINT n = d.size();
+	FINT n = n0? n0: d.size();
-	if(d1.size()!=n) laerror("inconsistent dimensions in diagonalize3");
+	if(d1.size()<n-1) laerror("inconsistent dimensions in diagonalize3");
 	if(v) {if(n!=v->nrows()||n!=v->ncols()) laerror("inconsistent dimensions in diagonalize3");}
 	d.copyonwrite();
 	d1.copyonwrite();

nonclass.h

@@ -184,7 +184,7 @@ extern void gdiagonalize(NRMat<std::complex<double> > &a, NRVec<double> &wr, NRV
                 NRMat<std::complex<double> > *b=NULL, NRVec<std::complex<double> > *beta=NULL);
 
 //diagonalization of symmetric real tridiagonal matrix
-extern void diagonalize3(NRVec<double> &d, NRVec<double> &d1, NRMat<double> *v, const bool corder=1);
+extern void diagonalize3(NRVec<double> &d, NRVec<double> &d1, NRMat<double> *v, const bool corder=1, int n0=0);
 
 //complex,real,imaginary parts of various entities
 template<typename T>