//   =====================< EXLEFT.CPP >=========================
//   * Examples of MatrixLeft 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, symm.obj					*
//   ************************************************************
//   * To execute: exleft <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 stdout				*
//   * If one puts Y/N equal to Y it prints messages as well	*
//   * The default is Y						*
//   ************************************************************
//   * exleft 0 pro.doc				                *
//   * Executes all examples Prints the result on PRO.DOC	*
//   * exleft 1 pro.doc						*
//   * Executes example 1 with messages and results on PRO.DOC	*
//   * exleft							*
//   * Executes everything and prints it to stdout	 	*
//   ============================================================

// Prototype examples
void Ex_MatrixLeft_1(void);
void Ex_MatrixLeft_2(void);
void Ex_MatrixLeft_3(void);
void Ex_MatrixLeft_4(void);
void Ex_MatrixLeft_5(void);
void Ex_MatrixLeft_6(void);
void Ex_MatrixLeft_7(void);


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


int main(int argc,char **argv)
{
const int NEXMATRIXLEFT = 7;
const int NMAX = NEXMATRIXLEFT+1;

int i;
void (*ExMatrixLeft[NMAX])(void);
char myfilename[81];
ExMatrixLeft[1]=Ex_MatrixLeft_1;
ExMatrixLeft[2]=Ex_MatrixLeft_2;
ExMatrixLeft[3]=Ex_MatrixLeft_3;
ExMatrixLeft[4]=Ex_MatrixLeft_4;
ExMatrixLeft[5]=Ex_MatrixLeft_5;
ExMatrixLeft[6]=Ex_MatrixLeft_6;
ExMatrixLeft[7]=Ex_MatrixLeft_7;

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

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

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

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

//   ===========================================================
//   ===============   INITIALISATIONS   =======================
//   ===========================================================
void Ex_MatrixLeft_1(void)
{
Message("\n***** MatrixLeft Initialisations **********");
   {
   MatrixLeft A;
	A.Print("A left without nothing");
   MatrixLeft B(2,2);
	B.SetValue(1,1,10.);
	B.SetValue(2,1,5.);	B.SetValue(2,2,15.);
	fprintf(bzzFileOut,"(2,1)=5. %f",B.GetValue(2,1));
	B.Print("MLeft B(2,2) (2,1)=5.");
   MatrixLeft C(3,3,
		1.,
		2.,3.,
		4.,5.,6.);
	C.Print("MatrixLeft C(3,3,1.,...)");
   float w[]={
	     1.,
	     2.,3.,
	     4.,5.,6.
	     };
   MatrixLeft D(3,3,w);	D.Print("MatrixLeft D(3,w)");
   MatrixLeft E=D;		E.Print("MatrixLeft E=D");
   MatrixLeft F("LEFT.DAT");	F.Print("From File");
   }
}

//   ============================================================
//   ==================   OPERATIONS   ==========================
//   ============================================================
void Ex_MatrixLeft_2(void)
{
   {
   Message("\n*************< Sum >*********************");
   MatrixLeft A(3,3,
		 1.,
		 2.,3.,
		 4.,5.,6.);
   MatrixLeft B(3,3,
		9.,
		8.,7.,
		6.,5.,4.);
   MatrixLeft C;
   Sum(A,B,&C);
   C.Print();
   }
   {
   Message("\n*************< Sum >*********************");
   MatrixLeft A(3,3,
		 1.,
		 2.,3.,
		 4.,5.,6.);
   MatrixLeft B(3,3,
		9.,
		8.,7.,
		6.,5.,4.);
   MatrixLeft C = A + B;
   C.Print();
   }
   {
   Message("\n*************< Sum >*********************");
   MatrixLeft A(3,3,
		 1.,
		 2.,3.,
		 4.,5.,6.);
   MatrixLeft B(3,3,
		 9.,
		 8.,7.,
		 6.,5.,4.);
   Sum(&A,B);
   A.Print();
   }
   {
   Message("\n*************< Sum >*********************");
   MatrixLeft A(3,3,
		 1.,
		 2.,3.,
		 4.,5.,6.);
   MatrixLeft B(3,3,
		 9.,
		 8.,7.,
		 6.,5.,4.);
   Sum(A,&B);
   B.Print();
   }
   {
   Message("\n*************< Sum >*********************");
   MatrixLeft A(3,3,
		 1.,
		 2.,3.,
		 4.,5.,6.);
   Sum(&A);
   A.Print();
   }
   {
   Message("\n*************< Sum >*********************");
   MatrixLeft A(3,3,
		 1.,
		 2.,3.,
		 4.,5.,6.);
   Sum(A,A,&A);
   A.Print();
   }
   {
   Message("\n******** MatrixLeft  Operators ******");
   float w[]={
	     1.,
	     2.,3.,
	     4.,5.,6.
	     };
   float q[]={
	     1.,0.,0.,
	     2.,3.,0.,
	     4.,5.,6.
	     };
   float z[]={
	      1.,0.,0.,
	      2.,3.,0.,
	      4.,5.,6.
	      };
   Matrix Q(3,3,q),T(3,3,z),U(3,3);
   MatrixLeft A(3,3,w),B(3,3,1.,2.,3.,4.,5.,6.);
   A.Print("A"); B.Print("B");
   MatrixLeft C=A+B;		C.Print("A+B");
   C=A-B;			C.Print("A-B");
   C=A*B;			C.Print("A*B");
   U=Q*T;                       U.Print("A*B");
   C=3.*A;                     	C.Print("3.*A");
   C=A/3.;                     	C.Print("A/3.");
   A*=3.;                       A.Print("A*=3.");
   A/=3.;                       A.Print("A/=3.");
   float e[]={1.,2.,3.};
   Vector b(3,e);
   Vector c=A*b;                c.Print("A*b");
   if(A == A)fprintf(bzzFileOut,"\nA == A");
   if(C != A)fprintf(bzzFileOut,"\nC != A");
   }
   {
   Message("\n*** TProduct *");
   MatrixLeft L(4,4,
	       1.,
	       2.,3.,
	       4.,5.,6.,
	       7.,8.,9.,10.);
   Matrix A(4,4,
	    1.,0.,0.,0.,
	    2.,3.,0.,0.,
	    4.,5.,6.,0.,
	    7.,8.,9.,10.);
   MatrixSymm S;
   TProduct(L,&S);
   S.Print("S");
   TProduct(&A);
   A.Print("verification with Matrix");
   }
   {
   Message("\n*** ProductT *");
   MatrixLeft L(4,4,
	       1.,
	       2.,3.,
	       4.,5.,6.,
	       7.,8.,9.,10.);
   Matrix A(4,4,
	    1.,0.,0.,0.,
	    2.,3.,0.,0.,
	    4.,5.,6.,0.,
	    7.,8.,9.,10.);
   MatrixSymm S;
   ProductT(L,&S);
   S.Print("S");
   ProductT(&A);
   A.Print("verification with Matrix");
   }
}

//   ===========================================================
//   ====================   NORMS   ============================
//   ===========================================================
void Ex_MatrixLeft_3(void)
{
    {
    Message("\n********** Matrix Norms **************");
    float w[]={
	       1.,0.,0.,
	       2.,3.,0.,
	       4.,5.,6.
	       };
    float q[]={1.,2.,3.,4.,5.,6.};
    Matrix A(3,3,w);		A.Print("A");
    MatrixLeft B(3,3,q);	B.Print("B");
    fprintf(bzzFileOut,"\nNormT %f",A.NormT());
    fprintf(bzzFileOut,"\nNormT Left %f",B.NormT());

    fprintf(bzzFileOut,"\nNormR %f",A.NormR());
    fprintf(bzzFileOut,"\nNormR Left %f",B.NormR());

    fprintf(bzzFileOut,"\nNormC %f",A.NormC());
    fprintf(bzzFileOut,"\nNormC Left %f",B.NormC());

    fprintf(bzzFileOut,"\nNormF %f",A.NormF());
    fprintf(bzzFileOut,"\nNormF Left %f",B.NormF());

    fprintf(bzzFileOut,"\nNormI %f",A.NormI());
    fprintf(bzzFileOut,"\nNormI Left %f",B.NormI());

    fprintf(bzzFileOut,"\nNorm1 %f",A.Norm1());
    fprintf(bzzFileOut,"\nNorm1 Left %f",B.Norm1());
    }
}

//   ===========================================================
//   ===============   Save and Recover   ======================
//   ===========================================================
void Ex_MatrixLeft_4(void)
{
   {
   Message("\n** WARNING: Recover and Save  uses drive C:\\TEMP **");
   TempFile f1("C:\\TEMP\\");
   MatrixLeft A(3,3,
		 1.,
		 2.,3.,
		 4.,5.,6.);
   A.Print();
   A.Save('*',f1.FileName());
   A(3,3)=100.; // modify A
   Recover(&A,'*',f1.FileName()); // recover previous A
   MatrixLeft B('*',f1.FileName()); // initialise B
   B.Print("B");
   A.Print("A");
   MatrixLeft C;
   Recover(&C,'*',f1.FileName()); // C is resized
   C.Print("C");
   }
}

//   ===========================================================
//   ===================   INVERSE  ============================
//   ===========================================================
void Ex_MatrixLeft_5(void)
{
    {
    Message("\n********** Inverse **************");
   MatrixLeft A(4,4,
		1.,
		2.,3.,
		4.,5.,6.,
		7.,8.,9.,10.);
   MatrixLeft B = A;
   Inverse(&B);
   B.Print("inverse");
   (A*B).Print("control");
    }
}

//   ============================================================
//   ===================   LLT LTL   ============================
//   ============================================================
void Ex_MatrixLeft_6(void)
{
   {
   Message("\n********** LLT and LTL **************");
   Matrix A(3,3,
	     11.,0.,0.,
	     4.1,55.,0.,
	     3.2,2.5,7.8);
   MatrixLeft L(3,3,
		 11.,
		 4.1,55.,
		 3.2,2.5,7.8);
   Matrix C;
   TProduct(A,A,&C);
   C.Print("ATA");
   MatrixSymm S;
   TProduct(L,&S);
   S.Print("LTL");
   ProductT(&A);
   A.Print("AAT");
   ProductT(L,&S);
   S.Print("LLT");
   }
}

//   ============================================================
//   ============   LEFT SYSTEM SOLUTIONS  ======================
//   ============================================================
void Ex_MatrixLeft_7(void)
{
   {
Message("\n*** Solution,condition,determinant\n");
   MatrixLeft L(3,3,
		1.,
		2.,3.,
		4.,5.,6.);
   Vector bx(3,1.,5.,15.);
   fprintf(bzzFileOut,"\n%f",L.ConditionNumber());
   fprintf(bzzFileOut,"\n%e",L.Determinant());
   Solve(L,&bx);
   bx.Print("Solution");
   }
   {
   Message("\n\Left system solution  ");
   MatrixLeft L(4,4,
	       1.,
	       2.,3.,
	       4.,5.,6.,
	       7.,8.,9.,10.);
   Vector b(4,1.,5.,15.,34.),x;
   Solve(L,b,&x);
   b.Print("Known terms");
   x.Print("x");
   }
   {
   Message("\n*** Solution BX\n");
   MatrixLeft L(3,3,
		 1.,
		 2.,3.,
		 4.,5.,6.);
   Matrix B(3,5,
	     1.,  1.,   2.,  1.,   5.,
	     5.,  5.,  10.,  5.,  13.,
	    15., 33.,  30., 15.,  31.);
   Solve(L,&B);
   B.Print();
   }
   {
   MatrixLeft L(4,4,
	       1.,
	       2.,3.,
	       4.,5.,6.,
	       7.,8.,9.,10.);
   Matrix B(4,2,1.,2.,
		5.,  10.,
		15., 30.,
		34., 68.),X;
   Solve(L,B,&X);
   B.Print("Known terms");
   X.Print("X");
   }
}