//   =====================< EXRIGHT.CPP >============================
//   * Examples of MatrixRight 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,	    *
//   * right.obj, left.obj, symm.obj			            *
//   ****************************************************************
//   * To execute: exright <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					            *
//   ****************************************************************
//   * exright 1 pro.doc					    *
//   * Executes example 1 with messages and results on PRO.DOC	    *
//   * exright						            *
//   * Executes everything and writes it on stdout		    *
//   ================================================================

// Prototype examples
void Ex_MatrixRight_1(void);
void Ex_MatrixRight_2(void);
void Ex_MatrixRight_3(void);
void Ex_MatrixRight_4(void);
void Ex_MatrixRight_5(void);
void Ex_MatrixRight_6(void);
void Ex_MatrixRight_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 "right.hpp"
#include "symm.hpp"

#define NEXMATRIXRIGHT 7
#define NMAX NEXMATRIXRIGHT+1

main(int argc,char **argv)
{
int i;
void (*ExMatrixRight[NMAX])(void);
char myfilename[81];

ExMatrixRight[1] = Ex_MatrixRight_1;
ExMatrixRight[2] = Ex_MatrixRight_2;
ExMatrixRight[3] = Ex_MatrixRight_3;
ExMatrixRight[4] = Ex_MatrixRight_4;
ExMatrixRight[5] = Ex_MatrixRight_5;
ExMatrixRight[6] = Ex_MatrixRight_6;
ExMatrixRight[7] = Ex_MatrixRight_7;

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

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

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

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

//   ================================================================
//   ===============   INITIALISATIONS   ============================
//   ================================================================
void Ex_MatrixRight_1(void)
{
Message("\n***** MatrixRight Initialisations **********");
   {
   MatrixRight A;
	A.Print("A right without anything");
   MatrixRight B(2,2);
	B.SetValue(1,1,10.);	B.SetValue(1,2,15.);
				B.SetValue(2,2,515.);
	fprintf(bzzFileOut,"\n(1,2)=15 %f",B.GetValue(1,2));
	B.Print("MRight B(2,2) (1,2)=15.");
   MatrixRight C(3,3,
		 1.,2.,3.,
		    4.,5.,
		       6.);
	C.Print("MatrixRight C(3,3,1.,...)");
   float w[]={
	     1.,2.,3.,
		4.,5.,
		   6.
	     };
   MatrixRight D(3,3,w);     D.Print("MatrixRight D(3,w)");
   MatrixRight E=D;		E.Print("MatrixRight E=D");
   MatrixRight F("RIGHT.DAT");	F.Print("From File");
   }
}

//   ================================================================
//   ==================   OPERATIONS   ==============================
//   ================================================================
void Ex_MatrixRight_2(void)
{
   {
   Message("\n*************< Sum >*********************");
   MatrixRight A(3,3,
		 1.,2.,3.,
		    4.,5.,
		       6.);
   MatrixRight B(3,3,
		9.,8.,7.,
		   6.,5.,
		      4.);
   MatrixRight C;
   Sum(A,B,&C);		 
   C.Print();
   }
   {
   Message("\n*************< Sum >*********************");
   MatrixRight A(3,3,
		 1.,2.,3.,
		    4.,5.,
		       6.);
   MatrixRight B(3,3,
		9.,8.,7.,
		   6.,5.,
		      4.);
   MatrixRight C = A + B;
   C.Print();
   }
   {
   Message("\n*************< Sum >*********************");
   MatrixRight A(3,3,
		 1.,1.,1.,
		    1.,1.,
		       1.);
   MatrixLeft B(3,3,
		1.,
		1.,1.,
		1.,1.,1.);
   Matrix C;
   Sum(A,B,&C);		 
   C.Print();
   }
   {
   Message("\n*************< Sum >*********************");
   MatrixRight A(3,3,
		 1.,1.,1.,
		    1.,1.,
		       1.);
   MatrixLeft B(3,3,
		1.,
		1.,1.,
		1.,1.,1.);
   Matrix C;
   Sum(B,A,&C);		 
   C.Print();
   }
   {
   Message("\n*************< Sum >*********************");
   MatrixRight A(3,3,
		 1.,2.,3.,
		    4.,5.,
		       6.);
   MatrixRight B(3,3,
		 9.,8.,7.,
		    6.,5.,
		       4.);

   Sum(&A,B);		 
   A.Print();
   }
   {
   Message("\n*************< Sum >*********************");
   MatrixRight A(3,3,
		 1.,2.,3.,
		    4.,5.,
		       6.);
   MatrixRight B(3,3,
		 9.,8.,7.,
		    6.,5.,
		       4.);

   Sum(A,&B);		 
   B.Print();
   }
   {
   Message("\n*************< Sum >*********************");
   MatrixRight A(3,3,
		 1.,2.,3.,
		    4.,5.,
		       6.);
   Sum(&A);		 
   A.Print();
   }
   {
   Message("\n*************< Sum >*********************");
   MatrixRight A(3,3,
		 1.,2.,3.,
		    4.,5.,
		       6.);
   Sum(A,A,&A);
   A.Print();
   }
   {
   Message("\n******** Operators + - MatrixRight  ******");
   float w[]={
	     1.,2.,3.,
		4.,5.,
		   6.
	     };
   MatrixRight A(3,3,w),B(3,3,1.,2.,3.,4.,5.,6.);
   MatrixRight C=A+B;		C.Print("A+B");
   A += B;                      A.Print("A += B");
   C=A-B;			C.Print("A-B");
   A -= B;                      A.Print("A -= B");
   C = -A;                      C.Print("C = -A");
   }
   {
   Message("\n******** Operators *  MatrixRight  ******");
   MatrixRight R1(3,3,
	     1.,2.,3.,
		4.,5.,
		   6.);
   Matrix M1(3,3,
	     1.,2.,3.,
	     0.,4.,5.,
	     0.,0.,6.);
   MatrixRight R2(3,3,
	     7.,2.,9.,
		1.,5.,
		   3.);
   Matrix M2(3,3,
	     7.,2.,9.,
	     0.,1.,5.,
	     0.,0.,3.);
   MatrixRight R = R1*R2;	R.Print("R1*R2");
   Matrix M = M1*M2;            M.Print("M1*M2");
   MatrixRight C=3.*R;         	C.Print("3.*A");
   C=R/3.;                     	C.Print("A/3.");
   R*=3.;                       R.Print("A*=3.");
   R/=3.;                       R.Print("A/=3.");
   if(R == R)fprintf(bzzFileOut,"\nR == R");
   if(R != R1)fprintf(bzzFileOut,"\nR != R1");
   }
}

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

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

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

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

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

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

//   ================================================================
//   ===============   Save and Recover   ===========================
//   ================================================================
void Ex_MatrixRight_4(void)
{
   {
   Message("\n** WARNING: Recover and Save  uses drive C:\\TEMP **");
   TempFile f1("C:\\TEMP\\");
   MatrixRight 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 
   MatrixRight B('*',f1.FileName()); // initialises B
   B.Print("B");
   A.Print("A");
   MatrixRight C;
   Recover(&C,'*',f1.FileName()); // C is resized
   C.Print("C");
   }
}

//   ================================================================
//   ===================   INVERSE  =================================
//   ================================================================
void Ex_MatrixRight_5(void)
{
    {
    Message("\n********** Inverse **************");
     MatrixRight R(3,3,
		 1.,6.,9.,
		    7.,5.,
		       6.);
     MatrixRight B=R;
     Inverse(&B);
     B.Print("Inverse");
     MatrixRight C=R*B;
     C.Print("control");
     float v[]={16.,12.,6.};
     Vector x(3,v);
     Vector y = B*x;
     y.Print("R-1 * x");
    }
    }

//   ================================================================
//   ===================   RRT RTR   ================================
//   ================================================================
void Ex_MatrixRight_6(void)
{
   {
   Message("\n********** RRT and RTR **************");
   Matrix A(3,3,
	     11.,2.1,3.2,
	     0.,4.1,55.,
	     0.,0.,6.5);
   MatrixRight R(3,3,
		 11.,2.1,3.2,
		     4.1,55.,
			 6.5);
   Matrix C;
   TProduct(A,A,&C);
   C.Print("ATA");
   MatrixSymm S;
   TProduct(R,&S);
   S.Print("RTR");
   ProductT(&A);
   A.Print("AAT");
   ProductT(R,&S);
   S.Print("RRT");
   }
}

//   ================================================================
//   ============   RIGHT SYSTEM SOLUTIONS  =========================
//   ================================================================
void Ex_MatrixRight_7(void)
{
   {
   Message("\nSystem solutions bx");
   float w[]={
	     1.,2.,3.,
		4.,5.,
		   6.
	     };
   MatrixRight A(3,3,w);    A.Print("MatrixRight A(3,w)");
   Vector b(3); b(1) = 6.; b(2) = 9.; b(3) = 6.;
   fprintf(bzzFileOut,"\ndet %e",A.Determinant());
   Solve(A,&b);
   b.Print("Solve(A,&b)");
   }
   {
   MatrixRight R(4,4,
	       1.,2.,3.,4.,
		  5.,6.,7.,
		     8.,9.,
			10.);
   Matrix b(4,2,10.,20.,
		18.,36.,
		17.,34.,
		10.,20.),x;
   Solve(R,b,&x);
   b.Print("Known terms");
   x.Print("Solution");
   }

   {
   Message("\nSystem solutions BX");   
   MatrixRight A(3,3,
		 1.,2.,3.,
		    4.,5.,
		       6.);
   Matrix B(3,2,
	    6.,7.,
	    9.,9.,
	    6.,6.);
   Solve(A,&B);

   B.Print();
   }
   {
   Message("\n ** TranposeSolve **");
   MatrixRight R(4,4,
	       1.,2.,3.,4.,
		  5.,6.,7.,
		     8.,9.,
			10.);
   Vector b(4,1.,7.,17.,30.),x;
   TransposeSolve(R,&b);
   b.Print("Known terms");
   x.Print("Solution");
   }
   {
   Message("\n ** TranposeSolve **");
   MatrixRight R(4,4,
	       1.,2.,3.,4.,
		  5.,6.,7.,
		     8.,9.,
			10.);
   Vector b(4,1.,7.,17.,30.),x;
   TransposeSolve(R,b,&x);
   b.Print("Known terms");
   x.Print("Solution");
   }

}