load "lapack"

func real[int,int]  arrayFromsparse(matrix & Asparse)
{
	real[int,int] Full(Asparse.n,Asparse.m);
	Full = 0.;
	int[int] I(1),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;	
}

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.;
	for(int i=0;i<sigma.n;++i)
	{
		res(i,i) = sigma(i);
	}
	return res;
}


real[int,int] C =  [ [5e-16, -1.863389981e-16, 0.894427191, 0.4472135955],
	 [2.916666667e-16, 1.490711985e-16, -5.123475383e-16, 0.1224361692],
	 [ -0,   0,   0,   0] ];


cout << "C before call to dgesdd " << C << endl;
real[int,int] U(1,1), V(1,1);
real[int] sigma;
dgesdd(C,U,sigma,V);
cout << " sigma " << sigma << endl;
cout << "C was modified " << C << endl;

cout << "Let's check that orginal C == U sigma V^T " << endl;

real[int,int] arrSigma = arrayFromRectDiag(U,sigma ,V);
real[int,int] arTmp = U*arrSigma;
real[int,int] VTT = V';
real[int,int] arTmp2 = arTmp*VTT;
cout << " U " << U << endl;
cout << " VT " << VTT << endl;
cout << " Sigma =  " << arrSigma << endl;			
cout << " U*  Sigma * V^T  = " << arTmp2 << endl;