src/LSLib/LSAux.cc
/* LU Decomposition *************************************-*-c++-*-
**							 	**
**  This routine is taken from "Numerical Recipies", W.H.Press,	**
**  B.P.Flannery, S.A.Teukolsky, & W.T.Vetterling, Cambridge	**
**  University Press, 1989.  See pages 31 through 36.  The code **
**  has been modified for GAMMA into C++.			**
**							 	**
**  The algorithm takes an input matrix A & overwrites it with  **
**  the LU decompositions of A' where A' is A with any needed	**
**  permutations. Any permutations used	are stored in the	**
**  integer array indx.  The function returns the integer d	**
**  which will be +1 if an even	number of permutations was used **
**  or -1 if an odd number of permutations was necessary.  Sub- **
**  sequently, d can be used in taking the determinant of the	**
**  input matrix A (or LU).					** 
**								**  
**  The LU decomposition formulates the equation		**
**								**  
**		     	        A = LU				**  
**								**  
**  where L is lower triangular and U is upper triangular.	**
**								**
**  [A11 A12 A13 A14]   [L11  0   0   0 ] [U11 U12 U12 U14]	**
**  |A21 A22 A23 A24]   [L21 L22  0   0 ] [ 0  U22 U23 U24|	**
**  |A31 A32 A33 A34] = [L31 L32 L33  0 ] [ 0   0  U33 U34|	**
**  [A41 A42 A43 A44]   [L41 L42 L43 L44] [ 0   0   0  U44]	**
**								**
**  Both L and U have elements on the diagonal, but those of L	**
**  are always set to be 1: <i|L|i> = 1.			**
**								**
**  [A11 A12 A13 A14]   [ 1   0   0   0 ] [U11 U12 U12 U14]	**
**  |A21 A22 A23 A24]   [L21  1   0   0 ] [ 0  U22 U23 U24|	**
**  |A31 A32 A33 A34] = [L31 L32  1   0 ] [ 0   0  U33 U34|	**
**  [A41 A42 A43 A44]   [L41 L42 L43  1 ] [ 0   0   0  U44]	**
**								**
**  Neglecting the diagonal of L, both L and U can be overlaid	**
**  for storage in a single array (here overwriting A).		**
**								**
**                      [U11 U12 U12 U14]			**
**  			[L21 U22 U23 U24|			**
**  			[L31 L32 U33 U34|			**
**  			[L41 L42 L43 U44]			**
**								**
**  Finally, the algorithm uses Crouts method to determine the	**
**  elements of L and U.  This sets the diagonal elements of L	**
**  to 1 and follows a specific order in computation of the L	**
**  and U elements which allows A to be overwritten as L and U	**
**  are determined.  For Crouts method to be stable, partial	**
**  pivoting is used (row interchanges).			**
**								**
*****************************************************************/


The algorithm takes an input matrix A & overwrites it with
the LU decompositions of A' where A' is A with any needed
permutations. Any permutations used   are stored in the
integer array indx.  The function returns the integer d
which will be +1 if an even   number of permutations was used
or -1 if an odd number of permutations was necessary.  Sub
sequently, d can be used in taking the determinant of the
input matrix A (or LU).

The LU decomposition formulates the equation

A = LU





Both L and U have elements on the diagonal, but those of L
are always set to be 1:  = 1.



Neglecting the diagonal of L, both L and U can be overlaid
for storage in a single array (here overwriting A).



Finally, the algorithm uses Crouts method to determine the
elements of L and U.  This sets the diagonal elements of L
to 1 and follows a specific order in computation of the L
and U elements which allows A to be overwritten as L and U
are determined.  For Crouts method to be stable, partial
pivoting is used (row interchanges).
The algorithm takes an input matrix A & overwrites it with the LU decompositions of A ' where A ' is A with any needed permutations . Any permutations used are stored in the integer array indx . The function returns the integer d which will be + 1 if an even number of permutations was used or - 1 if an odd number of permutations was necessary . Sub sequently d can be used in taking the determinant of the input matrix A ( or LU ) . The LU decomposition formulates the equation A = LU Both L and U have elements on the diagonal but those of L are always set to be 1 : = 1 . Neglecting the diagonal of L both L and U can be overlaid for storage in a single array ( here overwriting A ) . Finally the algorithm uses Crouts method to determine the elements of L and U . This sets the diagonal elements of L to 1 and follows a specific order in computation of the L and U elements which allows A to be overwritten as L and U are determined . For Crouts method to be stable partial pivoting is used ( row interchanges ) .
int LU_decomp(matrix& A, int* indx)
  {
  int n = A.rows();			
  const double tiny = 1.0e-20;		
  double *vv;				
  vv = new double[n];
  int d = 1;				
  double big, temp;
  int i;
  for(i=0; i<n; i++)			
    {					
    big=0.0;				
    for (int j=0; j<n; j++)		
      {
      temp = norm(A.get(i,j));
      if (temp > big)
        big=temp;
      }
    if(big == 0.0)			
      {
      std::cout << "\nSingular matrix input\n";
      exit(-1);
      }
    vv[i]=1.0/big;
    }
  complex Uij, Lij;			
  int imax=0;
  double dum;
  complex dumz;
  for (int j=0; j<n; j++)		
    {
		       i-1
		       ---
 <i|U|j> = <i|A|j> - \   <i|L|k><k|U|j>
 		       /
		       ---
		       k=0
 This is equation 2.3.12 on page 34 of the referenced text and its use is
 restricted to terms above the diagonal.
    for (i=0; i<j; i++)
      {
      Uij = A.get(i,j);
      for(int k=0; k<i; k++)
        Uij -= A.get(i,k)*A.get(k,j);
      A.put(Uij,i,j);
      }
 Next implement the equation below to compute the values of U (i>j, i<n)
 For the case when i=j, this equation is the same as thh previous equation
 when the scaling factor is neglected (and this is done).
			       j-1
	        1	       ---
 <i|L|j> = ------- <i|A|j> - \   <i|L|k><k|U|j>
 	     <j|U|j>	       /
			       ---
			       k=0
    big=0.0;				 Used to look for largest pivot
    for (i=j; i<n; i++)
      {
      Lij = A.get(i,j);
      for(int k=0; k<j; k++)
	Lij -= A.get(i,k)*A.get(k,j);
      A.put(Lij,i,j);
      dum = vv[i] * norm(Lij);		 For looking at the pivot
      if (dum >= big)			 Use it if better that big
        {
	big=dum;
	imax=i;
	}
      }
    if (j != imax)			 This part interchanges the
      {					 rows of A via permutation
      for(int k=0; k<n; k++)		 The sign of d is switched
        {				 and the permutation stored
        dumz = A.get(imax,k);		 in the array indx
        A.put(A(j,k),imax,k);
	A.put(dumz,j,k);
	}
      d = -d;
      vv[imax] = vv[j];
      }
    indx[j] = imax;
    if(norm(A.get(j,j)) == 0.0)		 If <j|A|j> is zero make it
      A.put(tiny,j,j);			 small instead!
    if (j != n-1)			 Divide by the pivot element
      {					 If <j|A|j> was zero it has 
      dumz = 1.0/A.get(j,j);		 since been set to tiny so it
      for(i=j+1; i<n; i++)		 won't blow up!
        A.put(dumz*A(i,j),i,j);
      }
    }					 Return to treat next column
  delete [] vv;
  return d;
  }
int LU_decomp ( matrix & A , int * indx ) { int n = A . rows ( ) ; const double tiny = 1.0e-20 ; double * vv ; vv = new double [ n ] ; int d = 1 ; double big , temp ; int i ; for ( i = 0 ; i < n ; i ++ ) { big = 0.0 ; for ( int j = 0 ; j < n ; j ++ ) { temp = norm ( A . get ( i , j ) ) ; if ( temp > big ) big = temp ; } if ( big == 0.0 ) { std :: cout << " \n \n " ; exit ( -1 ) ; } vv [ i ] = 1.0 / big ; } complex Uij , Lij ; int imax = 0 ; double dum ; complex dumz ; for ( int j = 0 ; j < n ; j ++ ) { for ( i = 0 ; i < j ; i ++ ) { Uij = A . get ( i , j ) ; for ( int k = 0 ; k < i ; k ++ ) Uij -= A . get ( i , k ) * A . get ( k , j ) ; A . put ( Uij , i , j ) ; } big = 0.0 ; for ( i = j ; i < n ; i ++ ) { Lij = A . get ( i , j ) ; for ( int k = 0 ; k < j ; k ++ ) Lij -= A . get ( i , k ) * A . get ( k , j ) ; A . put ( Lij , i , j ) ; dum = vv [ i ] * norm ( Lij ) ; if ( dum >= big ) { big = dum ; imax = i ; } } if ( j != imax ) { for ( int k = 0 ; k < n ; k ++ ) { dumz = A . get ( imax , k ) ; A . put ( A ( j , k ) , imax , k ) ; A . put ( dumz , j , k ) ; } d = - d ; vv [ imax ] = vv [ j ] ; } indx [ j ] = imax ; if ( norm ( A . get ( j , j ) ) == 0.0 ) A . put ( tiny , j , j ) ; if ( j != n - 1 ) { dumz = 1.0 / A . get ( j , j ) ; for ( i = j + 1 ; i < n ; i ++ ) A . put ( dumz * A ( i , j ) , i , j ) ; } } delete [ ] vv ; return d ; }
LU Decomposition *************************************-*-c++
This routine is taken from "Numerical Recipies", W.H.Press,
B.P.Flannery, S.A.Teukolsky, & W.T.Vetterling, Cambridge
University Press, 1989.
LU Decomposition ************************************* - * - c ++ This routine is taken from " Numerical Recipies " W . H . Press B . P . Flannery S . A . Teukolsky & W . T . Vetterling Cambridge University Press 1989 .
// Input		A     : Input matrix //			indx  : Integer array for permutations // Return		A     : LU Decomposition of input A // 			indx  : Index of Permutations // Note			      : Assumes A is real, non-singular //				and square // Get dimension of A // Set this to a small number // Vector for row scaling factors // Permutation counter (1 even, -1 odd) // For each row of A find the // the largest element for pivoting  // and store corresponding scaling // factor in array vv // Cannot have a zero row in A! // Element of U and L // Begin going through each A column // First implement equation below to compute the values of U (except <i|U|i>) // //		       i-1 //		       --- // <i|U|j> = <i|A|j> - \   <i|L|k><k|U|j> // 		       / //		       --- //		       k=0 // // This is equation 2.3.12 on page 34 of the referenced text and its use is // restricted to terms above the diagonal. // Next implement the equation below to compute the values of U (i>j, i<n) // For the case when i=j, this equation is the same as thh previous equation // when the scaling factor is neglected (and this is done). // //			       j-1 //	        1	       --- // <i|L|j> = ------- <i|A|j> - \   <i|L|k><k|U|j> // 	     <j|U|j>	       / //			       --- //			       k=0 // Used to look for largest pivot // For looking at the pivot // Use it if better that big // This part interchanges the // rows of A via permutation // The sign of d is switched // and the permutation stored // in the array indx // If <j|A|j> is zero make it // small instead! // Divide by the pivot element // If <j|A|j> was zero it has  // since been set to tiny so it // won't blow up! // Return to treat next column