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 #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 } else { // upper and lower triangular parts 00068 L(i,j) = U(j,i) = 0; 00069 L(j,i) = A(j,i); 00070 U(i,j) = A(i,j); 00071 } 00072 } 00073 } 00074 00075 p = p - 1; // Fortran counts from 1 00076 00077 return (info==0); 00078 } 00079 00080 // Slower than not using LAPACK when matrix size smaller than approx 20. 00081 bool lu(const cmat &X, cmat &L, cmat &U, ivec &p) 00082 { 00083 it_assert_debug(X.rows() == X.cols(), "lu: matrix is not quadratic"); 00084 //int m, n, lda, info; 00085 //m = n = lda = X.rows(); 00086 int m = X.rows(), info; 00087 00088 cmat A(X); 00089 L.set_size(m, m, false); 00090 U.set_size(m, m, false); 00091 p.set_size(m, false); 00092 00093 zgetrf_(&m, &m, A._data(), &m, p._data(), &info); 00094 00095 for (int i=0; i<m; i++) { 00096 for (int j=i; j<m; j++) { 00097 if (i == j) { // diagonal 00098 L(i,j) = 1; 00099 U(i,j) = A(i,j); 00100 } else { // upper and lower triangular parts 00101 L(i,j) = U(j,i) = 0; 00102 L(j,i) = A(j,i); 00103 U(i,j) = A(i,j); 00104 } 00105 } 00106 } 00107 00108 p = p - 1; // Fortran counts from 1 00109 00110 return (info==0); 00111 } 00112 00113 #else 00114 00115 bool lu(const mat &X, mat &L, mat &U, ivec &p) 00116 { 00117 it_error("LAPACK library is needed to use lu() function"); 00118 return false; 00119 } 00120 00121 bool lu(const cmat &X, cmat &L, cmat &U, ivec &p) 00122 { 00123 it_error("LAPACK library is needed to use lu() function"); 00124 return false; 00125 } 00126 00127 #endif // HAVE_LAPACK 00128 00129 00130 void interchange_permutations(vec &b, const ivec &p) 00131 { 00132 it_assert(b.size() == p.size(),"interchange_permutations(): dimension mismatch"); 00133 double temp; 00134 00135 for (int k=0; k<b.size(); k++) { 00136 temp=b(k); 00137 b(k)=b(p(k)); 00138 b(p(k))=temp; 00139 } 00140 } 00141 00142 bmat permutation_matrix(const ivec &p) 00143 { 00144 it_assert (p.size() > 0, "permutation_matrix(): vector must have nonzero size"); 00145 int n = p.size(), k; 00146 bmat P, identity; 00147 bvec row_k, row_pk; 00148 identity=eye_b(n); 00149 00150 00151 for (k=n-1; k>=0; k--) { 00152 // swap rows k and p(k) in identity 00153 row_k=identity.get_row(k); 00154 row_pk=identity.get_row(p(k)); 00155 identity.set_row(k, row_pk); 00156 identity.set_row(p(k), row_k); 00157 00158 if (k == n-1) { 00159 P=identity; 00160 } else { 00161 P*=identity; 00162 } 00163 00164 // swap back 00165 identity.set_row(k, row_k); 00166 identity.set_row(p(k), row_pk); 00167 } 00168 return P; 00169 } 00170 00171 } // namespace itpp
Generated on Sun Sep 14 18:52:25 2008 for IT++ by Doxygen 1.5.6