//   ======================< SYMM.CPP >=========================
//   * MatrixSymm class for symmetrical 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 "symm.hpp"

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

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

//   ********************< Initialize >*************************
//   * Purpose: Common function for initialising constructors  *
//   ***********************************************************
void MatrixSymm::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",
      MatrixSymm::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",
      MatrixSymm::ERROR,ERR_SPACE,ERR_FUNCTION);
   matrix[0] = new float[size]; 
   if(!matrix[0])
      Error("%s%s%sInitialize",
      MatrixSymm::ERROR,ERR_SPACE,ERR_FUNCTION);
   matrix[1] = matrix[0];
   for(int i = 2;i <= numRows;i++)
      matrix[i] = &matrix[i-1][i-1];
   }

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

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

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

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

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

//   ****************< Copy-Initializer >***********************
//   * Purpose: Defining and initialising MatrixSymm           *
//   *          Used in returns                                *
//   * Examples: MatrixSymm A=B;                               *
//   *           MatrixSymm A=B+2.*C;                          *
//   *           return A; (used implicitly)                   *
//   ***********************************************************
MatrixSymm::MatrixSymm(MatrixSymm &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 MatrixSymm*
//   *          Must then be assigned other values             *
//   * Example: MatrixSymm A(5,5);                             *
//   *          A(1,1)=3.; A(2,1)=7.;....                      *
//   ***********************************************************
MatrixSymm::MatrixSymm(int rows,int columns)
   {
   Initialize(rows,columns);
   if(numRows != 0)
      memset(matrix[0],0,size*sizeof(float));
   }

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

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

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

//   *********************< GetValue >**************************
//   * Purpose: Receiving values with control                  *
//   * Example: S.GetValue(1,2);                               *
//   ***********************************************************
float MatrixSymm::GetValue(int row,int col) const
   {
   if(row < 1 || row > numRows)
      Error("%s%s%sGetValue",
      MatrixSymm::ERROR,ERR_RANGE,ERR_FUNCTION);
   if(col < 1 || col > numRows)
      Error("%s%s%sGetValue",
      MatrixSymm::ERROR,ERR_RANGE,ERR_FUNCTION);
   if(row > col) return matrix[row][col];
   else return matrix[col][row];
   }

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

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

//   ************************< () >*****************************
//   * Purpose: Get and set with control                       *
//   * Example: A(1,2)=3.;//Set                                *
//   *          a=A(4,5); // Get                               *
//   ***********************************************************
float &MatrixSymm::operator ()(int row,int col)
   {
   if( (row < 1 || row > numRows)
	      ||
     (col < 1 || col > numColumns) )
       Error("%s%s%s()",
       MatrixSymm::ERROR,ERR_RANGE,ERR_OPERATOR);
   if(row > col)return matrix[row][col];
   else return matrix[col][row];
}

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

//   *************************< = >*****************************
//   * Purpose: Assignment                                     *
//   * Example: A=B;   // simple                               *
//   *          A=C+D;    // with operations                   *
//   *          A=B=C+D;  // multiple                          *
//   ***********************************************************
MatrixSymm &MatrixSymm::operator =(const MatrixSymm &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 MatrixSymms are identical      *
//   * Description: Returns 1 if they are identical            *
//   * Example: if(A==B)...                                    *
//   ***********************************************************
char operator == (const MatrixSymm &lval,
                const MatrixSymm &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: Check if two MatrixSymms are different         *
//   * Description: Returns 1 if they are different            *
//   * Example: if(A!=B)...                                    *
//   ***********************************************************
char operator != (const MatrixSymm &lval,
                const MatrixSymm &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   ***************************
//   ===========================================================

//   *************************< + >*****************************
//   * Purpose: Addition between matrices                      *
//   * Description: Uses Sum from utility.cpp                  *
//   * Example: C = A + B;                                     *
//   ***********************************************************
MatrixSymm operator + (const MatrixSymm &lval,
                     const MatrixSymm &rval)
   {
   if(lval.numRows != rval.numRows)
      Error("%s%s%s+",
      MatrixSymm::ERROR,ERR_CHECK_DIMENSION,ERR_OPERATOR);
   MatrixSymm result('*',lval.numRows,lval.numRows);
   Sum(rval.size,lval.matrix[1],rval.matrix[1],
            result.matrix[1]);
   return result;
   }

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

//   *************************< - >*****************************
//   * Purpose: Difference between matrices                    *
//   * Description: Uses Difference from utility.cpp           *
//   * Example:  A=B-C;                                        *
//   ***********************************************************
MatrixSymm operator - (const MatrixSymm &lval,
                     const MatrixSymm &rval)
   {
   if(lval.numRows != rval.numRows)
      Error("%s%s%s-",
      MatrixSymm::ERROR,ERR_CHECK_DIMENSION,ERR_OPERATOR);
   MatrixSymm result('*',lval.numRows,lval.numRows);
   Difference(rval.size,lval.matrix[1],rval.matrix[1],
            result.matrix[1]);
   return result;
   }

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

//   ***********************< Print >***************************
//   * Purpose: Printing                                       *
//   ***********************************************************
void MatrixSymm::Print(char *msg)
   {
   fprintf(bzzFileOut,"\nMatrixSymm No.%d",whoAmI);
   if(*msg)fprintf(bzzFileOut,"  %s",msg);
   fprintf(bzzFileOut,"\nrows %d cols %d\n",
          numRows,numRows);
   for(int row=1;row<=numRows;row++)
      {
      for(int col=1;col<=numRows;col++)
        {
        if(row>col)
           fprintf(bzzFileOut," %12.4e",matrix[row][col]);
        else
           fprintf(bzzFileOut," %12.4e",matrix[col][row]);
        }
      fprintf(bzzFileOut,"\n");
      }
   }
//   ************************< Save >***************************
//   * Purpose: Saving a MatrixSymm on formatted file          *
//   * Description: Can be read                                *
//   *              See TempFile in utility                    *
//   * Example: A.Save("Mat.Dat");                             *
//   ***********************************************************
void MatrixSymm::Save(char *filematrix) // formatted
   {
   FILE *bzz;
   if((bzz = fopen(filematrix,"w")) == NULL)
       Error("%s%s%s%sSave",
       MatrixSymm::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 <= row;column++)
        fprintf(bzz," %15.7e",matrix[row][column]);
      fprintf(bzz,"\n");
      }
   fclose(bzz);
   }

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


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

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

//   ***************< ChangeDimensions >************************
//   * Purpose: Change the dimensions of a MatrixSymm          *
//   * Example: ChangeDimensions(row,columns,&A);              *
//   ***********************************************************
void ChangeDimensions(int rows,int columns,
                    MatrixSymm *result,char zero)
   {
   if(rows != columns)
       Error("%s%s%sChangeDimensions",
       MatrixSymm::ERROR,ERR_CHECK_DIMENSION,ERR_FUNCTION);
   int who = result->whoAmI;
   if(rows != result->numRows ||
      columns != result->numColumns)
        {
        result->Deinitialize();
        result->Initialize(rows,columns);
        MatrixSymm::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(MatrixSymm *L,char *filematrix)
   {
   if((bzzFileIn = fopen(filematrix,"r")) == NULL)
       Error("%s%s%s%sRecover",
       MatrixSymm::ERROR,ERR_OPEN_FILE,filematrix,ERR_FUNCTION);
   int rows,columns;
   fscanf(bzzFileIn,"%d %d",&rows,&columns);
   ChangeDimensions(rows,columns,L,1);
   float *w = L->matrix[0];
   for(int i = 1;i < L->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(MatrixSymm *L,char,char *filematrix)
   {
   if((bzzFileIn=fopen(filematrix,"rb"))==NULL)
       Error("%s%s%s%sFILE",
       MatrixSymm::ERROR,ERR_OPEN_FILE,filematrix,ERR_CONSTRUCTOR);
   int rows,columns;
   if(fread(&rows,sizeof(int),1,bzzFileIn) != 1)
      Error("%sFILE",MatrixSymm::ERROR);
   if(fread(&columns,sizeof(int),1,bzzFileIn) != 1)
      Error("%sFILE",MatrixSymm::ERROR);
   ChangeDimensions(rows,columns,L,1);
   if(rows < 1 || columns < 1)
      {fclose(bzzFileIn);return;}
   if(fread(L->matrix[1],sizeof(float)*(L->size),
              1,bzzFileIn) != 1)
        Error("%s%s%sRecover",
        MatrixSymm::ERROR,ERR_READING_FILE,ERR_FUNCTION);
   fclose(bzzFileIn);
   }

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