//
//  MatrixSVD.h
//
//  Qdos - Class for performing a Singular Value Decomposition of a Numeric
//  Matrix
//
//  Janet Barnes
//  Sept 1999
//

#ifndef QDOS_MATRIXSVD_H
#define QDOS_MATRIXSVD_H

#include "Matrix.h"
#include <cfloat>
#ifdef HAVE_SSTREAM_H
#include <sstream>
#else
#include <strstream>
#endif

#define EPSILON 1e-12

using namespace std;

namespace QDOS{
	void svd_golub_reinsch(double * a, int m, int n, double * q, double * u, double * v);
	
	void svd_golub_reinsch(double_complex * a, int m, int n, double * q, double_complex * u, double_complex * v);

// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ class MatrixSVD
            
// Class for performing a Singular Value Decomposition of a rectangular
// matrix
	template<class T>
			class MatrixSVD
	{
		typedef typename promote_trait<T,double>::T_promote T_promote;
		typedef T value_type;
    
		size_t _nrows;                 // No. rows and cols of rank 2 matrix
		size_t _ncols; 
		size_t _mindim; 
		Matrix<T_promote> _u;          // _nrows*_nrows orthogonal matrix _u
		Matrix<T_promote> _v;          // _ncols*_ncols orthogonal matrix _v
		Matrix<double> _w;             // Rank 1 matrix of min(_nrows,_ncols) 
                                   // singular values
		double _zeroval;                                        
                                       
		public:
   // --- Constructors
   // Empty constructor
			MatrixSVD(){;}
   // Decompose Matrix m, of M rows and N columns, M>=N
			MatrixSVD(const Matrix<T>& A,double t=0):_zeroval(t)	        
			{
				try
				{
       // verify that matrix argument IS rank 2
					if (A.rank()!=2)
						throw MatrixDomainError();
					_nrows=A.dims()[0];
					_ncols=A.dims()[1];           
#ifdef DIAG
					cout << "Resizing u,v and w" << endl;
#endif
					_u.resize(Dim(_nrows,_nrows,ROW_MAJOR_ORDER));
					_v.resize(Dim(_ncols,_ncols,ROW_MAJOR_ORDER));
					_mindim=std::min(_nrows,_ncols);
					_w.resize(Dim(_mindim,ROW_MAJOR_ORDER));
					int m, n;
					m=_nrows;
					n=_ncols;
           // take copy of array and store in col-major-order
#ifdef DIAG 
					cout << "Taking copy of array and storing in row-major-order" << endl; 
#endif
					Matrix<T_promote> B(Dim(A.dims(),ROW_MAJOR_ORDER));
					B.copy(A);
           // call CLAPACK routine
#ifdef DIAG
					cout << "calling SVD routine" << endl;
#endif
					svd_golub_reinsch(B.ptr(), m, n, _w.ptr(), _u.ptr(), _v.ptr());
           // ... test result of info....
				}
				catch(Exception& e)
				{
					e.msg() << "[MatrixSVD<T>]"; 
					throw e;
				}
     
			}

   // --- Inquiries
			double& zeroval() { return _zeroval;}
			size_t get_rows() const { return _nrows;}
			size_t get_cols() const { return _ncols;}
			size_t get_mindim() const { return _mindim;}
			MatrixRes<T_promote> get_u() const { return MatrixRes<T_promote>(_u); }
			MatrixRes<T_promote> get_v() const { return MatrixRes<T_promote>(_v); }
			MatrixRes<double> get_w() const { return MatrixRes<double>(_w); }

   // --- Misc functions
			MatrixRes<T_promote> pinverse();
			MatrixRes<T_promote> solve(const Matrix<T>& B);
   
	};

// ++++++++++++++++++++++++++++++++++++++++++++++++++++++ class MatrixSVD::solve
// :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ----------------------------------------------------------------------- solve

	template<class T>
			inline MatrixRes<typename promote_trait<T,double>::T_promote> \
			MatrixSVD<T>::solve(const Matrix<T>& B) 
			{
				try
				{
					typedef typename promote_trait<T,double>::T_promote T_promote;
        // verify that B is rank 2 or 1
					if (B.rank()!=1 && B.rank()!=2) 
						throw MatrixRankError();
        // verify that B has correct no of rows
					if (B.dims()[0]!=get_rows())
						throw MatrixLengthError();
					bool zero_coeff;
					double sigi;
					double zero=0.;
					Matrix<T_promote> Ut;
					string str;
        // value against which to compare for determining zero values   
					if (zeroval()==0)
						zeroval()=get_mindim()*max(abs(get_w())) * DBL_EPSILON; // in <float.h> 
					zero_coeff=min(abs(get_w())) < zeroval();
					if (zero_coeff)
					{
						ostrstream ost;
						ost << "[MatrixSVD<T>::solve] AX=B : Linear dependency detected "\
								<< "among X" << endl;
						ost << "#  " << "Singular Value\t" << "Null Coefficients"<<endl;
						str+=ost.str();
					}
        // set Ut=transpose of U[0...m-1;0...n-1] i.e. m x n
					if (get_rows()>=get_cols()) 
					{
            // m > n
						Ut=get_u().col(slice(0,get_cols(),1)).shallowcopy(); 
					}
					else         
					{ 
						Ut.resize(Dim(get_rows(),get_cols()));
						Ut.col(slice(0,get_rows(),1))=get_u(); // n > m
					}
        // transpose to get n x m matrix
        //Ut.transpose(); 
					Matrix<T_promote> Vx=(!Ut)*B;         // Vx = n p (n m * m p)  B is rank 2
                                           // Vx = n (n m * m) B is rank 1                              
					Matrix<T_promote> V(get_v());                                   
        // Calculate U'B/sig
					for (int i=0;i<get_mindim();i++)
					{
						sigi=get_w()(i);
						if (sigi>zeroval())
							Vx[i]/=sigi; // Vx[i] will be a rank 1 or rank 0 Matrix
						else 
						{
							Vx[i]=zero;
							ostrstream ost;
							ost << i << setw(16) << sigi << "\t(";
							Matrix<T_promote> Vi=V.col(i);
							for (size_t j=0;j<Vi.size();j++)
								ost << setw(10) << Vi(j) << " ";  
							ost << ")" << endl;
							str+=ost.str();
						}  
					}
    
        // for cases where N > M i.e. mindim is M
					for (int i=get_mindim();i<get_cols();i++)
						Vx[i]=zero;
          
					if (str.size())
						cout << str << endl;
					return (V*Vx);
				}
				catch(Exception& e)
				{
					e.msg() << "[MatrixSVD<T>::solve]";
					throw e;
				}
            
			}


// +++++++++++++++++++++++++++++++++++++++++++++++++++ class MatrixSVD::pinverse
// :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// -------------------------------------------------------------------- pinverse

			template<class T>
					MatrixRes<typename promote_trait<T,double>::T_promote> MatrixSVD<T>::pinverse() 
			{
// return the pseudo-inverse of a singular-value-decomposed matrix
				try
				{
					size_t nrows=get_rows();
					Matrix<typename promote_trait<T,double>::T_promote> \
							I(Dim(nrows,nrows));
					I.unit_matrix();
					return solve(I);
				}
				catch(Exception& e)
				{
					e.msg() << "[MatrixSVD<T>::pinverse]";
					throw e;
				}
			}

}  // end namespace QDOS 
#endif
