//   ======================< EXQR.CPP >==============================
//   * Examples of class FactoredQR				    *
//   * 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, left.obj, right.obj, symm.obj, 		    *
//   * factored.obj, factQRLQ.obj				    *
//   ****************************************************************
//   * To execute: exQR <number> <file name > <Y/N>		    *
//   ****************************************************************
//   * exQR 0 pro.doc				                    *
//   * Executes all examples Prints the result on pro.doc	    *
//   * exQR 1 pro.doc						    *
//   * executes example 1 with messages and results on PRO.DOC	    *
//   * exQR							    *
//   * executes everything and sends it to stdout		    *
//   ================================================================

// Prototype examples
void Ex_QR_1(void);
void Ex_QR_2(void);
void Ex_QR_3(void);
void Ex_QR_4(void);
void Ex_QR_5(void);
void Ex_QR_6(void);
void Ex_QR_7(void);
void Ex_QR_8(void);
void Ex_QR_9(void);
void Ex_QR_10(void);
void Ex_QR_11(void);
void Ex_QR_12(void);
void Ex_QR_13(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 "right.hpp"
#include "factored.hpp"
#include "factqrlq.hpp"

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

ExFact[1]=Ex_QR_1;
ExFact[2]=Ex_QR_2;
ExFact[3]=Ex_QR_3;
ExFact[4]=Ex_QR_4;
ExFact[5]=Ex_QR_5;
ExFact[6]=Ex_QR_6;
ExFact[7]=Ex_QR_7;
ExFact[8]=Ex_QR_8;
ExFact[9]=Ex_QR_9;
ExFact[10]=Ex_QR_10;
ExFact[11]=Ex_QR_11;
ExFact[12]=Ex_QR_12;
ExFact[13]=Ex_QR_13;

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

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

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

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

void Ex_QR_1(void)
{
Message("\n\n\n FactoredQR examples **************\n\n");
   {
   FactoredQR A(3,3,
		60.,30.,20.,
		30.,20.,15.,
		20.,15.,12.);
   Vector bx(3,110.,65.,47.);
   Solve(A,&bx);
   A.Print("A");
   bx.Print("Solution found:");
   if(A.Singular())
      Message("\nSingular matrix");
   else Message("\nNon-singular matrix");
   Message("\n\nCondition is %f",
      A.ConditionNumber());
   }
   {
   Message("\n\n\n Commutation from Matrix **************");
   Matrix B(3,3,
		 60.,30.,20.,
		 30.,20.,15.,
		 20.,15.,12.);
   FactoredQR A;
   CommuteMatrixToFactored(&B,&A);
   B.Print("B after Commute");
   Vector bx(3,110.,65.,47.);
   Solve(&A,&bx);
   bx.Print("x e b");
   }
}

void Ex_QR_2(void)
{
   {
   Message("\nMatrix with a row of zeros\n");
   float s[]={60.,30.,20.,
	      0.,0.,0.,
	      20.,15.,12.};
   FactoredQR A(3,3,s);
   float w[]={110.,0.,47.};
   Vector bx(3,w);
   Solve(A,&bx);
   A.Print("A");
   bx.Print("Solution found:");
   if(A.Singular())
      Message("\nSingular matrix");
   else Message("\nNon-singular matrix");
   Message("\n\nCondition is %f",
      A.ConditionNumber());
   MatrixRight R;
   A.GetMatrixR(&R);
   R.Print("Matrix R");
   Matrix Q;
   A.GetMatrixQ(&Q);
   Q.Print("Matrix Q");
   (Q*R).Print("Check");
   }
   {
   Message("\nMatrix with a column of zeros\n");
   FactoredQR A(3,3,
		60.,0.,20.,
		30.,0.,15.,
		20.,0.,12.);
   float w[]={80.,45.,32.};
   Vector bx(3,w);
   Solve(A,&bx);
   A.Print("A");
   bx.Print("Solution found:");
   if(A.Singular())
      Message("\nSingular matrix");
   else Message("\nNon-singular matrix");
   Message("\n\nCondition is %f",
      A.ConditionNumber());
   MatrixRight R;
   A.GetMatrixR(&R);
   R.Print("Matrix R");
   Matrix Q;
   A.GetMatrixQ(&Q);
   Q.Print("Matrix Q");
   (Q*R).Print("Check");
   }
}

void Ex_QR_3(void)
{
   {
   Message("\n\nFour rows and three columns\n");
   float s[]={60.,30.,20.,
	      30.,20.,15.,
	      30.,20.,15.,
	      20.,15.,12.};
   FactoredQR A(4,3,s);
   float w[]={110.,65.,65.,47.};
   Vector bx(4,w);
   Solve(A,&bx);
   bx.Print("x e b");
   A.Print("A");
   }
}

void Ex_QR_4(void)
{
   {
   Message("\n** Regression problem **\n");
   float s[]={0.,0.,1.,
	      4.,0.,1.,
	      0.,4.,1.,
	      5.,6.,1.};
   FactoredQR A(4,3,s);
   float w[]={0.,-16.,-16.,-61.};
   Vector bx(4,w);
   Solve(A,&bx);
   A.Print("A");
   bx.Print("Solution found:");
   if(A.Singular())
      Message("\nSingular matrix");
   else Message("\nNon-singular matrix");
   MatrixRight R;
   A.GetMatrixR(&R);
   R.Print("Matrix R");
   Matrix Q;
   A.GetMatrixQ(&Q);
   Q.Print("Matrix Q");
   (Q*R).Print("Check");
   }
}

void Ex_QR_5(void)
{
   {
   Message("\n**Badly conditioned System**\n");
   float s[]={1.,2.,3.,
	      4.,5.,6.,
	      7.,8.,9.,
	      10.,11.,12.};
   FactoredQR A(4,3,s);
   float w[]={1.,0.,0.,0.};
   Vector bx(4,w);
   Solve(A,&bx);
   A.Print("A");
   bx.Print("Solution found:");
   if(A.Singular())
      Message("\nSingular matrix");
   else Message("\nNon-singular matrix");
   }
}

void Ex_QR_6(void)
{
   {
   Message("\n**Final column eliminated**\n");
   float s[]={1.,2.,
	      4.,5.,
	      7.,8.,
	      10.,11.};
   FactoredQR A(4,2,s);
   float w[]={1.,0.,0.,0.};
   Vector bx(4,w);
   Solve(A,&bx);
   A.Print("A");
   bx.Print("Solution found:");
   if(A.Singular())
      Message("\nSingular matrix");
   else Message("\nNon-singular matrix");
   }
}

void Ex_QR_7(void)
{
   {
   Message("\nSystem singular and incompatible\n");
   float s[]={1.,2.,-3.,
	      3.,-1.,2.,
	      5.,3.,-4.};
   FactoredQR 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("\nEstimated condition number%f",
      A.ConditionNumber());
   }
}

void Ex_QR_8(void)
{
Message("\n\n** Checking ConditionNumber **\n\n");
   {
   Message("\n** First example **\n");
   const float a = 50.;
   FactoredQR 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("\n\nEstimated condition number%f",
      A.ConditionNumber());
   }
   {
   Message("\n** Second example **\n");
   const float a=.01;
   FactoredQR 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("\n\nEstimated condition number%f",
      A.ConditionNumber());
   }
}

void Ex_QR_9(void)
{
   {
   Message("\n** Proving GetMatrixP,Q,R **");
   float s[]={0.,0.,1.,
	      4.,0.,1.,
	      0.,4.,1.,
	      5.,6.,1.};
   FactoredQR A(4,3,s);
   Matrix G(4,3,s);
   float w[]={0.,-16.,-16.,-61.};
   Vector bx(4,w);
   Solve(A,&bx);
   bx.Print("x e b");
   A.Print("A");
   Matrix P;
   A.GetMatrixP(&P);
   P.Print("P");
   Matrix RR = P*G;
   RR.Print("RR");
   MatrixRight R;
   A.GetMatrixR(&R);
   R.Print("R");
   Matrix Q;
   A.GetMatrixQ(&Q);
   Q.Print("Q");
   Matrix F = Q*R;
   F.Print("A");
   }
}

void Ex_QR_10(void)
{
  {
  Message("\n**** Calculating inverse\n");
   float s[]={5.,2.,1.,0.,
	      0.,4.,1.,0.,
	      3.,0.,10.,2.,
	      1.,2.,0.,1.};
   FactoredQR A(4,4,s);
   Matrix B(4,4,s);
   Matrix inv;
   A.Inverse(&inv);
   Matrix C;
   Product(inv,B,&C);
   inv.Print("Inv");
   C.Print("Check");
   }

  {
  Message("\n**** Calculating inverse\n");
   float s[]={5.,2.,1.,0.,
	      0.,4.,1.,0.,
	      3.,0.,10.,2.,
	      1.,2.,0.,1.};
   FactoredQR A(4,4,s);
   Matrix B(4,4,s);
   Matrix inv = A.Inverse();
   Matrix C;
   Product(inv,B,&C);
   inv.Print("Inv");
   C.Print("Test");
   }
}

void Ex_QR_11(void)
{
  {
   Message("\nProving HyperSolve\n");
   float s[]={6.,3.,2.,
	      0.,2.,1.,
	      6.01,3.,2.};
   FactoredQR A(3,3,s);
   float w[]={11.,3.,11.01};
   Vector bx(3,w);
   Solve(A,&bx);
   bx.Print("with Solve must be 1");
   Vector cx(3,w);
   HyperSolve(A,&cx);
   cx.Print("with HyperSolve must be 1");
   A.Print("A");
   Message("\n\nCondition is %f",
      A.ConditionNumber());
   }
}

void Ex_QR_12(void)
{
   {
   Message("\nCalculating linear model parameters \n");
   FactoredQR 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);
   Solve(F,y,&b);
   b.Print("Parameters");
   }
   {
   Message("\nCalculating variance y\n");
   FactoredQR 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 = .33,vy;
   MatrixRight R;
   F.GetMatrixR(&R);
   TransposeSolve(R,f,&x);
   vy = s2*(x%x);
   Message("Variance of y for x=3 is %e",vy);
   }
   {
   Message("\nCalculating linear model parameters \n");
   FactoredQR F(6,3,
	    1.,1.,1.,
	    1.,2.,4.,
	    1.,3.,9.,
	    1.,4.,16.,
	    1.,5.,25.,
	    1.,6.,36.);
   Vector y(6,5.51,8.98,13.4,19.1,24.4,33.2),b;
   y.Print();
   Solve(F,y,&b);
   b.Print("Parameters");
   Vector res;
   F.GetResiduals(&res);
   res.Print();
   float norm2 = res.Norm2();
   Message("\n S2 = %e",norm2*norm2/3.);
   }
}

void Ex_QR_13(void)
{
   {
   Message("\nDetermination of null space\n");
   Matrix F(4,2,
	    1.,5.,
	    2.,6.,
	    3.,7.,
	    4.,8.);
   FactoredQR A = F;
   Matrix P;
   A.GetMatrixP(&P);
   Transpose(&P);
   Matrix Q(4,2,1,3,P);
   P.Print("P");
   Q.Print("Q");
   MatrixRight R;
   A.GetMatrixR(&R);
   R.Print("R");
   (F%Q).Print("Check");
   }

}