src/Matrix/n_matrix.cc
// Input	nmx	: A normal matrix (complex Hermitian!)
//		U	: A second normal matrix
// Output	void    : Both matrices, nmx & U are altered
//			  Now nmx will be tri-diagonal and U
//			  will be the transfomation array
// Note			: It is assumed that array U is zeroed!
Input  nmx     : A normal matrix (complex Hermitian!)
U       : A second normal matrix
Output void    : Both matrices, nmx & U are altered
Now nmx will be tri-diagonal and U
will be the transfomation array
Note                   : It is assumed that array U is zeroed!
Input nmx : A normal matrix ( complex Hermitian! ) U : A second normal matrix Output void : Both matrices nmx & U are altered Now nmx will be tri - diagonal and U will be the transfomation array Note : It is assumed that array U is zeroed!
void n_matrix::cred(n_matrix& U)
  {
  int i,j,l,dim=rows();
  double s;			
  double sw;			
  double t;
  double csf;			
  complex* w;
  w = new complex[dim];
  for(i=0;i<dim;i++) U(i,i)=1.0;
  for(l=0; l<dim-2; l++) 		
    {
    csf = 0.0;				
    for(i=l+2; i<dim; i++) 		
      csf += AbsNorm((*this)(i,l));	
    if(csf > 0.0)                       
      {                                 
      csf += AbsNorm((*this)(l+1,l));	
      sw = 0.0;
      for(i=l+2; i<dim; i++) 
	sw += square_norm((*this)(i,l)/csf);
      s = sw + square_norm((*this)(l+1,l)/csf);
      s = csf*sqrt(s);
      if((*this)(l+1,l)!=0)			
	{					
	t = norm((*this)(l+1,l));		
	mul(w[l+1],(*this)(l+1,l),(1.0+s/t));   
	}
      else w[l+1] = -s;				 -s
      s = sw+square_norm(w[l+1]/csf);
      s = sqrt(2/s)/csf;
      for(i=l+2; i<dim; i++) 			
	mul(w[i],(*this)(i,l),s);		
      w[l+1] = w[l+1]*s;
	for(i=l; i<dim; i++) 
	  {
	  register double ssi=0, ssr=0;
	  for(j=l+1; j<dim; j++)		
	    {					
	    complex& t = (*this)(j,i);		
	    register double wr = Re(w[j]);
	    register double wi = Im(w[j]);
	    register double tr = Re(t);
	    register double ti = Im(t);
	    ssr += wr*tr+wi*ti;
	    ssi += wr*ti-wi*tr;
	    }
	  for(j=l+1; j<dim; j++) 		
	    {					
	    complex& t = (*this)(j,i);		
	    register double wr = Re(w[j]);
	    register double wi = Im(w[j]);
	    Re(t, Re(t) -ssr*wr+ssi*wi);	
	    Im(t, Im(t) -ssr*wi-ssi*wr);	
	    }
	  }
	  for(i=l; i<dim; i++) 
	    {
	    register double ssi=0, ssr=0;
	    for(j=l+1; j<dim; j++) 		
	      {					
	      complex& t = (*this)(i,j);	
	      register double wr = Re(w[j]);
 	      register double wi = Im(w[j]);
	      register double tr = Re(t);
	      register double ti = Im(t);	
	      ssr += wr*tr-wi*ti;		
	      ssi += wr*ti+wi*tr;		
	      }
	    for(j=l+1; j<dim; j++) 		
	      {					
	      complex& t = (*this)(i,j);	
	      register double wr = Re(w[j]);
	      register double wi = Im(w[j]);
	      Re(t,Re(t)-ssr*wr-ssi*wi);	
	      Im(t,Im(t)+ssr*wi-ssi*wr);	
	      }
	    }
       for(i=0; i<dim; i++) 			
	{					
	register double ssi=0, ssr=0;		
	for(j=l+1; j<dim; j++)			
	  { 					
	  complex& t = U(i,j);			
	  register double wr = Re(w[j]);	
	  register double wi = Im(w[j]);	
	  register double tr = Re(t);		
	  register double ti = Im(t);		
	  ssr += wr*tr-wi*ti;			
	  ssi += wr*ti+wi*tr;			
	  }
	for(j=l+1; j<dim; j++) 			
	  {					
	  complex& t = U(i,j);			
	  register double wr = Re(w[j]);	
	  register double wi = Im(w[j]);	
	  Re(t, Re(t)-ssr*wr-ssi*wi);		
	  Im(t, Im(t)+ssr*wi-ssi*wr);		
	  }
	}
      }
    }
  delete [] w;
  }
void n_matrix :: cred ( n_matrix & U ) { int i , j , l , dim = rows ( ) ; double s ; double sw ; double t ; double csf ; complex * w ; w = new complex [ dim ] ; for ( i = 0 ; i < dim ; i ++ ) U ( i , i ) = 1.0 ; for ( l = 0 ; l < dim - 2 ; l ++ ) { csf = 0.0 ; for ( i = l + 2 ; i < dim ; i ++ ) csf += AbsNorm ( ( * this ) ( i , l ) ) ; if ( csf > 0.0 ) { csf += AbsNorm ( ( * this ) ( l + 1 , l ) ) ; sw = 0.0 ; for ( i = l + 2 ; i < dim ; i ++ ) sw += square_norm ( ( * this ) ( i , l ) / csf ) ; s = sw + square_norm ( ( * this ) ( l + 1 , l ) / csf ) ; s = csf * sqrt ( s ) ; if ( ( * this ) ( l + 1 , l ) != 0 ) { t = norm ( ( * this ) ( l + 1 , l ) ) ; mul ( w [ l + 1 ] , ( * this ) ( l + 1 , l ) , ( 1.0 + s / t ) ) ; } else w [ l + 1 ] = - s ; s = sw + square_norm ( w [ l + 1 ] / csf ) ; s = sqrt ( 2 / s ) / csf ; for ( i = l + 2 ; i < dim ; i ++ ) mul ( w [ i ] , ( * this ) ( i , l ) , s ) ; w [ l + 1 ] = w [ l + 1 ] * s ; for ( i = l ; i < dim ; i ++ ) { register double ssi = 0 , ssr = 0 ; for ( j = l + 1 ; j < dim ; j ++ ) { complex & t = ( * this ) ( j , i ) ; register double wr = Re ( w [ j ] ) ; register double wi = Im ( w [ j ] ) ; register double tr = Re ( t ) ; register double ti = Im ( t ) ; ssr += wr * tr + wi * ti ; ssi += wr * ti - wi * tr ; } for ( j = l + 1 ; j < dim ; j ++ ) { complex & t = ( * this ) ( j , i ) ; register double wr = Re ( w [ j ] ) ; register double wi = Im ( w [ j ] ) ; Re ( t , Re ( t ) - ssr * wr + ssi * wi ) ; Im ( t , Im ( t ) - ssr * wi - ssi * wr ) ; } } for ( i = l ; i < dim ; i ++ ) { register double ssi = 0 , ssr = 0 ; for ( j = l + 1 ; j < dim ; j ++ ) { complex & t = ( * this ) ( i , j ) ; register double wr = Re ( w [ j ] ) ; register double wi = Im ( w [ j ] ) ; register double tr = Re ( t ) ; register double ti = Im ( t ) ; ssr += wr * tr - wi * ti ; ssi += wr * ti + wi * tr ; } for ( j = l + 1 ; j < dim ; j ++ ) { complex & t = ( * this ) ( i , j ) ; register double wr = Re ( w [ j ] ) ; register double wi = Im ( w [ j ] ) ; Re ( t , Re ( t ) - ssr * wr - ssi * wi ) ; Im ( t , Im ( t ) + ssr * wi - ssi * wr ) ; } } for ( i = 0 ; i < dim ; i ++ ) { register double ssi = 0 , ssr = 0 ; for ( j = l + 1 ; j < dim ; j ++ ) { complex & t = U ( i , j ) ; register double wr = Re ( w [ j ] ) ; register double wi = Im ( w [ j ] ) ; register double tr = Re ( t ) ; register double ti = Im ( t ) ; ssr += wr * tr - wi * ti ; ssi += wr * ti + wi * tr ; } for ( j = l + 1 ; j < dim ; j ++ ) { complex & t = U ( i , j ) ; register double wr = Re ( w [ j ] ) ; register double wi = Im ( w [ j ] ) ; Re ( t , Re ( t ) - ssr * wr - ssi * wi ) ; Im ( t , Im ( t ) + ssr * wi - ssi * wr ) ; } } } } delete [ ] w ; }
Input  nmx     : A normal matrix (complex Hermitian!)
Input nmx : A normal matrix ( complex Hermitian! )
// norm of vector a(l+1,l) to a(dim-1,l) // norm of vector a(l+2,l) to a(dim-1,l) // Column scaling factor //                 Begin With Transformation Array As Identity //                      Householder Reduction Section /* We Loop Through Columns Of The Input Hermitian Matrix. First, All Elements
   That Are Not Tridiagonal (in the lth Column) Are Checked To See If They Are
   Non-Zero. If So, Then We Do Some Math For The Transformation. Note That On
   Successive Columns We Use The Updated Array With Alterations From The 
   Previous Column Adjustment In Place.                                      */ // Loop rows of input array //      Check For Non-Zero Tridiagonal Elements Lower Part of Column l // This will be column scaling factor // Get scaling factor, it checks col // l for !0 elements off tri-diagonal /*            Next Generate The Vector W For Column l If Needed
 
  This section assumes that  1. nll=<l|mx|l>   2. nend=<rows_|mx|rows_>+1
                             3. wl = <l|W>                                   */ // Do nothing if row "tridiagonal" // else we must transform. // Adjust scaling factor (off-diag) // <l+1|W> = //                    1+s //    <l+1|mx|l> = ------------ //                 |<l+1|mx|l>| // <l+1|W> = -s // Get scaled column l //  w[i] = s*<i|nmx|l> //                        Transform Matrix From The Left //		         This Will Alter The Input Array /* For the current column (l), the diagonal values should already be set.
   The element directly below (tri-diagonal) is determined, as are others
   in below that (all zero).  In the transform, the elements in the lower
   right block are altered to an intermediate state.  This is depicted
   below for the transformation on the 1st and 3rd columns of a 5x5 array
   of <j|hmx|i> (with i>=l, j>l) Are Adjusted.  Note That The Vector w Is
   left untouched in these loops.

       t h h h h     t h h h h       t t 0 0 0     t t 0 0 0  
       h h h h h     t x x x x       t t t 0 0     t t t 0 0
       h h h h h --> 0 x x x x       0 t t x x --> 0 t t x x
       h h h h h     0 x x x x       0 0 x x x     0 0 t y y
       h h h h h     0 x x x x       0 0 x x x     0 0 0 y y
                 l=0                           l=2

   No diagonal elements of the array are altered, nor are any of elements
   above the row with this same column index. Only the lower triangular
   elements of the active column (l) are set to their final values.  The
   rest (in a square block to the lower right) are altered to intermediate
   values corrected in the next transformation (From The Right).  This
   intermediate result will likely not be Hermitian as was the input array

   The most direct code to perform this transformation is show below, but this
   is not used in favor of code which is faster (albeit more complicated)
 
               for(i=l; i<dim; i++)
                 {
                 ss = 0.0;
	         for(j=l+1; j<dim; j++) ss += conj(w[j],(*this)(j,i));
	         for(j=l+1; j<dim; j++) (*this)(j,i) -= ss*w[j];
                 }                                                           */ // This loop calculates values // ssi & ssr for this i // Set t as <j|mx|i> // Loop adjusts the matrix elements // below row l, column i & right // Set t as <j|mx|i> // This will subtract off // parts from element <j|mx|i> //                        Transform Matrix From The Right //			  This Will Alter The Input Array /* For the current column (l), the vertical elements are already set.
   The transformation adjusts all of the elements to the right of the
   diagonal in the row with the same index as the column. Only the first
   (tri-diagonal) will be non-zero.  In addition, the block to the 
   lower right of the diagonal will be altered to an intermediate state,
   except for the next diagonal which will be properly set.  This is depicted
   below for the transformation on the 1st and 3rd columns of a 5x5 array

       t h h h h     t t 0 0 0       t t 0 0 0     t t 0 0 0  
       t x x x x     t t y y y       t t t 0 0     t t t 0 0
       0 x x x x --> 0 y y y y       0 t t x x --> 0 t t t 0
       0 x x x x     0 y y y y       0 0 t x x     0 0 t t y
       0 x x x x     0 y y y y       0 0 0 x x     0 0 0 y y
                 l=0                           l=2

   Intermediate results (array states prior to full tri-diagonalization)
   following the right transformation should be Hermitian.

   The most direct code to perform this transformation is show below, but this
   is not used in favor of code which is faster (albeit more complicated)
 
               for(i=l; i<dim; i++)
                 {
                 ss = 0.0;
	         for(j=l+1; j<dim; j++) ss += w[j]*(*this)(i,j);
	         for(j=l+1; j<dim; j++) (*this)(i,j) -= conj(w[j],ss);
                 }                                                           */ // This loop calculates the // values ssi and ssr, this i // Set t = <i|mx|j> = mij // Now just add wj*mij to ss // ssr+=Re(wj)Re(mij)-Im(wj)Im(mij) // ssi+=Re(wj)Im(mij)+Im(wj)Re(mij) // This loop adjusts // Set t = <i|mx|j> // Alters <i|mx|j> by // subtractiong conj(wj)*ss //                        Accumulate Transformations //              (Maintains U As Both Hermitian And Unitary) /* We Are Through Altering The Input Array For This Column Adjustment.  Now
   We Are Altering The Transformation Array U Using The Vector w. The Vector
   w Is Left Unaltered.  For Current Column "l", We Will Adjust All Array U
   Elements To The Right Of The Column.  This Is Depicted Below For A 5x5
   Array On The l=2 (3rd) Column, Where u Are Completed U Elements And x,y
   Represent Elements Not Quite Finished.

                           u u u x x    u u u u y  
                           u u u x x    u u u u y
                           u u u x x -> u u u u y
                           u u u x x    u u u u y
                           u u u x x    u u u u y

   The most direct code to perform this transformation is show below, but this
   is not used in favor of code which is faster (albeit more complicated) 

                    for(i=0; i<dim; i++) 
	               {
	               ss = 0.0;
	               for(j=l+1; j<dim; j++) ss += w[j]*U(i,j);
	               for(j=l+1; j<dim; j++) U(i,j) -= conj(w[j],ss);
	               }                                                     */ // Loop all rows of U // We first need to determine // ssi & ssr for this row from // upper triangle row elements // of U. //   Get t = <i|U|j> //   Re<j|w> //   Im<j|w> //   Re<i|U|j> //   Im<i|U|j> //   Adjust ssr column j //   Adjust ssi column j // Having ssr & ssj for row i, // adjust U row i upper triang.  //   Get t = <i|U|j> //   Re<j|w> //   Im<j|w> //   Adjust Re(<i|U|j>) //   Adjust Im(<i|U|j>)