#include "libscl.h"

using namespace std;
using namespace scl;

int main()
{
  realmat A(5,5), B;
  INTEGER i,j;
/*
  for (j=1; j<=5; j++) for (i=1; i<=5; i++)  A(i,j)=pow(REAL(i),j);
  B = A(" ","1:2"); 
  cout << A << B;
  j = solve(A,B,1.e-10);
  cout << "return value = " << j << '\n';
  cout << A << B;
 */
  for (j=1; j<=5; j++) for (i=1; i<=5; i++)  A(i,j)=pow(REAL(i),j);
  for (i=1; i<=5; i++)  A(i,3)=A(i,1);
  B = cbind(A,A(" ","1"));
  cout << A << B;
  j = solve(A,B);
  cout << "return value = " << j << '\n';
  cout << A << B;
  
  realmat C;

  cout << starbox("Using solve/_");
  for (j=1; j<=5; j++) for (i=1; i<=5; i++) A(i,j)=pow(REAL(i),j);
  for (i=1; i<=5; i++)  A(i,3)=A(i,1);
  C=A+T(A); A = C;
  B.resize(5,5,0.0);
  for (i=1; i<=5; i++) B(i,i)=1.0;
  cout << A << B;
  j = solve(A,B,1.e-10);
  cout << "return value = " << j << '\n';
  cout << A << B << C*B*C;

  cout << starbox("Using psdsol/_");
  for (j=1; j<=5; j++) for (i=1; i<=5; i++)  A(i,j)=pow(REAL(i),j);
  for (i=1; i<=5; i++)  A(i,3)=A(i,1);
  C=A+T(A); A = C;
  B.resize(5,5,0.0);
  for (i=1; i<=5; i++) B(i,i)=1.0;
  cout << A << B;
  j = psdsol(A,B,1.e-10);
  cout << "return value = " << j << '\n';
  cout << A << B << C*B*C;


  realmat XX, Xy;
  realmat X(5,2);
  realmat y(5,1);
  realmat b(2,1);

  for (i=1; i<=5; i++) {
    X(i,1)=1.0;
    X(i,2)=i;
  }

  b(1,1)=1.0;
  b(2,1)=2.0;

  y = X*b;
  XX = T(X)*X;
  Xy = T(X)*y;

  cout << starbox("Solve normal eqs with solve/_") << XX << Xy;
  j=solve(XX,Xy);
  cout << "return value = " << j << '\n';
  cout << XX << Xy << b;

  y = X*b;
  XX = T(X)*X;
  Xy = T(X)*y;

  cout << starbox("Solve normal eqs with psdsol/_") << XX << Xy;
  j=psdsol(XX,Xy);
  cout << "return value = " << j << '\n';
  cout << XX << Xy << b;

  exit(0);
}












