From 05d04ceff80ef9eaf4888b5c0689a6d7e7897069 Mon Sep 17 00:00:00 2001
From: Chulwoo Jung <chulwoo@bnl.gov>
Date: Thu, 25 May 2017 12:30:47 -0400
Subject: [PATCH] Adding SimpleLanczos

---
 lib/algorithms/iterative/SimpleLanczos.h | 776 +++++++++++++++++++++++
 1 file changed, 776 insertions(+)
 create mode 100644 lib/algorithms/iterative/SimpleLanczos.h

diff --git a/lib/algorithms/iterative/SimpleLanczos.h b/lib/algorithms/iterative/SimpleLanczos.h
new file mode 100644
index 00000000..19ebbd26
--- /dev/null
+++ b/lib/algorithms/iterative/SimpleLanczos.h
@@ -0,0 +1,776 @@
+    /*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid 
+
+    Source file: ./lib/algorithms/iterative/ImplicitlyRestartedLanczos.h
+
+    Copyright (C) 2015
+
+Author: Peter Boyle <paboyle@ph.ed.ac.uk>
+Author: paboyle <paboyle@ph.ed.ac.uk>
+Author: Chulwoo Jung <chulwoo@bnl.gov>
+
+    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 2 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, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+    See the full license in the file "LICENSE" in the top level distribution directory
+    *************************************************************************************/
+    /*  END LEGAL */
+#ifndef GRID_LANC_H
+#define GRID_LANC_H
+
+#include <string.h> //memset
+
+#ifdef USE_LAPACK
+#ifdef USE_MKL
+#include<mkl_lapack.h>
+#else
+void LAPACK_dstegr(char *jobz, char *range, int *n, double *d, double *e,
+                   double *vl, double *vu, int *il, int *iu, double *abstol,
+                   int *m, double *w, double *z, int *ldz, int *isuppz,
+                   double *work, int *lwork, int *iwork, int *liwork,
+                   int *info);
+//#include <lapacke/lapacke.h>
+#endif
+#endif
+
+#include <Grid/algorithms/densematrix/DenseMatrix.h>
+#include <Grid/algorithms/iterative/EigenSort.h>
+
+// eliminate temorary vector in calc()
+#define MEM_SAVE
+
+namespace Grid {
+
+struct Bisection {
+
+#if 0
+static void get_eig2(int row_num,std::vector<RealD> &ALPHA,std::vector<RealD> &BETA, std::vector<RealD> & eig)
+{
+  int i,j;
+  std::vector<RealD> evec1(row_num+3);
+  std::vector<RealD> evec2(row_num+3);
+  RealD eps2;
+  ALPHA[1]=0.;
+  BETHA[1]=0.;
+  for(i=0;i<row_num-1;i++) {
+    ALPHA[i+1] = A[i*(row_num+1)].real();
+    BETHA[i+2] = A[i*(row_num+1)+1].real();
+  }
+  ALPHA[row_num] = A[(row_num-1)*(row_num+1)].real();
+  bisec(ALPHA,BETHA,row_num,1,row_num,1e-10,1e-10,evec1,eps2);
+  bisec(ALPHA,BETHA,row_num,1,row_num,1e-16,1e-16,evec2,eps2);
+
+  // Do we really need to sort here?
+  int begin=1;
+  int end = row_num;
+  int swapped=1;
+  while(swapped) {
+    swapped=0;
+    for(i=begin;i<end;i++){
+      if(mag(evec2[i])>mag(evec2[i+1]))	{
+	swap(evec2+i,evec2+i+1);
+	swapped=1;
+      }
+    }
+    end--;
+    for(i=end-1;i>=begin;i--){
+      if(mag(evec2[i])>mag(evec2[i+1]))	{
+	swap(evec2+i,evec2+i+1);
+	swapped=1;
+      }
+    }
+    begin++;
+  }
+
+  for(i=0;i<row_num;i++){
+    for(j=0;j<row_num;j++) {
+      if(i==j) H[i*row_num+j]=evec2[i+1];
+      else H[i*row_num+j]=0.;
+    }
+  }
+}
+#endif
+
+static void bisec(std::vector<RealD> &c,   
+		  std::vector<RealD> &b,
+		  int n,
+		  int m1,
+		  int m2,
+		  RealD eps1,
+		  RealD relfeh,
+		  std::vector<RealD> &x,
+		  RealD &eps2)
+{
+  std::vector<RealD> wu(n+2);
+
+  RealD h,q,x1,xu,x0,xmin,xmax; 
+  int i,a,k;
+
+  b[1]=0.0;
+  xmin=c[n]-fabs(b[n]);
+  xmax=c[n]+fabs(b[n]);
+  for(i=1;i<n;i++){
+    h=fabs(b[i])+fabs(b[i+1]);
+    if(c[i]+h>xmax) xmax= c[i]+h;
+    if(c[i]-h<xmin) xmin= c[i]-h;
+  }
+  xmax *=2.;
+
+  eps2=relfeh*((xmin+xmax)>0.0 ? xmax : -xmin);
+  if(eps1<=0.0) eps1=eps2;
+  eps2=0.5*eps1+7.0*(eps2);
+  x0=xmax;
+  for(i=m1;i<=m2;i++){
+    x[i]=xmax;
+    wu[i]=xmin;
+  }
+
+  for(k=m2;k>=m1;k--){
+    xu=xmin;
+    i=k;
+    do{
+      if(xu<wu[i]){
+	xu=wu[i];
+	i=m1-1;
+      }
+      i--;
+    }while(i>=m1);
+    if(x0>x[k]) x0=x[k];
+    while((x0-xu)>2*relfeh*(fabs(xu)+fabs(x0))+eps1){
+      x1=(xu+x0)/2;
+
+      a=0;
+      q=1.0;
+      for(i=1;i<=n;i++){
+	q=c[i]-x1-((q!=0.0)? b[i]*b[i]/q:fabs(b[i])/relfeh);
+	if(q<0) a++;
+      }
+//  	printf("x1=%0.14e a=%d\n",x1,a);
+      if(a<k){
+	if(a<m1){
+	  xu=x1;
+	  wu[m1]=x1;
+	}else {
+	  xu=x1;
+	  wu[a+1]=x1;
+	  if(x[a]>x1) x[a]=x1;
+	}
+      }else x0=x1;
+    }
+    printf("x0=%0.14e xu=%0.14e k=%d\n",x0,xu,k);
+    x[k]=(x0+xu)/2;
+  }
+}
+};
+
+/////////////////////////////////////////////////////////////
+// Implicitly restarted lanczos
+/////////////////////////////////////////////////////////////
+
+
+template<class Field> 
+    class SimpleLanczos {
+
+    const RealD small = 1.0e-16;
+public:       
+    int lock;
+    int get;
+    int Niter;
+    int converged;
+
+    int Nstop;   // Number of evecs checked for convergence
+    int Nk;      // Number of converged sought
+    int Np;      // Np -- Number of spare vecs in kryloc space
+    int Nm;      // Nm -- total number of vectors
+
+
+    RealD OrthoTime;
+
+    RealD eresid;
+
+    SortEigen<Field> _sort;
+
+    LinearOperatorBase<Field> &_Linop;
+
+    OperatorFunction<Field>   &_poly;
+
+    /////////////////////////
+    // Constructor
+    /////////////////////////
+    void init(void){};
+    void Abort(int ff, DenseVector<RealD> &evals,  DenseVector<DenseVector<RealD> > &evecs);
+
+    SimpleLanczos(
+				LinearOperatorBase<Field> &Linop, // op
+			       OperatorFunction<Field> & poly,   // polynmial
+			       int _Nstop, // sought vecs
+			       int _Nk, // sought vecs
+			       int _Nm, // spare vecs
+			       RealD _eresid, // resid in lmdue deficit 
+			       int _Niter) : // Max iterations
+      _Linop(Linop),
+      _poly(poly),
+      Nstop(_Nstop),
+      Nk(_Nk),
+      Nm(_Nm),
+      eresid(_eresid),
+      Niter(_Niter)
+    { 
+      Np = Nm-Nk; assert(Np>0);
+    };
+
+    /////////////////////////
+    // Sanity checked this routine (step) against Saad.
+    /////////////////////////
+    void RitzMatrix(DenseVector<Field>& evec,int k){
+
+      if(1) return;
+
+      GridBase *grid = evec[0]._grid;
+      Field w(grid);
+      std::cout<<GridLogMessage << "RitzMatrix "<<std::endl;
+      for(int i=0;i<k;i++){
+	_poly(_Linop,evec[i],w);
+	std::cout<<GridLogMessage << "["<<i<<"] ";
+	for(int j=0;j<k;j++){
+	  ComplexD in = innerProduct(evec[j],w);
+	  if ( fabs((double)i-j)>1 ) { 
+	    if (abs(in) >1.0e-9 )  { 
+	      std::cout<<GridLogMessage<<"oops"<<std::endl;
+	      abort();
+	    } else 
+	      std::cout<<GridLogMessage << " 0 ";
+	  } else { 
+	    std::cout<<GridLogMessage << " "<<in<<" ";
+	  }
+	}
+	std::cout<<GridLogMessage << std::endl;
+      }
+    }
+
+    void step(DenseVector<RealD>& lmd,
+	    DenseVector<RealD>& lme, 
+	    Field& last,
+	    Field& current, 
+	  	Field & next,
+		uint64_t k)
+    {
+	  if (lmd.size()<=k) lmd.resize(k+Nm);
+	  if (lme.size()<=k) lme.resize(k+Nm);
+      
+	
+      _poly(_Linop,current,next );   // 3. wk:=Avk−βkv_{k−1}
+      if(k>0){
+		next -= lme[k-1] * last;
+      }    
+//	std::cout<<GridLogMessage << "<last|next>" << innerProduct(last,next) <<std::endl;
+
+      ComplexD zalph = innerProduct(current,next); // 4. αk:=(wk,vk)
+      RealD     alph = real(zalph);
+
+      next = next - alph * current ;// 5. wk:=wk−αkvk
+//	std::cout<<GridLogMessage << "<current|next>" << innerProduct(current,next) <<std::endl;
+
+      RealD beta = normalise(next); // 6. βk+1 := ∥wk∥2. If βk+1 = 0 then Stop
+                                 // 7. vk+1 := wk/βk+1
+//	 norm=beta;
+
+	int interval = Nm/100+1;
+	if ((k%interval)==0) 
+	std::cout<<GridLogMessage << k << " : alpha = " << zalph << " beta "<<beta<<std::endl;
+      const RealD tiny = 1.0e-20;
+      if ( beta < tiny ) { 
+	std::cout<<GridLogMessage << " beta is tiny "<<beta<<std::endl;
+     }
+      lmd[k] = alph;
+      lme[k]  = beta;
+
+    }
+
+    void qr_decomp(DenseVector<RealD>& lmd,
+		   DenseVector<RealD>& lme,
+		   int Nk,
+		   int Nm,
+		   DenseVector<RealD>& Qt,
+		   RealD Dsh, 
+		   int kmin,
+		   int kmax)
+    {
+      int k = kmin-1;
+      RealD x;
+
+      RealD Fden = 1.0/hypot(lmd[k]-Dsh,lme[k]);
+      RealD c = ( lmd[k] -Dsh) *Fden;
+      RealD s = -lme[k] *Fden;
+      
+      RealD tmpa1 = lmd[k];
+      RealD tmpa2 = lmd[k+1];
+      RealD tmpb  = lme[k];
+
+      lmd[k]   = c*c*tmpa1 +s*s*tmpa2 -2.0*c*s*tmpb;
+      lmd[k+1] = s*s*tmpa1 +c*c*tmpa2 +2.0*c*s*tmpb;
+      lme[k]   = c*s*(tmpa1-tmpa2) +(c*c-s*s)*tmpb;
+      x        =-s*lme[k+1];
+      lme[k+1] = c*lme[k+1];
+      
+      for(int i=0; i<Nk; ++i){
+	RealD Qtmp1 = Qt[i+Nm*k  ];
+	RealD Qtmp2 = Qt[i+Nm*(k+1)];
+	Qt[i+Nm*k    ] = c*Qtmp1 - s*Qtmp2;
+	Qt[i+Nm*(k+1)] = s*Qtmp1 + c*Qtmp2; 
+      }
+
+      // Givens transformations
+      for(int k = kmin; k < kmax-1; ++k){
+
+	RealD Fden = 1.0/hypot(x,lme[k-1]);
+	RealD c = lme[k-1]*Fden;
+	RealD s = - x*Fden;
+	
+	RealD tmpa1 = lmd[k];
+	RealD tmpa2 = lmd[k+1];
+	RealD tmpb  = lme[k];
+
+	lmd[k]   = c*c*tmpa1 +s*s*tmpa2 -2.0*c*s*tmpb;
+	lmd[k+1] = s*s*tmpa1 +c*c*tmpa2 +2.0*c*s*tmpb;
+	lme[k]   = c*s*(tmpa1-tmpa2) +(c*c-s*s)*tmpb;
+	lme[k-1] = c*lme[k-1] -s*x;
+
+	if(k != kmax-2){
+	  x = -s*lme[k+1];
+	  lme[k+1] = c*lme[k+1];
+	}
+
+	for(int i=0; i<Nk; ++i){
+	  RealD Qtmp1 = Qt[i+Nm*k    ];
+	  RealD Qtmp2 = Qt[i+Nm*(k+1)];
+	  Qt[i+Nm*k    ] = c*Qtmp1 -s*Qtmp2;
+	  Qt[i+Nm*(k+1)] = s*Qtmp1 +c*Qtmp2;
+	}
+      }
+    }
+
+#ifdef USE_LAPACK
+#ifdef USE_MKL
+#define LAPACK_INT MKL_INT
+#else
+#define LAPACK_INT long long
+#endif
+    void diagonalize_lapack(DenseVector<RealD>& lmd,
+		     DenseVector<RealD>& lme, 
+		     int N1, // all
+		     int N2, // get
+		     GridBase *grid) {
+  const int size = Nm;
+  LAPACK_INT NN = N1;
+  double evals_tmp[NN];
+  double DD[NN];
+  double EE[NN];
+  for (int i = 0; i< NN; i++)
+    for (int j = i - 1; j <= i + 1; j++)
+      if ( j < NN && j >= 0 ) {
+        if (i==j) DD[i] = lmd[i];
+        if (i==j) evals_tmp[i] = lmd[i];
+        if (j==(i-1)) EE[j] = lme[j];
+      }
+  LAPACK_INT evals_found;
+  LAPACK_INT lwork = ( (18*NN) > (1+4*NN+NN*NN)? (18*NN):(1+4*NN+NN*NN)) ;
+  LAPACK_INT liwork =  3+NN*10 ;
+  LAPACK_INT iwork[liwork];
+  double work[lwork];
+  LAPACK_INT isuppz[2*NN];
+  char jobz = 'N'; // calculate evals only
+  char range = 'I'; // calculate il-th to iu-th evals
+  //    char range = 'A'; // calculate all evals
+  char uplo = 'U'; // refer to upper half of original matrix
+  char compz = 'I'; // Compute eigenvectors of tridiagonal matrix
+  int ifail[NN];
+  LAPACK_INT info;
+//  int total = QMP_get_number_of_nodes();
+//  int node = QMP_get_node_number();
+//  GridBase *grid = evec[0]._grid;
+  int total = grid->_Nprocessors;
+  int node = grid->_processor;
+  int interval = (NN/total)+1;
+  double vl = 0.0, vu = 0.0;
+  LAPACK_INT il = interval*node+1 , iu = interval*(node+1);
+  if (iu > NN)  iu=NN;
+  double tol = 0.0;
+    if (1) {
+      memset(evals_tmp,0,sizeof(double)*NN);
+      if ( il <= NN){
+        printf("total=%d node=%d il=%d iu=%d\n",total,node,il,iu);
+#ifdef USE_MKL
+        dstegr(&jobz, &range, &NN,
+#else
+        LAPACK_dstegr(&jobz, &range, &NN,
+#endif
+            (double*)DD, (double*)EE,
+            &vl, &vu, &il, &iu, // these four are ignored if second parameteris 'A'
+            &tol, // tolerance
+            &evals_found, evals_tmp, (double*)NULL, &NN,
+            isuppz,
+            work, &lwork, iwork, &liwork,
+            &info);
+        for (int i = iu-1; i>= il-1; i--){
+          printf("node=%d evals_found=%d evals_tmp[%d] = %g\n",node,evals_found, i - (il-1),evals_tmp[i - (il-1)]);
+          evals_tmp[i] = evals_tmp[i - (il-1)];
+          if (il>1) evals_tmp[i-(il-1)]=0.;
+        }
+      }
+      {
+         grid->GlobalSumVector(evals_tmp,NN);
+      }
+    } 
+// cheating a bit. It is better to sort instead of just reversing it, but the document of the routine says evals are sorted in increasing order. qr gives evals in decreasing order.
+}
+#undef LAPACK_INT 
+#endif
+
+
+    void diagonalize(DenseVector<RealD>& lmd,
+		     DenseVector<RealD>& lme, 
+		     int N2,
+		     int N1,
+		     GridBase *grid)
+    {
+
+#ifdef USE_LAPACK
+    const int check_lapack=0; // just use lapack if 0, check against lapack if 1
+
+    if(!check_lapack)
+	return diagonalize_lapack(lmd,lme,N2,N1,grid);
+
+//	diagonalize_lapack(lmd2,lme2,Nm2,Nm,Qt,grid);
+#endif
+    }
+
+#if 1
+    static RealD normalise(Field& v) 
+    {
+      RealD nn = norm2(v);
+      nn = sqrt(nn);
+      v = v * (1.0/nn);
+      return nn;
+    }
+
+    void orthogonalize(Field& w,
+		       DenseVector<Field>& evec,
+		       int k)
+    {
+      double t0=-usecond()/1e6;
+      typedef typename Field::scalar_type MyComplex;
+      MyComplex ip;
+
+      if ( 0 ) {
+	for(int j=0; j<k; ++j){
+	  normalise(evec[j]);
+	  for(int i=0;i<j;i++){
+	    ip = innerProduct(evec[i],evec[j]); // are the evecs normalised? ; this assumes so.
+	    evec[j] = evec[j] - ip *evec[i];
+	  }
+	}
+      }
+
+      for(int j=0; j<k; ++j){
+	ip = innerProduct(evec[j],w); // are the evecs normalised? ; this assumes so.
+	w = w - ip * evec[j];
+      }
+      normalise(w);
+      t0+=usecond()/1e6;
+      OrthoTime +=t0;
+    }
+
+    void setUnit_Qt(int Nm, DenseVector<RealD> &Qt) {
+      for(int i=0; i<Qt.size(); ++i) Qt[i] = 0.0;
+      for(int k=0; k<Nm; ++k) Qt[k + k*Nm] = 1.0;
+    }
+
+
+    void calc(
+		DenseVector<RealD>& eval,
+	      const Field& src,
+	      int& Nconv)
+      {
+
+	GridBase *grid = src._grid;
+//	assert(grid == src._grid);
+
+	std::cout<<GridLogMessage << " -- Nk = " << Nk << " Np = "<< Np << std::endl;
+	std::cout<<GridLogMessage << " -- Nm = " << Nm << std::endl;
+	std::cout<<GridLogMessage << " -- size of eval   = " << eval.size() << std::endl;
+	
+//	assert(c.size() && Nm == eval.size());
+	
+	DenseVector<RealD> lme(Nm);  
+	DenseVector<RealD> lmd(Nm);  
+
+	
+	Field current(grid);
+	Field last(grid);
+	Field next(grid);
+  
+	Nconv = 0;
+
+	RealD beta_k;
+  
+	// Set initial vector
+	// (uniform vector) Why not src??
+	//	evec[0] = 1.0;
+	current = src;
+	std:: cout<<GridLogMessage <<"norm2(src)= " << norm2(src)<<std::endl;
+	normalise(current);
+	std:: cout<<GridLogMessage <<"norm2(evec[0])= " << norm2(current) <<std::endl;
+	
+	// Initial Nk steps
+	OrthoTime=0.;
+	double t0=usecond()/1e6;
+	RealD norm; // sqrt norm of last vector
+
+	uint64_t iter=0;
+
+while(1){
+	std::vector <RealD> lme2(Nm);
+	std::vector <RealD> lmd2(Nm);
+	for(uint64_t k=0; k<Nm ; ++k,iter++){
+		step(lmd,lme,last,current,next,iter);
+		last=current;
+		current=next;
+	}
+	double t1=usecond()/1e6;
+	std::cout<<GridLogMessage <<"IRL::Initial steps: "<<t1-t0<< "seconds"<<std::endl; t0=t1;
+	std::cout<<GridLogMessage <<"IRL::Initial steps:OrthoTime "<<OrthoTime<< "seconds"<<std::endl;
+
+	  // getting eigenvalues
+	lmd2.resize(iter+2);
+	lme2.resize(iter+2);
+	for(uint64_t k=0; k<iter; ++k){
+		lmd2[k+1] = lmd[k];
+		lme2[k+2] = lme[k];
+	}
+	t1=usecond()/1e6;
+	std::cout<<GridLogMessage <<"IRL:: copy: "<<t1-t0<< "seconds"<<std::endl; t0=t1;
+{
+  int total = grid->_Nprocessors;
+  int node = grid->_processor;
+  int interval = (Nstop/total)+1;
+  int iu = (iter+1) - (interval*node+1);
+  int il = (iter+1) - (interval*(node+1));
+  RealD eps2;
+	  Bisection::bisec(lmd2,lme2,iter,il,iu,1e-16,1e-10, eval,eps2);
+//	  diagonalize(eval2,lme2,iter,Nk,grid);
+	for(int i=il;i<=iu;i++)
+	printf("eval[%d]=%0.14e\n",i,eval[i]);
+	t1=usecond()/1e6;
+	std::cout<<GridLogMessage <<"IRL:: diagonalize: "<<t1-t0<< "seconds"<<std::endl; t0=t1;
+}
+
+	for(uint64_t k=0; k<Nk; ++k){
+//	    eval[k] = eval2[k];
+    }
+}
+
+	}
+
+
+
+
+
+
+/**
+   There is some matrix Q such that for any vector y
+   Q.e_1 = y and Q is unitary.
+**/
+  template<class T>
+  static T orthQ(DenseMatrix<T> &Q, DenseVector<T> y){
+    int N = y.size();	//Matrix Size
+    Fill(Q,0.0);
+    T tau;
+    for(int i=0;i<N;i++){
+      Q[i][0]=y[i];
+    }
+    T sig = conj(y[0])*y[0];
+    T tau0 = fabs(sqrt(sig));
+    
+    for(int j=1;j<N;j++){
+      sig += conj(y[j])*y[j]; 
+      tau = abs(sqrt(sig) ); 	
+
+      if(abs(tau0) > 0.0){
+	
+	T gam = conj( (y[j]/tau)/tau0 );
+	for(int k=0;k<=j-1;k++){  
+	  Q[k][j]=-gam*y[k];
+	}
+	Q[j][j]=tau0/tau;
+      } else {
+	Q[j-1][j]=1.0;
+      }
+      tau0 = tau;
+    }
+    return tau;
+  }
+
+/**
+	There is some matrix Q such that for any vector y
+	Q.e_k = y and Q is unitary.
+**/
+  template< class T>
+  static T orthU(DenseMatrix<T> &Q, DenseVector<T> y){
+    T tau = orthQ(Q,y);
+    SL(Q);
+    return tau;
+  }
+
+
+/**
+	Wind up with a matrix with the first con rows untouched
+
+say con = 2
+	Q is such that Qdag H Q has {x, x, val, 0, 0, 0, 0, ...} as 1st colum
+	and the matrix is upper hessenberg
+	and with f and Q appropriately modidied with Q is the arnoldi factorization
+
+**/
+
+template<class T>
+static void Lock(DenseMatrix<T> &H, 	///Hess mtx	
+		 DenseMatrix<T> &Q, 	///Lock Transform
+		 T val, 		///value to be locked
+		 int con, 	///number already locked
+		 RealD small,
+		 int dfg,
+		 bool herm)
+{	
+  //ForceTridiagonal(H);
+
+  int M = H.dim;
+  DenseVector<T> vec; Resize(vec,M-con);
+
+  DenseMatrix<T> AH; Resize(AH,M-con,M-con);
+  AH = GetSubMtx(H,con, M, con, M);
+
+  DenseMatrix<T> QQ; Resize(QQ,M-con,M-con);
+
+  Unity(Q);   Unity(QQ);
+  
+  DenseVector<T> evals; Resize(evals,M-con);
+  DenseMatrix<T> evecs; Resize(evecs,M-con,M-con);
+
+  Wilkinson<T>(AH, evals, evecs, small);
+
+  int k=0;
+  RealD cold = abs( val - evals[k]); 
+  for(int i=1;i<M-con;i++){
+    RealD cnew = abs( val - evals[i]);
+    if( cnew < cold ){k = i; cold = cnew;}
+  }
+  vec = evecs[k];
+
+  ComplexD tau;
+  orthQ(QQ,vec);
+  //orthQM(QQ,AH,vec);
+
+  AH = Hermitian(QQ)*AH;
+  AH = AH*QQ;
+
+  for(int i=con;i<M;i++){
+    for(int j=con;j<M;j++){
+      Q[i][j]=QQ[i-con][j-con];
+      H[i][j]=AH[i-con][j-con];
+    }
+  }
+
+  for(int j = M-1; j>con+2; j--){
+
+    DenseMatrix<T> U; Resize(U,j-1-con,j-1-con);
+    DenseVector<T> z; Resize(z,j-1-con); 
+    T nm = norm(z); 
+    for(int k = con+0;k<j-1;k++){
+      z[k-con] = conj( H(j,k+1) );
+    }
+    normalise(z);
+
+    RealD tmp = 0;
+    for(int i=0;i<z.size()-1;i++){tmp = tmp + abs(z[i]);}
+
+    if(tmp < small/( (RealD)z.size()-1.0) ){ continue;}	
+
+    tau = orthU(U,z);
+
+    DenseMatrix<T> Hb; Resize(Hb,j-1-con,M);	
+	
+    for(int a = 0;a<M;a++){
+      for(int b = 0;b<j-1-con;b++){
+	T sum = 0;
+	for(int c = 0;c<j-1-con;c++){
+	  sum += H[a][con+1+c]*U[c][b];
+	}//sum += H(a,con+1+c)*U(c,b);}
+	Hb[b][a] = sum;
+      }
+    }
+	
+    for(int k=con+1;k<j;k++){
+      for(int l=0;l<M;l++){
+	H[l][k] = Hb[k-1-con][l];
+      }
+    }//H(Hb[k-1-con][l] , l,k);}}
+
+    DenseMatrix<T> Qb; Resize(Qb,M,M);	
+	
+    for(int a = 0;a<M;a++){
+      for(int b = 0;b<j-1-con;b++){
+	T sum = 0;
+	for(int c = 0;c<j-1-con;c++){
+	  sum += Q[a][con+1+c]*U[c][b];
+	}//sum += Q(a,con+1+c)*U(c,b);}
+	Qb[b][a] = sum;
+      }
+    }
+	
+    for(int k=con+1;k<j;k++){
+      for(int l=0;l<M;l++){
+	Q[l][k] = Qb[k-1-con][l];
+      }
+    }//Q(Qb[k-1-con][l] , l,k);}}
+
+    DenseMatrix<T> Hc; Resize(Hc,M,M);	
+	
+    for(int a = 0;a<j-1-con;a++){
+      for(int b = 0;b<M;b++){
+	T sum = 0;
+	for(int c = 0;c<j-1-con;c++){
+	  sum += conj( U[c][a] )*H[con+1+c][b];
+	}//sum += conj( U(c,a) )*H(con+1+c,b);}
+	Hc[b][a] = sum;
+      }
+    }
+
+    for(int k=0;k<M;k++){
+      for(int l=con+1;l<j;l++){
+	H[l][k] = Hc[k][l-1-con];
+      }
+    }//H(Hc[k][l-1-con] , l,k);}}
+
+  }
+}
+#endif
+
+
+ };
+
+}
+#endif
+