//   =====================< SPARSE.CPP >========================
//   * Class MatrixSparse                                      *
//   * Description: Chapter 13                                 *
//   ===========================================================

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

#include <alloc.h>

#include "utility.hpp"
#include "vector.hpp"
#include "sparse.hpp"

const char *const MatrixSparse::ERROR =
          "\n========>>> MatrixSparse error!!!!!!\n";

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

//   *************************< Initialize >********************
//   * Purpose: Common function for initialisation             *
//   * Description: Chapter 13                                 *
//   * Example: Initialize(rows,columns)                       *
//   ***********************************************************
void MatrixSparse::Initialize(int rows,int columns)
   {
   if(rows < 1 || columns < 1)
      {
      numRows = numColumns = 0;
      elementRow = 0;
      return;
      }
   numRows = rows;
   numColumns = columns;
   elementRow = new Element *[numRows + 1];
   if(!elementRow)
      Error("%s%s%sInitialize",
           MatrixSparse::ERROR,ERR_SPACE,ERR_FUNCTION);

   for(int row = 0;row <= numRows;row++)elementRow[row] = 0;
   return;
   }

//   *************************< Copy >**************************
//   * Purpose: Copying a sparse matrix                        *
//   * Description: Used in copy-initializer                   *
//   *              and in assignments                         *
//   ***********************************************************
void MatrixSparse::Copy(const MatrixSparse &rval)
   {
   for(int row = 1;row <= numRows;row++)
      {
      Element *elem,*elemRval;
      int first = 1;
      for(elemRval = rval.elementRow[row];
                     elemRval != 0;elemRval = elemRval->next)
        {
        Element *newElement = new Element;
        if(!newElement)
        Error("%s%s%sCopy",
           MatrixSparse::ERROR,ERR_SPACE,ERR_FUNCTION);
        newElement->column = elemRval->column;
        newElement->value = elemRval->value;
        newElement->next = 0;
        if(first == 1)
           {
           elem = elementRow[row] = newElement;
           first = 0;
           }
        else
           {
           elem->next = newElement;
           elem = elem->next; 
           }
        }
      }
   }

//   *******************< InsertElement >***********************
//   * Purpose: Inserting an element                           *
//   * Description: Used internally                            *
//   ***********************************************************
float &MatrixSparse::InsertElement
            (Element *elem,int row,int column,int first)
   {
   Element *newElement = new Element;
   if(!newElement)
      Error("%s%s%sInsertElement",
           MatrixSparse::ERROR,ERR_SPACE,ERR_FUNCTION);
   newElement->column = column;
   newElement->value = 0.;
   if(first == 1)
      {
      newElement->next = elem;
      elementRow[row] = newElement;
      }
   else
      {
      newElement->next = elem->next;  // inserted next
      elem->next = newElement;
      }
   return newElement->value;
   }

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

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

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

//   ******************< Copy-Initialize >**********************
//   * Purpose: For definition and initialisation              *
//   *          of MatrixSparse                                *
//   * Examples: MatrixSparse A=B;                             *
//   ***********************************************************
MatrixSparse::MatrixSparse(const MatrixSparse &rval)
   {
   Initialize(rval.numRows,rval.numColumns);
   Copy(rval);
   }

//   ********************< dimensions >*************************
//   * Purpose: Providing the dimensions of a matrix           *
//   * Examples: MatrixSparse A(120,350);                      *
//   ***********************************************************
MatrixSparse::MatrixSparse(int rows,int columns)
   {
   Initialize(rows,columns);
   }

//   ===========================================================
//   ********************   destructor   ***********************
//   ===========================================================
MatrixSparse::~MatrixSparse(void)
   {
   DeleteMatrix();
   }

//   ===========================================================
//   *****************    Access functions *********************
//   ===========================================================

//   ********************< operator () >************************
//   * Purpose: Receiving and assigning values with control    *
//   * Example: x = A(1,5); A(3,7) = 5.;                       *
//   ***********************************************************
float &MatrixSparse::operator () (int row,int column)
   {
   if(row < 1 || row > numRows || 
      column < 1 || column > numColumns)
        Error("%s%s%s()",
        MatrixSparse::ERROR,
        ERR_CHECK_DIMENSION,ERR_OPERATOR);
   Element *elem = elementRow[row];
   if(elem == 0 || column < elem->column)
      return InsertElement(elem,row,column,1);
   if(column == elem->column)
      return elem->value;
   for(;elem->next != 0;elem = elem->next)
      {
      if(column == elem->next->column)
        return elem->next->value;
      else if(column < elem->next->column)break;
      }
   return InsertElement(elem,row,column,0);
   }

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

//   *************************< = >*****************************
//   * Purpose: Assignment of a MatrixSparse                   *
//   * Example: A = B;                                         *
//   ***********************************************************
MatrixSparse &MatrixSparse::operator =
                         (const MatrixSparse &rval)
   {
   CleanMatrix();
   if(numColumns != rval.numColumns)
      numColumns = rval.numColumns;
   if(numRows != rval.numRows)
      {
      delete elementRow;
      Initialize(rval.numRows,rval.numColumns);
      }
   Copy(rval);
   return *this;
   }

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

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

//   *************************< + >*****************************
//   * Purpose: Addition between matrices                      *
//   * Example: A = B + C;                                     *
//   ***********************************************************
MatrixSparse operator +
      (const MatrixSparse &lval,const MatrixSparse &rval)
   {
   if(lval.numRows != rval.numRows || 
      lval.numColumns != rval.numColumns)
        Error("%s%s%s+",
        MatrixSparse::ERROR,ERR_CHECK_DIMENSION,
        ERR_OPERATOR);
   MatrixSparse result = lval;
   result += rval;
   return result;
   }

//   *************************< += >****************************
//   * Purpose: Adding a matrix which is also on the right     *
//   * Example: A += B; // A = A + B;                          *
//   ***********************************************************
MatrixSparse &MatrixSparse::operator +=
             (const MatrixSparse &rval)
   {
   if(numRows != rval.numRows
             ||
   numColumns != rval.numColumns)
       Error("%s%s%s+=",MatrixSparse::ERROR,
       ERR_CHECK_DIMENSION,ERR_OPERATOR);
   Element *elem;
   for(int row = 1;row <= numRows;row++)
      {
      elem = rval.elementRow[row];
      for(int column = 1;column <= numColumns;column++)
        {
        if(elem == 0)break;
        if(elem->column == column)
           {
           (*this)(row,column) += elem->value;
           elem = elem->next;
           }
        }
      }
   return *this;
   }

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

//   *************************< - >*****************************
//   * Purpose: Difference between matrices                    *
//   * Example: A = B - C;                                     *
//   ***********************************************************
MatrixSparse operator -
       (const MatrixSparse &lval,const MatrixSparse &rval)
   {
   if(lval.numRows != rval.numRows || 
      lval.numColumns != rval.numColumns)
        Error("%s%s%s+",
        MatrixSparse::ERROR,
        ERR_CHECK_DIMENSION,ERR_OPERATOR);
   MatrixSparse result = lval;
   result -= rval;
   return result;
   }

//   *************************< -= >****************************
//   * Purpose: Difference from a matrix also on the right     *
//   * Example: A -= B; // A = A - B;                          *
//   ***********************************************************
MatrixSparse &MatrixSparse::operator -=
                 (const MatrixSparse &rval)
   {
   if(numRows != rval.numRows
             ||
   numColumns != rval.numColumns)
       Error("%s%s%s+=",MatrixSparse::ERROR,
       ERR_CHECK_DIMENSION,ERR_OPERATOR);
   Element *elem;
   for(int row = 1;row <= numRows;row++)
      {
      elem = rval.elementRow[row];
      for(int column = 1;column <= numColumns;column++)
        {
        if(elem == 0)break;
        if(elem->column == column)
           {
           (*this)(row,column) -= elem->value;
           elem = elem->next;
           }
        }
      }
   return *this;
   }

//   ===========================================================
//   **********************   Product   ************************
//   ===========================================================
//   ***********************< Product >*************************
//   * Purpose: Product between sparse matrices and vectors    *
//   * Example: Product(A,x,&y); y = A*x;                      *
//   ***********************************************************
void Product(const MatrixSparse &lval,const Vector &rval,
            Vector *result)
   {
   if(lval.numColumns != rval.Size())
       Error("%s%s%s*",
       MatrixSparse::ERROR,
       ERR_CHECK_DIMENSION,ERR_OPERATOR);
   if(rval.WhoAmI() == result->WhoAmI())
      {
      Vector aux;
      Product(lval,rval,&aux); // call recursive
      Swap(&aux,result); // avoid copying
      }
   else
      {
      ChangeDimensions(lval.numRows,result);
      for(int i=1;i <= lval.numRows;i++)
        (*result)[i] = 0.; // initialise
      Element *elem;
      for(int row = 1;row <= lval.numRows;row++)
        {
        elem = lval.elementRow[row];
        for(int column = 1;
               column <= lval.numColumns;column++)
           {
           if(elem == 0)break;
           if(elem->column == column)
              {
              float *temp = &(*result)[row];
              *temp = Control(*temp +
                  elem->value * rval.GetValue(column));
              elem = elem->next;
              }
           }
        }
      }
   }

//   *************************< * >*****************************
//   * Purpose: Product between sparse matrices and vectors    *
//   * Example: y = A * x;                                     *
//   ***********************************************************
Vector operator *
       (const MatrixSparse &lval,const Vector &rval)
   {
   Vector result;
   Product(lval,rval,&result);
   return result;
   }

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

//   ***********************< Print >***************************
//   * Purpose: Printing out the elements of a matrix          *
//   * Description: Uses the file bzzFileOut global            *
//   *              in utility.hpp                             *
//   * Example: A.Print("Matrix A");                           *
//   ***********************************************************
void MatrixSparse::Print(char *msg)
   {
   Element *elem;
   if(*msg)fprintf(bzzFileOut,"\n%s",msg);
   fprintf(bzzFileOut,"\nrows %d columns %d\n",
          numRows,numColumns);
   for(int row = 1;row <= numRows;row++)
      {
      elem = elementRow[row];
      for(int column = 1;column <= numColumns;column++)
        {
        if(elem != 0 && elem->column == column)
           {
           fprintf(bzzFileOut," %12.4e",elem->value);
           elem=elem->next;
           }
        else
           fprintf(bzzFileOut,"       *     ");
        }
      fprintf(bzzFileOut,"\n");
      }
   }

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

//   ******************< DeleteElement >************************
//   * Purpose: Eliminating an element from a matrix           *
//   * Example: A.DeleteElement(i,j);                          *
//   ***********************************************************
void MatrixSparse::DeleteElement(int row,int column)
   {
   if(elementRow == 0)return;
   if(row < 1 || row > numRows)return;
   if(column < 1 || column > numColumns)return;
   Element *elem,*temp;
   elem = elementRow[row];
   if(elem == 0)return;
   if(column == elem->column)
      {
      elementRow[row] = elem->next;
      delete elem;
      return;
      }
   for(;elem->next != 0; elem = elem->next)
      {
      if(column == elem->next->column)
        {
        temp = elem->next;
        elem->next = elem->next->next;
        delete temp;
        return;
        }
      }
   }

//   ********************< DeleteRow >**************************
//   * Purpose: Eliminating a row from a matrix                *
//   * Example: A.DeleteRow(i);                                *
//   ***********************************************************
void MatrixSparse::DeleteRow(int row)
   {
   if(row < 1 || row > numRows)return;
   if(elementRow == 0)return;
   Element *elem,*temp;
   elem = elementRow[row];
   while(elem != 0)
      {
      temp = elem;
      elem = elem->next;
      delete temp;
      }
   elementRow[row] = 0;
   }

//   *******************< DeleteMatrix >************************
//   * Purpose: Eliminating a matrix of no use                 *
//   * Example: A.DeleteMatrix();                              *
//   ***********************************************************
void MatrixSparse::DeleteMatrix(void)
   {
   if(elementRow == 0)return;
   CleanMatrix();
   delete elementRow;
   elementRow = 0;
   numRows = numColumns = 0;
   }

//   *******************< CleanMatrix >*************************
//   * Purpose: Eliminating all coefficients of a matrix       *
//   * Example: A.CleanMatrix();                               *
//   ***********************************************************
void MatrixSparse::CleanMatrix(void)
   {
   if(elementRow == 0)return;
   Element *elem,*temp;
   for(int row =1; row <= numRows; row++)
      {
      elem = elementRow[row];
      while(elem != 0)
        {
        temp = elem;
        elem = elem->next;
        delete temp;
        }
      elementRow[row] = 0;
      }
   }

//   *******************< CleanMatrix >*************************
//   * Purpose: Eliminating the coefficients of a matrix which *
//   *          are smaller than eps                           *
//   * Example: A.CleanMatrix(1.e-20);                         *
//   ***********************************************************
void MatrixSparse::CleanMatrix(float eps)
   {
   if(elementRow == 0)return;
   eps = Abs(eps);
   Element *elem,*temp;
   for(int row =1; row <= numRows; row++)
      {
      elem = elementRow[row];
      for(int column = 1;elem != 0;column++)
        {
        if(elem->column == column)
           {
           if(Abs(elem->value) < eps)
              DeleteElement(row,column);
           elem=elem->next;
           }
        }
      }
   }

//   ===========================================================
//   ==========   System solving functions =====================
//   ===========================================================
//   ********************< SolveRight >*************************
//   * Purpose: Solving a system with a Right MatrixSparse     *
//   * Example: A.SolveRight(&bx);                             *
//   ***********************************************************
void MatrixSparse::SolveRight(Vector *bx)
   {
   if(numRows != bx->dimensions)
      Error("%s%s%sSolve",
      MatrixSparse::ERROR,
      ERR_CHECK_DIMENSION,ERR_FUNCTION);
   Element *elem;
   double sum;
   float temp,*x = bx->vector;
   for(int row = numRows;row >= 1;row--)
      {
      elem = elementRow[row];
      if(elem != 0 && elem->column < row)
        {
        Error("%sMatrix non Right",
               MatrixSparse::ERROR);
        }
      if(elem == 0 || (temp = elem->value) == 0.
                 || elem->column != row)
        {
        Message("%sSingular Matrix",
               MatrixSparse::ERROR);
        x[row]=0.;
        }
      else
        {
        sum = x[row];
        elem = elem->next;
        for(int column = row+1;
               column <= numColumns;column++)
           {
           if(elem != 0 && elem->column == column)
              {
              sum -= elem->value * x[column];
              elem = elem->next;
              }
           }
        x[row] = Control(sum/temp);
        }
      }
   }

//   ********************< SolveLeft >**************************
//   * Purpose: Solving a system with a Left MatrixSparse      *
//   * Example: A.SolveLeft(&bx);                              *
//   ***********************************************************
void MatrixSparse::SolveLeft(Vector *bx)
   {
   if(numRows != bx->dimensions)
      Error("%s%s%sSolve",
      MatrixSparse::ERROR,
      ERR_CHECK_DIMENSION,ERR_FUNCTION);
   Element *elem;
   double sum;
   float temp,*x = bx->vector;
   for(int row = 1;row <= numRows;row++)
      {
      elem = elementRow[row];
      if(elem == 0)
        {
        Message("%sSingular Matrix",
               MatrixSparse::ERROR);
        x[row]=0.;
        }
      else
        {
        sum = x[row];
        for(int column = 1;column < row;column++)
           {
           if(elem != 0
                  &&
           elem->column == column
                  && elem->column != row)
              {
              sum -= elem->value * x[column];
              elem = elem->next;
              }
           }
        if(elem->value == 0. || elem->column != row)
           {
           Message("%sSingular Matrix",
                  MatrixSparse::ERROR);
           x[row]=0.;
           }
         else
           x[row] = Control(sum/elem->value);
        }
      }
   }