src/Matrix/h_matrix.cc
//sosi- this function appears to be loosing memory!!!!!
sosi- this function appears to be loosing memory
sosi - this function appears to be loosing memory
void h_matrix::tqli(_matrix* (&U), _matrix* (&D), int newU)
  {
  if(!rows_) return;			
  D = new d_matrix(rows_,rows_);	
  double* odvals;
  odvals = new double[rows_];
  complex *hii = data; 			
  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 *dst = ((d_matrix*)D)->data;	
  complex *dii = dst;			
  int rl = rows_;			
  for(int ii=0; ii<rows_-1; ii++, dii++)
    {
    *dii = *hii;			
    odvals[ii] = zRe(*(hii+1));		
    hii += rl;				
    rl--;				
    }
  odvals[rows_-1] = 0.0;		
  *dii = *hii;				
  int i,l,m,iter;
  double dd,s,r,p,g,f,c,b,x;
  complex sxmki,cxmki,sxmkip1,cxmkip1;	
  complex *dmm, *dll, *diip1;		
  complex *u0i, *uki, *ukip1;		
  iter = 0;				
  for(l=0, dll=dst; l<rows_; l++,dll++)	
    {
    do
      {
      m = l;				
      dmm = dll; 			
      while(1)	 			
        {
        if(m == rows_-1) break;		
        dd = fabs(zRe(*dmm)) + fabs(zRe(*(dmm+1)));
        if(fabs(odvals[m])+dd == dd) break;
        m++;
        dmm++;				
        }
      if(m != l)
        {
        if(iter > 10*rows_) 		
          {
          HMxerror(19, 1);		
          HMxerror(20, 1);		
          HMxerror(55, 1);		
          HMxerror(53, 1);		
          HMxfatal(54);			
          }
        iter++;				
        g = (Re(*(dll+1))-Re(*dll))/(2.0*odvals[l]);
        r = sqrt(g*g + 1.0);
        x = odvals[l]/(g+sign(r,g));
        g = Re(*dmm) - Re(*dll) + x;
        s = 1.0;					
        c = 1.0;					
        p = 0.0;					
        for(i=m-1,dii=dmm-1,diip1=dmm; i>=l; i--,dii--,diip1--)
          {
          f = s*odvals[i];
          b = c*odvals[i];
          if(fabs(f) >= fabs(g))
            {
            c = g/f;
            r = sqrt(c*c+1.0);
            odvals[i+1] = f*r;
            s = 1.0/r;
            c = c*s;
            }
          else
            {
            s = f/g;
            r = sqrt(s*s + 1.0);
            odvals[i+1] = g*r;
            c = 1.0/r;
            s = s*c;
            }
          g = Re(*diip1) - p;
          r = (Re(*dii)-g)*s + 2.0*c*b;
          p = s*r;
          Re(*diip1, g+p);			
          g = c*r - b;
          u0i = u00 + i;			
          uki = u0i;				
          ukip1 = uki+1;			
          for(uki=u0i, ukip1=uki+1; uki<uend;
                       uki+=rows_,ukip1+=rows_)
            {
	    mul(sxmki, s, *uki);		
	    mul(cxmki, c, *uki);		
            mul(sxmkip1, s, *ukip1);		
            mul(cxmkip1, c, *ukip1);		
            sub(*uki, cxmki, sxmkip1);		
            ::add(*ukip1, sxmki, cxmkip1);
            }
          }
        Re((*dll), Re(*dll)-p);			
        odvals[l] = g;				
        odvals[m] = 0.0;			
        }
      } while (m!=l);				
    }
  delete [] odvals;
  }
void h_matrix :: tqli ( _matrix * ( & U ) , _matrix * ( & D ) , int newU ) { if ( ! rows_ ) return ; D = new d_matrix ( rows_ , rows_ ) ; double * odvals ; odvals = new double [ rows_ ] ; complex * hii = data ; 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 * dst = ( ( d_matrix * ) D ) -> data ; complex * dii = dst ; int rl = rows_ ; for ( int ii = 0 ; ii < rows_ - 1 ; ii ++ , dii ++ ) { * dii = * hii ; odvals [ ii ] = zRe ( * ( hii + 1 ) ) ; hii += rl ; rl -- ; } odvals [ rows_ - 1 ] = 0.0 ; * dii = * hii ; int i , l , m , iter ; double dd , s , r , p , g , f , c , b , x ; complex sxmki , cxmki , sxmkip1 , cxmkip1 ; complex * dmm , * dll , * diip1 ; complex * u0i , * uki , * ukip1 ; iter = 0 ; for ( l = 0 , dll = dst ; l < rows_ ; l ++ , dll ++ ) { do { m = l ; dmm = dll ; while ( 1 ) { if ( m == rows_ - 1 ) break ; dd = fabs ( zRe ( * dmm ) ) + fabs ( zRe ( * ( dmm + 1 ) ) ) ; if ( fabs ( odvals [ m ] ) + dd == dd ) break ; m ++ ; dmm ++ ; } if ( m != l ) { if ( iter > 10 * rows_ ) { HMxerror ( 19 , 1 ) ; HMxerror ( 20 , 1 ) ; HMxerror ( 55 , 1 ) ; HMxerror ( 53 , 1 ) ; HMxfatal ( 54 ) ; } iter ++ ; g = ( Re ( * ( dll + 1 ) ) - Re ( * dll ) ) / ( 2.0 * odvals [ l ] ) ; r = sqrt ( g * g + 1.0 ) ; x = odvals [ l ] / ( g + sign ( r , g ) ) ; g = Re ( * dmm ) - Re ( * dll ) + x ; s = 1.0 ; c = 1.0 ; p = 0.0 ; for ( i = m - 1 , dii = dmm - 1 , diip1 = dmm ; i >= l ; i -- , dii -- , diip1 -- ) { f = s * odvals [ i ] ; b = c * odvals [ i ] ; if ( fabs ( f ) >= fabs ( g ) ) { c = g / f ; r = sqrt ( c * c + 1.0 ) ; odvals [ i + 1 ] = f * r ; s = 1.0 / r ; c = c * s ; } else { s = f / g ; r = sqrt ( s * s + 1.0 ) ; odvals [ i + 1 ] = g * r ; c = 1.0 / r ; s = s * c ; } g = Re ( * diip1 ) - p ; r = ( Re ( * dii ) - g ) * s + 2.0 * c * b ; p = s * r ; Re ( * diip1 , g + p ) ; g = c * r - b ; u0i = u00 + i ; uki = u0i ; ukip1 = uki + 1 ; for ( uki = u0i , ukip1 = uki + 1 ; uki < uend ; uki += rows_ , ukip1 += rows_ ) { mul ( sxmki , s , * uki ) ; mul ( cxmki , c , * uki ) ; mul ( sxmkip1 , s , * ukip1 ) ; mul ( cxmkip1 , c , * ukip1 ) ; sub ( * uki , cxmki , sxmkip1 ) ; :: add ( * ukip1 , sxmki , cxmkip1 ) ; } } Re ( ( * dll ) , Re ( * dll ) - p ) ; odvals [ l ] = g ; odvals [ m ] = 0.0 ; } } while ( m != l ) ; } delete [ ] odvals ; }
sosi- this function appears to be loosing memory!!!!
sosi - this function appears to be loosing memory!!!!
// Input        hmx     : A tri-diagonal real symmetric array (this) //              U	: Pointer to mx to become eigenvectors mx //              D	: Pointer to mx to become diagonal mx //              newU    : Flag whether to begin with U=I // Output       void    : The matrices D and U are set to a diagonal //                        matrix of nmx eigenvalues and a matrix //                        the matrix of D eigenvectors respectively // Note			: The input array hmx is used only to supply //			  the diagonals & tri-diagonal elements // Note			: The pointer D must NOT point to any //			  allocated memory.  Memory will be allocated //			  herein and D set to be d_matrix* (diagonal) // Do nothing for 0x0 array //	       Copy Input Real Symmetric Tri-Diagonal Elements // Diagonals Go Into Diagonal Matrix D, Off-Diagonals Into Vector e[n] // Allocate space for eigenvalues //  double odvals[rows_];			// Array for off-diagonals // hii = <i|hmx|i> -> <0|hmx|0> // Allocate array U if from scratch // u00 = <0|U|0> // End of U: after <rows_|U|rows> // If U is to begin from scratch // we must set it to I // <i|U|i> <- <0|U|0> // Initialize U to I matrix // dst = <0|D|0> // dii = <i|D|i> -> <0|D|0> // Length of hmx row ii // Store diagonal of row ii // Store off-diagonal of row ii // Skip hii to next diagonal // Decrement row length // Last one set to zero // Last one <i|D|i> = <i|hmx|i> // Complex temporaries // Needed to get D elements // Needed to get U elements // Zero iteration count // Loop for eigenvalues // Start with row m // Pointer to <m|D|m> // Look for an eigenvalue // all eigenvalues found // A new eigenvalue found // Adjust to row <m|D|m> // Bail if too many iterations! //   Problems during diag. //   Too Many Iterations //   Unlikely Error //   Please Report Problem //   GAMMA WWW site // Increment iteration count // 				Form The Shift //                                        1/2 //     [ <l+1|hmx|l+1> - <l|hmx|l>       ] // r = | ------------------------- + 1.0 |        //     [       2*<l|hmx|l+1>             ] //                           This is The QL Step // Initialize s // Initialize c // Initialize p // <i+1|dmx|i+1> += g+p //          Accumulate Transforms In the Eigenvector Array U //   This Will Adjust Columns i and i+1 For This Interation in the Loop // Set <0|U|i> // Set <k|mkev|i> = <0|U|i> // Set <k|mkev|i+1> = <0|U|i+1> // sxmki = s*<k|U|i> (Need A Copy!) // cxmki = c*<k|U|i> // sxmkip1 = s*<k|U|i+1> // cxmkip1 = c*<k|U|i+1> // <k|U|i> = c*<k|U|i> - s*<k|U|i+1> //            add(*mkip1, sxmki, cxmkip1); // sosi - This is a little glitch.  GCC won't pick up on class complex "add" function //        instead it thinks it's a bad call to h_matrix add......  I don't even know //	  why the :: in front fixes things.... // <l|dmx|l> -= p // Set adjusted <l|hmx|l+1> // We've zeroed <m|hmx|m+1> // Loop till all eigenvalues found