[Toon-members] TooN LU.h

[Tom Drummond]

CVSROOT:        /sources/toon
Module name:    TooN
Changes by:     Tom Drummond <twd20>    06/07/09 19:53:01

Modified files:
        .              : LU.h 

Log message:
        added determinant()

CVSWeb URLs:
http://cvs.savannah.gnu.org/viewcvs/TooN/LU.h?cvsroot=toon&r1=1.10&r2=1.11

Patches:
Index: LU.h
===================================================================
RCS file: /sources/toon/TooN/LU.h,v
retrieving revision 1.10
retrieving revision 1.11
diff -u -b -r1.10 -r1.11
--- LU.h        17 Jun 2006 14:31:19 -0000      1.10
+++ LU.h        9 Jul 2006 19:53:01 -0000       1.11
@@ -45,9 +45,10 @@
     int lda = Size;
     int M = Size;
     int N = Size;
-    dgetrf_(&M,&N,my_lu.get_data_ptr(),&lda,my_IPIV,&my_info);
-    if(my_info < 0){
-      std::cerr << "error in LU, INFO was " << my_info << std::endl;
+    int info;
+    dgetrf_(&M,&N,my_lu.get_data_ptr(),&lda,my_IPIV,&info);
+    if(info != 0){
+      std::cerr << "error in LU, INFO was " << info << std::endl;
     }
   }
 
@@ -145,10 +146,11 @@
     int lda=Size;
     int lwork=-1;
     double size;
-    dgetri_(&N, Inverse.get_data_ptr(), &lda, my_IPIV, &size, &lwork, 
&my_info);
+    int info;
+    dgetri_(&N, Inverse.get_data_ptr(), &lda, my_IPIV, &size, &lwork, &info);
     lwork=int(size);
     double* WORK = new double[lwork];
-    dgetri_(&N, Inverse.get_data_ptr(), &lda, my_IPIV, WORK, &lwork, &my_info);
+    dgetri_(&N, Inverse.get_data_ptr(), &lda, my_IPIV, WORK, &lwork, &info);
     delete [] WORK;
     return Inverse;
   }
@@ -156,10 +158,26 @@
   Matrix<Size,Size,RowMajor>& get_lu(){return my_lu;}
   const Matrix<Size,Size,RowMajor>& get_lu()const {return my_lu;}
 
-  int get_info() const { return my_info; }
+  inline int get_sign() const {
+    int result=1;
+    for(int i=0; i<Size-1; i++){
+      if(my_IPIV[i] > i+1){
+       result=-result;
+      }
+    }
+    return result;
+  }
+
+  inline double determinant() const {
+    double result = get_sign();
+    for (int i=0; i<Size; i++){
+      result*=my_lu(i,i);
+    }
+    return result;
+  }
+
  private:
   Matrix<Size,Size,RowMajor> my_lu;
-  int my_info;
   int my_IPIV[Size];
 };
   
@@ -190,9 +208,10 @@
     int lda = my_lu.num_cols();
     int M = lda;
     int N = lda;
-    dgetrf_(&M,&N,my_lu.get_data_ptr(),&lda,my_IPIV,&my_info);
-    if(my_info < 0){
-      std::cerr << "error in LU, INFO was " << my_info << std::endl;
+    int info;
+    dgetrf_(&M,&N,my_lu.get_data_ptr(),&lda,my_IPIV,&info);
+    if(info != 0){
+      std::cerr << "error in LU, INFO was " << info << std::endl;
     }
   }
 
@@ -293,10 +312,11 @@
     int lda=my_lu.num_rows();;
     int lwork=-1;
     double size;
-    dgetri_(&N, Inverse.get_data_ptr(), &lda, my_IPIV, &size, &lwork, 
&my_info);
+    int info;
+    dgetri_(&N, Inverse.get_data_ptr(), &lda, my_IPIV, &size, &lwork, &info);
     lwork=int(size);
     double* WORK = new double[lwork];
-    dgetri_(&N, Inverse.get_data_ptr(), &lda, my_IPIV, WORK, &lwork, &my_info);
+    dgetri_(&N, Inverse.get_data_ptr(), &lda, my_IPIV, WORK, &lwork, &info);
     delete [] WORK;
     return Inverse;
   }
@@ -304,12 +324,9 @@
   Matrix<-1,-1,RowMajor>& get_lu(){return my_lu;}
   const Matrix<-1,-1,RowMajor>& get_lu()const {return my_lu;}
 
-  int get_info() const { return my_info; }
-
 
  private:
   Matrix<-1,-1,RowMajor> my_lu;
-  int my_info;
   int* my_IPIV;
 };