31 #ifndef TOON_INCLUDE_LU_H
32 #define TOON_INCLUDE_LU_H
68 template <
int Size=-1,
class Precision=
double>
74 template<
int S1,
int S2,
class Base>
81 template<
int S1,
int S2,
class Base>
89 int lda = m.num_rows();
96 std::cerr <<
"error in LU, INFO was " <<
my_info << std::endl;
102 template <
int Rows,
int NRHS,
class Base>
109 int M=rhs.num_cols();
110 int N=
my_lu.num_rows();
112 int lda=
my_lu.num_rows();
113 int ldb=rhs.num_cols();
114 trsm_(
"R",
"U",
"N",
"N",&M,&N,&alpha,&
my_lu[0][0],&lda,&result[0][0],&ldb);
115 trsm_(
"R",
"L",
"N",
"U",&M,&N,&alpha,&
my_lu[0][0],&lda,&result[0][0],&ldb);
118 for(
int i=N-1; i>=0; i--){
119 const int swaprow =
my_IPIV[i]-1;
120 for(
int j=0; j<NRHS; j++){
121 Precision temp = result[i][j];
122 result[i][j] = result[swaprow][j];
123 result[swaprow][j] = temp;
131 template <
int Rows,
class Base>
139 int N=
my_lu.num_rows();
141 int lda=
my_lu.num_rows();
143 trsm_(
"R",
"U",
"N",
"N",&M,&N,&alpha,&
my_lu[0][0],&lda,&result[0],&ldb);
144 trsm_(
"R",
"L",
"N",
"U",&M,&N,&alpha,&
my_lu[0][0],&lda,&result[0],&ldb);
147 for(
int i=N-1; i>=0; i--){
148 const int swaprow =
my_IPIV[i]-1;
149 Precision temp = result[i];
150 result[i] = result[swaprow];
151 result[swaprow] = temp;
160 int N =
my_lu.num_rows();
161 int lda=
my_lu.num_rows();
166 Precision* WORK =
new Precision[lwork];
182 for(
int i=0; i<
my_lu.num_rows()-1; i++){
194 for (
int i=0; i<
my_lu.num_rows(); i++){
LU(const Matrix< S1, S2, Precision, Base > &m)
const Matrix< Size, Size, Precision > & get_lu() const
void trsm_(const char *SIDE, const char *UPLO, const char *TRANSA, const char *DIAG, int *M, int *N, float *alpha, float *A, int *lda, float *B, int *ldb)
void compute(const Matrix< S1, S2, Precision, Base > &m)
Perform the LU decompsition of another matrix.
Everything lives inside this namespace.
Precision determinant() const
Calculate the determinant of the matrix.
Matrix< Size, NRHS, Precision > backsub(const Matrix< Rows, NRHS, Precision, Base > &rhs)
static void test(int s1, int s2)
void getri_(int *N, double *A, int *lda, int *IPIV, double *WORK, int *lwork, int *INFO)
Matrix< Size, Size, Precision > my_lu
int get_info() const
Get the LAPACK info.
Vector< Size, Precision > backsub(const Vector< Rows, Precision, Base > &rhs)
Vector< Size, int > my_IPIV
Matrix< Size, Size, Precision > get_inverse()
void getrf_(int *M, int *N, float *A, int *lda, int *IPIV, int *INFO)