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
Generated on Wed Jul 27 2011 16:27:04 for IT++ by Doxygen 1.7.4