58 template<
class vec_t=boost::numeric::ublas::vector<
double>,
59 class mat_t=boost::numeric::ublas::matrix<
double> >
class lanczos {
96 vec_t &eigen, vec_t &diag, vec_t &off_diag) {
113 for(i=1;i<size;i++) O2SCL_IX(w,i)=0.0;
115 for(i=0;i<size;i++) O2SCL_IX(v,i)=0;
119 for(i=0;i<size;i++) {
121 O2SCL_IX(w,i)=O2SCL_IX(v,i)/O2SCL_IX(off_diag,j-1);
122 O2SCL_IX(v,i)=-O2SCL_IX(off_diag,j-1)*t;
126 for(k=0;k<size;k++) O2SCL_IX(v,k)+=O2SCL_IX(prod,k);
127 O2SCL_IX(diag,j)=0.0;
128 O2SCL_IX(off_diag,j)=0.0;
130 for(k=0;k<size;k++) O2SCL_IX(diag,j)+=O2SCL_IX(w,k)*O2SCL_IX(v,k);
131 for(k=0;k<size;k++) O2SCL_IX(v,k)-=O2SCL_IX(diag,j)*O2SCL_IX(w,k);
132 for(k=0;k<size;k++) O2SCL_IX(off_diag,j)+=O2SCL_IX(v,k)*O2SCL_IX(v,k);
133 O2SCL_IX(off_diag,j)=sqrt(O2SCL_IX(off_diag,j));
136 if (j>=n_iter || O2SCL_IX(off_diag,j-1)==0.0) cont=
false;
139 for(k=0;k<size-1;k++) {
140 O2SCL_IX(b3,k+1)=O2SCL_IX(off_diag,k);
141 O2SCL_IX(eigen,k)=O2SCL_IX(diag,k);
143 O2SCL_IX(eigen,size-1)=O2SCL_IX(diag,size-1);
145 O2SCL_ERR2(
"Call to eigen_tdiag() in lanczos::",
176 double b,c,f,g,p,r,s,tst1,tst2;
180 for(
size_t ij=1;ij<n;ij++) {
181 O2SCL_IX(off_diag,ij-1)=O2SCL_IX(off_diag,ij);
183 O2SCL_IX(off_diag,n-1)=0.0;
190 while (done==
false && l<=((
int)n)) {
194 for(m=l;m<((int)n) && idone==
false;m++) {
195 tst1=fabs(O2SCL_IX(diag,m-1))+fabs(O2SCL_IX(diag,m));
196 tst2=tst1+fabs(O2SCL_IX(off_diag,m-1));
203 p=O2SCL_IX(diag,l-1);
205 if (m!=l && j==((
int)
td_iter)) {
209 O2SCL_ERR(
"No convergence in lanczos::eigen_tdiag()",
218 g=(O2SCL_IX(diag,l)-p)/(2.0*O2SCL_IX(off_diag,l-1));
221 g=O2SCL_IX(diag,m-1)-p+O2SCL_IX(off_diag,l-1)/
222 (g+(g>=0.0 ? fabs(r) : -fabs(r)));
228 for(
int ii=1;ii<=mml;ii++) {
231 f=s*O2SCL_IX(off_diag,i-1);
232 b=c*O2SCL_IX(off_diag,i-1);
234 O2SCL_IX(off_diag,i)=r;
240 O2SCL_IX(off_diag,m-1)=0.0;
247 g=O2SCL_IX(diag,i)-p;
248 r=(O2SCL_IX(diag,i-1)-g)*s+2.0*c*b;
250 O2SCL_IX(diag,i)=g+p;
257 O2SCL_IX(diag,l-1)-=p;
258 O2SCL_IX(off_diag,l-1)=g;
259 O2SCL_IX(off_diag,m-1)=0.0;
269 O2SCL_IX(diag,i-1)=p;
274 for(
int ii=2;ii<=l;ii++) {
276 if (p>=O2SCL_IX(diag,i-2)) {
280 O2SCL_IX(diag,i-1)=O2SCL_IX(diag,i-2);
284 if (skip==
false) i=1;
285 O2SCL_IX(diag,i-1)=p;
297 #ifndef DOXYGEN_INTERNAL
305 void product(
size_t n, mat_t &a, vec_t &w, vec_t &prod) {
308 O2SCL_IX(prod,i)=0.0;
310 O2SCL_IX(prod,i)+=O2SCL_IX2(a,i,j)*O2SCL_IX(w,j);