//   ========================< FACTORED.CPP >===================
//   * basic class for factorisations                          *
//   * Description: Chapter 15                                 *
//   * Scientific C++ Building Numerical Libraries             *
//   *                the Object-Oriented Way                  *
//   * by G. Buzzi-Ferraris                                    *
//   * Addison-Wesley (1993)                                   *
//   ===========================================================

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

#include <alloc.h>

#include "utility.hpp"
#include "vector.hpp"
#include "matrix.hpp"
#include "factored.hpp"

const char *const Factored::ERROR=
       "\n========>>> Factored error!!!!!!\n";
const char *const Factored::ERR_DESTROYED=
       "Matrix Destroyed";
int Factored::count = 0;

//   ===========================================================
//   ================   protected functions   ==================
//   ===========================================================

//   *********************< Initialize >************************
//   * Purpose: Common function for initialisation             *
//   ***********************************************************
void Factored::Initialize(int rows,int columns)
   {
   count++;
   whoAmI = count;
   singular = 0;
   norm = 0.;
   factorizationStatus = UNFACTORED;
   matrixStatus = AVAILABLE;
   factored = 0;
   if(rows < 1 || columns < 1)
      {
      numRows = numColumns = size = 0;
      matrix = 0; 
      return;
      }
   numRows = rows;
   numColumns = columns;
   size = numRows * numColumns + 1;
   matrix = new float *[numRows + 1]; 
   if(!matrix)
      Error("%s%s%sInitialize",
      Factored::ERROR,ERR_SPACE,ERR_FUNCTION);
   matrix[0]=new float[size]; 
   if(!matrix[0])
      Error("%s%s%sInitialize",
      Factored::ERROR,ERR_SPACE,ERR_FUNCTION);
   matrix[1] = matrix[0];
   for(int i = 2;i <= numRows;i++)
      matrix[i] = matrix[i-1]+numColumns;
   }

//   ****************< FactoredInitialize >*********************
//   * Purpose: Factored initialisation                        *
//   ***********************************************************
void Factored::FactoredInitialize(void)
   {
   factored = new float *[numRows + 1]; 
   if(!factored)
      Error("%s%s%sFactoredInizialization",
      Factored::ERROR,ERR_SPACE,ERR_FUNCTION);
   factored[0]=new float[size]; 
   if(!factored[0])
      Error("%s%s%sFactoredInizialization",
      Factored::ERROR,ERR_SPACE,ERR_FUNCTION);
   factored[1] = factored[0];
   for(int i = 2;i <= numRows;i++)
      factored[i] = factored[i-1] + numColumns;
   memcpy(factored[0],matrix[0],size*sizeof(float));
   }

//   ********************< Deinitialize >***********************
//   * Purpose: Common function for deinitialisation           *
//   ***********************************************************
void Factored::Deinitialize(void)
   {
   FactoredDeinitialize();
   if(matrix != 0)
      {
      delete matrix[0];
      delete matrix; 
      }
   }

//   ***************< FactoredDeinitialize >*********************
//   * Purpose: Common function for factored                    *
//   *          deinitialisation                                *
//   ***********************************************************
void Factored::FactoredDeinitialize(void)
   {
   if(factored != 0)
      {
      delete factored[0];
      delete factored; 
      }
   factored = 0;
   }

//   ***********************< PrepCopy >************************
//   * Purpose: Preparing an assignment                        *
//   ***********************************************************
void Factored::PrepCopy(int rRows,int rColumns)
   {
   int who = whoAmI;
   SpecificDeinitialize();
   FactoredDeinitialize();
   singular = 0;
   norm = 0.;
   factorizationStatus = UNFACTORED;
   matrixStatus = AVAILABLE;
   factored = 0;
   if(numRows != rRows || numColumns != rColumns || matrix == 0)
      {
      Deinitialize();
      Initialize(rRows,rColumns);
      count--;
      }
   whoAmI = who;
   }

//   ****************< private constructor >********************
//   * Purpose: Sizing a Factored without initialisation       *
//   *          Serves internally                              *
//   * Example: Factored A('*',4,5);                           *
//   ***********************************************************
Factored::Factored(char,int rows,int columns)
   {
   Initialize(rows,columns);
   }

//   *******************< PrepOnlySolve >***********************
//   * Purpose: Preparing the factorisation                    *
//   * Description: The factored matrix is superimposed on a   *
//   *              matrix                                     *
//   ***********************************************************
void Factored::PrepOnlySolve(void)
   {
   if(factorizationStatus == UNFACTORED && 
      matrixStatus == DESTROYED)
        Error("%s%s%sPrepOnlySolve",
        Factored::ERROR,
        Factored::ERR_DESTROYED,ERR_FUNCTION);
   else if(factorizationStatus == UNFACTORED &&
          matrixStatus == AVAILABLE)
      {
      Norm();
      factored = matrix;
      matrix = 0;
      factorizationStatus = FACTORED;
      matrixStatus = DESTROYED;
      SpecificInitialize();
      Factorization();
      }
   else if(factorizationStatus == FACTORED && 
          matrixStatus == MODIFIED)
      { // factorised and modified
      Norm();
      FactoredDeinitialize();
      factored = matrix;
      matrix = 0;
      matrixStatus = DESTROYED;
      Factorization();
      }
   else if(factorizationStatus == FACTORED && 
          matrixStatus == AVAILABLE)
      { //factorised with PrepSolve
      if(matrix != 0)  // free matrix
        {
        delete matrix[0];
        delete matrix;
        }
      matrix = 0;
      matrixStatus = DESTROYED;
      }
   }

//   *********************< PrepSolve >*************************
//   * Purpose: Preparing the factorisation                    *
//   * Description: Factored and matrix are independent        *
//   ***********************************************************
void Factored::PrepSolve(void)
   {
   if(factorizationStatus == UNFACTORED && 
      matrixStatus == DESTROYED)
        Error("%s%s%sPrepOnlySolve",
        Factored::ERROR,
        Factored::ERR_DESTROYED,ERR_FUNCTION);
   else if(factorizationStatus == UNFACTORED && 
          matrixStatus == AVAILABLE)
      {
      FactoredInitialize();
      SpecificInitialize();
      Factorization();
      factorizationStatus = FACTORED;
      }
   else if(factorizationStatus == FACTORED && 
          matrixStatus == MODIFIED)
      {
      memcpy(factored[0],matrix[0],size*sizeof(float));
      Factorization();
      matrixStatus = AVAILABLE;
      }
   }

//   ************************< Norm >***************************
//   * Purpose: Calculating the Condition Number norm          *
//   * Description: NormR of the matrix from the beginning     *
//   ***********************************************************
void Factored::Norm(void)
   {
   float aus;
   norm = 0.;
   float *w = matrix[0];
   for(int i = 1;i <= numRows;i++)
       {
       aus=0.;
       for(int j = 1;j <= numColumns;j++)aus += Abs(*++w);
       if(aus > norm)norm = aus;
       }
   }

//   ===========================================================
//   =================   public functions   ====================
//   ===========================================================

//   ===========================================================
//   ******************   constructors   ***********************
//   ===========================================================

//   **********************< default >**************************
//   * Purpose: To define a Factored which can then            *
//   *          receive an assignment                          *
//   * Examples: Factored A;                                   *
//   *           A = B;                                        *
//   ***********************************************************
Factored::Factored(void)
   {
   Initialize(0,0);
   }

//   ****************< Copy-Initializer >***********************
//   * Description: Objects of the Factored or derived type    *
//   *              must not be copied or initialised          *
//   ***********************************************************
Factored::Factored(const Factored &rval)
   {
   Error("\n%sUses copy-initializer constructor%s"
      "\n con Factored N. %d",Factored::ERROR,rval.whoAmI);
   }

//   *************< sizing constructor >************************
//   * Purpose: Sizing a Factored and initialising it to 0     *
//   *          They must then be assigned to other values     *
//   * Example: Factored A(4,5);                               *
//   *          A(1,1)=3.; A(1,2)=7.;....                      *
//   ***********************************************************
Factored::Factored(int rows,int columns)
   {
   Initialize(rows,columns);
   if(size != 0)
      memset(matrix[0],0,size*sizeof(float));
   }

//   ********************< constructor >************************
//   * Purpose: Initialising a Factored from array             *
//   * Example: float w[]={1.,2.,3.,4.,5.,6.};                 *
//   *          Factored A(2,3,w);                             *
//   ***********************************************************
Factored::Factored(int rows,int columns,float *initvalues)
   {
   Initialize(rows,columns);
   float *w=matrix[0];
   ++w;
   if(numRows != 0)
      memcpy(w,initvalues,(size-1)*sizeof(float));
   }

//   *******************< from Matrix >*************************
//   * Purpose: Defining and initialising from Matrix          *
//   * Examples: Factored A=B;  // B Matrix                    *
//   *           Factored A=B+2.*C; // B and C Matrix          *
//   ***********************************************************
Factored::Factored(const Matrix &rval)
   {
   Initialize(rval.numRows,rval.numColumns);
   if(numRows != 0)
      memcpy(matrix[0],rval.matrix[0],size*sizeof(float));
   }

//   ********************< from submatrix >*********************
//   * Purpose: From a matrix with an identical start          *
//   * Example: Matrix B(4,5,A);                               *
//   ***********************************************************
Factored::Factored(int rows,int columns,const Matrix &rval)
   {
   if(rows < 1 || rows > rval.Rows())
       Error("%s%s%sSubFactored",
       Factored::ERROR,ERR_CHECK_DIMENSION,ERR_CONSTRUCTOR);
   if(columns < 1 || columns > rval.Columns())
       Error("%s%s%sSubFactored",
       Factored::ERROR,ERR_CHECK_DIMENSION,ERR_CONSTRUCTOR);
   Initialize(rows,columns);
   for(int row=1;row<=numRows;row++)
      memcpy(matrix[row]+1,rval.matrix[row]+1,
            numColumns*sizeof(float));
   }

//   *****************< from internal submatrix >***************
//   * Purpose: From internal submatrix                        *
//   * Example: Matrix B(m,n,i,j,A);      B(m,n) from i,j      *
//   ***********************************************************
Factored::Factored(int rows,int columns,
                 int irow,int jcol,const Matrix &rval)
   {
   if(rows < 1 || irow < 1
              || 
     irow > rval.Rows()
              ||
     rows > (rval.Rows()-irow+1))
       Error("%s%s%sSubFactored",
       Factored::ERROR,ERR_RANGE,ERR_CONSTRUCTOR);
   if(columns < 1 || jcol < 1 || 
     jcol > rval.Columns() ||
     columns>(rval.Columns()-jcol+1))
       Error("%s%s%sSubFactored",
       Factored::ERROR,ERR_RANGE,ERR_CONSTRUCTOR);
   Initialize(rows,columns);
   for(int row=1;row<=numRows;row++)
      memcpy(matrix[row]+1,rval.matrix[row+irow-1]+jcol,
            numColumns*sizeof(float));
   }

//   ===========================================================
//   ********************   destructor   ***********************
//   ===========================================================
Factored::~Factored()
   {
   if(factored != 0)
      {
      delete factored[0];
      delete factored; 
      }
   if(matrix != 0)
      {
      delete matrix[0];
      delete matrix; 
      }
   }

//   ===========================================================
//   ***********   Non-modifying access functions   ************
//   ===========================================================

//   *********************< GetValue >**************************
//   * Purpose: Receiving values with control                  *
//   * Example: float xf = A.GetValue(2,3);                    *
//   ***********************************************************
float Factored::GetValue(int row,int col) const
   {
   if(matrixStatus == DESTROYED)
      Error("%s%s",Factored::ERROR,ERR_DESTROYED);
   if( (row < 1 || row > numRows)
              ||
       (col < 1 || col > numColumns) )
           Error("%s%s%sGetValue",
           Factored::ERROR,ERR_RANGE,ERR_FUNCTION);
   return matrix[row][col];
   }

//   **********************< GetRow >***************************
//   * Purpose: Receiving the values of a row                  *
//   * Example: Vector v = A.GetRow(3);                        *
//   ***********************************************************
Vector Factored::GetRow(int i) const
   {
   if(matrixStatus == DESTROYED)
      Error("%s%s",Factored::ERROR,ERR_DESTROYED);
   if(i < 1 || i > numRows)
      Error("%s%s%sGetRow",
      Factored::ERROR,ERR_RANGE,ERR_FUNCTION);
   Vector result('*',numColumns);
   float *w=matrix[i];
   float *r=result.vector;
   memcpy(r,w,(numColumns+1)*sizeof(float));
   return result;
   }

//   *********************< GetColumn >*************************
//   * Purpose: Receiving the values of a column               *
//   * Example: Vector v = A.GetVector(2);                     *
//   ***********************************************************
Vector Factored::GetColumn(int j) const
   {
   if(matrixStatus == DESTROYED)
      Error("%s%s",Factored::ERROR,ERR_DESTROYED);
   if(j < 1 || j > numColumns)
      Error("%s%s%sGetColumn",
      Factored::ERROR,ERR_RANGE,ERR_FUNCTION);
   Vector result('*',numRows);
   float *r = result.vector;
   for(int k = 1;k <= numRows;k++)
       *++r = matrix[k][j];
   return result;
   }

//   ===========================================================
//   **************    Modifying access functions **************
//   ===========================================================

//   ********************< operator () >************************
//   * Purpose: Receiving and assigning values with control    *
//   * Example: x = A(1,5); A(3,7) = 5.;                       *
//   ***********************************************************
float &Factored::operator () (int row,int col)
   {
   if(matrixStatus == DESTROYED)
      Error("%s%s",Factored::ERROR,ERR_DESTROYED);
   if(factorizationStatus == FACTORED)
      matrixStatus = MODIFIED;
   if( (row < 1 || row > numRows)
              ||
     (col < 1 || col > numColumns) )
       Error("%s",Factored::ERROR);
   return matrix[row][col];
   }

//   *********************< SetValue >**************************
//   * Purpose: Assigning values with control                  *
//   ***********************************************************
void Factored::SetValue(int row,int col,float val)
   {
   (*this)(row,col) = val;
   }

//   **********************< SetRow >***************************
//   * Purpose: Assigning a row with a Vector                  *
//   * Example: A.SetRow(3,v);                                 *
//   ***********************************************************
void Factored::SetRow(int i,const Vector &rval)
   {
   if(matrixStatus == DESTROYED)
      Error("%s%s",Factored::ERROR,ERR_DESTROYED);
   if(factorizationStatus == FACTORED)
      matrixStatus = MODIFIED;
   if( (i < 1 || i > numRows)
             ||
      (rval.dimensions != numColumns ) )
          Error("%s%s%sSetRow",
          Factored::ERROR,ERR_RANGE,ERR_FUNCTION);
   memcpy(matrix[i]+1,rval.vector+1,
         numColumns*sizeof(float));
   }

//   *********************< SetColumn >*************************
//   * Purpose: Assigning a column with a Vector               *
//   * Example: A.SetColumn(3,v);                              *
//   ***********************************************************
void Factored::SetColumn(int j,const Vector &rval)
   {
   if(matrixStatus == DESTROYED)
      Error("%s%s",Factored::ERROR,ERR_DESTROYED);
   if(factorizationStatus == FACTORED)
      matrixStatus = MODIFIED;
   if( (j < 1 || j > numColumns)
             ||
      (rval.dimensions != numRows) )
          Error("%s%s%sSetColumn",
          Factored::ERROR,ERR_RANGE,ERR_FUNCTION);
   float *r = rval.vector;
   for(int k = 1;k <= numRows;k++)
       matrix[k][j] = *++r;
   }

//   ===========================================================
//   ***************   assignment operators   ******************
//   ===========================================================

//   ******************< CopyFromMatrix >***********************
//   * Purpose: Assignment of a Matrix                         *
//   * Description: Used internally                            *
//   ***********************************************************
void Factored::CopyFromMatrix(const Matrix &rval)
   {
   PrepCopy(rval.numRows,rval.numColumns);
   if(numRows != 0)
      memcpy(matrix[0],rval.matrix[0],size*sizeof(float));
   }

//   ***************< CommuteMatrixToFactored >*****************
//   * Purpose: Transforming a Matrix in Factored              *
//   * Description: Destroys the Matrix                        *
//   ***********************************************************
void CommuteMatrixToFactored(Matrix *lval,Factored *rval)
   {
   Delete(rval);
   rval->numRows = lval->numRows;
   rval->numColumns = lval->numColumns;
   rval->size = lval->size;
   rval->matrix = lval->matrix;
   lval->numRows = 0;
   lval->numColumns = 0;
   lval->size = 0;
   lval->matrix = 0;
   }

//   ===========================================================
//   *              Non-modifying functions                    *
//   ===========================================================

//   ***********************< Print >***************************
//   * Purpose: Printing                                       *
//   ***********************************************************
void Factored::Print(char *msg)
   {
   fprintf(bzzFileOut,"\nFactored No.%d",whoAmI);
   if(*msg)fprintf(bzzFileOut,"  %s",msg);
   if(matrixStatus == DESTROYED)return;
   fprintf(bzzFileOut,"\nrows %d columns %d\n",
          numRows,numColumns);
   for(int row = 1;row <= numRows;row++)
      {
      for(int column = 1;column <= numColumns;column++)
        fprintf(bzzFileOut," %12.4e",matrix[row][column]);
      fprintf(bzzFileOut,"\n");
      }
   }

//   ===========================================================
//   *                 Modifying functions                             *
//   ===========================================================

//   **********************< Delete >***************************
//   * Purpose: Eliminating a factored of no use               *
//   *          without it leaving the purpose                 *
//   * Example: Delete(&A);                                    *
//   ***********************************************************
void Delete(Factored *result)
   {
   int who = result->whoAmI;
   result->SpecificDeinitialize();
   result->Deinitialize();
   result->Initialize(0,0);
   Factored::count--;
   result->whoAmI = who;
   }

//   ===========================================================
//   ==========   Functions for linear algebra   ===============
//   ===========================================================

//   ***********************< Solve >***************************
//   * Purpose: Solving a linear system for a vector           *
//   * Description: The matrix A is factorised                 *
//   *              the solution replaces bx                   *
//   * Example: Solve(&A,&bx);                                 *
//   ***********************************************************
void Solve(Factored *A,Vector *bx)
   {
   if(bx->Size() != A->numRows)
      Error("%s%s%sOnlySolve",
      Factored::ERROR,ERR_CHECK_DIMENSION,ERR_FUNCTION);
   A->PrepOnlySolve();
   A->Solution(bx);
   }

//   ***********************< Solve >***************************
//   * Purpose: Solving a linear system for a vector           *
//   * Description: The matrix A is factorised                 *
//   *              The solution replaces x                    *
//   * Example: Solve(&A,b,&x);                                *
//   ***********************************************************
void Solve(Factored *A,const Vector &b,Vector *x)
   {
   if(b.Size() != A->numRows)
      Error("%s%s%sOnlySolve",
      Factored::ERROR,ERR_CHECK_DIMENSION,ERR_FUNCTION);
   A->PrepOnlySolve();
   if(b.WhoAmI() == x->WhoAmI())
      A->Solution(x);
   else
      A->Solution(b,x);
   }

//   ***********************< Solve >***************************
//   * Purpose: Solving a linear system for a matrix           *
//   * Description: The matrix A is factorised                 *
//   *              The solution replaces BX                   *
//   * Example: Solve(&A,&BX);                                 *
//   ***********************************************************
void Solve(Factored *A,Matrix *BX)
   {
   if(BX->Rows() != A->numRows)
      Error("%s%s%sOnlySolve",
      Factored::ERROR,ERR_CHECK_DIMENSION,ERR_FUNCTION);
   A->PrepOnlySolve();
   A->Solution(BX);
   }

//   ***********************< Solve >***************************
//   * Purpose: Solving a linear system for a matrix           *
//   * Description: The matrix A is factorised                 *
//   *              The solution replaces X                    *
//   * Example: Solve(&A,B,&X);                                *
//   ***********************************************************
void Solve(Factored *A,const Matrix &B,Matrix *X)
   {
   if(B.Rows() != A->numRows)
      Error("%s%s%sOnlySolve",
      Factored::ERROR,ERR_CHECK_DIMENSION,ERR_FUNCTION);
   A->PrepOnlySolve();
   if(B.WhoAmI() == X->WhoAmI())
      A->Solution(X);
   else
      A->Solution(B,X);
   }

//   ***********************< Solve >***************************
//   * Purpose: Solving a linear system for a vector           *
//   * Description: The matrix A is factorised                 *
//   *              The solution replaces bx                   *
//   * Example: Solve(A,&bx);                                  *
//   ***********************************************************
void Solve(Factored &A,Vector *bx)
   {
   if(bx->Size() != A.numRows)
      Error("%s%s%sSolve",
      Factored::ERROR,ERR_CHECK_DIMENSION,ERR_FUNCTION);
   A.PrepSolve();
   A.Solution(bx);
   }

//   ***********************< Solve >***************************
//   * Purpose: Solving a linear system for a vector           *
//   * Description: The matrix A is saved                      *
//   *              The solution replaces x                    *
//   * Example: Solve(A,b,&x);                                 *
//   ***********************************************************
void Solve(Factored &A,const Vector &b,Vector *x)
   {
   if(b.Size() != A.numRows)
      Error("%s%s%sSolve",
      Factored::ERROR,ERR_CHECK_DIMENSION,ERR_FUNCTION);
   A.PrepSolve();
   if(b.WhoAmI() == x->WhoAmI())
      A.Solution(x);
   else
      A.Solution(b,x);
   }

//   ***********************< Solve >***************************
//   * Purpose: Solving a linear system for a matrix           *
//   * Description: The matrix A is saved                      *
//   *              The solution replaces BX                   *
//   * Example: Solve(A,&BX);                                  *
//   ***********************************************************
void Solve(Factored &A,Matrix *BX)
   {
   if(BX->Rows() != A.numRows)
      Error("%s%s%sSolve",
      Factored::ERROR,ERR_CHECK_DIMENSION,ERR_FUNCTION);
   A.PrepSolve();
   A.Solution(BX);
   }

//   ***********************< Solve >***************************
//   * Purpose: Solving a linear system for a matrix           *
//   * Description: The matrix A is saved                      *
//   *              The solution replaces X                    *
//   * Example: Solve(A,,B,&X);                                *
//   ***********************************************************
void Solve(Factored &A,const Matrix &B,Matrix *X)
   {
   if(B.Rows() != A.numRows)
      Error("%s%s%sSolve",
      Factored::ERROR,ERR_CHECK_DIMENSION,ERR_FUNCTION);
   A.PrepSolve();
   if(B.WhoAmI() == X->WhoAmI())
      A.Solution(X);
   else
      A.Solution(B,X);
   }

//   **************< TransposeOnlySolve >***********************
//   * Purpose: Solving a linear system for a vector           *
//   *          with A transposed                              *
//   * Description: The matrix A is factorised                 *
//   *              The solution replaces bx                   *
//   * Example: TransposeSolve(&A,&bx);                        *
//   ***********************************************************
void TransposeSolve(Factored *A,Vector *bx)
   {
   if(bx->Size() != A->numColumns)
      Error("%s%s%sTransposeOnlySolve",
      Factored::ERROR,ERR_CHECK_DIMENSION,ERR_FUNCTION);
   A->PrepOnlySolve();
   A->TransposeSolution(bx);
   }

//   **************< TransposeOnlySolve >***********************
//   * Purpose: Solving a linear system for a vector           *
//   *        with A transposed                                *
//   * Description: The matrix A is factorised                 *
//   *              The solution replaces x                    *
//   * Example: TransposeSolve(&A,b,&x);                       *
//   ***********************************************************
void TransposeSolve(Factored *A,const Vector &b,Vector *x)
   {
   if(b.Size() != A->numColumns)
      Error("%s%s%sTransposeOnlySolve",
      Factored::ERROR,ERR_CHECK_DIMENSION,ERR_FUNCTION);
   A->PrepOnlySolve();
   if(b.WhoAmI() == x->WhoAmI())
      A->TransposeSolution(x);
   else
      A->TransposeSolution(b,x);
   }

//   ******************< TransposeSolve >***********************
//   * Purpose: Solving a linear system for a vector           *
//   *          with A transposed                              *
//   * Description: The matrix A is saved                      *
//   *              the solution replaces bx                   *
//   * Example: TransposeSolve(A,&bx);                         *
//   ***********************************************************
void TransposeSolve(Factored &A,Vector *bx)
   {
   if(bx->Size() != A.numColumns)
      Error("%s%s%sTransposeSolve",
      Factored::ERROR,ERR_CHECK_DIMENSION,ERR_FUNCTION);
   A.PrepSolve();
   A.TransposeSolution(bx);
   }

//   ******************< TransposeSolve >***********************
//   * Purpose: Solving a linear system for a vector           *
//   *          with A transposed                              *
//   * Description: The matrix A is saved                      *
//   *              The solution replaces x                    *
//   * Example: TransposeSolve(A,b,&x);                        *
//   ***********************************************************
void TransposeSolve(Factored &A,const Vector &b,Vector *x)
   {
   if(b.Size() != A.numColumns)
      Error("%s%s%sTransposeSolve",
      Factored::ERROR,ERR_CHECK_DIMENSION,ERR_FUNCTION);
   A.PrepSolve();
   if(b.WhoAmI() == x->WhoAmI())
      A.TransposeSolution(x);
   else
      A.TransposeSolution(b,x);
   }

//   *********************< HyperSolve >************************
//   * Purpose: Solving a linear system for a vector           *
//   *          iteratively to improve the solution            *
//   * Description: The matrix A is saved                      *
//   *              The solution replaces bx                   *
//   * Example: HyperSolveSolve(A,&bx);                        *
//   ***********************************************************
void HyperSolve(Factored &A,Vector *bx)
   {
   const int NITER = 10;
   float **mat = A.matrix;
   Vector b = *bx;
   Vector x(A.numColumns);
   Vector r(A.numRows);
   Solve(A,bx);
   x = *bx;
   double delx;
   float normdxOld;
   for(int iter =1;iter <= NITER;iter++)
      {
      for(int i=1;i<=A.numRows;i++)
        {
        delx = b[i];
        for(int j = 1;j <= A.numColumns;j++)
           delx -= mat[i][j] * (*bx)[j];
        r[i] = Control(delx);
        }
      *bx = r;
      Solve(A,bx);
      float normdx = bx->NormI();
      float normx = x.NormI();
      if(iter == 1)
        {
        normdxOld = normx; 
        if(normdx > .5*normx)break;
        }
      else if(normdx >= normdxOld)break;
      normdxOld = normdx;
      x += *bx;
      *bx = x;
      if(normdx < 10.*MACH_EPS*normx)break;
      }
   }

//   *****************< ConditionNumber >***********************
//   * Purpose: Calculating the condition number               *
//   * Example: Float condition = A.ConditionNumber();         *
//   ***********************************************************
float Factored::ConditionNumber(void)
   {
   if(numRows != numColumns)
      Error("%s%s%sConditionNumber",
      Factored::ERROR,ERR_CHECK_DIMENSION,ERR_FUNCTION);
   PrepSolve();
   if(matrixStatus == AVAILABLE)
      Norm();
   return Condition();
   }

//   *******************< Determinant >*************************
//   * Purpose: Calculating the determinant of a matrix        *
//   * Example: Double determinant = A.Determinant();          *
//   ***********************************************************
double Factored::Determinant(void)
   {
   if(numRows != numColumns)
      Error("%s%s%sDeterminant",
      Factored::ERROR,ERR_CHECK_DIMENSION,ERR_FUNCTION);
   PrepSolve();
   return DeterminantEvaluation();
   }

//   **********************< Singular >*************************
//   * Purpose: Checking if the resulting matrix is singular   *
//   * Description: If 1 the matrix is singular, if 0 it is not*
//   * Example: char singular = A.Singular();                  *
//   ***********************************************************
char Factored::Singular(void)
   {
   PrepSolve();
   return singular;
   }

//   **********************< Inverse >**************************
//   * Purpose: Calculating the inverse of a matrix            *
//   * Example: Matrix inv;                                    *
//   *          A.Inverse(&inv);                               *
//   ***********************************************************
void Factored::Inverse(Matrix *inv)
   {
   if(numRows != numColumns)
      Error("%s%s%sDeterminant",
      Factored::ERROR,ERR_CHECK_DIMENSION,ERR_FUNCTION);
   PrepSolve();
   ChangeDimensions(numRows,numColumns,inv);
   for(int i = 1;i <=numRows;i++)
      {
      for(int j =1;j <= numColumns;j++)
        (*inv)[i][j] = 0.;
      (*inv)[i][i] = 1.;
      }
   Solve(*this,inv);
   }

//   **********************< Inverse >**************************
//   * Purpose: Calculating the inverse of a matrix            *
//   * Example: Matrix inv = A.Inverse();                      *
//   ***********************************************************
Matrix Factored::Inverse(void)
   {
   if(numRows != numColumns)
      Error("%s%s%sDeterminant",
      Factored::ERROR,ERR_CHECK_DIMENSION,ERR_FUNCTION);
   PrepSolve();
   Matrix inv(numRows,numColumns);
   for(int i = 1;i <=numRows;i++)
      inv[i][i] = 1.;
   Solve(*this,&inv);
   return inv;
   }