IT++ Logo Newcom Logo

lu.cpp

Go to the documentation of this file.
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
SourceForge Logo

Generated on Wed Mar 21 12:21:34 2007 for IT++ by Doxygen 1.4.7