IT++ Logo
lu.cpp
Go to the documentation of this file.
1
29#ifndef _MSC_VER
30# include <itpp/config.h>
31#else
32# include <itpp/config_msvc.h>
33#endif
34
35#if defined(HAVE_LAPACK)
36# include <itpp/base/algebra/lapack.h>
37#endif
38
40#include <itpp/base/specmat.h>
41
42
43namespace itpp
44{
45
46#if defined(HAVE_LAPACK)
47
48bool lu(const mat &X, mat &L, mat &U, ivec &p)
49{
50 it_assert_debug(X.rows() == X.cols(), "lu: matrix is not quadratic");
51 //int m, n, lda, info;
52 //m = n = lda = X.rows();
53 int m = X.rows(), info;
54
55 mat A(X);
56 L.set_size(m, m, false);
57 U.set_size(m, m, false);
58 p.set_size(m, false);
59
60 dgetrf_(&m, &m, A._data(), &m, p._data(), &info);
61
62 for (int i = 0; i < m; i++) {
63 for (int j = i; j < m; j++) {
64 if (i == j) { // diagonal
65 L(i, j) = 1;
66 U(i, j) = A(i, j);
67 }
68 else { // upper and lower triangular parts
69 L(i, j) = U(j, i) = 0;
70 L(j, i) = A(j, i);
71 U(i, j) = A(i, j);
72 }
73 }
74 }
75
76 p = p - 1; // Fortran counts from 1
77
78 return (info == 0);
79}
80
81// Slower than not using LAPACK when matrix size smaller than approx 20.
82bool lu(const cmat &X, cmat &L, cmat &U, ivec &p)
83{
84 it_assert_debug(X.rows() == X.cols(), "lu: matrix is not quadratic");
85 //int m, n, lda, info;
86 //m = n = lda = X.rows();
87 int m = X.rows(), info;
88
89 cmat A(X);
90 L.set_size(m, m, false);
91 U.set_size(m, m, false);
92 p.set_size(m, false);
93
94 zgetrf_(&m, &m, A._data(), &m, p._data(), &info);
95
96 for (int i = 0; i < m; i++) {
97 for (int j = i; j < m; j++) {
98 if (i == j) { // diagonal
99 L(i, j) = 1;
100 U(i, j) = A(i, j);
101 }
102 else { // upper and lower triangular parts
103 L(i, j) = U(j, i) = 0;
104 L(j, i) = A(j, i);
105 U(i, j) = A(i, j);
106 }
107 }
108 }
109
110 p = p - 1; // Fortran counts from 1
111
112 return (info == 0);
113}
114
115#else
116
117bool lu(const mat &X, mat &L, mat &U, ivec &p)
118{
119 it_error("LAPACK library is needed to use lu() function");
120 return false;
121}
122
123bool lu(const cmat &X, cmat &L, cmat &U, ivec &p)
124{
125 it_error("LAPACK library is needed to use lu() function");
126 return false;
127}
128
129#endif // HAVE_LAPACK
130
131
132void interchange_permutations(vec &b, const ivec &p)
133{
134 it_assert(b.size() == p.size(), "interchange_permutations(): dimension mismatch");
135 double temp;
136
137 for (int k = 0; k < b.size(); k++) {
138 temp = b(k);
139 b(k) = b(p(k));
140 b(p(k)) = temp;
141 }
142}
143
145{
146 it_assert(p.size() > 0, "permutation_matrix(): vector must have nonzero size");
147 int n = p.size(), k;
148 bmat P, identity;
149 bvec row_k, row_pk;
150 identity = eye_b(n);
151
152
153 for (k = n - 1; k >= 0; k--) {
154 // swap rows k and p(k) in identity
155 row_k = identity.get_row(k);
156 row_pk = identity.get_row(p(k));
157 identity.set_row(k, row_pk);
158 identity.set_row(p(k), row_k);
159
160 if (k == n - 1) {
161 P = identity;
162 }
163 else {
164 P *= identity;
165 }
166
167 // swap back
168 identity.set_row(k, row_k);
169 identity.set_row(p(k), row_pk);
170 }
171 return P;
172}
173
174} // namespace itpp
#define it_error(s)
Abort unconditionally.
Definition: itassert.h:126
#define it_assert_debug(t, s)
Abort if t is not true and NDEBUG is not defined.
Definition: itassert.h:107
#define it_assert(t, s)
Abort if t is not true.
Definition: itassert.h:94
bool lu(const mat &X, mat &L, mat &U, ivec &p)
LU factorisation of real matrix.
Definition: lu.cpp:117
bmat permutation_matrix(const ivec &p)
Make permutation matrix P from the interchange permutation vector p.
Definition: lu.cpp:144
void interchange_permutations(vec &b, const ivec &p)
Makes swapping of vector b according to the interchange permutation vector p.
Definition: lu.cpp:132
ITPP_EXPORT bmat eye_b(int size)
A Binary (size,size) unit matrix.
Definitions of LU factorisation functions.
Mat< bin > bmat
bin matrix
Definition: mat.h:508
itpp namespace
Definition: itmex.h:37
Definitions of special vectors and matrices.
SourceForge Logo

Generated on Tue Jan 24 2023 00:00:00 for IT++ by Doxygen 1.9.6