conversion constructor from vec3 and mat3 to nrvec and nrmat

la_random.h

@@ -3,23 +3,23 @@
 
 namespace LA {
 
-//RANDOM numbers defaulting to standard library but switchable to user's functions
-#ifndef RANDDOUBLE
-#define RANDDOUBLE randdouble
-#endif
-
-#ifndef RANDDOUBLESIGNED
-#define RANDDOUBLESIGNED randdoublesigned
-#endif
-
-#ifndef RANDINT32
-#define RANDINT32 randint32
-#endif
-
 extern double randdouble();
 extern double randdoublesigned();
 extern int randint32();
 
+//RANDOM numbers defaulting to standard library but switchable to user's functions
+#ifndef RANDDOUBLE
+#define RANDDOUBLE LA::randdouble
+#endif
+
+#ifndef RANDDOUBLESIGNED
+#define RANDDOUBLESIGNED LA::randdoublesigned
+#endif
+
+#ifndef RANDINT32
+#define RANDINT32 LA::randint32
+#endif
+
 #ifdef __GNUC__
 #define WEAK_SYMBOL __attribute__((weak)) 
 #else

la_traits.h

@@ -68,6 +68,13 @@ extern "C" {
 }
 #endif
 
+
+//forward declarations
+namespace LA_Vecmat3 {
+template<typename C> class Vec3;
+template<typename C> class Mat3;
+}
+
 namespace LA {
 
 

mat.cc

@@ -3411,6 +3411,8 @@ for(int j=0; j<nn; ++j) (*this)(j,i) *= f;
  ******************************************************************************/
 template class NRMat<double>;
 template class NRMat<std::complex<double> >;
+//template class NRMat<float>;
+//template class NRMat<std::complex<float> >;
 template class NRMat<long long>;
 template class NRMat<long>;
 template class NRMat<int>;

mat.h

@@ -22,6 +22,10 @@
 #define _LA_MAT_H_
 #include "la_traits.h"
 
+using namespace LA_Vecmat3;
+
+#include "vecmat3.h"
+
 
 namespace LA {
 
@@ -82,8 +86,9 @@ public:
 	inline NRMat(const T &a, const int n, const int m);
 
 	//! inlined constructor creating matrix of given size filled with data located at given memory location 
-	NRMat(const T *a, const int n, const int m);
+	inline NRMat(const T *a, const int n, const int m);
-
+	inline NRMat(const Mat3<T> &rhs) : NRMat(&rhs(0,0),3,3) {};
+	
 	//! inlined constructor creating matrix of given size from an array
 	template<int R, int C> inline NRMat(const T (&a)[R][C]);
 
@@ -615,7 +620,6 @@ NRMat<T>::NRMat(const T *a, const int n, const int m) : nn(n), mm(m), count(new
 		cublasSetVector(nm, sizeof(T), a, 1, v, 1);
 	}
 #endif
-
 }
 
 /***************************************************************************//**

smat.cc

@@ -872,7 +872,8 @@ NRSMat<std::complex<double> >  NRSMat<std::complex<double> >::inverse() {return
  ******************************************************************************/
 template class NRSMat<double>;
 template class NRSMat<std::complex<double> >;
-
+//template class NRSMat<float>;
+//template class NRSMat<std::complex<float> >;
 template class NRSMat<long long>;
 template class NRSMat<long>;
 template class NRSMat<int>;

t.cc

@@ -426,7 +426,7 @@ cout <<b;
 cout <<d;
 }
 
-if(1)
+if(0)
 {
 NRMat<double> a;
 cin >>a ;
@@ -2784,5 +2784,32 @@ NRSMat<char> adjperm = adj.permuted(p);
 cout <<"resorted graph = "<<adjperm<<endl;
 }
 
+if(1)
+{
+Mat3<double> a;
+a.randomize(10.);
+cout<<a<<endl;
+Mat3<double> ata = a.transpose()*a;
+Mat3<double> work(ata);
+Vec3<double> w;
+Mat3<double> v;
+work.eivec_sym(w,v,true);
+cout <<w<<endl<<sqrt(w[0])<<" "<<sqrt(w[1])<<" "<<sqrt(w[2])<<" "<<endl<<endl<<v<<endl;
+
+Mat3<double> aat = a*a.transpose();
+Mat3<double> work2(aat);
+Vec3<double> w2; 
+Mat3<double> v2;
+work2.eivec_sym(w2,v2,true);
+cout <<w2<<endl<<sqrt(w2[0])<<" "<<sqrt(w2[1])<<" "<<sqrt(w2[2])<<" "<<endl<<endl<<v2<<endl;
+
+
+NRMat<double> aa(a);
+cout <<aa.nrows()<<" "<<aa.ncols()<<endl;
+NRMat<double> uu(3,3),vv(3,3);
+NRVec<double> ss(3);
+singular_decomposition(aa,&uu,ss,&vv,0);
+cout <<uu<<ss<<vv<<endl;
+}
 
 }

vec.cc

@@ -27,7 +27,7 @@
 #include <errno.h>
 #include "vec.h"
 #include <unistd.h>
-
+#include "vecmat3.h"
 
 namespace LA {
 
@@ -904,6 +904,7 @@ void NRVec<T>::storesubvector(const NRVec<int> &selection, const NRVec &rhs)
 }
 
 
+
 /***************************************************************************//**
  * forced instantization in the corespoding object file
  ******************************************************************************/
@@ -1039,6 +1040,8 @@ INSTANTIZE_CONCAT(std::complex<double>)
 
 
 
+//template class NRVec<float>;
+//template class NRVec<std::complex<float> >;
 template class NRVec<double>;
 template class NRVec<std::complex<double> >;
 template class NRVec<char>;

vec.h

@@ -21,8 +21,11 @@
 #define _LA_VEC_H_
 
 #include "la_traits.h"
+#include "vecmat3.h"
 #include <list>
 
+using namespace LA_Vecmat3;
+
 namespace LA {
 
 /***************************************************************************//**
@@ -134,6 +137,7 @@ public:
 	
 	//! inlined constructor creating vector of given size filled with data located at given memory location 
         inline NRVec(const T *a, const int n);
+	inline NRVec(const Vec3<T> &rhs) : NRVec(&rhs[0],3) {};
 	
 	//! inlined constructor creating vector of given size filled with data located at given memory location 
 	inline NRVec(T *a, const int n, bool skeleton);

vecmat3.cc

@@ -17,6 +17,9 @@
 */
 
 #include "vecmat3.h"
+#ifndef QUAT_NO_RANDOM 
+#include "la_random.h"
+#endif
 
 namespace LA_Vecmat3 {
 
@@ -613,6 +616,24 @@ tmp=(q[2][1]+q[1][2])/2;
 q[2][1]=q[1][2]=tmp;
 }
 
+
+#ifndef QUAT_NO_RANDOM
+template <>
+void Vec3<double>::randomize(const double x)
+{
+for(int i=0;i<3;++i)  q[i]=x*RANDDOUBLESIGNED();
+}
+
+
+template <>
+void Mat3<double>::randomize(const double x, const bool symmetric)
+{
+if(symmetric)  for(int i=0;i<3;++i) for(int j=0; j<=i; ++j) q[j][i]=q[i][j]=x*RANDDOUBLESIGNED();
+else for(int i=0;i<3;++i) for(int j=0; j<3; ++j) q[i][j]=x*RANDDOUBLESIGNED();
+}
+
+#endif
+
 //eigensolver for 3x3 matrix by Joachim Kopp - analytic formula version,
 //might be unstable for ill-conditioned ones, then use other methods
 //cf. arxiv physics 0610206v3
@@ -651,7 +672,7 @@ q[2][1]=q[1][2]=tmp;
 #define SQR(x)      ((x)*(x))                        // x^2
 //
 template <typename T>
-void Mat3<T>::eival_sym(Vec3<T> &w) const
+void Mat3<T>::eival_sym(Vec3<T> &w, const bool sortdown) const
 {
 T  m, c1, c0;
   
@@ -685,10 +706,19 @@ T  m, c1, c0;
   w[0]  = w[1] + c;
   w[1] -= s;
 
+if(sortdown)
+	{
+        if(w[0]<w[1]) {T tmp=w[0]; w[0]=w[1]; w[1]=tmp;}
+        if(w[0]<w[2]) {T tmp=w[0]; w[0]=w[2]; w[2]=tmp;}
+        if(w[1]<w[2]) {T tmp=w[1]; w[1]=w[2]; w[2]=tmp;}
+	}
+else
 //sort in ascending order
-if(w[0]>w[1]) {T tmp=w[0]; w[0]=w[1]; w[1]=tmp;}
+	{
-if(w[0]>w[2]) {T tmp=w[0]; w[0]=w[2]; w[2]=tmp;}
+	if(w[0]>w[1]) {T tmp=w[0]; w[0]=w[1]; w[1]=tmp;}
-if(w[1]>w[2]) {T tmp=w[1]; w[1]=w[2]; w[2]=tmp;}
+	if(w[0]>w[2]) {T tmp=w[0]; w[0]=w[2]; w[2]=tmp;}
+	if(w[1]>w[2]) {T tmp=w[1]; w[1]=w[2]; w[2]=tmp;}
+	}
 }
 
 // Calculates the eigenvalues and normalized eigenvectors of a symmetric 3x3
@@ -717,7 +747,7 @@ if(w[1]>w[2]) {T tmp=w[1]; w[1]=w[2]; w[2]=tmp;}
 //   v1.0: First released version
 // ----------------------------------------------------------------------------
 template <typename T>
-void Mat3<T>::eivec_sym(Vec3<T> &w, Mat3 &v) const
+void Mat3<T>::eivec_sym(Vec3<T> &w, Mat3 &v, const bool sortdown) const
 {
   T norm;          // Squared norm or inverse norm of current eigenvector
   T n0, n1;        // Norm of first and second columns of A
@@ -729,7 +759,7 @@ void Mat3<T>::eivec_sym(Vec3<T> &w, Mat3 &v) const
   int i, j;             // Loop counters
 
   // Calculate eigenvalues
-  eival_sym(w);
+  eival_sym(w,sortdown);
 Mat3<T> A(*this); //scratch copy
 
   wmax = fabs(w[0]);

vecmat3.h

@@ -16,7 +16,7 @@
     along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */
 
-//this header defines simple classes for 3-dimensional vectors and matrices to describe rotations etc.
+//this header defines simple classes for 3-dimensional REAL-valued vectors and matrices to describe rotations etc.
 //the class is compatible with functions in quaternion.h used for SO(3) parametrization
 //it should be compilable separately from LA as well as being a part of LA
 
@@ -62,7 +62,7 @@ public:
 	//compiler generates default copy constructor and assignment operator
 	
 	//formal indexing
-	inline const T operator[](const int i) const {return q[i];};
+	inline const T& operator[](const int i) const {return q[i];};
 	inline T& operator[](const int i) {return q[i];};
 
 	//operations of Vec3s with scalars
@@ -89,6 +89,7 @@ public:
 	const Vec3 timesT(const Mat3<T> &rhs) const; //with transpose
 	Mat3<T> outer(const Vec3 &rhs) const; //tensor product
 	void inertia(Mat3<T> &itensor, const T weight) const; //contribution to inertia tensor
+	void randomize(const T x);
 
 	//C-style IO
 	int fprintf(FILE *f, const char *format) const {return ::fprintf(f,format,q[0],q[1],q[2]);};
@@ -130,7 +131,7 @@ public:
 	//formal indexing
 	inline const T* operator[](const int i) const {return q[i];};
 	inline T* operator[](const int i) {return q[i];};
-	inline const T operator()(const int i, const int j) const {return q[i][j];};
+	inline const T& operator()(const int i, const int j) const {return q[i][j];};
 	inline T& operator()(const int i, const int j) {return q[i][j];};
 	
 
@@ -145,6 +146,8 @@ public:
         const Mat3 operator*(const T rhs) const {return Mat3(*this) *= rhs;};
         const Mat3 operator/(const T rhs) const {return Mat3(*this) /= rhs;};
 
+	void randomize(const T x, const bool symmetric=false);
+
 	//Mat3 algebra
 	const Mat3 operator-() const {return *this * (T)-1;}; //unary minus
 	Mat3& operator+=(const Mat3 &rhs);
@@ -166,8 +169,8 @@ public:
         int fprintf(FILE *f, const char *format) const {int n= ::fprintf(f,format,q[0][0],q[0][1],q[0][2]); n+=::fprintf(f,format,q[1][0],q[1][1],q[1][2]); n+=::fprintf(f,format,q[2][0],q[2][1],q[2][2]); return n;};
         int fscanf(FILE *f, const char *format) const {return ::fscanf(f,format,q[0][0],q[0][1],q[0][2]) + ::fscanf(f,format,q[1][0],q[1][1],q[1][2]) + ::fscanf(f,format,q[2][0],q[2][1],q[2][2]);};
 	void symmetrize(); //average offdiagonal elements
-	void eival_sym(Vec3<T> &w) const; //only for real symmetric matrix, symmetry is not checked
+	void eival_sym(Vec3<T> &w, const bool sortdown=false) const; //only for real symmetric matrix, symmetry is not checked
-	void eivec_sym(Vec3<T> &w, Mat3 &v) const; //only for real symmetric matrix, symmetry is not checked
+	void eivec_sym(Vec3<T> &w, Mat3 &v, const bool sortdown=false) const; //only for real symmetric matrix, symmetry is not checked
 	T norm(const T scalar = 0) const;
 };