src/Matrix/n_matrix.cc
/*************************************************************************
**                                                                      **
** The routine tqli uses the QL algorithm to reduce an input tri-	**
** diagonal Hermitian array to a symmetric tridiagonal matrix.		**
**                                                                      **
**	Input:	nmx	: A symmetric tridiagonal matrix (this)		**
**		z	: Array of previous accumulated transformations	**
*	Output: a	: The tri-diagonalized matrix result		**
*		z	: Assumulated transformations.			**
**                                                                      **
*************************************************************************/
The routine tqli uses the QL algorithm to reduce an input tri
diagonal Hermitian array to a symmetric tridiagonal matrix.


The routine tqli uses the QL algorithm to reduce an input tri diagonal Hermitian array to a symmetric tridiagonal matrix .
void n_matrix::tqli(n_matrix& z, d_matrix& a)
  {
  int i,k,l,m,iter,n=rows();
  double dd;
  double s,r,p,g,f,c,b;
  double x;
  double *e;
  e = new double[n];
  complex xx;
  if(n>0)
    {
    for(i=0;i<n-1;i++) {	
      e[i]=zRe((*this)(i,i+1));
    };
    convert(&a);
    e[n-1] = 0.0;
    iter=0;
    for(l=0;l<n;l++)
      {						
      do {
        m=l;
        while (1) {		
          if (m==n-1) break;	
          dd=fabs(zRe(a(m,m)))+fabs(zRe(a(m+1,m+1)));
          if (fabs(e[m])+dd==dd) break;	
          m=m+1;
        };
        if (m!=l) {
          if (iter>10*n)
            {
            NMxerror(29,1);			
            NMxfatal(28);			
            }					
          iter=iter+1;
          g=(Re(a(l+1,l+1))-Re(a(l,l)))/(2.0*e[l]);
          r=sqrt(g*g+1.0);
          x=e[l];
          x=x/(g+sign(r,g));
          g=Re(a(m,m))-Re(a(l,l))+x;
          s=1.0;
          c=1.0;
          p=0.0;
          for(i=m-1; i>=l; i--) {
            f=s*e[i];
            b=c*e[i];
            if (fabs(f)>=fabs(g)) {
              c=g/f;
              r=sqrt(c*c+1.0);
              e[i+1]=f*r;
              s=1.0/r;
              c=c*s;
            } else {
              s=f/g;
              r=sqrt(s*s+1.0);
              e[i+1]=g*r;
              c=1.0/r;
              s=s*c;
            };
            g=Re(a(i+1,i+1))-p;
            r=(Re(a(i,i))-g)*s+2.0*c*b;
            p=s*r;
            Re(a(i+1,i+1),g+p);
            g=c*r-b;
            for (k=0;k<n;k++) {
              complex& zi =z(k,i);
              complex& zi1=z(k,i+1);
	      mul (xx, s, zi);
              sub( zi, c*zi, s*zi1);
              ::add( zi1, xx, c*zi1);
            };
          };
          Re(a(l,l),Re(a(l,l))-p);
          e[l]=g;
          e[m]=0.0;
        };
      } while (m!=l);		
    };
    }
  delete [] e;
  }
void n_matrix :: tqli ( n_matrix & z , d_matrix & a ) { int i , k , l , m , iter , n = rows ( ) ; double dd ; double s , r , p , g , f , c , b ; double x ; double * e ; e = new double [ n ] ; complex xx ; if ( n > 0 ) { for ( i = 0 ; i < n - 1 ; i ++ ) { e [ i ] = zRe ( ( * this ) ( i , i + 1 ) ) ; } ; convert ( & a ) ; e [ n - 1 ] = 0.0 ; iter = 0 ; for ( l = 0 ; l < n ; l ++ ) { do { m = l ; while ( 1 ) { if ( m == n - 1 ) break ; dd = fabs ( zRe ( a ( m , m ) ) ) + fabs ( zRe ( a ( m + 1 , m + 1 ) ) ) ; if ( fabs ( e [ m ] ) + dd == dd ) break ; m = m + 1 ; } ; if ( m != l ) { if ( iter > 10 * n ) { NMxerror ( 29 , 1 ) ; NMxfatal ( 28 ) ; } iter = iter + 1 ; g = ( Re ( a ( l + 1 , l + 1 ) ) - Re ( a ( l , l ) ) ) / ( 2.0 * e [ l ] ) ; r = sqrt ( g * g + 1.0 ) ; x = e [ l ] ; x = x / ( g + sign ( r , g ) ) ; g = Re ( a ( m , m ) ) - Re ( a ( l , l ) ) + x ; s = 1.0 ; c = 1.0 ; p = 0.0 ; for ( i = m - 1 ; i >= l ; i -- ) { f = s * e [ i ] ; b = c * e [ i ] ; if ( fabs ( f ) >= fabs ( g ) ) { c = g / f ; r = sqrt ( c * c + 1.0 ) ; e [ i + 1 ] = f * r ; s = 1.0 / r ; c = c * s ; } else { s = f / g ; r = sqrt ( s * s + 1.0 ) ; e [ i + 1 ] = g * r ; c = 1.0 / r ; s = s * c ; } ; g = Re ( a ( i + 1 , i + 1 ) ) - p ; r = ( Re ( a ( i , i ) ) - g ) * s + 2.0 * c * b ; p = s * r ; Re ( a ( i + 1 , i + 1 ) , g + p ) ; g = c * r - b ; for ( k = 0 ; k < n ; k ++ ) { complex & zi = z ( k , i ) ; complex & zi1 = z ( k , i + 1 ) ; mul ( xx , s , zi ) ; sub ( zi , c * zi , s * zi1 ) ; :: add ( zi1 , xx , c * zi1 ) ; } ; } ; Re ( a ( l , l ) , Re ( a ( l , l ) ) - p ) ; e [ l ] = g ; e [ m ] = 0.0 ; } ; } while ( m != l ) ; } ; } delete [ ] e ; }
The routine tqli uses the QL algorithm to reduce an input tri
diagonal Hermitian array to a symmetric tridiagonal matrix.
The routine tqli uses the QL algorithm to reduce an input tri diagonal Hermitian array to a symmetric tridiagonal matrix .
// copy to workspace // Loop for eigenvalues // look for an eigenvalue // all eigenvalues found // a new e.value found // Too many iterations // Unable to diagonalize // Almost surely bad input // form shift // QL step // accumulate transforms // all eigenvalues found