src/Matrix/h_matrix.cc
/*************************************************************************
**									**
** The routine rred will convert a general complex Hermitian tri-	**
** diagonal array into a real symmetric tri-diagonal array.  This 	**
** result may be subsequently used for a complete diagonalization.	**
** The tranformation is diagonal and unitary.				**
**                                                                      **
**      Input           hmx     : Input hermitian array (this)		**
**                                Assumed to be tri-diagonal		**
**                      mxev    : Any previously accumulated            **
**                                transformations in a normal array.    **
**                                If no previous transformations this   **
**                                should be an identity matrix (normal) **
**      Output          hmx     : Real symmetric tri-diagonal array     **
**                      mxev    : All accumulated transformations.      **
**                                                                      **
**  The routine loops over each of the diagonal elements of the tri-    **
**  diagonal array starting with the 2nd column (index 1).  It then     **
**  looks at the element directly left which may be complex.  If found  **
**  this is then made real (as is its symmetric partner), the next two  **
**  off-diagonals set appropriatly, & the transformation matrix fixed.  **
**  Below is a depiction of what the input tri-diagonal array is doing  **
**  during the looping over a 5x5 array.                                **
**                                                                      **
**  x x . . .    x ? . . .    x r . . .    x r . . .    x r . . .       **
**    x x . .      X x . .      x A x .    r x ? . .    r x r x .       **
**      x x . ->     x x . ->   B x x . -> . B X x . -> . r x C . ->    **
**        x x          x x          x x          x x        D   x       **
**          x            x            x            x            x       **
**                                                                      **
**   Initial    Check Row 1     Adjust    Check Row 2     Adjust        **
**                                                                      **
**  Above the x's are just the original hermitian tri-diagonal elements **
**  and the big X's are the diagonals used in the loop.  I used r's to  **
**  depict altered elements which are real and A,B,C,D to depict those  **
**  modified but not yet to the final form.  Note that elements on the  **
**  lower diagonal (B, D, ...) are used as storage locations.           **
**									**
*************************************************************************/
The routine rred will convert a general complex Hermitian tri
diagonal array into a real symmetric tri-diagonal array.  This
result may be subsequently used for a complete diagonalization.
The tranformation is diagonal and unitary.

Input           hmx     : Input hermitian array (this)
Assumed to be tri-diagonal
mxev    : Any previously accumulated
transformations in a normal array.
If no previous transformations this
should be an identity matrix (normal)
Output          hmx     : Real symmetric tri-diagonal array
mxev    : All accumulated transformations.

The routine loops over each of the diagonal elements of the tri
diagonal array starting with the 2nd column (index 1).  It then
looks at the element directly left which may be complex.  If found
this is then made real (as is its symmetric partner), the next two
off-diagonals set appropriatly, & the transformation matrix fixed.
Below is a depiction of what the input tri-diagonal array is doing
during the looping over a 5x5 array.

x x .

Initial    Check Row 1     Adjust    Check Row 2     Adjust

Above the x's are just the original hermitian tri-diagonal elements
and the big X's are the diagonals used in the loop.  I used r's to
depict altered elements which are real and A,B,C,D to depict those
modified but not yet to the final form.  Note that elements on the
lower diagonal (B, D, ...) are used as storage locations.
The routine rred will convert a general complex Hermitian tri diagonal array into a real symmetric tri - diagonal array . This result may be subsequently used for a complete diagonalization . The tranformation is diagonal and unitary . Input hmx : Input hermitian array ( this ) Assumed to be tri - diagonal mxev : Any previously accumulated transformations in a normal array . If no previous transformations this should be an identity matrix ( normal ) Output hmx : Real symmetric tri - diagonal array mxev : All accumulated transformations . The routine loops over each of the diagonal elements of the tri diagonal array starting with the 2nd column ( index 1 ) . It then looks at the element directly left which may be complex . If found this is then made real ( as is its symmetric partner ) the next two off - diagonals set appropriatly & the transformation matrix fixed . Below is a depiction of what the input tri - diagonal array is doing during the looping over a 5x5 array . x x . Initial Check Row 1 Adjust Check Row 2 Adjust Above the x ' s are just the original hermitian tri - diagonal elements and the big X ' s are the diagonals used in the loop . I used r ' s to depict altered elements which are real and A B C D to depict those modified but not yet to the final form . Note that elements on the lower diagonal ( B D ... ) are used as storage locations .
void h_matrix::rred(_matrix* (&U), int newU)
  {
  int i;
  double normod;			
  complex sf;				
  complex *hii = data + rows_;		
  complex *hup = data + 1;		
  complex *hrt = hii + 1;		
  if(newU)				
    U=new n_matrix(rows_,rows_,complex0);
  complex *U00 = ((n_matrix*)U)->data;	
  complex *Uend = U00 + rows_*rows_;	
  if(newU)				
    {					
    complex *Uii = U00;			
    for(; Uii<Uend; Uii+=rows_+1)	
    *Uii = complex1;
    }
  complex *Uji;				
  int rl = rows_ - 1;			
  for(i=1; i<rows_; i++)		
    {
    normod = norm(*hup);		
    if(normod!=0) 			
      {
      sf = (*hup)/normod;		
      (*hup) = normod;			
      if(i < rows_-1) (*hrt) *= sf;	
      for(Uji=U00+i;Uji<Uend;Uji+=rows_)
        *Uji = conj(sf, *Uji);
      }
    hup = hrt;				
    hii += rl;				
    hrt = hii + 1;			
    rl--;				
    }
  }
void h_matrix :: rred ( _matrix * ( & U ) , int newU ) { int i ; double normod ; complex sf ; complex * hii = data + rows_ ; complex * hup = data + 1 ; complex * hrt = hii + 1 ; if ( newU ) U = new n_matrix ( rows_ , rows_ , complex0 ) ; complex * U00 = ( ( n_matrix * ) U ) -> data ; complex * Uend = U00 + rows_ * rows_ ; if ( newU ) { complex * Uii = U00 ; for ( ; Uii < Uend ; Uii += rows_ + 1 ) * Uii = complex1 ; } complex * Uji ; int rl = rows_ - 1 ; for ( i = 1 ; i < rows_ ; i ++ ) { normod = norm ( * hup ) ; if ( normod != 0 ) { sf = ( * hup ) / normod ; ( * hup ) = normod ; if ( i < rows_ - 1 ) ( * hrt ) *= sf ; for ( Uji = U00 + i ; Uji < Uend ; Uji += rows_ ) * Uji = conj ( sf , * Uji ) ; } hup = hrt ; hii += rl ; hrt = hii + 1 ; rl -- ; } }
The routine rred will convert a general complex Hermitian tri
diagonal array into a real symmetric tri-diagonal array.
The routine rred will convert a general complex Hermitian tri diagonal array into a real symmetric tri - diagonal array .
// Input        hmx     : A tri-diagonal Hermitian array (this) //              U	: Previously accumulated transitions //		newU	: Flag whether to begin with U=I // Output	void    : The two input arrays will be modified. //			  The array hmx will return real symmetric  //			  tridiagonal and the array U will be //			  adjusted to contain the additional	 //			  transformation from hmx in to hmx out. // Note			: This routine assumes that U is an //			  n_matrix and USES ITS INTERNAL STRUCTURE //			  It also assumes U is square and its //			  dimension matches that of hmx // For norm of off-diagonal // Transformation scaling factor // Diagonal element: hii = <1|hmx|1> // Elem. above diag: hup = <0|hmx|1> // Elem. right diag: hrt = <1|hmx|2> // Allocate array U if from scratch // U00 = <0|U|0> // Uend = <rows_|U|rows_> + 1 // If U is to begin from scratch // we must set it to I // <i|U|i> <- <0|U|0> // Initialize U to I matrix // Uji  = <j|U|i> // Length of row containing hii // Loop rows of hmx (after row 0) // If off-diag. non-zero <i|hmx|i-1> // then we transform it to be real // 	Scaling factor // 	Set off-diagonal to new real // 	Adjust the next off diagonal // 	Accumulate this transform // Right becomes up on next row  // Jump to next row's diagonal // New element right of diagonal hii // New row length