src/Matrix/n_matrix.cc
/*
     this subroutine is a translation of a complex analogue of
     the algol procedure orthes, num. math. 12, 349-368(1968)
     by martin and wilkinson.
     handbook for auto. comp., vol.ii-linear algebra, 339-358(1971).

     given a complex general matrix, this subroutine
     reduces a submatrix situated in rows and columns
     low through igh to upper hessenberg form by
     unitary similarity transformations.

     on input
        low and igh are integers determined by the balancing
          subroutine  cbal.  if  cbal  has not been used,
          set low=1, igh=n.

        a contains the complex input matrix.

     on output

        a of the hessenberg matrix.  information
          about the unitary transformations used in the reduction
          is stored in the remaining triangles under the
          hessenberg matrix.

        ort contains further information about the
          transformations.  only elements low through igh are used.

*/
this subroutine is a translation of a complex analogue of
the algol procedure orthes, num. math.

given a complex general matrix, this subroutine
reduces a submatrix situated in rows and columns
low through igh to upper hessenberg form by
unitary similarity transformations.

on input
low and igh are integers determined by the balancing
subroutine  cbal.  if  cbal  has not been used,
set low=1, igh=n.

a contains the complex input matrix.

on output

a of the hessenberg matrix.  information
about the unitary transformations used in the reduction
is stored in the remaining triangles under the
hessenberg matrix.

ort contains further information about the
transformations.  only elements low through igh are used.
this subroutine is a translation of a complex analogue of the algol procedure orthes num . math . given a complex general matrix this subroutine reduces a submatrix situated in rows and columns low through igh to upper hessenberg form by unitary similarity transformations . on input low and igh are integers determined by the balancing subroutine cbal . if cbal has not been used set low = 1 igh = n . a contains the complex input matrix . on output a of the hessenberg matrix . information about the unitary transformations used in the reduction is stored in the remaining triangles under the hessenberg matrix . ort contains further information about the transformations . only elements low through igh are used .
complex* n_matrix::corth  (int low, int igh)
{
  int i,j,m;
  double f,g,h,sc;
  complex z;
  n_matrix& a=*this;
  int n=a.rows();
  complex *ort;
  ort=new complex[n];
  for (m=low+1; m<igh; m++) 
    {
      h=0;
      ort[m]=0;
      sc=0;
      for (i=m; i<=igh; i++) 
	sc+=AbsNorm( a(i,m-1) );
      if (sc>0.0) 
	{
	  for (i=igh;i>=m;i--) 
	    {
	      ort[i]= a(i,m-1) /sc;
	      h+=square_norm(ort[i]);
	    };
	  g=sqrt(h);
	  f=norm(ort[m]);
	  if (f==0.0)
	    {
	      ort[m]=g;
	      Re(a(m,m-1),sc);
	    }
	  else 
	    {
	      h+=f*g;
	      g/=f;
	      ort[m]*=(1.0+g);
	    };
	  for (j=m;j<n;j++)
	    {
	      z=0.0;
	      for (i=igh;i>=m;--i) 
		z += conj(ort[i],a(i,j));
	      z/=h;
	      for (i=m;i<=igh;i++) 
		a(i,j)-=z*ort[i];
	    };
	  for (i=0;i<=igh;i++) 
	    {
	      z=0.0;
	      for (j=igh;j>=m;--j) 
		z += ort[j]*a(i,j);
	      z/=h;
	      for (j=m;j<=igh;j++) 
		a(i,j)-= conj(ort[j],z);
	    };
	  ort[m]*=sc;
	  a(m,m-1) *= -g;
	};
    };
  return ort;
}
complex * n_matrix :: corth ( int low , int igh ) { int i , j , m ; double f , g , h , sc ; complex z ; n_matrix & a = * this ; int n = a . rows ( ) ; complex * ort ; ort = new complex [ n ] ; for ( m = low + 1 ; m < igh ; m ++ ) { h = 0 ; ort [ m ] = 0 ; sc = 0 ; for ( i = m ; i <= igh ; i ++ ) sc += AbsNorm ( a ( i , m - 1 ) ) ; if ( sc > 0.0 ) { for ( i = igh ; i >= m ; i -- ) { ort [ i ] = a ( i , m - 1 ) / sc ; h += square_norm ( ort [ i ] ) ; } ; g = sqrt ( h ) ; f = norm ( ort [ m ] ) ; if ( f == 0.0 ) { ort [ m ] = g ; Re ( a ( m , m - 1 ) , sc ) ; } else { h += f * g ; g /= f ; ort [ m ] *= ( 1.0 + g ) ; } ; for ( j = m ; j < n ; j ++ ) { z = 0.0 ; for ( i = igh ; i >= m ; -- i ) z += conj ( ort [ i ] , a ( i , j ) ) ; z /= h ; for ( i = m ; i <= igh ; i ++ ) a ( i , j ) -= z * ort [ i ] ; } ; for ( i = 0 ; i <= igh ; i ++ ) { z = 0.0 ; for ( j = igh ; j >= m ; -- j ) z += ort [ j ] * a ( i , j ) ; z /= h ; for ( j = m ; j <= igh ; j ++ ) a ( i , j ) -= conj ( ort [ j ] , z ) ; } ; ort [ m ] *= sc ; a ( m , m - 1 ) *= - g ; } ; } ; return ort ; }
this subroutine is a translation of a complex analogue of
the algol procedure orthes, num.
this subroutine is a translation of a complex analogue of the algol procedure orthes num .