//   =====================< EXGAUSS.CPP >========================
//   * Examples of FactoredGauss class				*
//   * Description: Chapter 15                                 *
//   * Scientific C++ Building Numerical Libraries             *
//   *                the Object-Oriented Way                  *
//   * by G. Buzzi-Ferraris                                    *
//   * Addison-Wesley (1993)                                   *
//   ************************************************************
//   * To be assembled with: utility.obj, vector.obj,		*
//   * matrix.obj, factored.obj, factplr.obj			*
//   ************************************************************
//   * To execute: exgauss <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 nothing is given then it writes on screen		*
//   * If one puts Y/N equal to Y it prints messages as well	*
//   * The default is Y						*
//   ************************************************************
//   * Examples:                                                *
//   * exgauss 0 pro.doc				        *
//   * Executes all examples Prints the result on PRO.DOC	*		*
//   * exgauss 1 pro.doc					*
//   * Executes example 1 with messages and results on PRO.DOC	*
//   * exgauss							*
//   * Executes everything and sends it to stdout		*
//   ============================================================

// Prototype examples
void Ex_FactGauss_1(void);
void Ex_FactGauss_2(void);
void Ex_FactGauss_3(void);
void Ex_FactGauss_4(void);
void Ex_FactGauss_5(void);
void Ex_FactGauss_6(void);
void Ex_FactGauss_7(void);
void Ex_FactGauss_8(void);
void Ex_FactGauss_9(void);

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <ctype.h>
#include <math.h>
#include "utility.hpp"
#include "vector.hpp"
#include "matrix.hpp"
#include "factored.hpp"
#include "factplr.hpp"

int main(int argc,char **argv)
{
const int NEXGAUSS = 9;
const int NMAX = NEXGAUSS+1;
int i;
void (*ExFact[NMAX])(void);
char myfilename[81];

ExFact[1]=Ex_FactGauss_1;
ExFact[2]=Ex_FactGauss_2;
ExFact[3]=Ex_FactGauss_3;
ExFact[4]=Ex_FactGauss_4;
ExFact[5]=Ex_FactGauss_5;
ExFact[6]=Ex_FactGauss_6;
ExFact[7]=Ex_FactGauss_7;
ExFact[8]=Ex_FactGauss_8;
ExFact[9]=Ex_FactGauss_9;

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

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

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

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

void Ex_FactGauss_1(void)
{
//   ============================================================
//   ===============  INITIALISATIONS ===========================
//   ============================================================
   Message("\n\n\n* Examples of initialisations *");
Message("\n Gauss is used, but it goes with everything\n");
   {
   Message("\n Direct initialisation\n");
   FactoredGauss A(3,3);
      A(1,1)=60.; A(1,2)=30.; A(1,3)=20.;
      A(2,1)=30.; A(2,2)=20.; A(2,3)=15.;
      A(3,1)=20.; A(3,2)=15.; A(3,3)=12.;
   A.Print("A");
   }
   {
   Message("\n Initialisation with constructor\n");
   FactoredGauss A(3,3,
		 60.,30.,20.,
		 30.,20.,15.,
		 20.,15.,12.);
   A.Print("A");
   }
   {
   Message("\n Initialisation from array\n");
   float s[]={60.,30.,20.,
	      30.,20.,15.,
	      20.,15.,12.};
   FactoredGauss A(3,3,s);
   A.Print("A");
   }
   {
   Message("\n\n\n Constructor from Matrix **************");
   Matrix B(3,3,
		 60.,30.,20.,
		 30.,20.,15.,
		 20.,15.,12.);
   FactoredGauss A = B;
   A.Print("A");
   B.Print("B");
   }
   {
   Message("\n\n\n Switching from Matrix **************");
   Matrix B(3,3,
		 60.,30.,20.,
		 30.,20.,15.,
		 20.,15.,12.);
   FactoredGauss A;
   CommuteMatrixToFactored(&B,&A);
   A.Print("A replaces B");
   B.Print("B is eliminated");
   }
   {
   Message("\n\n\n Assignment from Matrix **************");
   Matrix B(3,3,
		 60.,30.,20.,
		 30.,20.,15.,
		 20.,15.,12.);
   FactoredGauss A;
   A = B;
   A.Print("A");
   B.Print("B");
   }
   {
   Message("\n\n\n Constructor from SubMatrix *********");
   Matrix B(4,5,
		 60.,30.,20.,3.,2.,
		 30.,20.,15.,1.,5.,
		 20.,15.,12.,4.,5.,
		 1.,2.,3.,4.,5.);
   FactoredGauss A(3,3,B);
   A.Print("A");
   B.Print("B");
   }
   {
   Message("\n\n\n Constructor from internal SubMatrix *********");
   Matrix B(4,5,
		1.,2.,3.,4.,5.,
		1.,60.,30.,20.,1.,
		1.,30.,20.,15.,1.,
		1.,20.,15.,12.,1.);
   FactoredGauss A(3,3,2,2,B);
   A.Print("A");
   B.Print("B");
   }
}

void Ex_FactGauss_2(void)
{
//   ============================================================
//   ==================   SOLUTIONS   ===========================
//   ============================================================
Message("\n\n\n ******** Examples of solutions ******\n\n");
   {
   Message("\n\n\n Constructor from internal submatrix **");
   Matrix B(4,5,
		1.,2.,3.,4.,5.,
		1.,60.,30.,20.,1.,
		1.,30.,20.,15.,1.,
		1.,20.,15.,12.,1.);
   FactoredGauss A(3,3,2,2,B);
   Vector bx(3,110.,65.,47.);
   Message("\n Superimposes both A and b");
   Solve(&A,&bx);
   bx.Print("x e b");
   }
   {
   Message("\nTwo systems simultaneously\n");
   FactoredGauss A(3,3,
	      60.,30.,20.,
	      30.,20.,15.,
	      20.,15.,12.);
   Matrix BX(3,2,
	      110.,220.,
	      65.,130.,
	      47.,94.);
   Message("\nMatrix A saved, BX superimposed\n");
   Solve(A,&BX);
   A.Print("A");
   BX.Print("x e b first column 1 second 2");
   }
   {
   Message("\nMatrix badly conditioned\n");
   const int n=4;
   FactoredGauss A(n,n);
   Vector bx(n);
   for(int i = 1;i <= n;i++)
      {
      bx(i) = 0.;
      for(int j = 1;j <= n;j++)
	 {
	 A(i,j) =  1./(double)(i+j-1);
	 bx(i) +=  1./(double)(i+j-1);
	 }
      }
   Solve(A,&bx);
   bx.Print("Solution x e b");
   A.Print("A");
   }
}

void Ex_FactGauss_3(void)
{
//   ===========================================================
//   ===============   CONDITION-NUMBER   ======================
//   ===========================================================
Message("\n\n** Check ConditionNumber **\n\n");
   {
   Message("\n** First example **\n");
   const float a = 50.;
   FactoredGauss A(3,3,
		   11.+a,10.+a,14.+ a,
		   12.+a,11+a,-13.+a,
		   14.+a,13.+a,-66+a);
   Vector bx(3,3.*a+35,3.*a+10,3.*a-39.);
   Solve(A,&bx);
   bx.Print("x e b");
   A.Print("A");
   Message("\n\nTrue condition number%f",
      (3.*a+35.)*(1843.+166*a));
   Message("\nEstimated condition number%f",
      A.ConditionNumber());
   }
   {
   Message("\n** Second example **\n");
   const float a=.01;
   FactoredGauss A(2,2,
		   1.,1.+a,
		   1.-a,1.);
   float w[]={2.+a,2.-a};
   Vector bx(2,w);
   Solve(A,&bx);
   bx.Print("x e b");
   A.Print("A");
   Message("\n\nTrue condition number%f",
      (3.*a+35.)*(1843.+166*a));
   Message("\nEstimated condition number%f",
      A.ConditionNumber());
   }
}

void Ex_FactGauss_4(void)
{
//   ===========================================================
//   ==================   TRANSPOSE   ==========================
//   ===========================================================
Message("\n\n\n Example of TransposeSolve\n\n");
   {
   FactoredGauss A(3,3,
		   6.,4.,2.,
		   3.,2.,5.,
		   2.,1.,2.);
   Vector bx(3,11.,7.,9.);
   TransposeSolve(&A,&bx);
   bx.Print("Solution con TransposeSolve");
   A.Print("A");
   }
}

void Ex_FactGauss_5(void)
{
//   ===========================================================
//   ===============  ACCESS FUNCTIONS   =======================
//   ===========================================================
Message("\n\n\n** Examples of Set ***\n\n");
   {
   float s[]={3.,2.,1.,
	      6.,0.,1.,
	      0.,0.,1.};
   FactoredGauss A(3,3,s);
   float w[]={6.,7.,1.};
   Vector bx(3,w);
   Solve(A,&bx);
   bx.Print("x and b must be 1");
   A.Print("A");

   Message("\nModifies coefficient A(2,1) of A");
   A(2,1)=8.; bx(1)=6.;bx(2)=9.;bx(3)=1.;
   Solve(A,&bx);
   bx.Print("x and b must be 1");
   A.Print("A");

   Message("\nModify row 2 in A");
   Vector y(3,1.,1.,1.);
   A.SetRow(2,y);
   bx(1)=6.;bx(2)=3.;bx(3)=1.;
   Solve(A,&bx);
   bx.Print("x and b must be 1");
   A.Print("A");

   Message("\nModify column 3 in A");
   y *=2.;
   A.SetColumn(3,y); bx(1)=7.;bx(2)=4.;bx(3)=2.;
   Solve(A,&bx);
   bx.Print("x and b must be 1");
   A.Print("A");
   }
}
void Ex_FactGauss_6(void)
{
  {
//   ===========================================================
//   ====================  INVERSE   ===========================
//   ===========================================================
  Message("\n**** Calculate Inverse\n");
   float s[]={5.,2.,1.,0.,
	      0.,4.,1.,0.,
	      3.,0.,10.,2.,
	      1.,2.,0.,1.};
   FactoredGauss A(4,4,s);
   Matrix B(4,4,s);
   Matrix inv;
   A.Inverse(&inv);
   Matrix C;
   Product(inv,B,&C);
   inv.Print("Inverse");
   C.Print("Checks if identical to I");
   }
  {
  Message("\n**** Calculate Inverse\n");
   float s[]={5.,2.,1.,0.,
	      0.,4.,1.,0.,
	      3.,0.,10.,2.,
	      1.,2.,0.,1.};
   FactoredGauss A(4,4,s);
   Matrix B(4,4,s);
   Matrix inv = A.Inverse();
   Matrix C;
   Product(inv,B,&C);
   inv.Print("Inverse");
   C.Print("Check");
   }
}

void Ex_FactGauss_7(void)
{
//   ===========================================================
//   ==================   SOLUTIONS   ==========================
//   ===========================================================
   {
   Message("\nIncompatible and singular system\n");
   float s[]={1.,2.,-3.,
	      3.,-1.,2.,
	      5.,3.,-4.};
   FactoredGauss A(3,3,s);
   Vector bx(3,-1.,7.,2.);
   Solve(&A,&bx);
   bx.Print("Solution found:");
   if(A.Singular())
      Message("\nSingular matrix ");
   else Message("\nNon-singular matrix ");
   Message("\n\nCondition is %f",
      A.ConditionNumber());
   }
}

void Ex_FactGauss_8(void)
{
//   ==========================================================
//   ======================  HYPER  ===========================
//   ==========================================================
Message("\n ****** HyperSolve ****");
   {
   const float a = 50.;
   float s[]={11.+a,10.+a,14.+ a,
	      12.+a,11+a,-13.+a,
	      14.+a,13.+a,-66+a};
   FactoredGauss A(3,3,s);
   float w[]={3.*a+35,3.*a+10,3.*a-39.};
   Vector bx(3,w);
   Solve(A,&bx);
   bx.Print("Solution with Solve");
   Vector cx(3,w);
   HyperSolve(A,&cx);
   cx.Print("Solution with HyperSolve");
   A.Print("A");
   Message("\n\nTrue condition number%f",
      (3.*a+35.)*(1843.+166*a));
   Message("\n\nEstimated condition number%f",
      A.ConditionNumber());
   }
   {
   Message("\nHyperSolve Test\n");
   float s[]={6.,3.,2.,
	      0.,2.,1.,
	      6.01,3.,2.};
   FactoredGauss A(3,3,s);
   float w[]={11.,3.,11.01};
   Vector bx(3,w);
   Solve(A,&bx);
   bx.Print("with Solve it must be 1");
   Vector cx(3,w);
   HyperSolve(A,&cx);
   cx.Print("with HyperSolve it must be 1");
   A.Print("A");
   Message("\n\nCondition is %f",
      A.ConditionNumber());
   }
}

void Ex_FactGauss_9(void)
{
//   ==========================================================
//   =================  LINEAR  MODELS ========================
//   ==========================================================
   {
   Message("\nCalculating linear model parameters\n");
   Message("\nThis path is not useful (better to use QR)\n");
   Matrix F(6,3,
	    1.,1.,1.,
	    1.,2.,4.,
	    1.,3.,9.,
	    1.,4.,16.,
	    1.,5.,25.,
	    1.,6.,36.);
   Vector y(6),b;
   for(int i=1;i <=6;i++)
      y(i) = 3. + i*(2. + .5*i);
   Matrix FTF = F%F;
   Vector FTy = F%y;
   FactoredGauss commFTF;
   CommuteMatrixToFactored(&FTF,&commFTF);
   Solve(&commFTF,FTy,&b);
   b.Print("Parameters");
   }
   {
   Message("\nCalculating variance y\n");
   Message("\nBetter to use QR\n");
   Matrix F(6,3,
	    1.,1.,1.,
	    1.,2.,4.,
	    1.,3.,9.,
	    1.,4.,16.,
	    1.,5.,25.,
	    1.,6.,36.);
   Vector f(3,1.,3.,9.),x;
   float s2 = 1.e-4,vy;
   Matrix FTF = F%F;
   FactoredGauss commFTF;
   CommuteMatrixToFactored(&FTF,&commFTF);
   Solve(&commFTF,f,&x);
   vy = s2*(x%x);
   Message("The variance of y for x=3 is %e",vy);
   }
}