[Toon-members] TooN LU.h lapack.h test/lutest.cc

[Edward Rosten]

CVSROOT:        /cvsroot/toon
Module name:    TooN
Changes by:     Edward Rosten <edrosten>        09/02/26 21:50:27

Modified files:
        .              : LU.h lapack.h 
Added files:
        test           : lutest.cc 

Log message:
        Fixed LU to work in new TooN. Works for float and double. Only a single
        implementation is now required, covering static and dynamic, and only 
        a single instance of backsub is needed for all types of matrix, to its
        about 40% of the length of the old version.
        
        WARNING WARNING WARNING
        
        There's a bug somewhere so Matrix*Matrix doesn't work for some mixes
        of static and dynamic.

CVSWeb URLs:
http://cvs.savannah.gnu.org/viewcvs/TooN/LU.h?cvsroot=toon&r1=1.12&r2=1.13
http://cvs.savannah.gnu.org/viewcvs/TooN/lapack.h?cvsroot=toon&r1=1.7&r2=1.8
http://cvs.savannah.gnu.org/viewcvs/TooN/test/lutest.cc?cvsroot=toon&rev=1.1

Patches:
Index: LU.h
===================================================================
RCS file: /cvsroot/toon/TooN/LU.h,v
retrieving revision 1.12
retrieving revision 1.13
diff -u -b -r1.12 -r1.13
--- LU.h        10 Jul 2006 10:48:54 -0000      1.12
+++ LU.h        26 Feb 2009 21:50:13 -0000      1.13
@@ -1,5 +1,5 @@
 /*                       
-        Copyright (C) 2005 Tom Drummond
+        Copyright (C) 2005 Tom Drummond, 2009 Edward Rosten
 
      This library is free software; you can redistribute it and/or
      modify it under the terms of the GNU Lesser General Public
@@ -16,8 +16,8 @@
      Foundation, Inc.
      51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
 */
-#ifndef __LU_H
-#define __LU_H
+#ifndef TOON_INCLUDE_LU_H
+#define TOON_INCLUDE_LU_H
 
 #include <iostream>
 
@@ -25,47 +25,57 @@
 
 #include <TooN/TooN.h>
 
-#ifndef TOON_NO_NAMESPACE
 namespace TooN {
-#endif 
 
-template <int Size=-1>
+template <int Size=-1, class Precision=double>
 class LU {
   public:
 
-  template<class Accessor>
-  LU(const FixedMatrix<Size,Size,Accessor>& m){
+  template<int S1, int S2, class Base>
+  LU(const Matrix<S1,S2,Precision, Base>& m)
+  :my_lu(m.num_rows(),m.num_cols()),my_IPIV(m.num_rows()){
     compute(m);
   }
 
-  template<class Accessor>
-  void compute(const FixedMatrix<Size,Size,Accessor>& m){
+  template<int S1, int S2, class Base>
+  void compute(const Matrix<S1,S2,Precision,Base>& m){
+    //check for consistency with Size
+    SizeMismatch<Size, S1>::test(my_lu.num_rows(),m.num_rows());
+    SizeMismatch<Size, S2>::test(my_lu.num_rows(),m.num_cols());
+       
+    //Make a local copy. This is guaranteed contiguous
     my_lu=m;
-    int lda = Size;
-    int M = Size;
-    int N = Size;
-    dgetrf_(&M,&N,my_lu.get_data_ptr(),&lda,my_IPIV,&my_info);
+    int lda = m.num_rows();
+    int M = m.num_rows();
+    int N = m.num_rows();
+
+    getrf_(&M,&N,&my_lu[0][0],&lda,&my_IPIV[0],&my_info);
+
     if(my_info < 0){
       std::cerr << "error in LU, INFO was " << my_info << std::endl;
     }
   }
 
-  template <int NRHS, class Accessor>
-  Matrix<Size,NRHS,RowMajor> backsub(const FixedMatrix<Size,NRHS,Accessor>& 
rhs){
-    Matrix<Size,NRHS,RowMajor> result(rhs);
-    int M=NRHS;
-    int N=Size;
+  template <int Rows, int NRHS, class Base>
+  Matrix<Size,NRHS> backsub(const Matrix<Rows,NRHS,Precision,Base>& rhs){
+    //Check the number of rows is OK.
+    SizeMismatch<Size, Rows>::test(my_lu.num_rows(), rhs.num_rows());
+       
+    Matrix<Size, NRHS> result(rhs);
+
+    int M=rhs.num_cols();
+    int N=my_lu.num_rows();
     double alpha=1;
-    int lda=Size;
-    int ldb=NRHS;
-    
dtrsm_("R","U","N","N",&M,&N,&alpha,my_lu.get_data_ptr(),&lda,result.get_data_ptr(),&ldb);
-    
dtrsm_("R","L","N","U",&M,&N,&alpha,my_lu.get_data_ptr(),&lda,result.get_data_ptr(),&ldb);
+    int lda=my_lu.num_rows();
+    int ldb=rhs.num_cols();
+    trsm_("R","U","N","N",&M,&N,&alpha,&my_lu[0][0],&lda,&result[0][0],&ldb);
+    trsm_("R","L","N","U",&M,&N,&alpha,&my_lu[0][0],&lda,&result[0][0],&ldb);
 
     // now do the row swapping (lapack dlaswp.f only shuffles fortran rows = 
Rowmajor cols)
     for(int i=N-1; i>=0; i--){
       const int swaprow = my_IPIV[i]-1; // fortran arrays start at 1
       for(int j=0; j<NRHS; j++){
-       double temp = result[i][j];
+       Precision temp = result[i][j];
        result[i][j] = result[swaprow][j];
        result[swaprow][j] = temp;
       }
@@ -73,32 +83,9 @@
     return result;
   }
 
-  template <class Accessor>
-  Matrix<-1,-1,RowMajor> backsub(const DynamicMatrix<Accessor>& rhs){
-    Matrix<-1,-1,RowMajor> result(rhs);
-    assert(result.num_rows == my_lu.num_rows());
-    int M=result.num_cols();
-    int N=result.num_rows();
-    double alpha=1;
-    int lda=result.num_rows();
-    int ldb=result.num_cols();
-    
dtrsm_("R","U","N","N",&M,&N,&alpha,my_lu.get_data_ptr(),&lda,result.get_data_ptr(),&ldb);
-    
dtrsm_("R","L","N","U",&M,&N,&alpha,my_lu.get_data_ptr(),&lda,result.get_data_ptr(),&ldb);
 
-    // now do the row swapping (lapack dlaswp.f only shuffles fortran rows = 
Rowmajor cols)
-    for(int i=N-1; i>=0; i--){
-      const int swaprow = my_IPIV[i]-1; // fortran arrays start at 1
-      for(int j=0; j<result.num_cols(); j++){
-       double temp = result(i,j);
-       result(i,j) = result(swaprow,j);
-       result(swaprow,j) = temp;
-      }
-    }
-    return result;
-  }
-
-
-  template <class Accessor>
+  //FIXME fix for Vector
+  /*template <class Accessor>
   Vector<Size> backsub(const FixedVector<Size,Accessor>& rhs){
     Vector<Size> result(rhs);
     int M=1;
@@ -115,211 +102,24 @@
       result[swaprow]=temp;
     }
     return result;
-  }
-
-  template <class Accessor>
-  Vector<> backsub(const DynamicVector<Accessor>& rhs){
-    assert(rhs.size()==Size);
-    Vector<> result(rhs);
-    int M=1;
-    int N=Size;
-    double alpha=1;
-    int lda=Size;
-    int ldb=1;
-    
dtrsm_("R","U","N","N",&M,&N,&alpha,my_lu.get_data_ptr(),&lda,result.get_data_ptr(),&ldb);
-    
dtrsm_("R","L","N","U",&M,&N,&alpha,my_lu.get_data_ptr(),&lda,result.get_data_ptr(),&ldb);
-    for(int i=N-1; i>=0; i--){
-      const int swaprow = my_IPIV[i]-1; // fortran arrays start at 1
-      double temp = result[i];
-      result[i]=result[swaprow];
-      result[swaprow]=temp;
-    }
-    return result;
-  }
-
-
-  Matrix<Size,Size,RowMajor> get_inverse(){
-    Matrix<Size,Size,RowMajor> Inverse = my_lu;
-    int N = Size;
-    int lda=Size;
-    int lwork=-1;
-    double size;
-    dgetri_(&N, Inverse.get_data_ptr(), &lda, my_IPIV, &size, &lwork, 
&my_info);
-    lwork=int(size);
-    double* WORK = new double[lwork];
-    dgetri_(&N, Inverse.get_data_ptr(), &lda, my_IPIV, WORK, &lwork, &my_info);
-    delete [] WORK;
-    return Inverse;
-  }
-
-  Matrix<Size,Size,RowMajor>& get_lu(){return my_lu;}
-  const Matrix<Size,Size,RowMajor>& get_lu()const {return my_lu;}
-
-  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;
-  }
-
-  int get_info() const { return my_info; }
- private:
-  Matrix<Size,Size,RowMajor> my_lu;
-  int my_info;
-  int my_IPIV[Size];
-};
-  
-
-template <>
-class LU<> {
-  public:
-
-  LU(int size):my_lu(size,size){
-    my_IPIV = new int[size];
-  }
-
-  template<class Accessor>
-  LU(const DynamicMatrix<Accessor>& m) :
-    my_lu(m.num_rows(),m.num_cols())
-  {
-    my_IPIV = new int[m.num_rows()];
-    assert(m.num_rows() == m.num_cols());
-    compute(m);
-  }
-
-  ~LU(){delete[] my_IPIV;}
+  }*/
 
-
-  template<class Accessor>
-  void compute(const DynamicMatrix<Accessor>& m){
-    my_lu=m;
-    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;
-    }
-  }
-
-
-  template <int Rows, int Cols, class Accessor>
-  Matrix<Rows,Cols,RowMajor> backsub(const FixedMatrix<Rows,Cols,Accessor>& 
rhs){
-    assert(my_lu.num_rows() == rhs.num_rows());
-    Matrix<Rows,Cols,RowMajor> result(rhs);
-    int M=result.num_cols();
-    int N=result.num_rows();
-    double alpha=1;
-    int lda=my_lu.num_rows();
-    int ldb=result.num_cols();
-    
dtrsm_("R","U","N","N",&M,&N,&alpha,my_lu.get_data_ptr(),&lda,result.get_data_ptr(),&ldb);
-    
dtrsm_("R","L","N","U",&M,&N,&alpha,my_lu.get_data_ptr(),&lda,result.get_data_ptr(),&ldb);
-
-    // now do the row swapping (lapack dlaswp.f only shuffles fortran rows = 
Rowmajor cols)
-    for(int i=N-1; i>=0; i--){
-      const int swaprow = my_IPIV[i]-1; // fortran arrays start at 1
-      for(int j=0; j<result.num_cols(); j++){
-       double temp = result(i,j);
-       result(i,j) = result(swaprow,j);
-       result(swaprow,j) = temp;
-      }
-    }
-    return result;
-  }
-
-
-  template <class Accessor>
-  Matrix<-1,-1,RowMajor> backsub(const DynamicMatrix<Accessor>& rhs){
-    assert(my_lu.num_rows() == rhs.num_rows());
-    Matrix<-1,-1,RowMajor> result(rhs);
-    int M=result.num_cols();
-    int N=result.num_rows();
-    double alpha=1;
-    int lda=my_lu.num_rows();
-    int ldb=result.num_cols();
-    
dtrsm_("R","U","N","N",&M,&N,&alpha,my_lu.get_data_ptr(),&lda,result.get_data_ptr(),&ldb);
-    
dtrsm_("R","L","N","U",&M,&N,&alpha,my_lu.get_data_ptr(),&lda,result.get_data_ptr(),&ldb);
-
-    // now do the row swapping (lapack dlaswp.f only shuffles fortran rows = 
Rowmajor cols)
-    for(int i=N-1; i>=0; i--){
-      const int swaprow = my_IPIV[i]-1; // fortran arrays start at 1
-      for(int j=0; j<result.num_cols(); j++){
-       double temp = result(i,j);
-       result(i,j) = result(swaprow,j);
-       result(swaprow,j) = temp;
-      }
-    }
-    return result;
-  }
-
-  template <int Size, class Accessor>
-  Vector<Size> backsub(const FixedVector<Size,Accessor>& rhs){
-    assert(rhs.size() == my_lu.num_rows());
-    Vector<Size> result(rhs);
-    int M=1;
-    int N=result.size();
-    double alpha=1;
-    int lda=result.size();
-    int ldb=1;
-    
dtrsm_("R","U","N","N",&M,&N,&alpha,my_lu.get_data_ptr(),&lda,result.get_data_ptr(),&ldb);
-    
dtrsm_("R","L","N","U",&M,&N,&alpha,my_lu.get_data_ptr(),&lda,result.get_data_ptr(),&ldb);
-    for(int i=N-1; i>=0; i--){
-      const int swaprow = my_IPIV[i]-1; // fortran arrays start at 1
-      double temp = result[i];
-      result[i]=result[swaprow];
-      result[swaprow]=temp;
-    }
-    return result;
-  }
-
-  
-  template <class Accessor>
-  Vector<> backsub(const DynamicVector<Accessor>& rhs){
-    assert(rhs.size() == my_lu.num_rows());
-    Vector<> result(rhs);
-    int M=1;
-    int N=result.size();
-    double alpha=1;
-    int lda=result.size();
-    int ldb=1;
-    
dtrsm_("R","U","N","N",&M,&N,&alpha,my_lu.get_data_ptr(),&lda,result.get_data_ptr(),&ldb);
-    
dtrsm_("R","L","N","U",&M,&N,&alpha,my_lu.get_data_ptr(),&lda,result.get_data_ptr(),&ldb);
-    for(int i=N-1; i>=0; i--){
-      const int swaprow = my_IPIV[i]-1; // fortran arrays start at 1
-      double temp = result[i];
-      result[i]=result[swaprow];
-      result[swaprow]=temp;
-    }
-    return result;
-  }
-
-  Matrix<-1,-1,RowMajor> get_inverse(){
-    Matrix<-1,-1,RowMajor> Inverse(my_lu);
+  Matrix<Size,Size,Precision> get_inverse(){
+    Matrix<Size,Size,Precision> Inverse(my_lu);
     int N = my_lu.num_rows();
     int lda=my_lu.num_rows();
     int lwork=-1;
-    double size;
-    dgetri_(&N, Inverse.get_data_ptr(), &lda, my_IPIV, &size, &lwork, 
&my_info);
+    Precision size;
+    getri_(&N, &Inverse[0][0], &lda, &my_IPIV[0], &size, &lwork, &my_info);
     lwork=int(size);
-    double* WORK = new double[lwork];
-    dgetri_(&N, Inverse.get_data_ptr(), &lda, my_IPIV, WORK, &lwork, &my_info);
+    Precision* WORK = new Precision[lwork];
+    getri_(&N, &Inverse[0][0], &lda, &my_IPIV[0], WORK, &lwork, &my_info);
     delete [] WORK;
     return Inverse;
   }
 
-  Matrix<-1,-1,RowMajor>& get_lu(){return my_lu;}
-  const Matrix<-1,-1,RowMajor>& get_lu()const {return my_lu;}
+  Matrix<Size,Size,Precision>& get_lu(){return my_lu;}
+  const Matrix<Size,Size,Precision>& get_lu()const {return my_lu;}
 
   inline int get_sign() const {
     int result=1;
@@ -331,8 +131,8 @@
     return result;
   }
 
-  inline double determinant() const {
-    double result = get_sign();
+  inline Precision determinant() const {
+    Precision result = get_sign();
     for (int i=0; i<my_lu.num_rows(); i++){
       result*=my_lu(i,i);
     }
@@ -341,21 +141,14 @@
 
   int get_info() const { return my_info; }
 
-
  private:
-  Matrix<-1,-1,RowMajor> my_lu;
-  int my_info;
-  int* my_IPIV;
-};
-
 
+  Matrix<Size,Size,Precision> my_lu;
+  int my_info;
+  Vector<Size, int> my_IPIV;  //Convenient static-or-dynamic array of ints :-)
 
-#ifndef TOON_NO_NAMESPACE
+};
 }
-#endif 
-
-
-
 
 
 #endif

Index: lapack.h
===================================================================
RCS file: /cvsroot/toon/TooN/lapack.h,v
retrieving revision 1.7
retrieving revision 1.8
diff -u -b -r1.7 -r1.8
--- lapack.h    10 Oct 2007 11:33:05 -0000      1.7
+++ lapack.h    26 Feb 2009 21:50:19 -0000      1.8
@@ -17,24 +17,27 @@
      Foundation, Inc.
      51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
 */
-#ifndef __LAPACK_H
-#define __LAPACK_H
+#ifndef TOON_INCLUDE_LAPCK_H
+#define TOON_INCLUDE_LAPCK_H
 
 // LAPACK and BLAS routines
-
-#ifndef TOON_NO_NAMESPACE
 namespace TooN {
-#endif 
 
 extern "C" {
   // LU decomoposition of a general matrix
   void dgetrf_(int* M, int *N, double* A, int* lda, int* IPIV, int* INFO);
+  void sgetrf_(int* M, int *N, float* A, int* lda, int* IPIV, int* INFO);
+
+
 
   // generate inverse of a matrix given its LU decomposition
   void dgetri_(int* N, double* A, int* lda, int* IPIV, double* WORK, int* 
lwork, int* INFO);
+  void sgetri_(int* N, float* A, int* lda, int* IPIV, float* WORK, int* lwork, 
int* INFO);
 
   // inverse of a triangular matrix * a vector (BLAS level 2)
   void dtrsm_(char* SIDE, char* UPLO, char* TRANSA, char* DIAG, int* M, int* 
N, double* alpha, double* A, int* lda, double* B, int* ldb);
+  void strsm_(char* SIDE, char* UPLO, char* TRANSA, char* DIAG, int* M, int* 
N, float* alpha, float* A, int* lda, float* B, int* ldb);
+  
 
   // SVD of a general matrix of doubles
   void dgesvd_(const char* JOBU, const char* JOBVT, int* M, int *N, double* A, 
int* lda,
@@ -55,11 +58,29 @@
     void dpotri_(const char* UPLO, const int* N, double* A, const int* LDA, 
int* INFO);
 }
 
+void getrf_(int* M, int *N, float* A, int* lda, int* IPIV, int* INFO){
+  sgetrf_(M, N, A, lda, IPIV, INFO);
+}
 
-#ifndef TOON_NO_NAMESPACE
+void getrf_(int* M, int *N, double* A, int* lda, int* IPIV, int* INFO){
+  dgetrf_(M, N, A, lda, IPIV, INFO);
 }
-#endif 
 
+inline void trsm_(char* SIDE, char* UPLO, char* TRANSA, char* DIAG, int* M, 
int* N, float* alpha, float* A, int* lda, float* B, int* ldb)
+{ strsm_(SIDE, UPLO, TRANSA, DIAG, M, N, alpha, A, lda, B, ldb);
+}
+
+inline void trsm_(char* SIDE, char* UPLO, char* TRANSA, char* DIAG, int* M, 
int* N, double* alpha, double* A, int* lda, double* B, int* ldb) {
+  dtrsm_(SIDE, UPLO, TRANSA, DIAG, M, N, alpha, A, lda, B, ldb);
+}
 
+void getri_(int* N, double* A, int* lda, int* IPIV, double* WORK, int* lwork, 
int* INFO){
+       dgetri_(N, A, lda, IPIV, WORK, lwork, INFO);
+}
 
+void getri_(int* N, float* A, int* lda, int* IPIV, float* WORK, int* lwork, 
int* INFO){
+       sgetri_(N, A, lda, IPIV, WORK, lwork, INFO);
+}
+
+}
 #endif

Index: test/lutest.cc
===================================================================
RCS file: test/lutest.cc
diff -N test/lutest.cc
--- /dev/null   1 Jan 1970 00:00:00 -0000
+++ test/lutest.cc      26 Feb 2009 21:50:23 -0000      1.1
@@ -0,0 +1,36 @@
+#include <cstdlib>
+#include <iostream>
+#include <TooN/TooN.h>
+#include <TooN/LU.h>
+
+using namespace TooN;
+using namespace std;
+
+
+int main()
+{
+       Matrix<-1,5,float> m(5,5);
+
+       for(int i=0; i< m.num_rows(); i++)
+               for(int j=0; j< m.num_rows(); j++)
+                       m[i][j] = drand48();
+
+       
+       LU<5,float> mlu(m);
+
+       cout << mlu.get_inverse() << endl;
+       Matrix<5,5,float> inv = mlu.get_inverse();
+       Matrix<5,5,float> a = m*inv;
+
+       for(int i=0; i< m.num_rows(); i++)
+               for(int j=0; j< m.num_rows(); j++)
+               {
+                       if(round(a[i][j]) < 1e-10)
+                               a[i][j] = 0;
+               }
+
+       cout << a;
+
+
+
+}