From 3d284d544af76a58970737134b14453a43b33301 Mon Sep 17 00:00:00 2001
From: Jiri Pittner <jiri@pittnerovi.com>
Date: Tue, 9 Dec 2025 15:21:18 +0100
Subject: [PATCH] diagonalize3 wrapper for dsteqr

---
 nonclass.cc | 24 ++++++++++++++++++++++++
 nonclass.h  |  3 +++
 2 files changed, 27 insertions(+)

diff --git a/nonclass.cc b/nonclass.cc
index c9a3a48..1869614 100644
--- a/nonclass.cc
+++ b/nonclass.cc
@@ -944,7 +944,31 @@ extern "C" void FORNAME(zggev)(const char *JOBVL, const char *JOBVR, const FINT
                 std::complex<double> *VL, const FINT *LDVL,  std::complex<double> *VR, const FINT *LDVR,
                 std::complex<double> *WORK, const FINT *LWORK, double *RWORK, FINT *INFO );
 
+//tridiagonal
+extern "C" void FORNAME(dsteqr)(const char *compz, const  FINT *n, double *d, double *d1, double *z, const FINT *ldz,  double *work, FINT *info);
 
+void diagonalize3(NRVec<double> &d, NRVec<double> &d1, NRMat<double> *v, const bool corder)
+{
+	FINT n = d.size();
+	if(d1.size()!=n) laerror("inconsistent dimensions in diagonalize3");
+	if(v) {if(n!=v->nrows()||n!=v->ncols()) laerror("inconsistent dimensions in diagonalize3");}
+	d.copyonwrite();
+	d1.copyonwrite();
+	if(v) v->copyonwrite();
+	
+	FINT ldz=n;
+	FINT info;
+	int worksize=2*n-2;
+	if(worksize<1) worksize=1;
+        double *work = new double[worksize];
+	char job=v?'i':'n';
+	FORNAME(dsteqr)(&job,&n,&d[0],&d1[0],(v?&(*v)(0,0):NULL),&ldz,work,&info);
+	delete[] work;
+	if (v && corder) v->transposeme();
+
+	if (info < 0) laerror("illegal argument in dsteqr() fo diagonalize3");
+	if (info > 0) laerror("convergence problem in dsteqr() fo diagonalize3");
+}
 
 
 //statics for sorting
diff --git a/nonclass.h b/nonclass.h
index b09c7bf..3a3da84 100644
--- a/nonclass.h
+++ b/nonclass.h
@@ -183,6 +183,9 @@ extern void gdiagonalize(NRMat<std::complex<double> > &a, NRVec<double> &wr, NRV
                 NRMat<std::complex<double> > *vl=NULL, NRMat<std::complex<double> > *vr=NULL, const bool corder=1, int n=0, const int sorttype=0, const int biorthonormalize=0,
                 NRMat<std::complex<double> > *b=NULL, NRVec<std::complex<double> > *beta=NULL);
 
+//diagonalization of symmetric real tridiagonal matrix
+extern void diagonalize3(NRVec<double> &d, NRVec<double> &d1, NRMat<double> *v, const bool corder=1);
+
 //complex,real,imaginary parts of various entities
 template<typename T>
 extern const typename LA_traits<T>::realtype realpart(const T&);