Auxlib


Classes

class  auxlib
 wrapper for accessing external functions defined in ATLAS, LAPACK or BLAS libraries More...

Functions

template<typename eT >
static void auxlib::inv_noalias (Mat< eT > &out, const Mat< eT > &X)
 immediate matrix inverse
template<typename eT >
static void auxlib::inv_inplace (Mat< eT > &X)
 immediate inplace matrix inverse
template<typename eT >
static eT auxlib::det (const Mat< eT > &X)
 immediate determinant of a matrix using ATLAS or LAPACK
template<typename eT >
static void auxlib::lu (Mat< eT > &L, Mat< eT > &U, podarray< int > &ipiv, const Mat< eT > &X_orig)
 immediate LU decomposition of a matrix using ATLAS or LAPACK
template<typename eT >
static void auxlib::lu (Mat< eT > &L, Mat< eT > &U, Mat< eT > &P, const Mat< eT > &X)
template<typename eT >
static void auxlib::lu (Mat< eT > &L, Mat< eT > &U, const Mat< eT > &X)
template<typename eT >
static void auxlib::eig_sym (Col< eT > &eigval, const Mat< eT > &A)
 immediate eigenvalues of a symmetric real matrix using LAPACK
template<typename T >
static void auxlib::eig_sym (Col< T > &eigval, const Mat< std::complex< T > > &A)
 immediate eigenvalues of a hermitian complex matrix using LAPACK
template<typename eT >
static void auxlib::eig_sym (Col< eT > &eigval, Mat< eT > &eigvec, const Mat< eT > &A)
 immediate eigenvalues and eigenvectors of a symmetric real matrix using LAPACK
template<typename T >
static void auxlib::eig_sym (Col< T > &eigval, Mat< std::complex< T > > &eigvec, const Mat< std::complex< T > > &A)
 immediate eigenvalues and eigenvectors of a hermitian complex matrix using LAPACK
template<typename T >
void auxlib::eig_gen (Col< std::complex< T > > &eigval, Mat< T > &l_eigvec, Mat< T > &r_eigvec, const Mat< T > &A, const char side)
 Eigenvalues and eigenvectors of a general square real matrix using LAPACK. The argument 'side' specifies which eigenvectors should be calculated (see code for mode details).
template<typename T >
static void auxlib::eig_gen (Col< std::complex< T > > &eigval, Mat< std::complex< T > > &l_eigvec, Mat< std::complex< T > > &r_eigvec, const Mat< std::complex< T > > &A, const char side)
 Eigenvalues and eigenvectors of a general square complex matrix using LAPACK The argument 'side' specifies which eigenvectors should be calculated (see code for mode details).
template<typename eT >
static bool auxlib::chol (Mat< eT > &out, const Mat< eT > &X)
template<typename eT >
static bool auxlib::qr (Mat< eT > &Q, Mat< eT > &R, const Mat< eT > &X)
template<typename eT >
static bool auxlib::svd (Col< eT > &S, const Mat< eT > &X)
template<typename T >
static bool auxlib::svd (Col< T > &S, const Mat< std::complex< T > > &X)
template<typename eT >
static bool auxlib::svd (Mat< eT > &U, Col< eT > &S, Mat< eT > &V, const Mat< eT > &X)
template<typename T >
static bool auxlib::svd (Mat< std::complex< T > > &U, Col< T > &S, Mat< std::complex< T > > &V, const Mat< std::complex< T > > &X)
template<typename eT >
static bool auxlib::solve (Mat< eT > &out, const Mat< eT > &A, const Mat< eT > &B)
 Solve a system of linear equations Assumes that A.n_rows = A.n_cols and B.n_rows = A.n_rows.
template<typename eT >
static bool auxlib::solve_od (Mat< eT > &out, const Mat< eT > &A, const Mat< eT > &B)
 Solve an over-determined system. Assumes that A.n_rows > A.n_cols and B.n_rows = A.n_rows.
template<typename eT >
static bool auxlib::solve_ud (Mat< eT > &out, const Mat< eT > &A, const Mat< eT > &B)
 Solve an under-determined system. Assumes that A.n_rows < A.n_cols and B.n_rows = A.n_rows.

Function Documentation

template<typename eT >
void auxlib::inv_noalias ( Mat< eT > &  out,
const Mat< eT > &  X 
) [inline, static, inherited]

immediate matrix inverse

Definition at line 25 of file auxlib_meat.hpp.

References arma_check(), arma_stop(), Mat< eT >::at(), atlas::clapack_getrf(), atlas::clapack_getri(), Mat< eT >::colptr(), det(), lapack::getrf_(), lapack::getri_(), max(), podarray< T1 >::memptr(), Mat< eT >::memptr(), Mat< eT >::n_cols, Mat< eT >::n_rows, podarray< T1 >::set_size(), Mat< eT >::set_size(), and tmp_real().

Referenced by op_inv::apply().

00026   {
00027   arma_extra_debug_sigprint();
00028   
00029   switch(X.n_rows)
00030     {
00031     case 1:
00032       {
00033       out.set_size(1,1);
00034       out[0] = eT(1) / X[0];
00035       };
00036       break;
00037       
00038     case 2:
00039       {
00040       out.set_size(2,2);
00041       
00042       const eT a = X.at(0,0);
00043       const eT b = X.at(0,1);
00044       const eT c = X.at(1,0);
00045       const eT d = X.at(1,1);
00046       
00047       const eT k = eT(1) / (a*d - b*c);
00048       
00049       out.at(0,0) =  d*k;
00050       out.at(0,1) = -b*k;
00051       out.at(1,0) = -c*k;
00052       out.at(1,1) =  a*k;
00053       };
00054       break;
00055     
00056     case 3:
00057       {
00058       out.set_size(3,3);
00059       
00060       const eT* X_col0 = X.colptr(0);
00061       const eT a11 = X_col0[0];
00062       const eT a21 = X_col0[1];
00063       const eT a31 = X_col0[2];
00064       
00065       const eT* X_col1 = X.colptr(1);
00066       const eT a12 = X_col1[0];
00067       const eT a22 = X_col1[1];
00068       const eT a32 = X_col1[2];
00069       
00070       const eT* X_col2 = X.colptr(2);
00071       const eT a13 = X_col2[0];
00072       const eT a23 = X_col2[1];
00073       const eT a33 = X_col2[2];
00074       
00075       const eT k = eT(1) / ( a11*(a33*a22 - a32*a23) - a21*(a33*a12-a32*a13) + a31*(a23*a12 - a22*a13) );
00076       
00077       
00078       eT* out_col0 = out.colptr(0);
00079       out_col0[0] =  (a33*a22 - a32*a23) * k;
00080       out_col0[1] = -(a33*a21 - a31*a23) * k;
00081       out_col0[2] =  (a32*a21 - a31*a22) * k;
00082       
00083       eT* out_col1 = out.colptr(1);
00084       out_col1[0] = -(a33*a12 - a32*a13) * k;
00085       out_col1[1] =  (a33*a11 - a31*a13) * k;
00086       out_col1[2] = -(a32*a11 - a31*a12) * k;
00087       
00088       eT* out_col2 = out.colptr(2);
00089       out_col2[0] =  (a23*a12 - a22*a13) * k;
00090       out_col2[1] = -(a23*a11 - a21*a13) * k;
00091       out_col2[2] =  (a22*a11 - a21*a12) * k;
00092       };
00093       break;
00094       
00095     case 4:
00096       {
00097       out.set_size(4,4);
00098       
00099       out.at(0,0) = X.at(1,2)*X.at(2,3)*X.at(3,1) - X.at(1,3)*X.at(2,2)*X.at(3,1) + X.at(1,3)*X.at(2,1)*X.at(3,2) - X.at(1,1)*X.at(2,3)*X.at(3,2) - X.at(1,2)*X.at(2,1)*X.at(3,3) + X.at(1,1)*X.at(2,2)*X.at(3,3);
00100       out.at(1,0) = X.at(1,3)*X.at(2,2)*X.at(3,0) - X.at(1,2)*X.at(2,3)*X.at(3,0) - X.at(1,3)*X.at(2,0)*X.at(3,2) + X.at(1,0)*X.at(2,3)*X.at(3,2) + X.at(1,2)*X.at(2,0)*X.at(3,3) - X.at(1,0)*X.at(2,2)*X.at(3,3);
00101       out.at(2,0) = X.at(1,1)*X.at(2,3)*X.at(3,0) - X.at(1,3)*X.at(2,1)*X.at(3,0) + X.at(1,3)*X.at(2,0)*X.at(3,1) - X.at(1,0)*X.at(2,3)*X.at(3,1) - X.at(1,1)*X.at(2,0)*X.at(3,3) + X.at(1,0)*X.at(2,1)*X.at(3,3);
00102       out.at(3,0) = X.at(1,2)*X.at(2,1)*X.at(3,0) - X.at(1,1)*X.at(2,2)*X.at(3,0) - X.at(1,2)*X.at(2,0)*X.at(3,1) + X.at(1,0)*X.at(2,2)*X.at(3,1) + X.at(1,1)*X.at(2,0)*X.at(3,2) - X.at(1,0)*X.at(2,1)*X.at(3,2);
00103       
00104       out.at(0,1) = X.at(0,3)*X.at(2,2)*X.at(3,1) - X.at(0,2)*X.at(2,3)*X.at(3,1) - X.at(0,3)*X.at(2,1)*X.at(3,2) + X.at(0,1)*X.at(2,3)*X.at(3,2) + X.at(0,2)*X.at(2,1)*X.at(3,3) - X.at(0,1)*X.at(2,2)*X.at(3,3);
00105       out.at(1,1) = X.at(0,2)*X.at(2,3)*X.at(3,0) - X.at(0,3)*X.at(2,2)*X.at(3,0) + X.at(0,3)*X.at(2,0)*X.at(3,2) - X.at(0,0)*X.at(2,3)*X.at(3,2) - X.at(0,2)*X.at(2,0)*X.at(3,3) + X.at(0,0)*X.at(2,2)*X.at(3,3);
00106       out.at(2,1) = X.at(0,3)*X.at(2,1)*X.at(3,0) - X.at(0,1)*X.at(2,3)*X.at(3,0) - X.at(0,3)*X.at(2,0)*X.at(3,1) + X.at(0,0)*X.at(2,3)*X.at(3,1) + X.at(0,1)*X.at(2,0)*X.at(3,3) - X.at(0,0)*X.at(2,1)*X.at(3,3);
00107       out.at(3,1) = X.at(0,1)*X.at(2,2)*X.at(3,0) - X.at(0,2)*X.at(2,1)*X.at(3,0) + X.at(0,2)*X.at(2,0)*X.at(3,1) - X.at(0,0)*X.at(2,2)*X.at(3,1) - X.at(0,1)*X.at(2,0)*X.at(3,2) + X.at(0,0)*X.at(2,1)*X.at(3,2);
00108       
00109       out.at(0,2) = X.at(0,2)*X.at(1,3)*X.at(3,1) - X.at(0,3)*X.at(1,2)*X.at(3,1) + X.at(0,3)*X.at(1,1)*X.at(3,2) - X.at(0,1)*X.at(1,3)*X.at(3,2) - X.at(0,2)*X.at(1,1)*X.at(3,3) + X.at(0,1)*X.at(1,2)*X.at(3,3);
00110       out.at(1,2) = X.at(0,3)*X.at(1,2)*X.at(3,0) - X.at(0,2)*X.at(1,3)*X.at(3,0) - X.at(0,3)*X.at(1,0)*X.at(3,2) + X.at(0,0)*X.at(1,3)*X.at(3,2) + X.at(0,2)*X.at(1,0)*X.at(3,3) - X.at(0,0)*X.at(1,2)*X.at(3,3);
00111       out.at(2,2) = X.at(0,1)*X.at(1,3)*X.at(3,0) - X.at(0,3)*X.at(1,1)*X.at(3,0) + X.at(0,3)*X.at(1,0)*X.at(3,1) - X.at(0,0)*X.at(1,3)*X.at(3,1) - X.at(0,1)*X.at(1,0)*X.at(3,3) + X.at(0,0)*X.at(1,1)*X.at(3,3);
00112       out.at(3,2) = X.at(0,2)*X.at(1,1)*X.at(3,0) - X.at(0,1)*X.at(1,2)*X.at(3,0) - X.at(0,2)*X.at(1,0)*X.at(3,1) + X.at(0,0)*X.at(1,2)*X.at(3,1) + X.at(0,1)*X.at(1,0)*X.at(3,2) - X.at(0,0)*X.at(1,1)*X.at(3,2);
00113       
00114       out.at(0,3) = X.at(0,3)*X.at(1,2)*X.at(2,1) - X.at(0,2)*X.at(1,3)*X.at(2,1) - X.at(0,3)*X.at(1,1)*X.at(2,2) + X.at(0,1)*X.at(1,3)*X.at(2,2) + X.at(0,2)*X.at(1,1)*X.at(2,3) - X.at(0,1)*X.at(1,2)*X.at(2,3);
00115       out.at(1,3) = X.at(0,2)*X.at(1,3)*X.at(2,0) - X.at(0,3)*X.at(1,2)*X.at(2,0) + X.at(0,3)*X.at(1,0)*X.at(2,2) - X.at(0,0)*X.at(1,3)*X.at(2,2) - X.at(0,2)*X.at(1,0)*X.at(2,3) + X.at(0,0)*X.at(1,2)*X.at(2,3);
00116       out.at(2,3) = X.at(0,3)*X.at(1,1)*X.at(2,0) - X.at(0,1)*X.at(1,3)*X.at(2,0) - X.at(0,3)*X.at(1,0)*X.at(2,1) + X.at(0,0)*X.at(1,3)*X.at(2,1) + X.at(0,1)*X.at(1,0)*X.at(2,3) - X.at(0,0)*X.at(1,1)*X.at(2,3);
00117       out.at(3,3) = X.at(0,1)*X.at(1,2)*X.at(2,0) - X.at(0,2)*X.at(1,1)*X.at(2,0) + X.at(0,2)*X.at(1,0)*X.at(2,1) - X.at(0,0)*X.at(1,2)*X.at(2,1) - X.at(0,1)*X.at(1,0)*X.at(2,2) + X.at(0,0)*X.at(1,1)*X.at(2,2);      
00118       
00119       out /= det(X);
00120       };
00121       break;
00122       
00123     default:
00124       {
00125       #if defined(ARMA_USE_ATLAS)
00126         {
00127         out = X;
00128         podarray<int> ipiv(out.n_rows);
00129         
00130         int info = atlas::clapack_getrf(atlas::CblasColMajor, out.n_rows, out.n_cols, out.memptr(), out.n_rows, ipiv.memptr());
00131         
00132         if(info == 0)
00133           {
00134           info = atlas::clapack_getri(atlas::CblasColMajor, out.n_rows, out.memptr(), out.n_rows, ipiv.memptr());
00135           }
00136         
00137         arma_check( (info > 0), "auxlib::inv_noalias(): matrix appears to be singular" );
00138         }
00139       #elif defined(ARMA_USE_LAPACK)
00140         {
00141         out = X;
00142         
00143         int n_rows = out.n_rows;
00144         int n_cols = out.n_cols;
00145         int info   = 0;
00146         
00147         podarray<int> ipiv(out.n_rows);
00148         
00149         // 84 was empirically found -- it is the maximum value suggested by LAPACK (as provided by ATLAS v3.6)
00150         // based on tests with various matrix types on 32-bit and 64-bit machines
00151         //
00152         // the "work" array is deliberately long so that a secondary (time-consuming)
00153         // memory allocation is avoided, if possible
00154         
00155         int work_len = (std::max)(1, n_rows*84);
00156         podarray<eT> work(work_len);
00157         
00158         lapack::getrf_(&n_rows, &n_cols, out.memptr(), &n_rows, ipiv.memptr(), &info);
00159         
00160         if(info == 0)
00161           {
00162           // query for optimum size of work_len
00163           
00164           int work_len_tmp = -1;
00165           lapack::getri_(&n_rows, out.memptr(), &n_rows, ipiv.memptr(), work.memptr(), &work_len_tmp, &info);
00166           
00167           if(info == 0)
00168             {
00169             int proposed_work_len = static_cast<int>(auxlib::tmp_real(work[0]));
00170             
00171             // if necessary, allocate more memory
00172             if(work_len < proposed_work_len)
00173               {
00174               work_len = proposed_work_len;
00175               work.set_size(work_len);
00176               }
00177             }
00178           
00179           lapack::getri_(&n_rows, out.memptr(), &n_rows, ipiv.memptr(), work.memptr(), &work_len, &info);
00180           }
00181         
00182         arma_check( (info > 0), "auxlib::inv_noalias(): matrix appears to be singular" );
00183         }
00184       #else
00185         {
00186         arma_stop("auxlib::inv_noalias(): need ATLAS or LAPACK library");
00187         }
00188       #endif
00189       };
00190     }
00191   }

template<typename eT >
void auxlib::inv_inplace ( Mat< eT > &  X  )  [inline, static, inherited]

immediate inplace matrix inverse

Definition at line 199 of file auxlib_meat.hpp.

References arma_check(), arma_stop(), Mat< eT >::at(), atlas::clapack_getrf(), atlas::clapack_getri(), Mat< eT >::colptr(), det(), lapack::getrf_(), lapack::getri_(), max(), podarray< T1 >::memptr(), Mat< eT >::memptr(), Mat< eT >::n_cols, Mat< eT >::n_rows, podarray< T1 >::set_size(), and tmp_real().

Referenced by op_inv::apply().

00200   {
00201   arma_extra_debug_sigprint();
00202   
00203   // for more info, see:
00204   // http://www.dr-lex.34sp.com/random/matrix_inv.html
00205   // http://www.cvl.iis.u-tokyo.ac.jp/~miyazaki/tech/teche23.html
00206   // http://www.euclideanspace.com/maths/algebra/matrix/functions/inverse/fourD/index.htm
00207   // http://www.geometrictools.com//LibFoundation/Mathematics/Wm4Matrix4.inl
00208   
00209   switch(X.n_rows)
00210     {
00211     case 1:
00212       {
00213       X[0] = eT(1) / X[0];
00214       };
00215       break;
00216       
00217     case 2:
00218       {
00219       const eT a = X.at(0,0);
00220       const eT b = X.at(0,1);
00221       const eT c = X.at(1,0);
00222       const eT d = X.at(1,1);
00223       
00224       const eT k = eT(1) / (a*d - b*c);
00225       
00226       X.at(0,0) =  d*k;
00227       X.at(0,1) = -b*k;
00228       X.at(1,0) = -c*k;
00229       X.at(1,1) =  a*k;
00230       };
00231       break;
00232     
00233     case 3:
00234       {
00235       eT* X_col0 = X.colptr(0);
00236       eT* X_col1 = X.colptr(1);
00237       eT* X_col2 = X.colptr(2);
00238       
00239       const eT a11 = X_col0[0];
00240       const eT a21 = X_col0[1];
00241       const eT a31 = X_col0[2];
00242       
00243       const eT a12 = X_col1[0];
00244       const eT a22 = X_col1[1];
00245       const eT a32 = X_col1[2];
00246       
00247       const eT a13 = X_col2[0];
00248       const eT a23 = X_col2[1];
00249       const eT a33 = X_col2[2];
00250       
00251       const eT k = eT(1) / ( a11*(a33*a22 - a32*a23) - a21*(a33*a12-a32*a13) + a31*(a23*a12 - a22*a13) );
00252       
00253       X_col0[0] =  (a33*a22 - a32*a23) * k;
00254       X_col0[1] = -(a33*a21 - a31*a23) * k;
00255       X_col0[2] =  (a32*a21 - a31*a22) * k;
00256       
00257       X_col1[0] = -(a33*a12 - a32*a13) * k;
00258       X_col1[1] =  (a33*a11 - a31*a13) * k;
00259       X_col1[2] = -(a32*a11 - a31*a12) * k;
00260       
00261       X_col2[0] =  (a23*a12 - a22*a13) * k;
00262       X_col2[1] = -(a23*a11 - a21*a13) * k;
00263       X_col2[2] =  (a22*a11 - a21*a12) * k;
00264       };
00265       break;
00266       
00267     case 4:
00268       {
00269       const Mat<eT> A(X);
00270       
00271       X.at(0,0) = A.at(1,2)*A.at(2,3)*A.at(3,1) - A.at(1,3)*A.at(2,2)*A.at(3,1) + A.at(1,3)*A.at(2,1)*A.at(3,2) - A.at(1,1)*A.at(2,3)*A.at(3,2) - A.at(1,2)*A.at(2,1)*A.at(3,3) + A.at(1,1)*A.at(2,2)*A.at(3,3);
00272       X.at(1,0) = A.at(1,3)*A.at(2,2)*A.at(3,0) - A.at(1,2)*A.at(2,3)*A.at(3,0) - A.at(1,3)*A.at(2,0)*A.at(3,2) + A.at(1,0)*A.at(2,3)*A.at(3,2) + A.at(1,2)*A.at(2,0)*A.at(3,3) - A.at(1,0)*A.at(2,2)*A.at(3,3);
00273       X.at(2,0) = A.at(1,1)*A.at(2,3)*A.at(3,0) - A.at(1,3)*A.at(2,1)*A.at(3,0) + A.at(1,3)*A.at(2,0)*A.at(3,1) - A.at(1,0)*A.at(2,3)*A.at(3,1) - A.at(1,1)*A.at(2,0)*A.at(3,3) + A.at(1,0)*A.at(2,1)*A.at(3,3);
00274       X.at(3,0) = A.at(1,2)*A.at(2,1)*A.at(3,0) - A.at(1,1)*A.at(2,2)*A.at(3,0) - A.at(1,2)*A.at(2,0)*A.at(3,1) + A.at(1,0)*A.at(2,2)*A.at(3,1) + A.at(1,1)*A.at(2,0)*A.at(3,2) - A.at(1,0)*A.at(2,1)*A.at(3,2);
00275       
00276       X.at(0,1) = A.at(0,3)*A.at(2,2)*A.at(3,1) - A.at(0,2)*A.at(2,3)*A.at(3,1) - A.at(0,3)*A.at(2,1)*A.at(3,2) + A.at(0,1)*A.at(2,3)*A.at(3,2) + A.at(0,2)*A.at(2,1)*A.at(3,3) - A.at(0,1)*A.at(2,2)*A.at(3,3);
00277       X.at(1,1) = A.at(0,2)*A.at(2,3)*A.at(3,0) - A.at(0,3)*A.at(2,2)*A.at(3,0) + A.at(0,3)*A.at(2,0)*A.at(3,2) - A.at(0,0)*A.at(2,3)*A.at(3,2) - A.at(0,2)*A.at(2,0)*A.at(3,3) + A.at(0,0)*A.at(2,2)*A.at(3,3);
00278       X.at(2,1) = A.at(0,3)*A.at(2,1)*A.at(3,0) - A.at(0,1)*A.at(2,3)*A.at(3,0) - A.at(0,3)*A.at(2,0)*A.at(3,1) + A.at(0,0)*A.at(2,3)*A.at(3,1) + A.at(0,1)*A.at(2,0)*A.at(3,3) - A.at(0,0)*A.at(2,1)*A.at(3,3);
00279       X.at(3,1) = A.at(0,1)*A.at(2,2)*A.at(3,0) - A.at(0,2)*A.at(2,1)*A.at(3,0) + A.at(0,2)*A.at(2,0)*A.at(3,1) - A.at(0,0)*A.at(2,2)*A.at(3,1) - A.at(0,1)*A.at(2,0)*A.at(3,2) + A.at(0,0)*A.at(2,1)*A.at(3,2);
00280       
00281       X.at(0,2) = A.at(0,2)*A.at(1,3)*A.at(3,1) - A.at(0,3)*A.at(1,2)*A.at(3,1) + A.at(0,3)*A.at(1,1)*A.at(3,2) - A.at(0,1)*A.at(1,3)*A.at(3,2) - A.at(0,2)*A.at(1,1)*A.at(3,3) + A.at(0,1)*A.at(1,2)*A.at(3,3);
00282       X.at(1,2) = A.at(0,3)*A.at(1,2)*A.at(3,0) - A.at(0,2)*A.at(1,3)*A.at(3,0) - A.at(0,3)*A.at(1,0)*A.at(3,2) + A.at(0,0)*A.at(1,3)*A.at(3,2) + A.at(0,2)*A.at(1,0)*A.at(3,3) - A.at(0,0)*A.at(1,2)*A.at(3,3);
00283       X.at(2,2) = A.at(0,1)*A.at(1,3)*A.at(3,0) - A.at(0,3)*A.at(1,1)*A.at(3,0) + A.at(0,3)*A.at(1,0)*A.at(3,1) - A.at(0,0)*A.at(1,3)*A.at(3,1) - A.at(0,1)*A.at(1,0)*A.at(3,3) + A.at(0,0)*A.at(1,1)*A.at(3,3);
00284       X.at(3,2) = A.at(0,2)*A.at(1,1)*A.at(3,0) - A.at(0,1)*A.at(1,2)*A.at(3,0) - A.at(0,2)*A.at(1,0)*A.at(3,1) + A.at(0,0)*A.at(1,2)*A.at(3,1) + A.at(0,1)*A.at(1,0)*A.at(3,2) - A.at(0,0)*A.at(1,1)*A.at(3,2);
00285       
00286       X.at(0,3) = A.at(0,3)*A.at(1,2)*A.at(2,1) - A.at(0,2)*A.at(1,3)*A.at(2,1) - A.at(0,3)*A.at(1,1)*A.at(2,2) + A.at(0,1)*A.at(1,3)*A.at(2,2) + A.at(0,2)*A.at(1,1)*A.at(2,3) - A.at(0,1)*A.at(1,2)*A.at(2,3);
00287       X.at(1,3) = A.at(0,2)*A.at(1,3)*A.at(2,0) - A.at(0,3)*A.at(1,2)*A.at(2,0) + A.at(0,3)*A.at(1,0)*A.at(2,2) - A.at(0,0)*A.at(1,3)*A.at(2,2) - A.at(0,2)*A.at(1,0)*A.at(2,3) + A.at(0,0)*A.at(1,2)*A.at(2,3);
00288       X.at(2,3) = A.at(0,3)*A.at(1,1)*A.at(2,0) - A.at(0,1)*A.at(1,3)*A.at(2,0) - A.at(0,3)*A.at(1,0)*A.at(2,1) + A.at(0,0)*A.at(1,3)*A.at(2,1) + A.at(0,1)*A.at(1,0)*A.at(2,3) - A.at(0,0)*A.at(1,1)*A.at(2,3);
00289       X.at(3,3) = A.at(0,1)*A.at(1,2)*A.at(2,0) - A.at(0,2)*A.at(1,1)*A.at(2,0) + A.at(0,2)*A.at(1,0)*A.at(2,1) - A.at(0,0)*A.at(1,2)*A.at(2,1) - A.at(0,1)*A.at(1,0)*A.at(2,2) + A.at(0,0)*A.at(1,1)*A.at(2,2);      
00290       
00291       X /= det(A);
00292       };
00293       break;
00294       
00295     default:
00296       {
00297       #if defined(ARMA_USE_ATLAS)
00298         {
00299         Mat<eT>& out = X;
00300         podarray<int> ipiv(out.n_rows);
00301         
00302         int info = atlas::clapack_getrf(atlas::CblasColMajor, out.n_rows, out.n_cols, out.memptr(), out.n_rows, ipiv.memptr());
00303         
00304         if(info == 0)
00305           {
00306           info = atlas::clapack_getri(atlas::CblasColMajor, out.n_rows, out.memptr(), out.n_rows, ipiv.memptr());
00307           }
00308         
00309         arma_check( (info > 0), "auxlib::inv_inplace(): matrix appears to be singular" );
00310         }
00311       #elif defined(ARMA_USE_LAPACK)
00312         {
00313         Mat<eT>& out = X;
00314         
00315         int n_rows = out.n_rows;
00316         int n_cols = out.n_cols;
00317         int info   = 0;
00318         
00319         podarray<int> ipiv(out.n_rows);
00320         
00321         // 84 was empirically found -- it is the maximum value suggested by LAPACK (as provided by ATLAS v3.6)
00322         // based on tests with various matrix types on 32-bit and 64-bit machines
00323         //
00324         // the "work" array is deliberately long so that a secondary (time-consuming)
00325         // memory allocation is avoided, if possible
00326         
00327         int work_len = (std::max)(1, n_rows*84);
00328         podarray<eT> work(work_len);
00329         
00330         lapack::getrf_(&n_rows, &n_cols, out.memptr(), &n_rows, ipiv.memptr(), &info);
00331         
00332         if(info == 0)
00333           {
00334           // query for optimum size of work_len
00335           
00336           int work_len_tmp = -1;
00337           lapack::getri_(&n_rows, out.memptr(), &n_rows, ipiv.memptr(), work.memptr(), &work_len_tmp, &info);
00338           
00339           if(info == 0)
00340             {
00341             int proposed_work_len = static_cast<int>(auxlib::tmp_real(work[0]));
00342             
00343             // if necessary, allocate more memory
00344             if(work_len < proposed_work_len)
00345               {
00346               work_len = proposed_work_len;
00347               work.set_size(work_len);
00348               }
00349             }
00350           
00351           lapack::getri_(&n_rows, out.memptr(), &n_rows, ipiv.memptr(), work.memptr(), &work_len, &info);
00352           }
00353         
00354         arma_check( (info > 0), "auxlib::inv_noalias(): matrix appears to be singular" );
00355         }
00356       #else
00357         {
00358         arma_stop("auxlib::inv_inplace(): need ATLAS or LAPACK");
00359         }
00360       #endif
00361       }
00362     
00363     }
00364   
00365   }

template<typename eT >
eT auxlib::det ( const Mat< eT > &  X  )  [inline, static, inherited]

immediate determinant of a matrix using ATLAS or LAPACK

Definition at line 372 of file auxlib_meat.hpp.

References arma_stop(), Mat< eT >::at(), atlas::clapack_getrf(), Mat< eT >::colptr(), lapack::getrf_(), podarray< T1 >::mem, podarray< T1 >::memptr(), Mat< eT >::memptr(), Mat< eT >::n_cols, and Mat< eT >::n_rows.

Referenced by det(), inv_inplace(), and inv_noalias().

00373   {
00374   arma_extra_debug_sigprint();
00375   
00376   switch(X.n_rows)
00377     {
00378     case 0:
00379       return 0.0;
00380     
00381     case 1:
00382       return X[0];
00383     
00384     case 2:
00385       return (X.at(0,0)*X.at(1,1) - X.at(0,1)*X.at(1,0));
00386     
00387     case 3:
00388       {
00389       const eT* a_col0 = X.colptr(0);
00390       const eT  a11    = a_col0[0];
00391       const eT  a21    = a_col0[1];
00392       const eT  a31    = a_col0[2];
00393       
00394       const eT* a_col1 = X.colptr(1);
00395       const eT  a12    = a_col1[0];
00396       const eT  a22    = a_col1[1];
00397       const eT  a32    = a_col1[2];
00398       
00399       const eT* a_col2 = X.colptr(2);
00400       const eT  a13    = a_col2[0];
00401       const eT  a23    = a_col2[1];
00402       const eT  a33    = a_col2[2];
00403       
00404       return ( a11*(a33*a22 - a32*a23) - a21*(a33*a12-a32*a13) + a31*(a23*a12 - a22*a13) );
00405       
00406       // const double tmp1 = X.at(0,0) * X.at(1,1) * X.at(2,2);
00407       // const double tmp2 = X.at(0,1) * X.at(1,2) * X.at(2,0);
00408       // const double tmp3 = X.at(0,2) * X.at(1,0) * X.at(2,1);
00409       // const double tmp4 = X.at(2,0) * X.at(1,1) * X.at(0,2);
00410       // const double tmp5 = X.at(2,1) * X.at(1,2) * X.at(0,0);
00411       // const double tmp6 = X.at(2,2) * X.at(1,0) * X.at(0,1);
00412       // return (tmp1+tmp2+tmp3) - (tmp4+tmp5+tmp6);
00413       }
00414       
00415     case 4:
00416       {
00417       const eT val = \
00418           X.at(0,3) * X.at(1,2) * X.at(2,1) * X.at(3,0) \
00419         - X.at(0,2) * X.at(1,3) * X.at(2,1) * X.at(3,0) \
00420         - X.at(0,3) * X.at(1,1) * X.at(2,2) * X.at(3,0) \
00421         + X.at(0,1) * X.at(1,3) * X.at(2,2) * X.at(3,0) \
00422         + X.at(0,2) * X.at(1,1) * X.at(2,3) * X.at(3,0) \
00423         - X.at(0,1) * X.at(1,2) * X.at(2,3) * X.at(3,0) \
00424         - X.at(0,3) * X.at(1,2) * X.at(2,0) * X.at(3,1) \
00425         + X.at(0,2) * X.at(1,3) * X.at(2,0) * X.at(3,1) \
00426         + X.at(0,3) * X.at(1,0) * X.at(2,2) * X.at(3,1) \
00427         - X.at(0,0) * X.at(1,3) * X.at(2,2) * X.at(3,1) \
00428         - X.at(0,2) * X.at(1,0) * X.at(2,3) * X.at(3,1) \
00429         + X.at(0,0) * X.at(1,2) * X.at(2,3) * X.at(3,1) \
00430         + X.at(0,3) * X.at(1,1) * X.at(2,0) * X.at(3,2) \
00431         - X.at(0,1) * X.at(1,3) * X.at(2,0) * X.at(3,2) \
00432         - X.at(0,3) * X.at(1,0) * X.at(2,1) * X.at(3,2) \
00433         + X.at(0,0) * X.at(1,3) * X.at(2,1) * X.at(3,2) \
00434         + X.at(0,1) * X.at(1,0) * X.at(2,3) * X.at(3,2) \
00435         - X.at(0,0) * X.at(1,1) * X.at(2,3) * X.at(3,2) \
00436         - X.at(0,2) * X.at(1,1) * X.at(2,0) * X.at(3,3) \
00437         + X.at(0,1) * X.at(1,2) * X.at(2,0) * X.at(3,3) \
00438         + X.at(0,2) * X.at(1,0) * X.at(2,1) * X.at(3,3) \
00439         - X.at(0,0) * X.at(1,2) * X.at(2,1) * X.at(3,3) \
00440         - X.at(0,1) * X.at(1,0) * X.at(2,2) * X.at(3,3) \
00441         + X.at(0,0) * X.at(1,1) * X.at(2,2) * X.at(3,3) \
00442         ;
00443       
00444       return val;
00445       }
00446       
00447     default:  
00448       {
00449       #if defined(ARMA_USE_ATLAS)
00450         {
00451         Mat<eT> tmp = X;
00452         podarray<int> ipiv(tmp.n_rows);
00453         
00454         atlas::clapack_getrf(atlas::CblasColMajor, tmp.n_rows, tmp.n_cols, tmp.memptr(), tmp.n_rows, ipiv.memptr());
00455         
00456         // on output tmp appears to be L+U_alt, where U_alt is U with the main diagonal set to zero
00457         eT val = tmp.at(0,0);
00458         for(u32 i=1; i < tmp.n_rows; ++i)
00459           {
00460           val *= tmp.at(i,i);
00461           }
00462       
00463         eT sign = eT(1);
00464         for(u32 i=0; i < tmp.n_rows; ++i)
00465           {
00466           if(s32(i) != ipiv.mem[i])
00467             {
00468             sign *= eT(-1);
00469             }
00470           }
00471         
00472         return sign*val;
00473         }
00474       #elif defined(ARMA_USE_LAPACK)
00475         {
00476         Mat<eT> tmp = X;
00477         podarray<int> ipiv(tmp.n_rows);
00478         
00479         int info = 0;
00480         int n_rows = tmp.n_rows;
00481         int n_cols = tmp.n_cols;
00482         
00483         lapack::getrf_(&n_rows, &n_cols, tmp.memptr(), &n_rows, ipiv.memptr(), &info);
00484         
00485         // on output tmp appears to be L+U_alt, where U_alt is U with the main diagonal set to zero
00486         eT val = tmp.at(0,0);
00487         for(u32 i=1; i < tmp.n_rows; ++i)
00488           {
00489           val *= tmp.at(i,i);
00490           }
00491       
00492         eT sign = eT(1);
00493         for(u32 i=0; i < tmp.n_rows; ++i)
00494           {
00495           if(s32(i) != ipiv.mem[i])
00496             {
00497             sign *= eT(-1);
00498             }
00499           }
00500         
00501         return sign*val;
00502         }
00503       #else
00504         {
00505         arma_stop("auxlib::det(): need ATLAS or LAPACK library");
00506         return eT(0);
00507         }
00508       #endif
00509       }
00510     }
00511   }

template<typename eT >
void auxlib::lu ( Mat< eT > &  L,
Mat< eT > &  U,
podarray< int > &  ipiv,
const Mat< eT > &  X_orig 
) [inline, static, inherited]

immediate LU decomposition of a matrix using ATLAS or LAPACK

Definition at line 519 of file auxlib_meat.hpp.

References arma_stop(), Mat< eT >::at(), atlas::clapack_getrf(), lapack::getrf_(), podarray< T1 >::memptr(), Mat< eT >::memptr(), Mat< eT >::n_cols, Mat< eT >::n_rows, Mat< eT >::set_size(), and podarray< T1 >::set_size().

Referenced by lu(), and lu().

00520   {
00521   arma_extra_debug_sigprint();
00522   
00523   arma_debug_check( (&L == &U), "auxlib::lu(): L and U are the same object");
00524   
00525   unwrap_check< Mat<eT> > tmp1(X_orig, U);
00526   const Mat<eT>& X_tmp = tmp1.M;
00527   
00528   unwrap_check< Mat<eT> > tmp2(X_tmp, L);
00529   const Mat<eT>& X = tmp2.M;
00530   
00531   U = X;
00532   
00533   #if defined(ARMA_USE_ATLAS) || defined(ARMA_USE_LAPACK)
00534     {
00535     
00536     #if defined(ARMA_USE_ATLAS)
00537       {
00538       ipiv.set_size(U.n_rows);
00539     
00540       //int info = 
00541       atlas::clapack_getrf(atlas::CblasColMajor, U.n_rows, U.n_cols, U.memptr(), U.n_rows, ipiv.memptr());
00542       }
00543     #elif defined(ARMA_USE_LAPACK)
00544       {
00545       ipiv.set_size(U.n_rows);
00546       int info = 0;
00547       
00548       int n_rows = U.n_rows;
00549       int n_cols = U.n_cols;
00550       
00551       lapack::getrf_(&n_rows, &n_cols, U.memptr(), &n_rows, ipiv.memptr(), &info);
00552       }
00553     #endif
00554     
00555     
00556     L.set_size(U.n_rows, U.n_rows);
00557     
00558     for(u32 col=0; col<L.n_cols; ++col)
00559       {
00560       
00561       for(u32 row=0; row<col; ++row)
00562         {
00563         L.at(row,col) = eT(0);
00564         }
00565       
00566       L.at(col,col) = eT(1);
00567       
00568       for(u32 row=col+1; row<L.n_rows; ++row)
00569         {
00570         L.at(row,col) = U.at(row,col);
00571         U.at(row,col) = eT(0);
00572         }
00573       
00574       }
00575     }
00576   #else
00577     {
00578     arma_stop("auxlib::lu(): need ATLAS or LAPACK library");
00579     }
00580   #endif
00581   
00582   }

template<typename eT >
void auxlib::lu ( Mat< eT > &  L,
Mat< eT > &  U,
Mat< eT > &  P,
const Mat< eT > &  X 
) [inline, static, inherited]

Definition at line 589 of file auxlib_meat.hpp.

References eye(), lu(), podarray< T1 >::n_elem, Mat< eT >::set_size(), and Mat< eT >::swap_rows().

00590   {
00591   arma_extra_debug_sigprint();
00592   
00593   arma_debug_check( ( (&P == &L) || (&P == &U) ), "auxlib::lu(): P is the same object as L or U" );
00594   
00595   podarray<int> ipiv;
00596   auxlib::lu(L, U, ipiv, X);
00597   
00598   const u32 n = ipiv.n_elem;
00599   
00600   P.set_size(n,n);
00601   Mat<eT> ident = eye(n,n);
00602 
00603   for(u32 i=n; i>0; --i)
00604     {
00605     const u32 j = i-1;
00606     const u32 k = ipiv[j]-1;
00607     
00608     ident.swap_rows(j,k);
00609     
00610     if(i == n)
00611       P = ident;
00612     else
00613       P *= ident;
00614 
00615     ident.swap_rows(j,k);
00616     }
00617   
00618   }

template<typename eT >
void auxlib::lu ( Mat< eT > &  L,
Mat< eT > &  U,
const Mat< eT > &  X 
) [inline, static, inherited]

Definition at line 625 of file auxlib_meat.hpp.

References lu().

00626   {
00627   arma_extra_debug_sigprint();
00628   
00629   podarray<int> ipiv;
00630   auxlib::lu(L, U, ipiv, X);
00631   }

template<typename eT >
void auxlib::eig_sym ( Col< eT > &  eigval,
const Mat< eT > &  A 
) [inline, static, inherited]

immediate eigenvalues of a symmetric real matrix using LAPACK

Definition at line 639 of file auxlib_meat.hpp.

References arma_stop(), max(), podarray< T1 >::memptr(), Mat< eT >::memptr(), Mat< eT >::n_cols, Mat< eT >::n_rows, Col< eT >::set_size(), and lapack::syev_().

Referenced by eig_sym().

00640   {
00641   arma_extra_debug_sigprint();
00642   
00643   #if defined(ARMA_USE_LAPACK)
00644     {
00645     const unwrap_check<Mat<eT> > tmp(A_orig, eigval);
00646     const Mat<eT>& A = tmp.M;
00647     
00648     arma_debug_check( (A.n_rows != A.n_cols), "auxlib::eig_sym(): given matrix is not square");
00649     
00650     // rudimentary "better-than-nothing" test for symmetry
00651     //arma_debug_check( (A.at(A.n_rows-1, 0) != A.at(0, A.n_cols-1)), "auxlib::eig(): given matrix is not symmetric" );
00652 
00653     char jobz  = 'N';
00654     char uplo  = 'U';
00655     
00656     int n_rows = A.n_rows;
00657     int lwork  = (std::max)(1,3*n_rows-1);
00658     
00659     eigval.set_size(n_rows);
00660     podarray<eT> work(lwork);
00661   
00662     Mat<eT> A_copy = A;
00663     int info;
00664 
00665     arma_extra_debug_print("lapack::syev_()");
00666     lapack::syev_(&jobz, &uplo, &n_rows, A_copy.memptr(), &n_rows, eigval.memptr(), work.memptr(), &lwork, &info);
00667     }
00668   #else
00669     {
00670     arma_stop("auxlib::eig_sym(): need LAPACK library");
00671     }
00672   #endif
00673   }

template<typename T >
void auxlib::eig_sym ( Col< T > &  eigval,
const Mat< std::complex< T > > &  A 
) [inline, static, inherited]

immediate eigenvalues of a hermitian complex matrix using LAPACK

Definition at line 681 of file auxlib_meat.hpp.

References arma_stop(), lapack::heev_(), max(), podarray< T1 >::memptr(), Mat< eT >::memptr(), and Col< eT >::set_size().

00682   {
00683   arma_extra_debug_sigprint();
00684 
00685   typedef typename std::complex<T> eT;
00686   
00687   #if defined(ARMA_USE_LAPACK)
00688     {
00689     arma_debug_check( (A.n_rows != A.n_cols), "auxlib::eig_sym(): given matrix is not hermitian");
00690 
00691     char jobz  = 'N'; 
00692     char uplo  = 'U';
00693 
00694     int n_rows = A.n_rows;
00695     int lda    = A.n_rows;
00696     int lwork  = (std::max)(1, 2*n_rows - 1);  // TODO: automatically find best size of lwork
00697 
00698     eigval.set_size(n_rows);
00699 
00700     podarray<eT> work(lwork);
00701     podarray<T>  rwork( (std::max)(1, 3*n_rows - 2) );
00702   
00703     Mat<eT> A_copy = A;
00704     int info;
00705   
00706     arma_extra_debug_print("lapack::heev_()");
00707     lapack::heev_(&jobz, &uplo, &n_rows, A_copy.memptr(), &lda, eigval.memptr(), work.memptr(), &lwork, rwork.memptr(), &info);
00708     }
00709   #else
00710     {
00711     arma_stop("auxlib::eig_sym(): need LAPACK library");
00712     }
00713   #endif
00714   }

template<typename eT >
void auxlib::eig_sym ( Col< eT > &  eigval,
Mat< eT > &  eigvec,
const Mat< eT > &  A 
) [inline, static, inherited]

immediate eigenvalues and eigenvectors of a symmetric real matrix using LAPACK

Definition at line 722 of file auxlib_meat.hpp.

References arma_stop(), max(), podarray< T1 >::memptr(), Mat< eT >::memptr(), Mat< eT >::n_cols, Mat< eT >::n_rows, Col< eT >::set_size(), and lapack::syev_().

00723   {
00724   arma_extra_debug_sigprint();
00725   
00726   // TODO: check for aliasing
00727   
00728   #if defined(ARMA_USE_LAPACK)
00729     {
00730     const unwrap_check< Mat<eT> > tmp1(A_orig, eigval);
00731     const Mat<eT>& A_tmp = tmp1.M;
00732     
00733     const unwrap_check< Mat<eT> > tmp2(A_tmp, eigvec);
00734     const Mat<eT>& A = tmp2.M;
00735     
00736     arma_debug_check( (A.n_rows != A.n_cols), "auxlib::eig_sym(): given matrix is not square" );
00737     
00738     // rudimentary "better-than-nothing" test for symmetry
00739     //arma_debug_check( (A.at(A.n_rows-1, 0) != A.at(0, A.n_cols-1)), "auxlib::eig(): given matrix is not symmetric" );
00740     
00741     
00742     char jobz  = 'V';
00743     char uplo  = 'U';
00744     
00745     int n_rows = A.n_rows;
00746     int lwork  = (std::max)(1, 3*n_rows-1);
00747     
00748     eigval.set_size(n_rows);
00749     podarray<eT> work(lwork);
00750   
00751     eigvec = A;
00752     int info;
00753     
00754     arma_extra_debug_print("lapack::syev_()");
00755     lapack::syev_(&jobz, &uplo, &n_rows, eigvec.memptr(), &n_rows, eigval.memptr(), work.memptr(), &lwork, &info);
00756     }
00757   #else
00758     {
00759     arma_stop("auxlib::eig_sym(): need LAPACK library");
00760     }
00761   #endif
00762   
00763   }

template<typename T >
void auxlib::eig_sym ( Col< T > &  eigval,
Mat< std::complex< T > > &  eigvec,
const Mat< std::complex< T > > &  A 
) [inline, static, inherited]

immediate eigenvalues and eigenvectors of a hermitian complex matrix using LAPACK

Definition at line 771 of file auxlib_meat.hpp.

References arma_stop(), lapack::heev_(), max(), podarray< T1 >::memptr(), Mat< eT >::memptr(), Mat< eT >::n_cols, Mat< eT >::n_rows, and Col< eT >::set_size().

00772   {
00773   arma_extra_debug_sigprint();
00774   
00775   typedef typename std::complex<T> eT;
00776 
00777   #if defined(ARMA_USE_LAPACK)
00778     {
00779     const unwrap_check< Mat<eT> > tmp(A_orig, eigvec);
00780     const Mat<eT>& A = tmp.M;
00781     
00782     arma_debug_check( (A.n_rows != A.n_cols), "auxlib::eig_sym(): given matrix is not hermitian" );
00783     
00784     char jobz  = 'V';
00785     char uplo  = 'U';
00786 
00787     int n_rows = A.n_rows;
00788     int lda    = A.n_rows;
00789     int lwork  = (std::max)(1, 2*n_rows - 1);  // TODO: automatically find best size of lwork
00790     
00791     eigval.set_size(n_rows);
00792 
00793     podarray<eT> work(lwork);
00794     podarray<T>  rwork( (std::max)(1, 3*n_rows - 2) );
00795   
00796     eigvec = A;
00797     int info;
00798   
00799     arma_extra_debug_print("lapack::heev_()");
00800     lapack::heev_(&jobz, &uplo, &n_rows, eigvec.memptr(), &lda, eigval.memptr(), work.memptr(), &lwork, rwork.memptr(), &info);
00801     }
00802   #else
00803     {
00804     arma_stop("auxlib::eig_sym(): need LAPACK library");
00805     }
00806   #endif
00807   
00808   }

template<typename T >
void auxlib::eig_gen ( Col< std::complex< T > > &  eigval,
Mat< T > &  l_eigvec,
Mat< T > &  r_eigvec,
const Mat< T > &  A,
const char  side 
) [inline, inherited]

Eigenvalues and eigenvectors of a general square real matrix using LAPACK. The argument 'side' specifies which eigenvectors should be calculated (see code for mode details).

Definition at line 819 of file auxlib_meat.hpp.

References arma_stop(), lapack::geev_(), max(), podarray< T1 >::memptr(), Mat< eT >::memptr(), Mat< eT >::n_cols, Mat< eT >::n_rows, and Mat< eT >::set_size().

00826   {
00827   arma_extra_debug_sigprint();
00828 
00829   // TODO: check for aliasing
00830   
00831   #if defined(ARMA_USE_LAPACK)
00832     {
00833     arma_debug_check( (A.n_rows != A.n_cols), "auxlib::eig_gen(): given matrix is not square" );
00834     
00835     char jobvl;
00836     char jobvr;
00837 
00838     switch(side)
00839       {
00840       case 'l':  // left
00841         jobvl = 'V';
00842         jobvr = 'N';
00843         break;
00844         
00845       case 'r':  // right
00846         jobvl = 'N';
00847         jobvr = 'V';
00848         break;
00849         
00850       case 'b':  // both
00851         jobvl = 'V';
00852         jobvr = 'V';
00853         break;
00854         
00855       case 'n':  // neither
00856         jobvl = 'N';
00857         jobvr = 'N';
00858         break;
00859 
00860       default:
00861         arma_stop("auxlib::eig_gen(): parameter 'side' is invalid");
00862       }
00863 
00864        
00865     int n_rows = A.n_rows;
00866     int lda    = A.n_rows;
00867     int lwork  = (std::max)(1, 4*n_rows);  // TODO: automatically find best size of lwork
00868     
00869     eigval.set_size(n_rows);
00870     l_eigvec.set_size(n_rows, n_rows);
00871     r_eigvec.set_size(n_rows, n_rows);
00872     
00873     podarray<T> work(lwork);
00874     podarray<T> rwork( (std::max)(1, 3*n_rows) );
00875     
00876     podarray<T> wr(n_rows);
00877     podarray<T> wi(n_rows);
00878     
00879     Mat<T> A_copy = A;
00880     int info;
00881     
00882     arma_extra_debug_print("lapack::cx_geev_()");
00883     lapack::geev_(&jobvl, &jobvr, &n_rows, A_copy.memptr(), &lda, wr.memptr(), wi.memptr(), l_eigvec.memptr(), &n_rows, r_eigvec.memptr(), &n_rows, work.memptr(), &lwork, &info);
00884     
00885     
00886     eigval.set_size(n_rows);
00887     for(u32 i=0; i<u32(n_rows); ++i)
00888       {
00889       eigval[i] = std::complex<T>(wr[i], wi[i]);
00890       }
00891     
00892     }
00893   #else
00894     {
00895     arma_stop("auxlib::eig_gen(): need LAPACK library");
00896     }
00897   #endif
00898   
00899   }

template<typename T >
void auxlib::eig_gen ( Col< std::complex< T > > &  eigval,
Mat< std::complex< T > > &  l_eigvec,
Mat< std::complex< T > > &  r_eigvec,
const Mat< std::complex< T > > &  A,
const char  side 
) [inline, static, inherited]

Eigenvalues and eigenvectors of a general square complex matrix using LAPACK The argument 'side' specifies which eigenvectors should be calculated (see code for mode details).

Definition at line 912 of file auxlib_meat.hpp.

References arma_stop(), lapack::cx_geev_(), max(), podarray< T1 >::memptr(), Mat< eT >::memptr(), Mat< eT >::n_cols, Mat< eT >::n_rows, and Mat< eT >::set_size().

00919   {
00920   arma_extra_debug_sigprint();
00921 
00922   // TODO: check for aliasing
00923   
00924   typedef typename std::complex<T> eT;
00925 
00926   #if defined(ARMA_USE_LAPACK)
00927     {
00928     arma_debug_check( (A.n_rows != A.n_cols), "auxlib::eig_gen(): given matrix is not square" );
00929     
00930     char jobvl;
00931     char jobvr;
00932 
00933     switch(side)
00934       {
00935       case 'l':  // left
00936         jobvl = 'V';
00937         jobvr = 'N';
00938         break;
00939         
00940       case 'r':  // right
00941         jobvl = 'N';
00942         jobvr = 'V';
00943         break;
00944         
00945       case 'b':  // both
00946         jobvl = 'V';
00947         jobvr = 'V';
00948         break;
00949         
00950       case 'n':  // neither
00951         jobvl = 'N';
00952         jobvr = 'N';
00953         break;
00954 
00955       default:
00956         arma_stop("auxlib::eig_gen(): parameter 'side' is invalid");
00957       }
00958     
00959        
00960     int n_rows = A.n_rows;
00961     int lda    = A.n_rows;
00962     int lwork  = (std::max)(1, 4*n_rows);  // TODO: automatically find best size of lwork
00963     
00964     eigval.set_size(n_rows);
00965     l_eigvec.set_size(n_rows, n_rows);
00966     r_eigvec.set_size(n_rows, n_rows);
00967     
00968     podarray<eT> work(lwork);
00969     podarray<T>  rwork( (std::max)(1, 3*n_rows) );  // was 2,3
00970     
00971     Mat<eT> A_copy = A;
00972     int info;
00973     
00974     arma_extra_debug_print("lapack::cx_geev_()");
00975     lapack::cx_geev_(&jobvl, &jobvr, &n_rows, A_copy.memptr(), &lda, eigval.memptr(), l_eigvec.memptr(), &n_rows, r_eigvec.memptr(), &n_rows, work.memptr(), &lwork, rwork.memptr(), &info);
00976     }
00977   #else
00978     {
00979     arma_stop("auxlib::eig_gen(): need LAPACK library");
00980     }
00981   #endif
00982   
00983   }

template<typename eT >
bool auxlib::chol ( Mat< eT > &  out,
const Mat< eT > &  X 
) [inline, static, inherited]

Definition at line 990 of file auxlib_meat.hpp.

References arma_stop(), Mat< eT >::colptr(), Mat< eT >::memptr(), Mat< eT >::n_rows, and lapack::potrf_().

Referenced by chol().

00991   {
00992   arma_extra_debug_sigprint();
00993   
00994   // TODO: check for aliasing
00995   
00996   
00997   #if defined(ARMA_USE_LAPACK)
00998     {
00999     char uplo = 'U';
01000     int  n    = X.n_rows;
01001     int  lda  = X.n_rows;
01002     int  info;
01003     
01004     out = X;
01005     lapack::potrf_(&uplo, &n, out.memptr(), &lda, &info);
01006     
01007     for(u32 col=0; col<X.n_rows; ++col)
01008       {
01009       eT* colptr = out.colptr(col);
01010       for(u32 row=col+1; row<X.n_rows; ++row)
01011         {
01012         colptr[row] = eT(0);
01013         }
01014       }
01015     
01016     return (info == 0);
01017     }
01018   #else
01019     {
01020     arma_stop("auxlib::chol(): need LAPACK library");
01021     return false;
01022     }
01023   #endif
01024   }

template<typename eT >
bool auxlib::qr ( Mat< eT > &  Q,
Mat< eT > &  R,
const Mat< eT > &  X 
) [inline, static, inherited]

Definition at line 1031 of file auxlib_meat.hpp.

References arma_stop(), Mat< eT >::at(), lapack::geqrf_(), max(), Mat< eT >::mem, podarray< T1 >::memptr(), Mat< eT >::memptr(), min(), Mat< eT >::n_cols, Mat< eT >::n_elem, Mat< eT >::n_rows, lapack::orgqr_(), Mat< eT >::set_size(), podarray< T1 >::set_size(), tmp_real(), and lapack::ungqr_().

Referenced by qr().

01032   {
01033   arma_extra_debug_sigprint();
01034   
01035   // TODO: check for aliasing
01036   
01037   #if defined(ARMA_USE_LAPACK)
01038     {
01039     int m            = static_cast<int>(X.n_rows);
01040     int n            = static_cast<int>(X.n_cols);
01041     int work_len     = (std::max)(1,n);
01042     int work_len_tmp;
01043     int k            = (std::min)(m,n);
01044     int info;
01045     
01046     podarray<eT> tau(k);
01047     podarray<eT> work(work_len);
01048     
01049     R = X;
01050     
01051     // query for the optimum value of work_len
01052     work_len_tmp = -1;
01053     lapack::geqrf_(&m, &n, R.memptr(), &m, tau.memptr(), work.memptr(), &work_len_tmp, &info);
01054     
01055     if(info == 0)
01056       {
01057       work_len = static_cast<int>(auxlib::tmp_real(work[0]));
01058       work.set_size(work_len);
01059       }
01060     
01061     lapack::geqrf_(&m, &n, R.memptr(), &m, tau.memptr(), work.memptr(), &work_len, &info);
01062     
01063     Q.set_size(X.n_rows,X.n_rows);
01064     
01065           eT* Q_mem = Q.memptr();
01066     const eT* R_mem = R.mem;
01067     
01068     const u32 n_elem_copy = (std::min)(Q.n_elem, R.n_elem);
01069     for(u32 i=0; i < n_elem_copy; ++i)
01070       {
01071       Q_mem[i] = R_mem[i];
01072       }
01073     
01074     
01075     // construct R
01076     for(u32 row=0; row < R.n_rows; ++row)
01077       {
01078       const u32 n_elem_tmp = (std::min)(row, R.n_cols);
01079       for(u32 col=0; col < n_elem_tmp; ++col)
01080         {
01081         R.at(row,col) = eT(0);
01082         }
01083       }
01084       
01085     
01086     if( (is_float<eT>::value == true) || (is_double<eT>::value == true) )
01087       {
01088       // query for the optimum value of work_len
01089       work_len_tmp = -1;
01090       lapack::orgqr_(&m, &m, &k, Q.memptr(), &m, tau.memptr(), work.memptr(), &work_len_tmp, &info);
01091       
01092       if(info == 0)
01093         {
01094         work_len = static_cast<int>(auxlib::tmp_real(work[0]));
01095         work.set_size(work_len);
01096         }
01097       
01098       lapack::orgqr_(&m, &m, &k, Q.memptr(), &m, tau.memptr(), work.memptr(), &work_len, &info);
01099       }
01100     else
01101     if( (is_supported_complex_float<eT>::value == true) || (is_supported_complex_double<eT>::value == true) )
01102       {
01103       // query for the optimum value of work_len
01104       work_len_tmp = -1;
01105       lapack::ungqr_(&m, &m, &k, Q.memptr(), &m, tau.memptr(), work.memptr(), &work_len_tmp, &info);
01106       
01107       if(info == 0)
01108         {
01109         work_len = static_cast<int>(auxlib::tmp_real(work[0]));
01110         work.set_size(work_len);
01111         }
01112       
01113       lapack::ungqr_(&m, &m, &k, Q.memptr(), &m, tau.memptr(), work.memptr(), &work_len, &info);
01114       }
01115     
01116     return (info == 0);
01117     }
01118   #else
01119     {
01120     arma_stop("auxlib::qr(): need LAPACK library");
01121     return false;
01122     }
01123   #endif
01124   }

template<typename eT >
bool auxlib::svd ( Col< eT > &  S,
const Mat< eT > &  X 
) [inline, static, inherited]

Definition at line 1131 of file auxlib_meat.hpp.

References arma_stop(), podarray< T1 >::memptr(), Mat< eT >::memptr(), min(), Mat< eT >::n_cols, Mat< eT >::n_rows, podarray< T1 >::set_size(), and Col< eT >::set_size().

Referenced by svd().

01132   {
01133   arma_extra_debug_sigprint();
01134   
01135   // TODO: check for aliasing
01136   
01137   #if defined(ARMA_USE_LAPACK)
01138     {
01139     Mat<eT> A = X;
01140     
01141     char jobu  = 'N';
01142     char jobvt = 'N';
01143     
01144     int  m     = A.n_rows;
01145     int  n     = A.n_cols;
01146     int  lda   = A.n_rows;
01147     int  ldu   = A.n_rows;
01148     int  ldvt  = A.n_cols;
01149     int  lwork = 2;
01150     int  info;
01151     
01152     Mat<eT> U(1,1);
01153     Mat<eT> V(1,1);
01154     
01155     S.set_size( (std::min)(m, n) );
01156     
01157     podarray<eT> work(lwork);
01158   
01159   
01160     // let gesvd_() calculate the optimum size of the workspace
01161     int lwork_tmp = -1;
01162     
01163     lapack::gesvd_<eT>
01164       (
01165       &jobu, &jobvt,
01166       &m,&n,
01167       A.memptr(), &lda,
01168       S.memptr(),
01169       U.memptr(), &ldu,
01170       V.memptr(), &ldvt,
01171       work.memptr(), &lwork_tmp,
01172       &info
01173       );
01174     
01175     if(info == 0)
01176       {
01177       lwork = static_cast<int>(work[0]);
01178       work.set_size(lwork);
01179       
01180       lapack::gesvd_<eT>
01181         (
01182         &jobu, &jobvt,
01183         &m, &n,
01184         A.memptr(), &lda,
01185         S.memptr(),
01186         U.memptr(), &ldu,
01187         V.memptr(), &ldvt,
01188         work.memptr(), &lwork,
01189         &info
01190         );
01191       
01192       return (info == 0);
01193       }
01194     else
01195       {
01196       return false;
01197       }
01198     }
01199   #else
01200     {
01201     arma_stop("auxlib::svd(): need LAPACK library");
01202     return false;
01203     }
01204   #endif
01205   }

template<typename T >
bool auxlib::svd ( Col< T > &  S,
const Mat< std::complex< T > > &  X 
) [inline, static, inherited]

Definition at line 1212 of file auxlib_meat.hpp.

References arma_stop(), max(), podarray< T1 >::memptr(), Mat< eT >::memptr(), min(), Mat< eT >::n_cols, Mat< eT >::n_rows, real(), podarray< T1 >::set_size(), and Col< eT >::set_size().

01213   {
01214   arma_extra_debug_sigprint();
01215   
01216   typedef std::complex<T> eT;
01217   
01218   // TODO: check for aliasing
01219   
01220   #if defined(ARMA_USE_LAPACK)
01221     {
01222     Mat<eT> A = X;
01223     
01224     char jobu  = 'N';
01225     char jobvt = 'N';
01226     
01227     int  m     = A.n_rows;
01228     int  n     = A.n_cols;
01229     int  lda   = A.n_rows;
01230     int  ldu   = A.n_rows;
01231     int  ldvt  = A.n_cols;
01232     int  lwork = 2 * (std::min)(m,n) + (std::max)(m,n);
01233     int  info;
01234     
01235     Mat<eT> U(1,1);
01236     Mat<eT> V(1,1);
01237     
01238     S.set_size( (std::min)(m,n) );
01239     
01240     podarray<eT> work(lwork);
01241     podarray<T>  rwork( 5*(std::min)(m,n) );
01242   
01243     // let gesvd_() calculate the optimum size of the workspace
01244     int lwork_tmp = -1;
01245     
01246     lapack::cx_gesvd_<T>
01247       (
01248       &jobu, &jobvt,
01249       &m, &n,
01250       A.memptr(), &lda,
01251       S.memptr(),
01252       U.memptr(), &ldu,
01253       V.memptr(), &ldvt,
01254       work.memptr(), &lwork_tmp,
01255       rwork.memptr(),
01256       &info
01257       );
01258     
01259     if(info == 0)
01260       {
01261       int proposed_lwork = static_cast<int>(real(work[0]));
01262       if(proposed_lwork > lwork)
01263         {
01264         lwork = proposed_lwork;
01265         work.set_size(lwork);
01266         }
01267       
01268       lapack::cx_gesvd_<T>
01269         (
01270         &jobu, &jobvt,
01271         &m, &n,
01272         A.memptr(), &lda,
01273         S.memptr(),
01274         U.memptr(), &ldu,
01275         V.memptr(), &ldvt,
01276         work.memptr(), &lwork,
01277         rwork.memptr(),
01278         &info
01279         );
01280       
01281       return (info == 0);
01282       }
01283     else
01284       {
01285       return false;
01286       }
01287     }
01288   #else
01289     {
01290     arma_stop("auxlib::svd(): need LAPACK library");
01291     return false;
01292     }
01293   #endif
01294   }

template<typename eT >
bool auxlib::svd ( Mat< eT > &  U,
Col< eT > &  S,
Mat< eT > &  V,
const Mat< eT > &  X 
) [inline, static, inherited]

Definition at line 1301 of file auxlib_meat.hpp.

References op_trans::apply(), arma_stop(), podarray< T1 >::memptr(), Mat< eT >::memptr(), min(), Mat< eT >::n_cols, Mat< eT >::n_rows, podarray< T1 >::set_size(), Col< eT >::set_size(), and Mat< eT >::set_size().

01302   {
01303   arma_extra_debug_sigprint();
01304   
01305   // TODO: check for aliasing
01306 
01307   #if defined(ARMA_USE_LAPACK)
01308     {
01309     Mat<eT> A = X;
01310     
01311     char jobu  = 'A';
01312     char jobvt = 'A';
01313     
01314     int  m     = A.n_rows;
01315     int  n     = A.n_cols;
01316     int  lda   = A.n_rows;
01317     int  ldu   = A.n_rows;
01318     int  ldvt  = A.n_cols;
01319     int  lwork = 2;
01320     int  info;
01321     
01322     U.set_size(m,m);
01323     V.set_size(n,n);
01324     
01325     S.set_size( (std::min)(m,n) );
01326     podarray<eT> work(lwork);
01327   
01328     // let gesvd_() calculate the optimum size of the workspace
01329     int lwork_tmp = -1;
01330     
01331     lapack::gesvd_<eT>
01332       (
01333       &jobu, &jobvt,
01334       &m, &n,
01335       A.memptr(), &lda,
01336       S.memptr(),
01337       U.memptr(), &ldu,
01338       V.memptr(), &ldvt,
01339       work.memptr(), &lwork_tmp,
01340       &info
01341       );
01342     
01343     if(info == 0)
01344       {
01345       lwork = static_cast<int>(work[0]);
01346       work.set_size(lwork);
01347       
01348       lapack::gesvd_<eT>
01349         (
01350         &jobu, &jobvt,
01351         &m, &n,
01352         A.memptr(), &lda,
01353         S.memptr(),
01354         U.memptr(), &ldu,
01355         V.memptr(), &ldvt,
01356         work.memptr(), &lwork,
01357         &info
01358         );
01359       
01360       op_trans::apply(V,V);  // op_trans will work out that an in-place transpose can be done
01361       
01362       return (info == 0);
01363       }
01364     else
01365       {
01366       return false;
01367       }
01368     }
01369   #else
01370     {
01371     arma_stop("auxlib::svd(): need LAPACK library");
01372     return false;
01373     }
01374   #endif
01375   }

template<typename T >
bool auxlib::svd ( Mat< std::complex< T > > &  U,
Col< T > &  S,
Mat< std::complex< T > > &  V,
const Mat< std::complex< T > > &  X 
) [inline, static, inherited]

Definition at line 1382 of file auxlib_meat.hpp.

References op_htrans::apply(), arma_stop(), podarray< T1 >::memptr(), Mat< eT >::memptr(), min(), Mat< eT >::n_cols, Mat< eT >::n_rows, real(), podarray< T1 >::set_size(), and Col< eT >::set_size().

01383   {
01384   arma_extra_debug_sigprint();
01385   
01386   typedef std::complex<T> eT;
01387   
01388   // TODO: check for aliasing
01389   
01390   #if defined(ARMA_USE_LAPACK)
01391     {
01392     Mat<eT> A = X;
01393     
01394     char jobu  = 'A';
01395     char jobvt = 'A';
01396     
01397     int  m     = A.n_rows;
01398     int  n     = A.n_cols;
01399     int  lda   = A.n_rows;
01400     int  ldu   = A.n_rows;
01401     int  ldvt  = A.n_cols;
01402     int  lwork = 2;
01403     int  info;
01404     
01405     U.set_size(m,m);
01406     V.set_size(n,n);
01407     
01408     S.set_size( (std::min)(m,n) );
01409     
01410     podarray<eT> work(lwork);
01411     podarray<T>  rwork( 5*(std::min)(m,n) );
01412   
01413     // let gesvd_() calculate the optimum size of the workspace
01414     int lwork_tmp = -1;
01415     lapack::cx_gesvd_<T>
01416      (
01417      &jobu, &jobvt,
01418      &m, &n,
01419      A.memptr(), &lda,
01420      S.memptr(),
01421      U.memptr(), &ldu,
01422      V.memptr(), &ldvt,
01423      work.memptr(), &lwork_tmp,
01424      rwork.memptr(),
01425      &info
01426      );
01427     
01428     if(info == 0)
01429       {
01430       lwork = static_cast<int>(real(work[0]));
01431       work.set_size(lwork);
01432       
01433       lapack::cx_gesvd_<T>
01434         (
01435         &jobu, &jobvt,
01436         &m, &n,
01437         A.memptr(), &lda,
01438         S.memptr(),
01439         U.memptr(), &ldu,
01440         V.memptr(), &ldvt,
01441         work.memptr(), &lwork,
01442         rwork.memptr(),
01443         &info
01444         );
01445       
01446       op_htrans::apply(V,V);  // op_htrans will work out that an in-place transpose can be done
01447       
01448       return (info == 0);
01449       }
01450     else
01451       {
01452       return false;
01453       }
01454     }
01455   #else
01456     {
01457     arma_stop("auxlib::svd(): need LAPACK library");
01458     return false;
01459     }
01460   #endif
01461   
01462   }

template<typename eT >
bool auxlib::solve ( Mat< eT > &  out,
const Mat< eT > &  A,
const Mat< eT > &  B 
) [inline, static, inherited]

Solve a system of linear equations Assumes that A.n_rows = A.n_cols and B.n_rows = A.n_rows.

Definition at line 1471 of file auxlib_meat.hpp.

References arma_stop(), podarray< T1 >::memptr(), Mat< eT >::memptr(), Mat< eT >::n_cols, and Mat< eT >::n_rows.

Referenced by solve().

01472   {
01473   arma_extra_debug_sigprint();
01474   
01475   #if defined(ARMA_USE_LAPACK)
01476     {
01477     int n    = A.n_rows;
01478     int lda  = A.n_rows;
01479     int ldb  = A.n_rows;
01480     int nrhs = B.n_cols;
01481     int info;
01482     
01483     podarray<int> ipiv(n);
01484     
01485     out = B;
01486     Mat<eT> A_copy = A;
01487   
01488     lapack::gesv_<eT>(&n, &nrhs, A_copy.memptr(), &lda, ipiv.memptr(), out.memptr(), &ldb, &info);
01489   
01490     return (info == 0);
01491     }
01492   #else
01493     {
01494     arma_stop("auxlib::solve(): need LAPACK library");
01495     return false;
01496     }
01497   #endif
01498   }

template<typename eT >
bool auxlib::solve_od ( Mat< eT > &  out,
const Mat< eT > &  A,
const Mat< eT > &  B 
) [inline, static, inherited]

Solve an over-determined system. Assumes that A.n_rows > A.n_cols and B.n_rows = A.n_rows.

Definition at line 1508 of file auxlib_meat.hpp.

References arma_stop(), Mat< eT >::colptr(), syslib::copy_elem(), max(), podarray< T1 >::memptr(), Mat< eT >::memptr(), Mat< eT >::n_cols, Mat< eT >::n_rows, Mat< eT >::set_size(), and trans().

Referenced by solve().

01509   {
01510   arma_extra_debug_sigprint();
01511   
01512   #if defined(ARMA_USE_LAPACK)
01513     {
01514     char trans = 'N';
01515     
01516     int  m     = A.n_rows;
01517     int  n     = A.n_cols;
01518     int  lda   = A.n_rows;
01519     int  ldb   = A.n_rows;
01520     int  nrhs  = B.n_cols;
01521     int  lwork = n + (std::max)(n, nrhs);
01522     int  info;
01523     
01524     Mat<eT> A_copy = A;
01525     Mat<eT> tmp    = B;
01526     
01527     
01528     podarray<eT> work(lwork);
01529     
01530     arma_extra_debug_print("lapack::gels_()");
01531     
01532     // NOTE:
01533     // the dgels() function in the lapack library supplied by ATLAS 3.6
01534     // seems to have problems
01535     
01536     lapack::gels_<eT>
01537       (
01538       &trans, &m, &n, &nrhs,
01539       A_copy.memptr(), &lda,
01540       tmp.memptr(), &ldb,
01541       work.memptr(), &lwork,
01542       &info
01543       );
01544     
01545     arma_extra_debug_print("lapack::gels_() -- finished");
01546     
01547     out.set_size(A.n_cols, B.n_cols);
01548     
01549     for(u32 col=0; col<B.n_cols; ++col)
01550       {
01551       syslib::copy_elem( out.colptr(col), tmp.colptr(col), A.n_cols );
01552       }
01553     
01554     return (info == 0);
01555     }
01556   #else
01557     {
01558     arma_stop("auxlib::solve_od(): need LAPACK library");
01559     return false;
01560     }
01561   #endif
01562   }

template<typename eT >
bool auxlib::solve_ud ( Mat< eT > &  out,
const Mat< eT > &  A,
const Mat< eT > &  B 
) [inline, static, inherited]

Solve an under-determined system. Assumes that A.n_rows < A.n_cols and B.n_rows = A.n_rows.

Definition at line 1572 of file auxlib_meat.hpp.

References arma_stop(), Mat< eT >::colptr(), syslib::copy_elem(), max(), podarray< T1 >::memptr(), Mat< eT >::memptr(), Mat< eT >::n_cols, Mat< eT >::n_rows, Mat< eT >::set_size(), trans(), and Mat< eT >::zeros().

Referenced by solve().

01573   {
01574   arma_extra_debug_sigprint();
01575   
01576   #if defined(ARMA_USE_LAPACK)
01577     {
01578     char trans = 'N';
01579     
01580     int  m     = A.n_rows;
01581     int  n     = A.n_cols;
01582     int  lda   = A.n_rows;
01583     int  ldb   = A.n_cols;
01584     int  nrhs  = B.n_cols;
01585     int  lwork = m + (std::max)(m,nrhs);
01586     int  info;
01587     
01588     
01589     Mat<eT> A_copy = A;
01590     
01591     Mat<eT> tmp;
01592     tmp.zeros(A.n_cols, B.n_cols);
01593     
01594     for(u32 col=0; col<B.n_cols; ++col)
01595       {
01596       eT* tmp_colmem = tmp.colptr(col);
01597       
01598       syslib::copy_elem( tmp_colmem, B.colptr(col), B.n_rows );
01599       
01600       for(u32 row=B.n_rows; row<A.n_cols; ++row)
01601         {
01602         tmp_colmem[row] = eT(0);
01603         }
01604       }
01605     
01606     podarray<eT> work(lwork);
01607     
01608     arma_extra_debug_print("lapack::gels_()");
01609     
01610     // NOTE:
01611     // the dgels() function in the lapack library supplied by ATLAS 3.6
01612     // seems to have problems
01613     
01614     lapack::gels_<eT>
01615       (
01616       &trans, &m, &n, &nrhs,
01617       A_copy.memptr(), &lda,
01618       tmp.memptr(), &ldb,
01619       work.memptr(), &lwork,
01620       &info
01621       );
01622     
01623     arma_extra_debug_print("lapack::gels_() -- finished");
01624     
01625     out.set_size(A.n_cols, B.n_cols);
01626     
01627     for(u32 col=0; col<B.n_cols; ++col)
01628       {
01629       syslib::copy_elem( out.colptr(col), tmp.colptr(col), A.n_cols );
01630       }
01631   
01632     return (info == 0);
01633     }
01634   #else
01635     {
01636     arma_stop("auxlib::solve_ud(): need LAPACK library");
01637     return false;
01638     }
01639   #endif
01640   }