src/LSLib/LSAux.cc
/* LU Back Substitution *********************************-*-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 38.  The code **
**  has been modified for GAMMA into C++.			**
**							 	**
**  The algorithm takes an input matrix A' in its LU decomp-	**
**  osition format where A' is some original matrix A with any	**
**  needed row permutations to attain the LU form. The row	**
**  permutations used to relate A to A' are stored in the	**
**  integer array indx.  The function the proceeds to solve	**
**								**  
**			A|x> = |b>				**
**								**  
**  for |x> by considering the problem as			** 
**								**  
**              A|x> = (LU)|x> = L(U|x>) |b> 			**
**								**  
**  and first solving for the vector |y> where |y> = U|x>	**
**								**  
**			L|y> = |b>				**
**								**
**  followed by solving for |x> by				**
**								**
**			U|x> = |y>				**
**								**
**  Due to the triagular nature of the arrays L and U it is	**
**  relatively easy to solve the vector equations involving	**
**  them.  The only further complexity in this algorithm is	**
**  that the input |b> must be first permuted to |b'> so that	**
**  its element ordering match that of A'.  Following the	**
**  solution in which |x'> is obtained it is then un-permuted	**
**  to match the original A ordering.  The first equation	**
**  to be solved (actually involving L' not L) appears as	**
**								**
** 		 [ 1   0   0   0 ] [y1] = [b1]			**
** 		 [L21  1   0   0 ] [y2] = |b2|			**
** 		 [L31 L32  1   0 ] [y3] = |b3|			**
** 		 [L41 L42 L43  1 ] [y4] = [b4]			**
**								**
**  because the diagonal elements of the matrix L are always 1.	**
**  The first element of |y> (actually |y'>) is given by	**
**								**
**			y1 = b1					**
**								**
**  and then subseqent elements of this vector are found by	**
**								**
**			      [	     ---	  ]		**
**		         1    |      \	 	  |		**
**	         yi = ------- | b  - /  <i|L|j>y  |		**
**		      <i|L|i> [  i   ---	j ]		**
**								**
*****************************************************************/
LU Back Substitution 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 38.  The code
has been modified for GAMMA into C++.

The algorithm takes an input matrix A' in its LU decomp
osition format where A' is some original matrix A with any
needed row permutations to attain the LU form. The row
permutations used to relate A to A' are stored in the
integer array indx.  The function the proceeds to solve



for |x> by considering the problem as



and first solving for the vector |y> where |y> = U|x>



followed by solving for |x> by



Due to the triagular nature of the arrays L and U it is
relatively easy to solve the vector equations involving
them.  The only further complexity in this algorithm is
that the input |b> must be first permuted to |b'> so that
its element ordering match that of A'.  Following the
solution in which |x'> is obtained it is then un-permuted
to match the original A ordering.  The first equation
to be solved (actually involving L' not L) appears as



because the diagonal elements of the matrix L are always 1.
The first element of |y> (actually |y'>) is given by

y1 = b1

and then subseqent elements of this vector are found by


LU Back Substitution 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 38 . The code has been modified for GAMMA into C ++ . The algorithm takes an input matrix A ' in its LU decomp osition format where A ' is some original matrix A with any needed row permutations to attain the LU form . The row permutations used to relate A to A ' are stored in the integer array indx . The function the proceeds to solve for |x > by considering the problem as and first solving for the vector |y > where |y > = U|x > followed by solving for |x > by Due to the triagular nature of the arrays L and U it is relatively easy to solve the vector equations involving them . The only further complexity in this algorithm is that the input |b > must be first permuted to |b ' > so that its element ordering match that of A ' . Following the solution in which |x ' > is obtained it is then un - permuted to match the original A ordering . The first equation to be solved ( actually involving L ' not L ) appears as because the diagonal elements of the matrix L are always 1 . The first element of |y > ( actually |y ' > ) is given by y1 = b1 and then subseqent elements of this vector are found by
void LU_backsub(matrix &ALU, int* indx, matrix& b)
  {
  int n = ALU.rows();			
  int ii=0, ip=0;
  complex sum=0;
  int i;
  for(i=0; i<n; i++)			
    {					
    ip = indx[i];			
    sum = b.get(ip,0);			
    b.put(b.get(i,0),ip,0);		
    if (ii)
      for (int j=0; j<i; j++)
        sum -= ALU.get(i,j)*b.get(j,0);
    else if(norm(sum))
      ii=1;
    b.put(sum,i,0);
    }
  for(i=n-1; i>=0; i--)			
    {					
    sum = b.get(i,0);			
    for(int j=i+1; j<n; j++)
      sum -= ALU.get(i,j)*b.get(j,0);
    b.put(sum/ALU.get(i,i),i,0);	
    }
  return;
  }
void LU_backsub ( matrix & ALU , int * indx , matrix & b ) { int n = ALU . rows ( ) ; int ii = 0 , ip = 0 ; complex sum = 0 ; int i ; for ( i = 0 ; i < n ; i ++ ) { ip = indx [ i ] ; sum = b . get ( ip , 0 ) ; b . put ( b . get ( i , 0 ) , ip , 0 ) ; if ( ii ) for ( int j = 0 ; j < i ; j ++ ) sum -= ALU . get ( i , j ) * b . get ( j , 0 ) ; else if ( norm ( sum ) ) ii = 1 ; b . put ( sum , i , 0 ) ; } for ( i = n - 1 ; i >= 0 ; i -- ) { sum = b . get ( i , 0 ) ; for ( int j = i + 1 ; j < n ; j ++ ) sum -= ALU . get ( i , j ) * b . get ( j , 0 ) ; b . put ( sum / ALU . get ( i , i ) , i , 0 ) ; } return ; }
LU Back Substitution *********************************-*-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 Back Substitution ********************************* - * - 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		ALU   : Input matrix A in LU form //			indx  : Integer array for permutations //			b     : Input column vector b // Return		      : // Note			      : // Get dimension of array //		First Solve L|y> = |b> For |y> // Works forwards, |y(0)> // solved first because L // is lower triangular // A is input permuted! Get // intial b's in correct order //		Now Solve U|x> = |y> For |x> // Works backwards, |x(n-1)> // solved first because U // is upper-triangular // This is now |x(i)>