Ravelin
LinAlg.h
1 /****************************************************************************
2  * Copyright 2013 Evan Drumwright
3  * This library is distributed under the terms of the Apache V2.0
4  * License (obtainable from http://www.apache.org/licenses/LICENSE-2.0).
5  ****************************************************************************/
6 
7 #ifndef LINALG
8 #error This class is not to be included by the user directly. Use LinAlgf.h or LinAlgd.h instead.
9 #endif
10 
12 
15 class LINALG
16 {
17  public:
18  enum SVD { eSVD1, eSVD2 };
19 
20  private:
21  static REAL log2(REAL x);
22  static void lartg_(REAL* F, REAL* G, REAL* CS, REAL* SN, REAL* R);
23  static void gtsv_(INTEGER* N, INTEGER* NRHS, REAL* DL, REAL* D, REAL* DU, REAL* B, INTEGER* LDB, INTEGER* INFO);
24  static void trtrs_(char* UPLO, char* TRANS, INTEGER* N, INTEGER* NRHS, REAL* AP, INTEGER* LDA, REAL* B, INTEGER* LDB, INTEGER* INFO);
25  void ormqr_(char* SIDE, char* TRANS, INTEGER* M, INTEGER* N, INTEGER* K, REAL* A, INTEGER* LDA, REAL* TAU, REAL* C, INTEGER* LDC, INTEGER* INFO);
26  static void sptrf_(char* UPLO, INTEGER* N, REAL* AP, INTEGER* IPIV, INTEGER* INFO);
27  static void sptrs_(char* UPLO, INTEGER* N, INTEGER* NRHS, REAL* AP, INTEGER* IPIV, REAL* B, INTEGER* LDB, INTEGER* INFO);
28  void gelsd_(INTEGER* M, INTEGER* N, INTEGER* NRHS, REAL* A, INTEGER*
29 LDA, REAL* B, INTEGER* LDB, REAL* RCOND, INTEGER* INFO);
30  void sysv_(char* UPLO, INTEGER* N, INTEGER* NRHS, REAL* A, INTEGER* LDA, INTEGER* IPIV, REAL* B, INTEGER* LDB, INTEGER* INFO);
31  void gesdd_(char* JOBZ, INTEGER* M, INTEGER* N, REAL* A, INTEGER* LDA, REAL* S, REAL* U, INTEGER* LDU, REAL* V, INTEGER* LDVT, INTEGER* INFO);
32  void gesvd_(char* JOBU, char* JOBV, INTEGER* M, INTEGER* N, REAL* A, INTEGER* LDA, REAL* S, REAL* U, INTEGER* LDU, REAL* V, INTEGER* LDVT, INTEGER* INFO);
33  void syevd_(char* JOBZ, char* UPLO, INTEGER* N, REAL* A, INTEGER* LDA, REAL* EVALS, INTEGER* INFO);
34  static void gesv_(INTEGER* N, INTEGER* NRHS, REAL* A, INTEGER* LDA, INTEGER* IPIV, REAL* X, INTEGER* LDX, INTEGER* INFO);
35  static void potrs_(char* UPLO, INTEGER* N, INTEGER* NRHS, REAL* A, INTEGER* LDA, REAL* B, INTEGER* LDB, INTEGER* INFO);
36  static void potri_(char* UPLO, INTEGER* N, REAL* A, INTEGER* LDA, INTEGER* INFO);
37  static void potrf_(char* UPLO, INTEGER* N, REAL* A, INTEGER* LDA, INTEGER* INFO);
38  static void posv_(char* UPLO, INTEGER* N, INTEGER* NRHS, REAL* A, INTEGER* LDA, REAL* B, INTEGER* LDB, INTEGER* INFO);
39  void sytri_(char* UPLO, INTEGER* N, REAL* A, INTEGER* LDA, INTEGER* IPIV, INTEGER* INFO);
40  void sytrf_(char* UPLO, INTEGER* N, REAL* A, INTEGER* LDA, INTEGER* IPIV, INTEGER* INFO);
41  void geqp3_(INTEGER* M, INTEGER* N, REAL* A, INTEGER* LDA, INTEGER* JPVT, REAL* TAU, INTEGER* INFO);
42  void geqrf_(INTEGER* M, INTEGER* N, REAL* A, INTEGER* LDA, REAL* TAU, INTEGER* INFO);
43  static void getrf_(INTEGER* M, INTEGER* N, REAL* A, INTEGER* LDA, INTEGER* IPIV, INTEGER* INFO);
44  void getri_(INTEGER* N, REAL* A, INTEGER* LDA, INTEGER* IPIV, INTEGER* INFO);
45  static void getrs_(char* TRANS, INTEGER* N, INTEGER* NRHS, REAL* A, INTEGER* LDA, INTEGER* IPIV, REAL* B, INTEGER* LDB, INTEGER* INFO);
46  void orgqr_(INTEGER* M, INTEGER* N, INTEGER* K, REAL* A, INTEGER* LDA, REAL* TAU, INTEGER* INFO);
47 
48  public:
49  void compress();
50  void free_memory();
51  static void factor_LDL(MATRIXN& M, std::vector<int>& IPIV);
52  MATRIXN& pseudo_invert(MATRIXN& A, REAL tol=(REAL) -1.0);
53  static void givens(REAL a, REAL b, REAL& c, REAL& s);
54  static MATRIX2 givens(REAL c, REAL s);
55  static void householder(REAL alpha, const VECTORN& x, REAL& tau, VECTORN& v);
56  static VECTORN& solve_sparse_direct(const SPARSEMATRIXN& A, const VECTORN& b, Transposition trans, VECTORN& x);
57  static MATRIXN& solve_sparse_direct(const SPARSEMATRIXN& A, const MATRIXN& B, Transposition trans, MATRIXN& X);
58  void update_QR_rank1(MATRIXN& Q, MATRIXN& R, const VECTORN& u, const VECTORN& v);
59  void update_QR_delete_cols(MATRIXN& Q, MATRIXN& R, unsigned k, unsigned p);
60  void update_QR_insert_cols(MATRIXN& Q, MATRIXN& R, MATRIXN& U, unsigned k);
61  void update_QR_insert_rows(MATRIXN& Q, MATRIXN& R, MATRIXN& U, unsigned k);
62  void update_QR_delete_rows(MATRIXN& Q, MATRIXN& R, unsigned k, unsigned p);
63 
65  FastThreadable<MATRIXN> workM, workM2;
66 
68  FastThreadable<MATRIXN> U;
69 
71  FastThreadable<std::vector<INTEGER> > pivwork;
72 
74  FastThreadable<MATRIXN> V;
75 
77  FastThreadable<VECTORN> S;
78 
80  FastThreadable<VECTORN> workv, workv2;
81 
83  FastThreadable<std::vector<INTEGER> > iworkv;
84 
85  // include templated routines here...
86  #include "LinAlg.inl"
87 
88 }; // end class
89 
static void givens(REAL a, REAL b, REAL &c, REAL &s)
Less robust least squares solver (solves Ax = b)
Definition: LinAlg.cpp:184
A general 2x2 matrix.
Definition: Matrix2.h:16
FastThreadable< MATRIXN > U
work matrix (for SVD)
Definition: LinAlg.h:68
void compress()
Compresses all memory.
Definition: LinAlg.cpp:17
FastThreadable< VECTORN > S
work vector (for SVD/eigenvalues)
Definition: LinAlg.h:77
Linear algebra routines.
Definition: LinAlg.h:15
void free_memory()
Frees all allocated memory.
Definition: LinAlg.cpp:2
A generic, possibly non-square matrix.
Definition: MatrixN.h:18
void update_QR_delete_cols(MATRIXN &Q, MATRIXN &R, unsigned k, unsigned p)
Updates a QR factorization by deleting p columns starting at column idx k.
Definition: LinAlg.cpp:236
FastThreadable< VECTORN > workv
work vectors (LAPACK routines)
Definition: LinAlg.h:80
A sparse matrix.
Definition: SparseMatrixN.h:2
FastThreadable< std::vector< INTEGER > > pivwork
work STL integer vector (pivoting)
Definition: LinAlg.h:71
void update_QR_insert_rows(MATRIXN &Q, MATRIXN &R, MATRIXN &U, unsigned k)
Updates a QR factorization by inserting a block of rows, starting at index k.
Definition: LinAlg.cpp:501
void update_QR_delete_rows(MATRIXN &Q, MATRIXN &R, unsigned k, unsigned p)
Updates a QR factorization by deleting a block of rows.
Definition: LinAlg.cpp:714
A generic N-dimensional floating point vector.
Definition: VectorN.h:16
void update_QR_insert_cols(MATRIXN &Q, MATRIXN &R, MATRIXN &U, unsigned k)
Updates a QR factorization by inserting one or more columns at column idx k.
Definition: LinAlg.cpp:380
static void householder(REAL alpha, const VECTORN &x, REAL &tau, VECTORN &v)
Computes a householder vector.
Definition: LinAlg.cpp:203
FastThreadable< std::vector< INTEGER > > iworkv
work STL integer vector (LAPACK routines)
Definition: LinAlg.h:83
FastThreadable< MATRIXN > V
work matrix (for SVD)
Definition: LinAlg.h:74
FastThreadable< MATRIXN > workM
work matrices
Definition: LinAlg.h:65
void update_QR_rank1(MATRIXN &Q, MATRIXN &R, const VECTORN &u, const VECTORN &v)
Updates a QR factorization by a rank-1 update.
Definition: LinAlg.cpp:223