// Gnurbs - A curve and surface library
// Copyright (C) 2008-2026 Eric Bechet
//
// See the LICENSE file for contributions and license information.
// Please report all bugs and problems to <bechet@cadxfem.org>.
//

#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();

// 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();
}