[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Toon-members] TooN QR_Lapack.h lapack.h
|
From: |
Edward Rosten |
|
Subject: |
[Toon-members] TooN QR_Lapack.h lapack.h |
|
Date: |
Wed, 07 Apr 2010 13:09:59 +0000 |
CVSROOT: /cvsroot/toon
Module name: TooN
Changes by: Edward Rosten <edrosten> 10/04/07 13:09:59
Modified files:
. : QR_Lapack.h lapack.h
Log message:
Optional column pivoting in QR decomposition.
CVSWeb URLs:
http://cvs.savannah.gnu.org/viewcvs/TooN/QR_Lapack.h?cvsroot=toon&r1=1.1&r2=1.2
http://cvs.savannah.gnu.org/viewcvs/TooN/lapack.h?cvsroot=toon&r1=1.17&r2=1.18
Patches:
Index: QR_Lapack.h
===================================================================
RCS file: /cvsroot/toon/TooN/QR_Lapack.h,v
retrieving revision 1.1
retrieving revision 1.2
diff -u -b -r1.1 -r1.2
--- QR_Lapack.h 6 Apr 2010 16:58:55 -0000 1.1
+++ QR_Lapack.h 7 Apr 2010 13:09:59 -0000 1.2
@@ -11,9 +11,20 @@
/**
Performs %QR decomposition.
-***WARNING*** this will only work if the number of columns is greater than
address@hidden this will only work if the number of columns is greater than
the number of rows!
+The QR decomposition operates on a matrix A. It can be performed with
+or without column pivoting. In general:
+\f[
+AP = QR
+\f]
+Where \f$P\f$ is a permutation matrix constructed to permute the columns
+of A. In practise, \f$P\f$ is stored as a vector of integer elements.
+
+With column pivoting, the elements of the leading diagonal of \f$R\f$ will
+be sorted from largest in magnitude to smallest in magnitude.
+
@ingroup gDecomps
*/
template<int Rows=Dynamic, int Cols=Rows, class Precision=double>
@@ -25,14 +36,19 @@
public:
/// Construct the %QR decomposition of a matrix. This
initialises the class, and
/// performs the decomposition immediately.
+ /// @param m The matrix to decompose
+ /// @param p Whether or not to perform pivoting
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())
+ QR_Lapack(const Matrix<R,C,P,B>& m, bool p=0)
+ :copy(m),tau(square_size()), Q(square_size(), square_size()),
do_pivoting(p), pivot(Zeros(square_size()))
{
+ //pivot is set to all zeros, which means all columns
are free columns
+ //and can take part in column pivoting.
+
compute();
}
- ///Return L
+ ///Return R
const Matrix<Rows, Cols, Precision, ColMajor>& get_R()
{
return copy;
@@ -44,6 +60,12 @@
return Q;
}
+ ///Return the permutation vector. The definition is that column
\f$i\f$ of A is
+ ///column \f$P(i)\f$ of \f$QR\f$.
+ const Vector<Cols, int>& get_P()
+ {
+ return pivot;
+ }
private:
@@ -58,14 +80,22 @@
Precision size;
+ //Set up the pivot vector
+ if(do_pivoting)
+ pivot = Zeros;
+ else
+ for(int i=0; i < pivot.size(); i++)
+ pivot[i] = i+1;
+
+
//Compute the working space
- geqrf_(&M, &N, copy.get_data_ptr(), &lda,
tau.get_data_ptr(), &size, &LWORK, &INFO);
+ geqp3_(&M, &N, copy.get_data_ptr(), &lda,
pivot.get_data_ptr(), 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);
+ geqp3_(&M, &N, copy.get_data_ptr(), &lda,
pivot.get_data_ptr(), tau.get_data_ptr(), work, &LWORK, &INFO);
if(INFO < 0)
@@ -93,19 +123,23 @@
for(int c=0; c<r; c++)
copy[r][c] = 0;
+ //Now fix the pivot matrix.
+ //We need to go from FORTRAN to C numbering.
+ for(int i=0; i < pivot.size(); i++)
+ pivot[i]--;
}
Matrix<Rows, Cols, Precision, ColMajor> copy;
Matrix<square_Size, square_Size, Precision, ColMajor> Q;
Vector<square_Size, Precision> tau;
+ Vector<Cols, int> pivot;
+ bool do_pivoting;
int square_size()
{
return std::min(copy.num_rows(), copy.num_cols());
}
-
-
};
}
Index: lapack.h
===================================================================
RCS file: /cvsroot/toon/TooN/lapack.h,v
retrieving revision 1.17
retrieving revision 1.18
diff -u -b -r1.17 -r1.18
--- lapack.h 7 Apr 2010 09:31:50 -0000 1.17
+++ lapack.h 7 Apr 2010 13:09:59 -0000 1.18
@@ -76,9 +76,9 @@
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);
+ // Computes a QR decomposition of a general rectangular matrix
with column pivoting
+ void sgeqp3_(int* M, int* N, float* A, int* LDA, int* JPVT,
float* TAU, float* WORK, int* LWORK, int* INFO );
+ void dgeqp3_(int* M, int* N, double* A, int* LDA, int* JPVT,
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 );
@@ -161,14 +161,14 @@
}
//QR decomposition
- inline void geqrf_(int *m, int *n, float *a, int *lda, float *tau,
float *work, int *lwork, int *info)
+ inline void geqp3_(int* M, int* N, float* A, int* LDA, int* JPVT,
float* TAU, float* WORK, int* LWORK, int* INFO )
{
- sgeqrf_(m, n, a, lda, tau, work, lwork, info);
+ sgeqp3_(M, N, A, LDA, JPVT, TAU, WORK, LWORK, INFO);
}
- inline void geqrf_(int *m, int *n, double *a, int *lda, double *tau,
double *work, int *lwork, int *info)
+ inline void geqp3_(int* M, int* N, double* A, int* LDA, int* JPVT,
double* TAU, double* WORK, int* LWORK, int* INFO )
{
- dgeqrf_(m, n, a, lda, tau, work, lwork, info);
+ dgeqp3_(M, N, A, LDA, JPVT, TAU, WORK, LWORK, INFO);
}
inline void orgqr_(int* M,int* N,int* K, float* A, int* LDA, float*
TAU, float* WORK, int* LWORK, int* INFO )
| [Prev in Thread] |
Current Thread |
[Next in Thread] |
- [Toon-members] TooN QR_Lapack.h lapack.h,
Edward Rosten <=