From 837963193f0a78e8eb519aabdaba84792d642ffa Mon Sep 17 00:00:00 2001
From: Jiri Pittner <jiri@pittnerovi.com>
Date: Sat, 5 Sep 2026 18:15:39 +0200
Subject: [PATCH] adjoint matrix

---
 mat.cc     | 19 +++++++++++++++++--
 mat.h      |  3 +++
 nonclass.h | 17 +++++++++++++----
 t.cc       |  8 ++++++++
 4 files changed, 41 insertions(+), 6 deletions(-)

diff --git a/mat.cc b/mat.cc
index 80b6c0e..40c1762 100644
--- a/mat.cc
+++ b/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);
 }
 
 
diff --git a/mat.h b/mat.h
index 5c20355..0285e3f 100644
--- a/mat.h
+++ b/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
diff --git a/nonclass.h b/nonclass.h
index ac30e98..21d15cf 100644
--- a/nonclass.h
+++ b/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];
-	v.diagmultr(w);
+	if(do_absadjoint) //absadjoint matrix - missing det(U)*det(V) to be full adjoint - a sign or complex phase
+		{
+		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
 }
diff --git a/t.cc b/t.cc
index 6f96253..169a746 100644
--- a/t.cc
+++ b/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;
+
 }