///////////////////////////////////////////////////////////////////////
//linss.h
#ifndef LINSS_H
#define LINSS_H

#include "vector.h"
#include "matrix.h"

class Linss{
  int n;
  Vector x, y;
  Matrix lhsMat(double lam) const;

public:
  Linss(const Vector& X, const Vector& Y);
  Vector smooth(double lam, double* GCV) const;
}; 
#endif

//linss.cpp
#include "linss.h"

Linss::Linss(const Vector& X, const Vector& Y){
  x = Vector(X);
  y = Vector(Y);
  n = x.getnRows();
}

Vector Linss::smooth(double lam, double* GCV) const{
  Matrix G(n, n, 0.);
  Matrix H = lhsMat(lam);
  Vector fit = H.bandChol(y, G, 1);
 
  //compute the GCV criterion
  double temp = 1/(G[n - 1][n - 1]*G[n - 1][n - 1]);
  double tr = temp;
  *GCV = (y[n - 1] - fit[n - 1])*(y[n - 1] - fit[n - 1]);
  for(int i = (n - 2); i >= 0; i--){    
    temp = (1 + G[i + 1][i]*G[i + 1][i]*temp)/(G[i][i]*G[i][i]);
    *GCV += (y[i] - fit[i])*(y[i] - fit[i]);
    tr += temp;
  }
  *GCV /= (1. - tr/((double)n))*(1. - tr/((double)n))*(double)n;

  return fit;
}

Matrix Linss::lhsMat(double lam) const{
  double** pH = new double*[n];
  pH[0] = new double[n];
  pH[n - 1] = new double[n];
  pH[0][0] = 1. + (lam/(x[1] - x[0]));
  pH[0][1] = -lam/(x[1] - x[0]);
  pH[n - 1][n - 2] = -lam/(x[n - 1] - x[n - 2]);
  pH[n - 1][n - 1] = 1. + (lam/(x[n - 1] - x[n - 2]));
  for(int i = 1; i < (n - 1); i++){
    pH[i] = new double[n];
    pH[i][i - 1] = -lam/(x[i] - x[i - 1]);
    pH[i][i] = 1 + (lam/(x[i] - x[i - 1])) + (lam/(x[i + 1] - x[i]));
    pH[i][i + 1] = -lam/(x[i + 1] - x[i]);
  }
  Matrix H(n, n, pH);

  for(int i = 0; i < n; i++)
    delete[] pH[i];
  delete[] pH;

  return H;
}

// linssDriver.cpp 
#include <cstdlib>
#include <iostream>
#include <fstream>
#include "vector.h"
#include "linss.h"

int main(int argc, char* argv[]){
  double lam = atof(argv[1]);
  int n = atoi(argv[2]);
  std::ifstream inFile(argv[3]);
  double* py = new double[n];
  double* px = new double[n];

  for(int i = 0; i < n; i++){
    inFile >> px[i];
    inFile >> py[i];
  }
  inFile.close();
  Vector y(n, py);
  Vector x(n, px);
  Linss S(x, y);
  double gcv = 0.;
  Vector fit = S.smooth(lam, &gcv);
  std::cout << "The value of GCV at lambda = " 
            << lam << " is " << gcv << std::endl;

  return 0;  
}
///////////////////////////////////////////////////////////////////////
//linssDriver.cc
#include "R.h" 
#include "vector.hh"
#include "linss.hh"

extern "C" {

  void linss(int* n, double* lam, double* px, double* py, 
             double* pfit, double* GCV){

    Vector y(*n, py);
    Vector x(*n, px);

    Linss S(x, y);
    Vector fit = S.smooth(*lam, GCV);
    for(int i = 0; i < *n; i++)
      pfit[i] = fit[i];
  }
}
///////////////////////////////////////////////////////////////////////
//inheritEx.cpp
#include <iostream>

using std::cout; using std::endl;

struct A{
  int iA;
  A(int ia){iA = ia;}
  void output(){f();}

  void f(){cout << iA << endl;}
  //virtual void f(){cout << iA << endl;}
};

struct B : public A {
  int iB;
  B(int ia, int ib):A(ia){
    iB = ib;
  }
  void f(){cout << iA << " " << iB << endl;}
};

int main(){
  A aObj(2);
  B bObj(3, 5);
  A* pA;
  pA = &aObj;
  pA->f();
  pA = &bObj;
  pA->f();
  bObj.output();
  return 0;
}
///////////////////////////////////////////////////////////////////////
//simpleMat.h
#ifndef MATRIX_H
#define MATRIX_H
#include <iostream>
#include "vector.h"

class Matrix{
  int nRows, nCols;
  double** pA; 
  Vector gauss(const Vector& RHS) const;

 public:
  //constructors and destructors
  Matrix(int nrows = 0, int ncols = 0, double a = 0.);
  Matrix(int nrows, int ncols, const double* const* pa);
  Matrix(const Matrix& A);
  virtual ~Matrix();

  //overloaded operators
  Matrix& operator=(const Matrix& B);
  const double* operator[](int i) const {return pA[i];}

  //solution of a linear system
  virtual Vector solve(const Vector& RHS) const{
    std::cout << "Solution via gauss from class Matrix" << std::endl;
    return gauss(RHS);
  }

  //utilities and accessors
  void printMatrix() const;
  int getNRows() const {return nRows;}
  int getNCols() const {return nCols;}
};
#endif

//simpleMatrix.cpp
#include "simpleMatrix.h"
#include <cstdlib>
#include <cmath>
#include <limits>

using std::cout; using std::endl;

//diagonal matrix constructor
Matrix:: Matrix(int nrows, int ncols, double a){
  nRows = nrows;
  ncols = ncols;
  //set pA to null pointer in default case
  if(nRows == 0 || nCols == 0){
    pA = 0;
    return;
  }

  pA = new double*[nRows];

  for(int i = 0; i < nRows; i++){
    pA[i] = new double[nCols];
    for(int j = 0; j < nCols; j++){
      if(i == j)
        pA[i][j] = a;
      else
        pA[i][j] = 0.;
    }
  }
}

//constructor from pointer-to-pointer 
Matrix::Matrix(int nrows, int ncols, const double* const* pa){
  nRows = nrows;
  nCols = ncols;
  pA = new double*[nRows];
  
  for(int i = 0; i < nRows; i++){
    pA[i] = new double[nCols];
  }

  for(int i = 0; i < nRows; i++)
    for(int j = 0; j < nCols; j++){
      pA[i][j] = pa[i][j];
    }
}

//copy constructor
Matrix::Matrix(const Matrix& A){
  nRows = A.nRows;
  nCols = A.nCols;
  pA = new double*[nRows];
      
  for(int i = 0; i < nRows; i++){
    pA[i] = new double[nCols];
  }
  
  for(int i = 0; i < nRows; i++)
    for(int j = 0; j < nCols; j++)
      pA[i][j] = A.pA[i][j];
}    
      
//destructor
Matrix::~Matrix(){   
  if(pA != 0){
    for(int i = 0; i < nRows; i++)
      if(pA[i] != 0)
        delete[] pA[i];
    delete[] pA;
  }
}
    
//Overloaded = operators:

Matrix& Matrix::operator=(const Matrix& A){

  if(this == &A)//avoid self assignment
    return *this;

  if(pA != 0){
    for(int i = 0; i < nRows; i++)//release current memory
      delete[] pA[i];
    delete[] pA;
  }
    
  //define/redefine object's members
  nRows = A.nRows; 
  nCols = A.nCols;
  pA = new double*[nRows];
      
  for(int i = 0; i < nRows; i++){
    pA[i] = new double[nCols];
  }

  for(int i = 0; i < nRows; i++)
    for(int j = 0; j < nCols; j++)
      pA[i][j] = A.pA[i][j];

  return *this;
}    


Vector Matrix::gauss(const Vector& RHS) const{
  double multiplier = 0;
  Matrix G(nRows, nRows, pA);
  Vector h(nRows, RHS.pA);
  
  for(int j = 0; j < nCols; j++){//column to be swept
    
    for(int jj = j + 1; jj < nCols; jj++){//current operation column

      if(G[j][j] == 0){
        cout << "Oops! Division by 0!" << endl;
        exit(1);
      }

      for(int i = j + 1; i < nRows; i++){//work down rows of jjth column
        multiplier = G[i][j]/G[j][j];
        G.pA[i][jj] = G.pA[i][jj] - multiplier*G.pA[j][jj];
      }
    }
    
    for(int i = j + 1; i < nRows; i++){//do the same to the RHS
      multiplier = G[i][j]/G[j][j];
      h.pA[i] = h.pA[i] - multiplier*h.pA[j];
    }
  }
  //now backsolve

  double temp;//temporary storage
  
  Vector b(nRows, 0.);//solution vector
  
  //initialize recursion on last row of all right-hand-sides
  if(G.pA[nRows - 1][nRows - 1] != 0)
    b.pA[nRows - 1] = h.pA[nRows - 1]/G.pA[nRows - 1][nRows - 1];
  else{
    cout << "Singular system!" << endl;
    exit(1);
  }
  
  //now work through the remaining rows
  for(int i = (nRows - 2); i >= 0; i--){
        if(G.pA[i][i] != 0){
      temp = h.pA[i];
      for(int k = (i + 1); k < nRows; k++)
        temp -= b.pA[k]*G.pA[i][k]; 
      b.pA[i] = temp/G.pA[i][i];
    }
    else{
      cout << "Singular system!" << endl;
      exit(1);
    }
  }
  return b;
}

// Utilities:

void Matrix::printMatrix() const{
  for(int i = 0; i < nRows; i++){
    for(int j = 0; j < nCols; j++)
      cout << " " << pA[i][j] << "  ";
    cout << endl;
  }
}

//pdBand.h
#ifndef PDBAND_H
#define PDBAND_H
#include <iostream>
#include "simpleMatrix.h";
#include "vector.h"

class pdBand : public Matrix{
  int q;//the bandwidth
  Vector bandBack(const Vector& RHS) const;
  Vector bandForward(const Vector& RHS) const;
  Vector bandChol(const Vector& RHS) const;

 public:
 
  pdBand(int nrows = 0, int Q = -1, double a = 0.);
  pdBand(int nrows, int Q, const double* const* pa);
  pdBand(const pdBand& A);
 
  ~pdBand(){};

  pdBand& operator=(const pdBand& B);
 
  Vector solve(const Vector& RHS) const {
    std::cout << "Solution via banded Cholesky from class pdBand" 
	      << std::endl;  
    return bandChol(RHS);
  }
 
  Matrix bandToFull() const;
};
#endif

//pdBand.cpp
#include <cstdlib>
#include <cmath>
#include "pdBand.h"

using std::cout; using std::endl;
 
pdBand::pdBand(int nrows, int Q, const double* const* pa) 
  : Matrix(nrows, Q + 1, pa) {
  q = Q;
}
 
pdBand::pdBand(int nrows, int Q, double a) 
  : Matrix(nrows, Q + 1, a) {
  q = Q;
}
 
pdBand::pdBand(const pdBand& A) 
  : Matrix(A) {
  q = A.q;
}

pdBand& pdBand::operator=(const pdBand& A){

  if(this == &A)//avoid self assignment
    return *this;
  this->Matrix::operator=(A);
  q = A.q;

  return *this;
}
   
//back-solve a banded  upper triangular system 
Vector pdBand::bandBack(const Vector& RHS) const{
  double temp;//temporary storage
  int up;
  int nrows = this->getNRows();
  Vector b(nrows, 0.);//solution vector

  //initialize recursion on last row of right-hand-side
  if((*this)[nrows - 1][0] != 0)
    b.pA[nrows - 1] = RHS.pA[nrows - 1]/(*this)[nrows - 1][0];
  else{
    cout << "Singular system!" << endl;
    exit(1);
  }
  
  //now work through the remaining rows
  for(int i = (nrows - 2); i >= 0; i--){
    if((*this)[i][0] != 0){
      temp = RHS.pA[i];
      up = std::min(nrows, q + 1);
      for(int k = 1; k < up; k++)
        temp -= b.pA[i+k]*(*this)[i][k]; 
      b.pA[i] = temp/(*this)[i][0];
    }
    else{
      cout << "Singular system!" << endl;
      exit(1);
    }
  }
  return b;
}

//forward-solve a banded lower triangular system 
Vector pdBand::bandForward(const Vector& RHS) const{
  double temp;//temporary storage
  int low;
  int nrows = this->getNRows();
  Vector b(nrows, 0.);//solution vector

  //initialize recursion on first row of right-hand-side
  if(*this[0][0] != 0){
    b.pA[0] = RHS.pA[0]/(*this)[0][0];
  }
  else{
    cout << "Singular system!" << endl;
    exit(1);
  }
  
  //now work through the remaining rows
  for(int i = 1; i < nrows; i++){
    if((*this)[i][0] != 0){
      temp = RHS.pA[i];
      low = std::max(0, i - q);
      for(int k = low; k < i; k++)
        temp -= b.pA[k]*(*this)[k][i - k]; 
      b.pA[i] = temp/(*this)[i][0];
    }
    else{
      cout << "Singular system!" << endl;
      exit(1);
    }
  }
  return b;
}

//solve a banded system via Cholesky factorization
Vector pdBand::bandChol(const Vector& RHS) const{
  double temp = 0;
  int low, up;
  int nrows = this->getNRows();
  
  //temporary storage
  double** pG = new double*[nrows];
  for(int i = 0; i < nrows; i++){
    pG[i] = new double[q + 1];
    for(int j = 0; j < q + 1; j++){
      pG[i][j] = (*this)[i][j];
    }
  }
  for(int i = 0; i < nrows; i++){
    if((*this)[i][0] == 0){
      cout << "Singular system!" << endl;
      exit(1);
    }
    temp = (*this)[i][0];
    
    low = std::max(0, i - q);
    for(int k = low; k < i; k++)
      temp -= pG[k][i - k]*pG[k][i - k];
      
    pG[i][0] = sqrt(temp);

    up = std::min(i + q, nrows - 1);
      
    for(int j = (i + 1); j <= up; j++){
      temp = (*this)[i][j - i];
      low = std::max(0, j - q);
      for(int k = low; k < i; k++)
        temp -= pG[k][i - k]*pG[k][j - k];
      pG[i][j - i] = temp/pG[i][0];
    }
  }
  pdBand G(nrows, q, pG);
  Vector h = G.bandForward(RHS);
  Vector b = G.bandBack(h);

  for(int i = 0; i < nrows; i++)
    delete[] pG[i];
  delete[] pG;
  return b;
}    

Matrix pdBand::bandToFull() const{
  int nrows = this->getNRows();
  double** pTemp = new double*[nrows];
  for(int i = 0; i < nrows; i++){
    pTemp[i] = new double[nrows];
    for(int j = i; j < nrows; j++){
      if(j <= std::min(nrows, i + q))
	pTemp[i][j] = (*this)[i][j - i];
      else pTemp[i][j] = 0.;
    }
  }
  for(int i = 1; i < nrows; i ++)
    for(int j = 0; j < i; j++)
      pTemp[i][j] = pTemp[j][i];
  Matrix A(nrows, nrows, pTemp);
  for(int i = 0; i < nrows; i++)
    delete[] pTemp[i];
  delete[] pTemp;
  return A;
}

//pdBandDriver.cpp 
#include <iostream>
#include <fstream>
#include "simpleMatrix.h"
#include "vector.h"
#include "pdBand.h"

int main(int argc, char* argv[]){
  double* pc = new double[6];
  double** pA = new double*[6];

  std::ifstream fIn1, fIn2;
  fIn1.open(argv[1]); fIn2.open(argv[2]);

  for(int i = 0; i < 6; i++){
    pA[i] = new double[3];
    for(int j = 0; j < 3; j++){
      fIn1 >> pA[i][j];
    }
    fIn2 >> pc[i];
  }       
  fIn1.close(); fIn2.close();
  
  pdBand A(6, 2, pA);
  Vector c(6, pc);
  Matrix* pMat = &A;
  pMat->solve(c).printVec();

  Matrix B = A.bandToFull();
  pMat = &B;
  pMat->solve(c).printVec();
  return 0;
}
///////////////////////////////////////////////////////////////////////
//tntEx1.cpp 
#include "tnt.h"
#include <iostream>
#include <fstream>

using std::cout; using std::endl;

int main(int argc, char** argv){
  int m = atoi(argv[2]);
  int n = atoi(argv[3]);
  double* pc = new double[n];
  double* pA = new double[m*n];
  std::ifstream fIn;
  fIn.open(argv[1]);

  for(int i = 0; i < m; i++){
    for(int j = 0; j < n; j++)
      fIn >> pA[n*i + j];//row dominant form
    fIn >> pc[i];
  }
  fIn.close();
  TNT::Array1D<double> c(n, pc);
  cout << "The vector c" << endl;
  cout << c << endl;
  TNT::Array2D<double> A(m, n, pA);
  cout << "The matrix A" << endl;
  cout << A << endl;
  
  TNT::Array2D<double> B(m, n);
  for(int i = 0; i < m; i++){
    for(int j = 0; j < n; j++)
      B[i][j] = (double)(i +1)*(j + 1);
  }
  cout << "The matrix B" << endl;
  cout << B << endl;

  cout << "The sum of A and B" << endl;
  cout << A + B << endl;
  cout << "The product of A and B" << endl;
  cout << matmult(A, B) << endl;
  return 0;
}
///////////////////////////////////////////////////////////////////////
//tntEx2.cpp 
#include <iostream>
#include <fstream>
#include <algorithm>
#include "tnt.h"
#include "jama_cholesky.h"
#include "jama_eig.h"
#include "jama_svd.h"

using std::cout; using std::endl;

int main(int argc, char** argv){
  int m = atoi(argv[1]);
  int n = atoi(argv[2]);
  double* pA1 = new double[n*n];
  double* pc = new double[n];
  double* pA2 = new double[n*m];
  std::ifstream fIn1, fIn2;
  fIn1.open(argv[3]);
  fIn2.open(argv[4]);

  for(int i = 0; i < n; i++){
    for(int j = 0; j < n; j++)
      fIn1 >> pA1[n*i + j];//row dominant form
    fIn1 >> pc[i];
  }
 for(int i = 0; i < n; i++)
    for(int j = 0; j < m; j++)
      fIn2 >> pA2[m*i + j];

  fIn1.close(); fIn2.close();
  TNT::Array2D<double> A1(n, n, pA1);
  cout << "The matrix A1" << endl;
  cout << A1 << endl;
  TNT::Array1D<double> c(n, pc);
  cout << "The right-hand side" << endl;
  cout << c << endl;
  JAMA::Cholesky<double> Ch(A1);
  cout << "The solution by Cholesky" << endl;
  cout << Ch.solve(c) << endl;
  
  JAMA::Eigenvalue<double> A1Eig(A1);
  TNT::Array2D<double> V1(n, n);
  TNT::Array1D<double> lambda(n);
  A1Eig.getRealEigenvalues(lambda);
  A1Eig.getV(V1);
  cout << "The eigenvalues for A1" << endl;
  cout << lambda << endl;
  cout << "The eigenvectors for A1" << endl;
  cout << V1 << endl;

  TNT::Array2D<double> A2(n, m, pA2);
  cout << "The matrix A2" << endl;  
  cout << A2 << endl;

  int dMin = min(m, n);
  JAMA::SVD<double> A2Svd(A2);
  TNT::Array2D<double> V2(m, m);
  TNT::Array2D<double> U(n, n);
  TNT::Array1D<double> s(dMin);
  A2Svd.getSingularValues(s);
  A2Svd.getU(U);
  A2Svd.getV(V2);

  cout << "The singular values for A2" << endl;
  cout << s << endl;
  cout << "The left singular vectors for A2" << endl;
  cout << U << endl;
  cout << "The right singular vectors for A2" << endl;
  cout << V2 << endl;

  return 0;
}
///////////////////////////////////////////////////////////////////////
//Additional drivers for the choesky/gauss, eigenvalue, svd and qr examples
//cholGaussEx.cpp 
#include "matrix.h"
#include "vector.h"
#include <iostream>
#include <fstream>

using std::cout; using std::endl; using std::ifstream;


int main(int argc, char** argv){
  double* pc = new double[5];
  double** pA = new double*[5];
  ifstream fIn;
  fIn.open("exampleMatrices.txt");
  for(int i = 0; i < 5; i++){
    pA[i] = new double[5];
    for(int j = 0; j < 5; j++)
      fIn >> pA[i][j];
    fIn >> pc[i];
  }
  fIn.close();
  cout << "Here is A" << endl;
  Matrix A(5, 5, pA);
  A.printMatrix();
    
  Vector c(5, pc);
  cout << "Here is c" << endl;
  c.printVec();
    
  Vector b = A.gauss(c);

  cout << "Gauss solution" << endl;
  b.printVec();
  Matrix G(5, 5, 0.);
  b = A.cholesky(c);
  cout << "Cholesky solution" << endl;
  b.printVec();
    return 0;
}

//////////

//eigSvdEx.cpp 
#include "matrix.h"
#include "vector.h"
#include <iostream>
#include <fstream>

using std::cout; using std::endl; using std::ifstream;

int main(int argc, char** argv){
  double* pc = new double[5];
  double** pA = new double*[5];
  ifstream fIn;
  fIn.open("exampleMatrices.txt");
  for(int i = 0; i < 5; i++){
    pA[i] = new double[5];
    for(int j = 0; j < 5; j++)
      fIn >> pA[i][j];
    fIn >> pc[i];
  }
 
  cout << "Here is A" << endl;
  Matrix A(5, 5, pA);
  A.printMatrix();
    
  Vector v(5, pc);
  cout << "Here is v" << endl;
  v.printVec();
  Matrix U(5, 5, 0.);
  Vector lam = A.eigen(5, v, U, 100, .0000000001);
  cout << "Eigenvalues and vectors for A" << endl;
  lam.printVec();
  cout << endl;
  U.printMatrix();

  for(int i = 0; i < 5; i++){
    pA[i] = new double[4];
    for(int j = 0; j < 4; j++)
      fIn >> pA[i][j];
  }
  fIn.close();
  
  A = Matrix(5, 4, pA);
  cout << "Matrix for the svd" << endl;
  A.printMatrix();
  cout << "Vector for the svd" << endl;
  v = Vector(4, pc);
  v.printVec();
  U = Matrix(5, 4, 0.);
  Matrix V(4, 4, 0.);
  lam = A.SVD(4, v, U, V, 100, .0000000001);
  cout << endl;
  cout << "Singular values for A" << endl;
  lam.printVec();
  cout << endl;
  cout << "Right singular vectors" << endl;
  V.printMatrix();
  cout << endl;
  cout << "Left singular vectors" << endl;
  U.printMatrix();
}

//////////

//svdLS.cpp 
#include "matrix.h"
#include "vector.h"
#include <iostream>
#include <fstream>

using std::cout; using std::endl; using std::ifstream;

int main(int argc, char** argv){

  double** pA = new double*[6];
  ifstream fIn;
  double* pV = new double[3];
  double* pC = new double[6];

  ifstream fin;
  fIn.open("exampleMatrices.txt");
  for(int i = 0; i < 15; i++)
    fIn.ignore(255, '\n');
 
  for(int i = 0; i < 6; i++){
    pA[i] = new double[3];
    for(int j = 0; j < 3; j++)
      fIn >> pA[i][j];
    fIn >> pC[i];
  }
  for(int i = 0; i < 3; i++)
    fIn >> pV[i];
  fin.close();

  Matrix A(6, 3, pA);
  cout << "Matrix for the svd" << endl;
  A.printMatrix();
  cout << "Vector to start power method recursion" << endl;
  Vector v(3, pV);
  v.printVec();
  cout << "Right hand side vector" << endl;
  Vector c(6, pC);
  c.printVec();

  Matrix U = Matrix(6, 2, 0.);
  Matrix V(3, 2, 0.);
  Vector lam = A.SVD(2, v, U, V, 100, .0000000001);

  cout << "Singular values for A" << endl;
  lam.printVec();
  cout << endl;
  cout << "Left singular vectors" << endl;
  V.printMatrix();
  cout << endl;
  cout << "Right singular vectors" << endl;
  U.printMatrix();

  double** pLamInv = new double*[2];
  for(int i = 0; i < 2; i++){
    pLamInv[i] = new double[2];
    for(int j = 0; j <2; j++)
      pLamInv[i][j] = 0.;
  }
  for(int i = 0; i < 2; i++)
    pLamInv[i][i] = 1./lam[i];
  Matrix LamInv(2, 2, pLamInv);
  Matrix Ainv = V*LamInv*U.trans();
  cout << "Moore-Penrose inverse of A" << endl;
  Ainv.printMatrix();

  cout << "Least square solution" << endl;
  Vector b = Ainv*c;
  b.printVec();//least squares solution
}

//////////

//qrEx.cpp
#include "matrix.h"
#include "vector.h"
#include <iostream>
#include <fstream>

using std::cout; using std::endl; using std::ifstream;

int main(int argc, char** argv){
  double** pA = new double*[5];
  double* pc = new double[5];


  std::ifstream fIn;
  fIn.open("exampleMatrices.txt");
  for(int i = 0; i < 10; i++)
    fIn.ignore(255, '\n');
 
  for(int i = 0; i < 5; i++){
    pA[i] = new double[4];
    for(int j = 0; j < 4; j++){
      fIn >> pA[i][j];
    }
    fIn >> pc[i];
  }
  fIn.close();
  Matrix A(5, 4, pA);
  cout << "The A matrix" << endl;
  A.printMatrix();
  Vector c(5, pc);
  cout << "The y vector" << endl;
  c.printVec();

  Matrix Q(A.getnRows(), A.getnCols());
  Matrix R(A.getnCols(), A.getnCols());
  A.QR(Q, R);
  c = Q.trans()*c;
  Vector b = R.backward(c);
  cout << "Least squares solution" << endl;
  b.printVec();
}

//////////
//exampleMatrices.txt
/*
2.81837298736488 -3.14525367334217 -0.208262669293409 -4.19317243555944 -1.08278242271503 -1.68669331074241
-3.14525367334217 5.42465301722828 1.92991794589062 5.51100748630534 0.145300484178578 0.837787044494525
-0.208262669293409 1.92991794589062 2.10967002326924 1.91872220821717 -1.52998036566192 0.153373117836515
-4.19317243555944 5.51100748630534 1.91872220821717 8.02392851378211 -0.21456814263257 -1.13813693701195
-1.08278242271503 0.145300484178578 -1.52998036566192 -0.21456814263257 3.1624021689251 1.25381492106993
-0.560475646552213 1.71506498688328 1.22408179743946 1.78691313680308
-0.23017748948328 0.460916205989202 0.359813827057364 0.49785047822924
1.55870831414912 -1.26506123460653 0.400771450594052 -1.96661715662964
0.070508391424576 -0.686852851893526 0.110682715945120 0.701355901563685
0.129287735160946 -0.445661970099958 -0.555841134754075 -0.472791407727934
-0.560475646552213 1.71506498688328 1.22408179743946 1.78691313680308 -1.68669331074241
-0.23017748948328 0.460916205989202 0.359813827057364 0.49785047822924 0.837787044494525
1.55870831414912 -1.26506123460653 0.400771450594052 -1.96661715662964 0.153373117836515
0.070508391424576 -0.686852851893526 0.110682715945120 0.701355901563685 -1.13813693701195
0.129287735160946 -0.445661970099958 -0.555841134754075 -0.472791407727934 1.25381492106993
1 1 0 -0.560475646552213
1 1 0 -0.23017748948328
1 1 0 1.55870831414912
1 0 1 0.070508391424576
1 0 1 0.129287735160946
1 0 1 1.71506498688328
0.460916205989202
-1.26506123460653
-0.686852851893526
*/