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 #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
SourceForge Logo

Generated on Sat Apr 19 10:43:50 2008 for IT++ by Doxygen 1.5.5