//   =====================< MATRIX.CPP >========================
//   * Matrix Class:                                           *
//   * for operations between matrices and vectors             *
//   * Description: Chapter 12                                 *
//   * 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 <stdarg.h>
#include <math.h>

#include <alloc.h>

#include "utility.hpp"
#include "vector.hpp"
#include "matrix.hpp"
#include "right.hpp"
#include "left.hpp"
#include "symm.hpp"

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

//   ===========================================================
//   ================= Private functions =======================
//   ===========================================================

//   **********************< Initialize >***********************
//   * Purpose: Common function for initialisation             *
//   * Description: Chapter 12                                 *
//   * Example: Initialize(rows,columns)                       *
//   ***********************************************************
void Matrix::Initialize(int rows,int columns)
   {
   count++;
   whoAmI = count;
   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",
      Matrix::ERROR,ERR_SPACE,ERR_FUNCTION);
   matrix[0]=new float[size]; 
   if(!matrix[0])
      Error("%s%s%sInitialize",
      Matrix::ERROR,ERR_SPACE,ERR_FUNCTION);
   matrix[1] = matrix[0];
   for(int i = 2;i <= numRows;i++)
      matrix[i] = matrix[i-1]+numColumns;
   }

//   **********************< Deinitialize >*********************
//   * Purpose: Common function for deinitialisation           *
//   * Description: Chapter 12                                 *
//   * Example: Deinitialize();                                *
//   ***********************************************************
void Matrix::Deinitialize(void)
   {
   if(matrix == 0) return;
   delete matrix[0];
   delete matrix;
   }

//   ***********************< PrepCopy >************************
//   * Purpose: Preparing an assignment                        *
//   ***********************************************************
void Matrix::PrepCopy(int rRows,int rColumns)
   {
   int who = whoAmI;
   if(numRows != rRows || numColumns != rColumns)
      {
      Deinitialize();
      Initialize(rRows,rColumns);
      count--;
      }
   whoAmI = who;
   }

//   *********************< CopyRight >*************************
//   * Purpose: Copying a MatrixRight                          *
//   ***********************************************************
void Matrix::CopyRight(const MatrixRight &rval)
   {
   if(numRows != 0)
      {
      for(int row = 1;row <= numRows;row++)
         {
         for(int col = 1;col <= numColumns;col++)
            {
            if(col < row) 
               matrix[row][col]=0.;
            else  
               matrix[row][col] = rval.matrix[row][col];
            }
         }
      }

   }

//   *********************< CopyLeft >**************************
//   * Purpose: Copying a MatrixLeft                           *
//   ***********************************************************
void Matrix::CopyLeft(const MatrixLeft &rval)
   {
   if(numRows != 0)
      {
      for(int row = 1;row <= numRows;row++)
         {
         for(int col = 1;col <= numColumns;col++)
            {
            if(col > row) 
               matrix[row][col] = 0.;
            else  
               matrix[row][col] = rval.matrix[row][col];
            }
         }
      }
   }

//   *********************< CopySymm >**************************
//   * Purpose: Copying a MatrixSymm                           *
//   ***********************************************************
void Matrix::CopySymm(const MatrixSymm &rval)
   {
   if(numRows != 0)
      {
      for(int row = 1;row <= numRows;row++)
         {
         for(int col = 1;col <= numColumns;col++)
            {
            if(col > row)
               matrix[row][col] = rval.matrix[col][row];
            else  
               matrix[row][col] = rval.matrix[row][col];
            }
         }
      }
   }

//   ****************< private constructor >********************
//   * Purpose: Set size of  Matrix without initialisation     *
//   *          Serving internally                             *
//   * Example: Matrix A('*',4,5);                             *
//   ***********************************************************
Matrix::Matrix(char,int rows,int columns)
   {
   Initialize(rows,columns);
   }

//   ===========================================================
//   ================== Public functions =======================
//   ===========================================================

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

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

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

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

//   **************< sizing and initialising >******************
//   * Purpose: Initialising a Matrix                          *
//   * Example: Matrix A(2,3,1.,2.,3.,4.,5.,6.);               *
//   ***********************************************************
Matrix::Matrix(int rows,int columns,double a11,...)
   {
   Initialize(rows,columns);
   float *w=matrix[0] + 1;
   va_list puntaList;
   va_start(puntaList,a11);
   int i;
   *w = Control(a11);
   for(i = 2;i < size;i++)
      *++w = Control(va_arg(puntaList,double));
   va_end(puntaList);
   }

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

//   ********************< ASCII FILE >*************************
//   * Purpose: Initialising a Matrix from ASCII FILE          *
//   * Example: Matrix A("MAT.DAT");                           *
//   ***********************************************************
Matrix::Matrix(char *filematrix)
   {
   if((bzzFileIn=fopen(filematrix,"r"))==NULL)
       Error("%s%s%s%sFILE",
       Matrix::ERROR,ERR_OPEN_FILE,filematrix,ERR_CONSTRUCTOR);
   int rows,columns;
   fscanf(bzzFileIn,"%d %d",&rows,&columns);
   if(rows < 1 || columns < 1)
       Error("%s%s%sFILE",
       Matrix::ERROR,ERR_CHECK_DIMENSION,ERR_CONSTRUCTOR);
   Initialize(rows,columns);
   float *w = matrix[0];
   for(int i = 1;i < size;i++)fscanf(bzzFileIn,"%f",++w);
   fclose(bzzFileIn);
   }

//   ********************< from binary FILE >*******************
//   * Purpose: Initialising a Matrix from binary FILE         *
//   *          obtained with Save('*',file name)              *
//   * Example: Matrix A('*',"MAT.BIN");                       *
//   ***********************************************************
Matrix::Matrix(char,char *filematrix)
   {
   if((bzzFileIn=fopen(filematrix,"rb"))==NULL)
       Error("%s%s%s%sFILE",
       Matrix::ERROR,ERR_OPEN_FILE,filematrix,ERR_CONSTRUCTOR);
   int rows,columns;
   if(fread(&rows,sizeof(int),1,bzzFileIn) != 1)
      Error("%s%s%sFILE",
      Matrix::ERROR,ERR_READING_FILE,ERR_CONSTRUCTOR);
   if(fread(&columns,sizeof(int),1,bzzFileIn) != 1)
      Error("%s%s%sFILE",
      Matrix::ERROR,ERR_READING_FILE,ERR_CONSTRUCTOR);
   if(rows < 1 || columns < 1)
       Error("%s%s%sFILE",
       Matrix::ERROR,ERR_CHECK_DIMENSION,ERR_CONSTRUCTOR);
   Initialize(rows,columns);
   if(fread(matrix[1],sizeof(float)*size,1,bzzFileIn) != 1)
      Error("%s%s%sFILE",
      Matrix::ERROR,ERR_READING_FILE,ERR_CONSTRUCTOR);
   fclose(bzzFileIn);
   }

//   ********************<  from submatrix >********************
//   * Purpose: From a submatrix with identical beginning      *
//   * Example: Matrix B(4,5,A);                               *
//   ***********************************************************
Matrix::Matrix(int rows,int columns,const Matrix &rval)
   {
   if(rows < 1 || rows > rval.numRows)
       Error("%s%s%sSubMatrix",
       Matrix::ERROR,ERR_CHECK_DIMENSION,ERR_CONSTRUCTOR);
   if(columns < 1 || columns > rval.numColumns)
       Error("%s%s%sSubMatrix",
       Matrix::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      *
//   ***********************************************************
Matrix::Matrix(int rows,int columns,int irow,int jcol,
        const Matrix &rval)
   {
   if(rows < 1 || irow < 1
               ||
     irow>rval.numRows
               ||
     rows>(rval.numRows-irow+1))
        Error("%s%s%sSubMatrix",
        Matrix::ERROR,ERR_RANGE,ERR_CONSTRUCTOR);
   if(columns < 1 || jcol < 1 || 
      jcol>rval.numColumns ||
      columns>(rval.numColumns-jcol+1))
         Error("%s%s%sSubMatrix",
         Matrix::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));
   }

//   *********************< from Vector >***********************
//   * Purpose: Vector constructor                             *
//   * Example: Matrix A(y);                                   *
//   *  else:   Matrix A = y;                                  *
//   ***********************************************************
Matrix::Matrix(const Vector &rval)
   {
   Initialize(rval.dimensions,1);
   if(numRows != 0)
      memcpy(matrix[0],rval.vector,size*sizeof(float));
   }

//   **********************< from Right >***********************
//   * Purpose: From MatrixRight                               *
//   * Example: Matrix A(R);                                   *
//   *  else:   Matrix A = R;                                  *
//   ***********************************************************
Matrix::Matrix(const MatrixRight &rval)
   {
   Initialize(rval.numRows,rval.numColumns);
   CopyRight(rval);
   }

//   ************************< from Left >**********************
//   * Purpose: From MatrixLeft                                *
//   * Example: Matrix A(L);                                   *
//   *  else:   Matrix A = L;                                  *
//   ***********************************************************
Matrix::Matrix(const MatrixLeft &rval)
   {
   Initialize(rval.numRows,rval.numColumns);
   CopyLeft(rval);
   }


//   ************************< from Symm >**********************
//   * Purpose: From MatrixSymm                                *
//   * Example: Matrix A(S);                                   *
//   *  else:   Matrix A = S;                                  *
//   ***********************************************************
Matrix::Matrix(const MatrixSymm &rval)
   {
   Initialize(rval.numRows,rval.numColumns);
   CopySymm(rval);
   }

//   ===========================================================
//   ********************   destructor   ***********************
//   ===========================================================
Matrix::~Matrix(void)
   {
   Deinitialize();
   }

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

//   *********************< GetValue >**************************
//   * Purpose: Receiving values with control                  *
//   * Description: Must be used when the Matrix               *
//   *              is defined const                           *
//   * Example: float x = A.GetValue(3,5);  // A const         *
//   ***********************************************************
float Matrix::GetValue(int row,int col) const
   {
   if( (row < 1 || row > numRows)
                ||
       (col < 1 || col > numColumns) )
            Error("%s%s%sGetValue",
            Matrix::ERROR,ERR_RANGE,ERR_FUNCTION);
   return matrix[row][col];
   }

//   **********************< GetRow >***************************
//   * Purpose: Returning a row of a Matrix                    *
//   * Description: Places the row in a Vector                 *
//   * Example: Vector x = A.GetRow(3);                        *
//   ***********************************************************
Vector Matrix::GetRow(int i) const
   {
   if(i < 1 || i > numRows)
      Error("%s%s%sGetRow",
      Matrix::ERROR,ERR_RANGE,ERR_FUNCTION);
   Vector result('*',numColumns);
   memcpy(result.vector+1,matrix[i]+1,
          numColumns*sizeof(float));
   return result;
   }

//   *********************< GetColumn >*************************
//   * Purpose: Returning a column of a Matrix                 *
//   * Description: Places the column in a Vector              *
//   * Example: Vector v = GetColumn(j);                       *
//   ***********************************************************
Vector Matrix::GetColumn(int j) const
   {
   if(j < 1 || j > numColumns)
      Error("%s%s%sGetColumn",
      Matrix::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    *
//   * Description: As with A[i][j] but with control           *
//   * Example: x = A(1,5); A(3,7) = 5.;                       *
//   ***********************************************************
float &Matrix::operator () (int row,int col)
   {
   if( (row < 1 || row > numRows)
                ||
     (col < 1 || col > numColumns) )
        Error("%s%s%s()",
        Matrix::ERROR,ERR_RANGE,ERR_OPERATOR);
   return matrix[row][col];
   }

//   *********************< SetValue >**************************
//   * Purpose: Assigning values with control                  *
//   * Example:  A.SetValue(3,4,5.);                           *
//   ***********************************************************
void Matrix::SetValue(int row,int col,float val)
   {
   if( (row < 1 || row > numRows)
                ||
       (col < 1 || col > numColumns) )
            Error("%s%s%sSetValue",
            Matrix::ERROR,ERR_RANGE,ERR_FUNCTION);
   matrix[row][col] = val;
   }

//   **********************< SetRow >***************************
//   * Purpose: Setting a row of a matrix equal to a Vector    *
//   * Description: Uses memcpy                                *
//   * Example: A.SetRow(3,v);                                 *
//   ***********************************************************
void Matrix::SetRow(int i,const Vector &rval)
   {
   if( (i < 1 || i > numRows)
              ||
      (rval.dimensions != numColumns ) )
           Error("%s%s%sSetRow",
           Matrix::ERROR,ERR_RANGE,ERR_FUNCTION);
   memcpy(matrix[i]+1,rval.vector+1,
          numColumns*sizeof(float));
   }

//   *********************< SetColumn >*************************
//   * Purpose: Setting a column of a matrix equal to a Vector *
//   * Example: A.SetColumn(2,v);                              *
//   ***********************************************************
void Matrix::SetColumn(int j,const Vector &rval)
   {
   if( (j < 1 || j > numColumns)
              ||
      (rval.dimensions != numRows) )
           Error("%s%s%sSetColumn",
           Matrix::ERROR,ERR_RANGE,ERR_FUNCTION);
   float *r = rval.vector;
   for(int k = 1;k <= numRows;k++)
       matrix[k][j] = *++r;
   }

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

//   *************************< = >*****************************
//   * Purpose: Assignment of a Matrix                         *
//   * Description: Uses memcpy                                *
//   * Example: A=B;      // simple                            *
//   *          A=C+D;    // with operations                   *
//   *          A=B=C+D;  // multiple                          *
//   ***********************************************************
Matrix &Matrix::operator =(const Matrix &rval)
   {
   PrepCopy(rval.numRows,rval.numColumns);
   if(numRows != 0)
      memcpy(matrix[0],rval.matrix[0],size*sizeof(float));
   return *this;
   }

//   *************************< = >*****************************
//   * Purpose: Assignment of a Matrix from Vector             *
//   * Description: Uses memcpy                                *
//   * Example: A = v;         // simple                       *
//   *          A = 2.*v;      //  with operations             *
//   *          A = B = 2.*v;  // multiple                     *
//   ***********************************************************
Matrix &Matrix::operator =(const Vector &rval)
   {
   PrepCopy(rval.dimensions,1);
   if(numRows != 0)
      memcpy(matrix[0],rval.vector,size*sizeof(float));
   return *this;
   }
//   *************************< = >*****************************
//   * Purpose: Assignment of a Matrix from MatrixRight        *
//   * Example: A = R;          // simple                      *
//   *          A = 2.*R;       // with operations             *
//   *          A = B = 2.*R;   // multiple                    *
//   ***********************************************************
Matrix &Matrix::operator =(const MatrixRight &rval)
   {
   PrepCopy(rval.numRows,rval.numColumns);
   CopyRight(rval);
   return *this;
   }

//   *************************< = >*****************************
//   * Purpose: Assignment of a Matrix from MatrixLeft         *
//   * Example: A = L;         // simple                       *
//   *          A = 2.*L;      // with operations              *
//   *          A = B = 2.*L;  // multiple                     *
//   ***********************************************************
Matrix &Matrix::operator =(const MatrixLeft &rval)
   {
   PrepCopy(rval.numRows,rval.numColumns);
   CopyLeft(rval);
   return *this;
   }

//   *************************< = >*****************************
//   * Purpose: Assignment of a Matrix from MatrixSymm         *
//   * Example: A = S;          // simple                      *
//   *          A = 2.*S;       // with operations             *
//   *          A = B = 2.*S;   // multiple                    *
//   ***********************************************************
Matrix &Matrix::operator =(const MatrixSymm &rval)
   {
   PrepCopy(rval.numRows,rval.numColumns);
   CopySymm(rval);
   return *this;
   }

//   ===========================================================
//   **********   operators for composing matrices   ***********
//   ===========================================================

//   ************************< && >*****************************
//   * Purpose: Constructs one Matrix from two Matrices        *
//   *          inserting the second matrix beneath the first  *
//   *          They must have the same number of columns      *
//   * Description: Uses memcpy                                *
//   * Example: A = B&&C;                                      *
//   ***********************************************************
Matrix operator && (const Matrix &lval,const Matrix &rval)
   {
   if(lval.numColumns != rval.numColumns)
       Error("%s%s%s&&",
       Matrix::ERROR,ERR_CHECK_DIMENSION,ERR_OPERATOR);
   Matrix result('*',lval.numRows+rval.numRows,
                 lval.numColumns);
   if(lval.numRows != 0)
      {
      memcpy(result.matrix[0],lval.matrix[0],
             lval.size*sizeof(float));
      float *r = &result[lval.numRows][lval.numColumns] + 1;
      memcpy(r,rval.matrix[0]+1,
             (rval.size-1)*sizeof(float));
      }
   return result;
   }

//   ************************< || >*****************************
//   * Purpose: Constructs one Matrix from two Matrices        *
//   *          inserting the second matrix to the side        *
//   *          of the first                                   *
//   *          They must have the same number of rows         *
//   * Description: Uses memcpy                                *
//   * Example: A = B||C;                                      *
//   ***********************************************************
Matrix operator || (const Matrix &lval,const Matrix &rval)
   {
   if(lval.numRows != rval.numRows)
       Error("%s%s%s||",
       Matrix::ERROR,ERR_CHECK_DIMENSION,ERR_OPERATOR);
   Matrix result('*',lval.numRows,
                 lval.numColumns+rval.numColumns);
   float **s=lval.matrix;
   float **d=rval.matrix;
   float **r=result.matrix;
   if(result.numRows != 0)
      {
      for(int i=1;i<=lval.numRows;i++)
         {
         memcpy(r[i]+1,s[i]+1,
                (lval.numColumns)*sizeof(float));
         memcpy(r[i]+lval.numColumns+1,d[i]+1,
                (rval.numColumns)*sizeof(float));
         }
      }
   return result;
   }

//   ===========================================================
//   ******************   test operators ***********************
//   ===========================================================

//   ************************< == >*****************************
//   * Purpose: If it is == return 1 otherwise 0               *
//   * Example: if( A == B )                                   *
//   ***********************************************************
char operator == (const Matrix &lval,const Matrix &rval)
   {
   if(lval.whoAmI == rval.whoAmI)return 1;
   char ch = 1;
   if(   lval.numRows != rval.numRows
                      ||
      lval.numColumns != rval.numColumns) ch = 0;
   else
      {
      if(memcmp(lval.matrix[0]+1,rval.matrix[0]+1,
           (rval.size-1)*sizeof(float)) == 0) ch = 1;
      else ch = 0;
      }
   return ch;
   }

//   ************************< != >*****************************
//   * Purpose: If it is != return 1 otherwise 0               *
//   * Example: if( A != B)                                    *
//   ***********************************************************

char operator != (const Matrix &lval,const Matrix &rval)
   {
   if(lval.whoAmI == rval.whoAmI) return 0;
   char ch=0;
   if(lval.numRows != rval.numRows
             ||
      lval.numColumns != rval.numColumns) ch = 1;
   else
      ch = memcmp(lval.matrix[0]+1,rval.matrix[0]+1,
                (rval.size-1)*sizeof(float));
   return ch;
   }

//   ===========================================================
//   ====================   OPERATIONS   =======================
//   ===========================================================

//   ===========================================================
//   ***********************   Sum   ***************************
//   ===========================================================

//   ************************< Sum >****************************
//   * Purpose: Addition between matrices                      *
//   * Description: Uses Sum from utility.cpp                  *
//   * Example: Sum(A,B,&C); C = A + B;                        *
//   ***********************************************************
void Sum(const Matrix &lval,const Matrix &rval,
         Matrix *result)
   {
   if(lval.numColumns != rval.numColumns
                      ||
         lval.numRows != rval.numRows)
       Error("%s%s%sSum",
       Matrix::ERROR,ERR_CHECK_DIMENSION,ERR_FUNCTION);
   if(result->whoAmI == lval.whoAmI)
      (*result) += rval;
   else if(result->whoAmI == rval.whoAmI)
      Sum(lval,result);
   else
      {
      ChangeDimensions(lval.numRows,lval.numColumns,result,1);
      Sum(rval.size-1,lval.matrix[1]+1,rval.matrix[1]+1,
          result->matrix[1]+1);
      }
   }

//   *************************< + >*****************************
//   * Purpose: Addition of matrices                           *
//   * Example: A = B + C;                                     *
//   ***********************************************************
Matrix operator + (const Matrix &lval,const Matrix &rval)
   {
   Matrix result('*',lval.numRows,lval.numColumns);
   Sum(lval,rval,&result);
   return result;
   }

//   ************************< Sum >****************************
//   * Purpose: Addition of matrices                           *
//   * Description: The result goes in the left matrix         *
//   * Example: Sum(&A,B); A += B;                             *
//   ***********************************************************
void Sum(Matrix *lvalAndResult,const Matrix &rval)
   {
   (*lvalAndResult) += rval;
   }

//   *************************< += >****************************
//   * Purpose: Adds a matrix which is also on the right       *
//   * Description: Uses Sum from utility.cpp                  *
//   * Example: A += B; // A = A + B;                          *
//   ***********************************************************
Matrix &Matrix::operator += (const Matrix &rval)
   {
   if(numColumns != rval.numColumns ||
      numRows != rval.numRows)
         Error("%s%s%s+=",
         Matrix::ERROR,ERR_CHECK_DIMENSION,ERR_OPERATOR);
   if(whoAmI == rval.whoAmI)
      Sum(this);
   else
      Sum(rval.size-1,matrix[1]+1,rval.matrix[1]+1);
   return *this;
   }

//   ************************< Sum >****************************
//   * Purpose: Addition between matrices                      *
//   * Description: The result goes in the right matrix        *
//   * Example: Sum(B,&A); A = B + A;                          *
//   ***********************************************************
void Sum(const Matrix &lval,Matrix *rvalAndResult)
   {
   (*rvalAndResult) += lval;
   }

//   ************************< Sum >****************************
//   * Purpose: Addition of two identical matrices             *
//   * Description: The result replaces the matrix             *
//   * Example: Sum(&A); A = A + A;                            *
//   ***********************************************************
void Sum(Matrix *lvalRvalAndResult)
   {
   Sum(lvalRvalAndResult->size-1,
       lvalRvalAndResult->matrix[1]+1);
   }

//   ===========================================================
//   ********************   Difference   ***********************
//   ===========================================================

//   **********************< Difference >***********************
//   * Purpose: Difference between matrices                    *
//   * Description: Uses Difference from utility.cpp           *
//   * Example: Difference(B,C,&A); A = B - C;                 *
//   ***********************************************************
void Difference(const Matrix &lval,const Matrix &rval,
                Matrix *result)
   {
   if(lval.numColumns != rval.numColumns
                      ||
         lval.numRows != rval.numRows)
       Error("%s%s%sSum",
       Matrix::ERROR,ERR_CHECK_DIMENSION,ERR_FUNCTION);
   if(result->whoAmI == lval.whoAmI)
      (*result) -= rval;
   else if(result->whoAmI == rval.whoAmI)
      Difference(lval,result);
   else
      {
      ChangeDimensions(lval.numRows,lval.numColumns,result,1);
      Difference(rval.size-1,lval.matrix[1]+1,
                 rval.matrix[1]+1,result->matrix[1]+1);
      }
   }

//   *************************< - >*****************************
//   * Purpose: Difference between matrices                    *
//   * Description: Uses Difference                            *
//   * Example: A = B - C;                                     *
//   ***********************************************************
Matrix operator - (const Matrix &lval,const Matrix &rval)
   {
   Matrix result('*',lval.numRows,lval.numColumns);
   Difference(lval,rval,&result);
   return result;
   }

//   **********************< Difference >***********************
//   * Purpose: Difference between matrices                    *
//   * Description: The result goes in the left matrix         *
//   * Example: Difference(&A,B); A = A - B;                   *
//   ***********************************************************
void Difference(Matrix *lvalAndResult,const Matrix &rval)
   {
   (*lvalAndResult) -= rval;
   }

//   *************************< -= >****************************
//   * Purpose: Difference from a matrix also on the right     *
//   * Description: Uses Sum from utility.cpp                  *
//   * Example: A -= B; // A = A - B;                          *
//   ***********************************************************
Matrix &Matrix::operator -= (const Matrix &rval)
   {
   if(numColumns != rval.numColumns ||
      numRows != rval.numRows)
         Error("%s%s%s+=",
         Matrix::ERROR,ERR_CHECK_DIMENSION,ERR_OPERATOR);
   if(whoAmI == rval.whoAmI)
      Difference(this);
   else
      Difference(rval.size-1,matrix[1]+1,rval.matrix[1]+1);
   return *this;
   }

//   **********************< Difference >***********************
//   * Purpose: Difference between matrices                    *
//   * Description: The result goes in the right matrix        *
//   * Example: Difference(A,&B); B= A - B;                    *
//   ***********************************************************
void Difference(const Matrix &lval,Matrix *rvalAndResult)
   {
   if(lval.numColumns != rvalAndResult->numColumns
                ||
      lval.numRows != rvalAndResult->numRows)
       Error("%s%s%sDiffernece",
       Matrix::ERROR,ERR_CHECK_DIMENSION,ERR_FUNCTION);
   Difference(lval.size-1,lval.matrix[1]+1,
              rvalAndResult->matrix[1]+1,1);
   }

//   **********************< Difference >***********************
//   * Purpose: Difference between identical matrices          *
//   * Description:  Null matrix                               *
//   * Example: Difference(&A); A = A - A;                     *
//   ***********************************************************
void Difference(Matrix *lvalRvalAndResult)
   {
   float *w = lvalRvalAndResult->matrix[0];
   for(int i = 1;i < lvalRvalAndResult->size;i++)
       *++w = 0.;
   }

//   ===========================================================
//   ***********************   Minus   *************************
//   ===========================================================

//   ************************< Minus >**************************
//   * Purpose: Unary minus                                    *
//   * Example: Minus(A,&B); B = -A;                           *
//   ***********************************************************
void Minus (const Matrix &rval,Matrix *result)
   {
   ChangeDimensions(rval.numRows,rval.numColumns,result,1);
   float *w=rval.matrix[1];
   float *r=result->matrix[1];
   for(int i = 1;i < result->size;i++)
   r[i] = -w[i];
   }    

//   *************************< - >*****************************
//   * Purpose: Unary minus                                    *
//   * Example: A = -B;                                        *
//   ***********************************************************
Matrix operator -(const Matrix &rval) //unary minus
   {
   Matrix result('*',rval.numRows,rval.numColumns);
   Minus(rval,&result);
   return result;
   }

//   ************************< Minus >**************************
//   * Purpose: Unary minus                                    *
//   * Example: Minus(&A); A = -A;                             *
//   ***********************************************************
void Minus (Matrix *rvalAndResult)
   {
   float *w=rvalAndResult->matrix[1];
   for(int i = 1;i < rvalAndResult->size;i++)
   w[i] = -w[i];
   }    

//   ===========================================================
//   **********************   Product   ************************
//   ===========================================================

//   ***********************< Product >*************************
//   * Purpose: Product between matrices                       *
//   * Description: Pointers are used to eliminate the         *
//                  use of double indices                      *
//   * Example: Product(A,B,&C); C = A * B;                    *
//   ***********************************************************
void Product(const Matrix &lval,const Matrix &rval,
             Matrix *result)
   {
   if(lval.numColumns != rval.numRows)
       Error("%s%s%sProduct",
       Matrix::ERROR,ERR_CHECK_DIMENSION,ERR_OPERATOR);
   if(result->whoAmI == lval.whoAmI)
      (*result) *= rval;
   else if(result->whoAmI == rval.whoAmI)
      Product(lval,result);
   else
      {
      ChangeDimensions(lval.numRows,rval.numColumns,
                       result,1);
      float *v = result->matrix[1];
      float *b;
      int nc = rval.numColumns;
      for(int row = 1;row <= lval.numRows;row++)
         {
         for(int column = 1;column <= rval.numColumns;
                 column++)
            {
            b = rval.matrix[1] + column;
            float *w = lval.matrix[row];
            double sum=0.;
            for(int i=1;i<=lval.numColumns;i++)
               {
               sum += *++w * (*b);
               b += nc;
               }
            *++v = Control(sum);
            }
         }
      }
   }
//   *************************< * >*****************************
//   * Purpose: Product between matrices                       *
//   * Description: Uses Product(A,B,&C)                       *
//   * Example: C = A * B;                                     *
//   ***********************************************************
Matrix operator * (const Matrix &lval,const Matrix &rval)
   {
   Matrix result('*',lval.numRows,rval.numColumns);
   Product(lval,rval,&result);
   return result;
   }

//   ***********************< Product >*************************
//   * Purpose: Product between matrices                       *
//   * Description: The result goes in the left matrix         *
//   * Example:  Product(&A,B); A *= B;                        *
//   ***********************************************************
void Product(Matrix *lvalAndResult,const Matrix &rval)
   {
   (*lvalAndResult) *= rval;
   }

//   ************************< *= >*****************************
//   * Purpose: The product of one matrix with another,        *
//   *          stored replacing the first                     *
//   * Description: Always favour this alternative             *
//   *               when possible!                            *
//   * Example:  A *= B; // equivalent to  A = A * B           *
//   ***********************************************************
Matrix &Matrix::operator *= (const Matrix &rval)
   {
   if(numColumns != rval.numRows)
       Error("%s%s%s*=",
       Matrix::ERROR,ERR_CHECK_DIMENSION,ERR_OPERATOR);
   if(whoAmI == rval.whoAmI)
      Product(this); // A = A*A;
   else
      {
      // if rval is not square *this changes dimensions
      if(rval.numRows == rval.numColumns)
         {
         float *b;
         int nc = rval.numColumns;
         float *v = matrix[1];
         float *w = new float[numColumns+1];
         for(int row = 1;row <= numRows;row++)
            {
            memcpy(w,matrix[row],
                    (numColumns+1)*sizeof(float));
            for(int column = 1;column <= numColumns;
                    column++)
               {
               b = rval.matrix[1] + column;
               double sum=0.;
               for(int i=1;i<=numColumns;i++)
                  {
                  sum += w[i] * (*b);
                  b += nc;
                  }
               *++v = Control(sum);
               }
            }
         delete w;
         }
      else
         {
         Matrix result;
         Product(*this,rval,&result);
         Swap(this,&result);// avoids copy
         }
      }
   return *this;
   }


//   ***********************< Product >*************************
//   * Purpose: Product between matrices                       *
//   * Description: The result goes in the right matrix        *
//   * Example: Product(B,&A); A = B * A;                      *
//   ***********************************************************
void Product(const Matrix &lval,Matrix *rvalAndResult)
   {
   Matrix C;
   Product(lval,*rvalAndResult,&C);
   Swap(rvalAndResult,&C);
   }

//   ***********************< Product >*************************
//   * Purpose: Product between identical matrices             *
//   * Description: The result replaces the matrix             *
//   * Example: Product(&A); A = A * A;                        *
//   ***********************************************************
void Product(Matrix *lvalRvalAndResult)
   {
   Matrix C;
   Product(*lvalRvalAndResult,*lvalRvalAndResult,&C);
   Swap(lvalRvalAndResult,&C);
   }


//   ***********************< Product >*************************
//   * Purpose: Product of a float with a Matrix               *
//   * Description: Uses Product from utility.cpp              *
//   * Example: A = 3.5 * B;                                   *
//   ***********************************************************
void Product(float lval,const Matrix &rval,Matrix *result)
   {
   if(rval.whoAmI == result->whoAmI)
      Product(result->size - 1,lval,result->matrix[1]+1);
   else
      {
      ChangeDimensions(rval.numRows,rval.numColumns,
                       result,1);
      Product(rval.size - 1,lval,rval.matrix[1]+1,
           result->matrix[1]+1);
      }
   }

//   *************************< * >*****************************
//   * Purpose: Product of a float with a Matrix               *
//   * Description: Uses Product from utility.cpp              *
//   * Example: A = 3.5 * B;                                   *
//   ***********************************************************
Matrix operator * (float lval,const Matrix &rval)
   {
   Matrix result('*',rval.numRows,rval.numColumns);
   Product(lval,rval,&result);
   return result;
   }

//   ***********************< Product >*************************
//   * Purpose: Product of a float with a Matrix               *
//   * Description: Uses Product from utility.cpp              *
//   * Example: Product(&A,3.); A = 3.* A;                     *
//   ***********************************************************
void Product(float lval,Matrix *rvalAndResult)
   {
   Product(rvalAndResult->size - 1,lval,
     rvalAndResult->matrix[1]+1);
   }

//   ************************< *= >*****************************
//   * Purpose: Multiplies and modifies a Matrix by a float    *
//   * Example: A *= 3.; // A =  3.*A;                         *
//   ***********************************************************
Matrix &Matrix::operator *= (float lval)
   {
   Product(size - 1,lval,matrix[1] + 1);//in UTILITY.CPP
   return *this;
   }

//   ===========================================================
//   *********************   TProduct   ************************
//   ===========================================================

//   **********************< TProduct >*************************
//   * Purpose: Product of a transposed matrix with a          *
//   *          normal matrix                                  *
//   * Example: TProduct(A,B,&C);  // C = ATB                 *
//   ***********************************************************
void TProduct(const Matrix &lval,const Matrix &rval,
              Matrix *result)
   {
   if(lval.numRows!=rval.numRows)
       Error("%s%s%s %",
       Matrix::ERROR,ERR_CHECK_DIMENSION,ERR_OPERATOR);
   if(result->whoAmI == lval.whoAmI)
      (*result) %= rval;
   else if(result->whoAmI == rval.whoAmI)
      TProduct(lval,result);
   else
      {
      ChangeDimensions(lval.numColumns,
                       rval.numColumns,result,1);
      if(lval.whoAmI != rval.whoAmI)
         {
         float *p,*b;
         int np = lval.numColumns;
         int nc = rval.numColumns;
         float *v = result->matrix[1];
         for(int cc = 1;cc <= lval.numColumns;cc++)
            {
            for(int column = 1;column <= rval.numColumns;column++)
               {
               double sum=0.;
               p = lval.matrix[1] + cc; 
               b = rval.matrix[1] + column; 
               for(int i = 1;i <= lval.numRows;i++)
                  {
                  sum += (*p) * (*b);
                  p += np;
                  b += nc;
                  }
               *++v = Control(sum);
               }
            }
         }
      else  // identical matrices ATA
         {
         int n = lval.numColumns;
         int m = lval.numRows;
         int i,j,k;
         float **a = lval.matrix;
         float **aTa = result->matrix;
         for(i = 1;i <=n;i++)
            {
            for(j = i;j <= n;j++)
               {
               double sum = 0.;
               for(k = 1;k <= m;k++)
                  sum += a[k][i]*a[k][j];        
               aTa[i][j] = Control(sum);
               }
            }
         for(i = 1;i <= n;i++)
            for(j = 1;j < i;j++) 
               aTa[i][j] = aTa[j][i];
         }
      }
   }

//   *************************< % >*****************************
//   * Purpose: Product of a transposed matrix with a          *
//   *          normal matrix                                  *
//   * Example: C = A % B; // C = ATB                          *
//   ***********************************************************
Matrix operator % (const Matrix &lval,const Matrix &rval)
   {
   Matrix result('*',lval.numColumns,rval.numColumns);
   TProduct(lval,rval,&result);
   return result;
   }

//   **********************< TProduct >*************************
//   * Purpose: Executes A = ATB                               *
//   * Description: The result goes in the left matrix         *
//   * Example: TProduct(&A,B); A %= B;                        *
//   ***********************************************************
void TProduct(Matrix *lvalAndResult,const Matrix &rval)
   {
   (*lvalAndResult) %= rval;
   }

//   ************************< %= >*****************************
//   * Purpose: Executes A = ATB                               *
//   * Example: A %= B; A %= A;                                *
//   ***********************************************************
Matrix &Matrix::operator %=(const Matrix &rval)
   {
   if(numRows != rval.numRows)
       Error("%s%s %=",
       Matrix::ERROR,ERR_CHECK_DIMENSION,ERR_OPERATOR);
   if(*this == rval) TProduct(this);  // A=ATA
   else if(numRows == numColumns &&
           rval.numRows == rval.numColumns)
      { // to optimise if square
      Matrix result;
      TProduct(*this,rval,&result);
      Swap(this,&result);
      }
   else
      { // general case
      Matrix result;
      TProduct(*this,rval,&result);
      Swap(this,&result);
      }
   return *this;
   }

//   ***********************< TProduct >************************
//   * Purpose: Calculates A = BTA                             *
//   * Example: TProduc(B,&A); A = BTA;                        *
//   ***********************************************************
void TProduct(const Matrix &lval,Matrix *rvalAndResult)
   {
   Matrix C;
   TProduct(lval,*rvalAndResult,&C);
   Swap(rvalAndResult,&C);
   }

//   **********************< TProduct >*************************
//   * Purpose: Calculates ATA and puts it in A                *
//   * Example: TProduct(&A);                                  *
//   ***********************************************************
void TProduct(Matrix *lvalRvalAndResult)
   {
   if(lvalRvalAndResult->numRows == 
         lvalRvalAndResult->numColumns)
      { // square matrices
      int i,j,k,n = lvalRvalAndResult->numRows;
      Vector s(n);
      float **aTa = lvalRvalAndResult->matrix;
      for(j = 1;j <= n;j++)
         {
         for(i = j;i <= n;i++)
            {
            double sum = 0.;
            for(k = 1;k <= n;k++) 
               sum += aTa[k][i]*aTa[k][j];
            s[i] = Control(sum);
            }
         for(i = j;i <= n;i++) aTa[i][j] = s[i];
         }
      for(j = 1;j <= n;j++)
         for(i = j+1;i <= n;i++) 
            aTa[j][i] = aTa[i][j];
      }
   else // rectangular matrices
      {
      Matrix result;
      TProduct(*lvalRvalAndResult,*lvalRvalAndResult,
               &result);
      Swap(lvalRvalAndResult,&result);
      }
   }

//   ===========================================================
//   *********************   ProductT   ************************
//   ===========================================================

//   **********************< ProductT *>************************
//   * Purpose: Executes AATand puts it in A                   *
//   * Description: Symmetrical matrix                         *
//   * Example: ProductT(&A); A = AAT;                         *
//   ***********************************************************
void ProductT (Matrix *lvalRvalAndResult)
   {
   if(lvalRvalAndResult->numRows == 
      lvalRvalAndResult->numColumns)
      { // square matrices 
      int i,j,k,n = lvalRvalAndResult->numRows;
      Vector s(n);
      float **a = lvalRvalAndResult->matrix;
      for(i = 1;i <= n;i++)
         {
         for(j = i;j <= n;j++)
            {
            double sum = 0.;
            float *ai = a[i];
            float *aj = a[j];
            for(k = 1;k <= n;k++) sum += (*++ai) * (*++aj);
            s[j] = Control(sum);
            }
         memcpy(a[i]+i,&s[0]+i,(n-i+1)*sizeof(float));
         }
      for(i = 1;i <= n;i++)
      for(j = i+1;j <= n;j++) a[j][i] = a[i][j];
      }
   else // rectangular matrices 
      {
      int n = lvalRvalAndResult->numRows;
      int m = lvalRvalAndResult->numColumns;
      int i,j,k;
      Matrix result('*',n,n);
      float **a = lvalRvalAndResult->matrix;
      float **aaT = result.matrix;
      for(i = 1;i <= n;i++)
         {
         for(j = i;j <= n;j++)
            {
            double sum = 0.;
            float *ai = a[i];
            float *aj = a[j];
            for(k = 1;k <= m;k++) sum += (*++ai) * (*++aj);
            aaT[i][j] = Control(sum);
            }
         }
      for(i = 1;i <= n;i++)
         for(j = 1;j < i;j++) 
            aaT[i][j] = aaT[j][i];
      Swap(lvalRvalAndResult,&result);
      }
   }

//   ===========================================================
//   **********************   Division   ***********************
//   ===========================================================

//   ***********************< Division >************************
//   * Purpose: Dividing the coefficients of a matrix          *
//   *          by a float                                     *
//   * Description: Uses Division from utility.cpp             *
//   * Example: Division(A,3.,&B); B = A/3.;                   *
//   ***********************************************************
void Division(const Matrix &lval,float rval,Matrix *result)
   {
   if(lval.whoAmI == result->whoAmI)
      Division(result->size - 1,result->matrix[1]+1,rval);
   else
      {
      ChangeDimensions(lval.numRows,lval.numColumns,
                       result,1);
      Division(lval.size - 1,lval.matrix[1]+1,
               rval,result->matrix[1]+1);
      }
   }

//   *************************< / >*****************************
//   * Purpose: Dividing the coefficients of a matrix          *
//   *          by a float                                     *
//   * Description: Uses Division from utility.cpp             *
//   * Example: A = B/3.;                                      *
//   ***********************************************************
Matrix operator / (const Matrix &lval,float rval)
   {
   Matrix result('*',lval.numRows,lval.numColumns);
   Division(lval,rval,&result);
   return result;
   }

//   ***********************< Division >************************
//   * Purpose: Dividing the coefficients of a matrix          *
//   *          by a float                                     *
//   * Description: Uses Division from utility.cpp             *
//   * Example: Division(&A,3.); A = A/3.;                     *
//   ***********************************************************
void Division(Matrix *lvalAndResult,float rval)
   {
   Division(lvalAndResult->size - 1,
            lvalAndResult->matrix[1]+1,rval);
   }

//   ************************< /= >*****************************
//   * Purpose: Dividing and modifying a Matrix by a float     *
//   * Description: Uses Division from utility.cpp             *
//   * Example: A /= 2.; // A = A/2.;                          *
//   ***********************************************************
Matrix &Matrix::operator /= (float rval)
   {
   Division(size - 1,matrix[1] + 1,rval);//in UTILITY.CPP
   return *this;
   }

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

//   ***********************< Print >***************************
//   * Purpose: Printing out the elements of a matrix          *
//   * Description: Uses the file bzzFileOut global            *
//   *              in utility.hpp (Chapter 8)                 *
//   * Example: A.Print("Matrix A");                           *
//   ***********************************************************
void Matrix::Print(char *msg)
   {
   fprintf(bzzFileOut,"\nMatrix No.%d",whoAmI);
   if(*msg)fprintf(bzzFileOut,"  %s",msg);
   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");
      }
   }

//   ************************< Save >***************************
//   * Purpose: Saving  a Matrix on formatted file             *
//   * Description: Can be read.                               *
//   *              See TempFile in utility                    *
//   * Example: A.Save("Mat.Dat");                             *
//   ***********************************************************
void Matrix::Save(char *filematrix) // formatted
   {
   FILE *bzz;
   if((bzz = fopen(filematrix,"w")) == NULL)
       Error("%s%s%s%sSave",
       Matrix::ERROR,ERR_OPEN_FILE,filematrix,ERR_FUNCTION);

   fprintf(bzz," %d %d\n",numRows,numColumns);
   for(int row = 1;row <= numRows;row++)
      {
      for(int column = 1;column <= numColumns;column++)
         fprintf(bzz," %15.7e",matrix[row][column]);
      fprintf(bzz,"\n");
      }
   fclose(bzz);
   }

//   ************************< Save >***************************
//   * Purpose: Saving a Matrix on an unformatted file         *
//   * Description: Can be read                                *
//   * Example: A.Save('*',MAT.BIN");                          *
//   *          See also TempFile from utility.cpp             *
//   ***********************************************************
void Matrix::Save(char,char *filematrix) // unformatted
   {
   FILE *bzz;
   if((bzz = fopen(filematrix,"wb")) == NULL)
       Error("%s%s%s%sSave",
       Matrix::ERROR,ERR_OPEN_FILE,filematrix,ERR_FUNCTION);
   if(fwrite(&numRows,sizeof(int),1,bzz) != 1)
       Error("%s%s%sSave",
       Matrix::ERROR,ERR_WRITING_FILE,ERR_FUNCTION);
   if(fwrite(&numColumns,sizeof(int),1,bzz) != 1)
       Error("%s%s%sSave",
       Matrix::ERROR,ERR_WRITING_FILE,ERR_FUNCTION);
   if(fwrite(matrix[1],sizeof(float)*size,1,bzz) != 1)
       Error("%s%s%sSave",
       Matrix::ERROR,ERR_WRITING_FILE,ERR_FUNCTION);
   fclose(bzz);
   }

//   ===========================================================
//   *********************   Norms   ***************************
//   ===========================================================

//   ***********************< NormT >***************************
//   * Purpose: Total matrix NORM  (n times the max of aij)    *
//   * Example: float normT = A.NormT ();                      *
//   ***********************************************************
float Matrix::NormT(void)
   {
   float norm=0.,aus;
   int n=(numRows > numColumns) ? numRows : numColumns;
   float *w = matrix[0];
   for(int i = 1;i < size;i++)
      if((aus = Abs(*++w)) > norm) norm = aus;
   return (float)n*norm;
   }

//   ***********************< NormR >***************************
//   * Purpose: Max of the sums of the abs values              *
//   *          of coeff. rows                                 *
//   * Example: float normR = A.NormR ();                      *
//   ***********************************************************
float Matrix::NormR(void)
   {
   float norm = 0.;
   float *w = matrix[0];
   for(int i = 1;i <= numRows;i++)
       {
       double sum = 0.;
       for(int j = 1;j <= numColumns;j++)sum += Abs(*++w);
       if(sum > norm) norm = Control(sum);
       }
   return norm;
   }

//   ***********************< NormC >***************************
//   * Purpose: Norm C                                         *
//   * Example: float normC = A.NormC ();                      *
//   ***********************************************************
float Matrix::NormC(void)
   {
   float norm = 0.;
   for(int j = 1;j <= numColumns;j++)
       {
       double sum = 0.;
       for(int i = 1;i <= numRows;i++)
          sum += Abs(matrix[i][j]);
       if(sum > norm) norm = Control(sum);
       }
   return norm;
   }

//   ***********************< NormF >***************************
//   * Purpose: Frobenius' norm (the square root of the        *
//   *          sums of squares of all the coefficients)       *
//   * Description: Uses the SqrtSumSqr from utility.cpp       *
//   * Example: float normF = A.NormF ();                      *
//   ***********************************************************
float Matrix::NormF(void)
   {
   return SqrtSumSqr(size-1,matrix[0]+1);
   }

//   ===========================================================
//   *                       Max Min                           *
//   ===========================================================

//   ************************< Max >****************************
//   * Purpose: Maximum of a  Matrix                           *
//   * Description: Uses the function described in UTILITY.CPP *
//   * Example: float xmax = A.Max();                          *
//   *  else:    float xmax = A.Max(&imax,&jmax);              *
//   ***********************************************************
float Matrix::Max(int *imax,int *jmax)
   {
   int iaux,im,jm;
   float xm;
   xm = ::Max(size-1,matrix[0] + 1,&iaux);
   if(imax == 0 && jmax == 0)
      return xm;
   im = iaux/numColumns + 1;
   jm = iaux - numColumns*(im - 1) + 1;
   if(imax != 0)*imax = im;
   if(jmax != 0)*jmax = jm;
   return xm;
   }

//   **********************< MaxAbs >***************************
//   * Purpose: Absolute maximum of a Matrix                   *
//   * Description: Uses the function defined in UTILITY.CPP   *
//   * Example: float xmax = A.MaxAbs();                       *
//   *  else:   float xmax = A.MaxAbs(&imax,&jmax);            *
//   ***********************************************************
float Matrix::MaxAbs(int *imax,int *jmax)
   {
   int iaux,im,jm;
   float xm;
   xm = ::MaxAbs(size-1,matrix[0] + 1,&iaux);
   if(imax == 0 && jmax == 0)
      return xm;
   im = iaux/numColumns + 1;
   jm = iaux - numColumns*(im - 1) + 1;
   if(imax != 0)*imax = im;
   if(jmax != 0)*jmax = jm;
   return xm;
   }

//   ************************< Min >****************************
//   * Purpose: Minimum of a Matrix                            *
//   * Description: Uses the function defined in UTILITY.CPP   *
//   * Example: float xmin = A.Min();                          *
//   *  else:   float xmin = A.Min(&imin,&jmin);               *
//   ***********************************************************
float Matrix::Min(int *imin,int *jmin)
   {
   int iaux,im,jm;
   float xm;
   xm = ::Min(size-1,matrix[0] + 1,&iaux);
   if(imin == 0 && jmin == 0)
      return xm;
   im = iaux/numColumns + 1;
   jm = iaux - numColumns*(im - 1) + 1;
   if(imin != 0)*imin = im;
   if(jmin != 0)*jmin = jm;
   return xm;
   }

//   **********************< MinAbs >***************************
//   * Purpose: Absolute minimum of a Matrix                   *
//   * Description: Uses the function defined in UTILITY.CPP   *
//   * Example: float xmin = A.MinAbs();                       *
//   *  else:   float xmin = A.MinAbs(&imin,&jmin);            *
//   ***********************************************************
float Matrix::MinAbs(int *imin,int *jmin)
   {
   int iaux,im,jm;
   float xm;
   xm = ::MinAbs(size-1,matrix[0] + 1,&iaux);
   if(imin == 0 && jmin == 0)
      return xm;
   im = iaux/numColumns + 1;
   jm = iaux - numColumns*(im - 1) + 1;
   if(imin != 0)*imin = im;
   if(jmin != 0)*jmin = jm;
   return xm;
   }

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

//   **********************< Delete >***************************
//   * Purpose: Eliminating a matrix of no use                 *
//   *          without leaving the purpose                    *
//   * Example: Delete(&A);                                    *
//   ***********************************************************
void Delete(Matrix *result)
   {
   result->Deinitialize();
   result->size = result->numRows = result->numColumns = 0;
   result->matrix = 0;
   }

//   ***************< ChangeDimensions >************************
//   * Purpose: Changing the dimensions of a Matrix            *
//   * Example: ChangeDimensions(row,columns,&A);              *
//   ***********************************************************
void ChangeDimensions(int rows,int columns,
                      Matrix *result,char zero)
   {
   result->PrepCopy(rows,columns);
   if(zero == 0 && result->size != 0)
      memset(result->matrix[0],0,
             result->size*sizeof(float));
   }

//   **********************< Recover >**************************
//   * Purpose: Recovering a matrix stored with                *
//   *          formatted Save(file name)                      *
//   * Description: It is not necessary to have the same       *
//   *              dimensions as the saved matrix             *
//   * Example: Recover(&x,file name);                         *
//   ***********************************************************
void Recover(Matrix *A,char *filematrix)
   {
   if((bzzFileIn=fopen(filematrix,"r"))==NULL)
       Error("%s%s%s%sRecover",
       Matrix::ERROR,ERR_OPEN_FILE,filematrix,ERR_FUNCTION);
   int rows,columns;
   fscanf(bzzFileIn,"%d %d",&rows,&columns);
   ChangeDimensions(rows,columns,A,1);
   float *w = A->matrix[0];
   for(int i = 1;i < A->size;i++)fscanf(bzzFileIn,"%f",++w);
   fclose(bzzFileIn);
   }

//   **********************< Recover >**************************
//   * Purpose: Recovering a matrix stored by                  *
//   *          unformatted  Save('*', file name)              *
//   * Description: It is not necessary to have the same       *
//   *              dimensions as the saved matrix             *
//   * Example: Recover(&A,'*',file name);                     *
//   ***********************************************************
void Recover(Matrix *A,char,char *filematrix)
   {
   if((bzzFileIn = fopen(filematrix,"rb")) == NULL)
       Error("%s%s%s%sRecover *",
       Matrix::ERROR,ERR_OPEN_FILE,filematrix,ERR_FUNCTION);
   int rows,columns;
   if(fread(&rows,sizeof(int),1,bzzFileIn) != 1)
      Error("%s%s%sRecover",
      Matrix::ERROR,ERR_READING_FILE,ERR_FUNCTION);
   if(fread(&columns,sizeof(int),1,bzzFileIn) != 1)
      Error("%s%s%sRecover",
      Matrix::ERROR,ERR_READING_FILE,ERR_FUNCTION);
   ChangeDimensions(rows,columns,A,1);
   if(rows < 1 || columns < 1)
       {fclose(bzzFileIn);return;}
   if(fread(A->matrix[1],sizeof(float)*(A->size),
          1,bzzFileIn) != 1)
         Error("%s%s%sRecover",
         Matrix::ERROR,ERR_READING_FILE,ERR_FUNCTION);
   fclose(bzzFileIn);
   }

//   *********************< Transpose >*************************
//   * Purpose: Modifying a matrix into its transpose          *
//   * Description: If the matrix is square it replaces itself *
//   * Example: Transpose(&A);                                 *
//   ***********************************************************
void Transpose(Matrix *A)
   {
   int i,j;
   float **a=A->matrix;
   if(A->numRows == A->numColumns)
      {
      for(i = 1;i <= A->numRows;i++)
      for(j = i+1;j <= A->numColumns;j++)
         Swap(&a[i][j],&a[j][i]); // in utility.cpp
      }
   else
      {
      Matrix B('*',A->numColumns,A->numRows);
      for(i = 1;i <= A->numRows;i++)
         for(j = 1;j <= A->numColumns;j++)
            B[j][i] = a[i][j];
      Swap(A,&B); // avoid copying
      }
   }

//   *********************< SumRankOne >************************
//   * Purpose: Sum at matrix one of rank one                  *
//   * Example: SumRankOne(u,vt,&A); A = A + u*vT              *
//   ***********************************************************
void SumRankOne(const Vector &u,
          const Vector &vT,Matrix *result)
   {
   if(u.dimensions != vT.dimensions ||
      u.dimensions != result->numRows ||
      u.dimensions != result->numColumns)
         Error("%s%s%sSumRankOne",
         Matrix::ERROR,ERR_CHECK_DIMENSION,ERR_FUNCTION);
   float *r = result->matrix[1] + 1;
   float *left = u.vector + 1;
   float *right = vT.vector;
   for(int i = 1;i <= u.dimensions;i++)
      {
      for(int j = 1;j <= u.dimensions;j++)
         *r++ = Control(*r + (*left) * right[j]);
      left++;
      }
   }

//   *********************< SumRankOne >************************
//   * Purpose: Sum at matrix one of rank one                  *
//   * Example: SumRankOne(a,u,vt,&A); A = A + a*u*vT          *
//   ***********************************************************
void SumRankOne(float alfa,const Vector &u,
          const Vector &vT,Matrix *result)
   {
   Vector au;
   Product(alfa,u,&au);
   SumRankOne(au,vT,result);
   }

//   *********************< SumRankOne >************************
//   * Purpose: Sum at matrix one of rank one                  *
//   * Example: SumRankOne(u,vt,b,&A); A = A + u*vT/b          *
//   ***********************************************************
void SumRankOne(const Vector &u,const Vector &vT,
                float beta,Matrix *result)
   {
   Vector au;
   Division(u,beta,&au);
   SumRankOne(au,vT,result);
   }
 
//   ************************< Swap >***************************
//   * Purpose: Swapping the elements of any two Matrices      *
//   * Description: Provides an efficient method of copying    *
//   *              when a Matrix remains in the purpose and   *
//   *              another one leaves, they swap              *
//   * Example: Swap(&A,&B);                                   *
//   ***********************************************************
void Swap(Matrix *lval,Matrix *rval)
   {
   float **temp = lval->matrix;
   lval->matrix = rval->matrix;
   rval->matrix = temp;
   Swap(&lval->numColumns,&rval->numColumns);
   Swap(&lval->numRows,&rval->numRows);
   Swap(&lval->size,&rval->size);
   }