//   ======================< EXLQ.CPP >============================
//   * Examples of class FactoredLQ				  *
//   * 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: exLQ <number> <file name > <Y/N>		  *
//   * Example: 						  *
//   * 	EXLQ 1 pro.doc						  *
//   * 	Executes example 1 with  messages and results on pro.doc  *
//   ==============================================================

// Prototype examples
void Ex_LQ_1(void);
void Ex_LQ_2(void);
void Ex_LQ_3(void);
void Ex_LQ_4(void);
void Ex_LQ_5(void);
void Ex_LQ_6(void);
void Ex_LQ_7(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 "factqrlq.hpp"

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

   ExFact[1]=Ex_LQ_1;
   ExFact[2]=Ex_LQ_2;
   ExFact[3]=Ex_LQ_3;
   ExFact[4]=Ex_LQ_4;
   ExFact[5]=Ex_LQ_5;
   ExFact[6]=Ex_LQ_6;
   ExFact[7]=Ex_LQ_7;

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

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

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

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

void Ex_LQ_1(void)
{
Message("\n\n\n*** FactoredLQ  examples ************");
   {
   Message("\nFirst system\n");
   float s[]={60.,30.,20.,
	      30.,20.,15.,
	      20.,15.,12.};
   FactoredLQ A(3,3,s);
   float w[]={110.,65.,47.};
   Vector bx(3,w);
   Solve(A,&bx);
   bx.Print("Solution in bx");
   if(A.Singular() == 1)
      Message("\nSingular matrix\n");
   else
      Message("\nNon-singular matrix\n");
   Message("\nCondition number: %e",
      A.ConditionNumber());
   Message("\nDeterminant: %e",
      A.Determinant());
   }
}

void Ex_LQ_2(void)
{
   {
   Message("\nRow of zeros\n");
   FactoredLQ A(3,3,
	      60.,30.,20.,
	      0.,0.,0.,
	      20.,15.,12.);
   Vector bx(3,110.,0.,47.);
   Solve(A,&bx);
   bx.Print("Solution");
   A.Print("A");
   if(A.Singular() == 1)
      Message("\nSingular matrix\n");
   else
      Message("\nNon-singular matrix\n");
   Message("\nCondition number: %e",
      A.ConditionNumber());
   Message("\nDeterminant: %e",
      A.Determinant());
   }
}

void Ex_LQ_3(void)
{
   {
   Message("\nColumn of zeros\n");
   FactoredLQ A(3,3,
	      60.,0.,20.,
	      30.,0.,15.,
	      20.,0.,12.);
   Vector bx(3,80.,45.,37.);
   Solve(A,&bx);
   bx.Print("Solution");
   A.Print("A");
   if(A.Singular() == 1)
      Message("\nSingular matrix\n");
   else
      Message("\nNon-singular matrix\n");
   Message("\nCondition number: %e",
      A.ConditionNumber());
   Message("\nDeterminant: %e",
      A.Determinant());
   }
}

void Ex_LQ_4(void)
{
   {
   Message("\n**Checking the minimum norms of x\n");
   FactoredLQ A(1,2,1.,1.);
   Vector b(1,1.),x;
   Solve(&A,b,&x);
   x.Print("Solution x");
   }
}

void Ex_LQ_5(void)
{
   {
   Message("\n**Example with rectangular matrix\n");
   FactoredLQ A(3,4,
	      1.,4.,7.,10.,
	      2.,5.,8.,11.,
	      3.,6.,9.,12.);
   Vector b(3,1.,0.,0.),x;
   Solve(&A,b,&x);
   x.Print("Solution x");
   }
}

void Ex_LQ_6(void)
{
Message("\n*** Checking HyperSolve ****\n");
   {
   const float a = 50.;
   float s[]={11.+a,10.+a,14.+ a,
	      12.+a,11+a,-13.+a,
	      14.+a,13.+a,-66+a};
   FactoredLQ 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");
   }
}

void Ex_LQ_7(void)
{
   Message("\n**Checking the Condition number\n");
   {
   const float a=.01;
   FactoredLQ A(2,2,
		1.,1.+a,
		1.-a,1.);
   float w[]={2.+a,2.-a};
   Vector bx(2,w);
   Message("\nThe true condition number is: %f",
	   (2.+a)*(2.+a)/(a*a));
   Message("\nThat estimated is: %f",A.ConditionNumber());
   }
   {
   Message("\n**Checking the condition number\n");
   const float a = 50.;
   float s[]={11.+a,10.+a,14.+ a,
	      12.+a,11+a,-13.+a,
	      14.+a,13.+a,-66+a};
   FactoredLQ A(3,3,s);
   float w[]={3.*a+35,3.*a+10,3.*a-39.};
   Vector bx(3,w);
   Message("\nThe true condition number is: %f",
	   (3.*a+35.)*(1843.+166*a));
   Message("\nThat estimated is: %f",A.ConditionNumber());
   }
   {
   Message("\n**Checking the condition number\n");
   const float a = 50.;
   float s[]={11.+a,10.+a,14.+ a,
	      12.+a,11+a,-13.+a,
	      14.+a,13.+a,-66+a};
   FactoredLQ A(3,3,s);
   Message("\nThe true condition number is: %f",
	   (3.*a+35.)*(1843.+166*a));
   Message("\nThat estimated is: %f",A.ConditionNumber());
   }
}