IT++ Logo
lu.cpp
Go to the documentation of this file.
00001 
00029 #ifndef _MSC_VER
00030 #  include <itpp/config.h>
00031 #else
00032 #  include <itpp/config_msvc.h>
00033 #endif
00034 
00035 #if defined(HAVE_LAPACK)
00036 #  include <itpp/base/algebra/lapack.h>
00037 #endif
00038 
00039 #include <itpp/base/algebra/lu.h>
00040 #include <itpp/base/specmat.h>
00041 
00042 
00043 namespace itpp
00044 {
00045 
00046 #if defined(HAVE_LAPACK)
00047 
00048 bool lu(const mat &X, mat &L, mat &U, ivec &p)
00049 {
00050   it_assert_debug(X.rows() == X.cols(), "lu: matrix is not quadratic");
00051   //int m, n, lda, info;
00052   //m = n = lda = X.rows();
00053   int m = X.rows(), info;
00054 
00055   mat A(X);
00056   L.set_size(m, m, false);
00057   U.set_size(m, m, false);
00058   p.set_size(m, false);
00059 
00060   dgetrf_(&m, &m, A._data(), &m, p._data(), &info);
00061 
00062   for (int i = 0; i < m; i++) {
00063     for (int j = i; j < m; j++) {
00064       if (i == j) { // diagonal
00065         L(i, j) = 1;
00066         U(i, j) = A(i, j);
00067       }
00068       else { // upper and lower triangular parts
00069         L(i, j) = U(j, i) = 0;
00070         L(j, i) = A(j, i);
00071         U(i, j) = A(i, j);
00072       }
00073     }
00074   }
00075 
00076   p = p - 1; // Fortran counts from 1
00077 
00078   return (info == 0);
00079 }
00080 
00081 // Slower than not using LAPACK when matrix size smaller than approx 20.
00082 bool lu(const cmat &X, cmat &L, cmat &U, ivec &p)
00083 {
00084   it_assert_debug(X.rows() == X.cols(), "lu: matrix is not quadratic");
00085   //int m, n, lda, info;
00086   //m = n = lda = X.rows();
00087   int m = X.rows(), info;
00088 
00089   cmat A(X);
00090   L.set_size(m, m, false);
00091   U.set_size(m, m, false);
00092   p.set_size(m, false);
00093 
00094   zgetrf_(&m, &m, A._data(), &m, p._data(), &info);
00095 
00096   for (int i = 0; i < m; i++) {
00097     for (int j = i; j < m; j++) {
00098       if (i == j) { // diagonal
00099         L(i, j) = 1;
00100         U(i, j) = A(i, j);
00101       }
00102       else { // upper and lower triangular parts
00103         L(i, j) = U(j, i) = 0;
00104         L(j, i) = A(j, i);
00105         U(i, j) = A(i, j);
00106       }
00107     }
00108   }
00109 
00110   p = p - 1; // Fortran counts from 1
00111 
00112   return (info == 0);
00113 }
00114 
00115 #else
00116 
00117 bool lu(const mat &X, mat &L, mat &U, ivec &p)
00118 {
00119   it_error("LAPACK library is needed to use lu() function");
00120   return false;
00121 }
00122 
00123 bool lu(const cmat &X, cmat &L, cmat &U, ivec &p)
00124 {
00125   it_error("LAPACK library is needed to use lu() function");
00126   return false;
00127 }
00128 
00129 #endif // HAVE_LAPACK
00130 
00131 
00132 void interchange_permutations(vec &b, const ivec &p)
00133 {
00134   it_assert(b.size() == p.size(), "interchange_permutations(): dimension mismatch");
00135   double temp;
00136 
00137   for (int k = 0; k < b.size(); k++) {
00138     temp = b(k);
00139     b(k) = b(p(k));
00140     b(p(k)) = temp;
00141   }
00142 }
00143 
00144 bmat permutation_matrix(const ivec &p)
00145 {
00146   it_assert(p.size() > 0, "permutation_matrix(): vector must have nonzero size");
00147   int n = p.size(), k;
00148   bmat P, identity;
00149   bvec row_k, row_pk;
00150   identity = eye_b(n);
00151 
00152 
00153   for (k = n - 1; k >= 0; k--) {
00154     // swap rows k and p(k) in identity
00155     row_k = identity.get_row(k);
00156     row_pk = identity.get_row(p(k));
00157     identity.set_row(k, row_pk);
00158     identity.set_row(p(k), row_k);
00159 
00160     if (k == n - 1) {
00161       P = identity;
00162     }
00163     else {
00164       P *= identity;
00165     }
00166 
00167     // swap back
00168     identity.set_row(k, row_k);
00169     identity.set_row(p(k), row_pk);
00170   }
00171   return P;
00172 }
00173 
00174 } // namespace itpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
SourceForge Logo

Generated on Wed Jul 27 2011 16:27:04 for IT++ by Doxygen 1.7.4