Bayes++ Bayesian Filtering Classes  Release 2014.5 - Copyright (c) 2003,2004,2005,2006,2011,2012,2014 Michael Stevens
uLAPACK.hpp
Go to the documentation of this file.
1 /*
2  * Bayes++ the Bayesian Filtering Library
3  * Copyright (c) 2002 Michael Stevens
4  * See accompanying Bayes++.htm for terms and conditions of use.
5  *
6  * $Id$
7  */
8 
9 /*
10  * uBLAS to LAPACK Interface
11  * Very basic we only functions for information_root_filter supported
12  */
13 
14 /* Filter Matrix Namespace */
15 namespace Bayesian_filter_matrix
16 {
17 namespace LAPACK {
18 
19 namespace rawLAPACK {
20  // LAPACK routines as C callable functions
21  extern "C" {
22 
23  void dgeqrf_(
24  const int & m,
25  const int & n,
26  double da[],
27  const int & lda,
28  double dtau[],
29  double dwork[],
30  const int& ldwork,
31  int& info);
32  void sgeqrf_(
33  const int & m,
34  const int & n,
35  float da[],
36  const int & lda,
37  float dtau[],
38  float dwork[],
39  const int& ldwork,
40  int& info);
41 
42  void dgetrs_(
43  const char& transa,
44  const int& n,
45  const int& nrhs,
46  const double da[],
47  const int& lda,
48  int ipivot[],
49  double db[],
50  const int& ldb,
51  int& info);
52  void sgetrs_(
53  const char& transa,
54  const int& n,
55  const int& nrhs,
56  const float da[],
57  const int& lda,
58  int ipivot[],
59  float db[],
60  const int& ldb,
61  int& info);
62 
63  void dgetrf_(
64  const int& m,
65  const int& n,
66  double da[],
67  const int& lda,
68  int ipivot[],
69  int& info);
70  void sgetrf_(
71  const int& m,
72  const int& n,
73  float da[],
74  const int& lda,
75  int ipivot[],
76  int& info);
77 
78  }// extern "C"
79 
80  // Type overloads for C++
81  inline void geqrf( const int & m, const int & n, double da[], const int & lda, double dtau[], double dwork[], const int& ldwork, int& info)
82  {
83  dgeqrf_(m,n,da,lda,dtau,dwork,ldwork,info);
84  }
85  inline void geqrf( const int & m, const int & n, float da[], const int & lda, float dtau[], float dwork[], const int& ldwork, int& info)
86  {
87  sgeqrf_(m,n,da,lda,dtau,dwork,ldwork,info);
88  }
89  inline void getrs( const char& transa, const int& n, const int& nrhs, const double da[], const int& lda, int ipivot[], double db[], const int& ldb, int& info)
90  {
91  dgetrs_(transa,n,nrhs,da,lda,ipivot,db,ldb,info);
92  }
93  inline void getrs( const char& transa, const int& n, const int& nrhs, const float da[], const int& lda, int ipivot[], float db[], const int& ldb, int& info)
94  {
95  sgetrs_(transa,n,nrhs,da,lda,ipivot,db,ldb,info);
96  }
97  inline void getrf( const int& m, const int& n, double da[], const int& lda, int ipivot[], int& info)
98  {
99  dgetrf_(m,n,da,lda,ipivot,info);
100  }
101  inline void getrf( const int& m, const int& n, float da[], const int& lda, int ipivot[], int& info)
102  {
103  sgetrf_(m,n,da,lda,ipivot,info);
104  }
105 }// namespace rawLAPACK
106 
107 
108 /* Support types */
109 typedef boost::numeric::ublas::vector<int> pivot_t;
112 
113 /* LAPACK Interface*/
114 
115 
116 
117 // QR Factorization of a MxN General Matrix A.
118 // a (IN/OUT - matrix(M,N)) On entry, the coefficient matrix A. On exit , the upper triangle and diagonal is the min(M,N) by N upper triangular matrix R. The lower triangle, together with the tau vector, is the orthogonal matrix Q as a product of min(M,N) elementary reflectors.
119 // tau (OUT - vector (min(M,N))) Vector of the same numerical type as A. The scalar factors of the elementary reflectors.
120 // info (OUT - int)
121 // 0 : function completed normally
122 // < 0 : The ith argument, where i = abs(return value) had an illegal value.
123 int geqrf (matrix_t& a, vector_t& tau)
124 {
125  int _m = int(a.size1());
126  int _n = int(a.size2());
127  int _lda = int(a.size1());
128  int _info;
129 
130  // make_sure tau's size is greater than or equal to min(m,n)
131  if (int(tau.size()) < (_n<_m ? _n : _m) )
132  return -104;
133 
134  int ldwork = _n*_n;
135  vector_t dwork(ldwork);
136  rawLAPACK::geqrf (_m, _n, a.data().begin(), _lda, tau.data().begin(), dwork.data().begin(), ldwork, _info);
137 
138  return _info;
139 }
140 
141 // LU factorization of a general matrix A.
142 // Computes an LU factorization of a general M-by-N matrix A using
143 // partial pivoting with row interchanges. Factorization has the form
144 // A = P*L*U.
145 // a (IN/OUT - matrix(M,N)) On entry, the coefficient matrix A to be factored. On exit, the factors L and U from the factorization A = P*L*U.
146 // ipivot (OUT - vector(min(M,N))) Integer vector. The row i of A was interchanged with row IPIV(i).
147 // info (OUT - int)
148 // 0 : successful exit
149 // < 0 : If INFO = -i, then the i-th argument had an illegal value.
150 // > 0 : If INFO = i, then U(i,i) is exactly zero. The factorization has been completed, but the factor U is exactly singular, and division by zero will occur if it is used to solve a system of equations.
151 int getrf (matrix_t& a, pivot_t& ipivot)
152 {
153  matrix_t::value_type* _a = a.data().begin();
154  int _m = int(a.size1());
155  int _n = int(a.size2());
156  int _lda = _m; // minor size
157  int _info;
158 
159  rawLAPACK::getrf (_m, _n, _a, _lda, ipivot.data().begin(), _info);
160 
161  return _info;
162 }
163 
164 // Solution to a system using LU factorization
165 // Solves a system of linear equations A*X = B with a general NxN
166 // matrix A using the LU factorization computed by GETRF.
167 // transa (IN - char) 'T' for the transpose of A, 'N' otherwise.
168 // a (IN - matrix(M,N)) The factors L and U from the factorization A = P*L*U as computed by GETRF.
169 // ipivot (IN - vector(min(M,N))) Integer vector. The pivot indices from GETRF; row i of A was interchanged with row IPIV(i).
170 // b (IN/OUT - matrix(ldb,NRHS)) Matrix of same numerical type as A. On entry, the right hand side matrix B. On exit, the solution matrix X.
171 //
172 // info (OUT - int)
173 // 0 : function completed normally
174 // < 0 : The ith argument, where i = abs(return value) had an illegal value.
175 // > 0 : if INFO = i, U(i,i) is exactly zero; the matrix is singular and its inverse could not be computed.
176 int getrs (char transa, matrix_t& a,
177  pivot_t& ipivot, matrix_t& b)
178 {
179  matrix_t::value_type* _a = a.data().begin();
180  int a_n = int(a.size1());
181  int _lda = a_n;
182  int p_n = int(ipivot.size());
183 
184  matrix_t::value_type* _b = b.data().begin();
185  int b_n = int(b.size1());
186  int _ldb = b_n;
187  int _nrhs = int(b.size2()); /* B's size2 is the # of vectors on rhs */
188 
189  if (a_n != b_n) /*Test to see if AX=B has correct dimensions */
190  return -101;
191  if (p_n < a_n) /*Check to see if ipivot is big enough */
192  return -102;
193 
194  int _info;
195  rawLAPACK::getrs (transa, a_n, _nrhs, _a, _lda, ipivot.data().begin(),
196  _b, _ldb, _info);
197 
198  return _info;
199 }
200 
201 }//namespace LAPACK
202 }//namespace Bayesian_filter_matrix
FMVec< detail::BaseDenseVector > DenseVec
Definition: uBLASmatrix.hpp:333
void sgetrs_(const char &transa, const int &n, const int &nrhs, const float da[], const int &lda, int ipivot[], float db[], const int &ldb, int &info)
void getrf(const int &m, const int &n, double da[], const int &lda, int ipivot[], int &info)
Definition: uLAPACK.hpp:97
void sgeqrf_(const int &m, const int &n, float da[], const int &lda, float dtau[], float dwork[], const int &ldwork, int &info)
boost::numeric::ublas::vector< int > pivot_t
Definition: uLAPACK.hpp:109
void getrs(const char &transa, const int &n, const int &nrhs, const double da[], const int &lda, int ipivot[], double db[], const int &ldb, int &info)
Definition: uLAPACK.hpp:89
FMMatrix< detail::BaseDenseColMatrix > DenseColMatrix
Definition: uBLASmatrix.hpp:336
Bayesian_filter_matrix::DenseColMatrix matrix_t
Definition: uLAPACK.hpp:111
Bayesian_filter_matrix::DenseVec vector_t
Definition: uLAPACK.hpp:110
Definition: matSup.cpp:33
void dgetrf_(const int &m, const int &n, double da[], const int &lda, int ipivot[], int &info)
void geqrf(const int &m, const int &n, double da[], const int &lda, double dtau[], double dwork[], const int &ldwork, int &info)
Definition: uLAPACK.hpp:81
void dgeqrf_(const int &m, const int &n, double da[], const int &lda, double dtau[], double dwork[], const int &ldwork, int &info)
void sgetrf_(const int &m, const int &n, float da[], const int &lda, int ipivot[], int &info)
void dgetrs_(const char &transa, const int &n, const int &nrhs, const double da[], const int &lda, int ipivot[], double db[], const int &ldb, int &info)