NRMat svdinverse()

mat.cc

@@ -3495,21 +3495,32 @@ copyonwrite();
 }
 
 
-
+//cannot be in mat.h - would be undefined for intereger types
 template<>
-NRMat<double> NRMat<double>::inverse()
+NRMat<double> NRMat<double>::inverse() const
 {
-NRMat<double> tmp(*this);
+return calcinverse(*this);
-return calcinverse(tmp);
 }
 
 template<>
-NRMat<std::complex<double> > NRMat<std::complex<double> >::inverse()
+NRMat<std::complex<double> > NRMat<std::complex<double> >::inverse() const
 {
-NRMat<std::complex<double> >  tmp(*this);
+return calcinverse(*this);
-return calcinverse(tmp);
 }       
 
+template<>
+NRMat<double> NRMat<double>::svdinverse(const double thr) const
+{
+return calcsvdinverse(*this,thr);
+}
+
+template<>
+NRMat<std::complex<double> > NRMat<std::complex<double> >::svdinverse(const double thr) const
+{
+return calcsvdinverse(*this,thr);
+}
+
+
 
 
 

mat.h

@@ -173,7 +173,10 @@ public:
 	void identity() {*this = (T)1;}
 
 	//! inverse matrix
-	NRMat inverse();
+	NRMat inverse() const;
+
+	//! pseudo-svd-inverse matrix
+        NRMat svdinverse(const LA_traits<T>::normtype thr=0) const;
 
 	//! add scalar value to the diagonal elements
 	NRMat & operator+=(const T &a);

nonclass.h

@@ -249,13 +249,11 @@ extern const typename LA_traits<T>::complextype complexmatrix (const T&, const T
 extern void cholesky(NRMat<double> &a, bool upper=1);
 extern void cholesky(NRMat<std::complex<double> > &a, bool upper=1);
 
-//inverse by means of linear solve, preserving rhs intact
+//inverse by means of linear solve, pass by value to preserve argument intact
 template<typename T>
 const NRMat<T> calcinverse(NRMat<T> a, T *det=NULL)
 {
-#ifdef DEBUG
 	if(a.nrows()!=a.ncols()) laerror("inverse() for non-square matrix");
-#endif
 	NRMat<T> result(a.nrows(),a.nrows());
 	result = (T)1.;
 	a.copyonwrite();
@@ -264,6 +262,24 @@ const NRMat<T> calcinverse(NRMat<T> a, T *det=NULL)
 	return result;
 }
 
+
+//inverse by means of SVD , pass by value to preserve argument intact
+template<typename T>
+const NRMat<T> calcsvdinverse(NRMat<T> a, double thr=0)
+{
+        if(a.nrows()!=a.ncols()) laerror("svdinverse() for non-square matrix");
+	int n=a.nrows();
+        a.copyonwrite();
+	NRMat<T> u(n,n),v(n,n);
+	NRVec<double> w(n);
+	singular_decomposition(a,&u,w,&v,true);
+	for(int i=0; i<n; ++i) w[i] = (w[i]<thr)? 0. : 1./w[i];
+	v.diagmultr(w);
+	u.transposeme();
+        return v*u; //could use gemm instead of separate transpose too
+}
+
+
 //several matrix norms
 template<class MAT>
 typename LA_traits<MAT>::normtype MatrixNorm(const MAT  &A, const char norm);

t.cc

@@ -4798,7 +4798,7 @@ cout<<"part\n"<<tt.subtensor1(1);
 cout<<"part\n"<<tt.subtensor1(2);
 }
 
-if(1)
+if(0)
 {
 int n;
 cin >>n;
@@ -4809,4 +4809,18 @@ if(p!=q) laerror("inverseme failed");
 else cout <<p<<"OK\n";
 }
 
+if(1)
+{
+int n;
+cin >>n;
+NRMat<double> a(n,n);
+a.randomize(1.);
+NRMat<double> b=a.inverse();
+NRMat<double> c=a.svdinverse(1e-14);
+cout<< "inverses diff = "<<(b-c).norm()<<endl;
+cout<< "inverse error = "<<(a*b).norm(1.)<<endl;
+cout<< "svdinverse error = "<<(a*c).norm(1.)<<endl;
+}
+
+
 }//main