//   =====================< EXMATR1.CPP >===========================
//   * Examples of Matrix class 				   *
//   * Description: Chapter 12                                     *
//   * Scientific C++ Building Numerical Libraries                 *
//   *                the Object-Oriented Way                      *
//   * by G. Buzzi-Ferraris                                        *
//   * Addison-Wesley (1993)                                       *
//   ***************************************************************
//   * To be linked with: utility.obj, vector.obj, matrix.obj,	   *
//   * left.obj, right.obj, symm.obj				   *
//   ***************************************************************
//   * To execute: exmatr1 <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						   *
//   ***************************************************************
//   * exmatr1 0 pro.doc				           *
//   * Executes all examples prints the result on pro.doc	   *
//   * exmatr1 1 pro.doc					   *
//   * Executes example 1 with messages and results on PRO.DOC	   *
//   * exmatr1							   *
//   * Executes everything and sends it to stdout		   *
//   ===============================================================


// Prototype examples
void Ex_Matrix_1(void);
void Ex_Matrix_2(void);
void Ex_Matrix_3(void);
void Ex_Matrix_4(void);
void Ex_Matrix_5(void);
void Ex_Matrix_6(void);
void Ex_Matrix_7(void);
void Ex_Matrix_8(void);

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <ctype.h>
#include "utility.hpp"
#include "vector.hpp"
#include "matrix.hpp"
#include "right.hpp"
#include "left.hpp"
#include "symm.hpp"


int main(int argc,char **argv)
{
const int NEXMATRIX = 8;
const int NMAX = NEXMATRIX+1;
int i;
void (*ExMatrix[NMAX])(void);
char myfilename[81];

ExMatrix[1]=Ex_Matrix_1;
ExMatrix[2]=Ex_Matrix_2;
ExMatrix[3]=Ex_Matrix_3;
ExMatrix[4]=Ex_Matrix_4;
ExMatrix[5]=Ex_Matrix_5;
ExMatrix[6]=Ex_Matrix_6;
ExMatrix[7]=Ex_Matrix_7;
ExMatrix[8]=Ex_Matrix_8;

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

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

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

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

void Ex_Matrix_1(void)
{
//   ===========================================================
//   ===============   INITIALISATIONS   =======================
//   ===========================================================

   Message("\n\n\n* Examples of Matrix initialisations *");
    {
    Matrix A; 		A.Print("A without anything");
    Matrix B(3,2);
		B(1,1)=11.;
		B(2,1)=21.; B(2,2)=22.;
			    B(3,2)=32.;
    B.Print("B(3,2) with assignment to coefficients");
    ChangeDimensions(2,5,&B);
    B.Print("ChangeDimensions(2,5,&B)");
    Matrix C(2,2,
	     1.,2.,
	     3.,4.);
	C.Print("C(2,2,1.,...)");
    A = B = C;
	A.Print("Multiple assignment A=B=C");
    float w[] = {1.,2.,3.,4.,5.,6.};
    Matrix D(2,3,w);	 	D.Print("D(2,3,w)");
    char *fileName="MATR.DAT";
    Matrix E(fileName);		E.Print("From File");
    Vector v(6,w);
    Matrix F(v);		F.Print("F(Vector v)");

    MatrixRight G(3,3,w);
    Matrix H(G);		H.Print("From MatrixRight");
    MatrixLeft I(3,3,w);
    Matrix L(I);		L.Print("From MatrixLeft");
    MatrixSymm M(3,3,w);
    Matrix N(M);		N.Print("From MatrixSymm");

    }
}

void Ex_Matrix_2(void)
{
//   ============================================================
//   ====================   NORMS   =============================
//   ============================================================
Message("\n********** Matrix Norms **************");
    {
    float w[]={1.,-2.,-3.,4.};
    Matrix A(2,2,w);
    Message("\nNormT %f",A.NormT());
    Message("\nNormR %f",A.NormR());
    Message("\nNormC %f",A.NormC());
    Message("\nNormF %f",A.NormF());
    Message("\nNormI %f",A.NormI());
    Message("\nNorm1 %f",A.Norm1());
    }
}

void Ex_Matrix_3(void)
{
//   ===========================================================
//   ===================   MAX-MIN   ===========================
//   ===========================================================
   {
   Message("\n\n********* Max Min **********");
   Matrix A(3,5,
	   0.,-10.,1.,30.,50.,
	   -20.,20.,1.,2.,12.,
	   7.,-120.,45.,22.,-1000.);
   A.Print("Matrix");
   float xf;
   int i,j;
   xf=A.Max();            Message("\nMax %f",xf);
   xf=A.Max(&i,&j);       Message("\nMax %d %d %f",i,j,xf);
   xf=A.MaxAbs();         Message("\nMaxAbs %f",xf);
   xf=A.MaxAbs(&i,&j);    Message("\nMaxAbs %d %d %f",i,j,xf);
   xf=A.Min();            Message("\nMin %f",xf);
   xf=A.Min(&i,&j);       Message("\nMin %d %d %f",i,j,xf);
   xf=A.MinAbs();         Message("\nMinAbs %f",xf);
   xf=A.MinAbs(&i,&j);    Message("\nMinAbs %d %d %f",i,j,xf);
   }
}
void Ex_Matrix_4(void)
{
//   ===========================================================
//   ================    SUBMATRICES  ==========================
//   ===========================================================

Message("\nExamples of construction of matrices from submatrices"
	"\n and vice versa");
     {
     Matrix A(2,2,3.,3.,3.,3.),
	    B(2,3,5.,5.,5.,5.,5.,5.),
	    C(5,2,8.,8.,8.,8.,8.,8.,8.,8.,8.,8.),
	    D(2,3,1.,1.,1.,1.,1.,1.),
	    E(3,3,9.,9.,9.,9.,9.,9.,9.,9.,9.);
     Matrix X;
     X=(A||B)&&(C||(D&&E));
     X.Print("X=(A||B)&&(C||(D&&E))");
     Matrix Y(2,3,3,3,X);
     Y.Print("\nY should be identical a d");

//   ============================================================
//   ================   ASSIGNMENTS   ===========================
//   ============================================================

Message("\n**** Assignment examples ********************");
     X(1,2)=12.;
	Message("\nX(1,2)=12.    %f",X(1,2));
     X[1][2]=15.;
	Message("\nX[1][2]=15.    %f",X[1][2]);
     X.SetValue(2,3,23.);
       Message("\nSetGet 23.    %f",X.GetValue(2,3));
     Vector b,c;
     b=X.GetRow(3);  // reads row 3
	b.Print("Row 3 of X");
     c=X.GetColumn(1); //reads column 1
	c.Print("column 1 of X");

     b(4)=333.; //changes element 4 of the vector
     X.SetRow(2,b);	X.Print("X with row 2 changed");
     Message("\nX 2,3=%f 3,4=%f",X(2,3),X.GetValue(3,4));
     }
}

void Ex_Matrix_5(void)
{
//   ===========================================================
//   ===============   Save and Recover   ======================
//   ===========================================================

Message("\n****** save and read ******");
   {
   Message("\n** WARNING: Recover and Save  uses drive C:\\TEMP **");
   TempFile f1("C:\\TEMP\\"),f2("C:\\TEMP\\");
   Matrix A(3,3,1.,2.,3.,4.,5.,6.,7.,8.,9.);
   A.Save(f2.FileName());
   Matrix B(f2.FileName());
   B.Print("B from Formatted File ");
   A.Save('*',f1.FileName());
   Matrix C('*',f1.FileName());
   C.Print("C from unformatted File ");
   }

Message("\n****** save and recover ******");
   {
   TempFile f1("C:\\TEMP\\"),f2("C:\\TEMP\\");
   Matrix A(10,10);
   int i,j;
   for(i=1; i <= 10;i++)
      for(j=1;j <=10;j++)
      A(i,j) = i+j;
   A.Save(f2.FileName());
   A(10,10)=5.;
   Recover(&A,f2.FileName());
   Message("\n%f",A(10,10));
   A.Save('*',f1.FileName());
   A(10,10)=5.;
   Recover(&A,'*',f1.FileName());
   Message("\n%f",A(10,10));
   }
}

void Ex_Matrix_6(void)
{
//   ===========================================================
//   ====================   Transpose   ========================
//   ===========================================================
Message("\n******** Modifying functions ********");
     {
     float a[]={1.,2.,3.,4.,5.,6.,7.,8.,9.,10.};
     Matrix A(2,5,a);	A.Print("A");
     Transpose(&A);	A.Print("A Transpose");
     Transpose(&A);	A.Print("A Transpose");
     Vector w(5,9.,9.,9.,9.,9.);	Matrix B=w;
     Transpose(&B);
     Matrix C=A&&B;	C.Print("C=A && Vector 1");
     A.Print("A");
     }
}

void Ex_Matrix_7(void)
{
//   ==========================================================
//   ======================   Swap   ==========================
//   ==========================================================
Message("\n******** Swap Matrices ********");
   {
   float wa[]={1.,2.,3.,4.,5.,6.,7.,8.};
   float wb[]={100.,200.,300.,400.,500.,600.};
   Matrix A(2,4,wa);
      {
      Matrix B(3,2,wb);
	A.Print("A before Swap");
	B.Print("B before Swap");
      Swap(&A,&B);
	A.Print("A after Swap");
	B.Print("B after Swap");
      }
   A.Print("A");
   }
}

void Ex_Matrix_8(void)
{
//   ============================================================
//   ====================   SumRankOne   ========================
//   ============================================================
   {
   Message("\n***********< SumRankOne >*************\n");
   float w[]=
      {1.,2.,3.,
       4.,5.,6.,
       7.,8.,9.};
   Matrix A(3.,3.,w);
   float x[]=
      {1.,2.,3.};
   Vector u(3,x);
   float y[]=
      {1.,2.,3.};
   Vector vT(3,y);
   SumRankOne(u,vT,&A);
   A.Print("SumRankOne(u,vT,&A)");
   }
   {
   Message("\n***********< SumRankOne >*************\n");
   float w[]=
      {1.,2.,3.,
       4.,5.,6.,
       7.,8.,9.};
   Matrix A(3.,3.,w);
   float x[]=
      {.5,1.,1.5};
   Vector u(3,x);
   float y[]=
      {1.,2.,3.};
   Vector vT(3,y);
   SumRankOne(2.,u,vT,&A);
   A.Print("SumRankOne(2.,u,vT,&A)");
   }
   {
   Message("\n***********< SumRankOne >*************\n");
   float w[]=
      {1.,2.,3.,
       4.,5.,6.,
       7.,8.,9.};
   Matrix A(3.,3.,w);
   float x[]=
      {3.,6.,9.};
   Vector u(3,x);
   float y[]=
      {1.,2.,3.};
   Vector vT(3,y);
   SumRankOne(u,vT,3.,&A);
   A.Print("SumRankOne(u,vT,3.,&A)");
   }

}