#include "libscl.h"

using namespace std;
using namespace scl;

int main()
{
  REAL zero = 0;
  REAL one = 1;

  INTEGER r=5, c=4;
  realmat A(r,c);
  realmat Aplus0, Aplus1;

  INTEGER i,j;
  INTEGER rank;

  for (j=1; j<=c; j++) {
    for (i=1; i<=r; i++) {
      A(i,j)=j+i;
    }
  }

  realmat U,S,V;

  cout << "Input matrix" << A << '\n';

  rank = svd(A,U,S,V);

  cout << "rank = " << rank << '\n';

//  cout << "Singular value decomposition\n" << U << S << V << '\n';

  for (i=1; i<=S.nrow(); i++)
    {if (i<=rank) {S[i] = one/S[i];} else {S[i]=zero;}}

  Aplus0 = V * diag(S) * T(U);

  cout << "Aplus from svd\n" << Aplus0 << '\n';

  Aplus1 = ginv(A);

  cout << "Aplus from ginv\n" << Aplus1 << '\n';

  cout << "difference\n" << Aplus1 - Aplus0 << '\n';


  A = T(A);

  cout << "Input matrix" << A << '\n';

  rank = svd(A,U,S,V);

  cout << "rank = " << rank << '\n';

//  cout << "Singular value decomposition\n" << U << S << V << '\n';

  for (i=1; i<=S.nrow(); i++)
    {if (i<=rank) {S[i] = one/S[i];} else {S[i]=zero;}}

  Aplus0 = V * diag(S) * T(U);

  cout << "Aplus from svd\n" << Aplus0 << '\n';

  Aplus1 = ginv(A);

  cout << "Aplus from ginv\n" << Aplus1 << '\n';

  cout << "difference\n" << Aplus1 - Aplus0 << '\n';


  A.resize(r,c,zero);

  cout << "Input matrix" << A << '\n';

  Aplus1 = ginv(A);

  cout << "Aplus from ginv\n" << Aplus1 << '\n';

  A.resize(r,r);

  static INT_32BIT ix; ix=983874;

  for (i=1; i<=r; i++){
  for (j=1; j<=r; j++){
    A(i,j) = (i+j)*unsk(&ix);
  }
  }

  cout << "Input matrix" << A << '\n';
  Aplus0 = inv(A);
  cout << "Aplus from inv\n" << Aplus0 << '\n';
  Aplus1 = ginv(A);
  cout << "Aplus from ginv\n" << Aplus1 << '\n';
  cout << "difference\n" << Aplus1 - Aplus0 << '\n';


  realmat D;
  cout << "Input matrix" << D << '\n';
  Aplus1 = ginv(D);

  return 0;

}


