IT++ Logo

lu.cpp

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

Generated on Wed Mar 2 2011 22:04:50 for IT++ by Doxygen 1.7.3