//   ======================< MATRIX.HPP >=======================
//   * 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)                                   *
//   ===========================================================
//   ****** Constructors by Matrix:                            *
//   * Matrix A; // default                                    *
//   * Matrix A = B; // copy-initializer                       *
//   * Matrix A(m,n); // sizes and places at 0                 *
//   * Matrix A(2,3,1.,2.,3.,4.,5.,6.); // matrix 2X3          *
//   * float x[6]={1.,2.,3.,4.,5.,6.};                         *
//   * Matrix A(2,3,x); // from array                          *
//   * Matrix A("MAT.DAT"); // Formatted  file                 *
//   * Matrix A('*',MAT.BIN"); // Binary file                  *
//   * Matrix A(m,n,B); // submatrix of B                      *
//   * Matrix A(m,n,i,j,B);// submatrix of B from i,j          *
//   * Matrix A = v; // from Vector                            *
//   * Matrix A = R; // from MatrixRight                       *
//   * Matrix A = L; // from MatrixLeft                        *
//   * Matrix A = R; // from MatrixSymm                        *
//   ***********************************************************
//   ***** Access functions:                                   *
//   * i = A.Rows(); // numRows                                *
//   * i = A.Columns(); // numColumns                          *
//   * who = A.WhoAmI();                                       *
//   * xf = A.GetValue(i,j);                                   *
//   * v = GetRow(i);                                          *
//   * v = GetColumn(j);                                       *
//   * xf = A(i,j);                                            *
//   * A(i,j) = xf;                                            *
//   * xf = A[i][j];                                           *
//   * A[i][j] = xf;                                           *
//   * A.SetValue(i,j,xf);                                     *
//   * A.SetRow(i,v);                                          *
//   * A.SetColumn(j,v);                                       *
//   ***********************************************************
//   ***** Assignments:                                        *
//   * A = B;  // Matrix                                       *
//   * A = v;  // Vector                                       *
//   * A = R;  // MatrixRight                                  *
//   * A = L;  // MatrixLeft                                   *
//   * A = S;  // MatrixSymm                                   *
//   ***********************************************************
//   *****  Operators for composing matrices:                  *
//   * A = B&&C; // Adds C beneath B                           *
//   * A = B||C; // Adds C onto the side of B                  *
//   ***********************************************************
//   *****  Operators for tests:                               *
//   * if(A == B)                                              *
//   * if(A != B)                                              *
//   ***********************************************************
//   ***** Implemented operations :                            *
//   * Sum(A,B,&C); // C = A + B;                              *
//   * C = A + B;   // C = A + B;                              *
//   * Sum(A,B,&A); // A = A + B;                              *
//   * Sum(&A,B);   // A = A + B;                              *
//   * A += B;      // A = A + B;                              *
//   * Sum(A,B,&B); // B = A + B;                              *
//   * Sum(A,&B);   // B = A + B;                              *
//   * Sum(A,A,&A); // A = A + A;                              *
//   * Sum(&A);     // A = A + A;                              *
//   * Difference(A,B,&C); // C = A - B;                       *
//   * C = A - B;          // C = A - B;                       *
//   * Difference(A,B,&A); // A = A - B;                       *
//   * Difference(&A,B);   // A = A - B;                       *
//   * A -= B;             // A = A - B;                       *
//   * Difference(A,B,&B); // B = A - B;                       *
//   * Difference(A,&B);   // B = A - B;                       *
//   * Difference(A,A,&A); // A = A - A;                       *
//   * Difference(&A);     // A = A - A;                       *
//   * Minus(A,&B);       // B = -A;                           *
//   * B = -A;            // B = -A;                           *
//   * Minus(&A);         // A = -A;                           *
//   * Product(A,B,&C);   // C = A*B;                          *
//   * C = A*B;                                                *
//   * Product(A,B,&A);   // A = A*B;                          *
//   * Product(&A,B);     // A = A*B;                          *
//   * A *= B;            // A = A*B;                          *
//   * Product(A,B,&B);   // B = A*B;                          *
//   * Product(A,&B);     // B = A*B;                          *
//   * Product(A,A,&A);   // A = A*A;                          *
//   * Product(&A);       // A = A*A;                          *
//   * Product(A,x,&y);   // y = A*x;                          *
//   * y = A*x;                                                *
//   * Product(A,x,&x);   // x = A*x;                          *
//   * Product(A,&x);     // x = A*x;                          *
//   * Product(3.,A,&B);  // B = 3.*A;                         *
//   * B = 3.*A;                                               *
//   * Product(3.,&A);    // A = 3.*A;                         *
//   * A *= 3.;           // A = 3.*A;                         *
//   * TProduct(A,B,&C);  // C = ATB;                          *
//   * C = A%B;           // C = ATB;                          *
//   * TProduct(A,B,&A);  // A = ATB;                          *
//   * TProduct(&A,B);    // A = ATB;                          *
//   * A %= B;            // A = ATB;                          *
//   * TProduct(A,B,&B);  // B = ATB;                          *
//   * TProduct(A,&B);    // B = ATB;                          *
//   * TProduct(A,A,&A);  // A = ATA;                          *
//   * TProduct(&A);      // A = ATA;                          *
//   * A %= A;            // A = ATA;                          *
//   * y = A%x;           // y = ATx;                          *
//   * ProductT(x,y,&A);  // A = xyT;                          *
//   * A = x->*y;         // A = xyT;                          *
//   * ProductT(&A);      // A = AAT                           *
//   * Division(B,3.,&A); // A = B/3.;                         *
//   * A = B/3.;                                               *
//   * Division(&A,3.);   // A = A/3.;                         *
//   * A /= 3.;           // A = A/3.                          *
//   ***********************************************************
//   * Other functions:                                        *
//   * A.Print("Comment");                                     *
//   * A.Save("MAT.DAT");  // formatted                        *
//   * A.Save('*',"MAT.DAT");  // unformatted                  *
//   * xf = A.Max(&imax,&jmax);                                *
//   * xf = A.MaxAbs(&imax,&jmax);                             *
//   * xf = A.Min(&imin,&jmin);                                *
//   * xf = A.MinAbs(&imin,&jmin);                             *
//   * xf = A.NormT();                                         *
//   * xf = A.NormR();                                         *
//   * xf = A.NormC();                                         *
//   * xf = A.NormF();                                         *
//   * xf = A.NormI();                                         *
//   * xf = A.Norm1();                                         *
//   * Delete(&A);                                             *
//   * ChangeDimensions(newr,newc,&A);                         *
//   * Recover(&A,"MAT.DAT");                                  *
//   * Recover(&A,'*',"MAT.DAT");                              *
//   * Transpose(&A);                                          *
//   * SumRankOne(u,vT,&A);      // A = A + uvT;               *
//   * SumRankOne(3.,u,vT,&A);   // A = A + 3.*uvT;            *
//   * SumRankOne(u,vT,3.,&A);   // A = A + uvT/3.;            *
//   * Swap(&A,&B);                                            *
//   ***********************************************************


#ifndef MATRIX_HPP
#define MATRIX_HPP

// preventive statements 
class Vector;
class Matrix;
class MatrixLeft;
class MatrixRight;
class MatrixSymm;

class MatrixSparse;
class Factored;
class FactoredPLR;
class FactoredGauss;
class FactoredQRLQ;
class FactoredQR;
class FactoredLQ;
class FactoredSVD;
class FactoredSymm;

//   ===========================================================
//   ===================   class Matrix   ======================
//   ===========================================================
class Matrix
   {
   friend class Vector;

   friend class MatrixRight;
   friend class MatrixLeft;
   friend class MatrixSymm;

   friend class MatrixSparse;

   friend class Factored;
   friend class FactoredPLR;
   friend class FactoredGauss;
   friend class FactoredQRLQ;
   friend class FactoredQR;
   friend class FactoredLQ;
   friend class FactoredSVD;
   friend class FactoredSymm;

private:
   static const char *const ERROR;
   static int count; // for whoAmI

   float **matrix;
   int numRows,numColumns;
   int size;
   int whoAmI;

   // initialisation constructors
   void Initialize(int mrows,int ncolumns);

   // deinitialisation
   void Deinitialize(void);

   // prepares assignments
   void PrepCopy(int rRows,int rColumns);

   // for assigning and initialising from MatrixRight
   void CopyRight(const MatrixRight & rval);

   // for assigning and initialising from MatrixLeft
   void CopyLeft(const MatrixLeft & rval);

   // for assigning and initialising from MatrixSymm
   void CopySymm(const MatrixSymm & rval);

   // private constructor Matrix A type('*',3,5);
   Matrix(char,int mrows,int columns);

public:
//   ===========================================================
//   *******************   constructors   **********************
//   ===========================================================
   // default // Matrix A;
   Matrix(void);

   // copy-initializer // Matrix A = B;
   Matrix(const Matrix &rval);

   // sizes and initialises at 0 //Matrix A(3,5);
   Matrix(int rows,int columns);

   // sizes and initialises
   // Matrix A(2,3,1.,2.,3.,4.,5.,6.);
   Matrix(int rows,int columns,double a11,...);

   // initialises from array
   // float w[6]={1.,2.,3.,4.,5.6}; Matrix A(3,5,w);
   Matrix(int rows,int columns,float *initvalues);

   // initialises from formatted File 
   // Matrix A("MAT.DAT");
   Matrix(char *filematrix);

   // initialises from binary File 
   // Matrix A('*',"MAT.BIN"); See Save
   Matrix(char,char *filematrix);

   // makes a submatrix of rows,columns
   Matrix(int rows,int columns,const Matrix &rval); 

   // likewise starting from irow,jcol
   Matrix(int rows,int columns,int irow,int jcol,
	  const Matrix &rval);

   // from Vector
   Matrix(const Vector &rval);

   // from MatrixRight
   Matrix(const MatrixRight &rval);

   // from MatrixLeft
   Matrix(const MatrixLeft &rval);

   // from MatrixSymm
   Matrix(const MatrixSymm &rval);

//   ===========================================================
//   ********************   destructor   ***********************
//   ===========================================================
   ~Matrix(void);  

//   ===========================================================
//   ********   Non-modifying access functions   ***************
//   ===========================================================
   // number of rows
   int Rows(void) const
      {return numRows;}

   // number of columns
   int Columns(void) const
      {return numColumns;}

   int WhoAmI(void) const {return whoAmI;}

   // receives the values with control 
   float GetValue(int row,int col) const;

   //takes the row i of a matrix
   Vector GetRow(int i) const;  

   //takes the column j of a matrix
   Vector GetColumn(int j) const; 

//   ===========================================================
//   ************   Modifying access functions   ***************
//   ===========================================================
   // assigns and receives vector values with control
   float &operator () (int row,int col);

   // assigns and receives vector values without control
   float *operator [] (int r)
      {return matrix[r];}

   // assigns values with control 
   void SetValue(int row,int column,float val);

   //substitutes column i with rval
   void SetRow(int i,const Vector &rval);

   // substitutes column j with rval
   void SetColumn(int j,const Vector &rval); 

//   ===========================================================
//   **************  assignment operators **********************
//   ===========================================================
   Matrix &operator =
	(const Matrix &rval);

   Matrix &operator =
	(const Vector &rval);

   Matrix &operator =
	(const MatrixLeft &rval);

   Matrix &operator =
	(const MatrixRight &rval);

   Matrix &operator =
	(const MatrixSymm &rval);

   // transforms a Matrix in Factored
   // Matrix gets destroyed
   friend void CommuteMatrixToFactored
	(Matrix *lval,Factored *rval);

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

      //adds a matrix beneath another
   friend Matrix operator &&
	(const Matrix &lval,const Matrix &rval);

     //add a matrix to the side of another
   friend Matrix operator ||
	(const Matrix &lval,const Matrix &rval);

//   ===========================================================
//   *****************    test operators   *********************
//   ===========================================================
   friend char operator ==
	(const Matrix &lval,const Matrix &rval);

   friend char operator !=
	(const Matrix &lval,const Matrix &rval);

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

//   ===========================================================
//   ***********************   Sum   ***************************
//   ===========================================================
    // Sum(A,B,&C); C = A + B;
   friend void Sum             
      (const Matrix &lval,const Matrix &rval,Matrix *result);

    // C = A + B;
   friend Matrix operator +   
      (const Matrix &lval,const Matrix &rval);

    // Sum(&A,B); A = A + B;
   friend void Sum 
      (Matrix *lvalAndResult,const Matrix &rval);

    // A += B; A = A + B;
   Matrix &operator +=   
      (const Matrix &rval);

 // Sum(B,&A); A = B + A;
   friend void Sum     
      (const Matrix &lval,Matrix *rvalAndResult);

 // Sum(&A); A = A + A;
   friend void Sum   
      (Matrix *lvalRvalAndResult);

//   ===========================================================
//   ********************   Difference   ***********************
//   ===========================================================
   // Difference(A,B,&C); C = A - B;
   friend void Difference
      (const Matrix &lval,const Matrix &rval,Matrix *result);

   // C = A - B;
   friend Matrix operator -
      (const Matrix &lval,const Matrix &rval);

   // Difference(&A,B); A = A - B;
   friend void Difference
      (Matrix *lvalAndResult,const Matrix &rval);

   // A -= B; A = A - B;
   Matrix &operator -=
      (const Matrix &rval);

   // Difference(B,&A); A = B - A;
   friend void Difference
      (const Matrix &lval,Matrix *rvalAndResult);

   // Difference(&A); A = A - A;
   friend void Difference
      (Matrix *lvalRvalAndResult);

//   ===========================================================
//   ***********************   Minus   *************************
//   ===========================================================
   friend void Minus
      (const Matrix &rval,Matrix *result);

   friend Matrix operator -
      (const Matrix &rval);

   friend void Minus
      (Matrix *rvalAndResult);

//   ===========================================================
//   **********************   Product   ************************
//   ===========================================================
   friend void Product // C=A*B;
      (const Matrix &lval,const Matrix &rval,Matrix *result);

   friend Matrix operator *     // C=A*B;
      (const Matrix &lval,const Matrix &rval);

   friend void Product       // A = A * B;
      (Matrix *lvalAndResult,const Matrix &rval);

   Matrix &operator *=        // A = A * B;   
	(const Matrix &rval);

   friend void Product       // B = A * B;
      (const Matrix &lval,Matrix *rvalAndResult);

   friend void Product       // A = A * A;
      (Matrix *lvalRvalAndResult);

   friend void Product     // y =A*x;
      (const Matrix &lval,const Vector &rval,Vector *result);

   friend Vector operator *    // y = A*x;
      (const Matrix &lval,const Vector &rval);

   friend void Product       // A = 3.*B;
      (float lval,const Matrix &rval,Matrix *result);

   friend Matrix operator *      // A = 3.*B;
      (float lval,const Matrix &rval);

   friend void Product       // A = 3.* A;
      (float lval,Matrix *rvalAndResult);

   Matrix &operator *=    // A = 3.*A;
	(float rval);

//   ===========================================================
//   *********************   TProduct   ************************
//   ===========================================================
   friend void TProduct // C=ATB;
      (const Matrix &lval,const Matrix &rval,Matrix *result);

   friend Matrix operator %
	(const Matrix &lval,const Matrix &rval); // ATB

   // TProduct(A,x,&y); y =ATx; y =A%x;
   friend void TProduct(const Matrix &lval,
	 const Vector &rval,Vector *result);

   friend Vector operator %
	(const Matrix &lval,const Vector &rval);

   friend void TProduct       // A = ATB;
      (Matrix *lvalAndResult,const Matrix &rval);

   Matrix &operator %=
	(const Matrix &rval);

   friend void TProduct       // B = ATB;
      (const Matrix &lval,Matrix *rvalAndResult);

   friend void TProduct       // A = A % A;
      (Matrix *lvalRvalAndResult);


   // TProduct(A,&x); x = ATx;
   friend void TProduct
      (const Matrix &lval,Vector *rvalAndresult);


//   ===========================================================
//   *********************   ProductT   ************************
//   ===========================================================
   // ProductT(x,y,&A); A = xyT; A = x->*y;
   friend void ProductT(const Vector &lval,
	const Vector &rval,Matrix *result);

   // A = x->*y; A = xyT;
   friend Matrix operator ->*
	(const Vector &lval,const Vector &rval);  


   friend void ProductT       // A = AAT ;
      (Matrix *lvalRvalAndResult);

//   ===========================================================
//   **********************   Division   ***********************
//   ===========================================================
   friend void Division  // A =B/3.;
      (const Matrix &lval,float rval,Matrix *result);

   friend Matrix operator /      // A= B/3.;
	(const Matrix &lval,float rval);

   friend void Division          // A = A/3.;
      (Matrix *lvalAndResult,float rval);

   Matrix &operator /=            // A /= 3.;
	(float rval);

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

//   ***********************< Print >***************************
   void Print(char *msg="");

//   ************************< Save >***************************
   void Save(char *filematrix); // formatted
   void Save(char,char *filematrix);// binary

//   *********************< Max and Min >***********************
   //to have the position Max(&im,&jc)
   float Max(int *imax=0,int *jmax=0);
   float MaxAbs(int *imax=0,int *jmax=0);

   float Min(int *imin=0,int *jmin=0);
   float MinAbs(int *imin=0,int *jmin=0);

//   ***********************< Norms >***************************
   float NormT(void);
   float NormR(void);
   float NormC(void);
   float NormF(void);
   float NormI(void){return NormR();}
   float Norm1(void){return NormC();}

//   ===========================================================
//   =============   Modifying Functions =======================
//   ===========================================================
   friend void Delete(Matrix *A); // eliminate Matrix
   friend void ChangeDimensions(int rows,int columns,
	       Matrix *result,char zero = 0);

   // recovery of formatted Save 
   friend void Recover(Matrix *A,char *filevector);

   // recovery of binary Save
   friend void Recover(Matrix *A,char,char *filevector);

   friend void Transpose(Matrix *A);
   friend char Inverse(Matrix *A);  // Gauss Jordan

   friend void SumRankOne(const Vector &u,
	  const Vector &vT,Matrix *result);
   friend void SumRankOne(float product,const Vector &u,
	  const Vector &vT,Matrix *result);
   friend void SumRankOne(const Vector &u,
	  const Vector &vT,float divisor,Matrix *result);

   friend void Swap(Matrix *lval,Matrix *rval);
   friend void Solve(const MatrixRight &R,Matrix *BX);
   friend void Solve(const MatrixLeft &R,Matrix *BX);
};

#endif MATRIX_HPP