toon-members
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[Toon-members] TooN Makefile.in lapack.h regressions/regressio...


From: Edward Rosten
Subject: [Toon-members] TooN Makefile.in lapack.h regressions/regressio...
Date: Tue, 06 Apr 2010 16:58:55 +0000

CVSROOT:        /cvsroot/toon
Module name:    TooN
Changes by:     Edward Rosten <edrosten>        10/04/06 16:58:55

Modified files:
        .              : Makefile.in lapack.h 
        regressions    : regression.h 
Added files:
        .              : QR_Lapack.h 
        regressions    : qr.cc qr.txt 

Log message:
        Added preliminary LAPACK-based QR decomposition.
        
        Only woks for cols >= rows.

CVSWeb URLs:
http://cvs.savannah.gnu.org/viewcvs/TooN/Makefile.in?cvsroot=toon&r1=1.25&r2=1.26
http://cvs.savannah.gnu.org/viewcvs/TooN/lapack.h?cvsroot=toon&r1=1.15&r2=1.16
http://cvs.savannah.gnu.org/viewcvs/TooN/QR_Lapack.h?cvsroot=toon&rev=1.1
http://cvs.savannah.gnu.org/viewcvs/TooN/regressions/regression.h?cvsroot=toon&r1=1.2&r2=1.3
http://cvs.savannah.gnu.org/viewcvs/TooN/regressions/qr.cc?cvsroot=toon&rev=1.1
http://cvs.savannah.gnu.org/viewcvs/TooN/regressions/qr.txt?cvsroot=toon&rev=1.1

Patches:
Index: Makefile.in
===================================================================
RCS file: /cvsroot/toon/TooN/Makefile.in,v
retrieving revision 1.25
retrieving revision 1.26
diff -u -b -r1.25 -r1.26
--- Makefile.in 1 Feb 2010 19:35:17 -0000       1.25
+++ Makefile.in 6 Apr 2010 16:58:54 -0000       1.26
@@ -36,7 +36,7 @@
        doxygen 
 
 
-TESTS=lu slice vector_resize gauss_jordan eigen-sqrt determinant chol_toon 
chol_lapack simplex sym_eigen fill so3 complex
+TESTS=lu slice vector_resize gauss_jordan eigen-sqrt determinant chol_toon 
chol_lapack simplex sym_eigen fill so3 complex qr
 
 
 TEST_RESULT=$(TESTS:%=regressions/%.result)

Index: lapack.h
===================================================================
RCS file: /cvsroot/toon/TooN/lapack.h,v
retrieving revision 1.15
retrieving revision 1.16
diff -u -b -r1.15 -r1.16
--- lapack.h    1 Feb 2010 19:41:40 -0000       1.15
+++ lapack.h    6 Apr 2010 16:58:55 -0000       1.16
@@ -1,6 +1,6 @@
 // -*- c++ -*-
 
-// Copyright (C) 2005,2009 Tom Drummond (address@hidden)
+// Copyright (C) 2005,2009,2010 Tom Drummond (address@hidden), E. Rosten
 //
 // This file is part of the TooN Library.  This library is free
 // software; you can redistribute it and/or modify it under the
@@ -75,6 +75,14 @@
                // Cholesky inverse given decomposition
                void dpotri_(const char* UPLO, const int* N, double* A, const 
int* LDA, int* INFO);
                void spotri_(const char* UPLO, const int* N, float* A, const 
int* LDA, int* INFO);
+               
+               // Computes a QR factorization of a general rectangular matrix.
+               void sgeqrf_(int *m, int *n, float *a, int *lda, float *tau, 
float *work, int *lwork, int *info);
+               void dgeqrf_(int *m, int *n, double *a, int *lda, double *tau, 
double *work, int *lwork, int *info);
+               
+               //Reconstruct Q from a QR decomposition
+               void sorgqr_(int* M,int* N,int* K, float* A, int* LDA, float* 
TAU, float* WORK, int* LWORK, int* INFO );
+               void dorgqr_(int* M,int* N,int* K, double* A, int* LDA, double* 
TAU, double* WORK, int* LWORK, int* INFO );
        }
 
 
@@ -152,6 +160,28 @@
                ssyev_(JOBZ, UPLO, N, A, lda, W, WORK, LWORK, INFO);
        }
 
+       //QR decomposition
+       void geqrf_(int *m, int *n, float *a, int *lda, float *tau, float 
*work, int *lwork, int *info)
+       {
+               sgeqrf_(m, n, a, lda, tau, work, lwork, info);
+       }
+
+       void geqrf_(int *m, int *n, double *a, int *lda, double *tau, double 
*work, int *lwork, int *info)
+       {
+               dgeqrf_(m, n, a, lda, tau, work, lwork, info);
+       }
+       
+       void orgqr_(int* M,int* N,int* K, float* A, int* LDA, float* TAU, 
float* WORK, int* LWORK, int* INFO )
+       {
+               sorgqr_(M, N, K, A, LDA, TAU, WORK, LWORK, INFO);
+       }
+
+       void orgqr_(int* M,int* N,int* K, double* A, int* LDA, double* TAU, 
double* WORK, int* LWORK, int* INFO )
+       {
+               dorgqr_(M, N, K, A, LDA, TAU, WORK, LWORK, INFO);
+       }
+
+       //Non symmetric (general) eigen decomposition
        inline void geev_(const char* JOBVL, const char* JOBVR, int* N, double* 
A, int* lda, double* WR, double* WI, double* VL, int* LDVL, double* VR, int* 
LDVR , double* WORK, int* LWORK, int* INFO){
                dgeev_(JOBVL, JOBVR, N,  A,  lda,  WR,  WI,  VL,  LDVL,  VR,  
LDVR ,  WORK,  LWORK,  INFO);
        }

Index: regressions/regression.h
===================================================================
RCS file: /cvsroot/toon/TooN/regressions/regression.h,v
retrieving revision 1.2
retrieving revision 1.3
diff -u -b -r1.2 -r1.3
--- regressions/regression.h    15 Oct 2009 14:51:50 -0000      1.2
+++ regressions/regression.h    6 Apr 2010 16:58:55 -0000       1.3
@@ -1,6 +1,7 @@
 #include <TooN/TooN.h>
 #include <TooN/LU.h>
 #include <TooN/SVD.h>
+#include <TooN/QR_Lapack.h>
 #include <TooN/helpers.h>
 #include <TooN/determinant.h>
 #include <iomanip>

Index: QR_Lapack.h
===================================================================
RCS file: QR_Lapack.h
diff -N QR_Lapack.h
--- /dev/null   1 Jan 1970 00:00:00 -0000
+++ QR_Lapack.h 6 Apr 2010 16:58:55 -0000       1.1
@@ -0,0 +1,114 @@
+#ifndef TOON_INCLUDE_QR_LAPACK_H
+#define TOON_INCLUDE_QR_LAPACK_H
+
+
+#include <TooN/TooN.h>
+#include <TooN/lapack.h>
+#include <utility>
+
+namespace TooN{
+
+/**
+Performs %QR decomposition.
+
+***WARNING*** this will only work if the number of columns is greater than 
+the number of rows!
+
address@hidden gDecomps
+*/
+template<int Rows=Dynamic, int Cols=Rows, class Precision=double>
+class QR_Lapack{
+
+       private:
+               static const int square_Size = (Rows>=0 && 
Cols>=0)?(Rows<Cols?Rows:Cols):Dynamic;
+
+       public: 
+               /// Construct the %QR decomposition of a matrix. This 
initialises the class, and
+               /// performs the decomposition immediately.
+               template<int R, int C, class P, class B> 
+               QR_Lapack(const Matrix<R,C,P,B>& m)
+               :copy(m),tau(square_size()), Q(square_size(), square_size())
+               {
+                       compute();
+               }
+               
+               ///Return L
+               const Matrix<Rows, Cols, Precision, ColMajor>& get_R()
+               {
+                       return copy;
+               }
+               
+               ///Return Q
+               const Matrix<square_Size, square_Size, Precision, ColMajor>& 
get_Q()
+               {
+                       return Q;
+               }       
+
+
+       private:
+
+               void compute()
+               {       
+                       int M = copy.num_rows();
+                       int N = copy.num_cols();
+                       
+                       int LWORK=-1;
+                       int INFO;
+                       int lda = M;
+
+                       Precision size;
+                       
+                       //Compute the working space
+                       geqrf_(&M, &N, copy.get_data_ptr(), &lda, 
tau.get_data_ptr(), &size, &LWORK, &INFO);
+
+                       LWORK = (int) size;
+
+                       Precision* work = new Precision[LWORK];
+
+                       geqrf_(&M, &N, copy.get_data_ptr(), &lda, 
tau.get_data_ptr(), work, &LWORK, &INFO);
+
+
+                       if(INFO < 0)
+                               std::cerr << "error in QR, INFO was " << INFO 
<< std::endl;
+
+                       //The upper "triangle+" of copy is R
+                       //The lower right and tau contain enough information to 
reconstruct Q
+                       
+                       //LAPACK provides a handy function to do the 
reconstruction
+                       Q = copy.template slice<0,0,square_Size, 
square_Size>(0,0,square_size(), square_size());
+                       
+                       int K = square_size();
+                       M=K;
+                       N=K;
+                       lda = K;
+                       orgqr_(&M, &N, &K, Q.get_data_ptr(), &lda, 
tau.get_data_ptr(), work, &LWORK, &INFO);
+
+                       if(INFO < 0)
+                               std::cerr << "error in QR, INFO was " << INFO 
<< std::endl;
+
+                       delete [] work;
+                       
+                       //Now zero out the lower triangle
+                       for(int r=1; r < square_size(); r++)
+                               for(int c=0; c<r; c++)
+                                       copy[r][c] = 0;
+
+               }
+
+               Matrix<Rows, Cols, Precision, ColMajor> copy;
+               Matrix<square_Size, square_Size, Precision, ColMajor> Q;
+               Vector<square_Size, Precision> tau;
+
+
+               int square_size()
+               {
+                       return std::min(copy.num_rows(), copy.num_cols());      
+               }
+
+
+};
+
+}
+
+
+#endif

Index: regressions/qr.cc
===================================================================
RCS file: regressions/qr.cc
diff -N regressions/qr.cc
--- /dev/null   1 Jan 1970 00:00:00 -0000
+++ regressions/qr.cc   6 Apr 2010 16:58:55 -0000       1.1
@@ -0,0 +1,33 @@
+#include "regressions/regression.h"
+using namespace TooN;
+using namespace std;
+
+int main()
+{
+       Matrix<3,4> m;
+       
+       m = Data(5.4388593399963903e-01,
+9.9370462412085203e-01,
+1.0969746452319418e-01,
+4.4837291206649532e-01,
+7.2104662057981139e-01,
+2.1867663239963386e-01,
+6.3591370975105699e-02,
+3.6581617683817125e-01,
+5.2249530577710213e-01,
+1.0579827325022817e-01,
+4.0457999585762583e-01,
+7.6350464084881342e-01);
+       
+       cout << setprecision(20) << scientific;
+       cout << m << endl;
+
+       QR_Lapack<3, 4> q(m);
+
+       cout << q.get_R() << endl;
+       cout << q.get_Q() << endl;
+
+       cout << q.get_Q() * q.get_R() - m << endl;
+
+}
+

Index: regressions/qr.txt
===================================================================
RCS file: regressions/qr.txt
diff -N regressions/qr.txt
--- /dev/null   1 Jan 1970 00:00:00 -0000
+++ regressions/qr.txt  6 Apr 2010 16:58:55 -0000       1.1
@@ -0,0 +1,23 @@
+#The matrix
+5.4388593399963903e-01 9.9370462412085203e-01 1.0969746452319418e-01 
4.4837291206649532e-01
+7.2104662057981139e-01 2.1867663239963386e-01 6.3591370975105699e-02 
3.6581617683817125e-01
+5.2249530577710213e-01 1.0579827325022817e-01 4.0457999585762583e-01 
7.6350464084881342e-01
+
+#QR decomp computed using MATLAB
+
+#> t 1e-10
+
+#R
+-1.0434181725517979e+00 -7.2206631564721513e-01 -3.0371945598870015e-01 
-8.6883845111442815e-01
+0.0000000000000000e+00 7.2462532386550971e-01 -7.3953941785256813e-02 
-2.9029928014361439e-02
+0.0000000000000000e+00 0.0000000000000000e+00 2.8643965469504806e-01 
4.0258674748893997e-01
+
+#Q
+-5.2125403630790146e-01 8.5192253468958445e-01 5.0221753461968818e-02
+-6.9104280483864855e-01 -3.8682349403614080e-01 -6.1059595997877780e-01
+-5.0075350374555994e-01 -3.5298099006850986e-01 7.9034824548220517e-01
+
+#Difference
+0 0 0 0 
+0 0 0 0 
+0 0 0 0 




reply via email to

[Prev in Thread] Current Thread [Next in Thread]