#include <iostream>
#include "linear_algebra.h"

int main(void)
{

//  Solving AX=B

  Square_Matrix A(5);
  // A(L,C) refers to line L and column C of the matrix
  A(0,0) =1;
  A(0,1) =0.0;
  A(0,2) =0;
  A(0,3) =0;
  A(0,4) =0;
  A(1,0) =0;
  A(1,1) =1;
  A(1,2) =1;
  A(1,3) =1;
  A(1,4) =1;
  A(2,0) =0;
  A(2,1) =0;
  A(2,2) =0.5;
  A(2,3) =1;
  A(2,4) =1;
  A(3,0) =0;
  A(3,1) =0;
  A(3,2) =0;
  A(3,3) =0.25;
  A(3,4) =1;
  A(4,0) =0;
  A(4,1) =1;
  A(4,2) =0;
  A(4,3) =0;
  A(4,4) =0.5;
  Vector B(5),X(5);
  B(0) =1;
  B(1) =2;
  B(2) =-1;
  B(3) =0;
  B(4) =1;
  LU_Matrix LU(A);    // this does transfom A into an L*U (rather costly)
  LU.Solve_Linear_System(B,X);    // this does the back substitution to solve L*U*X=B (cheap)
  std::cout << "Solving AX=B" << std::endl;
  std::cout << "A=" ;
  A.Display();
  std::cout << "B=" ;
  B.Display();
  std::cout << "X=" ;
  X.Display();
  A.Mult(X,B);    // check if solution is ok
  std::cout << "AX=" ;
  B.Display();
  
  //eigenvalues and eigenvectors
  
  Square_Matrix AA(3),CC(3);
  Vector VV(3);

  AA(0,0) =2.0;
  AA(0,1) =0.0;
  AA(0,2) =0.0;
  AA(1,0) =0.0;
  AA(1,1) =3.0;
  AA(1,2) =4.0;
  AA(2,0) =0.0;
  AA(2,1) =4.0;
  AA(2,2) =9.0;
  std::cout << "AA=" ;
  AA.Display();
  int err=CC.Eigen_Vectors(AA,VV);
  if (!err)
  {
    std::cout << "Eigenvalues of AA : " ;
    VV.Display();
    std::cout << "Eigenvectors of AA : " ;
    CC.Display();
  }
  else std::cout << "Problem finding eigenvectors of AA" << std::endl;

  // manipulating rectangular matrices

  Rect_Matrix J(2,3);
  Rect_Matrix JT(3,2);
  J(0,0) =1;
  J(0,1) =0.0;
  J(0,2) =0;
  J(1,0) =0;
  J(1,1) =1;
  J(1,2) =1;
  for (int i=0;i<J.SizeL();++i)
    for (int j=0;j<J.SizeC();++j)
      JT(j,i) =J(i,j);
  A.Resize(2);    // A is square
  A.Mult(J,JT);   // A=J*JT
  std::cout << "J=" ;
  J.Display();
  std::cout << "JT=";
  JT.Display();
  std::cout << "A=J*JT=" ;
  A.Display();
  A.Square(J);   // A=J*JT
  std::cout << "A=square(J) = J*JT" ;
  A.Display();
  std::cout << "Solving AX=B" << std::endl;
  LU_Matrix LU2(A);
  LU2.Solve_Linear_System(B,X);    // B and X are too big, but that's OK. Only the first indices are used.
  std::cout << "B=";
  B.Display();
  std::cout << "X=";
  X.Display();
}