From ea6547f9f602ee082cc184a5712746a0cd47adc5 Mon Sep 17 00:00:00 2001
From: jiri <jiri>
Date: Mon, 19 Oct 2009 19:38:57 +0000
Subject: [PATCH] *** empty log message ***

---
 mat.cc | 77 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 77 insertions(+)

diff --git a/mat.cc b/mat.cc
index 98e6757..013f89c 100644
--- a/mat.cc
+++ b/mat.cc
@@ -1167,7 +1167,80 @@ for (int i=0; i< nn; i++) v[i][i] = r[i];
 #endif
 }
 
+template<>
+void NRMat<double>::orthonormalize(const bool rowcol, const NRSMat<double> *metric) //modified Gram-Schmidt
+{
+if(metric) //general metric
+{
+if(rowcol) //vectors are rows
+	{
+	if((*metric).nrows() != mm) laerror("incompatible metric in orthonormalize");
+	        for(int j=0; j<nn; ++j)
+                {
+                for(int i=0; i<j; ++i)
+                        {
+			NRVec<double> tmp = *metric * (*this).row(i);
+                        double fact = cblas_ddot(mm,(*this)[j],1,tmp,1);
+                        cblas_daxpy(mm,-fact,(*this)[i],1,(*this)[j],1);
+                        }
+		NRVec<double> tmp = *metric * (*this).row(j);
+                double norm = cblas_ddot(mm,(*this)[j],1,tmp,1);
+                if(norm<=0.) laerror("zero vector in orthonormalize or nonpositive metric");
+                cblas_dscal(mm,1./sqrt(norm),(*this)[j],1);
+                }
+	}
+else	//vectors are columns
+	{
+	if((*metric).nrows() != nn) laerror("incompatible metric in orthonormalize");
+        for(int j=0; j<mm; ++j)
+                {
+                for(int i=0; i<j; ++i)
+                        {
+			NRVec<double> tmp = *metric * (*this).column(i);
+                        double fact = cblas_ddot(nn,&(*this)[0][j],mm,tmp,1);
+                        cblas_daxpy(nn,-fact,&(*this)[0][i],mm,&(*this)[0][j],mm);
+                        }
+		NRVec<double> tmp = *metric * (*this).column(j);
+		double norm = cblas_ddot(nn,&(*this)[0][j],mm,tmp,1);
+                if(norm<=0.) laerror("zero vector in orthonormalize or nonpositive metric");
+                cblas_dscal(nn,1./sqrt(norm),&(*this)[0][j],mm);
+                }
+	}
+}
 
+else //unit metric
+
+{
+if(rowcol) //vectors are rows
+	{
+	for(int j=0; j<nn; ++j)
+		{
+		for(int i=0; i<j; ++i)
+			{
+			double fact = cblas_ddot(mm,(*this)[j],1,(*this)[i],1);
+			cblas_daxpy(mm,-fact,(*this)[i],1,(*this)[j],1);
+			}
+		double norm = cblas_dnrm2(mm,(*this)[j],1);
+		if(norm==0.) laerror("zero vector in orthonormalize");
+		cblas_dscal(mm,1./norm,(*this)[j],1);
+		}
+	}
+else //vectors are columns
+	{
+        for(int j=0; j<mm; ++j)
+                {
+                for(int i=0; i<j; ++i)
+                        {
+                        double fact = cblas_ddot(nn,&(*this)[0][j],mm,&(*this)[0][i],mm);
+                        cblas_daxpy(nn,-fact,&(*this)[0][i],mm,&(*this)[0][j],mm);
+                        }
+                double norm = cblas_dnrm2(nn,&(*this)[0][j],mm);
+                if(norm==0.) laerror("zero vector in orthonormalize");
+                cblas_dscal(nn,1./norm,&(*this)[0][j],mm);
+                }
+	}
+}
+}
 
 
 
@@ -1177,10 +1250,14 @@ for (int i=0; i< nn; i++) v[i][i] = r[i];
 //// forced instantization in the corresponding object file
 template class NRMat<double>;
 template class NRMat<complex<double> >;
+template class NRMat<long long>;
+template class NRMat<long>;
 template class NRMat<int>;
 template class NRMat<short>;
 template class NRMat<char>;
 template class NRMat<unsigned char>;
+template class NRMat<unsigned short>;
 template class NRMat<unsigned int>;
 template class NRMat<unsigned long>;
+template class NRMat<unsigned long long>;