//   **********************< RIGHT.CPP >************************
//   * square nxn, triangular, Right or Upper matrices         *
//   * 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 "utility.hpp"
#include "vector.hpp"
#include "matrix.hpp"
#include "left.hpp"
#include "right.hpp"
#include "symm.hpp"

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

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

//   ********************< Initialize >*************************
//   * Purpose: Common function for initialising constructors  *
//   ***********************************************************
void MatrixRight::Initialize(int rows,int columns)
   {
   count++;
   whoAmI = count;
   if(rows < 1 || columns < 1)
      {
      numRows = numColumns = size = 0;
      matrix = 0; 
      return;
      }
   if(rows != columns)
      Error("%s%s%sInitialize",
      MatrixRight::ERROR,ERR_CHECK_DIMENSION,ERR_FUNCTION);
   numRows = numColumns = rows;
   size = numRows*(numRows+1)/2+1;
   matrix = new float *[numRows+1]; 
   if(!matrix)
      Error("%s%s%sInitialize",
      MatrixRight::ERROR,ERR_SPACE,ERR_FUNCTION);
   matrix[0] = new float[size]; 
   if(!matrix[0])
      Error("%s%s%sInitialize",
      MatrixRight::ERROR,ERR_SPACE,ERR_FUNCTION);
   matrix[1] = matrix[0];
   for(int i = 2;i <= numRows;i++)
      matrix[i] = &matrix[i-1][numRows-i+1];
   }

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

//   ****************< private constructor >********************
//   * Purpose: Sets size of a MatrixRight without             *
//   *          initialisation. Serving internally             *
//   * Example: MatrixRight R('*',4,5);                        *
//   ***********************************************************
MatrixRight::MatrixRight(char,int rows,int columns)
   {
   Initialize(rows,columns);
   }

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

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

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

//   ****************< Copy-Initializer >***********************
//   * Purpose: Definition and initialisation of MatrixRigh    *
//   *          Used in returns                                *
//   * Examples: MatrixRight A=B;                              *
//   *           MatrixRight A=B+2.*C;                         *
//   *           return A; (used implicitly)                   *
//   ***********************************************************
MatrixRight::MatrixRight(MatrixRight &rval)
   {
   Initialize(rval.numRows,rval.numRows);
   if(numRows != 0)
      memcpy(matrix[0],rval.matrix[0],size*sizeof(float));
   }

//   *************< sizing constructor >************************
//   * Purpose: Setting size and initialising to 0             *
//   *          a MatrixRight                                  *
//   *          Must be able to assign it other values         *
//   * Description: Uses the function memset                   *
//   * Example: MatrixRight A(5,5);                            *
//   *          A(1,1)=3.; A(1,2)=7.;....                      *
//   ***********************************************************
MatrixRight::MatrixRight(int rows,int columns)
   {
   Initialize(rows,columns);
   if(numRows != 0)
      memset(matrix[0],0,size*sizeof(float));
   }

//   *************< sizing and initialising  >******************
//   * Purpose: Initialising a MatrixRight                     *
//   * Example: MatrixRight A(2,2,1.,2.,3.);                   *
//   ***********************************************************
MatrixRight::MatrixRight(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 MatrixRight from an array       *
//   * Example: float w[]={1.,2.,3.,                           *
//   *                        4.,5.,                           *
//   *                           6.};                          *
//   *          MatrixRight A(3,3,w);                          *
//   ***********************************************************
MatrixRight::MatrixRight(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 MatrixRight from FILE           *
//   * Example: MatrixRight A("MAT.DAT");                      *
//   ***********************************************************
MatrixRight::MatrixRight(char *filematrix)
   {
   if((bzzFileIn=fopen(filematrix,"r"))==NULL)
       Error("%s%s%s%sConstructor",
       MatrixRight::ERROR,ERR_OPEN_FILE,filematrix,ERR_FUNCTION);
   int rows,columns;
   fscanf(bzzFileIn,"%d %d",&rows,&columns);
   if(rows < 1 || columns < 1)
       Error("%s%s%sFILE",
       MatrixRight::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);
   }
//   *******************< binary FILE >*************************
//   * Purpose: Initialising a MatrixRight from                *
//   *          unformatted FILE                               *
//   * Example: MatrixRight A('*',"MAT.DAT");                  *
//   ***********************************************************
MatrixRight::MatrixRight(char,char *filematrix)
   {
   if((bzzFileIn=fopen(filematrix,"rb"))==NULL)
       Error("%s%s%s%sFILE",
       MatrixRight::ERROR,ERR_OPEN_FILE,filematrix,ERR_CONSTRUCTOR);
   int rows,columns;
   if(fread(&rows,sizeof(int),1,bzzFileIn) != 1)
      Error("%sFILE",MatrixRight::ERROR);
   if(fread(&columns,sizeof(int),1,bzzFileIn) != 1)
      Error("%sFILE",MatrixRight::ERROR);
   if(rows < 1 || columns < 1)
       Error("%s%s%sFILE",
       MatrixRight::ERROR,ERR_CHECK_DIMENSION,
       ERR_CONSTRUCTOR);
   Initialize(rows,columns);
   if(fread(matrix[1],sizeof(float)*size,1,bzzFileIn) != 1)
      Error("%sFILE",MatrixRight::ERROR);
   fclose(bzzFileIn);
   }

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

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

//   *********************< GetValue >**************************
//   * Purpose: Receiving values with control                  *
//   * Example: x = R.GetValue(i,j);                           *
//   ***********************************************************
float MatrixRight::GetValue(int row,int column) const
   {
   if(row < 1 || row > numRows)
      Error("%s%s%sGetValue",
      MatrixRight::ERROR,ERR_RANGE,ERR_FUNCTION);
   if(column < 1 || column > numRows)
      Error("%s%s%sGetValue",
      MatrixRight::ERROR,ERR_RANGE,ERR_FUNCTION);
   if(column < row)return 0.;
   return matrix[row][column];
   }

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

//   *********************< SetValue >**************************
//   * Purpose: Assigning values with control                  *
//   ***********************************************************
void MatrixRight::SetValue(int row,int column,float val)
   {
   if(row < 1 || row > numRows)
      Error("%s%s%sSetValue",
      MatrixRight::ERROR,ERR_RANGE,ERR_FUNCTION);
   if(column < 1 || column > numRows)
      Error("%s%s%sSetValue",
      MatrixRight::ERROR,ERR_RANGE,ERR_FUNCTION);
   if(column < row)
      Error("%s\ncolumn<row%sSetValue",
      MatrixRight::ERROR,ERR_FUNCTION);
   matrix[row][column]=val;
   }

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

//   ===========================================================
//   ***************   assignment operator   *******************
//   ===========================================================

//   *************************< = >*****************************
//   * Purpose: Assignment                                     *
//   * Example: A=B;   // simple                               *
//   *          A=C+D;    // with operations                   *
//   *          A=B=C+D;  // multiple                          *
//   ***********************************************************
MatrixRight &MatrixRight::operator =(const MatrixRight &rval)
   {
   int who = whoAmI;
   if(numRows != rval.numRows)
      {
      Deinitialize();
      Initialize(rval.numRows,rval.numRows);
      count--;
      }
   whoAmI = who;
   if(numRows != 0)
      memcpy(matrix[0],rval.matrix[0],size*sizeof(float));
   return *this;
}

//   ===========================================================
//   ***************   operators for tests   *******************
//   ===========================================================

//   ************************< == >*****************************
//   * Purpose: Checking if two MatrixRights are identical     *
//   * Description: Returns 1 if they are identical            *
//   * Example: if(A==B)...                                    *
//   ***********************************************************
char operator == (const MatrixRight &lval,
                const MatrixRight &rval)
   {
   if(lval.whoAmI == rval.whoAmI)return 1;
   char ch = 1;
   if(lval.numRows != rval.numRows)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: Checking if two MatrixRights are different     *
//   * Description: Returns 1 if they are different            *
//   * Example: if(A!=B)...                                    *
//   ***********************************************************
char operator != (const MatrixRight &lval,
                const MatrixRight &rval)
   {
   if(lval.whoAmI == rval.whoAmI) return 0;
   char ch = 0;
   if(lval.numRows != rval.numRows) 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 MatrixRight &lval,const MatrixRight &rval,
        MatrixRight *result)
   {
   if(lval.numRows != rval.numRows)
      Error("%s%s%s+",
      MatrixRight::ERROR,ERR_CHECK_DIMENSION,ERR_OPERATOR);
   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,lval.matrix[1],rval.matrix[1],
            result->matrix[1]);
      }
   }

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

//   ************************< Sum >****************************
//   * Purpose: Addition between Left and Right matrices       *
//   * Example: Sum(L,R,&C); C = L + R;                        *
//   ***********************************************************
void Sum(MatrixLeft &lval,MatrixRight &rval,Matrix *result)
   {
   if(lval.numRows != rval.numRows)
      Error("%s%s%s+",
      MatrixRight::ERROR,ERR_CHECK_DIMENSION,ERR_OPERATOR);
   ChangeDimensions(lval.numRows,lval.numColumns,
                  result,1);
   for(int i = 1;i <= lval.numRows;i++)
      {
      for(int j = 1;j <= lval.numColumns;j++)
        {
        if(i < j)
           (*result)[i][j] = lval[i][j];
        else if(i > j)
           (*result)[i][j] = rval[i][j];
        else
           (*result)[i][j] = Control(lval[i][j] 
                                  + rval[i][j]);
        }
      } 
   }

//   ************************< Sum >****************************
//   * Purpose: Addition between Left and Right matrices       *
//   * Example: Sum(R,L,&C); C = R + L;                        *
//   ***********************************************************
void Sum(MatrixRight &lval,MatrixLeft &rval,Matrix *result)
   {
   if(lval.numRows != rval.numRows)
      Error("%s%s%s+",
      MatrixRight::ERROR,ERR_CHECK_DIMENSION,ERR_OPERATOR);
   ChangeDimensions(lval.numRows,lval.numColumns,
                  result,1);
   for(int i = 1;i <= lval.numRows;i++)
      {
      for(int j = 1;j <= lval.numColumns;j++)
        {
        if(i < j)
           (*result)[i][j] = rval[i][j];
        else if(i > j)
           (*result)[i][j] = lval[i][j];
        else
           (*result)[i][j] = Control(lval[i][j] 
                                  + rval[i][j]);
        }
      } 
   }

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


//   ************************< += >*****************************
//   * Purpose: Addition of two Right matrices                 *
//   * Example: A += B;                                        *
//   ***********************************************************
MatrixRight &MatrixRight::operator +=
            (const MatrixRight &rval)
   {
  if(numRows != rval.numRows)
      Error("%s%s%s+",
      MatrixRight::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 MatrixRight &lval,MatrixRight *rvalAndResult)
   {
   (*rvalAndResult) += lval;
   }

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

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

//   *************************< - >*****************************
//   * Purpose: Difference between Right matrices              *
//   ***********************************************************
MatrixRight operator -
       (const MatrixRight &lval,const MatrixRight &rval)
   {
   if(lval.numRows != rval.numRows)
      Error("%s%s%s-",
      MatrixRight::ERROR,ERR_CHECK_DIMENSION,ERR_OPERATOR);
   MatrixRight result('*',lval.numRows,lval.numRows);
   Difference(rval.size,lval.matrix[1],rval.matrix[1],
            result.matrix[1]);
   return result;
   }

//   ************************< -= >*****************************
//   * Purpose: Difference between two Right matrices          *
//   * Example: A-=B;                                          *
//   ***********************************************************
MatrixRight &MatrixRight::operator -=
            (const MatrixRight &rval)
   {
   if(numRows != rval.numRows)
      Error("%s%s%s-",
      MatrixRight::ERROR,ERR_CHECK_DIMENSION,ERR_OPERATOR);
   Difference(rval.size,matrix[1],rval.matrix[1]);
   return *this;
   }

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

//   *************************< - >*****************************
//   * Purpose: Unary minus                                    *
//   * Example: A=-B;                                          *
//   ***********************************************************
MatrixRight operator -(const MatrixRight &rval) //unary minus
   {
   MatrixRight result('*',rval.numRows,rval.numRows);
   float *w = rval.matrix[1];
   float *r = result.matrix[1];
   for(int i = 1;i < result.size;i++)
   r[i] = -w[i];
   return result;
   }

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

//   *************************< * >*****************************
//   * Purpose: Product of two Right matrices                  *
//   * Description: It produces a matrixRight                  *
//   * Example: C = A*B;                                       *
//   ***********************************************************
MatrixRight operator *
       (const MatrixRight &lval,const MatrixRight &rval)
   {
   if(lval.numRows != rval.numRows)
      Error("%s%s%s*",
      MatrixRight::ERROR,ERR_CHECK_DIMENSION,ERR_OPERATOR);
   MatrixRight result('*',lval.numRows,lval.numRows);
   float **b=rval.matrix;
   float *v=result.matrix[1];
   for(int row=1;row<=lval.numRows;row++)
       {
       for(int column=row;column<=rval.numRows;column++)
         {
         float *w = lval.matrix[row] + row;
         double sum = 0.;
         for(int i = row,k = 1;k <= column-row+1;i++,k++)
            {
            sum += *w++ * b[i][column];
            }
         *++v = Control(sum); 
         }
      }
   return result;
   }

//   *************************< * >*****************************
//   * Purpose: Product of a Right matrix with a vector        *
//   * Example: y = R*x;                                       *
//   ***********************************************************
Vector operator *
       (const MatrixRight &lval,const Vector &rval)
   {
   if(lval.numColumns != rval.dimensions)
       Error("%s%s%s*",
       MatrixRight::ERROR,ERR_CHECK_DIMENSION,ERR_OPERATOR);
    Vector result('*',lval.numRows);
    float *r = result.vector;
    for(int i = 1;i <= lval.numRows;i++)
       r[i] = Dot(rval.dimensions-i+1,
             lval.matrix[i] + i,rval.vector + i);
    return result;
   }

//   *************************< * >*****************************
//   * Purpose: Product of a float with a Right matrix         *
//   * Example: C = 3.*B;                                      *
//   ***********************************************************
MatrixRight operator * (float lval,const MatrixRight &rval)
   {
   MatrixRight result('*',rval.numRows,rval.numColumns);
   Product(rval.size - 1,lval,rval.matrix[1]+1,
          result.matrix[1]+1);
   return result;
   }

//   ************************< *= >*****************************
//   * Purpose: Product of a float with a Right matrix         *
//   * Example: R *= 3.;                                       *
//   ***********************************************************
MatrixRight &MatrixRight::operator *= (float rval)
   {
   Product(size - 1,rval,matrix[1] + 1);//in UTILITY.CPP
   return *this;
   }

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

//   **********************< TProduct >*************************
//   * Purpose: Product of a transposed right matrix           *
//   *          multiplied by itself                           *
//   * Description: The result is MatrixSymm                   *
//   * Example: TProduct(R,&S);                                *
//   ***********************************************************
void TProduct(const MatrixRight &lvalAndRval,
                  MatrixSymm *result)
   {
   float **r = lvalAndRval.matrix;
   int n = lvalAndRval.numRows;
   ChangeDimensions(n,n,result,1);
   float **s = result->matrix;
   for(int i=1;i <= n;i++)
      {
      for(int j=1;j <= i;j++)
        {
        double sum = 0.;
        for(int k=1;k <= i && k <=j;k++)
           sum += r[k][i]*r[k][j];
        s[i][j] = Control(sum);
        }
      }
   }
//   ===========================================================
//   *********************   ProductT   ************************
//   ===========================================================

//   **********************< ProductT >*************************
//   * Purpose: Product of a Right matrix                      *
//   *          with its transpose                             *
//   * Description: The result is MatrixSymm                   *
//   * Example: ProductT(R,&S);                                *
//   ***********************************************************
void ProductT(const MatrixRight &lvalAndRval,
                  MatrixSymm *result)
   {
   float **r = lvalAndRval.matrix;
   int n = lvalAndRval.numRows;
   ChangeDimensions(n,n,result,1);
   float **s = result->matrix;
   for(int i=1;i <= n;i++)
      {
      for(int j=1;j <= i;j++)
        {
        double sum = 0.;
        for(int k=i;k <= n && k>=j;k++)
           sum += r[i][k]*r[j][k];
        s[i][j] = Control(sum);
        }
      }
   }

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

//   *************************< / >*****************************
//   * Purpose: Dividing by a float                            *
//   * Example: C = B/3.;                                      *
//   ***********************************************************
MatrixRight operator / (const MatrixRight &lval,float rval)
   {
   MatrixRight result('*',lval.numRows,lval.numColumns);
   Division(lval.size - 1,lval.matrix[1]+1,rval,
           result.matrix[1]+1);
   return result;
   }

//   ************************< /= >*****************************
//   * Purpose: Dividing by a float                            *
//   * Example: R /= 3.;                                       *
//   ***********************************************************
MatrixRight &MatrixRight::operator /= (float rval)
   {
   Division(size - 1,matrix[1] + 1,rval);//in UTILITY.CPP
   return *this;
   }

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

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

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

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

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

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

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

//   ***********************< NormR >***************************
//   * Purpose: Max of the sums of the abs values              *
//   *          of coeff. rows                                 *
//   ***********************************************************
float MatrixRight::NormR(void)
   {
   float norm=0.,aus;
   float *w=matrix[0];
   for(int i=1;i<=numRows;i++)
      {
      aus=0.;
      for(int j=i;j<=numRows;j++)aus+=Abs(*++w);
      if(aus>norm)norm=aus;
      }
   return norm;
   }

//   ***********************< NormC >***************************
//   * Purpose: Norm C                                         *
//   ***********************************************************
float MatrixRight::NormC(void)
   {
   float norm=0.,aus;
   norm=0.;
   for(int j=1;j<=numRows;j++)
      {
      aus=0.;
      for(int i=1;i<=j;i++)aus+=Abs(matrix[i][j]);
      if(aus>norm)norm=aus;
      }
   return norm;
   }

//   ***********************< NormF >***************************
//   * Purpose: Calculating Frobenius' norm                    *
//   * Description: Uses the functions SqrtSumSqr              *
//   * Example: float normF = R.NormF();                        *
//   ***********************************************************
float MatrixRight::NormF(void)
   {
   return SqrtSumSqr(size-1,matrix[0]+1);
   }

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

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

//   ***************< ChangeDimensions >************************
//   * Purpose: Changing the dimensions of a MatrixRight       *
//   * Example: ChangeDimensions(row,columns,&A);              *
//   ***********************************************************
void ChangeDimensions(int rows,int columns,
                    MatrixRight *result,char zero)
   {
   if(rows != columns)
       Error("%s%s%sChangeDimensions",
       MatrixRight::ERROR,ERR_CHECK_DIMENSION,ERR_FUNCTION);
   int who = result->whoAmI;
   if(rows != result->numRows ||
      columns != result->numColumns)
        {
        result->Deinitialize();
        result->Initialize(rows,columns);
        MatrixRight::count--;
        }
   result->whoAmI = who;
   if(zero == 0 && result->size != 0)
      memset(result->matrix[0],0,
            result->size*sizeof(float));
   }

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

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

//   **********************< Inverse >**************************
//   * Purpose: Executing the inverse of aMatrixRight          *
//   *          The inverse will also be a MatrixRight         *
//   * Example: Inverse(&R);                                   *
//   ***********************************************************
char Inverse(MatrixRight *R)
   {
   int row,column,k,ok=1;
   float **r = R->matrix;
   float *v = R->matrix[1];
   float *w;
   int n = R->numRows;
   double sum;
   for(row = 1;row <= n;row++)
      {
      if(r[row][row] == 0.)
         {
         Message
         ("\nSingular matrix in MatrixRight::Inverse\n");
         ok = 0;
         r[row][row] = 1.;
         for(column = row+1;column <= n;column++)
            r[row][column] = 0.;
         }
      else r[row][row] = 1./r[row][row];
      }
   for(row = 1;row <= n-1;row++)
      {
      v++;
      for(column = row+1;column <= n;column++)
        {
        w = R->matrix[row] + row;
        sum = 0.;
        for(k = row;k <= column-1;k++)
           sum += *w++ * r[k][column];
        *++v = Control(-sum * r[column][column]);
        }
      }
   return ok;
   }

//   ************************< Swap >***************************
//   * Purpose: Swapping the elements of any two Matrices      *
//   * Description: Provides an efficient method of copying    *
//   * Example: Swap(&A,&B);                                   *
//   ***********************************************************
void Swap(MatrixRight *lval,MatrixRight *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);
   }

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

//   ***********************< Solve >***************************
//   * Purpose: Solving a Right system R*bx = bx               *
//   * Description: Given the Vector bx it returns             *
//   *              the solution in bx                         *
//   * Example: Solve(R,&bx);                                  *
//   ***********************************************************
void Solve(const MatrixRight &R,Vector *bx)
   {
   if(R.numRows != bx->Size())
      Error("%s%s%sSolve",
      MatrixRight::ERROR,ERR_CHECK_DIMENSION,ERR_FUNCTION);
   double sum;
   float *r,tmp,*x = &(*bx)[0];
   float **mat = R.matrix;
   for(int i = R.numRows;i >= 1;i--)
     {
     r = mat[i] + i;
     tmp = *r;
     if(tmp == 0.)
       {
       Message("%sSingular Matrix",MatrixRight::ERROR);
       x[i]=0.;
       }
     else
       {
       sum = x[i];
       for(int j = i+1;j <= R.numColumns;j++)
          sum -= *++r * x[j];
       x[i]=Control(sum/tmp);
       }
     }
  }

//   ***********************< Solve >***************************
//   * Purpose: Solving a Right system R*b = x                 *
//   * Description: Given the Vectors b and x  it returns      *
//   *              the solution in x                          *
//   * Example:  Solve(R,b,&x);                                *
//   ***********************************************************
void Solve(const MatrixRight &R,const Vector &b,Vector *x)
   {
   if(b.WhoAmI() == x->WhoAmI())
      Solve(R,x);
   else
      {
      *x = b;
      Solve(R,x);
      }
   }

//   ***********************< Solve >***************************
//   * Purpose: Solving a Right system for a Matrix BX         *
//   * Description: Substitutes the solution in BX             *
//   * Example: R.Solve(&BX);                                  *
//   ***********************************************************
void Solve(const MatrixRight &R,Matrix *BX)
   {
   if(R.numRows != BX->Rows())
      Error("%s%s%sSolve",
      MatrixRight::ERROR,ERR_CHECK_DIMENSION,ERR_FUNCTION);
   double sum;
   float tmp,*r,**bx = BX->matrix;
   float **mat = R.matrix;
   for(int k = 1;k <= BX->Columns();k++)
      {
      for(int i = R.numRows;i >= 1;i--)
        {
        r = mat[i] + i;
        tmp = *r;
        if(tmp == 0.)
           {
           Message("%sSingular Matrix",MatrixRight::ERROR);
           bx[i][k] = 0.;
           }
        else
           {
           sum = bx[i][k];
           for(int j = i+1;j <= R.numColumns;j++)
           sum -= *++r * bx[j][k];
           bx[i][k] = Control(sum/tmp);
           }
        }
      }
   }

//   ***********************< Solve >***************************
//   * Purpose: Solving a Right system for a Matrix B          *
//   * Description: Puts the solution in X                     *
//   * Example: Solve(R,B,&X);                                 *
//   ***********************************************************
void Solve(const MatrixRight &R,const Matrix &B,Matrix *X)
   {
   if(B.WhoAmI() == X->WhoAmI())
      Solve(R,X);
   else
      {
      *X = B;
      Solve(R,X);
      }
   }

//   ******************< TransposeSolve >***********************
//   * Purpose: Solving a Right system RT*bx = bx              *
//   * Description: Given Vector bx it returns                 *
//   *              the solution in bx                         *
//   * Example: TransposeSolve(R,&bx);                         *
//   ***********************************************************
void TransposeSolve(const MatrixRight &R,Vector *bx)
   {
   if(R.numRows != bx->Size())
      Error("%s%s%sTransposeSolve",
      MatrixRight::ERROR,ERR_CHECK_DIMENSION,ERR_FUNCTION);
   double sum;
   float tmp,*x = &(*bx)[0];
   float **mat = R.matrix;
   int n = R.numRows;
   for(int i = 1;i <= n;i++)
     {
     tmp = mat[i][i];
     if(tmp == 0.)
       {
       Message("%sSingular Matrix",MatrixRight::ERROR);
       x[i]=0.;
       }
     else
       {
       sum = x[i];
       for(int j = 1;j < i;j++)
          sum -= mat[j][i] * x[j];
       x[i]=Control(sum/tmp);
       }
     }
  }

//   ******************< TransposeSolve >***********************
//   * Purpose: Solving a Right system RT*b = x                *
//   * Description: Given the Vectors b and x it returns       *
//   *              the solution in x                          *
//   * Example: TransposeSolve(R,b,&x);                        *
//   ***********************************************************
void TransposeSolve(const MatrixRight &R,
                  const Vector &b,Vector *x)
   {
   if(b.WhoAmI() == x->WhoAmI())
      TransposeSolve(R,x);
   else
      {
      *x = b;
      TransposeSolve(R,x);
      }
   }


//   ********************< Determinant >************************
//   * Purpose: Calculating the determinant                    *
//   * Example: R.Determinant();                               *
//   ***********************************************************
double MatrixRight::Determinant(void)
   {
   double determinant = 1.;
   for(int i = 1;i <= numRows;i++)
      determinant *= matrix[i][i];
   return determinant;
   } 

//   ******************< ConditionNumber >**********************
//   * Purpose: Calculating the condition number               *
//   * Example: R.ConditionNumber();                           *
//   ***********************************************************
float MatrixRight::ConditionNumber(void)
   {
   double max = TINY_FLOAT;
   double min = BIG_FLOAT;
   for(int i = 1;i <= numRows;i++)
      {
      double aux = Abs(matrix[i][i]);
      if(max < aux) max = aux;
      if(min > aux) min = aux;
      }
   if(min == 0.)return BIG_FLOAT;
   else return Control(max/min);
   }