adjoint matrix

mat.cc

@@ -3511,13 +3511,28 @@ return calcinverse(*this);
 template<>
 NRMat<double> NRMat<double>::svdinverse(const double thr) const
 {
-return calcsvdinverse(*this,thr);
+return calcsvdinverse(*this,thr,false);
 }
 
 template<>
 NRMat<std::complex<double> > NRMat<std::complex<double> >::svdinverse(const double thr) const
 {
-return calcsvdinverse(*this,thr);
+return calcsvdinverse(*this,thr,false);
+}
+
+
+
+
+template<>
+NRMat<double> NRMat<double>::absadjoint(const double thr) const
+{
+return calcsvdinverse(*this,thr,true);
+}
+
+template<>
+NRMat<std::complex<double> > NRMat<std::complex<double> >::absadjoint(const double thr) const
+{
+return calcsvdinverse(*this,thr,true);
 }
 
 

mat.h

@@ -178,6 +178,9 @@ public:
 	//! pseudo-svd-inverse matrix
         NRMat svdinverse(const LA_traits<T>::normtype thr=0) const;
 
+	//! adjoint matrix (transposed cofactor matrix)
+	NRMat absadjoint(const LA_traits<T>::normtype thr=0) const;
+
 	//! add scalar value to the diagonal elements
 	NRMat & operator+=(const T &a);
 	//! subtract scalar value to the diagonal elements

nonclass.h

@@ -265,16 +265,25 @@ const NRMat<T> calcinverse(NRMat<T> a, T *det=NULL)
 
 //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)
+const NRMat<T> calcsvdinverse(NRMat<T> a, double thr=0, bool do_absadjoint=false)
 {
         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);
+	NRVec<double> w(n),s(n);
 	singular_decomposition(a,&u,w,&v,true);
-	for(int i=0; i<n; ++i) w[i] = (w[i]<thr)? 0. : 1./w[i];
+	if(do_absadjoint) //absadjoint matrix - missing det(U)*det(V) to be full adjoint - a sign or complex phase
-	v.diagmultr(w);
+		{
+		for(int i=0; i<n; ++i)
+			{
+			s[i] = 1;
+			for(int j=0; j<n; ++j) if(j!=i) s[i] *= w[j];
+			}
+		}
+	else //inverse
+		for(int i=0; i<n; ++i) s[i] = (w[i]<thr)? 0. : 1./w[i];
+	v.diagmultr(s);
 	u.transposeme();
         return v*u; //could use gemm instead of separate transpose too
 }

t.cc

@@ -4820,6 +4820,14 @@ 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;
+
+NRMat<double> aa(a);
+double det=determinant_destroy(aa);
+
+NRMat<double> adj=a.absadjoint();
+cout <<"det = "<<det<<endl;
+cout <<"adjoint error = "<<(a*adj).norm(abs(det))<<endl;
+
 }