//  svd_test.edp
//
//  Discussion:
//
//    Demonstrate the interface between FreeFem++ and LAPACK by calling DGESVD.
//
//  Location:
//
//    http://people.sc.fsu.edu/~jburkardt/freefem_src/svd_test/svd_test.edp
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    21 July 2015
//
//  Author:
//
//    John Burkardt
//
cout << "\n";
cout << "svd_test:\n";
cout << "  FreeFem++ version\n";
cout << "  Demonstrate the interface to LAPACK by accessing\n";
cout << "  the DGESVD routine which computes the singular value decomposition.\n";

load "lapack"
//
//  Here's a function which we don't actually use in this example.
//  But it demonstrates how one can take a sparse matrix "Asparse"
//  as stored in FreeFem++ and copy it into a dense matrix "Full"
//  suitable for working with LAPACK.
//
func real[int,int]  arrayFromsparse ( matrix & Asparse )
{
  real[int,int] Full(Asparse.n,Asparse.m);
  Full = 0.0;
  int[int] I(1);
  int[int] J(1);
  real[int] C(1);
  [I,J,C] = Asparse;
  for ( int i = 0; i < I.n; ++i )
  {
    Full ( I(i), J(i) ) = C(i);
  }
  return Full;	
}
//
//  Here's a function we will need.   It wants to turn the vector sigma into an array.
//  The dimensions of this array are determined by the column size of U and the row
//  size of V.
//
func real[int,int] arrayFromRectDiag ( real[int,int] & U, 
  real[int] & sigma ,real[int,int] & V )
{
  real[int,int] res(U.n,V.m);
  res = 0.0;
  for ( int i = 0; i < sigma.n; ++i )
  {
    res(i,i) = sigma(i);
  }
  return res;
}
//
//  Here is a small dense example matrix C for which we might want the
//  singular value decomposition.
//
real[int,int] C =  [ 
  [  5.0e-16,         -1.863389981e-16,  0.894427191,     0.4472135955 ],
  [  2.916666667e-16,  1.490711985e-16, -5.123475383e-16, 0.1224361692 ],
  [ -0.0,              0.0,              0.0,             0.0 ] ];

cout << "\n";
cout << "  The C array:" << C << endl;
//
//  Declare arrays for SVD decomposition C = U * SIGMA * V'
//
real[int,int] U(1,1);
real[int,int] V(1,1);
real[int] sigma;
//
//  Call the LAPACK routine DGESVD for the SVD decomposition.
//  "dgesdd" is a special Freefem++-to-LAPACK interface.
//
dgesdd ( C, U, sigma, V );
//
//  The most interesting item is the vector SIGMA of singular values.
// 
cout << "  Singular value vector sigma:" << sigma << endl;
//
//  By the way, the call to dgesdd overwrote the original information in C.
//
cout << "\n";
cout << "Let's check that original C == U * sigma * V'" << endl;
//
//  Call "arrayFromRectDiag" to build a matrix from the vector sigma.
//  The arrays U and V are only needed for dimension information.
//
real[int,int] S = arrayFromRectDiag ( U, sigma, V );
//
//  Postmultiply U by SIGMA.
//
real[int,int] US = U * S;
//
//  Transpose V.
//
real[int,int] VT = V';
//
//  Postmultiply U*SIGMA by V'
//
real[int,int] USVT = US * VT;
//
//  We should have recovered C.
//
cout << "\n";
cout << " U " << U << endl;
cout << " V' " << VT << endl;
cout << " S =  " << S << endl;			
cout << " U * S * V'  = " << USVT << endl;
//
//  Terminate.
//
cout << "\n";
cout << "svd_test:\n";
cout << "  Normal end of execution.\n";

exit ( 0 );