//   ================================================================
//   =============   Examples of MatrixSparse   =====================
//   * Description: Chapter 13                                      *
//   * Scientific C++ Building Numerical Libraries                  *
//   *                the Object-Oriented Way                       *
//   * by G. Buzzi-Ferraris                                         *
//   * Addison-Wesley (1993)                                        *
//   ****************************************************************
//   * To be compiled with:					    *
//   * UTILITY.OBJ, VECTOR.OBJ, MATRIX.OBJ, SPARSE.OBJ	            *
//   ****************************************************************
//   * To execute: exsparse <number> <file name> <Y/N>		    *
//   * If <number> equals zero it executes all examples		    *
//   * otherwise it executes only that example			    *
//   * It writes the results on < file name>			    *
//   * If not used it writes on screen				    *
//   * If one puts Y/N equal to Y it prints messages as well	    *
//   * The default is Y						    *
//   ****************************************************************
//   * Use: exsparse 0 pro.doc					    *
//   * Description:	Executes all examples			    *
//   * 	Prints the result on pro.doc			            *
//   * Use:	exsparse 2 pro.doc N				    *
//   * Description:	Executes example 2			    *
//   * 	Prints the result on pro.doc			            *
//   * 	Does not print from Message			            *
//   * Use:	exsparse					    *
//   * Description:	Executes all examples			    *
//   * 	Prints the result on stdout			            *
//   ================================================================

//   ================================================================
//   *                         PROTOTYPES                           *
//   ================================================================
void Ex_MatrixSparse1(void);
void Ex_MatrixSparse2(void);

//   ================================================================
//   *                       include				    *
//   ================================================================
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <ctype.h>

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

#define NEXMATRIX_SPARSE 2
#define NMAX NEXMATRIX_SPARSE+1

main(int argc,char **argv)
{
int i;
void (*ExMatrixSparse[NMAX])(void);
char myfilename[81];
ExMatrixSparse[1]=Ex_MatrixSparse1;
ExMatrixSparse[2]=Ex_MatrixSparse2;

if(argc>=3)
   {
   strcpy(myfilename,argv[2]);

   if((bzzFileOut=fopen(myfilename,"w"))==NULL)
      Error("fopen Failed in EXMATRIX_SPARSE");
      if(argc==4)bzzYesNo=toupper(argv[3][0]);
   }
else
   {
   bzzFileOut=stdout;
   }

if(argc==1||atoi(argv[1])<=0||atoi(argv[1])>NEXMATRIX_SPARSE)
   {     // all examples
   for(i=1;i<=NEXMATRIX_SPARSE;i++)
   (*ExMatrixSparse[i])();
   }
else
   {
   (*ExMatrixSparse[atoi(argv[1])])();
   }
return 1;
}

//   ================================================================
//   ====================   EXAMPLES   ==============================
//   ================================================================

void Ex_MatrixSparse1(void)
{
Message("\n\n*******  MatrixSparse example.**************");
    {
    MatrixSparse A(5,5),B(5,5);

    A(2,2)=22.;	A(3,3)=33.;
    A(4,4)=44.;	A(1,1)=11.;
    A(1,5)=15.;	A(1,3)=13.;

    B(2,2)=22.;	B(3,3)=33.;
    B(3,5)=35.;	B(5,2)=52.;
    B(5,5)=55.;

    A.Print("A 2,2 3,3 4,4 1,1 1,5 1,3");
    B.Print("A 2,2 3,3 4,4 1,1 1,5 1,3");

Message("\n ******* Copy-Inizializer **************");
   MatrixSparse C=A+B;
   C.Print("MatrixSparse C = A + B;");

Message("\n ******* operator = **************");
   MatrixSparse D;
   D = A + B - A;
   D.Print("D = A + B - A");

Message("\n ******* CleanMatrix **************");
   D.CleanMatrix(1.e-20);
   D.Print("D.CleanMatrix(1.e-20)");

Message("\n ******* CleanMatrix **************");
   D.CleanMatrix();
   D.Print("D.CleanMatrix()");

Message("\n ******* operator += **************");
   A += B;
   A.Print("A += B");

Message("\n ******* operator -= **************");
   A -= B;
   A.Print("A -= B");

Message("\n ******* DeleteElement 1,3 1,5 3,3 **********");
   C.DeleteElement(1,3);
   C.DeleteElement(1,5);
   C.DeleteElement(3,3);
   C.Print("DeleteElement");

Message("\n ******* DeleteRow 5 **************");
   C.DeleteRow(5);
   C.Print("C without row 5");

Message("\n ******* DeleteMatrix **************");
   C.DeleteMatrix();
   C.Print("void");
   }
   {
   MatrixSparse A(7,5);
   Vector x(5);
   Message("\n ******* y = A*x **************");
   A(1,1)=1.; A(1,4)=3.; A(1,5)=3.;
   A(2,3)=5.; A(2,4)=1.; A(2,5)=2.;
   A(4,1)=4.; A(4,2)=4.; A(4,5)=1.;
   A(7,1)=2.; A(7,4)=1.; A(7,5)=3.;
   A.Print("A");
   x(1) = 1.; x(2) = 2.; x(3) =3.; x(4) = 4.; x(5) =5.;
   Vector y = A*x;
   y.Print("y = A*x");
   Product(A,x,&y);
   y.Print("Product(A,x,&y)");
   }
}
void Ex_MatrixSparse2(void)
{
  Message("\n\n*******  Solve Sparse example.***********");
  Message("\n\n*******  SolveRight example**************");
    {
    MatrixSparse A(1000,1000);

    A(1,1)=1.;	A(1,2)=2.;
    A(2,2)=2.;	A(2,3)=3.;
    A(3,3)=5.;
    for(int i = 4;i <= 1000;i++)A(i,i)=3.;
    for(i = 4;i <= 900;i++)A(i,i+2)=1.;
    Vector bx(1000);
    bx(1) = 3.; bx(2) = 5.; bx(3) = 5.;
    for(i = 4;i <= 1000;i++)bx(i) = 3.;
    for(i = 4;i <= 900;i++)bx(i) += 1.;
    A.SolveRight(&bx);
    bx.Print("bx SolveRight");
    }
Message("\n\n******* SolveLeft example**************");
    {
    MatrixSparse A(4,4);

    A(1,1)=1.;
    A(2,1)=1.;	A(2,2)=1.;
			  A(3,3)=1.;
    A(4,1)=1.;	A(4,2)=1.;          A(4,4)=1.;

    Vector bx(4);
    bx(1)=1.;bx(2)=2.;bx(3)=1.;bx(4)=3.;
    A.SolveLeft(&bx);
    bx.Print("LEFT");
   }
}