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 
39 #include <itpp/base/algebra/lu.h>
40 #include <itpp/base/specmat.h>
41 
42 
43 namespace itpp
44 {
45 
46 #if defined(HAVE_LAPACK)
47 
48 bool 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.
82 bool 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 
117 bool 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 
123 bool 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 
132 void 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 
144 bmat permutation_matrix(const ivec &p)
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
SourceForge Logo

Generated on Sat Jul 6 2013 10:54:19 for IT++ by Doxygen 1.8.2