/*
  CLAPACK demo to show how one might compute SVD.
  This code is not definitive, as there is more than one way to skin this cat.
  It is geared to a matrix A with less rows than columns (M<N).

  For Cholesky decomposition, check out spotrf.

  J.K. Johnstone, February 2007
*/

#include <iostream>
using namespace std;

// to avoid C++ incompatibilities (after all, it is expecting C)
extern "C" {               // make sure these are visible in directory lookup
#include "f2c.h"           // our solution: copy them to /usr/include
#include "clapack.h"
}

int main (int argc, char **argv)
{
  // computing the SVD A = U \sigma V^t, A is an MxN matrix
  char JOBU = 'A'; // 'A' ==> compute all M columns of U 
                   // (see manpages for sgesvd for other options)
  char JOBVT = 'A'; // 'A' ==> compute all N rows of V^t
  integer M=2, N=2;
  float *A; // matrix that will undergo SVD, as a 1d array in column-major order
  A = new float[M*N]; // M*N 
  // GvL p. 71 example, in column-major order
  A[0] = .96; A[1] = 2.28; A[2] = 1.72; A[3] = .96;
  cout << "A: "; for (int k=0; k<4; k++) cout << A[k] << ",";
  cout << endl;
  integer LDA=M; // 'leading dimension of A'
  float *S;  // singular values (nonincreasing)
  S = new float[M]; // min(M,N)
  float *U;  // U in the SVD
  U = new float[M*M]; // M*M
  integer LDU = M; // leading dimension of U (see manpages for restrictions)
  float *VT;  // V^t in the SVD
  VT = new float[N*N]; // N*N
  integer LDVT = N; // leading dimension of V^t
  integer LWORK=10;  // dimension of WORK
  float *WORK;  // workspace (man page recommends at least 5*min(M,N); 
                // see man page for more)
  WORK = new float[LWORK];
  integer INFO;  // 0 iff successful exit (o.w. indicates source of problems)

  sgesvd_ (&JOBU, &JOBVT, &M, &N, A, &LDA, S, U, &LDU, VT, &LDVT, 
	   WORK, &LWORK, &INFO);

  // note: A's contents are destroyed upon exit, since JOBU = JOBVT = 'A'
  cout << "INFO = " << INFO << endl;
  cout << "WORK[0] = " << WORK[0] << endl;
  if (INFO>=0)
    {
      int i,j;
      cout << "U = " << endl;
      for (i=0; i<M; i++)
	{
	  for (j=0; j<M; j++)
	    cout << U[j*M + i] << " "; // column-major order!
	  cout << endl;
	}
      cout << "singular values = ";
      for (i=0; i<2; i++)
	cout << S[i] << " ";
      cout << endl;
      
      cout << "V^t = " << endl;
      for (i=0; i<N; i++)
	{
	  for (j=0; j<N; j++)
	    cout << VT[j*N + i] << " "; // column-major order!
	  cout << endl;
	}
    }
  return 0;
}