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
Generated on Sun Jul 26 08:54:24 2009 for IT++ by Doxygen 1.5.9