//   =====================< FACTQRLQ.HPP >======================
//   * Class FactoredQRLQ derived from Factored                *
//   * Class FactoredQR derived from FactoredQRLQ              *
//   * Class FactoredLQ derived from FactoredQRLQ              *
//   * Description: Chapter 15                                 *
//   * Scientific C++ Building Numerical Libraries             *
//   *                the Object-Oriented Way                  *
//   * by G. Buzzi-Ferraris                                    *
//   * Addison-Wesley (1993)                                   *
//   ===========================================================
//   * FactoredQR permits linear system solutions              *
//   * and solutions to related problems                       *
//   * by numRows >= numColumns                                *
//   * The  solution minimises the residual sum of squares.    *
//   * It is used for linear regression                        *
//   ***********************************************************
//   * Other functions specified for FactoredQR:               *
//   * void GetMatrixQ(Matrix *Q);                             *
//   *    calculates the matrix Q of the factorisation A = QR  *
//   * void GetMatrixP(Matrix *P);                             *
//   *    calculates the matrix P of the factorisation PA = R  *
//   * void GetMatrixR(MatrixRight *R);                        *
//   *    calculates the matrix R of the factorisation A = QR  *
//   * void GetResiduals(Vector *residuals);                   *
//   *    provides the residuals (only if numRows > numColumns)*
//   ===========================================================
//   * FactoredLQ permits linear system solution               *
//   * and solutions to related problems                       *
//   * by numRows <= numColumns                                *
//   * The solution minimises the solution norm                *
//   ***********************************************************
//   * Other specific functions for FactoredLQ:                *
//   * void GetMatrixQ(Matrix *Q);                             *
//   *    calculates the matrix Q of the factorisation A = LQ  *
//   * void GetMatrixP(Matrix *P);                             *
//   *    calculates the matrix P of the factorisation AP = L  *
//   * void GetMatrixL(MatrixLeft *L);                         *
//   *    calculates the matrix L of the factorisation A = LQ  *
//   ***********************************************************


#ifndef FACTOREDQRLQ_HPP
#define FACTOREDQRLQ_HPP

//   ===========================================================
//   ==========   prototype functions for QR and LQ   ==========
//   ===========================================================
   char QRFactorization
       (int m,int n,float **a,float *d);
   void QRSolution
       (int m,int n,float **a,float *d,float *b,float *x);
   float QRCondition
       (int n,float **a,float *d);

   char LQFactorization
       (int m,int n,float **a,float *d);
   void LQSolution
       (int m,int n,float **a,float *d,float *b,float *x);
   float LQCondition
       (int n,float **a,float *d);

//   ===========================================================
//   ===============   class FactoredQRLQ   ====================
//   ===========================================================
class Vector;
class Matrix;
class FactoredQRLQ : public Factored
   {
   friend void Householder
       (int j, int k,float *aux,float *d);
   friend void HouseholderApplyLeft
       (int ri, int rf,int ci,int cf,float **a,float *aux);
   friend void HouseholderApplyRight
       (int ri, int rf,int ci,int cf,float **a,float *aux);
   friend void Givens
       (float *x1,float *x2,float *c,float *s,float *w);
   friend void GivensApply
       (float *y1,float *y2,float c,float s,float w);
   friend void Hessenberg
       (int m,int n,float **a);
   friend void HessenbergP
       (int m,int n,float **a,float **p);
   friend void Bidiagonal
       (int m,int n,float **a,float *d,float *ds);
   friend void BidiagonalR
       (int m,int n,float **a,float *d,
       float *ds,float **an);

protected:
   int *indx;
   float *dqr;  //per LQ e QR
   float *bqr;  //per LQ e QR

   // initialisation of constructors
   void FurtherInit(void);

   //  initialisation of special vectors
   virtual void SpecificInitialize(void);

   //  deinitialisation of special vectors
   virtual void SpecificDeinitialize(void);

   // constructor A('*',3,5);
   FactoredQRLQ(char ch,int rows,int columns)
   : Factored(ch,rows,columns){FurtherInit();}

//   ===========================================================
//   ***********  Functions for linear algebra    **************
//   ===========================================================
   virtual void    Factorization(void)=0; 
   virtual void    Solution(Vector *bx)=0;
   virtual void    Solution(const Vector &b,Vector *x)=0;
   virtual void    TransposeSolution(Vector *bx)=0;
   virtual void    TransposeSolution
                 (const Vector &b,Vector *x)=0;
   virtual void    Solution(Matrix *BX)=0;
   virtual void    Solution(const Matrix &B,Matrix *X)=0;
   virtual float   Condition(void)=0;
   virtual double  DeterminantEvaluation(void);

public:
//   ===========================================================
//   ******************   constructors   ***********************
//   ===========================================================
   // default constructor type Factored A;
   FactoredQRLQ(void)
   : Factored(){FurtherInit();}

   // copy constructor
   FactoredQRLQ(const Factored &rval)
   : Factored(rval){FurtherInit();}

   // constructor A(3,5);
   FactoredQRLQ(int rows,int columns)
   : Factored(rows,columns){FurtherInit();}

   // constructor A(3,5,w);
   FactoredQRLQ(int rows,int columns,float *initvalues)
   : Factored(rows,columns,initvalues){FurtherInit();}

   // constructor from Matrix
   FactoredQRLQ(const Matrix &rval)
   : Factored(rval){FurtherInit();}

   // make a submatrix with rows,columns
   FactoredQRLQ(int rows,int columns,const Matrix &rval)
   : Factored(rows,columns,rval){FurtherInit();}

   // as above, commencing from irow,jcol
   FactoredQRLQ(int rows,int columns,
              int irow,int jcol,const Matrix &rval)
   : Factored(rows,columns,irow,jcol,rval){FurtherInit();}

//   ===========================================================
//   ********************   destructor   ***********************
//   ===========================================================
   ~FactoredQRLQ(void);

   };

//   ===========================================================
//   ================   class FactoredQR   =====================
//   ===========================================================
class FactoredQR : public FactoredQRLQ
   {
   friend void Householder
       (int j, int k,float *aux,float *d);
   friend void HouseholderApplyLeft
       (int ri, int rf,int ci,int cf,float **a,float *aux);
   friend void HouseholderApplyRight
       (int ri, int rf,int ci,int cf,float **a,float *aux);
   friend void Givens
       (float *x1,float *x2,float *c,float *s,float *w);
   friend void GivensApply
       (float *y1,float *y2,float c,float s,float w);
   friend void Hessenberg
       (int m,int n,float **a);
   friend void HessenbergP
       (int m,int n,float **a,float **p);
   friend void Bidiagonal
       (int m,int n,float **a,float *d,float *ds);
   friend void BidiagonalR
       (int m,int n,float **a,float *d,
       float *ds,float **an);

private:
//   ===========================================================
//   ************   Functions for linear algebra ***************
//   ===========================================================
   virtual void    Factorization(void);
   virtual void    Solution(Vector *bx);
   virtual void    Solution(const Vector &b,Vector *x);
   virtual void    TransposeSolution(Vector *bx);
   virtual void    TransposeSolution
                 (const Vector &b,Vector *x);
   virtual void    Solution(Matrix *BX);
   virtual void    Solution(const Matrix &B,Matrix *X);
   virtual float   Condition(void);

public:
//   ===========================================================
//   ******************   constructors   ***********************
//   ===========================================================
   // default constructor type Factored A;
   FactoredQR(void)
   : FactoredQRLQ(){}

   // copy constructor
   FactoredQR(const Factored &rval)
   : FactoredQRLQ(rval){}

   // constructor A(3,5);
   FactoredQR(int rows,int columns)
   : FactoredQRLQ(rows,columns){}

   // sized and initialised
   //  A(2,3,1.,2.,3.,4.,5.6.);
   FactoredQR(int rows,int columns,double a11,...);

   // constructor A(3,5,w);
   FactoredQR(int rows,int columns,float *initvalues)
   : FactoredQRLQ(rows,columns,initvalues){}

   // constructor from Matrix
   FactoredQR(const Matrix &rval)
   : FactoredQRLQ(rval){}

   // make a submatrix with rows,columns
   FactoredQR(int rows,int columns,const Matrix &rval)
   : FactoredQRLQ(rows,columns,rval){}

   //as above, commencing from irow,jcol
   FactoredQR(int rows,int columns,
             int irow,int jcol,const Matrix &rval)
   : FactoredQRLQ(rows,columns,irow,jcol,rval){}

//   ===========================================================
//   ********************   destructor   ***********************
//   ===========================================================
   ~FactoredQR(void);  

//   ===========================================================
//   ***************   assignment operators  *******************
//   ===========================================================
   void operator = (const FactoredQR &rval);
   FactoredQR &operator = (const Matrix &rval);

//   ===========================================================
//   ======   Functions which do not modify the matrix  ========
//   ===========================================================
   void GetMatrixQ(Matrix *Q); // A = QR
   void GetMatrixP(Matrix *P);  // PA = R
   void GetMatrixR(MatrixRight *R);
   void GetResiduals(Vector *residuals); // nRows > nColumns

   };

//   ===========================================================
//   ================   class FactoredLQ   =====================
//   ===========================================================
class FactoredLQ : public FactoredQRLQ
   {
   friend void Householder
       (int j, int k,float *aux,float *d);
   friend void HouseholderApplyLeft
       (int ri, int rf,int ci,int cf,float **a,float *aux);
   friend void HouseholderApplyRight
       (int ri, int rf,int ci,int cf,float **a,float *aux);
   friend void Givens
       (float *x1,float *x2,float *c,float *s,float *w);
   friend void GivensApply
       (float *y1,float *y2,float c,float s,float w);
   friend void Hessenberg
       (int m,int n,float **a);
   friend void HessenbergP
       (int m,int n,float **a,float **p);
   friend void Bidiagonal
       (int m,int n,float **a,float *d,float *ds);
   friend void BidiagonalR
       (int m,int n,float **a,float *d,
       float *ds,float **an);

private:
//   ===========================================================
//   ************   Functions for linear algebra ***************
//   ===========================================================
   virtual void    Factorization(void); 
   virtual void    Solution(Vector *bx);
   virtual void    Solution(const Vector &b,Vector *x);
   virtual void    TransposeSolution(Vector *bx);
   virtual void    TransposeSolution
                 (const Vector &b,Vector *x);
   virtual void    Solution(Matrix *BX);
   virtual void    Solution(const Matrix &B,Matrix *X);
   virtual float   Condition(void);

public:
//   ===========================================================
//   ******************   constructors   ***********************
//   ===========================================================
   // default constructor type Factored A;
   FactoredLQ(void)
   : FactoredQRLQ(){}

   // copy constructor
   FactoredLQ(const Factored &rval)
   : FactoredQRLQ(rval){}

   // constructor A(3,5);
   FactoredLQ(int rows,int columns)
   : FactoredQRLQ(rows,columns){}

   // sized and initialised
   //  A(2,3,1.,2.,3.,4.,5.6.);
   FactoredLQ(int rows,int columns,double a11,...);

   // constructor A(3,5,w);
   FactoredLQ(int rows,int columns,float *initvalues)
   : FactoredQRLQ(rows,columns,initvalues){}

   // constructor from Matrix
   FactoredLQ(const Matrix &rval)
   : FactoredQRLQ(rval){}

   // makes a submatrix with rows,columns
   FactoredLQ(int rows,int columns,const Matrix &rval)
   : FactoredQRLQ(rows,columns,rval){}

   //as above, commencing from irow,jcol
   FactoredLQ(int rows,int columns,
             int irow,int jcol,const Matrix &rval)
   : FactoredQRLQ(rows,columns,irow,jcol,rval){}

//   ===========================================================
//   ********************   destructor   ***********************
//   ===========================================================
   ~FactoredLQ(void);

//   ===========================================================
//   ***************   assignment operators  *******************
//   ===========================================================
   void operator = (const FactoredLQ &rval);
   FactoredLQ &operator = (const Matrix &rval);

//   ===========================================================
//   ======   Functions which do not modify the matrix  ========
//   ===========================================================
   void GetMatrixQ(Matrix *Q); // A = LQ
   void GetMatrixP(Matrix *P);  // AP = L
   void GetMatrixL(MatrixLeft *L);


   };


#endif FACTOREDQRLQ_HPP