[Top][All Lists]
[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
[Prev in Thread] |
Current Thread |
[Next in Thread] |
- [Toon-members] TooN Makefile.in lapack.h regressions/regressio...,
Edward Rosten <=