max@0: // Copyright (C) 2008-2012 NICTA (www.nicta.com.au)
max@0: // Copyright (C) 2008-2012 Conrad Sanderson
max@0: // Copyright (C) 2009 Edmund Highcock
max@0: // Copyright (C) 2011 James Sanders
max@0: // Copyright (C) 2011 Stanislav Funiak
max@0: // 
max@0: // This file is part of the Armadillo C++ library.
max@0: // It is provided without any warranty of fitness
max@0: // for any purpose. You can redistribute this file
max@0: // and/or modify it under the terms of the GNU
max@0: // Lesser General Public License (LGPL) as published
max@0: // by the Free Software Foundation, either version 3
max@0: // of the License or (at your option) any later version.
max@0: // (see http://www.opensource.org/licenses for more info)
max@0: 
max@0: 
max@0: //! \addtogroup auxlib
max@0: //! @{
max@0: 
max@0: 
max@0: 
max@0: //! immediate matrix inverse
max@0: template<typename eT, typename T1>
max@0: inline
max@0: bool
max@0: auxlib::inv(Mat<eT>& out, const Base<eT,T1>& X, const bool slow)
max@0:   {
max@0:   arma_extra_debug_sigprint();
max@0:   
max@0:   out = X.get_ref();
max@0:   
max@0:   arma_debug_check( (out.is_square() == false), "inv(): given matrix is not square" );
max@0:   
max@0:   bool status = false;
max@0:   
max@0:   const uword N = out.n_rows;
max@0:   
max@0:   if( (N <= 4) && (slow == false) )
max@0:     {
max@0:     status = auxlib::inv_inplace_tinymat(out, N);
max@0:     }
max@0:     
max@0:   if( (N > 4) || (status == false) )
max@0:     {
max@0:     status = auxlib::inv_inplace_lapack(out);
max@0:     }
max@0:   
max@0:   return status;
max@0:   }
max@0: 
max@0: 
max@0: 
max@0: template<typename eT>
max@0: inline
max@0: bool
max@0: auxlib::inv(Mat<eT>& out, const Mat<eT>& X, const bool slow)
max@0:   {
max@0:   arma_extra_debug_sigprint();
max@0:   
max@0:   arma_debug_check( (X.is_square() == false), "inv(): given matrix is not square" );
max@0:   
max@0:   bool status = false;
max@0:   
max@0:   const uword N = X.n_rows;
max@0:   
max@0:   if( (N <= 4) && (slow == false) )
max@0:     {
max@0:     status = (&out != &X) ? auxlib::inv_noalias_tinymat(out, X, N) : auxlib::inv_inplace_tinymat(out, N);
max@0:     }
max@0:   
max@0:   if( (N > 4) || (status == false) )
max@0:     {
max@0:     out = X;
max@0:     status = auxlib::inv_inplace_lapack(out);
max@0:     }
max@0:   
max@0:   return status;
max@0:   }
max@0: 
max@0: 
max@0: 
max@0: template<typename eT>
max@0: inline
max@0: bool
max@0: auxlib::inv_noalias_tinymat(Mat<eT>& out, const Mat<eT>& X, const uword N)
max@0:   {
max@0:   arma_extra_debug_sigprint();
max@0:   
max@0:   bool det_ok = true;
max@0:   
max@0:   out.set_size(N,N);
max@0:   
max@0:   switch(N)
max@0:     {
max@0:     case 1:
max@0:       {
max@0:       out[0] = eT(1) / X[0];
max@0:       };
max@0:       break;
max@0:       
max@0:     case 2:
max@0:       {
max@0:       const eT* Xm = X.memptr();
max@0:       
max@0:       const eT a = Xm[pos<0,0>::n2];
max@0:       const eT b = Xm[pos<0,1>::n2];
max@0:       const eT c = Xm[pos<1,0>::n2];
max@0:       const eT d = Xm[pos<1,1>::n2];
max@0:       
max@0:       const eT tmp_det = (a*d - b*c);
max@0:       
max@0:       if(tmp_det != eT(0))
max@0:         {
max@0:         eT* outm = out.memptr();
max@0:         
max@0:         outm[pos<0,0>::n2] =  d / tmp_det;
max@0:         outm[pos<0,1>::n2] = -b / tmp_det;
max@0:         outm[pos<1,0>::n2] = -c / tmp_det;
max@0:         outm[pos<1,1>::n2] =  a / tmp_det;
max@0:         }
max@0:       else
max@0:         {
max@0:         det_ok = false;
max@0:         }
max@0:       };
max@0:       break;
max@0:     
max@0:     case 3:
max@0:       {
max@0:       const eT* X_col0 = X.colptr(0);
max@0:       const eT a11 = X_col0[0];
max@0:       const eT a21 = X_col0[1];
max@0:       const eT a31 = X_col0[2];
max@0:       
max@0:       const eT* X_col1 = X.colptr(1);
max@0:       const eT a12 = X_col1[0];
max@0:       const eT a22 = X_col1[1];
max@0:       const eT a32 = X_col1[2];
max@0:       
max@0:       const eT* X_col2 = X.colptr(2);
max@0:       const eT a13 = X_col2[0];
max@0:       const eT a23 = X_col2[1];
max@0:       const eT a33 = X_col2[2];
max@0:       
max@0:       const eT tmp_det = a11*(a33*a22 - a32*a23) - a21*(a33*a12-a32*a13) + a31*(a23*a12 - a22*a13);
max@0:       
max@0:       if(tmp_det != eT(0))
max@0:         {
max@0:         eT* out_col0 = out.colptr(0);
max@0:         out_col0[0] =  (a33*a22 - a32*a23) / tmp_det;
max@0:         out_col0[1] = -(a33*a21 - a31*a23) / tmp_det;
max@0:         out_col0[2] =  (a32*a21 - a31*a22) / tmp_det;
max@0:         
max@0:         eT* out_col1 = out.colptr(1);
max@0:         out_col1[0] = -(a33*a12 - a32*a13) / tmp_det;
max@0:         out_col1[1] =  (a33*a11 - a31*a13) / tmp_det;
max@0:         out_col1[2] = -(a32*a11 - a31*a12) / tmp_det;
max@0:         
max@0:         eT* out_col2 = out.colptr(2);
max@0:         out_col2[0] =  (a23*a12 - a22*a13) / tmp_det;
max@0:         out_col2[1] = -(a23*a11 - a21*a13) / tmp_det;
max@0:         out_col2[2] =  (a22*a11 - a21*a12) / tmp_det;
max@0:         }
max@0:       else
max@0:         {
max@0:         det_ok = false;
max@0:         }
max@0:       };
max@0:       break;
max@0:       
max@0:     case 4:
max@0:       {
max@0:       const eT tmp_det = det(X);
max@0:       
max@0:       if(tmp_det != eT(0))
max@0:         {
max@0:         const eT* Xm   = X.memptr();
max@0:               eT* outm = out.memptr();
max@0:         
max@0:         outm[pos<0,0>::n4] = ( Xm[pos<1,2>::n4]*Xm[pos<2,3>::n4]*Xm[pos<3,1>::n4] - Xm[pos<1,3>::n4]*Xm[pos<2,2>::n4]*Xm[pos<3,1>::n4] + Xm[pos<1,3>::n4]*Xm[pos<2,1>::n4]*Xm[pos<3,2>::n4] - Xm[pos<1,1>::n4]*Xm[pos<2,3>::n4]*Xm[pos<3,2>::n4] - Xm[pos<1,2>::n4]*Xm[pos<2,1>::n4]*Xm[pos<3,3>::n4] + Xm[pos<1,1>::n4]*Xm[pos<2,2>::n4]*Xm[pos<3,3>::n4] ) / tmp_det;
max@0:         outm[pos<1,0>::n4] = ( Xm[pos<1,3>::n4]*Xm[pos<2,2>::n4]*Xm[pos<3,0>::n4] - Xm[pos<1,2>::n4]*Xm[pos<2,3>::n4]*Xm[pos<3,0>::n4] - Xm[pos<1,3>::n4]*Xm[pos<2,0>::n4]*Xm[pos<3,2>::n4] + Xm[pos<1,0>::n4]*Xm[pos<2,3>::n4]*Xm[pos<3,2>::n4] + Xm[pos<1,2>::n4]*Xm[pos<2,0>::n4]*Xm[pos<3,3>::n4] - Xm[pos<1,0>::n4]*Xm[pos<2,2>::n4]*Xm[pos<3,3>::n4] ) / tmp_det;
max@0:         outm[pos<2,0>::n4] = ( Xm[pos<1,1>::n4]*Xm[pos<2,3>::n4]*Xm[pos<3,0>::n4] - Xm[pos<1,3>::n4]*Xm[pos<2,1>::n4]*Xm[pos<3,0>::n4] + Xm[pos<1,3>::n4]*Xm[pos<2,0>::n4]*Xm[pos<3,1>::n4] - Xm[pos<1,0>::n4]*Xm[pos<2,3>::n4]*Xm[pos<3,1>::n4] - Xm[pos<1,1>::n4]*Xm[pos<2,0>::n4]*Xm[pos<3,3>::n4] + Xm[pos<1,0>::n4]*Xm[pos<2,1>::n4]*Xm[pos<3,3>::n4] ) / tmp_det;
max@0:         outm[pos<3,0>::n4] = ( Xm[pos<1,2>::n4]*Xm[pos<2,1>::n4]*Xm[pos<3,0>::n4] - Xm[pos<1,1>::n4]*Xm[pos<2,2>::n4]*Xm[pos<3,0>::n4] - Xm[pos<1,2>::n4]*Xm[pos<2,0>::n4]*Xm[pos<3,1>::n4] + Xm[pos<1,0>::n4]*Xm[pos<2,2>::n4]*Xm[pos<3,1>::n4] + Xm[pos<1,1>::n4]*Xm[pos<2,0>::n4]*Xm[pos<3,2>::n4] - Xm[pos<1,0>::n4]*Xm[pos<2,1>::n4]*Xm[pos<3,2>::n4] ) / tmp_det;
max@0:         
max@0:         outm[pos<0,1>::n4] = ( Xm[pos<0,3>::n4]*Xm[pos<2,2>::n4]*Xm[pos<3,1>::n4] - Xm[pos<0,2>::n4]*Xm[pos<2,3>::n4]*Xm[pos<3,1>::n4] - Xm[pos<0,3>::n4]*Xm[pos<2,1>::n4]*Xm[pos<3,2>::n4] + Xm[pos<0,1>::n4]*Xm[pos<2,3>::n4]*Xm[pos<3,2>::n4] + Xm[pos<0,2>::n4]*Xm[pos<2,1>::n4]*Xm[pos<3,3>::n4] - Xm[pos<0,1>::n4]*Xm[pos<2,2>::n4]*Xm[pos<3,3>::n4] ) / tmp_det;
max@0:         outm[pos<1,1>::n4] = ( Xm[pos<0,2>::n4]*Xm[pos<2,3>::n4]*Xm[pos<3,0>::n4] - Xm[pos<0,3>::n4]*Xm[pos<2,2>::n4]*Xm[pos<3,0>::n4] + Xm[pos<0,3>::n4]*Xm[pos<2,0>::n4]*Xm[pos<3,2>::n4] - Xm[pos<0,0>::n4]*Xm[pos<2,3>::n4]*Xm[pos<3,2>::n4] - Xm[pos<0,2>::n4]*Xm[pos<2,0>::n4]*Xm[pos<3,3>::n4] + Xm[pos<0,0>::n4]*Xm[pos<2,2>::n4]*Xm[pos<3,3>::n4] ) / tmp_det;
max@0:         outm[pos<2,1>::n4] = ( Xm[pos<0,3>::n4]*Xm[pos<2,1>::n4]*Xm[pos<3,0>::n4] - Xm[pos<0,1>::n4]*Xm[pos<2,3>::n4]*Xm[pos<3,0>::n4] - Xm[pos<0,3>::n4]*Xm[pos<2,0>::n4]*Xm[pos<3,1>::n4] + Xm[pos<0,0>::n4]*Xm[pos<2,3>::n4]*Xm[pos<3,1>::n4] + Xm[pos<0,1>::n4]*Xm[pos<2,0>::n4]*Xm[pos<3,3>::n4] - Xm[pos<0,0>::n4]*Xm[pos<2,1>::n4]*Xm[pos<3,3>::n4] ) / tmp_det;
max@0:         outm[pos<3,1>::n4] = ( Xm[pos<0,1>::n4]*Xm[pos<2,2>::n4]*Xm[pos<3,0>::n4] - Xm[pos<0,2>::n4]*Xm[pos<2,1>::n4]*Xm[pos<3,0>::n4] + Xm[pos<0,2>::n4]*Xm[pos<2,0>::n4]*Xm[pos<3,1>::n4] - Xm[pos<0,0>::n4]*Xm[pos<2,2>::n4]*Xm[pos<3,1>::n4] - Xm[pos<0,1>::n4]*Xm[pos<2,0>::n4]*Xm[pos<3,2>::n4] + Xm[pos<0,0>::n4]*Xm[pos<2,1>::n4]*Xm[pos<3,2>::n4] ) / tmp_det;
max@0:         
max@0:         outm[pos<0,2>::n4] = ( Xm[pos<0,2>::n4]*Xm[pos<1,3>::n4]*Xm[pos<3,1>::n4] - Xm[pos<0,3>::n4]*Xm[pos<1,2>::n4]*Xm[pos<3,1>::n4] + Xm[pos<0,3>::n4]*Xm[pos<1,1>::n4]*Xm[pos<3,2>::n4] - Xm[pos<0,1>::n4]*Xm[pos<1,3>::n4]*Xm[pos<3,2>::n4] - Xm[pos<0,2>::n4]*Xm[pos<1,1>::n4]*Xm[pos<3,3>::n4] + Xm[pos<0,1>::n4]*Xm[pos<1,2>::n4]*Xm[pos<3,3>::n4] ) / tmp_det;
max@0:         outm[pos<1,2>::n4] = ( Xm[pos<0,3>::n4]*Xm[pos<1,2>::n4]*Xm[pos<3,0>::n4] - Xm[pos<0,2>::n4]*Xm[pos<1,3>::n4]*Xm[pos<3,0>::n4] - Xm[pos<0,3>::n4]*Xm[pos<1,0>::n4]*Xm[pos<3,2>::n4] + Xm[pos<0,0>::n4]*Xm[pos<1,3>::n4]*Xm[pos<3,2>::n4] + Xm[pos<0,2>::n4]*Xm[pos<1,0>::n4]*Xm[pos<3,3>::n4] - Xm[pos<0,0>::n4]*Xm[pos<1,2>::n4]*Xm[pos<3,3>::n4] ) / tmp_det;
max@0:         outm[pos<2,2>::n4] = ( Xm[pos<0,1>::n4]*Xm[pos<1,3>::n4]*Xm[pos<3,0>::n4] - Xm[pos<0,3>::n4]*Xm[pos<1,1>::n4]*Xm[pos<3,0>::n4] + Xm[pos<0,3>::n4]*Xm[pos<1,0>::n4]*Xm[pos<3,1>::n4] - Xm[pos<0,0>::n4]*Xm[pos<1,3>::n4]*Xm[pos<3,1>::n4] - Xm[pos<0,1>::n4]*Xm[pos<1,0>::n4]*Xm[pos<3,3>::n4] + Xm[pos<0,0>::n4]*Xm[pos<1,1>::n4]*Xm[pos<3,3>::n4] ) / tmp_det;
max@0:         outm[pos<3,2>::n4] = ( Xm[pos<0,2>::n4]*Xm[pos<1,1>::n4]*Xm[pos<3,0>::n4] - Xm[pos<0,1>::n4]*Xm[pos<1,2>::n4]*Xm[pos<3,0>::n4] - Xm[pos<0,2>::n4]*Xm[pos<1,0>::n4]*Xm[pos<3,1>::n4] + Xm[pos<0,0>::n4]*Xm[pos<1,2>::n4]*Xm[pos<3,1>::n4] + Xm[pos<0,1>::n4]*Xm[pos<1,0>::n4]*Xm[pos<3,2>::n4] - Xm[pos<0,0>::n4]*Xm[pos<1,1>::n4]*Xm[pos<3,2>::n4] ) / tmp_det;
max@0:         
max@0:         outm[pos<0,3>::n4] = ( Xm[pos<0,3>::n4]*Xm[pos<1,2>::n4]*Xm[pos<2,1>::n4] - Xm[pos<0,2>::n4]*Xm[pos<1,3>::n4]*Xm[pos<2,1>::n4] - Xm[pos<0,3>::n4]*Xm[pos<1,1>::n4]*Xm[pos<2,2>::n4] + Xm[pos<0,1>::n4]*Xm[pos<1,3>::n4]*Xm[pos<2,2>::n4] + Xm[pos<0,2>::n4]*Xm[pos<1,1>::n4]*Xm[pos<2,3>::n4] - Xm[pos<0,1>::n4]*Xm[pos<1,2>::n4]*Xm[pos<2,3>::n4] ) / tmp_det;
max@0:         outm[pos<1,3>::n4] = ( Xm[pos<0,2>::n4]*Xm[pos<1,3>::n4]*Xm[pos<2,0>::n4] - Xm[pos<0,3>::n4]*Xm[pos<1,2>::n4]*Xm[pos<2,0>::n4] + Xm[pos<0,3>::n4]*Xm[pos<1,0>::n4]*Xm[pos<2,2>::n4] - Xm[pos<0,0>::n4]*Xm[pos<1,3>::n4]*Xm[pos<2,2>::n4] - Xm[pos<0,2>::n4]*Xm[pos<1,0>::n4]*Xm[pos<2,3>::n4] + Xm[pos<0,0>::n4]*Xm[pos<1,2>::n4]*Xm[pos<2,3>::n4] ) / tmp_det;
max@0:         outm[pos<2,3>::n4] = ( Xm[pos<0,3>::n4]*Xm[pos<1,1>::n4]*Xm[pos<2,0>::n4] - Xm[pos<0,1>::n4]*Xm[pos<1,3>::n4]*Xm[pos<2,0>::n4] - Xm[pos<0,3>::n4]*Xm[pos<1,0>::n4]*Xm[pos<2,1>::n4] + Xm[pos<0,0>::n4]*Xm[pos<1,3>::n4]*Xm[pos<2,1>::n4] + Xm[pos<0,1>::n4]*Xm[pos<1,0>::n4]*Xm[pos<2,3>::n4] - Xm[pos<0,0>::n4]*Xm[pos<1,1>::n4]*Xm[pos<2,3>::n4] ) / tmp_det;
max@0:         outm[pos<3,3>::n4] = ( Xm[pos<0,1>::n4]*Xm[pos<1,2>::n4]*Xm[pos<2,0>::n4] - Xm[pos<0,2>::n4]*Xm[pos<1,1>::n4]*Xm[pos<2,0>::n4] + Xm[pos<0,2>::n4]*Xm[pos<1,0>::n4]*Xm[pos<2,1>::n4] - Xm[pos<0,0>::n4]*Xm[pos<1,2>::n4]*Xm[pos<2,1>::n4] - Xm[pos<0,1>::n4]*Xm[pos<1,0>::n4]*Xm[pos<2,2>::n4] + Xm[pos<0,0>::n4]*Xm[pos<1,1>::n4]*Xm[pos<2,2>::n4] ) / tmp_det;
max@0:         }
max@0:       else
max@0:         {
max@0:         det_ok = false;
max@0:         }
max@0:       };
max@0:       break;
max@0:     
max@0:     default:
max@0:       ;
max@0:     }
max@0:   
max@0:   return det_ok;
max@0:   }
max@0: 
max@0: 
max@0: 
max@0: template<typename eT>
max@0: inline
max@0: bool
max@0: auxlib::inv_inplace_tinymat(Mat<eT>& X, const uword N)
max@0:   {
max@0:   arma_extra_debug_sigprint();
max@0:   
max@0:   bool det_ok = true;
max@0:   
max@0:   // for more info, see:
max@0:   // http://www.dr-lex.34sp.com/random/matrix_inv.html
max@0:   // http://www.cvl.iis.u-tokyo.ac.jp/~miyazaki/tech/teche23.html
max@0:   // http://www.euclideanspace.com/maths/algebra/matrix/functions/inverse/fourD/index.htm
max@0:   // http://www.geometrictools.com//LibFoundation/Mathematics/Wm4Matrix4.inl
max@0:   
max@0:   switch(N)
max@0:     {
max@0:     case 1:
max@0:       {
max@0:       X[0] = eT(1) / X[0];
max@0:       };
max@0:       break;
max@0:       
max@0:     case 2:
max@0:       {
max@0:       const eT a = X[pos<0,0>::n2];
max@0:       const eT b = X[pos<0,1>::n2];
max@0:       const eT c = X[pos<1,0>::n2];
max@0:       const eT d = X[pos<1,1>::n2];
max@0:       
max@0:       const eT tmp_det = (a*d - b*c);
max@0:       
max@0:       if(tmp_det != eT(0))
max@0:         {
max@0:         X[pos<0,0>::n2] =  d / tmp_det;
max@0:         X[pos<0,1>::n2] = -b / tmp_det;
max@0:         X[pos<1,0>::n2] = -c / tmp_det;
max@0:         X[pos<1,1>::n2] =  a / tmp_det;
max@0:         }
max@0:       else
max@0:         {
max@0:         det_ok = false;
max@0:         }
max@0:       };
max@0:       break;
max@0:     
max@0:     case 3:
max@0:       {
max@0:       eT* X_col0 = X.colptr(0);
max@0:       eT* X_col1 = X.colptr(1);
max@0:       eT* X_col2 = X.colptr(2);
max@0:       
max@0:       const eT a11 = X_col0[0];
max@0:       const eT a21 = X_col0[1];
max@0:       const eT a31 = X_col0[2];
max@0:       
max@0:       const eT a12 = X_col1[0];
max@0:       const eT a22 = X_col1[1];
max@0:       const eT a32 = X_col1[2];
max@0:       
max@0:       const eT a13 = X_col2[0];
max@0:       const eT a23 = X_col2[1];
max@0:       const eT a33 = X_col2[2];
max@0:       
max@0:       const eT tmp_det = a11*(a33*a22 - a32*a23) - a21*(a33*a12-a32*a13) + a31*(a23*a12 - a22*a13);
max@0:       
max@0:       if(tmp_det != eT(0))
max@0:         {
max@0:         X_col0[0] =  (a33*a22 - a32*a23) / tmp_det;
max@0:         X_col0[1] = -(a33*a21 - a31*a23) / tmp_det;
max@0:         X_col0[2] =  (a32*a21 - a31*a22) / tmp_det;
max@0:         
max@0:         X_col1[0] = -(a33*a12 - a32*a13) / tmp_det;
max@0:         X_col1[1] =  (a33*a11 - a31*a13) / tmp_det;
max@0:         X_col1[2] = -(a32*a11 - a31*a12) / tmp_det;
max@0:         
max@0:         X_col2[0] =  (a23*a12 - a22*a13) / tmp_det;
max@0:         X_col2[1] = -(a23*a11 - a21*a13) / tmp_det;
max@0:         X_col2[2] =  (a22*a11 - a21*a12) / tmp_det;
max@0:         }
max@0:       else
max@0:         {
max@0:         det_ok = false;
max@0:         }
max@0:       };
max@0:       break;
max@0:       
max@0:     case 4:
max@0:       {
max@0:       const eT tmp_det = det(X);
max@0:       
max@0:       if(tmp_det != eT(0))
max@0:         {
max@0:         const Mat<eT> A(X);
max@0:         
max@0:         const eT* Am = A.memptr();
max@0:               eT* Xm = X.memptr();
max@0:         
max@0:         Xm[pos<0,0>::n4] = ( Am[pos<1,2>::n4]*Am[pos<2,3>::n4]*Am[pos<3,1>::n4] - Am[pos<1,3>::n4]*Am[pos<2,2>::n4]*Am[pos<3,1>::n4] + Am[pos<1,3>::n4]*Am[pos<2,1>::n4]*Am[pos<3,2>::n4] - Am[pos<1,1>::n4]*Am[pos<2,3>::n4]*Am[pos<3,2>::n4] - Am[pos<1,2>::n4]*Am[pos<2,1>::n4]*Am[pos<3,3>::n4] + Am[pos<1,1>::n4]*Am[pos<2,2>::n4]*Am[pos<3,3>::n4] ) / tmp_det;
max@0:         Xm[pos<1,0>::n4] = ( Am[pos<1,3>::n4]*Am[pos<2,2>::n4]*Am[pos<3,0>::n4] - Am[pos<1,2>::n4]*Am[pos<2,3>::n4]*Am[pos<3,0>::n4] - Am[pos<1,3>::n4]*Am[pos<2,0>::n4]*Am[pos<3,2>::n4] + Am[pos<1,0>::n4]*Am[pos<2,3>::n4]*Am[pos<3,2>::n4] + Am[pos<1,2>::n4]*Am[pos<2,0>::n4]*Am[pos<3,3>::n4] - Am[pos<1,0>::n4]*Am[pos<2,2>::n4]*Am[pos<3,3>::n4] ) / tmp_det;
max@0:         Xm[pos<2,0>::n4] = ( Am[pos<1,1>::n4]*Am[pos<2,3>::n4]*Am[pos<3,0>::n4] - Am[pos<1,3>::n4]*Am[pos<2,1>::n4]*Am[pos<3,0>::n4] + Am[pos<1,3>::n4]*Am[pos<2,0>::n4]*Am[pos<3,1>::n4] - Am[pos<1,0>::n4]*Am[pos<2,3>::n4]*Am[pos<3,1>::n4] - Am[pos<1,1>::n4]*Am[pos<2,0>::n4]*Am[pos<3,3>::n4] + Am[pos<1,0>::n4]*Am[pos<2,1>::n4]*Am[pos<3,3>::n4] ) / tmp_det;
max@0:         Xm[pos<3,0>::n4] = ( Am[pos<1,2>::n4]*Am[pos<2,1>::n4]*Am[pos<3,0>::n4] - Am[pos<1,1>::n4]*Am[pos<2,2>::n4]*Am[pos<3,0>::n4] - Am[pos<1,2>::n4]*Am[pos<2,0>::n4]*Am[pos<3,1>::n4] + Am[pos<1,0>::n4]*Am[pos<2,2>::n4]*Am[pos<3,1>::n4] + Am[pos<1,1>::n4]*Am[pos<2,0>::n4]*Am[pos<3,2>::n4] - Am[pos<1,0>::n4]*Am[pos<2,1>::n4]*Am[pos<3,2>::n4] ) / tmp_det;
max@0:         
max@0:         Xm[pos<0,1>::n4] = ( Am[pos<0,3>::n4]*Am[pos<2,2>::n4]*Am[pos<3,1>::n4] - Am[pos<0,2>::n4]*Am[pos<2,3>::n4]*Am[pos<3,1>::n4] - Am[pos<0,3>::n4]*Am[pos<2,1>::n4]*Am[pos<3,2>::n4] + Am[pos<0,1>::n4]*Am[pos<2,3>::n4]*Am[pos<3,2>::n4] + Am[pos<0,2>::n4]*Am[pos<2,1>::n4]*Am[pos<3,3>::n4] - Am[pos<0,1>::n4]*Am[pos<2,2>::n4]*Am[pos<3,3>::n4] ) / tmp_det;
max@0:         Xm[pos<1,1>::n4] = ( Am[pos<0,2>::n4]*Am[pos<2,3>::n4]*Am[pos<3,0>::n4] - Am[pos<0,3>::n4]*Am[pos<2,2>::n4]*Am[pos<3,0>::n4] + Am[pos<0,3>::n4]*Am[pos<2,0>::n4]*Am[pos<3,2>::n4] - Am[pos<0,0>::n4]*Am[pos<2,3>::n4]*Am[pos<3,2>::n4] - Am[pos<0,2>::n4]*Am[pos<2,0>::n4]*Am[pos<3,3>::n4] + Am[pos<0,0>::n4]*Am[pos<2,2>::n4]*Am[pos<3,3>::n4] ) / tmp_det;
max@0:         Xm[pos<2,1>::n4] = ( Am[pos<0,3>::n4]*Am[pos<2,1>::n4]*Am[pos<3,0>::n4] - Am[pos<0,1>::n4]*Am[pos<2,3>::n4]*Am[pos<3,0>::n4] - Am[pos<0,3>::n4]*Am[pos<2,0>::n4]*Am[pos<3,1>::n4] + Am[pos<0,0>::n4]*Am[pos<2,3>::n4]*Am[pos<3,1>::n4] + Am[pos<0,1>::n4]*Am[pos<2,0>::n4]*Am[pos<3,3>::n4] - Am[pos<0,0>::n4]*Am[pos<2,1>::n4]*Am[pos<3,3>::n4] ) / tmp_det;
max@0:         Xm[pos<3,1>::n4] = ( Am[pos<0,1>::n4]*Am[pos<2,2>::n4]*Am[pos<3,0>::n4] - Am[pos<0,2>::n4]*Am[pos<2,1>::n4]*Am[pos<3,0>::n4] + Am[pos<0,2>::n4]*Am[pos<2,0>::n4]*Am[pos<3,1>::n4] - Am[pos<0,0>::n4]*Am[pos<2,2>::n4]*Am[pos<3,1>::n4] - Am[pos<0,1>::n4]*Am[pos<2,0>::n4]*Am[pos<3,2>::n4] + Am[pos<0,0>::n4]*Am[pos<2,1>::n4]*Am[pos<3,2>::n4] ) / tmp_det;
max@0:         
max@0:         Xm[pos<0,2>::n4] = ( Am[pos<0,2>::n4]*Am[pos<1,3>::n4]*Am[pos<3,1>::n4] - Am[pos<0,3>::n4]*Am[pos<1,2>::n4]*Am[pos<3,1>::n4] + Am[pos<0,3>::n4]*Am[pos<1,1>::n4]*Am[pos<3,2>::n4] - Am[pos<0,1>::n4]*Am[pos<1,3>::n4]*Am[pos<3,2>::n4] - Am[pos<0,2>::n4]*Am[pos<1,1>::n4]*Am[pos<3,3>::n4] + Am[pos<0,1>::n4]*Am[pos<1,2>::n4]*Am[pos<3,3>::n4] ) / tmp_det;
max@0:         Xm[pos<1,2>::n4] = ( Am[pos<0,3>::n4]*Am[pos<1,2>::n4]*Am[pos<3,0>::n4] - Am[pos<0,2>::n4]*Am[pos<1,3>::n4]*Am[pos<3,0>::n4] - Am[pos<0,3>::n4]*Am[pos<1,0>::n4]*Am[pos<3,2>::n4] + Am[pos<0,0>::n4]*Am[pos<1,3>::n4]*Am[pos<3,2>::n4] + Am[pos<0,2>::n4]*Am[pos<1,0>::n4]*Am[pos<3,3>::n4] - Am[pos<0,0>::n4]*Am[pos<1,2>::n4]*Am[pos<3,3>::n4] ) / tmp_det;
max@0:         Xm[pos<2,2>::n4] = ( Am[pos<0,1>::n4]*Am[pos<1,3>::n4]*Am[pos<3,0>::n4] - Am[pos<0,3>::n4]*Am[pos<1,1>::n4]*Am[pos<3,0>::n4] + Am[pos<0,3>::n4]*Am[pos<1,0>::n4]*Am[pos<3,1>::n4] - Am[pos<0,0>::n4]*Am[pos<1,3>::n4]*Am[pos<3,1>::n4] - Am[pos<0,1>::n4]*Am[pos<1,0>::n4]*Am[pos<3,3>::n4] + Am[pos<0,0>::n4]*Am[pos<1,1>::n4]*Am[pos<3,3>::n4] ) / tmp_det;
max@0:         Xm[pos<3,2>::n4] = ( Am[pos<0,2>::n4]*Am[pos<1,1>::n4]*Am[pos<3,0>::n4] - Am[pos<0,1>::n4]*Am[pos<1,2>::n4]*Am[pos<3,0>::n4] - Am[pos<0,2>::n4]*Am[pos<1,0>::n4]*Am[pos<3,1>::n4] + Am[pos<0,0>::n4]*Am[pos<1,2>::n4]*Am[pos<3,1>::n4] + Am[pos<0,1>::n4]*Am[pos<1,0>::n4]*Am[pos<3,2>::n4] - Am[pos<0,0>::n4]*Am[pos<1,1>::n4]*Am[pos<3,2>::n4] ) / tmp_det;
max@0:         
max@0:         Xm[pos<0,3>::n4] = ( Am[pos<0,3>::n4]*Am[pos<1,2>::n4]*Am[pos<2,1>::n4] - Am[pos<0,2>::n4]*Am[pos<1,3>::n4]*Am[pos<2,1>::n4] - Am[pos<0,3>::n4]*Am[pos<1,1>::n4]*Am[pos<2,2>::n4] + Am[pos<0,1>::n4]*Am[pos<1,3>::n4]*Am[pos<2,2>::n4] + Am[pos<0,2>::n4]*Am[pos<1,1>::n4]*Am[pos<2,3>::n4] - Am[pos<0,1>::n4]*Am[pos<1,2>::n4]*Am[pos<2,3>::n4] ) / tmp_det;
max@0:         Xm[pos<1,3>::n4] = ( Am[pos<0,2>::n4]*Am[pos<1,3>::n4]*Am[pos<2,0>::n4] - Am[pos<0,3>::n4]*Am[pos<1,2>::n4]*Am[pos<2,0>::n4] + Am[pos<0,3>::n4]*Am[pos<1,0>::n4]*Am[pos<2,2>::n4] - Am[pos<0,0>::n4]*Am[pos<1,3>::n4]*Am[pos<2,2>::n4] - Am[pos<0,2>::n4]*Am[pos<1,0>::n4]*Am[pos<2,3>::n4] + Am[pos<0,0>::n4]*Am[pos<1,2>::n4]*Am[pos<2,3>::n4] ) / tmp_det;
max@0:         Xm[pos<2,3>::n4] = ( Am[pos<0,3>::n4]*Am[pos<1,1>::n4]*Am[pos<2,0>::n4] - Am[pos<0,1>::n4]*Am[pos<1,3>::n4]*Am[pos<2,0>::n4] - Am[pos<0,3>::n4]*Am[pos<1,0>::n4]*Am[pos<2,1>::n4] + Am[pos<0,0>::n4]*Am[pos<1,3>::n4]*Am[pos<2,1>::n4] + Am[pos<0,1>::n4]*Am[pos<1,0>::n4]*Am[pos<2,3>::n4] - Am[pos<0,0>::n4]*Am[pos<1,1>::n4]*Am[pos<2,3>::n4] ) / tmp_det;
max@0:         Xm[pos<3,3>::n4] = ( Am[pos<0,1>::n4]*Am[pos<1,2>::n4]*Am[pos<2,0>::n4] - Am[pos<0,2>::n4]*Am[pos<1,1>::n4]*Am[pos<2,0>::n4] + Am[pos<0,2>::n4]*Am[pos<1,0>::n4]*Am[pos<2,1>::n4] - Am[pos<0,0>::n4]*Am[pos<1,2>::n4]*Am[pos<2,1>::n4] - Am[pos<0,1>::n4]*Am[pos<1,0>::n4]*Am[pos<2,2>::n4] + Am[pos<0,0>::n4]*Am[pos<1,1>::n4]*Am[pos<2,2>::n4] ) / tmp_det;
max@0:         }
max@0:       else
max@0:         {
max@0:         det_ok = false;
max@0:         }
max@0:       };
max@0:       break;
max@0:       
max@0:     default:
max@0:       ;
max@0:     }
max@0:   
max@0:   return det_ok;
max@0:   }
max@0: 
max@0: 
max@0: 
max@0: template<typename eT>
max@0: inline
max@0: bool
max@0: auxlib::inv_inplace_lapack(Mat<eT>& out)
max@0:   {
max@0:   arma_extra_debug_sigprint();
max@0: 
max@0:   if(out.is_empty())
max@0:     {
max@0:     return true;
max@0:     }
max@0:   
max@0:   #if defined(ARMA_USE_ATLAS)
max@0:     {
max@0:     podarray<int> ipiv(out.n_rows);
max@0:     
max@0:     int info = atlas::clapack_getrf(atlas::CblasColMajor, out.n_rows, out.n_cols, out.memptr(), out.n_rows, ipiv.memptr());
max@0:     
max@0:     if(info == 0)
max@0:       {
max@0:       info = atlas::clapack_getri(atlas::CblasColMajor, out.n_rows, out.memptr(), out.n_rows, ipiv.memptr());
max@0:       }
max@0:     
max@0:     return (info == 0);
max@0:     }
max@0:   #elif defined(ARMA_USE_LAPACK)
max@0:     {
max@0:     blas_int n_rows = out.n_rows;
max@0:     blas_int n_cols = out.n_cols;
max@0:     blas_int info   = 0;
max@0:     
max@0:     podarray<blas_int> ipiv(out.n_rows);
max@0:     
max@0:     // 84 was empirically found -- it is the maximum value suggested by LAPACK (as provided by ATLAS v3.6)
max@0:     // based on tests with various matrix types on 32-bit and 64-bit machines
max@0:     //
max@0:     // the "work" array is deliberately long so that a secondary (time-consuming)
max@0:     // memory allocation is avoided, if possible
max@0:     
max@0:     blas_int work_len = (std::max)(blas_int(1), n_rows*84);
max@0:     podarray<eT> work( static_cast<uword>(work_len) );
max@0:     
max@0:     lapack::getrf(&n_rows, &n_cols, out.memptr(), &n_rows, ipiv.memptr(), &info);
max@0:     
max@0:     if(info == 0)
max@0:       {
max@0:       // query for optimum size of work_len
max@0:       
max@0:       blas_int work_len_tmp = -1;
max@0:       lapack::getri(&n_rows, out.memptr(), &n_rows, ipiv.memptr(), work.memptr(), &work_len_tmp, &info);
max@0:       
max@0:       if(info == 0)
max@0:         {
max@0:         blas_int proposed_work_len = static_cast<blas_int>(access::tmp_real(work[0]));
max@0:         
max@0:         // if necessary, allocate more memory
max@0:         if(work_len < proposed_work_len)
max@0:           {
max@0:           work_len = proposed_work_len;
max@0:           work.set_size( static_cast<uword>(work_len) );
max@0:           }
max@0:         }
max@0:       
max@0:       lapack::getri(&n_rows, out.memptr(), &n_rows, ipiv.memptr(), work.memptr(), &work_len, &info);
max@0:       }
max@0:     
max@0:     return (info == 0);
max@0:     }
max@0:   #else
max@0:     {
max@0:     arma_ignore(out);
max@0:     arma_stop("inv(): use of ATLAS or LAPACK needs to be enabled");
max@0:     return false;
max@0:     }
max@0:   #endif
max@0:   }
max@0: 
max@0: 
max@0: 
max@0: template<typename eT, typename T1>
max@0: inline
max@0: bool
max@0: auxlib::inv_tr(Mat<eT>& out, const Base<eT,T1>& X, const uword layout)
max@0:   {
max@0:   arma_extra_debug_sigprint();
max@0:   
max@0:   out = X.get_ref();
max@0:   
max@0:   arma_debug_check( (out.is_square() == false), "inv(): given matrix is not square" );
max@0:   
max@0:   if(out.is_empty())
max@0:     {
max@0:     return true;
max@0:     }
max@0:   
max@0:   bool status;
max@0:   
max@0:   #if defined(ARMA_USE_LAPACK)
max@0:     {
max@0:     char     uplo = (layout == 0) ? 'U' : 'L';
max@0:     char     diag = 'N';
max@0:     blas_int n    = blas_int(out.n_rows);
max@0:     blas_int info = 0;
max@0:     
max@0:     lapack::trtri(&uplo, &diag, &n, out.memptr(), &n, &info);
max@0:     
max@0:     status = (info == 0);
max@0:     }
max@0:   #else
max@0:     {
max@0:     arma_ignore(layout);
max@0:     arma_stop("inv(): use of LAPACK needs to be enabled");
max@0:     status = false;
max@0:     }
max@0:   #endif
max@0:   
max@0:   
max@0:   if(status == true)
max@0:     {
max@0:     if(layout == 0)
max@0:       {
max@0:       // upper triangular
max@0:       out = trimatu(out);
max@0:       }
max@0:     else
max@0:       {
max@0:       // lower triangular      
max@0:       out = trimatl(out);
max@0:       }
max@0:     }
max@0:   
max@0:   return status;
max@0:   }
max@0: 
max@0: 
max@0: 
max@0: template<typename eT, typename T1>
max@0: inline
max@0: bool
max@0: auxlib::inv_sym(Mat<eT>& out, const Base<eT,T1>& X, const uword layout)
max@0:   {
max@0:   arma_extra_debug_sigprint();
max@0:   
max@0:   out = X.get_ref();
max@0:   
max@0:   arma_debug_check( (out.is_square() == false), "inv(): given matrix is not square" );
max@0:   
max@0:   if(out.is_empty())
max@0:     {
max@0:     return true;
max@0:     }
max@0:   
max@0:   bool status;
max@0:   
max@0:   #if defined(ARMA_USE_LAPACK)
max@0:     {
max@0:     char     uplo  = (layout == 0) ? 'U' : 'L';
max@0:     blas_int n     = blas_int(out.n_rows);
max@0:     blas_int lwork = n*n; // TODO: use lwork = -1 to determine optimal size
max@0:     blas_int info  = 0;
max@0:     
max@0:     podarray<blas_int> ipiv;
max@0:     ipiv.set_size(out.n_rows);
max@0:     
max@0:     podarray<eT> work;
max@0:     work.set_size( uword(lwork) );
max@0:     
max@0:     lapack::sytrf(&uplo, &n, out.memptr(), &n, ipiv.memptr(), work.memptr(), &lwork, &info);
max@0:     
max@0:     status = (info == 0);
max@0:     
max@0:     if(status == true)
max@0:       {
max@0:       lapack::sytri(&uplo, &n, out.memptr(), &n, ipiv.memptr(), work.memptr(), &info);
max@0:       
max@0:       out = (layout == 0) ? symmatu(out) : symmatl(out);
max@0:       
max@0:       status = (info == 0);
max@0:       }
max@0:     }
max@0:   #else
max@0:     {
max@0:     arma_ignore(layout);
max@0:     arma_stop("inv(): use of LAPACK needs to be enabled");
max@0:     status = false;
max@0:     }
max@0:   #endif
max@0:   
max@0:   return status;
max@0:   }
max@0: 
max@0: 
max@0: 
max@0: template<typename eT, typename T1>
max@0: inline
max@0: bool
max@0: auxlib::inv_sympd(Mat<eT>& out, const Base<eT,T1>& X, const uword layout)
max@0:   {
max@0:   arma_extra_debug_sigprint();
max@0:   
max@0:   out = X.get_ref();
max@0:   
max@0:   arma_debug_check( (out.is_square() == false), "inv(): given matrix is not square" );
max@0:   
max@0:   if(out.is_empty())
max@0:     {
max@0:     return true;
max@0:     }
max@0:   
max@0:   bool status;
max@0:   
max@0:   #if defined(ARMA_USE_LAPACK)
max@0:     {
max@0:     char     uplo = (layout == 0) ? 'U' : 'L';
max@0:     blas_int n    = blas_int(out.n_rows);
max@0:     blas_int info = 0;
max@0:     
max@0:     lapack::potrf(&uplo, &n, out.memptr(), &n, &info);
max@0:     
max@0:     status = (info == 0);
max@0:     
max@0:     if(status == true)
max@0:       {
max@0:       lapack::potri(&uplo, &n, out.memptr(), &n, &info);
max@0:     
max@0:       out = (layout == 0) ? symmatu(out) : symmatl(out);
max@0:     
max@0:       status = (info == 0);
max@0:       }
max@0:     }
max@0:   #else
max@0:     {
max@0:     arma_ignore(layout);
max@0:     arma_stop("inv(): use of LAPACK needs to be enabled");
max@0:     status = false;
max@0:     }
max@0:   #endif
max@0:   
max@0:   return status;
max@0:   }
max@0: 
max@0: 
max@0: 
max@0: template<typename eT, typename T1>
max@0: inline
max@0: eT
max@0: auxlib::det(const Base<eT,T1>& X, const bool slow)
max@0:   {
max@0:   const unwrap<T1>   tmp(X.get_ref());
max@0:   const Mat<eT>& A = tmp.M;
max@0:   
max@0:   arma_debug_check( (A.is_square() == false), "det(): matrix is not square" );
max@0:   
max@0:   const bool make_copy = (is_Mat<T1>::value == true) ? true : false;
max@0:   
max@0:   if(slow == false)
max@0:     {
max@0:     const uword N = A.n_rows;
max@0:     
max@0:     switch(N)
max@0:       {
max@0:       case 0:
max@0:       case 1:
max@0:       case 2:
max@0:         return auxlib::det_tinymat(A, N);
max@0:         break;
max@0:       
max@0:       case 3:
max@0:       case 4:
max@0:         {
max@0:         const eT tmp_det = auxlib::det_tinymat(A, N);
max@0:         return (tmp_det != eT(0)) ? tmp_det : auxlib::det_lapack(A, make_copy);
max@0:         }
max@0:         break;
max@0:       
max@0:       default:
max@0:         return auxlib::det_lapack(A, make_copy);
max@0:       }
max@0:     }
max@0:   else
max@0:     {
max@0:     return auxlib::det_lapack(A, make_copy);
max@0:     }
max@0:   }
max@0: 
max@0: 
max@0: 
max@0: template<typename eT>
max@0: inline
max@0: eT
max@0: auxlib::det_tinymat(const Mat<eT>& X, const uword N)
max@0:   {
max@0:   arma_extra_debug_sigprint();
max@0:   
max@0:   switch(N)
max@0:     {
max@0:     case 0:
max@0:       return eT(1);
max@0:       break;
max@0:     
max@0:     case 1:
max@0:       return X[0];
max@0:       break;
max@0:     
max@0:     case 2:
max@0:       {
max@0:       const eT* Xm = X.memptr();
max@0:       
max@0:       return ( Xm[pos<0,0>::n2]*Xm[pos<1,1>::n2] - Xm[pos<0,1>::n2]*Xm[pos<1,0>::n2] );
max@0:       }
max@0:       break;
max@0:     
max@0:     case 3:
max@0:       {
max@0:       // const double tmp1 = X.at(0,0) * X.at(1,1) * X.at(2,2);
max@0:       // const double tmp2 = X.at(0,1) * X.at(1,2) * X.at(2,0);
max@0:       // const double tmp3 = X.at(0,2) * X.at(1,0) * X.at(2,1);
max@0:       // const double tmp4 = X.at(2,0) * X.at(1,1) * X.at(0,2);
max@0:       // const double tmp5 = X.at(2,1) * X.at(1,2) * X.at(0,0);
max@0:       // const double tmp6 = X.at(2,2) * X.at(1,0) * X.at(0,1);
max@0:       // return (tmp1+tmp2+tmp3) - (tmp4+tmp5+tmp6);
max@0:       
max@0:       const eT* a_col0 = X.colptr(0);
max@0:       const eT  a11    = a_col0[0];
max@0:       const eT  a21    = a_col0[1];
max@0:       const eT  a31    = a_col0[2];
max@0:       
max@0:       const eT* a_col1 = X.colptr(1);
max@0:       const eT  a12    = a_col1[0];
max@0:       const eT  a22    = a_col1[1];
max@0:       const eT  a32    = a_col1[2];
max@0:       
max@0:       const eT* a_col2 = X.colptr(2);
max@0:       const eT  a13    = a_col2[0];
max@0:       const eT  a23    = a_col2[1];
max@0:       const eT  a33    = a_col2[2];
max@0:       
max@0:       return ( a11*(a33*a22 - a32*a23) - a21*(a33*a12-a32*a13) + a31*(a23*a12 - a22*a13) );
max@0:       }
max@0:       break;
max@0:       
max@0:     case 4:
max@0:       {
max@0:       const eT* Xm = X.memptr();
max@0:       
max@0:       const eT val = \
max@0:           Xm[pos<0,3>::n4] * Xm[pos<1,2>::n4] * Xm[pos<2,1>::n4] * Xm[pos<3,0>::n4] \
max@0:         - Xm[pos<0,2>::n4] * Xm[pos<1,3>::n4] * Xm[pos<2,1>::n4] * Xm[pos<3,0>::n4] \
max@0:         - Xm[pos<0,3>::n4] * Xm[pos<1,1>::n4] * Xm[pos<2,2>::n4] * Xm[pos<3,0>::n4] \
max@0:         + Xm[pos<0,1>::n4] * Xm[pos<1,3>::n4] * Xm[pos<2,2>::n4] * Xm[pos<3,0>::n4] \
max@0:         + Xm[pos<0,2>::n4] * Xm[pos<1,1>::n4] * Xm[pos<2,3>::n4] * Xm[pos<3,0>::n4] \
max@0:         - Xm[pos<0,1>::n4] * Xm[pos<1,2>::n4] * Xm[pos<2,3>::n4] * Xm[pos<3,0>::n4] \
max@0:         - Xm[pos<0,3>::n4] * Xm[pos<1,2>::n4] * Xm[pos<2,0>::n4] * Xm[pos<3,1>::n4] \
max@0:         + Xm[pos<0,2>::n4] * Xm[pos<1,3>::n4] * Xm[pos<2,0>::n4] * Xm[pos<3,1>::n4] \
max@0:         + Xm[pos<0,3>::n4] * Xm[pos<1,0>::n4] * Xm[pos<2,2>::n4] * Xm[pos<3,1>::n4] \
max@0:         - Xm[pos<0,0>::n4] * Xm[pos<1,3>::n4] * Xm[pos<2,2>::n4] * Xm[pos<3,1>::n4] \
max@0:         - Xm[pos<0,2>::n4] * Xm[pos<1,0>::n4] * Xm[pos<2,3>::n4] * Xm[pos<3,1>::n4] \
max@0:         + Xm[pos<0,0>::n4] * Xm[pos<1,2>::n4] * Xm[pos<2,3>::n4] * Xm[pos<3,1>::n4] \
max@0:         + Xm[pos<0,3>::n4] * Xm[pos<1,1>::n4] * Xm[pos<2,0>::n4] * Xm[pos<3,2>::n4] \
max@0:         - Xm[pos<0,1>::n4] * Xm[pos<1,3>::n4] * Xm[pos<2,0>::n4] * Xm[pos<3,2>::n4] \
max@0:         - Xm[pos<0,3>::n4] * Xm[pos<1,0>::n4] * Xm[pos<2,1>::n4] * Xm[pos<3,2>::n4] \
max@0:         + Xm[pos<0,0>::n4] * Xm[pos<1,3>::n4] * Xm[pos<2,1>::n4] * Xm[pos<3,2>::n4] \
max@0:         + Xm[pos<0,1>::n4] * Xm[pos<1,0>::n4] * Xm[pos<2,3>::n4] * Xm[pos<3,2>::n4] \
max@0:         - Xm[pos<0,0>::n4] * Xm[pos<1,1>::n4] * Xm[pos<2,3>::n4] * Xm[pos<3,2>::n4] \
max@0:         - Xm[pos<0,2>::n4] * Xm[pos<1,1>::n4] * Xm[pos<2,0>::n4] * Xm[pos<3,3>::n4] \
max@0:         + Xm[pos<0,1>::n4] * Xm[pos<1,2>::n4] * Xm[pos<2,0>::n4] * Xm[pos<3,3>::n4] \
max@0:         + Xm[pos<0,2>::n4] * Xm[pos<1,0>::n4] * Xm[pos<2,1>::n4] * Xm[pos<3,3>::n4] \
max@0:         - Xm[pos<0,0>::n4] * Xm[pos<1,2>::n4] * Xm[pos<2,1>::n4] * Xm[pos<3,3>::n4] \
max@0:         - Xm[pos<0,1>::n4] * Xm[pos<1,0>::n4] * Xm[pos<2,2>::n4] * Xm[pos<3,3>::n4] \
max@0:         + Xm[pos<0,0>::n4] * Xm[pos<1,1>::n4] * Xm[pos<2,2>::n4] * Xm[pos<3,3>::n4] \
max@0:         ;
max@0:       
max@0:       return val;
max@0:       }
max@0:       break;
max@0:     
max@0:     default:
max@0:       return eT(0);
max@0:       ;
max@0:     }
max@0:   }
max@0: 
max@0: 
max@0: 
max@0: //! immediate determinant of a matrix using ATLAS or LAPACK
max@0: template<typename eT>
max@0: inline
max@0: eT
max@0: auxlib::det_lapack(const Mat<eT>& X, const bool make_copy)
max@0:   {
max@0:   arma_extra_debug_sigprint();
max@0:   
max@0:   Mat<eT> X_copy;
max@0:   
max@0:   if(make_copy == true)
max@0:     {
max@0:     X_copy = X;
max@0:     }
max@0:   
max@0:   Mat<eT>& tmp = (make_copy == true) ? X_copy : const_cast< Mat<eT>& >(X);
max@0:   
max@0:   if(tmp.is_empty())
max@0:     {
max@0:     return eT(1);
max@0:     }
max@0:   
max@0:   
max@0:   #if defined(ARMA_USE_ATLAS)
max@0:     {
max@0:     podarray<int> ipiv(tmp.n_rows);
max@0:     
max@0:     //const int info =
max@0:     atlas::clapack_getrf(atlas::CblasColMajor, tmp.n_rows, tmp.n_cols, tmp.memptr(), tmp.n_rows, ipiv.memptr());
max@0:     
max@0:     // on output tmp appears to be L+U_alt, where U_alt is U with the main diagonal set to zero
max@0:     eT val = tmp.at(0,0);
max@0:     for(uword i=1; i < tmp.n_rows; ++i)
max@0:       {
max@0:       val *= tmp.at(i,i);
max@0:       }
max@0:     
max@0:     int sign = +1;
max@0:     for(uword i=0; i < tmp.n_rows; ++i)
max@0:       {
max@0:       if( int(i) != ipiv.mem[i] )  // NOTE: no adjustment required, as the clapack version of getrf() assumes counting from 0
max@0:         {
max@0:         sign *= -1;
max@0:         }
max@0:       }
max@0:     
max@0:     return ( (sign < 0) ? -val : val );
max@0:     }
max@0:   #elif defined(ARMA_USE_LAPACK)
max@0:     {
max@0:     podarray<blas_int> ipiv(tmp.n_rows);
max@0:     
max@0:     blas_int info   = 0;
max@0:     blas_int n_rows = blas_int(tmp.n_rows);
max@0:     blas_int n_cols = blas_int(tmp.n_cols);
max@0:     
max@0:     lapack::getrf(&n_rows, &n_cols, tmp.memptr(), &n_rows, ipiv.memptr(), &info);
max@0:     
max@0:     // on output tmp appears to be L+U_alt, where U_alt is U with the main diagonal set to zero
max@0:     eT val = tmp.at(0,0);
max@0:     for(uword i=1; i < tmp.n_rows; ++i)
max@0:       {
max@0:       val *= tmp.at(i,i);
max@0:       }
max@0:     
max@0:     blas_int sign = +1;
max@0:     for(uword i=0; i < tmp.n_rows; ++i)
max@0:       {
max@0:       if( blas_int(i) != (ipiv.mem[i] - 1) )  // NOTE: adjustment of -1 is required as Fortran counts from 1
max@0:         {
max@0:         sign *= -1;
max@0:         }
max@0:       }
max@0:     
max@0:     return ( (sign < 0) ? -val : val );
max@0:     }
max@0:   #else
max@0:     {
max@0:     arma_ignore(X);
max@0:     arma_ignore(make_copy);
max@0:     arma_ignore(tmp);
max@0:     arma_stop("det(): use of ATLAS or LAPACK needs to be enabled");
max@0:     return eT(0);
max@0:     }
max@0:   #endif
max@0:   }
max@0: 
max@0: 
max@0: 
max@0: //! immediate log determinant of a matrix using ATLAS or LAPACK
max@0: template<typename eT, typename T1>
max@0: inline
max@0: bool
max@0: auxlib::log_det(eT& out_val, typename get_pod_type<eT>::result& out_sign, const Base<eT,T1>& X)
max@0:   {
max@0:   arma_extra_debug_sigprint();
max@0:   
max@0:   typedef typename get_pod_type<eT>::result T;
max@0:   
max@0:   #if defined(ARMA_USE_ATLAS)
max@0:     {
max@0:     Mat<eT> tmp(X.get_ref());
max@0:     arma_debug_check( (tmp.is_square() == false), "log_det(): given matrix is not square" );
max@0:     
max@0:     if(tmp.is_empty())
max@0:       {
max@0:       out_val  = eT(0);
max@0:       out_sign =  T(1);
max@0:       return true;
max@0:       }
max@0:     
max@0:     podarray<int> ipiv(tmp.n_rows);
max@0:     
max@0:     const int info = atlas::clapack_getrf(atlas::CblasColMajor, tmp.n_rows, tmp.n_cols, tmp.memptr(), tmp.n_rows, ipiv.memptr());
max@0:     
max@0:     // on output tmp appears to be L+U_alt, where U_alt is U with the main diagonal set to zero
max@0:     
max@0:     sword sign = (is_complex<eT>::value == false) ? ( (access::tmp_real( tmp.at(0,0) ) < T(0)) ? -1 : +1 ) : +1;
max@0:     eT   val = (is_complex<eT>::value == false) ? std::log( (access::tmp_real( tmp.at(0,0) ) < T(0)) ? tmp.at(0,0)*T(-1) : tmp.at(0,0) ) : std::log( tmp.at(0,0) );
max@0:     
max@0:     for(uword i=1; i < tmp.n_rows; ++i)
max@0:       {
max@0:       const eT x = tmp.at(i,i);
max@0:       
max@0:       sign *= (is_complex<eT>::value == false) ? ( (access::tmp_real(x) < T(0)) ? -1 : +1 ) : +1;
max@0:       val  += (is_complex<eT>::value == false) ? std::log( (access::tmp_real(x) < T(0)) ? x*T(-1) : x ) : std::log(x);
max@0:       }
max@0:     
max@0:     for(uword i=0; i < tmp.n_rows; ++i)
max@0:       {
max@0:       if( int(i) != ipiv.mem[i] )  // NOTE: no adjustment required, as the clapack version of getrf() assumes counting from 0
max@0:         {
max@0:         sign *= -1;
max@0:         }
max@0:       }
max@0:     
max@0:     out_val  = val;
max@0:     out_sign = T(sign);
max@0:     
max@0:     return (info == 0);
max@0:     }
max@0:   #elif defined(ARMA_USE_LAPACK)
max@0:     {
max@0:     Mat<eT> tmp(X.get_ref());
max@0:     arma_debug_check( (tmp.is_square() == false), "log_det(): given matrix is not square" );
max@0:     
max@0:     if(tmp.is_empty())
max@0:       {
max@0:       out_val  = eT(0);
max@0:       out_sign =  T(1);
max@0:       return true;
max@0:       }
max@0:     
max@0:     podarray<blas_int> ipiv(tmp.n_rows);
max@0:     
max@0:     blas_int info   = 0;
max@0:     blas_int n_rows = blas_int(tmp.n_rows);
max@0:     blas_int n_cols = blas_int(tmp.n_cols);
max@0:     
max@0:     lapack::getrf(&n_rows, &n_cols, tmp.memptr(), &n_rows, ipiv.memptr(), &info);
max@0:     
max@0:     // on output tmp appears to be L+U_alt, where U_alt is U with the main diagonal set to zero
max@0:     
max@0:     sword sign = (is_complex<eT>::value == false) ? ( (access::tmp_real( tmp.at(0,0) ) < T(0)) ? -1 : +1 ) : +1;
max@0:     eT   val = (is_complex<eT>::value == false) ? std::log( (access::tmp_real( tmp.at(0,0) ) < T(0)) ? tmp.at(0,0)*T(-1) : tmp.at(0,0) ) : std::log( tmp.at(0,0) );
max@0:     
max@0:     for(uword i=1; i < tmp.n_rows; ++i)
max@0:       {
max@0:       const eT x = tmp.at(i,i);
max@0:       
max@0:       sign *= (is_complex<eT>::value == false) ? ( (access::tmp_real(x) < T(0)) ? -1 : +1 ) : +1;
max@0:       val  += (is_complex<eT>::value == false) ? std::log( (access::tmp_real(x) < T(0)) ? x*T(-1) : x ) : std::log(x);
max@0:       }
max@0:     
max@0:     for(uword i=0; i < tmp.n_rows; ++i)
max@0:       {
max@0:       if( blas_int(i) != (ipiv.mem[i] - 1) )  // NOTE: adjustment of -1 is required as Fortran counts from 1
max@0:         {
max@0:         sign *= -1;
max@0:         }
max@0:       }
max@0:     
max@0:     out_val  = val;
max@0:     out_sign = T(sign);
max@0:     
max@0:     return (info == 0);
max@0:     }
max@0:   #else
max@0:     {
max@0:     out_val  = eT(0);
max@0:     out_sign =  T(0);
max@0:     
max@0:     arma_stop("log_det(): use of ATLAS or LAPACK needs to be enabled");
max@0:     
max@0:     return false;
max@0:     }
max@0:   #endif
max@0:   }
max@0: 
max@0: 
max@0: 
max@0: //! immediate LU decomposition of a matrix using ATLAS or LAPACK
max@0: template<typename eT, typename T1>
max@0: inline
max@0: bool
max@0: auxlib::lu(Mat<eT>& L, Mat<eT>& U, podarray<blas_int>& ipiv, const Base<eT,T1>& X)
max@0:   {
max@0:   arma_extra_debug_sigprint();
max@0:   
max@0:   U = X.get_ref();
max@0:   
max@0:   const uword U_n_rows = U.n_rows;
max@0:   const uword U_n_cols = U.n_cols;
max@0:   
max@0:   if(U.is_empty())
max@0:     {
max@0:     L.set_size(U_n_rows, 0);
max@0:     U.set_size(0, U_n_cols);
max@0:     ipiv.reset();
max@0:     return true;
max@0:     }
max@0:   
max@0:   #if defined(ARMA_USE_ATLAS) || defined(ARMA_USE_LAPACK)
max@0:     {
max@0:     bool status;
max@0:     
max@0:     #if defined(ARMA_USE_ATLAS)
max@0:       {
max@0:       ipiv.set_size( (std::min)(U_n_rows, U_n_cols) );
max@0:       
max@0:       int info = atlas::clapack_getrf(atlas::CblasColMajor, U_n_rows, U_n_cols, U.memptr(), U_n_rows, ipiv.memptr());
max@0:       
max@0:       status = (info == 0);
max@0:       }
max@0:     #elif defined(ARMA_USE_LAPACK)
max@0:       {
max@0:       ipiv.set_size( (std::min)(U_n_rows, U_n_cols) );
max@0:       
max@0:       blas_int info = 0;
max@0:       
max@0:       blas_int n_rows = U_n_rows;
max@0:       blas_int n_cols = U_n_cols;
max@0:       
max@0:       
max@0:       lapack::getrf(&n_rows, &n_cols, U.memptr(), &n_rows, ipiv.memptr(), &info);
max@0:       
max@0:       // take into account that Fortran counts from 1
max@0:       arrayops::inplace_minus(ipiv.memptr(), blas_int(1), ipiv.n_elem);
max@0:       
max@0:       status = (info == 0);
max@0:       }
max@0:     #endif
max@0:     
max@0:     L.copy_size(U);
max@0:     
max@0:     for(uword col=0; col < U_n_cols; ++col)
max@0:       {
max@0:       for(uword row=0; (row < col) && (row < U_n_rows); ++row)
max@0:         {
max@0:         L.at(row,col) = eT(0);
max@0:         }
max@0:       
max@0:       if( L.in_range(col,col) == true )
max@0:         {
max@0:         L.at(col,col) = eT(1);
max@0:         }
max@0:       
max@0:       for(uword row = (col+1); row < U_n_rows; ++row)
max@0:         {
max@0:         L.at(row,col) = U.at(row,col);
max@0:         U.at(row,col) = eT(0);
max@0:         }
max@0:       }
max@0:     
max@0:     return status;
max@0:     }
max@0:   #else
max@0:     {
max@0:     arma_stop("lu(): use of ATLAS or LAPACK needs to be enabled");
max@0:     
max@0:     return false;
max@0:     }
max@0:   #endif
max@0:   }
max@0: 
max@0: 
max@0: 
max@0: template<typename eT, typename T1>
max@0: inline
max@0: bool
max@0: auxlib::lu(Mat<eT>& L, Mat<eT>& U, Mat<eT>& P, const Base<eT,T1>& X)
max@0:   {
max@0:   arma_extra_debug_sigprint();
max@0:   
max@0:   podarray<blas_int> ipiv1;
max@0:   const bool status = auxlib::lu(L, U, ipiv1, X);
max@0:   
max@0:   if(status == true)
max@0:     {
max@0:     if(U.is_empty())
max@0:       {
max@0:       // L and U have been already set to the correct empty matrices
max@0:       P.eye(L.n_rows, L.n_rows);
max@0:       return true;
max@0:       }
max@0:     
max@0:     const uword n      = ipiv1.n_elem;
max@0:     const uword P_rows = U.n_rows;
max@0:     
max@0:     podarray<blas_int> ipiv2(P_rows);
max@0:     
max@0:     const blas_int* ipiv1_mem = ipiv1.memptr();
max@0:           blas_int* ipiv2_mem = ipiv2.memptr();
max@0:     
max@0:     for(uword i=0; i<P_rows; ++i)
max@0:       {
max@0:       ipiv2_mem[i] = blas_int(i);
max@0:       }
max@0:     
max@0:     for(uword i=0; i<n; ++i)
max@0:       {
max@0:       const uword k = static_cast<uword>(ipiv1_mem[i]);
max@0:       
max@0:       if( ipiv2_mem[i] != ipiv2_mem[k] )
max@0:         {
max@0:         std::swap( ipiv2_mem[i], ipiv2_mem[k] );
max@0:         }
max@0:       }
max@0:     
max@0:     P.zeros(P_rows, P_rows);
max@0:     
max@0:     for(uword row=0; row<P_rows; ++row)
max@0:       {
max@0:       P.at(row, static_cast<uword>(ipiv2_mem[row])) = eT(1);
max@0:       }
max@0:     
max@0:     if(L.n_cols > U.n_rows)
max@0:       {
max@0:       L.shed_cols(U.n_rows, L.n_cols-1);
max@0:       }
max@0:       
max@0:     if(U.n_rows > L.n_cols)
max@0:       {
max@0:       U.shed_rows(L.n_cols, U.n_rows-1);
max@0:       }
max@0:     }
max@0:   
max@0:   return status;
max@0:   }
max@0: 
max@0: 
max@0: 
max@0: template<typename eT, typename T1>
max@0: inline
max@0: bool
max@0: auxlib::lu(Mat<eT>& L, Mat<eT>& U, const Base<eT,T1>& X)
max@0:   {
max@0:   arma_extra_debug_sigprint();
max@0:   
max@0:   podarray<blas_int> ipiv1;
max@0:   const bool status = auxlib::lu(L, U, ipiv1, X);
max@0:   
max@0:   if(status == true)
max@0:     {
max@0:     if(U.is_empty())
max@0:       {
max@0:       // L and U have been already set to the correct empty matrices
max@0:       return true;
max@0:       }
max@0:     
max@0:     const uword n      = ipiv1.n_elem;
max@0:     const uword P_rows = U.n_rows;
max@0:     
max@0:     podarray<blas_int> ipiv2(P_rows);
max@0:     
max@0:     const blas_int* ipiv1_mem = ipiv1.memptr();
max@0:           blas_int* ipiv2_mem = ipiv2.memptr();
max@0:     
max@0:     for(uword i=0; i<P_rows; ++i)
max@0:       {
max@0:       ipiv2_mem[i] = blas_int(i);
max@0:       }
max@0:     
max@0:     for(uword i=0; i<n; ++i)
max@0:       {
max@0:       const uword k = static_cast<uword>(ipiv1_mem[i]);
max@0:       
max@0:       if( ipiv2_mem[i] != ipiv2_mem[k] )
max@0:         {
max@0:         std::swap( ipiv2_mem[i], ipiv2_mem[k] );
max@0:         L.swap_rows( static_cast<uword>(ipiv2_mem[i]), static_cast<uword>(ipiv2_mem[k]) );
max@0:         }
max@0:       }
max@0:     
max@0:     if(L.n_cols > U.n_rows)
max@0:       {
max@0:       L.shed_cols(U.n_rows, L.n_cols-1);
max@0:       }
max@0:       
max@0:     if(U.n_rows > L.n_cols)
max@0:       {
max@0:       U.shed_rows(L.n_cols, U.n_rows-1);
max@0:       }
max@0:     }
max@0:   
max@0:   return status;
max@0:   }
max@0: 
max@0: 
max@0: 
max@0: //! immediate eigenvalues of a symmetric real matrix using LAPACK
max@0: template<typename eT, typename T1>
max@0: inline
max@0: bool
max@0: auxlib::eig_sym(Col<eT>& eigval, const Base<eT,T1>& X)
max@0:   {
max@0:   arma_extra_debug_sigprint();
max@0:   
max@0:   #if defined(ARMA_USE_LAPACK)
max@0:     {
max@0:     Mat<eT> A(X.get_ref());
max@0:     
max@0:     arma_debug_check( (A.is_square() == false), "eig_sym(): given matrix is not square");
max@0:     
max@0:     if(A.is_empty())
max@0:       {
max@0:       eigval.reset();
max@0:       return true;
max@0:       }
max@0:     
max@0:     // rudimentary "better-than-nothing" test for symmetry
max@0:     //arma_debug_check( (A.at(A.n_rows-1, 0) != A.at(0, A.n_cols-1)), "auxlib::eig(): given matrix is not symmetric" );
max@0:     
max@0:     char jobz  = 'N';
max@0:     char uplo  = 'U';
max@0:     
max@0:     blas_int n_rows = A.n_rows;
max@0:     blas_int lwork  = (std::max)(blas_int(1), 3*n_rows-1);
max@0:     
max@0:     eigval.set_size( static_cast<uword>(n_rows) );
max@0:     podarray<eT> work( static_cast<uword>(lwork) );
max@0:     
max@0:     blas_int info;
max@0:     
max@0:     arma_extra_debug_print("lapack::syev()");
max@0:     lapack::syev(&jobz, &uplo, &n_rows, A.memptr(), &n_rows, eigval.memptr(), work.memptr(), &lwork, &info);
max@0:     
max@0:     return (info == 0);
max@0:     }
max@0:   #else
max@0:     {
max@0:     arma_ignore(eigval);
max@0:     arma_ignore(X);
max@0:     arma_stop("eig_sym(): use of LAPACK needs to be enabled");
max@0:     return false;
max@0:     }
max@0:   #endif
max@0:   }
max@0: 
max@0: 
max@0: 
max@0: //! immediate eigenvalues of a hermitian complex matrix using LAPACK
max@0: template<typename T, typename T1>
max@0: inline
max@0: bool
max@0: auxlib::eig_sym(Col<T>& eigval, const Base<std::complex<T>,T1>& X)
max@0:   {
max@0:   arma_extra_debug_sigprint();
max@0:   
max@0:   typedef typename std::complex<T> eT;
max@0:   
max@0:   #if defined(ARMA_USE_LAPACK)
max@0:     {
max@0:     Mat<eT> A(X.get_ref());
max@0:     arma_debug_check( (A.is_square() == false), "eig_sym(): given matrix is not hermitian");
max@0:     
max@0:     if(A.is_empty())
max@0:       {
max@0:       eigval.reset();
max@0:       return true;
max@0:       }
max@0:     
max@0:     char jobz  = 'N'; 
max@0:     char uplo  = 'U';
max@0:     
max@0:     blas_int n_rows = A.n_rows;
max@0:     blas_int lda    = A.n_rows;
max@0:     blas_int lwork  = (std::max)(blas_int(1), 2*n_rows - 1);  // TODO: automatically find best size of lwork
max@0:     
max@0:     eigval.set_size( static_cast<uword>(n_rows) );
max@0:     
max@0:     podarray<eT> work( static_cast<uword>(lwork) );
max@0:     podarray<T>  rwork( static_cast<uword>((std::max)(blas_int(1), 3*n_rows - 2)) );
max@0:     
max@0:     blas_int info;
max@0:     
max@0:     arma_extra_debug_print("lapack::heev()");
max@0:     lapack::heev(&jobz, &uplo, &n_rows, A.memptr(), &lda, eigval.memptr(), work.memptr(), &lwork, rwork.memptr(), &info);
max@0:     
max@0:     return (info == 0);
max@0:     }
max@0:   #else
max@0:     {
max@0:     arma_ignore(eigval);
max@0:     arma_ignore(X);
max@0:     arma_stop("eig_sym(): use of LAPACK needs to be enabled");
max@0:     return false;
max@0:     }
max@0:   #endif
max@0:   }
max@0: 
max@0: 
max@0: 
max@0: //! immediate eigenvalues and eigenvectors of a symmetric real matrix using LAPACK
max@0: template<typename eT, typename T1>
max@0: inline
max@0: bool
max@0: auxlib::eig_sym(Col<eT>& eigval, Mat<eT>& eigvec, const Base<eT,T1>& X)
max@0:   {
max@0:   arma_extra_debug_sigprint();
max@0:   
max@0:   #if defined(ARMA_USE_LAPACK)
max@0:     {
max@0:     eigvec = X.get_ref();
max@0:     
max@0:     arma_debug_check( (eigvec.is_square() == false), "eig_sym(): given matrix is not square" );
max@0:     
max@0:     if(eigvec.is_empty())
max@0:       {
max@0:       eigval.reset();
max@0:       eigvec.reset();
max@0:       return true;
max@0:       }
max@0:     
max@0:     // rudimentary "better-than-nothing" test for symmetry
max@0:     //arma_debug_check( (A.at(A.n_rows-1, 0) != A.at(0, A.n_cols-1)), "auxlib::eig(): given matrix is not symmetric" );
max@0:     
max@0:     char jobz  = 'V';
max@0:     char uplo  = 'U';
max@0:     
max@0:     blas_int n_rows = eigvec.n_rows;
max@0:     blas_int lwork  = (std::max)(blas_int(1), 3*n_rows-1);
max@0:     
max@0:     eigval.set_size( static_cast<uword>(n_rows) );
max@0:     podarray<eT> work( static_cast<uword>(lwork) );
max@0:     
max@0:     blas_int info;
max@0:     
max@0:     arma_extra_debug_print("lapack::syev()");
max@0:     lapack::syev(&jobz, &uplo, &n_rows, eigvec.memptr(), &n_rows, eigval.memptr(), work.memptr(), &lwork, &info);
max@0:     
max@0:     return (info == 0);
max@0:     }
max@0:   #else
max@0:     {
max@0:     arma_ignore(eigval);
max@0:     arma_ignore(eigvec);
max@0:     arma_stop("eig_sym(): use of LAPACK needs to be enabled");
max@0:     
max@0:     return false;
max@0:     }
max@0:   #endif
max@0:   }
max@0: 
max@0: 
max@0: 
max@0: //! immediate eigenvalues and eigenvectors of a hermitian complex matrix using LAPACK
max@0: template<typename T, typename T1>
max@0: inline
max@0: bool
max@0: auxlib::eig_sym(Col<T>& eigval, Mat< std::complex<T> >& eigvec, const Base<std::complex<T>,T1>& X)
max@0:   {
max@0:   arma_extra_debug_sigprint();
max@0:   
max@0:   typedef typename std::complex<T> eT;
max@0:   
max@0:   #if defined(ARMA_USE_LAPACK)
max@0:     {
max@0:     eigvec = X.get_ref();
max@0:     
max@0:     arma_debug_check( (eigvec.is_square() == false), "eig_sym(): given matrix is not hermitian" );
max@0:     
max@0:     if(eigvec.is_empty())
max@0:       {
max@0:       eigval.reset();
max@0:       eigvec.reset();
max@0:       return true;
max@0:       }
max@0:     
max@0:     char jobz  = 'V';
max@0:     char uplo  = 'U';
max@0:     
max@0:     blas_int n_rows = eigvec.n_rows;
max@0:     blas_int lda    = eigvec.n_rows;
max@0:     blas_int lwork  = (std::max)(blas_int(1), 2*n_rows - 1);  // TODO: automatically find best size of lwork
max@0:     
max@0:     eigval.set_size( static_cast<uword>(n_rows) );
max@0:     
max@0:     podarray<eT> work( static_cast<uword>(lwork) );
max@0:     podarray<T>  rwork( static_cast<uword>((std::max)(blas_int(1), 3*n_rows - 2)) );
max@0:     
max@0:     blas_int info;
max@0:     
max@0:     arma_extra_debug_print("lapack::heev()");
max@0:     lapack::heev(&jobz, &uplo, &n_rows, eigvec.memptr(), &lda, eigval.memptr(), work.memptr(), &lwork, rwork.memptr(), &info);
max@0:     
max@0:     return (info == 0);
max@0:     }
max@0:   #else
max@0:     {
max@0:     arma_ignore(eigval);
max@0:     arma_ignore(eigvec);
max@0:     arma_ignore(X);
max@0:     arma_stop("eig_sym(): use of LAPACK needs to be enabled");
max@0:     return false;
max@0:     }
max@0:   #endif
max@0:   }
max@0: 
max@0: 
max@0: 
max@0: //! Eigenvalues and eigenvectors of a general square real matrix using LAPACK.
max@0: //! The argument 'side' specifies which eigenvectors should be calculated
max@0: //! (see code for mode details).
max@0: template<typename T, typename T1>
max@0: inline
max@0: bool
max@0: auxlib::eig_gen
max@0:   (
max@0:   Col< std::complex<T> >& eigval,
max@0:   Mat<T>& l_eigvec,
max@0:   Mat<T>& r_eigvec,
max@0:   const Base<T,T1>& X,
max@0:   const char side
max@0:   )
max@0:   {
max@0:   arma_extra_debug_sigprint();
max@0:   
max@0:   #if defined(ARMA_USE_LAPACK)
max@0:     {
max@0:     char jobvl;
max@0:     char jobvr;
max@0:     
max@0:     switch(side)
max@0:       {
max@0:       case 'l':  // left
max@0:         jobvl = 'V';
max@0:         jobvr = 'N';
max@0:         break;
max@0:         
max@0:       case 'r':  // right
max@0:         jobvl = 'N';
max@0:         jobvr = 'V';
max@0:         break;
max@0:         
max@0:       case 'b':  // both
max@0:         jobvl = 'V';
max@0:         jobvr = 'V';
max@0:         break;
max@0:         
max@0:       case 'n':  // neither
max@0:         jobvl = 'N';
max@0:         jobvr = 'N';
max@0:         break;
max@0:       
max@0:       default:
max@0:         arma_stop("eig_gen(): parameter 'side' is invalid");
max@0:         return false;
max@0:       }
max@0:     
max@0:     Mat<T> A(X.get_ref());
max@0:     arma_debug_check( (A.is_square() == false), "eig_gen(): given matrix is not square" );
max@0:     
max@0:     if(A.is_empty())
max@0:       {
max@0:       eigval.reset();
max@0:       l_eigvec.reset();
max@0:       r_eigvec.reset();
max@0:       return true;
max@0:       }
max@0:     
max@0:     uword A_n_rows = A.n_rows;
max@0:     
max@0:     blas_int n_rows = A_n_rows;
max@0:     blas_int lda    = A_n_rows;
max@0:     blas_int lwork  = (std::max)(blas_int(1), 4*n_rows);  // TODO: automatically find best size of lwork
max@0:     
max@0:     eigval.set_size(A_n_rows);
max@0:     l_eigvec.set_size(A_n_rows, A_n_rows);
max@0:     r_eigvec.set_size(A_n_rows, A_n_rows);
max@0:     
max@0:     podarray<T> work( static_cast<uword>(lwork) );
max@0:     podarray<T> rwork( static_cast<uword>((std::max)(blas_int(1), 3*n_rows)) );
max@0:     
max@0:     podarray<T> wr(A_n_rows);
max@0:     podarray<T> wi(A_n_rows);
max@0:     
max@0:     Mat<T> A_copy = A;
max@0:     blas_int info;
max@0:     
max@0:     arma_extra_debug_print("lapack::geev()");
max@0:     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);
max@0:     
max@0:     
max@0:     eigval.set_size(A_n_rows);
max@0:     for(uword i=0; i<A_n_rows; ++i)
max@0:       {
max@0:       eigval[i] = std::complex<T>(wr[i], wi[i]);
max@0:       }
max@0:     
max@0:     return (info == 0);
max@0:     }
max@0:   #else
max@0:     {
max@0:     arma_ignore(eigval);
max@0:     arma_ignore(l_eigvec);
max@0:     arma_ignore(r_eigvec);
max@0:     arma_ignore(X);
max@0:     arma_ignore(side);
max@0:     arma_stop("eig_gen(): use of LAPACK needs to be enabled");
max@0:     return false;
max@0:     }
max@0:   #endif
max@0:   }
max@0: 
max@0: 
max@0: 
max@0: 
max@0: 
max@0: //! Eigenvalues and eigenvectors of a general square complex matrix using LAPACK
max@0: //! The argument 'side' specifies which eigenvectors should be calculated
max@0: //! (see code for mode details).
max@0: template<typename T, typename T1>
max@0: inline
max@0: bool
max@0: auxlib::eig_gen
max@0:   (
max@0:   Col< std::complex<T> >& eigval,
max@0:   Mat< std::complex<T> >& l_eigvec, 
max@0:   Mat< std::complex<T> >& r_eigvec, 
max@0:   const Base< std::complex<T>, T1 >& X, 
max@0:   const char side
max@0:   )
max@0:   {
max@0:   arma_extra_debug_sigprint();
max@0:   
max@0:   typedef typename std::complex<T> eT;
max@0:   
max@0:   #if defined(ARMA_USE_LAPACK)
max@0:     {
max@0:     char jobvl;
max@0:     char jobvr;
max@0:     
max@0:     switch(side)
max@0:       {
max@0:       case 'l':  // left
max@0:         jobvl = 'V';
max@0:         jobvr = 'N';
max@0:         break;
max@0:         
max@0:       case 'r':  // right
max@0:         jobvl = 'N';
max@0:         jobvr = 'V';
max@0:         break;
max@0:         
max@0:       case 'b':  // both
max@0:         jobvl = 'V';
max@0:         jobvr = 'V';
max@0:         break;
max@0:         
max@0:       case 'n':  // neither
max@0:         jobvl = 'N';
max@0:         jobvr = 'N';
max@0:         break;
max@0:       
max@0:       default:
max@0:         arma_stop("eig_gen(): parameter 'side' is invalid");
max@0:         return false;
max@0:       }
max@0:     
max@0:     Mat<eT> A(X.get_ref());
max@0:     arma_debug_check( (A.is_square() == false), "eig_gen(): given matrix is not square" );
max@0:     
max@0:     if(A.is_empty())
max@0:       {
max@0:       eigval.reset();
max@0:       l_eigvec.reset();
max@0:       r_eigvec.reset();
max@0:       return true;
max@0:       }
max@0:     
max@0:     uword A_n_rows = A.n_rows;
max@0:        
max@0:     blas_int n_rows = A_n_rows;
max@0:     blas_int lda    = A_n_rows;
max@0:     blas_int lwork  = (std::max)(blas_int(1), 4*n_rows);  // TODO: automatically find best size of lwork
max@0:     
max@0:     eigval.set_size(A_n_rows);
max@0:     l_eigvec.set_size(A_n_rows, A_n_rows);
max@0:     r_eigvec.set_size(A_n_rows, A_n_rows);
max@0:     
max@0:     podarray<eT> work( static_cast<uword>(lwork) );
max@0:     podarray<T>  rwork( static_cast<uword>((std::max)(blas_int(1), 3*n_rows)) );  // was 2,3
max@0:     
max@0:     blas_int info;
max@0:     
max@0:     arma_extra_debug_print("lapack::cx_geev()");
max@0:     lapack::cx_geev(&jobvl, &jobvr, &n_rows, A.memptr(), &lda, eigval.memptr(), l_eigvec.memptr(), &n_rows, r_eigvec.memptr(), &n_rows, work.memptr(), &lwork, rwork.memptr(), &info);
max@0:     
max@0:     return (info == 0);
max@0:     }
max@0:   #else
max@0:     {
max@0:     arma_ignore(eigval);
max@0:     arma_ignore(l_eigvec);
max@0:     arma_ignore(r_eigvec);
max@0:     arma_ignore(X);
max@0:     arma_ignore(side);
max@0:     arma_stop("eig_gen(): use of LAPACK needs to be enabled");
max@0:     return false;
max@0:     }
max@0:   #endif
max@0:   }
max@0: 
max@0: 
max@0: 
max@0: template<typename eT, typename T1>
max@0: inline
max@0: bool
max@0: auxlib::chol(Mat<eT>& out, const Base<eT,T1>& X)
max@0:   {
max@0:   arma_extra_debug_sigprint();
max@0:   
max@0:   #if defined(ARMA_USE_LAPACK)
max@0:     {
max@0:     out = X.get_ref();
max@0:     
max@0:     arma_debug_check( (out.is_square() == false), "chol(): given matrix is not square" );
max@0:     
max@0:     if(out.is_empty())
max@0:       {
max@0:       return true;
max@0:       }
max@0:     
max@0:     const uword out_n_rows = out.n_rows;
max@0:     
max@0:     char      uplo = 'U';
max@0:     blas_int  n    = out_n_rows;
max@0:     blas_int  info;
max@0:     
max@0:     lapack::potrf(&uplo, &n, out.memptr(), &n, &info);
max@0:     
max@0:     for(uword col=0; col<out_n_rows; ++col)
max@0:       {
max@0:       eT* colptr = out.colptr(col);
max@0:       
max@0:       for(uword row=(col+1); row < out_n_rows; ++row)
max@0:         {
max@0:         colptr[row] = eT(0);
max@0:         }
max@0:       }
max@0:     
max@0:     return (info == 0);
max@0:     }
max@0:   #else
max@0:     {
max@0:     arma_ignore(out);
max@0:     arma_stop("chol(): use of LAPACK needs to be enabled");
max@0:     return false;
max@0:     }
max@0:   #endif
max@0:   }
max@0: 
max@0: 
max@0: 
max@0: template<typename eT, typename T1>
max@0: inline
max@0: bool 
max@0: auxlib::qr(Mat<eT>& Q, Mat<eT>& R, const Base<eT,T1>& X)
max@0:   {
max@0:   arma_extra_debug_sigprint();
max@0:   
max@0:   #if defined(ARMA_USE_LAPACK)
max@0:     {
max@0:     R = X.get_ref();
max@0:     
max@0:     const uword R_n_rows = R.n_rows;
max@0:     const uword R_n_cols = R.n_cols;
max@0:     
max@0:     if(R.is_empty())
max@0:       {
max@0:       Q.eye(R_n_rows, R_n_rows);
max@0:       return true;
max@0:       }
max@0:     
max@0:     blas_int m            = static_cast<blas_int>(R_n_rows);
max@0:     blas_int n            = static_cast<blas_int>(R_n_cols);
max@0:     blas_int work_len     = (std::max)(blas_int(1),n);
max@0:     blas_int work_len_tmp;
max@0:     blas_int k            = (std::min)(m,n);
max@0:     blas_int info;
max@0:     
max@0:     podarray<eT> tau( static_cast<uword>(k) );
max@0:     podarray<eT> work( static_cast<uword>(work_len) );
max@0:     
max@0:     // query for the optimum value of work_len
max@0:     work_len_tmp = -1;
max@0:     lapack::geqrf(&m, &n, R.memptr(), &m, tau.memptr(), work.memptr(), &work_len_tmp, &info);
max@0:     
max@0:     if(info == 0)
max@0:       {
max@0:       work_len = static_cast<blas_int>(access::tmp_real(work[0]));
max@0:       work.set_size( static_cast<uword>(work_len) );
max@0:       }
max@0:     
max@0:     lapack::geqrf(&m, &n, R.memptr(), &m, tau.memptr(), work.memptr(), &work_len, &info);
max@0:     
max@0:     Q.set_size(R_n_rows, R_n_rows);
max@0:     
max@0:     arrayops::copy( Q.memptr(), R.memptr(), (std::min)(Q.n_elem, R.n_elem) );
max@0:     
max@0:     //
max@0:     // construct R
max@0:     
max@0:     for(uword col=0; col < R_n_cols; ++col)
max@0:       {
max@0:       for(uword row=(col+1); row < R_n_rows; ++row)
max@0:         {
max@0:         R.at(row,col) = eT(0);
max@0:         }
max@0:       }
max@0:     
max@0:     
max@0:     if( (is_float<eT>::value == true) || (is_double<eT>::value == true) )
max@0:       {
max@0:       // query for the optimum value of work_len
max@0:       work_len_tmp = -1;
max@0:       lapack::orgqr(&m, &m, &k, Q.memptr(), &m, tau.memptr(), work.memptr(), &work_len_tmp, &info);
max@0:       
max@0:       if(info == 0)
max@0:         {
max@0:         work_len = static_cast<blas_int>(access::tmp_real(work[0]));
max@0:         work.set_size( static_cast<uword>(work_len) );
max@0:         }
max@0:       
max@0:       lapack::orgqr(&m, &m, &k, Q.memptr(), &m, tau.memptr(), work.memptr(), &work_len, &info);
max@0:       }
max@0:     else
max@0:     if( (is_supported_complex_float<eT>::value == true) || (is_supported_complex_double<eT>::value == true) )
max@0:       {
max@0:       // query for the optimum value of work_len
max@0:       work_len_tmp = -1;
max@0:       lapack::ungqr(&m, &m, &k, Q.memptr(), &m, tau.memptr(), work.memptr(), &work_len_tmp, &info);
max@0:       
max@0:       if(info == 0)
max@0:         {
max@0:         work_len = static_cast<blas_int>(access::tmp_real(work[0]));
max@0:         work.set_size( static_cast<uword>(work_len) );
max@0:         }
max@0:       
max@0:       lapack::ungqr(&m, &m, &k, Q.memptr(), &m, tau.memptr(), work.memptr(), &work_len, &info);
max@0:       }
max@0:     
max@0:     return (info == 0);
max@0:     }
max@0:   #else
max@0:     {
max@0:     arma_ignore(Q);
max@0:     arma_ignore(R);
max@0:     arma_ignore(X);
max@0:     arma_stop("qr(): use of LAPACK needs to be enabled");
max@0:     return false;
max@0:     }
max@0:   #endif
max@0:   }
max@0: 
max@0: 
max@0: 
max@0: template<typename eT, typename T1>
max@0: inline
max@0: bool
max@0: auxlib::svd(Col<eT>& S, const Base<eT,T1>& X, uword& X_n_rows, uword& X_n_cols)
max@0:   {
max@0:   arma_extra_debug_sigprint();
max@0:   
max@0:   #if defined(ARMA_USE_LAPACK)
max@0:     {
max@0:     Mat<eT> A(X.get_ref());
max@0:     
max@0:     X_n_rows = A.n_rows;
max@0:     X_n_cols = A.n_cols;
max@0:     
max@0:     if(A.is_empty())
max@0:       {
max@0:       S.reset();
max@0:       return true;
max@0:       }
max@0:     
max@0:     Mat<eT> U(1, 1);
max@0:     Mat<eT> V(1, A.n_cols);
max@0:     
max@0:     char jobu  = 'N';
max@0:     char jobvt = 'N';
max@0:     
max@0:     blas_int  m     = A.n_rows;
max@0:     blas_int  n     = A.n_cols;
max@0:     blas_int  lda   = A.n_rows;
max@0:     blas_int  ldu   = U.n_rows;
max@0:     blas_int  ldvt  = V.n_rows;
max@0:     blas_int  lwork = 2 * (std::max)(blas_int(1), (std::max)( (3*(std::min)(m,n) + (std::max)(m,n)), 5*(std::min)(m,n) ) );
max@0:     blas_int  info;
max@0:     
max@0:     S.set_size( static_cast<uword>((std::min)(m, n)) );
max@0:     
max@0:     podarray<eT> work( static_cast<uword>(lwork) );
max@0:   
max@0:   
max@0:     // let gesvd_() calculate the optimum size of the workspace
max@0:     blas_int lwork_tmp = -1;
max@0:     
max@0:     lapack::gesvd<eT>
max@0:       (
max@0:       &jobu, &jobvt,
max@0:       &m,&n,
max@0:       A.memptr(), &lda,
max@0:       S.memptr(),
max@0:       U.memptr(), &ldu,
max@0:       V.memptr(), &ldvt,
max@0:       work.memptr(), &lwork_tmp,
max@0:       &info
max@0:       );
max@0:     
max@0:     if(info == 0)
max@0:       {
max@0:       blas_int proposed_lwork = static_cast<blas_int>(work[0]);
max@0:       
max@0:       if(proposed_lwork > lwork)
max@0:         {
max@0:         lwork = proposed_lwork;
max@0:         work.set_size( static_cast<uword>(lwork) );
max@0:         }
max@0:       
max@0:       lapack::gesvd<eT>
max@0:         (
max@0:         &jobu, &jobvt,
max@0:         &m, &n,
max@0:         A.memptr(), &lda,
max@0:         S.memptr(),
max@0:         U.memptr(), &ldu,
max@0:         V.memptr(), &ldvt,
max@0:         work.memptr(), &lwork,
max@0:         &info
max@0:         );
max@0:       }
max@0:     
max@0:     return (info == 0);
max@0:     }
max@0:   #else
max@0:     {
max@0:     arma_ignore(S);
max@0:     arma_ignore(X);
max@0:     arma_ignore(X_n_rows);
max@0:     arma_ignore(X_n_cols);
max@0:     arma_stop("svd(): use of LAPACK needs to be enabled");
max@0:     return false;
max@0:     }
max@0:   #endif
max@0:   }
max@0: 
max@0: 
max@0: 
max@0: template<typename T, typename T1>
max@0: inline
max@0: bool
max@0: auxlib::svd(Col<T>& S, const Base<std::complex<T>, T1>& X, uword& X_n_rows, uword& X_n_cols)
max@0:   {
max@0:   arma_extra_debug_sigprint();
max@0:   
max@0:   typedef std::complex<T> eT;
max@0:   
max@0:   #if defined(ARMA_USE_LAPACK)
max@0:     {
max@0:     Mat<eT> A(X.get_ref());
max@0:     
max@0:     X_n_rows = A.n_rows;
max@0:     X_n_cols = A.n_cols;
max@0:     
max@0:     if(A.is_empty())
max@0:       {
max@0:       S.reset();
max@0:       return true;
max@0:       }
max@0:     
max@0:     Mat<eT> U(1, 1);
max@0:     Mat<eT> V(1, A.n_cols);
max@0:     
max@0:     char jobu  = 'N';
max@0:     char jobvt = 'N';
max@0:     
max@0:     blas_int  m     = A.n_rows;
max@0:     blas_int  n     = A.n_cols;
max@0:     blas_int  lda   = A.n_rows;
max@0:     blas_int  ldu   = U.n_rows;
max@0:     blas_int  ldvt  = V.n_rows;
max@0:     blas_int  lwork = 2 * (std::max)(blas_int(1), 2*(std::min)(m,n)+(std::max)(m,n) );
max@0:     blas_int  info;
max@0:     
max@0:     S.set_size( static_cast<uword>((std::min)(m,n)) );
max@0:     
max@0:     podarray<eT> work( static_cast<uword>(lwork) );
max@0:     podarray<T>  rwork( static_cast<uword>(5*(std::min)(m,n)) );
max@0:     
max@0:     // let gesvd_() calculate the optimum size of the workspace
max@0:     blas_int lwork_tmp = -1;
max@0:     
max@0:     lapack::cx_gesvd<T>
max@0:       (
max@0:       &jobu, &jobvt,
max@0:       &m, &n,
max@0:       A.memptr(), &lda,
max@0:       S.memptr(),
max@0:       U.memptr(), &ldu,
max@0:       V.memptr(), &ldvt,
max@0:       work.memptr(), &lwork_tmp,
max@0:       rwork.memptr(),
max@0:       &info
max@0:       );
max@0:     
max@0:     if(info == 0)
max@0:       {
max@0:       blas_int proposed_lwork = static_cast<blas_int>(real(work[0]));
max@0:       if(proposed_lwork > lwork)
max@0:         {
max@0:         lwork = proposed_lwork;
max@0:         work.set_size( static_cast<uword>(lwork) );
max@0:         }
max@0:       
max@0:       lapack::cx_gesvd<T>
max@0:         (
max@0:         &jobu, &jobvt,
max@0:         &m, &n,
max@0:         A.memptr(), &lda,
max@0:         S.memptr(),
max@0:         U.memptr(), &ldu,
max@0:         V.memptr(), &ldvt,
max@0:         work.memptr(), &lwork,
max@0:         rwork.memptr(),
max@0:         &info
max@0:         );
max@0:       }
max@0:         
max@0:     return (info == 0);
max@0:     }
max@0:   #else
max@0:     {
max@0:     arma_ignore(S);
max@0:     arma_ignore(X);
max@0:     arma_ignore(X_n_rows);
max@0:     arma_ignore(X_n_cols);
max@0: 
max@0:     arma_stop("svd(): use of LAPACK needs to be enabled");
max@0:     return false;
max@0:     }
max@0:   #endif
max@0:   }
max@0: 
max@0: 
max@0: 
max@0: template<typename eT, typename T1>
max@0: inline
max@0: bool
max@0: auxlib::svd(Col<eT>& S, const Base<eT,T1>& X)
max@0:   {
max@0:   arma_extra_debug_sigprint();
max@0:   
max@0:   uword junk;
max@0:   return auxlib::svd(S, X, junk, junk);
max@0:   }
max@0: 
max@0: 
max@0: 
max@0: template<typename T, typename T1>
max@0: inline
max@0: bool
max@0: auxlib::svd(Col<T>& S, const Base<std::complex<T>, T1>& X)
max@0:   {
max@0:   arma_extra_debug_sigprint();
max@0:   
max@0:   uword junk;
max@0:   return auxlib::svd(S, X, junk, junk);
max@0:   }
max@0: 
max@0: 
max@0: 
max@0: template<typename eT, typename T1>
max@0: inline
max@0: bool
max@0: auxlib::svd(Mat<eT>& U, Col<eT>& S, Mat<eT>& V, const Base<eT,T1>& X)
max@0:   {
max@0:   arma_extra_debug_sigprint();
max@0:   
max@0:   #if defined(ARMA_USE_LAPACK)
max@0:     {
max@0:     Mat<eT> A(X.get_ref());
max@0:     
max@0:     if(A.is_empty())
max@0:       {
max@0:       U.eye(A.n_rows, A.n_rows);
max@0:       S.reset();
max@0:       V.eye(A.n_cols, A.n_cols);
max@0:       return true;
max@0:       }
max@0:     
max@0:     U.set_size(A.n_rows, A.n_rows);
max@0:     V.set_size(A.n_cols, A.n_cols);
max@0:     
max@0:     char jobu  = 'A';
max@0:     char jobvt = 'A';
max@0:     
max@0:     blas_int  m     = A.n_rows;
max@0:     blas_int  n     = A.n_cols;
max@0:     blas_int  lda   = A.n_rows;
max@0:     blas_int  ldu   = U.n_rows;
max@0:     blas_int  ldvt  = V.n_rows;
max@0:     blas_int  lwork = 2 * (std::max)(blas_int(1), (std::max)( (3*(std::min)(m,n) + (std::max)(m,n)), 5*(std::min)(m,n) ) );
max@0:     blas_int  info;
max@0:     
max@0:     
max@0:     S.set_size( static_cast<uword>((std::min)(m,n)) );
max@0:     podarray<eT> work( static_cast<uword>(lwork) );
max@0:   
max@0:     // let gesvd_() calculate the optimum size of the workspace
max@0:     blas_int lwork_tmp = -1;
max@0:     
max@0:     lapack::gesvd<eT>
max@0:       (
max@0:       &jobu, &jobvt,
max@0:       &m, &n,
max@0:       A.memptr(), &lda,
max@0:       S.memptr(),
max@0:       U.memptr(), &ldu,
max@0:       V.memptr(), &ldvt,
max@0:       work.memptr(), &lwork_tmp,
max@0:       &info
max@0:       );
max@0:     
max@0:     if(info == 0)
max@0:       {
max@0:       blas_int proposed_lwork = static_cast<blas_int>(work[0]);
max@0:       if(proposed_lwork > lwork)
max@0:         {
max@0:         lwork = proposed_lwork;
max@0:         work.set_size( static_cast<uword>(lwork) );
max@0:         }
max@0:       
max@0:       lapack::gesvd<eT>
max@0:         (
max@0:         &jobu, &jobvt,
max@0:         &m, &n,
max@0:         A.memptr(), &lda,
max@0:         S.memptr(),
max@0:         U.memptr(), &ldu,
max@0:         V.memptr(), &ldvt,
max@0:         work.memptr(), &lwork,
max@0:         &info
max@0:         );
max@0:       
max@0:       op_strans::apply(V,V);  // op_strans will work out that an in-place transpose can be done
max@0:       }
max@0:     
max@0:     return (info == 0);
max@0:     }
max@0:   #else
max@0:     {
max@0:     arma_ignore(U);
max@0:     arma_ignore(S);
max@0:     arma_ignore(V);
max@0:     arma_ignore(X);
max@0:     arma_stop("svd(): use of LAPACK needs to be enabled");
max@0:     return false;
max@0:     }
max@0:   #endif
max@0:   }
max@0: 
max@0: 
max@0: 
max@0: template<typename T, typename T1>
max@0: inline
max@0: bool
max@0: auxlib::svd(Mat< std::complex<T> >& U, Col<T>& S, Mat< std::complex<T> >& V, const Base< std::complex<T>, T1>& X)
max@0:   {
max@0:   arma_extra_debug_sigprint();
max@0:   
max@0:   typedef std::complex<T> eT;
max@0:   
max@0:   #if defined(ARMA_USE_LAPACK)
max@0:     {
max@0:     Mat<eT> A(X.get_ref());
max@0:     
max@0:     if(A.is_empty())
max@0:       {
max@0:       U.eye(A.n_rows, A.n_rows);
max@0:       S.reset();
max@0:       V.eye(A.n_cols, A.n_cols);
max@0:       return true;
max@0:       }
max@0:     
max@0:     U.set_size(A.n_rows, A.n_rows);
max@0:     V.set_size(A.n_cols, A.n_cols);
max@0:     
max@0:     char jobu  = 'A';
max@0:     char jobvt = 'A';
max@0:     
max@0:     blas_int  m     = A.n_rows;
max@0:     blas_int  n     = A.n_cols;
max@0:     blas_int  lda   = A.n_rows;
max@0:     blas_int  ldu   = U.n_rows;
max@0:     blas_int  ldvt  = V.n_rows;
max@0:     blas_int  lwork = 2 * (std::max)(blas_int(1), 2*(std::min)(m,n)+(std::max)(m,n) );
max@0:     blas_int  info;
max@0:     
max@0:     S.set_size( static_cast<uword>((std::min)(m,n)) );
max@0:     
max@0:     podarray<eT> work( static_cast<uword>(lwork) );
max@0:     podarray<T>  rwork( static_cast<uword>(5*(std::min)(m,n)) );
max@0:     
max@0:     // let gesvd_() calculate the optimum size of the workspace
max@0:     blas_int lwork_tmp = -1;
max@0:     lapack::cx_gesvd<T>
max@0:      (
max@0:      &jobu, &jobvt,
max@0:      &m, &n,
max@0:      A.memptr(), &lda,
max@0:      S.memptr(),
max@0:      U.memptr(), &ldu,
max@0:      V.memptr(), &ldvt,
max@0:      work.memptr(), &lwork_tmp,
max@0:      rwork.memptr(),
max@0:      &info
max@0:      );
max@0:     
max@0:     if(info == 0)
max@0:       {
max@0:       blas_int proposed_lwork = static_cast<blas_int>(real(work[0]));
max@0:       if(proposed_lwork > lwork)
max@0:         {
max@0:         lwork = proposed_lwork;
max@0:         work.set_size( static_cast<uword>(lwork) );
max@0:         }
max@0:       
max@0:       lapack::cx_gesvd<T>
max@0:         (
max@0:         &jobu, &jobvt,
max@0:         &m, &n,
max@0:         A.memptr(), &lda,
max@0:         S.memptr(),
max@0:         U.memptr(), &ldu,
max@0:         V.memptr(), &ldvt,
max@0:         work.memptr(), &lwork,
max@0:         rwork.memptr(),
max@0:         &info
max@0:         );
max@0:       
max@0:       op_htrans::apply(V,V);  // op_htrans will work out that an in-place transpose can be done
max@0:       }
max@0:     
max@0:     return (info == 0);
max@0:     }
max@0:   #else
max@0:     {
max@0:     arma_ignore(U);
max@0:     arma_ignore(S);
max@0:     arma_ignore(V);
max@0:     arma_ignore(X);
max@0:     arma_stop("svd(): use of LAPACK needs to be enabled");
max@0:     return false;
max@0:     }
max@0:   #endif
max@0:   
max@0:   }
max@0: 
max@0: 
max@0: 
max@0: template<typename eT, typename T1>
max@0: inline
max@0: bool
max@0: auxlib::svd_econ(Mat<eT>& U, Col<eT>& S, Mat<eT>& V, const Base<eT,T1>& X, const char mode)
max@0:   {
max@0:   arma_extra_debug_sigprint();
max@0:   
max@0:   #if defined(ARMA_USE_LAPACK)
max@0:     {
max@0:     Mat<eT> A(X.get_ref());
max@0:     
max@0:     blas_int m    = A.n_rows;
max@0:     blas_int n    = A.n_cols;
max@0:     blas_int lda  = A.n_rows;
max@0:     
max@0:     S.set_size( static_cast<uword>((std::min)(m,n)) );
max@0:     
max@0:     blas_int ldu  = 0;
max@0:     blas_int ldvt = 0;
max@0:     
max@0:     char jobu;
max@0:     char jobvt;
max@0:     
max@0:     switch(mode)
max@0:       {
max@0:       case 'l':
max@0:         jobu  = 'S';
max@0:         jobvt = 'N';
max@0:         
max@0:         ldu  = m;
max@0:         ldvt = 1;
max@0:         
max@0:         U.set_size( static_cast<uword>(ldu), static_cast<uword>((std::min)(m,n)) );
max@0:         V.reset();
max@0:         
max@0:         break;
max@0:       
max@0:       
max@0:       case 'r':
max@0:         jobu  = 'N';
max@0:         jobvt = 'S';
max@0:         
max@0:         ldu = 1;
max@0:         ldvt = (std::min)(m,n);
max@0:         
max@0:         U.reset();
max@0:         V.set_size( static_cast<uword>(ldvt), static_cast<uword>(n) );
max@0:         
max@0:         break;
max@0:       
max@0:       
max@0:       case 'b':
max@0:         jobu  = 'S';
max@0:         jobvt = 'S';
max@0:         
max@0:         ldu  = m;
max@0:         ldvt = (std::min)(m,n);
max@0:         
max@0:         U.set_size( static_cast<uword>(ldu),  static_cast<uword>((std::min)(m,n)) );
max@0:         V.set_size( static_cast<uword>(ldvt), static_cast<uword>(n)               );
max@0:         
max@0:         break;
max@0:       
max@0:       
max@0:       default:
max@0:         U.reset();
max@0:         S.reset();
max@0:         V.reset();
max@0:         return false;
max@0:       }
max@0:     
max@0:     
max@0:     if(A.is_empty())
max@0:       {
max@0:       U.eye();
max@0:       S.reset();
max@0:       V.eye();
max@0:       return true;
max@0:       }
max@0:     
max@0:     
max@0:     blas_int lwork = 2 * (std::max)(blas_int(1), (std::max)( (3*(std::min)(m,n) + (std::max)(m,n)), 5*(std::min)(m,n) ) );
max@0:     blas_int info  = 0;
max@0:     
max@0:     
max@0:     podarray<eT> work( static_cast<uword>(lwork) );
max@0:     
max@0:     // let gesvd_() calculate the optimum size of the workspace
max@0:     blas_int lwork_tmp = -1;
max@0:     
max@0:     lapack::gesvd<eT>
max@0:       (
max@0:       &jobu, &jobvt,
max@0:       &m, &n,
max@0:       A.memptr(), &lda,
max@0:       S.memptr(),
max@0:       U.memptr(), &ldu,
max@0:       V.memptr(), &ldvt,
max@0:       work.memptr(), &lwork_tmp,
max@0:       &info
max@0:       );
max@0:     
max@0:     if(info == 0)
max@0:       {
max@0:       blas_int proposed_lwork = static_cast<blas_int>(work[0]);
max@0:       if(proposed_lwork > lwork)
max@0:         {
max@0:         lwork = proposed_lwork;
max@0:         work.set_size( static_cast<uword>(lwork) );
max@0:         }
max@0:       
max@0:       lapack::gesvd<eT>
max@0:         (
max@0:         &jobu, &jobvt,
max@0:         &m, &n,
max@0:         A.memptr(), &lda,
max@0:         S.memptr(),
max@0:         U.memptr(), &ldu,
max@0:         V.memptr(), &ldvt,
max@0:         work.memptr(), &lwork,
max@0:         &info
max@0:         );
max@0:       
max@0:       op_strans::apply(V,V);  // op_strans will work out that an in-place transpose can be done
max@0:       }
max@0:     
max@0:     return (info == 0);
max@0:     }
max@0:   #else
max@0:     {
max@0:     arma_ignore(U);
max@0:     arma_ignore(S);
max@0:     arma_ignore(V);
max@0:     arma_ignore(X);
max@0:     arma_ignore(mode);
max@0:     arma_stop("svd(): use of LAPACK needs to be enabled");
max@0:     return false;
max@0:     }
max@0:   #endif
max@0:   }
max@0: 
max@0: 
max@0: 
max@0: template<typename T, typename T1>
max@0: inline
max@0: bool
max@0: auxlib::svd_econ(Mat< std::complex<T> >& U, Col<T>& S, Mat< std::complex<T> >& V, const Base< std::complex<T>, T1>& X, const char mode)
max@0:   {
max@0:   arma_extra_debug_sigprint();
max@0:   
max@0:   typedef std::complex<T> eT;
max@0:   
max@0:   #if defined(ARMA_USE_LAPACK)
max@0:     {
max@0:     Mat<eT> A(X.get_ref());
max@0:     
max@0:     blas_int m    = A.n_rows;
max@0:     blas_int n    = A.n_cols;
max@0:     blas_int lda  = A.n_rows;
max@0:     
max@0:     S.set_size( static_cast<uword>((std::min)(m,n)) );
max@0:     
max@0:     blas_int ldu  = 0;
max@0:     blas_int ldvt = 0;
max@0:     
max@0:     char jobu;
max@0:     char jobvt;
max@0:     
max@0:     switch(mode)
max@0:       {
max@0:       case 'l':
max@0:         jobu  = 'S';
max@0:         jobvt = 'N';
max@0:         
max@0:         ldu  = m;
max@0:         ldvt = 1;
max@0:         
max@0:         U.set_size( static_cast<uword>(ldu), static_cast<uword>((std::min)(m,n)) );
max@0:         V.reset();
max@0:         
max@0:         break;
max@0:       
max@0:       
max@0:       case 'r':
max@0:         jobu  = 'N';
max@0:         jobvt = 'S';
max@0:         
max@0:         ldu  = 1;
max@0:         ldvt = (std::min)(m,n);
max@0:         
max@0:         U.reset();
max@0:         V.set_size( static_cast<uword>(ldvt), static_cast<uword>(n) );
max@0:         
max@0:         break;
max@0:       
max@0:       
max@0:       case 'b':
max@0:         jobu  = 'S';
max@0:         jobvt = 'S';
max@0:         
max@0:         ldu  = m;
max@0:         ldvt = (std::min)(m,n);
max@0:         
max@0:         U.set_size( static_cast<uword>(ldu),  static_cast<uword>((std::min)(m,n)) );
max@0:         V.set_size( static_cast<uword>(ldvt), static_cast<uword>(n)               );
max@0:         
max@0:         break;
max@0:       
max@0:       
max@0:       default:
max@0:         U.reset();
max@0:         S.reset();
max@0:         V.reset();
max@0:         return false;
max@0:       }
max@0:     
max@0:     
max@0:     if(A.is_empty())
max@0:       {
max@0:       U.eye();
max@0:       S.reset();
max@0:       V.eye();
max@0:       return true;
max@0:       }
max@0:     
max@0:     
max@0:     blas_int lwork = 2 * (std::max)(blas_int(1), (std::max)( (3*(std::min)(m,n) + (std::max)(m,n)), 5*(std::min)(m,n) ) );
max@0:     blas_int info  = 0;
max@0:     
max@0:     
max@0:     podarray<eT>  work( static_cast<uword>(lwork) );
max@0:     podarray<T>  rwork( static_cast<uword>(5*(std::min)(m,n)) );
max@0:     
max@0:     // let gesvd_() calculate the optimum size of the workspace
max@0:     blas_int lwork_tmp = -1;
max@0:     
max@0:     lapack::cx_gesvd<T>
max@0:       (
max@0:       &jobu, &jobvt,
max@0:       &m, &n,
max@0:       A.memptr(), &lda,
max@0:       S.memptr(),
max@0:       U.memptr(), &ldu,
max@0:       V.memptr(), &ldvt,
max@0:       work.memptr(), &lwork_tmp,
max@0:       rwork.memptr(),
max@0:       &info
max@0:       );
max@0:     
max@0:     if(info == 0)
max@0:       {
max@0:       blas_int proposed_lwork = static_cast<blas_int>(real(work[0]));
max@0:       if(proposed_lwork > lwork)
max@0:         {
max@0:         lwork = proposed_lwork;
max@0:         work.set_size( static_cast<uword>(lwork) );
max@0:         }
max@0:       
max@0:       lapack::cx_gesvd<T>
max@0:         (
max@0:         &jobu, &jobvt,
max@0:         &m, &n,
max@0:         A.memptr(), &lda,
max@0:         S.memptr(),
max@0:         U.memptr(), &ldu,
max@0:         V.memptr(), &ldvt,
max@0:         work.memptr(), &lwork,
max@0:         rwork.memptr(),
max@0:         &info
max@0:         );
max@0:       
max@0:       op_htrans::apply(V,V);  // op_strans will work out that an in-place transpose can be done
max@0:       }
max@0:     
max@0:     return (info == 0);
max@0:     }
max@0:   #else
max@0:     {
max@0:     arma_ignore(U);
max@0:     arma_ignore(S);
max@0:     arma_ignore(V);
max@0:     arma_ignore(X);
max@0:     arma_ignore(mode);
max@0:     arma_stop("svd(): use of LAPACK needs to be enabled");
max@0:     return false;
max@0:     }
max@0:   #endif
max@0:   }
max@0: 
max@0: 
max@0: 
max@0: //! Solve a system of linear equations.
max@0: //! Assumes that A.n_rows = A.n_cols and B.n_rows = A.n_rows
max@0: template<typename eT>
max@0: inline
max@0: bool
max@0: auxlib::solve(Mat<eT>& out, Mat<eT>& A, const Mat<eT>& B, const bool slow)
max@0:   {
max@0:   arma_extra_debug_sigprint();
max@0:   
max@0:   if(A.is_empty() || B.is_empty())
max@0:     {
max@0:     out.zeros(A.n_cols, B.n_cols);
max@0:     return true;
max@0:     }
max@0:   else
max@0:     {
max@0:     const uword A_n_rows = A.n_rows;
max@0:     
max@0:     bool status = false;
max@0:     
max@0:     if( (A_n_rows <= 4) && (slow == false) )
max@0:       {
max@0:       Mat<eT> A_inv;
max@0:       
max@0:       status = auxlib::inv_noalias_tinymat(A_inv, A, A_n_rows);
max@0:       
max@0:       if(status == true)
max@0:         {
max@0:         out.set_size(A_n_rows, B.n_cols);
max@0:         
max@0:         gemm_emul<false,false,false,false>::apply(out, A_inv, B);
max@0:         
max@0:         return true;
max@0:         }
max@0:       }
max@0:     
max@0:     if( (A_n_rows > 4) || (status == false) )
max@0:       {
max@0:       #if defined(ARMA_USE_ATLAS)
max@0:         {
max@0:         podarray<int> ipiv(A_n_rows);
max@0:         
max@0:         out = B;
max@0:         
max@0:         int info = atlas::clapack_gesv<eT>(atlas::CblasColMajor, A_n_rows, B.n_cols, A.memptr(), A_n_rows, ipiv.memptr(), out.memptr(), A_n_rows);
max@0:         
max@0:         return (info == 0);
max@0:         }
max@0:       #elif defined(ARMA_USE_LAPACK)
max@0:         {
max@0:         blas_int n    = A_n_rows;
max@0:         blas_int lda  = A_n_rows;
max@0:         blas_int ldb  = A_n_rows;
max@0:         blas_int nrhs = B.n_cols;
max@0:         blas_int info;
max@0:         
max@0:         podarray<blas_int> ipiv(A_n_rows);
max@0:         
max@0:         out = B;
max@0:         
max@0:         lapack::gesv<eT>(&n, &nrhs, A.memptr(), &lda, ipiv.memptr(), out.memptr(), &ldb, &info);
max@0:         
max@0:         return (info == 0);
max@0:         }
max@0:       #else
max@0:         {
max@0:         arma_stop("solve(): use of ATLAS or LAPACK needs to be enabled");
max@0:         return false;
max@0:         }
max@0:       #endif
max@0:       }
max@0:     }
max@0:   
max@0:   return true;
max@0:   }
max@0: 
max@0: 
max@0: 
max@0: //! Solve an over-determined system.
max@0: //! Assumes that A.n_rows > A.n_cols and B.n_rows = A.n_rows
max@0: template<typename eT>
max@0: inline
max@0: bool
max@0: auxlib::solve_od(Mat<eT>& out, Mat<eT>& A, const Mat<eT>& B)
max@0:   {
max@0:   arma_extra_debug_sigprint();
max@0:   
max@0:   #if defined(ARMA_USE_LAPACK)
max@0:     {
max@0:     if(A.is_empty() || B.is_empty())
max@0:       {
max@0:       out.zeros(A.n_cols, B.n_cols);
max@0:       return true;
max@0:       }
max@0:     
max@0:     char trans = 'N';
max@0:     
max@0:     blas_int  m     = A.n_rows;
max@0:     blas_int  n     = A.n_cols;
max@0:     blas_int  lda   = A.n_rows;
max@0:     blas_int  ldb   = A.n_rows;
max@0:     blas_int  nrhs  = B.n_cols;
max@0:     blas_int  lwork = n + (std::max)(n, nrhs);
max@0:     blas_int  info;
max@0:     
max@0:     Mat<eT> tmp = B;
max@0:     
max@0:     podarray<eT> work( static_cast<uword>(lwork) );
max@0:     
max@0:     arma_extra_debug_print("lapack::gels()");
max@0:     
max@0:     // NOTE: the dgels() function in the lapack library supplied by ATLAS 3.6 seems to have problems
max@0:     
max@0:     lapack::gels<eT>
max@0:       (
max@0:       &trans, &m, &n, &nrhs,
max@0:       A.memptr(), &lda,
max@0:       tmp.memptr(), &ldb,
max@0:       work.memptr(), &lwork,
max@0:       &info
max@0:       );
max@0:     
max@0:     arma_extra_debug_print("lapack::gels() -- finished");
max@0:     
max@0:     out.set_size(A.n_cols, B.n_cols);
max@0:     
max@0:     for(uword col=0; col<B.n_cols; ++col)
max@0:       {
max@0:       arrayops::copy( out.colptr(col), tmp.colptr(col), A.n_cols );
max@0:       }
max@0:     
max@0:     return (info == 0);
max@0:     }
max@0:   #else
max@0:     {
max@0:     arma_ignore(out);
max@0:     arma_ignore(A);
max@0:     arma_ignore(B);
max@0:     arma_stop("solve(): use of LAPACK needs to be enabled");
max@0:     return false;
max@0:     }
max@0:   #endif
max@0:   }
max@0: 
max@0: 
max@0: 
max@0: //! Solve an under-determined system.
max@0: //! Assumes that A.n_rows < A.n_cols and B.n_rows = A.n_rows
max@0: template<typename eT>
max@0: inline
max@0: bool
max@0: auxlib::solve_ud(Mat<eT>& out, Mat<eT>& A, const Mat<eT>& B)
max@0:   {
max@0:   arma_extra_debug_sigprint();
max@0:   
max@0:   #if defined(ARMA_USE_LAPACK)
max@0:     {
max@0:     if(A.is_empty() || B.is_empty())
max@0:       {
max@0:       out.zeros(A.n_cols, B.n_cols);
max@0:       return true;
max@0:       }
max@0:     
max@0:     char trans = 'N';
max@0:     
max@0:     blas_int  m     = A.n_rows;
max@0:     blas_int  n     = A.n_cols;
max@0:     blas_int  lda   = A.n_rows;
max@0:     blas_int  ldb   = A.n_cols;
max@0:     blas_int  nrhs  = B.n_cols;
max@0:     blas_int  lwork = m + (std::max)(m,nrhs);
max@0:     blas_int  info;
max@0:     
max@0:     
max@0:     Mat<eT> tmp;
max@0:     tmp.zeros(A.n_cols, B.n_cols);
max@0:     
max@0:     for(uword col=0; col<B.n_cols; ++col)
max@0:       {
max@0:       eT* tmp_colmem = tmp.colptr(col);
max@0:       
max@0:       arrayops::copy( tmp_colmem, B.colptr(col), B.n_rows );
max@0:       
max@0:       for(uword row=B.n_rows; row<A.n_cols; ++row)
max@0:         {
max@0:         tmp_colmem[row] = eT(0);
max@0:         }
max@0:       }
max@0:     
max@0:     podarray<eT> work( static_cast<uword>(lwork) );
max@0:     
max@0:     arma_extra_debug_print("lapack::gels()");
max@0:     
max@0:     // NOTE: the dgels() function in the lapack library supplied by ATLAS 3.6 seems to have problems
max@0:     
max@0:     lapack::gels<eT>
max@0:       (
max@0:       &trans, &m, &n, &nrhs,
max@0:       A.memptr(), &lda,
max@0:       tmp.memptr(), &ldb,
max@0:       work.memptr(), &lwork,
max@0:       &info
max@0:       );
max@0:     
max@0:     arma_extra_debug_print("lapack::gels() -- finished");
max@0:     
max@0:     out.set_size(A.n_cols, B.n_cols);
max@0:     
max@0:     for(uword col=0; col<B.n_cols; ++col)
max@0:       {
max@0:       arrayops::copy( out.colptr(col), tmp.colptr(col), A.n_cols );
max@0:       }
max@0:     
max@0:     return (info == 0);
max@0:     }
max@0:   #else
max@0:     {
max@0:     arma_ignore(out);
max@0:     arma_ignore(A);
max@0:     arma_ignore(B);
max@0:     arma_stop("solve(): use of LAPACK needs to be enabled");
max@0:     return false;
max@0:     }
max@0:   #endif
max@0:   }
max@0: 
max@0: 
max@0: 
max@0: //
max@0: // solve_tr
max@0: 
max@0: template<typename eT>
max@0: inline
max@0: bool
max@0: auxlib::solve_tr(Mat<eT>& out, const Mat<eT>& A, const Mat<eT>& B, const uword layout)
max@0:   {
max@0:   arma_extra_debug_sigprint();
max@0:   
max@0:   #if defined(ARMA_USE_LAPACK)
max@0:     {
max@0:     if(A.is_empty() || B.is_empty())
max@0:       {
max@0:       out.zeros(A.n_cols, B.n_cols);
max@0:       return true;
max@0:       }
max@0:     
max@0:     out = B;
max@0:     
max@0:     char     uplo  = (layout == 0) ? 'U' : 'L';
max@0:     char     trans = 'N';
max@0:     char     diag  = 'N';
max@0:     blas_int n     = blas_int(A.n_rows);
max@0:     blas_int nrhs  = blas_int(B.n_cols);
max@0:     blas_int info  = 0;
max@0:     
max@0:     lapack::trtrs<eT>(&uplo, &trans, &diag, &n, &nrhs, A.memptr(), &n, out.memptr(), &n, &info);
max@0:     
max@0:     return (info == 0);
max@0:     }
max@0:   #else
max@0:     {
max@0:     arma_ignore(out);
max@0:     arma_ignore(A);
max@0:     arma_ignore(B);
max@0:     arma_ignore(layout);
max@0:     arma_stop("solve(): use of LAPACK needs to be enabled");
max@0:     return false;
max@0:     }
max@0:   #endif
max@0:   }
max@0: 
max@0: 
max@0: 
max@0: //
max@0: // Schur decomposition
max@0: 
max@0: template<typename eT>
max@0: inline
max@0: bool
max@0: auxlib::schur_dec(Mat<eT>& Z, Mat<eT>& T, const Mat<eT>& A)
max@0:   {
max@0:   arma_extra_debug_sigprint();
max@0:   
max@0:   #if defined(ARMA_USE_LAPACK)
max@0:     {
max@0:     arma_debug_check( (A.is_square() == false), "schur_dec(): given matrix is not square" );
max@0:     
max@0:     if(A.is_empty())
max@0:       {
max@0:       Z.reset();
max@0:       T.reset();
max@0:       return true;
max@0:       }
max@0:     
max@0:     const uword A_n_rows = A.n_rows;
max@0:     
max@0:     char    jobvs    = 'V';                // get Schur vectors (Z)
max@0:     char     sort    = 'N';                // do not sort eigenvalues/vectors
max@0:     blas_int* select = 0;                  // pointer to sorting function
max@0:     blas_int    n    = blas_int(A_n_rows);
max@0:     blas_int sdim    = 0;                  // output for sorting
max@0:     
max@0:     blas_int lwork = 3 * n;                // workspace must be at least 3 * n (if set to -1, optimal size is output in work(0) and nothing else is done
max@0:     
max@0:     podarray<eT>       work( static_cast<uword>(lwork) );
max@0:     podarray<blas_int> bwork(A_n_rows);
max@0:     
max@0:     blas_int info = 0;
max@0:     
max@0:     Z.set_size(A_n_rows, A_n_rows);
max@0:     T = A;
max@0:     
max@0:     podarray<eT> wr(A_n_rows);             // output for eigenvalues
max@0:     podarray<eT> wi(A_n_rows);             // output for eigenvalues
max@0:     
max@0:     lapack::gees(&jobvs, &sort, select, &n, T.memptr(), &n, &sdim, wr.memptr(), wi.memptr(), Z.memptr(), &n, work.memptr(), &lwork, bwork.memptr(), &info);
max@0:     
max@0:     return (info == 0);
max@0:     }
max@0:   #else
max@0:     {
max@0:     arma_ignore(Z);
max@0:     arma_ignore(T);
max@0:     arma_stop("schur_dec(): use of LAPACK needs to be enabled");
max@0:     return false;
max@0:     }
max@0:   #endif
max@0:   }
max@0: 
max@0: 
max@0: 
max@0: template<typename cT>
max@0: inline
max@0: bool
max@0: auxlib::schur_dec(Mat<std::complex<cT> >& Z, Mat<std::complex<cT> >& T, const Mat<std::complex<cT> >& A)
max@0:   {
max@0:   arma_extra_debug_sigprint();
max@0:   
max@0:   #if defined(ARMA_USE_LAPACK)
max@0:     {
max@0:     arma_debug_check( (A.is_square() == false), "schur_dec(): matrix A is not square" );
max@0:     
max@0:     if(A.is_empty())
max@0:       {
max@0:       Z.reset();
max@0:       T.reset();
max@0:       return true;
max@0:       }
max@0:     
max@0:     typedef std::complex<cT> eT;
max@0:     
max@0:     const uword A_n_rows = A.n_rows;
max@0:     
max@0:     char        jobvs = 'V';                // get Schur vectors (Z)
max@0:     char         sort = 'N';                // do not sort eigenvalues/vectors
max@0:     blas_int*  select = 0;                  // pointer to sorting function
max@0:     blas_int        n = blas_int(A_n_rows);
max@0:     blas_int     sdim = 0;                  // output for sorting
max@0:     
max@0:     blas_int lwork = 3 * n;                 // workspace must be at least 3 * n (if set to -1, optimal size is output in work(0) and nothing else is done
max@0:     
max@0:     podarray<eT>       work( static_cast<uword>(lwork) );
max@0:     podarray<blas_int> bwork(A_n_rows);
max@0:     
max@0:     blas_int info = 0;
max@0:     
max@0:     Z.set_size(A_n_rows, A_n_rows);
max@0:     T = A;
max@0:     
max@0:     podarray<eT>     w(A_n_rows);           // output for eigenvalues
max@0:     podarray<cT> rwork(A_n_rows);
max@0:     
max@0:     lapack::cx_gees(&jobvs, &sort, select, &n, T.memptr(), &n, &sdim, w.memptr(), Z.memptr(), &n, work.memptr(), &lwork, rwork.memptr(), bwork.memptr(), &info);
max@0:     
max@0:     return (info == 0);
max@0:     }
max@0:   #else
max@0:     {
max@0:     arma_ignore(Z);
max@0:     arma_ignore(T);
max@0:     arma_stop("schur_dec(): use of LAPACK needs to be enabled");
max@0:     return false;
max@0:     }
max@0:   #endif
max@0:   }
max@0: 
max@0: 
max@0: 
max@0: //
max@0: // syl (solution of the Sylvester equation AX + XB = C)
max@0: 
max@0: template<typename eT>
max@0: inline
max@0: bool
max@0: auxlib::syl(Mat<eT>& X, const Mat<eT>& A, const Mat<eT>& B, const Mat<eT>& C)
max@0:   {
max@0:   arma_extra_debug_sigprint();
max@0:   
max@0:   arma_debug_check
max@0:     (
max@0:     (A.is_square() == false) || (B.is_square() == false),
max@0:     "syl(): given matrix is not square"
max@0:     );
max@0:     
max@0:   arma_debug_check
max@0:     (
max@0:     (C.n_rows != A.n_rows) || (C.n_cols != B.n_cols),
max@0:     "syl(): matrices are not conformant"
max@0:     );
max@0:   
max@0:   if(A.is_empty() || B.is_empty() || C.is_empty())
max@0:     {
max@0:     X.reset();
max@0:     return true;
max@0:     }
max@0:   
max@0:   #if defined(ARMA_USE_LAPACK)
max@0:     {
max@0:     Mat<eT> Z1, Z2, T1, T2;
max@0:     
max@0:     const bool status_sd1 = auxlib::schur_dec(Z1, T1, A);
max@0:     const bool status_sd2 = auxlib::schur_dec(Z2, T2, B);
max@0:     
max@0:     if( (status_sd1 == false) || (status_sd2 == false) )
max@0:       {
max@0:       return false;
max@0:       }
max@0:     
max@0:     char     trana = 'N';
max@0:     char     tranb = 'N';
max@0:     blas_int  isgn = +1;
max@0:     blas_int     m = blas_int(T1.n_rows);
max@0:     blas_int     n = blas_int(T2.n_cols);
max@0:     
max@0:     eT       scale = eT(0);
max@0:     blas_int  info = 0;
max@0:     
max@0:     Mat<eT> Y = trans(Z1) * C * Z2;
max@0:     
max@0:     lapack::trsyl<eT>(&trana, &tranb, &isgn, &m, &n, T1.memptr(), &m, T2.memptr(), &n, Y.memptr(), &m, &scale, &info);
max@0:     
max@0:     //Y /= scale;
max@0:     Y /= (-scale);
max@0:     
max@0:     X = Z1 * Y * trans(Z2);
max@0:     
max@0:     return (info >= 0);
max@0:     }
max@0:   #else
max@0:     {
max@0:     arma_stop("syl(): use of LAPACK needs to be enabled");
max@0:     return false;
max@0:     }
max@0:   #endif
max@0:   }
max@0:   
max@0:   
max@0:   
max@0: //
max@0: // lyap (solution of the continuous Lyapunov equation AX + XA^H + Q = 0)
max@0: 
max@0: template<typename eT>
max@0: inline
max@0: bool
max@0: auxlib::lyap(Mat<eT>& X, const Mat<eT>& A, const Mat<eT>& Q)
max@0:   {
max@0:   arma_extra_debug_sigprint();
max@0:   
max@0:   arma_debug_check( (A.is_square() == false), "lyap(): matrix A is not square");
max@0:   arma_debug_check( (Q.is_square() == false), "lyap(): matrix Q is not square");
max@0:   arma_debug_check( (A.n_rows != Q.n_rows),   "lyap(): matrices A and Q have different dimensions");
max@0:   
max@0:   Mat<eT> htransA;
max@0:   op_htrans::apply_noalias(htransA, A);
max@0:   
max@0:   const Mat<eT> mQ = -Q;
max@0:   
max@0:   return auxlib::syl(X, A, htransA, mQ);
max@0:   }
max@0: 
max@0: 
max@0: 
max@0: //
max@0: // dlyap (solution of the discrete Lyapunov equation AXA^H - X + Q = 0)
max@0: 
max@0: template<typename eT>
max@0: inline
max@0: bool
max@0: auxlib::dlyap(Mat<eT>& X, const Mat<eT>& A, const Mat<eT>& Q)
max@0:   {
max@0:   arma_extra_debug_sigprint();
max@0:   
max@0:   arma_debug_check( (A.is_square() == false), "dlyap(): matrix A is not square");
max@0:   arma_debug_check( (Q.is_square() == false), "dlyap(): matrix Q is not square");
max@0:   arma_debug_check( (A.n_rows != Q.n_rows),   "dlyap(): matrices A and Q have different dimensions");
max@0:   
max@0:   const Col<eT> vecQ = reshape(Q, Q.n_elem, 1);
max@0:   
max@0:   const Mat<eT> M = eye< Mat<eT> >(Q.n_elem, Q.n_elem) - kron(conj(A), A);
max@0:   
max@0:   Col<eT> vecX;
max@0:   
max@0:   const bool status = solve(vecX, M, vecQ);
max@0:   
max@0:   if(status == true)
max@0:     {
max@0:     X = reshape(vecX, Q.n_rows, Q.n_cols);
max@0:     return true;
max@0:     }
max@0:   else
max@0:     {
max@0:     X.reset();
max@0:     return false;
max@0:     }
max@0:   }
max@0: 
max@0: 
max@0: 
max@0: //! @}