00001 00033 #ifndef _MSC_VER 00034 # include <itpp/config.h> 00035 #else 00036 # include <itpp/config_msvc.h> 00037 #endif 00038 00039 #if defined(HAVE_LAPACK) 00040 # include <itpp/base/lapack.h> 00041 #endif 00042 00043 #include <itpp/base/lu.h> 00044 #include <itpp/base/specmat.h> 00045 00046 00047 namespace itpp { 00048 00049 #if defined(HAVE_LAPACK) 00050 00051 bool lu(const mat &X, mat &L, mat &U, ivec &p) 00052 { 00053 it_assert1(X.rows() == X.cols(), "lu: matrix is not quadratic"); 00054 //int m, n, lda, info; 00055 //m = n = lda = X.rows(); 00056 int m = X.rows(), info; 00057 00058 mat A(X); 00059 L.set_size(m, m, false); 00060 U.set_size(m, m, false); 00061 p.set_size(m, false); 00062 00063 dgetrf_(&m, &m, A._data(), &m, p._data(), &info); 00064 00065 for (int i=0; i<m; i++) { 00066 for (int j=i; j<m; j++) { 00067 if (i == j) { // diagonal 00068 L(i,j) = 1; 00069 U(i,j) = A(i,j); 00070 } else { // upper and lower triangular parts 00071 L(i,j) = U(j,i) = 0; 00072 L(j,i) = A(j,i); 00073 U(i,j) = A(i,j); 00074 } 00075 } 00076 } 00077 00078 p = p - 1; // Fortran counts from 1 00079 00080 return (info==0); 00081 } 00082 00083 // Slower than not using LAPACK when matrix size smaller than approx 20. 00084 bool lu(const cmat &X, cmat &L, cmat &U, ivec &p) 00085 { 00086 it_assert1(X.rows() == X.cols(), "lu: matrix is not quadratic"); 00087 //int m, n, lda, info; 00088 //m = n = lda = X.rows(); 00089 int m = X.rows(), info; 00090 00091 cmat A(X); 00092 L.set_size(m, m, false); 00093 U.set_size(m, m, false); 00094 p.set_size(m, false); 00095 00096 zgetrf_(&m, &m, A._data(), &m, p._data(), &info); 00097 00098 for (int i=0; i<m; i++) { 00099 for (int j=i; j<m; j++) { 00100 if (i == j) { // diagonal 00101 L(i,j) = 1; 00102 U(i,j) = A(i,j); 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 } 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 Thu Apr 19 14:14:57 2007 for IT++ by Doxygen 1.5.1