//   =====================< EXMATR2.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 a 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						  *
//   **************************************************************
//   * Examples:                                                  *
//   * exmatr1 0 pro.doc				          *
//   * Executes all examples Prints the result on pro.doc         *
//   * exmatr2 1 pro.doc					  *
//   * Executes example 1 with messages and results on PRO.DOC 	  *
//   * exmatr2							  *
//   * 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);

#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 = 7;
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;

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   =========================
//   ===========================================================

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

//   ===========================================================
//   =====================   Sum   =============================
//   ===========================================================
void Ex_Matrix_1(void)
{
   {
   Message("\n***********< Sum >*************\n");
   float w[]=
      {1.,2.,3.,4.,5.,6.,
       7.,8.,9.,10.,11.,12.};
   Matrix A(2,6,w);
   Matrix B(2,6,w);
   Matrix C;
   C = A + B;
   C.Print("C = A + B");
   }
   {
   Message("\n***********< Sum >*************\n");
   float w[]=
      {1.,2.,3.,4.,5.,6.,
       7.,8.,9.,10.,11.,12.};
   Matrix A(2,6,w);
   Matrix B(2,6,w);
   Matrix C;
   Sum(A,B,&C);
   C.Print("Sum(A,B,&C)");
   }
   {
   Message("\n***********< Sum >*************\n");
   float w[]=
      {1.,2.,3.,4.,5.,6.,
       7.,8.,9.,10.,11.,12.};
   Matrix A(2,6,w);
   Matrix B(2,6,w);
   A +=B;
   A.Print("A += B");
   }
   {
   Message("\n***********< Sum >*************\n");
   float w[]=
      {1.,2.,3.,4.,5.,6.,
       7.,8.,9.,10.,11.,12.};
   Matrix A(2,6,w);
   Matrix B(2,6,w);
   Sum(&A,B);
   A.Print("Sum(&A,B)");
   }
   {
   Message("\n***********< Sum >*************\n");
   float w[]=
      {1.,2.,3.,4.,5.,6.,
       7.,8.,9.,10.,11.,12.};
   Matrix A(2,6,w);
   Matrix B(2,6,w);
   Sum(A,&B);
   B.Print("Sum(A,&B)");
   }
   {
   Message("\n*************< Sum >*********************");
   Matrix A(2,6,
	    1.,2.,3.,4.,5.,6.,
	    7.,8.,9.,10.,11.,12.);
   Sum(&A);
   A.Print("Sum(&A)");
   }
   {
   Message("\n*************< Sum >*********************");
   Matrix A(2,6,
	    1.,2.,3.,4.,5.,6.,
	    7.,8.,9.,10.,11.,12.);
   Sum(A,A,&A);
   A.Print("Sum(A,A,&A)");
   }
}

//   ===========================================================
//   ==================   Difference   =========================
//   ===========================================================
void Ex_Matrix_2(void)
{
   {
   Message("\n***********< Difference >*************\n");
   float w[]=
      {1.,2.,3.,4.,5.,6.,
       7.,8.,9.,10.,11.,12.};
   Matrix A(2,6,w);
   Matrix B(2,6,w);
   Matrix C;
   C = A - B;
   C.Print("C = A - B");
   }
   {
   Message("\n***********< Difference >*************\n");
   float w[]=
      {1.,2.,3.,4.,5.,6.,
       7.,8.,9.,10.,11.,12.};
   Matrix A(2,6,w);
   Matrix B(2,6,w);
   Matrix C;
   Difference(A,B,&C);
   C.Print("Difference(A,B,&C)");
   }
   {
   Message("\n***********< Difference >*************\n");
   float w[]=
      {1.,2.,3.,4.,5.,6.,
       7.,8.,9.,10.,11.,12.};
   Matrix A(2,6,w);
   Matrix B(2,6,w);
   A -=B;
   A.Print("A -= B");
   }
   {
   Message("\n***********< Difference >*************\n");
   float w[]=
      {1.,2.,3.,4.,5.,6.,
       7.,8.,9.,10.,11.,12.};
   Matrix A(2,6,w);
   Matrix B(2,6,w);
   Difference(&A,B);
   A.Print("Difference(&A,B); A = A - B;");
   }
   {
   Message("\n***********< Difference >*************\n");
   float w[]=
      {1.,2.,3.,4.,5.,6.,
       7.,8.,9.,10.,11.,12.};
   Matrix A(2,6,w);
   Matrix B(2,6,w);
   Difference(A,&B);
   B.Print("Difference(A,&B); B = A - B;");
   }
}

//   ==========================================================
//   ====================   Minus   ===========================
//   ==========================================================
void Ex_Matrix_3(void)
{
   {
   Message("\n***********< Minus >*************\n");
   float a[]=
      {1.,3.,
       2.,-1.};
   float b[]=
      {2.,0.,-4.,
       3.,-2.,6.};
   Matrix A(2,2,a);
   Matrix B(2,3,b);
   Minus(A,&B);
   B.Print("Minus(A,&B); B = -A");
   }
   {
   Message("\n***********< Minus >*************\n");
   float a[]=
      {1.,3.,
       2.,-1.};
   float b[]=
      {2.,0.,-4.,
       3.,-2.,6.};
   Matrix A(2,2,a);
   Matrix B(2,3,b);
   B = - A;
   B.Print("B = -A");
   }
   {
   Message("\n***********< Minus >*************\n");
   float a[]=
      {1.,3.,
       2.,-1.};
   Matrix A(2,2,a);
   Minus(&A);
   A.Print("Minus(&A); A = -A");
   }
}

//   ==========================================================
//   ===================   Product   ==========================
//   ==========================================================
void Ex_Matrix_4(void)
{
   {
   Message("\n***********< Product >*************\n");
   float a[]=
      {9.,3.,6.,
       3.,-30.,15.};
   Matrix A(2,3,a);
   Matrix B;
   Product(3.,A,&B);
   B.Print("Product(3.,A,&B); B = 3.A");
   }
   {
   Message("\n***********< Product >*************\n");
   float a[]=
      {9.,3.,6.,
       3.,-30.,15.};
   Matrix A(2,3,a);
   Matrix B;
   B = 3.*A;
   B.Print("B =3.* A;");
   }
   {
   Message("\n***********< Product >*************\n");
   float a[]=
      {9.,3.,6.,
       3.,-30.,15.};
   Matrix A(2,3,a);
   Product(3.,&A);
   A.Print("Product(3.,&A)");
   }
   {
   Message("\n***********< Product >*************\n");
   float a[]=
      {9.,3.,6.,
       3.,-30.,15.};
   Matrix A(2,3,a);
   A *=3.;
   A.Print("A *=3.;");
   }
   {
   Message("\n***********< Product >*************\n");
    float a[]={1.,-1.,1.,-3.,2.,-1.,-2.,1.,0.};
	Matrix A(3,3,a);
    float b[]={1.,2.,3.,2.,4.,6.,1.,2.,3.};
	Matrix B(3,3,b);
    Matrix C;
    Product(A,B,&C);
	C.Print("Product(A,B,&C); C = AB = 0");
   }
   {
   Message("\n***********< Product >*************\n");
    float a[]={1.,-1.,1.,-3.,2.,-1.,-2.,1.,0.};
	Matrix A(3,3,a);
    float b[]={1.,2.,3.,2.,4.,6.,1.,2.,3.};
	Matrix B(3,3,b);
    Matrix C = A*B;
	C.Print("C = A*B must be = 0");
   }
   {
   Message("\n***********< Product >*************\n");
    float a[]={1.,-1.,1.,-3.,2.,-1.,-2.,1.,0.};
	Matrix A(3,3,a);
    float b[]={1.,2.,3.,2.,4.,6.,1.,2.,3.};
	Matrix B(3,3,b);
    Product(&A,B);
	A.Print("Product(&A,B); A = AB = 0");
   }
   {
   Message("\n***********< Product >*************\n");
    float a[]={1.,-1.,1.,-3.,2.,-1.,-2.,1.,0.};
	Matrix A(3,3,a);
    float b[]={1.,2.,3.,2.,4.,6.,1.,2.,3.};
	Matrix B(3,3,b);
    A *=B;
	A.Print("A *= B; A = AB = 0");
   }
   {
   Message("\n***********< Product >*************\n");
    float a[]={1.,-1.,1.,-3.,2.,-1.,-2.,1.,0.};
	Matrix A(3,3,a);
    float b[]={1.,2.,3.,2.,4.,6.,1.,2.,3.};
	Matrix B(3,3,b);
    Product(A,B,&B);
	B.Print("Product(A,B,&B); B = AB = 0");
   }
   {
   Message("\n***********< Product >*************\n");
    float a[]={1.,-1.,1.,-3.,2.,-1.,-2.,1.,0.};
	Matrix A(3,3,a);
    Matrix C;
    Product(A,A,&C);
	C.Print("Product(A,A,&C); C = AA ");
    Product(A,A,&A);
	A.Print("Product(A,A,&A); A = AA ");
   }
}

//   ===========================================================
//   ==================   TProduct   ===========================
//   ===========================================================
void Ex_Matrix_5(void)
{
   {
   Message("\n***********< TProduct >*************\n");
	Matrix A(3,3,
		 1.,-3.,-2.,
		 -1.,2., 1.,
		 1.,-1.,0.);
    float b[]={1.,2.,3.,2.,4.,6.,1.,2.,3.};
	Matrix B(3,3,b);
    Matrix C;
    TProduct(A,B,&C);
	C.Print("TProduct(A,B,&C); C = ATB = 0");
   }
   {
   Message("\n***********< TProduct >*************\n");
	Matrix A(3,3,
		 1.,-3.,-2.,
		 -1.,2., 1.,
		 1.,-1.,0.);
    float b[]={1.,2.,3.,2.,4.,6.,1.,2.,3.};
	Matrix B(3,3,b);
    Matrix C =A%B;
	C.Print("C = A%B; C = ATB = 0");
   }
   {
   Message("\n***********< TProduct >*************\n");
	Matrix A(3,3,
		 1.,-3.,-2.,
		 -1.,2., 1.,
		 1.,-1.,0.);
    float b[]={1.,2.,3.,2.,4.,6.,1.,2.,3.};
	Matrix B(3,3,b);
     A %= B;
	A.Print("A %=B; A = ATB = 0");
   }
   {
   Message("\n***********< TProduct >*************\n");
	Matrix A(3,3,
		 1.,-3.,-2.,
		 -1.,2., 1.,
		 1.,-1.,0.);
    float b[]={1.,2.,3.,2.,4.,6.,1.,2.,3.};
	Matrix B(3,3,b);
     TProduct(&A,B);
	A.Print("TProduct(&A,B); A = ATB = 0");
   }
   {
   Message("\n***********< TProduct >*************\n");
	Matrix A(3,3,
		 1.,-3.,-2.,
		 -1.,2., 1.,
		 1.,-1.,0.);
    float b[]={1.,2.,3.,2.,4.,6.,1.,2.,3.};
	Matrix B(3,3,b);
     TProduct(A,&B);
	B.Print("TProduct(A,&B); B = ATB = 0");
   }
   {
   Message("\n***********< TProduct >*************\n");
	Matrix A(5,2,
		 1.,2.,
		 3.,4.,
		 5.,6.,
		 7.,8.,
		 9.,10.);
     Matrix C;
     TProduct(A,A,&C);
     C.Print("TProduct(A,A,&C);");
   }
   {
   Message("\n***********< TProduct >*************\n");
	Matrix A(5,2,
		 1.,2.,
		 3.,4.,
		 5.,6.,
		 7.,8.,
		 9.,10.);
     TProduct(&A);
	A.Print("TProduct(&A);");
   }
}

//   ===========================================================
//   =================   ProductT   ============================
//   ===========================================================
void Ex_Matrix_6(void)
{
   {
   Message("\n***********< ProductT >*************\n");
   float a[]=
      {9.,3.,6.,
       3.,-30.,15.,
       1.,2.,3.};
   Matrix A(3,3,a);
   ProductT(&A);
   A.Print("ProductT(&A)");
   }
   {
   Message("\n***********< ProductT >*************\n");
	Matrix A(2.5,
		 1.,2.,3.,4.,5.,
		 6.,7.,8.,9.,10.);
     ProductT(&A);
	A.Print("ProductT(&A);");
   }
}

//   ============================================================
//   ==================   Division   ============================
//   ============================================================
void Ex_Matrix_7(void)
{
   {
   Message("\n***********< Division >*************\n");
   float a[]=
      {9.,3.,6.,
       3.,-30.,15.};
   Matrix A(2,3,a);
   Matrix B;
   Division(A,3.,&B);
   B.Print("Division(A,3.,&B)");
   }
   {
   Message("\n***********< Division >*************\n");
   float a[]=
      {9.,3.,6.,
       3.,-30.,15.};
   Matrix A(2,3,a);
   Matrix B;
   B = A/3.;
   B.Print("B = A/3.;");
   }
   {
   Message("\n***********< Division >*************\n");
   float a[]=
      {9.,3.,6.,
       3.,-30.,15.};
   Matrix A(2,3,a);
   Division(&A,3.);
   A.Print("Division(&A,3.)");
   }
   {
   Message("\n***********< Division >*************\n");
   float a[]=
      {9.,3.,6.,
       3.,-30.,15.};
   Matrix A(2,3,a);
   A /=3.;
   A.Print("A /=3.;");
   }
}