annotate armadillo-3.900.4/include/armadillo_bits/auxlib_meat.hpp @ 84:55a047986812 tip

Update library URI so as not to be document-local

// Copyright (C) 2008-2013 NICTA (www.nicta.com.au)
// Copyright (C) 2008-2013 Conrad Sanderson
// Copyright (C) 2009 Edmund Highcock
// Copyright (C) 2011 James Sanders
// Copyright (C) 2011 Stanislav Funiak
// Copyright (C) 2012 Eric Jon Sundstrom
// Copyright (C) 2012 Michael McNeil Forbes
// 
// This Source Code Form is subject to the terms of the Mozilla Public
// License, v. 2.0. If a copy of the MPL was not distributed with this
// file, You can obtain one at http://mozilla.org/MPL/2.0/.


//! \addtogroup auxlib
//! @{



//! immediate matrix inverse
template<typename eT, typename T1>
inline
bool
auxlib::inv(Mat<eT>& out, const Base<eT,T1>& X, const bool slow)
  {
  arma_extra_debug_sigprint();
  
  out = X.get_ref();
  
  arma_debug_check( (out.is_square() == false), "inv(): given matrix is not square" );
  
  bool status = false;
  
  const uword N = out.n_rows;
  
  if( (N <= 4) && (slow == false) )
    {
    status = auxlib::inv_inplace_tinymat(out, N);
    }
    
  if( (N > 4) || (status == false) )
    {
    status = auxlib::inv_inplace_lapack(out);
    }
  
  return status;
  }



template<typename eT>
inline
bool
auxlib::inv(Mat<eT>& out, const Mat<eT>& X, const bool slow)
  {
  arma_extra_debug_sigprint();
  
  arma_debug_check( (X.is_square() == false), "inv(): given matrix is not square" );
  
  bool status = false;
  
  const uword N = X.n_rows;
  
  if( (N <= 4) && (slow == false) )
    {
    status = (&out != &X) ? auxlib::inv_noalias_tinymat(out, X, N) : auxlib::inv_inplace_tinymat(out, N);
    }
  
  if( (N > 4) || (status == false) )
    {
    out = X;
    status = auxlib::inv_inplace_lapack(out);
    }
  
  return status;
  }



template<typename eT>
inline
bool
auxlib::inv_noalias_tinymat(Mat<eT>& out, const Mat<eT>& X, const uword N)
  {
  arma_extra_debug_sigprint();
  
  bool det_ok = true;
  
  out.set_size(N,N);
  
  switch(N)
    {
    case 1:
      {
      out[0] = eT(1) / X[0];
      };
      break;
      
    case 2:
      {
      const eT* Xm = X.memptr();
      
      const eT a = Xm[pos<0,0>::n2];
      const eT b = Xm[pos<0,1>::n2];
      const eT c = Xm[pos<1,0>::n2];
      const eT d = Xm[pos<1,1>::n2];
      
      const eT tmp_det = (a*d - b*c);
      
      if(tmp_det != eT(0))
        {
        eT* outm = out.memptr();
        
        outm[pos<0,0>::n2] =  d / tmp_det;
        outm[pos<0,1>::n2] = -b / tmp_det;
        outm[pos<1,0>::n2] = -c / tmp_det;
        outm[pos<1,1>::n2] =  a / tmp_det;
        }
      else
        {
        det_ok = false;
        }
      };
      break;
    
    case 3:
      {
      const eT* X_col0 = X.colptr(0);
      const eT a11 = X_col0[0];
      const eT a21 = X_col0[1];
      const eT a31 = X_col0[2];
      
      const eT* X_col1 = X.colptr(1);
      const eT a12 = X_col1[0];
      const eT a22 = X_col1[1];
      const eT a32 = X_col1[2];
      
      const eT* X_col2 = X.colptr(2);
      const eT a13 = X_col2[0];
      const eT a23 = X_col2[1];
      const eT a33 = X_col2[2];
      
      const eT tmp_det = a11*(a33*a22 - a32*a23) - a21*(a33*a12-a32*a13) + a31*(a23*a12 - a22*a13);
      
      if(tmp_det != eT(0))
        {
        eT* out_col0 = out.colptr(0);
        out_col0[0] =  (a33*a22 - a32*a23) / tmp_det;
        out_col0[1] = -(a33*a21 - a31*a23) / tmp_det;
        out_col0[2] =  (a32*a21 - a31*a22) / tmp_det;
        
        eT* out_col1 = out.colptr(1);
        out_col1[0] = -(a33*a12 - a32*a13) / tmp_det;
        out_col1[1] =  (a33*a11 - a31*a13) / tmp_det;
        out_col1[2] = -(a32*a11 - a31*a12) / tmp_det;
        
        eT* out_col2 = out.colptr(2);
        out_col2[0] =  (a23*a12 - a22*a13) / tmp_det;
        out_col2[1] = -(a23*a11 - a21*a13) / tmp_det;
        out_col2[2] =  (a22*a11 - a21*a12) / tmp_det;
        }
      else
        {
        det_ok = false;
        }
      };
      break;
      
    case 4:
      {
      const eT tmp_det = det(X);
      
      if(tmp_det != eT(0))
        {
        const eT* Xm   = X.memptr();
              eT* outm = out.memptr();
        
        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;
        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;
        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;
        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;
        
        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;
        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;
        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;
        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;
        
        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;
        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;
        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;
        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;
        
        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;
        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;
        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;
        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;
        }
      else
        {
        det_ok = false;
        }
      };
      break;
    
    default:
      ;
    }
  
  return det_ok;
  }



template<typename eT>
inline
bool
auxlib::inv_inplace_tinymat(Mat<eT>& X, const uword N)
  {
  arma_extra_debug_sigprint();
  
  bool det_ok = true;
  
  // for more info, see:
  // http://www.dr-lex.34sp.com/random/matrix_inv.html
  // http://www.cvl.iis.u-tokyo.ac.jp/~miyazaki/tech/teche23.html
  // http://www.euclideanspace.com/maths/algebra/matrix/functions/inverse/fourD/index.htm
  // http://www.geometrictools.com//LibFoundation/Mathematics/Wm4Matrix4.inl
  
  switch(N)
    {
    case 1:
      {
      X[0] = eT(1) / X[0];
      };
      break;
      
    case 2:
      {
      const eT a = X[pos<0,0>::n2];
      const eT b = X[pos<0,1>::n2];
      const eT c = X[pos<1,0>::n2];
      const eT d = X[pos<1,1>::n2];
      
      const eT tmp_det = (a*d - b*c);
      
      if(tmp_det != eT(0))
        {
        X[pos<0,0>::n2] =  d / tmp_det;
        X[pos<0,1>::n2] = -b / tmp_det;
        X[pos<1,0>::n2] = -c / tmp_det;
        X[pos<1,1>::n2] =  a / tmp_det;
        }
      else
        {
        det_ok = false;
        }
      };
      break;
    
    case 3:
      {
      eT* X_col0 = X.colptr(0);
      eT* X_col1 = X.colptr(1);
      eT* X_col2 = X.colptr(2);
      
      const eT a11 = X_col0[0];
      const eT a21 = X_col0[1];
      const eT a31 = X_col0[2];
      
      const eT a12 = X_col1[0];
      const eT a22 = X_col1[1];
      const eT a32 = X_col1[2];
      
      const eT a13 = X_col2[0];
      const eT a23 = X_col2[1];
      const eT a33 = X_col2[2];
      
      const eT tmp_det = a11*(a33*a22 - a32*a23) - a21*(a33*a12-a32*a13) + a31*(a23*a12 - a22*a13);
      
      if(tmp_det != eT(0))
        {
        X_col0[0] =  (a33*a22 - a32*a23) / tmp_det;
        X_col0[1] = -(a33*a21 - a31*a23) / tmp_det;
        X_col0[2] =  (a32*a21 - a31*a22) / tmp_det;
        
        X_col1[0] = -(a33*a12 - a32*a13) / tmp_det;
        X_col1[1] =  (a33*a11 - a31*a13) / tmp_det;
        X_col1[2] = -(a32*a11 - a31*a12) / tmp_det;
        
        X_col2[0] =  (a23*a12 - a22*a13) / tmp_det;
        X_col2[1] = -(a23*a11 - a21*a13) / tmp_det;
        X_col2[2] =  (a22*a11 - a21*a12) / tmp_det;
        }
      else
        {
        det_ok = false;
        }
      };
      break;
      
    case 4:
      {
      const eT tmp_det = det(X);
      
      if(tmp_det != eT(0))
        {
        const Mat<eT> A(X);
        
        const eT* Am = A.memptr();
              eT* Xm = X.memptr();
        
        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;
        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;
        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;
        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;
        
        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;
        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;
        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;
        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;
        
        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;
        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;
        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;
        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;
        
        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;
        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;
        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;
        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;
        }
      else
        {
        det_ok = false;
        }
      };
      break;
      
    default:
      ;
    }
  
  return det_ok;
  }



template<typename eT>
inline
bool
auxlib::inv_inplace_lapack(Mat<eT>& out)
  {
  arma_extra_debug_sigprint();

  if(out.is_empty())
    {
    return true;
    }
  
  #if defined(ARMA_USE_ATLAS)
    {
    podarray<int> ipiv(out.n_rows);
    
    int info = atlas::clapack_getrf(atlas::CblasColMajor, out.n_rows, out.n_cols, out.memptr(), out.n_rows, ipiv.memptr());
    
    if(info == 0)
      {
      info = atlas::clapack_getri(atlas::CblasColMajor, out.n_rows, out.memptr(), out.n_rows, ipiv.memptr());
      }
    
    return (info == 0);
    }
  #elif defined(ARMA_USE_LAPACK)
    {
    blas_int n_rows    = out.n_rows;
    blas_int n_cols    = out.n_cols;
    blas_int lwork     = 0;
    blas_int lwork_min = (std::max)(blas_int(1), n_rows);
    blas_int info      = 0;
    
    podarray<blas_int> ipiv(out.n_rows);
    
    eT        work_query[2];
    blas_int lwork_query = -1;
    
    lapack::getri(&n_rows, out.memptr(), &n_rows, ipiv.memptr(), &work_query[0], &lwork_query, &info);
    
    if(info == 0)
      {
      const blas_int lwork_proposed = static_cast<blas_int>( access::tmp_real(work_query[0]) );
      
      lwork = (lwork_proposed > lwork_min) ? lwork_proposed : lwork_min;
      }
    else
      {
      return false;
      }
    
    podarray<eT> work( static_cast<uword>(lwork) );
    
    lapack::getrf(&n_rows, &n_cols, out.memptr(), &n_rows, ipiv.memptr(), &info);
    
    if(info == 0)
      {
      lapack::getri(&n_rows, out.memptr(), &n_rows, ipiv.memptr(), work.memptr(), &lwork, &info);
      }
    
    return (info == 0);
    }
  #else
    {
    arma_stop("inv(): use of ATLAS or LAPACK needs to be enabled");
    return false;
    }
  #endif
  }



template<typename eT, typename T1>
inline
bool
auxlib::inv_tr(Mat<eT>& out, const Base<eT,T1>& X, const uword layout)
  {
  arma_extra_debug_sigprint();
  
  out = X.get_ref();
  
  arma_debug_check( (out.is_square() == false), "inv(): given matrix is not square" );
  
  if(out.is_empty())
    {
    return true;
    }
  
  bool status;
  
  #if defined(ARMA_USE_LAPACK)
    {
    char     uplo = (layout == 0) ? 'U' : 'L';
    char     diag = 'N';
    blas_int n    = blas_int(out.n_rows);
    blas_int info = 0;
    
    lapack::trtri(&uplo, &diag, &n, out.memptr(), &n, &info);
    
    status = (info == 0);
    }
  #else
    {
    arma_ignore(layout);
    arma_stop("inv(): use of LAPACK needs to be enabled");
    status = false;
    }
  #endif
  
  
  if(status == true)
    {
    if(layout == 0)
      {
      // upper triangular
      out = trimatu(out);
      }
    else
      {
      // lower triangular      
      out = trimatl(out);
      }
    }
  
  return status;
  }



template<typename eT, typename T1>
inline
bool
auxlib::inv_sym(Mat<eT>& out, const Base<eT,T1>& X, const uword layout)
  {
  arma_extra_debug_sigprint();
  
  out = X.get_ref();
  
  arma_debug_check( (out.is_square() == false), "inv(): given matrix is not square" );
  
  if(out.is_empty())
    {
    return true;
    }
  
  bool status;
  
  #if defined(ARMA_USE_LAPACK)
    {
    char     uplo  = (layout == 0) ? 'U' : 'L';
    blas_int n     = blas_int(out.n_rows);
    blas_int lwork = 3 * (n*n); // TODO: use lwork = -1 to determine optimal size
    blas_int info  = 0;
    
    podarray<blas_int> ipiv;
    ipiv.set_size(out.n_rows);
    
    podarray<eT> work;
    work.set_size( uword(lwork) );
    
    lapack::sytrf(&uplo, &n, out.memptr(), &n, ipiv.memptr(), work.memptr(), &lwork, &info);
    
    status = (info == 0);
    
    if(status == true)
      {
      lapack::sytri(&uplo, &n, out.memptr(), &n, ipiv.memptr(), work.memptr(), &info);
      
      out = (layout == 0) ? symmatu(out) : symmatl(out);
      
      status = (info == 0);
      }
    }
  #else
    {
    arma_ignore(layout);
    arma_stop("inv(): use of LAPACK needs to be enabled");
    status = false;
    }
  #endif
  
  return status;
  }



template<typename eT, typename T1>
inline
bool
auxlib::inv_sympd(Mat<eT>& out, const Base<eT,T1>& X, const uword layout)
  {
  arma_extra_debug_sigprint();
  
  out = X.get_ref();
  
  arma_debug_check( (out.is_square() == false), "inv(): given matrix is not square" );
  
  if(out.is_empty())
    {
    return true;
    }
  
  bool status;
  
  #if defined(ARMA_USE_LAPACK)
    {
    char     uplo = (layout == 0) ? 'U' : 'L';
    blas_int n    = blas_int(out.n_rows);
    blas_int info = 0;
    
    lapack::potrf(&uplo, &n, out.memptr(), &n, &info);
    
    status = (info == 0);
    
    if(status == true)
      {
      lapack::potri(&uplo, &n, out.memptr(), &n, &info);
    
      out = (layout == 0) ? symmatu(out) : symmatl(out);
    
      status = (info == 0);
      }
    }
  #else
    {
    arma_ignore(layout);
    arma_stop("inv(): use of LAPACK needs to be enabled");
    status = false;
    }
  #endif
  
  return status;
  }



template<typename eT, typename T1>
inline
eT
auxlib::det(const Base<eT,T1>& X, const bool slow)
  {
  const unwrap<T1>   tmp(X.get_ref());
  const Mat<eT>& A = tmp.M;
  
  arma_debug_check( (A.is_square() == false), "det(): matrix is not square" );
  
  const bool make_copy = (is_Mat<T1>::value == true) ? true : false;
  
  if(slow == false)
    {
    const uword N = A.n_rows;
    
    switch(N)
      {
      case 0:
      case 1:
      case 2:
        return auxlib::det_tinymat(A, N);
        break;
      
      case 3:
      case 4:
        {
        const eT tmp_det = auxlib::det_tinymat(A, N);
        return (tmp_det != eT(0)) ? tmp_det : auxlib::det_lapack(A, make_copy);
        }
        break;
      
      default:
        return auxlib::det_lapack(A, make_copy);
      }
    }
  
  return auxlib::det_lapack(A, make_copy);
  }



template<typename eT>
inline
eT
auxlib::det_tinymat(const Mat<eT>& X, const uword N)
  {
  arma_extra_debug_sigprint();
  
  switch(N)
    {
    case 0:
      return eT(1);
      break;
    
    case 1:
      return X[0];
      break;
    
    case 2:
      {
      const eT* Xm = X.memptr();
      
      return ( Xm[pos<0,0>::n2]*Xm[pos<1,1>::n2] - Xm[pos<0,1>::n2]*Xm[pos<1,0>::n2] );
      }
      break;
    
    case 3:
      {
      // const double tmp1 = X.at(0,0) * X.at(1,1) * X.at(2,2);
      // const double tmp2 = X.at(0,1) * X.at(1,2) * X.at(2,0);
      // const double tmp3 = X.at(0,2) * X.at(1,0) * X.at(2,1);
      // const double tmp4 = X.at(2,0) * X.at(1,1) * X.at(0,2);
      // const double tmp5 = X.at(2,1) * X.at(1,2) * X.at(0,0);
      // const double tmp6 = X.at(2,2) * X.at(1,0) * X.at(0,1);
      // return (tmp1+tmp2+tmp3) - (tmp4+tmp5+tmp6);
      
      const eT* a_col0 = X.colptr(0);
      const eT  a11    = a_col0[0];
      const eT  a21    = a_col0[1];
      const eT  a31    = a_col0[2];
      
      const eT* a_col1 = X.colptr(1);
      const eT  a12    = a_col1[0];
      const eT  a22    = a_col1[1];
      const eT  a32    = a_col1[2];
      
      const eT* a_col2 = X.colptr(2);
      const eT  a13    = a_col2[0];
      const eT  a23    = a_col2[1];
      const eT  a33    = a_col2[2];
      
      return ( a11*(a33*a22 - a32*a23) - a21*(a33*a12-a32*a13) + a31*(a23*a12 - a22*a13) );
      }
      break;
      
    case 4:
      {
      const eT* Xm = X.memptr();
      
      const eT val = \
          Xm[pos<0,3>::n4] * Xm[pos<1,2>::n4] * Xm[pos<2,1>::n4] * Xm[pos<3,0>::n4] \
        - Xm[pos<0,2>::n4] * Xm[pos<1,3>::n4] * Xm[pos<2,1>::n4] * Xm[pos<3,0>::n4] \
        - Xm[pos<0,3>::n4] * Xm[pos<1,1>::n4] * Xm[pos<2,2>::n4] * Xm[pos<3,0>::n4] \
        + Xm[pos<0,1>::n4] * Xm[pos<1,3>::n4] * Xm[pos<2,2>::n4] * Xm[pos<3,0>::n4] \
        + Xm[pos<0,2>::n4] * Xm[pos<1,1>::n4] * Xm[pos<2,3>::n4] * Xm[pos<3,0>::n4] \
        - Xm[pos<0,1>::n4] * Xm[pos<1,2>::n4] * Xm[pos<2,3>::n4] * Xm[pos<3,0>::n4] \
        - Xm[pos<0,3>::n4] * Xm[pos<1,2>::n4] * Xm[pos<2,0>::n4] * Xm[pos<3,1>::n4] \
        + Xm[pos<0,2>::n4] * Xm[pos<1,3>::n4] * Xm[pos<2,0>::n4] * Xm[pos<3,1>::n4] \
        + Xm[pos<0,3>::n4] * Xm[pos<1,0>::n4] * Xm[pos<2,2>::n4] * Xm[pos<3,1>::n4] \
        - Xm[pos<0,0>::n4] * Xm[pos<1,3>::n4] * Xm[pos<2,2>::n4] * Xm[pos<3,1>::n4] \
        - Xm[pos<0,2>::n4] * Xm[pos<1,0>::n4] * Xm[pos<2,3>::n4] * Xm[pos<3,1>::n4] \
        + Xm[pos<0,0>::n4] * Xm[pos<1,2>::n4] * Xm[pos<2,3>::n4] * Xm[pos<3,1>::n4] \
        + Xm[pos<0,3>::n4] * Xm[pos<1,1>::n4] * Xm[pos<2,0>::n4] * Xm[pos<3,2>::n4] \
        - Xm[pos<0,1>::n4] * Xm[pos<1,3>::n4] * Xm[pos<2,0>::n4] * Xm[pos<3,2>::n4] \
        - Xm[pos<0,3>::n4] * Xm[pos<1,0>::n4] * Xm[pos<2,1>::n4] * Xm[pos<3,2>::n4] \
        + Xm[pos<0,0>::n4] * Xm[pos<1,3>::n4] * Xm[pos<2,1>::n4] * Xm[pos<3,2>::n4] \
        + Xm[pos<0,1>::n4] * Xm[pos<1,0>::n4] * Xm[pos<2,3>::n4] * Xm[pos<3,2>::n4] \
        - Xm[pos<0,0>::n4] * Xm[pos<1,1>::n4] * Xm[pos<2,3>::n4] * Xm[pos<3,2>::n4] \
        - Xm[pos<0,2>::n4] * Xm[pos<1,1>::n4] * Xm[pos<2,0>::n4] * Xm[pos<3,3>::n4] \
        + Xm[pos<0,1>::n4] * Xm[pos<1,2>::n4] * Xm[pos<2,0>::n4] * Xm[pos<3,3>::n4] \
        + Xm[pos<0,2>::n4] * Xm[pos<1,0>::n4] * Xm[pos<2,1>::n4] * Xm[pos<3,3>::n4] \
        - Xm[pos<0,0>::n4] * Xm[pos<1,2>::n4] * Xm[pos<2,1>::n4] * Xm[pos<3,3>::n4] \
        - Xm[pos<0,1>::n4] * Xm[pos<1,0>::n4] * Xm[pos<2,2>::n4] * Xm[pos<3,3>::n4] \
        + Xm[pos<0,0>::n4] * Xm[pos<1,1>::n4] * Xm[pos<2,2>::n4] * Xm[pos<3,3>::n4] \
        ;
      
      return val;
      }
      break;
    
    default:
      return eT(0);
      ;
    }
  }



//! immediate determinant of a matrix using ATLAS or LAPACK
template<typename eT>
inline
eT
auxlib::det_lapack(const Mat<eT>& X, const bool make_copy)
  {
  arma_extra_debug_sigprint();
  
  Mat<eT> X_copy;
  
  if(make_copy == true)
    {
    X_copy = X;
    }
  
  Mat<eT>& tmp = (make_copy == true) ? X_copy : const_cast< Mat<eT>& >(X);
  
  if(tmp.is_empty())
    {
    return eT(1);
    }
  
  
  #if defined(ARMA_USE_ATLAS)
    {
    podarray<int> ipiv(tmp.n_rows);
    
    //const int info =
    atlas::clapack_getrf(atlas::CblasColMajor, tmp.n_rows, tmp.n_cols, tmp.memptr(), tmp.n_rows, ipiv.memptr());
    
    // on output tmp appears to be L+U_alt, where U_alt is U with the main diagonal set to zero
    eT val = tmp.at(0,0);
    for(uword i=1; i < tmp.n_rows; ++i)
      {
      val *= tmp.at(i,i);
      }
    
    int sign = +1;
    for(uword i=0; i < tmp.n_rows; ++i)
      {
      if( int(i) != ipiv.mem[i] )  // NOTE: no adjustment required, as the clapack version of getrf() assumes counting from 0
        {
        sign *= -1;
        }
      }
    
    return ( (sign < 0) ? -val : val );
    }
  #elif defined(ARMA_USE_LAPACK)
    {
    podarray<blas_int> ipiv(tmp.n_rows);
    
    blas_int info   = 0;
    blas_int n_rows = blas_int(tmp.n_rows);
    blas_int n_cols = blas_int(tmp.n_cols);
    
    lapack::getrf(&n_rows, &n_cols, tmp.memptr(), &n_rows, ipiv.memptr(), &info);
    
    // on output tmp appears to be L+U_alt, where U_alt is U with the main diagonal set to zero
    eT val = tmp.at(0,0);
    for(uword i=1; i < tmp.n_rows; ++i)
      {
      val *= tmp.at(i,i);
      }
    
    blas_int sign = +1;
    for(uword i=0; i < tmp.n_rows; ++i)
      {
      if( blas_int(i) != (ipiv.mem[i] - 1) )  // NOTE: adjustment of -1 is required as Fortran counts from 1
        {
        sign *= -1;
        }
      }
    
    return ( (sign < 0) ? -val : val );
    }
  #else
    {
    arma_ignore(X);
    arma_ignore(make_copy);
    arma_ignore(tmp);
    arma_stop("det(): use of ATLAS or LAPACK needs to be enabled");
    return eT(0);
    }
  #endif
  }



//! immediate log determinant of a matrix using ATLAS or LAPACK
template<typename eT, typename T1>
inline
bool
auxlib::log_det(eT& out_val, typename get_pod_type<eT>::result& out_sign, const Base<eT,T1>& X)
  {
  arma_extra_debug_sigprint();
  
  typedef typename get_pod_type<eT>::result T;
  
  #if defined(ARMA_USE_ATLAS)
    {
    Mat<eT> tmp(X.get_ref());
    arma_debug_check( (tmp.is_square() == false), "log_det(): given matrix is not square" );
    
    if(tmp.is_empty())
      {
      out_val  = eT(0);
      out_sign =  T(1);
      return true;
      }
    
    podarray<int> ipiv(tmp.n_rows);
    
    const int info = atlas::clapack_getrf(atlas::CblasColMajor, tmp.n_rows, tmp.n_cols, tmp.memptr(), tmp.n_rows, ipiv.memptr());
    
    // on output tmp appears to be L+U_alt, where U_alt is U with the main diagonal set to zero
    
    sword sign = (is_complex<eT>::value == false) ? ( (access::tmp_real( tmp.at(0,0) ) < T(0)) ? -1 : +1 ) : +1;
    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) );
    
    for(uword i=1; i < tmp.n_rows; ++i)
      {
      const eT x = tmp.at(i,i);
      
      sign *= (is_complex<eT>::value == false) ? ( (access::tmp_real(x) < T(0)) ? -1 : +1 ) : +1;
      val  += (is_complex<eT>::value == false) ? std::log( (access::tmp_real(x) < T(0)) ? x*T(-1) : x ) : std::log(x);
      }
    
    for(uword i=0; i < tmp.n_rows; ++i)
      {
      if( int(i) != ipiv.mem[i] )  // NOTE: no adjustment required, as the clapack version of getrf() assumes counting from 0
        {
        sign *= -1;
        }
      }
    
    out_val  = val;
    out_sign = T(sign);
    
    return (info == 0);
    }
  #elif defined(ARMA_USE_LAPACK)
    {
    Mat<eT> tmp(X.get_ref());
    arma_debug_check( (tmp.is_square() == false), "log_det(): given matrix is not square" );
    
    if(tmp.is_empty())
      {
      out_val  = eT(0);
      out_sign =  T(1);
      return true;
      }
    
    podarray<blas_int> ipiv(tmp.n_rows);
    
    blas_int info   = 0;
    blas_int n_rows = blas_int(tmp.n_rows);
    blas_int n_cols = blas_int(tmp.n_cols);
    
    lapack::getrf(&n_rows, &n_cols, tmp.memptr(), &n_rows, ipiv.memptr(), &info);
    
    // on output tmp appears to be L+U_alt, where U_alt is U with the main diagonal set to zero
    
    sword sign = (is_complex<eT>::value == false) ? ( (access::tmp_real( tmp.at(0,0) ) < T(0)) ? -1 : +1 ) : +1;
    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) );
    
    for(uword i=1; i < tmp.n_rows; ++i)
      {
      const eT x = tmp.at(i,i);
      
      sign *= (is_complex<eT>::value == false) ? ( (access::tmp_real(x) < T(0)) ? -1 : +1 ) : +1;
      val  += (is_complex<eT>::value == false) ? std::log( (access::tmp_real(x) < T(0)) ? x*T(-1) : x ) : std::log(x);
      }
    
    for(uword i=0; i < tmp.n_rows; ++i)
      {
      if( blas_int(i) != (ipiv.mem[i] - 1) )  // NOTE: adjustment of -1 is required as Fortran counts from 1
        {
        sign *= -1;
        }
      }
    
    out_val  = val;
    out_sign = T(sign);
    
    return (info == 0);
    }
  #else
    {
    arma_ignore(X);
    
    out_val  = eT(0);
    out_sign =  T(0);
    
    arma_stop("log_det(): use of ATLAS or LAPACK needs to be enabled");
    
    return false;
    }
  #endif
  }



//! immediate LU decomposition of a matrix using ATLAS or LAPACK
template<typename eT, typename T1>
inline
bool
auxlib::lu(Mat<eT>& L, Mat<eT>& U, podarray<blas_int>& ipiv, const Base<eT,T1>& X)
  {
  arma_extra_debug_sigprint();
  
  U = X.get_ref();
  
  const uword U_n_rows = U.n_rows;
  const uword U_n_cols = U.n_cols;
  
  if(U.is_empty())
    {
    L.set_size(U_n_rows, 0);
    U.set_size(0, U_n_cols);
    ipiv.reset();
    return true;
    }
  
  #if defined(ARMA_USE_ATLAS) || defined(ARMA_USE_LAPACK)
    {
    bool status;
    
    #if defined(ARMA_USE_ATLAS)
      {
      ipiv.set_size( (std::min)(U_n_rows, U_n_cols) );
      
      int info = atlas::clapack_getrf(atlas::CblasColMajor, U_n_rows, U_n_cols, U.memptr(), U_n_rows, ipiv.memptr());
      
      status = (info == 0);
      }
    #elif defined(ARMA_USE_LAPACK)
      {
      ipiv.set_size( (std::min)(U_n_rows, U_n_cols) );
      
      blas_int info = 0;
      
      blas_int n_rows = U_n_rows;
      blas_int n_cols = U_n_cols;
      
      
      lapack::getrf(&n_rows, &n_cols, U.memptr(), &n_rows, ipiv.memptr(), &info);
      
      // take into account that Fortran counts from 1
      arrayops::inplace_minus(ipiv.memptr(), blas_int(1), ipiv.n_elem);
      
      status = (info == 0);
      }
    #endif
    
    L.copy_size(U);
    
    for(uword col=0; col < U_n_cols; ++col)
      {
      for(uword row=0; (row < col) && (row < U_n_rows); ++row)
        {
        L.at(row,col) = eT(0);
        }
      
      if( L.in_range(col,col) == true )
        {
        L.at(col,col) = eT(1);
        }
      
      for(uword row = (col+1); row < U_n_rows; ++row)
        {
        L.at(row,col) = U.at(row,col);
        U.at(row,col) = eT(0);
        }
      }
    
    return status;
    }
  #else
    {
    arma_stop("lu(): use of ATLAS or LAPACK needs to be enabled");
    
    return false;
    }
  #endif
  }



template<typename eT, typename T1>
inline
bool
auxlib::lu(Mat<eT>& L, Mat<eT>& U, Mat<eT>& P, const Base<eT,T1>& X)
  {
  arma_extra_debug_sigprint();
  
  podarray<blas_int> ipiv1;
  const bool status = auxlib::lu(L, U, ipiv1, X);
  
  if(status == true)
    {
    if(U.is_empty())
      {
      // L and U have been already set to the correct empty matrices
      P.eye(L.n_rows, L.n_rows);
      return true;
      }
    
    const uword n      = ipiv1.n_elem;
    const uword P_rows = U.n_rows;
    
    podarray<blas_int> ipiv2(P_rows);
    
    const blas_int* ipiv1_mem = ipiv1.memptr();
          blas_int* ipiv2_mem = ipiv2.memptr();
    
    for(uword i=0; i<P_rows; ++i)
      {
      ipiv2_mem[i] = blas_int(i);
      }
    
    for(uword i=0; i<n; ++i)
      {
      const uword k = static_cast<uword>(ipiv1_mem[i]);
      
      if( ipiv2_mem[i] != ipiv2_mem[k] )
        {
        std::swap( ipiv2_mem[i], ipiv2_mem[k] );
        }
      }
    
    P.zeros(P_rows, P_rows);
    
    for(uword row=0; row<P_rows; ++row)
      {
      P.at(row, static_cast<uword>(ipiv2_mem[row])) = eT(1);
      }
    
    if(L.n_cols > U.n_rows)
      {
      L.shed_cols(U.n_rows, L.n_cols-1);
      }
      
    if(U.n_rows > L.n_cols)
      {
      U.shed_rows(L.n_cols, U.n_rows-1);
      }
    }
  
  return status;
  }



template<typename eT, typename T1>
inline
bool
auxlib::lu(Mat<eT>& L, Mat<eT>& U, const Base<eT,T1>& X)
  {
  arma_extra_debug_sigprint();
  
  podarray<blas_int> ipiv1;
  const bool status = auxlib::lu(L, U, ipiv1, X);
  
  if(status == true)
    {
    if(U.is_empty())
      {
      // L and U have been already set to the correct empty matrices
      return true;
      }
    
    const uword n      = ipiv1.n_elem;
    const uword P_rows = U.n_rows;
    
    podarray<blas_int> ipiv2(P_rows);
    
    const blas_int* ipiv1_mem = ipiv1.memptr();
          blas_int* ipiv2_mem = ipiv2.memptr();
    
    for(uword i=0; i<P_rows; ++i)
      {
      ipiv2_mem[i] = blas_int(i);
      }
    
    for(uword i=0; i<n; ++i)
      {
      const uword k = static_cast<uword>(ipiv1_mem[i]);
      
      if( ipiv2_mem[i] != ipiv2_mem[k] )
        {
        std::swap( ipiv2_mem[i], ipiv2_mem[k] );
        L.swap_rows( static_cast<uword>(ipiv2_mem[i]), static_cast<uword>(ipiv2_mem[k]) );
        }
      }
    
    if(L.n_cols > U.n_rows)
      {
      L.shed_cols(U.n_rows, L.n_cols-1);
      }
      
    if(U.n_rows > L.n_cols)
      {
      U.shed_rows(L.n_cols, U.n_rows-1);
      }
    }
  
  return status;
  }



//! immediate eigenvalues of a symmetric real matrix using LAPACK
template<typename eT, typename T1>
inline
bool
auxlib::eig_sym(Col<eT>& eigval, const Base<eT,T1>& X)
  {
  arma_extra_debug_sigprint();
  
  #if defined(ARMA_USE_LAPACK)
    {
    Mat<eT> A(X.get_ref());
    
    arma_debug_check( (A.is_square() == false), "eig_sym(): given matrix is not square");
    
    if(A.is_empty())
      {
      eigval.reset();
      return true;
      }
    
    eigval.set_size(A.n_rows);
    
    char jobz  = 'N';
    char uplo  = 'U';
    
    blas_int N     = blas_int(A.n_rows);
    blas_int lwork = 3 * ( (std::max)(blas_int(1), 3*N-1) );
    blas_int info  = 0;
    
    podarray<eT> work( static_cast<uword>(lwork) );
    
    lapack::syev(&jobz, &uplo, &N, A.memptr(), &N, eigval.memptr(), work.memptr(), &lwork, &info);
    
    return (info == 0);
    }
  #else
    {
    arma_ignore(eigval);
    arma_ignore(X);
    arma_stop("eig_sym(): use of LAPACK needs to be enabled");
    return false;
    }
  #endif
  }



//! immediate eigenvalues of a hermitian complex matrix using LAPACK
template<typename T, typename T1>
inline
bool
auxlib::eig_sym(Col<T>& eigval, const Base<std::complex<T>,T1>& X)
  {
  arma_extra_debug_sigprint();
  
  #if defined(ARMA_USE_LAPACK)
    {
    typedef typename std::complex<T> eT;
    
    Mat<eT> A(X.get_ref());
    
    arma_debug_check( (A.is_square() == false), "eig_sym(): given matrix is not square");
    
    if(A.is_empty())
      {
      eigval.reset();
      return true;
      }
    
    eigval.set_size(A.n_rows);
    
    char jobz  = 'N'; 
    char uplo  = 'U';
    
    blas_int N     = blas_int(A.n_rows);
    blas_int lwork = 3 * ( (std::max)(blas_int(1), 2*N-1) );
    blas_int info  = 0;
    
    podarray<eT>  work( static_cast<uword>(lwork) );
    podarray<T>  rwork( static_cast<uword>( (std::max)(blas_int(1), 3*N-2) ) );
    
    arma_extra_debug_print("lapack::heev()");
    lapack::heev(&jobz, &uplo, &N, A.memptr(), &N, eigval.memptr(), work.memptr(), &lwork, rwork.memptr(), &info);
    
    return (info == 0);
    }
  #else
    {
    arma_ignore(eigval);
    arma_ignore(X);
    arma_stop("eig_sym(): use of LAPACK needs to be enabled");
    return false;
    }
  #endif
  }



//! immediate eigenvalues and eigenvectors of a symmetric real matrix using LAPACK
template<typename eT, typename T1>
inline
bool
auxlib::eig_sym(Col<eT>& eigval, Mat<eT>& eigvec, const Base<eT,T1>& X)
  {
  arma_extra_debug_sigprint();
  
  #if defined(ARMA_USE_LAPACK)
    {
    eigvec = X.get_ref();
    
    arma_debug_check( (eigvec.is_square() == false), "eig_sym(): given matrix is not square" );
    
    if(eigvec.is_empty())
      {
      eigval.reset();
      eigvec.reset();
      return true;
      }
    
    eigval.set_size(eigvec.n_rows);
    
    char jobz  = 'V';
    char uplo  = 'U';
    
    blas_int N     = blas_int(eigvec.n_rows);
    blas_int lwork = 3 * ( (std::max)(blas_int(1), 3*N-1) );
    blas_int info  = 0;
    
    podarray<eT> work( static_cast<uword>(lwork) );
    
    lapack::syev(&jobz, &uplo, &N, eigvec.memptr(), &N, eigval.memptr(), work.memptr(), &lwork, &info);
    
    return (info == 0);
    }
  #else
    {
    arma_ignore(eigval);
    arma_ignore(eigvec);
    arma_ignore(X);
    arma_stop("eig_sym(): use of LAPACK needs to be enabled");
    return false;
    }
  #endif
  }



//! immediate eigenvalues and eigenvectors of a hermitian complex matrix using LAPACK
template<typename T, typename T1>
inline
bool
auxlib::eig_sym(Col<T>& eigval, Mat< std::complex<T> >& eigvec, const Base<std::complex<T>,T1>& X)
  {
  arma_extra_debug_sigprint();
  
  #if defined(ARMA_USE_LAPACK)
    {
    typedef typename std::complex<T> eT;
    
    eigvec = X.get_ref();
    
    arma_debug_check( (eigvec.is_square() == false), "eig_sym(): given matrix is not square" );
    
    if(eigvec.is_empty())
      {
      eigval.reset();
      eigvec.reset();
      return true;
      }
    
    eigval.set_size(eigvec.n_rows);
    
    char jobz  = 'V';
    char uplo  = 'U';
    
    blas_int N     = blas_int(eigvec.n_rows);
    blas_int lwork = 3 * ( (std::max)(blas_int(1), 2*N-1) );
    blas_int info  = 0;
    
    podarray<eT>  work( static_cast<uword>(lwork) );
    podarray<T>  rwork( static_cast<uword>((std::max)(blas_int(1), 3*N-2)) );
    
    arma_extra_debug_print("lapack::heev()");
    lapack::heev(&jobz, &uplo, &N, eigvec.memptr(), &N, eigval.memptr(), work.memptr(), &lwork, rwork.memptr(), &info);
    
    return (info == 0);
    }
  #else
    {
    arma_ignore(eigval);
    arma_ignore(eigvec);
    arma_ignore(X);
    arma_stop("eig_sym(): use of LAPACK needs to be enabled");
    return false;
    }
  #endif
  }



//! immediate eigenvalues and eigenvectors of a symmetric real matrix using LAPACK (divide and conquer algorithm)
template<typename eT, typename T1>
inline
bool
auxlib::eig_sym_dc(Col<eT>& eigval, Mat<eT>& eigvec, const Base<eT,T1>& X)
  {
  arma_extra_debug_sigprint();
  
  #if defined(ARMA_USE_LAPACK)
    {
    eigvec = X.get_ref();
    
    arma_debug_check( (eigvec.is_square() == false), "eig_sym(): given matrix is not square" );
    
    if(eigvec.is_empty())
      {
      eigval.reset();
      eigvec.reset();
      return true;
      }
    
    eigval.set_size(eigvec.n_rows);
    
    char jobz = 'V';
    char uplo = 'U';
    
    blas_int N      = blas_int(eigvec.n_rows);
    blas_int lwork  = 3 * (1 + 6*N + 2*(N*N));
    blas_int liwork = 3 * (3 + 5*N + 2);
    blas_int info   = 0;
    
    podarray<eT>        work( static_cast<uword>( lwork) );
    podarray<blas_int> iwork( static_cast<uword>(liwork) ); 
    
    arma_extra_debug_print("lapack::syevd()");
    lapack::syevd(&jobz, &uplo, &N, eigvec.memptr(), &N, eigval.memptr(), work.memptr(), &lwork, iwork.memptr(), &liwork, &info);
    
    return (info == 0);
    }
  #else
    {
    arma_ignore(eigval);
    arma_ignore(eigvec);
    arma_ignore(X);
    arma_stop("eig_sym(): use of LAPACK needs to be enabled");
    return false;
    }
  #endif
  }



//! immediate eigenvalues and eigenvectors of a hermitian complex matrix using LAPACK (divide and conquer algorithm)
template<typename T, typename T1>
inline
bool
auxlib::eig_sym_dc(Col<T>& eigval, Mat< std::complex<T> >& eigvec, const Base<std::complex<T>,T1>& X)
  {
  arma_extra_debug_sigprint();
  
  #if defined(ARMA_USE_LAPACK)
    {
    typedef typename std::complex<T> eT;
    
    eigvec = X.get_ref();
    
    arma_debug_check( (eigvec.is_square() == false), "eig_sym(): given matrix is not square" );
    
    if(eigvec.is_empty())
      {
      eigval.reset();
      eigvec.reset();
      return true;
      }
    
    eigval.set_size(eigvec.n_rows);
    
    char jobz  = 'V';
    char uplo  = 'U';
    
    blas_int N      = blas_int(eigvec.n_rows);
    blas_int lwork  = 3 * (2*N + N*N);
    blas_int lrwork = 3 * (1 + 5*N + 2*(N*N));
    blas_int liwork = 3 * (3 + 5*N);
    blas_int info   = 0;
    
    podarray<eT>        work( static_cast<uword>(lwork)  );
    podarray<T>        rwork( static_cast<uword>(lrwork) );
    podarray<blas_int> iwork( static_cast<uword>(liwork) ); 
    
    arma_extra_debug_print("lapack::heevd()");
    lapack::heevd(&jobz, &uplo, &N, eigvec.memptr(), &N, eigval.memptr(), work.memptr(), &lwork, rwork.memptr(), &lrwork, iwork.memptr(), &liwork, &info);
    
    return (info == 0);
    }
  #else
    {
    arma_ignore(eigval);
    arma_ignore(eigvec);
    arma_ignore(X);
    arma_stop("eig_sym(): use of LAPACK needs to be enabled");
    return false;
    }
  #endif
  }



//! Eigenvalues and eigenvectors of a general square real matrix using LAPACK.
//! The argument 'side' specifies which eigenvectors should be calculated
//! (see code for mode details).
template<typename T, typename T1>
inline
bool
auxlib::eig_gen
  (
  Col< std::complex<T> >&   eigval,
  Mat<T>&                 l_eigvec,
  Mat<T>&                 r_eigvec,
  const Base<T,T1>&       X,
  const char              side
  )
  {
  arma_extra_debug_sigprint();
  
  #if defined(ARMA_USE_LAPACK)
    {
    char jobvl;
    char jobvr;
    
    switch(side)
      {
      case 'l':  // left
        jobvl = 'V';
        jobvr = 'N';
        break;
        
      case 'r':  // right
        jobvl = 'N';
        jobvr = 'V';
        break;
        
      case 'b':  // both
        jobvl = 'V';
        jobvr = 'V';
        break;
        
      case 'n':  // neither
        jobvl = 'N';
        jobvr = 'N';
        break;
      
      default:
        arma_stop("eig_gen(): parameter 'side' is invalid");
        return false;
      }
    
    Mat<T> A(X.get_ref());
    arma_debug_check( (A.is_square() == false), "eig_gen(): given matrix is not square" );
    
    if(A.is_empty())
      {
      eigval.reset();
      l_eigvec.reset();
      r_eigvec.reset();
      return true;
      }
    
    const uword A_n_rows = A.n_rows;
    
    eigval.set_size(A_n_rows);
    
    l_eigvec.set_size(A_n_rows, A_n_rows);
    r_eigvec.set_size(A_n_rows, A_n_rows);
    
    blas_int N     = blas_int(A_n_rows);
    blas_int lwork = 3 * ( (std::max)(blas_int(1), 4*N) );
    blas_int info  = 0;
    
    podarray<T> work( static_cast<uword>(lwork) );
    
    podarray<T> wr(A_n_rows);
    podarray<T> wi(A_n_rows);
    
    arma_extra_debug_print("lapack::geev()");
    lapack::geev(&jobvl, &jobvr, &N, A.memptr(), &N, wr.memptr(), wi.memptr(), l_eigvec.memptr(), &N, r_eigvec.memptr(), &N, work.memptr(), &lwork, &info);
    
    eigval.set_size(A_n_rows);
    for(uword i=0; i<A_n_rows; ++i)
      {
      eigval[i] = std::complex<T>(wr[i], wi[i]);
      }
    
    return (info == 0);
    }
  #else
    {
    arma_ignore(eigval);
    arma_ignore(l_eigvec);
    arma_ignore(r_eigvec);
    arma_ignore(X);
    arma_ignore(side);
    arma_stop("eig_gen(): use of LAPACK needs to be enabled");
    return false;
    }
  #endif
  }





//! Eigenvalues and eigenvectors of a general square complex matrix using LAPACK
//! The argument 'side' specifies which eigenvectors should be calculated
//! (see code for mode details).
template<typename T, typename T1>
inline
bool
auxlib::eig_gen
  (
  Col< std::complex<T> >&              eigval,
  Mat< std::complex<T> >&            l_eigvec, 
  Mat< std::complex<T> >&            r_eigvec, 
  const Base< std::complex<T>, T1 >& X, 
  const char                         side
  )
  {
  arma_extra_debug_sigprint();
  
  #if defined(ARMA_USE_LAPACK)
    {
    typedef typename std::complex<T> eT;
    
    char jobvl;
    char jobvr;
    
    switch(side)
      {
      case 'l':  // left
        jobvl = 'V';
        jobvr = 'N';
        break;
        
      case 'r':  // right
        jobvl = 'N';
        jobvr = 'V';
        break;
        
      case 'b':  // both
        jobvl = 'V';
        jobvr = 'V';
        break;
        
      case 'n':  // neither
        jobvl = 'N';
        jobvr = 'N';
        break;
      
      default:
        arma_stop("eig_gen(): parameter 'side' is invalid");
        return false;
      }
    
    Mat<eT> A(X.get_ref());
    arma_debug_check( (A.is_square() == false), "eig_gen(): given matrix is not square" );
    
    if(A.is_empty())
      {
      eigval.reset();
      l_eigvec.reset();
      r_eigvec.reset();
      return true;
      }
    
    const uword A_n_rows = A.n_rows;
    
    eigval.set_size(A_n_rows);
    
    l_eigvec.set_size(A_n_rows, A_n_rows);
    r_eigvec.set_size(A_n_rows, A_n_rows);
    
    blas_int N     = blas_int(A_n_rows);
    blas_int lwork = 3 * ( (std::max)(blas_int(1), 2*N) );
    blas_int info  = 0;
    
    podarray<eT>  work( static_cast<uword>(lwork) );
    podarray<T>  rwork( static_cast<uword>(2*N)   );
    
    arma_extra_debug_print("lapack::cx_geev()");
    lapack::cx_geev(&jobvl, &jobvr, &N, A.memptr(), &N, eigval.memptr(), l_eigvec.memptr(), &N, r_eigvec.memptr(), &N, work.memptr(), &lwork, rwork.memptr(), &info);
    
    return (info == 0);
    }
  #else
    {
    arma_ignore(eigval);
    arma_ignore(l_eigvec);
    arma_ignore(r_eigvec);
    arma_ignore(X);
    arma_ignore(side);
    arma_stop("eig_gen(): use of LAPACK needs to be enabled");
    return false;
    }
  #endif
  }



template<typename eT, typename T1>
inline
bool
auxlib::chol(Mat<eT>& out, const Base<eT,T1>& X)
  {
  arma_extra_debug_sigprint();
  
  #if defined(ARMA_USE_LAPACK)
    {
    out = X.get_ref();
    
    arma_debug_check( (out.is_square() == false), "chol(): given matrix is not square" );
    
    if(out.is_empty())
      {
      return true;
      }
    
    const uword out_n_rows = out.n_rows;
    
    char      uplo = 'U';
    blas_int  n    = out_n_rows;
    blas_int  info = 0;
    
    lapack::potrf(&uplo, &n, out.memptr(), &n, &info);
    
    for(uword col=0; col<out_n_rows; ++col)
      {
      eT* colptr = out.colptr(col);
      
      for(uword row=(col+1); row < out_n_rows; ++row)
        {
        colptr[row] = eT(0);
        }
      }
    
    return (info == 0);
    }
  #else
    {
    arma_ignore(out);
    arma_ignore(X);
    
    arma_stop("chol(): use of LAPACK needs to be enabled");
    return false;
    }
  #endif
  }



template<typename eT, typename T1>
inline
bool
auxlib::qr(Mat<eT>& Q, Mat<eT>& R, const Base<eT,T1>& X)
  {
  arma_extra_debug_sigprint();
  
  #if defined(ARMA_USE_LAPACK)
    {
    R = X.get_ref();
    
    const uword R_n_rows = R.n_rows;
    const uword R_n_cols = R.n_cols;
    
    if(R.is_empty())
      {
      Q.eye(R_n_rows, R_n_rows);
      return true;
      }
    
    blas_int m         = static_cast<blas_int>(R_n_rows);
    blas_int n         = static_cast<blas_int>(R_n_cols);
    blas_int lwork     = 0;
    blas_int lwork_min = (std::max)(blas_int(1), (std::max)(m,n));  // take into account requirements of geqrf() _and_ orgqr()/ungqr()
    blas_int k         = (std::min)(m,n);
    blas_int info      = 0;
    
    podarray<eT> tau( static_cast<uword>(k) );
    
    eT        work_query[2];
    blas_int lwork_query = -1;
    
    lapack::geqrf(&m, &n, R.memptr(), &m, tau.memptr(), &work_query[0], &lwork_query, &info);
    
    if(info == 0)
      {
      const blas_int lwork_proposed = static_cast<blas_int>( access::tmp_real(work_query[0]) );
      
      lwork = (lwork_proposed > lwork_min) ? lwork_proposed : lwork_min;
      }
    else
      {
      return false;
      }
    
    podarray<eT> work( static_cast<uword>(lwork) );
    
    lapack::geqrf(&m, &n, R.memptr(), &m, tau.memptr(), work.memptr(), &lwork, &info);
    
    Q.set_size(R_n_rows, R_n_rows);
    
    arrayops::copy( Q.memptr(), R.memptr(), (std::min)(Q.n_elem, R.n_elem) );
    
    //
    // construct R
    
    for(uword col=0; col < R_n_cols; ++col)
      {
      for(uword row=(col+1); row < R_n_rows; ++row)
        {
        R.at(row,col) = eT(0);
        }
      }
    
    
    if( (is_float<eT>::value == true) || (is_double<eT>::value == true) )
      {
      lapack::orgqr(&m, &m, &k, Q.memptr(), &m, tau.memptr(), work.memptr(), &lwork, &info);
      }
    else
    if( (is_supported_complex_float<eT>::value == true) || (is_supported_complex_double<eT>::value == true) )
      {
      lapack::ungqr(&m, &m, &k, Q.memptr(), &m, tau.memptr(), work.memptr(), &lwork, &info);
      }
    
    return (info == 0);
    }
  #else
    {
    arma_ignore(Q);
    arma_ignore(R);
    arma_ignore(X);
    arma_stop("qr(): use of LAPACK needs to be enabled");
    return false;
    }
  #endif
  }



template<typename eT, typename T1>
inline
bool 
auxlib::qr_econ(Mat<eT>& Q, Mat<eT>& R, const Base<eT,T1>& X)
  {
  arma_extra_debug_sigprint();
  
  // This function implements a memory-efficient QR for a non-square X that has dimensions m x n.
  // This basically discards the basis for the null-space.
  // 
  // if m <= n: (use standard routine)
  //     Q[m,m]*R[m,n] = X[m,n]
  //     geqrf Needs A[m,n]: Uses R
  //     orgqr Needs A[m,m]: Uses Q
  // otherwise: (memory-efficient routine)
  //     Q[m,n]*R[n,n] = X[m,n]
  //     geqrf Needs A[m,n]: Uses Q
  //     geqrf Needs A[m,n]: Uses Q
  
  #if defined(ARMA_USE_LAPACK)
    {
    if(is_Mat<T1>::value == true)
      {
      const unwrap<T1>   tmp(X.get_ref());
      const Mat<eT>& M = tmp.M;
      
      if(M.n_rows < M.n_cols)
        {
        return auxlib::qr(Q, R, X);
        }
      }
    
    Q = X.get_ref();
    
    const uword Q_n_rows = Q.n_rows;
    const uword Q_n_cols = Q.n_cols;
    
    if( Q_n_rows <= Q_n_cols )
      {
      return auxlib::qr(Q, R, Q);
      }
    
    if(Q.is_empty())
      {
      Q.set_size(Q_n_rows, 0       );
      R.set_size(0,        Q_n_cols);
      return true;
      }
    
    blas_int m         = static_cast<blas_int>(Q_n_rows);
    blas_int n         = static_cast<blas_int>(Q_n_cols);
    blas_int lwork     = 0;
    blas_int lwork_min = (std::max)(blas_int(1), (std::max)(m,n));  // take into account requirements of geqrf() _and_ orgqr()/ungqr()
    blas_int k         = (std::min)(m,n);
    blas_int info      = 0;
    
    podarray<eT> tau( static_cast<uword>(k) );
    
    eT        work_query[2];
    blas_int lwork_query = -1;
    
    lapack::geqrf(&m, &n, Q.memptr(), &m, tau.memptr(), &work_query[0], &lwork_query, &info);
    
    if(info == 0)
      {
      const blas_int lwork_proposed = static_cast<blas_int>( access::tmp_real(work_query[0]) );
      
      lwork = (lwork_proposed > lwork_min) ? lwork_proposed : lwork_min;
      }
    else
      {
      return false;
      }
    
    podarray<eT> work( static_cast<uword>(lwork) );
    
    lapack::geqrf(&m, &n, Q.memptr(), &m, tau.memptr(), work.memptr(), &lwork, &info);
    
    // Q now has the elements on and above the diagonal of the array
    // contain the min(M,N)-by-N upper trapezoidal matrix Q
    // (Q is upper triangular if m >= n);
    // the elements below the diagonal, with the array TAU,
    // represent the orthogonal matrix Q as a product of min(m,n) elementary reflectors.
    
    R.set_size(Q_n_cols, Q_n_cols);
    
    //
    // construct R
    
    for(uword col=0; col < Q_n_cols; ++col)
      {
      for(uword row=0; row <= col; ++row)
        {
        R.at(row,col) = Q.at(row,col);
        }
      
      for(uword row=(col+1); row < Q_n_cols; ++row)
        {
        R.at(row,col) = eT(0);
        }
      }
    
    if( (is_float<eT>::value == true) || (is_double<eT>::value == true) )
      {
      lapack::orgqr(&m, &n, &k, Q.memptr(), &m, tau.memptr(), work.memptr(), &lwork, &info);
      }
    else
    if( (is_supported_complex_float<eT>::value == true) || (is_supported_complex_double<eT>::value == true) )
      {
      lapack::ungqr(&m, &n, &k, Q.memptr(), &m, tau.memptr(), work.memptr(), &lwork, &info);
      }
    
    return (info == 0);
    }
  #else
    {
    arma_ignore(Q);
    arma_ignore(R);
    arma_ignore(X);
    arma_stop("qr_econ(): use of LAPACK needs to be enabled");
    return false;
    }
  #endif
  }



template<typename eT, typename T1>
inline
bool
auxlib::svd(Col<eT>& S, const Base<eT,T1>& X, uword& X_n_rows, uword& X_n_cols)
  {
  arma_extra_debug_sigprint();
  
  #if defined(ARMA_USE_LAPACK)
    {
    Mat<eT> A(X.get_ref());
    
    X_n_rows = A.n_rows;
    X_n_cols = A.n_cols;
    
    if(A.is_empty())
      {
      S.reset();
      return true;
      }
    
    Mat<eT> U(1, 1);
    Mat<eT> V(1, A.n_cols);
    
    char jobu  = 'N';
    char jobvt = 'N';
    
    blas_int m          = A.n_rows;
    blas_int n          = A.n_cols;
    blas_int min_mn     = (std::min)(m,n);
    blas_int lda        = A.n_rows;
    blas_int ldu        = U.n_rows;
    blas_int ldvt       = V.n_rows;
    blas_int lwork      = 0;
    blas_int lwork_min  = (std::max)( blas_int(1), (std::max)( (3*min_mn + (std::max)(m,n)), 5*min_mn ) );
    blas_int info   = 0;
    
    S.set_size( static_cast<uword>(min_mn) );
    
    eT        work_query[2];
    blas_int lwork_query = -1;
    
    lapack::gesvd<eT>
      (
      &jobu, &jobvt, &m, &n, A.memptr(), &lda, S.memptr(), U.memptr(), &ldu, V.memptr(), &ldvt, &work_query[0], &lwork_query, &info
      );
    
    if(info == 0)
      {
      const blas_int lwork_proposed = static_cast<blas_int>( work_query[0] );
      
      lwork = (lwork_proposed > lwork_min) ? lwork_proposed : lwork_min;
      
      podarray<eT> work( static_cast<uword>(lwork) );
      
      lapack::gesvd<eT>
        (
        &jobu, &jobvt, &m, &n, A.memptr(), &lda, S.memptr(), U.memptr(), &ldu, V.memptr(), &ldvt, work.memptr(), &lwork, &info
        );
      }
    
    return (info == 0);
    }
  #else
    {
    arma_ignore(S);
    arma_ignore(X);
    arma_ignore(X_n_rows);
    arma_ignore(X_n_cols);
    arma_stop("svd(): use of LAPACK needs to be enabled");
    return false;
    }
  #endif
  }



template<typename T, typename T1>
inline
bool
auxlib::svd(Col<T>& S, const Base<std::complex<T>, T1>& X, uword& X_n_rows, uword& X_n_cols)
  {
  arma_extra_debug_sigprint();
  
  #if defined(ARMA_USE_LAPACK)
    {
    typedef std::complex<T> eT;
    
    Mat<eT> A(X.get_ref());
    
    X_n_rows = A.n_rows;
    X_n_cols = A.n_cols;
    
    if(A.is_empty())
      {
      S.reset();
      return true;
      }
    
    Mat<eT> U(1, 1);
    Mat<eT> V(1, A.n_cols);
    
    char jobu  = 'N';
    char jobvt = 'N';
    
    blas_int  m      = A.n_rows;
    blas_int  n      = A.n_cols;
    blas_int  min_mn = (std::min)(m,n);
    blas_int  lda    = A.n_rows;
    blas_int  ldu    = U.n_rows;
    blas_int  ldvt   = V.n_rows;
    blas_int  lwork  = 3 * ( (std::max)(blas_int(1), 2*min_mn+(std::max)(m,n) ) );
    blas_int  info   = 0;
    
    S.set_size( static_cast<uword>(min_mn) );
    
    podarray<eT>   work( static_cast<uword>(lwork   ) );
    podarray< T>  rwork( static_cast<uword>(5*min_mn) );
    
    // let gesvd_() calculate the optimum size of the workspace
    blas_int lwork_tmp = -1;
    
    lapack::cx_gesvd<T>
      (
      &jobu, &jobvt,
      &m, &n,
      A.memptr(), &lda,
      S.memptr(),
      U.memptr(), &ldu,
      V.memptr(), &ldvt,
      work.memptr(), &lwork_tmp,
      rwork.memptr(),
      &info
      );
    
    if(info == 0)
      {
      blas_int proposed_lwork = static_cast<blas_int>(real(work[0]));
      if(proposed_lwork > lwork)
        {
        lwork = proposed_lwork;
        work.set_size( static_cast<uword>(lwork) );
        }
      
      lapack::cx_gesvd<T>
        (
        &jobu, &jobvt,
        &m, &n,
        A.memptr(), &lda,
        S.memptr(),
        U.memptr(), &ldu,
        V.memptr(), &ldvt,
        work.memptr(), &lwork,
        rwork.memptr(),
        &info
        );
      }
        
    return (info == 0);
    }
  #else
    {
    arma_ignore(S);
    arma_ignore(X);
    arma_ignore(X_n_rows);
    arma_ignore(X_n_cols);

    arma_stop("svd(): use of LAPACK needs to be enabled");
    return false;
    }
  #endif
  }



template<typename eT, typename T1>
inline
bool
auxlib::svd(Col<eT>& S, const Base<eT,T1>& X)
  {
  arma_extra_debug_sigprint();
  
  uword junk;
  return auxlib::svd(S, X, junk, junk);
  }



template<typename T, typename T1>
inline
bool
auxlib::svd(Col<T>& S, const Base<std::complex<T>, T1>& X)
  {
  arma_extra_debug_sigprint();
  
  uword junk;
  return auxlib::svd(S, X, junk, junk);
  }



template<typename eT, typename T1>
inline
bool
auxlib::svd(Mat<eT>& U, Col<eT>& S, Mat<eT>& V, const Base<eT,T1>& X)
  {
  arma_extra_debug_sigprint();
  
  #if defined(ARMA_USE_LAPACK)
    {
    Mat<eT> A(X.get_ref());
    
    if(A.is_empty())
      {
      U.eye(A.n_rows, A.n_rows);
      S.reset();
      V.eye(A.n_cols, A.n_cols);
      return true;
      }
    
    U.set_size(A.n_rows, A.n_rows);
    V.set_size(A.n_cols, A.n_cols);
    
    char jobu  = 'A';
    char jobvt = 'A';
    
    blas_int  m          = blas_int(A.n_rows);
    blas_int  n          = blas_int(A.n_cols);
    blas_int  min_mn     = (std::min)(m,n);
    blas_int  lda        = blas_int(A.n_rows);
    blas_int  ldu        = blas_int(U.n_rows);
    blas_int  ldvt       = blas_int(V.n_rows);
    blas_int  lwork_min  = (std::max)( blas_int(1), (std::max)( (3*min_mn + (std::max)(m,n)), 5*min_mn ) );
    blas_int  lwork      = 0;
    blas_int  info       = 0;
    
    S.set_size( static_cast<uword>(min_mn) );
    
    // let gesvd_() calculate the optimum size of the workspace
    eT        work_query[2];
    blas_int lwork_query = -1;
    
    lapack::gesvd<eT>
      (
      &jobu, &jobvt, &m, &n, A.memptr(), &lda, S.memptr(), U.memptr(), &ldu, V.memptr(), &ldvt, &work_query[0], &lwork_query, &info
      );
    
    if(info == 0)
      {
      const blas_int lwork_proposed = static_cast<blas_int>( work_query[0] );
      
      lwork = (lwork_proposed > lwork_min) ? lwork_proposed : lwork_min;
      
      podarray<eT> work( static_cast<uword>(lwork) );
      
      lapack::gesvd<eT>
        (
        &jobu, &jobvt, &m, &n, A.memptr(), &lda, S.memptr(), U.memptr(), &ldu, V.memptr(), &ldvt, work.memptr(), &lwork, &info
        );
      
      op_strans::apply(V,V);  // op_strans will work out that an in-place transpose can be done
      }
    
    return (info == 0);
    }
  #else
    {
    arma_ignore(U);
    arma_ignore(S);
    arma_ignore(V);
    arma_ignore(X);
    arma_stop("svd(): use of LAPACK needs to be enabled");
    return false;
    }
  #endif
  }



template<typename T, typename T1>
inline
bool
auxlib::svd(Mat< std::complex<T> >& U, Col<T>& S, Mat< std::complex<T> >& V, const Base< std::complex<T>, T1>& X)
  {
  arma_extra_debug_sigprint();
  
  #if defined(ARMA_USE_LAPACK)
    {
    typedef std::complex<T> eT;
    
    Mat<eT> A(X.get_ref());
    
    if(A.is_empty())
      {
      U.eye(A.n_rows, A.n_rows);
      S.reset();
      V.eye(A.n_cols, A.n_cols);
      return true;
      }
    
    U.set_size(A.n_rows, A.n_rows);
    V.set_size(A.n_cols, A.n_cols);
    
    char jobu  = 'A';
    char jobvt = 'A';
    
    blas_int  m      = blas_int(A.n_rows);
    blas_int  n      = blas_int(A.n_cols);
    blas_int  min_mn = (std::min)(m,n);
    blas_int  lda    = blas_int(A.n_rows);
    blas_int  ldu    = blas_int(U.n_rows);
    blas_int  ldvt   = blas_int(V.n_rows);
    blas_int  lwork  = 3 * ( (std::max)(blas_int(1), 2*min_mn + (std::max)(m,n) ) );
    blas_int  info   = 0;
    
    S.set_size( static_cast<uword>(min_mn) );
    
    podarray<eT>  work( static_cast<uword>(lwork   ) );
    podarray<T>  rwork( static_cast<uword>(5*min_mn) );
    
    // let gesvd_() calculate the optimum size of the workspace
    blas_int lwork_tmp = -1;
    lapack::cx_gesvd<T>
     (
     &jobu, &jobvt,
     &m, &n,
     A.memptr(), &lda,
     S.memptr(),
     U.memptr(), &ldu,
     V.memptr(), &ldvt,
     work.memptr(), &lwork_tmp,
     rwork.memptr(),
     &info
     );
    
    if(info == 0)
      {
      blas_int proposed_lwork = static_cast<blas_int>(real(work[0]));
      
      if(proposed_lwork > lwork)
        {
        lwork = proposed_lwork;
        work.set_size( static_cast<uword>(lwork) );
        }
      
      lapack::cx_gesvd<T>
        (
        &jobu, &jobvt,
        &m, &n,
        A.memptr(), &lda,
        S.memptr(),
        U.memptr(), &ldu,
        V.memptr(), &ldvt,
        work.memptr(), &lwork,
        rwork.memptr(),
        &info
        );
      
      op_htrans::apply(V,V);  // op_htrans will work out that an in-place transpose can be done
      }
    
    return (info == 0);
    }
  #else
    {
    arma_ignore(U);
    arma_ignore(S);
    arma_ignore(V);
    arma_ignore(X);
    arma_stop("svd(): use of LAPACK needs to be enabled");
    return false;
    }
  #endif
  }



template<typename eT, typename T1>
inline
bool
auxlib::svd_econ(Mat<eT>& U, Col<eT>& S, Mat<eT>& V, const Base<eT,T1>& X, const char mode)
  {
  arma_extra_debug_sigprint();
  
  #if defined(ARMA_USE_LAPACK)
    {
    Mat<eT> A(X.get_ref());
    
    blas_int m      = blas_int(A.n_rows);
    blas_int n      = blas_int(A.n_cols);
    blas_int min_mn = (std::min)(m,n);
    blas_int lda    = blas_int(A.n_rows);
    
    S.set_size( static_cast<uword>(min_mn) );
    
    blas_int ldu  = 0;
    blas_int ldvt = 0;
    
    char jobu;
    char jobvt;
    
    switch(mode)
      {
      case 'l':
        jobu  = 'S';
        jobvt = 'N';
        
        ldu  = m;
        ldvt = 1;
        
        U.set_size( static_cast<uword>(ldu), static_cast<uword>(min_mn) );
        V.reset();
        
        break;
      
      
      case 'r':
        jobu  = 'N';
        jobvt = 'S';
        
        ldu = 1;
        ldvt = (std::min)(m,n);
        
        U.reset();
        V.set_size( static_cast<uword>(ldvt), static_cast<uword>(n) );
        
        break;
      
      
      case 'b':
        jobu  = 'S';
        jobvt = 'S';
        
        ldu  = m;
        ldvt = (std::min)(m,n);
        
        U.set_size( static_cast<uword>(ldu),  static_cast<uword>(min_mn) );
        V.set_size( static_cast<uword>(ldvt), static_cast<uword>(n     ) );
        
        break;
      
      
      default:
        U.reset();
        S.reset();
        V.reset();
        return false;
      }
    
    
    if(A.is_empty())
      {
      U.eye();
      S.reset();
      V.eye();
      return true;
      }
    
    
    blas_int lwork = 3 * ( (std::max)(blas_int(1), (std::max)( (3*min_mn + (std::max)(m,n)), 5*min_mn ) ) );
    blas_int info  = 0;
    
    
    podarray<eT> work( static_cast<uword>(lwork) );
    
    // let gesvd_() calculate the optimum size of the workspace
    blas_int lwork_tmp = -1;
    
    lapack::gesvd<eT>
      (
      &jobu, &jobvt,
      &m, &n,
      A.memptr(), &lda,
      S.memptr(),
      U.memptr(), &ldu,
      V.memptr(), &ldvt,
      work.memptr(), &lwork_tmp,
      &info
      );
    
    if(info == 0)
      {
      blas_int proposed_lwork = static_cast<blas_int>(work[0]);
      if(proposed_lwork > lwork)
        {
        lwork = proposed_lwork;
        work.set_size( static_cast<uword>(lwork) );
        }
      
      lapack::gesvd<eT>
        (
        &jobu, &jobvt,
        &m, &n,
        A.memptr(), &lda,
        S.memptr(),
        U.memptr(), &ldu,
        V.memptr(), &ldvt,
        work.memptr(), &lwork,
        &info
        );
      
      op_strans::apply(V,V);  // op_strans will work out that an in-place transpose can be done
      }
    
    return (info == 0);
    }
  #else
    {
    arma_ignore(U);
    arma_ignore(S);
    arma_ignore(V);
    arma_ignore(X);
    arma_ignore(mode);
    arma_stop("svd(): use of LAPACK needs to be enabled");
    return false;
    }
  #endif
  }



template<typename T, typename T1>
inline
bool
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)
  {
  arma_extra_debug_sigprint();
  
  #if defined(ARMA_USE_LAPACK)
    {
    typedef std::complex<T> eT;
    
    Mat<eT> A(X.get_ref());
    
    blas_int m      = blas_int(A.n_rows);
    blas_int n      = blas_int(A.n_cols);
    blas_int min_mn = (std::min)(m,n);
    blas_int lda    = blas_int(A.n_rows);
    
    S.set_size( static_cast<uword>(min_mn) );
    
    blas_int ldu  = 0;
    blas_int ldvt = 0;
    
    char jobu;
    char jobvt;
    
    switch(mode)
      {
      case 'l':
        jobu  = 'S';
        jobvt = 'N';
        
        ldu  = m;
        ldvt = 1;
        
        U.set_size( static_cast<uword>(ldu), static_cast<uword>(min_mn) );
        V.reset();
        
        break;
      
      
      case 'r':
        jobu  = 'N';
        jobvt = 'S';
        
        ldu  = 1;
        ldvt = (std::min)(m,n);
        
        U.reset();
        V.set_size( static_cast<uword>(ldvt), static_cast<uword>(n) );
        
        break;
      
      
      case 'b':
        jobu  = 'S';
        jobvt = 'S';
        
        ldu  = m;
        ldvt = (std::min)(m,n);
        
        U.set_size( static_cast<uword>(ldu),  static_cast<uword>(min_mn) );
        V.set_size( static_cast<uword>(ldvt), static_cast<uword>(n)      );
        
        break;
      
      
      default:
        U.reset();
        S.reset();
        V.reset();
        return false;
      }
    
    
    if(A.is_empty())
      {
      U.eye();
      S.reset();
      V.eye();
      return true;
      }
    
    
    blas_int lwork = 3 * ( (std::max)(blas_int(1), (std::max)( (3*min_mn + (std::max)(m,n)), 5*min_mn ) ) );
    blas_int info  = 0;
    
    
    podarray<eT>  work( static_cast<uword>(lwork   ) );
    podarray<T>  rwork( static_cast<uword>(5*min_mn) );
    
    // let gesvd_() calculate the optimum size of the workspace
    blas_int lwork_tmp = -1;
    
    lapack::cx_gesvd<T>
      (
      &jobu, &jobvt,
      &m, &n,
      A.memptr(), &lda,
      S.memptr(),
      U.memptr(), &ldu,
      V.memptr(), &ldvt,
      work.memptr(), &lwork_tmp,
      rwork.memptr(),
      &info
      );
    
    if(info == 0)
      {
      blas_int proposed_lwork = static_cast<blas_int>(real(work[0]));
      if(proposed_lwork > lwork)
        {
        lwork = proposed_lwork;
        work.set_size( static_cast<uword>(lwork) );
        }
      
      lapack::cx_gesvd<T>
        (
        &jobu, &jobvt,
        &m, &n,
        A.memptr(), &lda,
        S.memptr(),
        U.memptr(), &ldu,
        V.memptr(), &ldvt,
        work.memptr(), &lwork,
        rwork.memptr(),
        &info
        );
      
      op_htrans::apply(V,V);  // op_strans will work out that an in-place transpose can be done
      }
    
    return (info == 0);
    }
  #else
    {
    arma_ignore(U);
    arma_ignore(S);
    arma_ignore(V);
    arma_ignore(X);
    arma_ignore(mode);
    arma_stop("svd(): use of LAPACK needs to be enabled");
    return false;
    }
  #endif
  }



template<typename eT, typename T1>
inline
bool
auxlib::svd_dc(Mat<eT>& U, Col<eT>& S, Mat<eT>& V, const Base<eT,T1>& X)
  {
  arma_extra_debug_sigprint();
  
  #if defined(ARMA_USE_LAPACK)
    {
    Mat<eT> A(X.get_ref());
    
    if(A.is_empty())
      {
      U.eye(A.n_rows, A.n_rows);
      S.reset();
      V.eye(A.n_cols, A.n_cols);
      return true;
      }
    
    U.set_size(A.n_rows, A.n_rows);
    V.set_size(A.n_cols, A.n_cols);
    
    char jobz = 'A';
    
    blas_int  m      = blas_int(A.n_rows);
    blas_int  n      = blas_int(A.n_cols);
    blas_int  min_mn = (std::min)(m,n);
    blas_int  lda    = blas_int(A.n_rows);
    blas_int  ldu    = blas_int(U.n_rows);
    blas_int  ldvt   = blas_int(V.n_rows);
    blas_int  lwork  = 3 * ( 3*min_mn*min_mn + (std::max)( (std::max)(m,n), 4*min_mn*min_mn + 4*min_mn ) );
    blas_int  info   = 0;
    
    S.set_size( static_cast<uword>(min_mn) );
    
    podarray<eT>        work( static_cast<uword>(lwork   ) );
    podarray<blas_int> iwork( static_cast<uword>(8*min_mn) );
    
    lapack::gesdd<eT>
      (
      &jobz, &m, &n,
      A.memptr(), &lda, S.memptr(), U.memptr(), &ldu, V.memptr(), &ldvt,
      work.memptr(), &lwork, iwork.memptr(), &info
      );
    
    op_strans::apply(V,V);  // op_strans will work out that an in-place transpose can be done
    
    return (info == 0);
    }
  #else
    {
    arma_ignore(U);
    arma_ignore(S);
    arma_ignore(V);
    arma_ignore(X);
    arma_stop("svd(): use of LAPACK needs to be enabled");
    return false;
    }
  #endif
  }



template<typename T, typename T1>
inline
bool
auxlib::svd_dc(Mat< std::complex<T> >& U, Col<T>& S, Mat< std::complex<T> >& V, const Base< std::complex<T>, T1>& X)
  {
  arma_extra_debug_sigprint();
  
  #if defined(ARMA_USE_LAPACK)
    {
    typedef std::complex<T> eT;
    
    Mat<eT> A(X.get_ref());
    
    if(A.is_empty())
      {
      U.eye(A.n_rows, A.n_rows);
      S.reset();
      V.eye(A.n_cols, A.n_cols);
      return true;
      }
    
    U.set_size(A.n_rows, A.n_rows);
    V.set_size(A.n_cols, A.n_cols);
    
    char jobz = 'A';
    
    blas_int  m      = blas_int(A.n_rows);
    blas_int  n      = blas_int(A.n_cols);
    blas_int  min_mn = (std::min)(m,n);
    blas_int  lda    = blas_int(A.n_rows);
    blas_int  ldu    = blas_int(U.n_rows);
    blas_int  ldvt   = blas_int(V.n_rows);
    blas_int  lwork  = 3 * (min_mn*min_mn + 2*min_mn + (std::max)(m,n));
    blas_int  info   = 0;
    
    S.set_size( static_cast<uword>(min_mn) );
    
    podarray<eT>        work( static_cast<uword>(lwork                     ) );
    podarray<T>        rwork( static_cast<uword>(5*min_mn*min_mn + 7*min_mn) );
    podarray<blas_int> iwork( static_cast<uword>(8*min_mn                  ) );
    
    lapack::cx_gesdd<T>
      (
      &jobz, &m, &n,
      A.memptr(), &lda, S.memptr(), U.memptr(), &ldu, V.memptr(), &ldvt,
      work.memptr(), &lwork, rwork.memptr(), iwork.memptr(), &info
      );
    
    op_htrans::apply(V,V);  // op_htrans will work out that an in-place transpose can be done
    
    return (info == 0);
    }
  #else
    {
    arma_ignore(U);
    arma_ignore(S);
    arma_ignore(V);
    arma_ignore(X);
    arma_stop("svd(): use of LAPACK needs to be enabled");
    return false;
    }
  #endif
  }



//! Solve a system of linear equations.
//! Assumes that A.n_rows = A.n_cols and B.n_rows = A.n_rows
template<typename eT, typename T1>
inline
bool
auxlib::solve(Mat<eT>& out, Mat<eT>& A, const Base<eT,T1>& X, const bool slow)
  {
  arma_extra_debug_sigprint();
  
  bool status = false;
  
  const uword A_n_rows = A.n_rows;
  
  if( (A_n_rows <= 4) && (slow == false) )
    {
    Mat<eT> A_inv;
    
    status = auxlib::inv_noalias_tinymat(A_inv, A, A_n_rows);
    
    if(status == true)
      {
      const unwrap_check<T1> Y( X.get_ref(), out );
      const Mat<eT>& B     = Y.M;
      
      const uword B_n_rows = B.n_rows;
      const uword B_n_cols = B.n_cols;
      
      arma_debug_check( (A_n_rows != B_n_rows), "solve(): number of rows in the given objects must be the same" );
      
      if(A.is_empty() || B.is_empty())
        {
        out.zeros(A.n_cols, B_n_cols);
        return true;
        }
      
      out.set_size(A_n_rows, B_n_cols);
      
      gemm_emul<false,false,false,false>::apply(out, A_inv, B);
      
      return true;
      }
    }
  
  if( (A_n_rows > 4) || (status == false) )
    {
    out = X.get_ref();
    
    const uword B_n_rows = out.n_rows;
    const uword B_n_cols = out.n_cols;
      
    arma_debug_check( (A_n_rows != B_n_rows), "solve(): number of rows in the given objects must be the same" );
      
    if(A.is_empty() || out.is_empty())
      {
      out.zeros(A.n_cols, B_n_cols);
      return true;
      }
    
    #if defined(ARMA_USE_ATLAS)
      {
      podarray<int> ipiv(A_n_rows + 2);  // +2 for paranoia: old versions of Atlas might be trashing memory
      
      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);
      
      return (info == 0);
      }
    #elif defined(ARMA_USE_LAPACK)
      {
      blas_int n    = blas_int(A_n_rows);  // assuming A is square
      blas_int lda  = blas_int(A_n_rows);
      blas_int ldb  = blas_int(A_n_rows);
      blas_int nrhs = blas_int(B_n_cols);
      blas_int info = 0;
      
      podarray<blas_int> ipiv(A_n_rows + 2);  // +2 for paranoia: some versions of Lapack might be trashing memory
      
      arma_extra_debug_print("lapack::gesv()");
      lapack::gesv<eT>(&n, &nrhs, A.memptr(), &lda, ipiv.memptr(), out.memptr(), &ldb, &info);
      
      arma_extra_debug_print("lapack::gesv() -- finished");
      
      return (info == 0);
      }
    #else
      {
      arma_stop("solve(): use of ATLAS or LAPACK needs to be enabled");
      return false;
      }
    #endif
    }
  
  return true;
  }



//! Solve an over-determined system.
//! Assumes that A.n_rows > A.n_cols and B.n_rows = A.n_rows
template<typename eT, typename T1>
inline
bool
auxlib::solve_od(Mat<eT>& out, Mat<eT>& A, const Base<eT,T1>& X)
  {
  arma_extra_debug_sigprint();
  
  #if defined(ARMA_USE_LAPACK)
    {
    Mat<eT> tmp = X.get_ref();
    
    const uword A_n_rows = A.n_rows;
    const uword A_n_cols = A.n_cols;
    
    const uword B_n_rows = tmp.n_rows;
    const uword B_n_cols = tmp.n_cols;
      
    arma_debug_check( (A_n_rows != B_n_rows), "solve(): number of rows in the given objects must be the same" );
    
    out.set_size(A_n_cols, B_n_cols);
    
    if(A.is_empty() || tmp.is_empty())
      {
      out.zeros();
      return true;
      }
    
    char trans = 'N';
    
    blas_int  m     = blas_int(A_n_rows);
    blas_int  n     = blas_int(A_n_cols);
    blas_int  lda   = blas_int(A_n_rows);
    blas_int  ldb   = blas_int(A_n_rows);
    blas_int  nrhs  = blas_int(B_n_cols);
    blas_int  lwork = 3 * ( (std::max)(blas_int(1), n + (std::max)(n, nrhs)) );
    blas_int  info  = 0;
    
    podarray<eT> work( static_cast<uword>(lwork) );
    
    // NOTE: the dgels() function in the lapack library supplied by ATLAS 3.6 seems to have problems
    arma_extra_debug_print("lapack::gels()");
    lapack::gels<eT>( &trans, &m, &n, &nrhs, A.memptr(), &lda, tmp.memptr(), &ldb, work.memptr(), &lwork, &info );
    
    arma_extra_debug_print("lapack::gels() -- finished");
    
    for(uword col=0; col<B_n_cols; ++col)
      {
      arrayops::copy( out.colptr(col), tmp.colptr(col), A_n_cols );
      }
    
    return (info == 0);
    }
  #else
    {
    arma_ignore(out);
    arma_ignore(A);
    arma_ignore(X);
    arma_stop("solve(): use of LAPACK needs to be enabled");
    return false;
    }
  #endif
  }



//! Solve an under-determined system.
//! Assumes that A.n_rows < A.n_cols and B.n_rows = A.n_rows
template<typename eT, typename T1>
inline
bool
auxlib::solve_ud(Mat<eT>& out, Mat<eT>& A, const Base<eT,T1>& X)
  {
  arma_extra_debug_sigprint();
  
  // TODO: this function provides the same results as Octave 3.4.2.
  // TODO: however, these results are different than Matlab 7.12.0.635.
  // TODO: figure out whether both Octave and Matlab are correct, or only one of them
  
  #if defined(ARMA_USE_LAPACK)
    {
    const unwrap<T1>   Y( X.get_ref() );
    const Mat<eT>& B = Y.M;
    
    const uword A_n_rows = A.n_rows;
    const uword A_n_cols = A.n_cols;
    
    const uword B_n_rows = B.n_rows;
    const uword B_n_cols = B.n_cols;
    
    arma_debug_check( (A_n_rows != B_n_rows), "solve(): number of rows in the given objects must be the same" );
    
    // B could be an alias of "out", hence we need to check whether B is empty before setting the size of "out"
    if(A.is_empty() || B.is_empty())
      {
      out.zeros(A_n_cols, B_n_cols);
      return true;
      }
    
    char trans = 'N';
    
    blas_int  m     = blas_int(A_n_rows);
    blas_int  n     = blas_int(A_n_cols);
    blas_int  lda   = blas_int(A_n_rows);
    blas_int  ldb   = blas_int(A_n_cols);
    blas_int  nrhs  = blas_int(B_n_cols);
    blas_int  lwork = 3 * ( (std::max)(blas_int(1), m + (std::max)(m,nrhs)) );
    blas_int  info  = 0;
    
    Mat<eT> tmp(A_n_cols, B_n_cols);
    tmp.zeros();
    
    for(uword col=0; col<B_n_cols; ++col)
      {
      eT* tmp_colmem = tmp.colptr(col);
      
      arrayops::copy( tmp_colmem, B.colptr(col), B_n_rows );
      
      for(uword row=B_n_rows; row<A_n_cols; ++row)
        {
        tmp_colmem[row] = eT(0);
        }
      }
    
    podarray<eT> work( static_cast<uword>(lwork) );
    
    // NOTE: the dgels() function in the lapack library supplied by ATLAS 3.6 seems to have problems
    arma_extra_debug_print("lapack::gels()");
    lapack::gels<eT>( &trans, &m, &n, &nrhs, A.memptr(), &lda, tmp.memptr(), &ldb, work.memptr(), &lwork, &info );
    
    arma_extra_debug_print("lapack::gels() -- finished");
    
    out.set_size(A_n_cols, B_n_cols);
    
    for(uword col=0; col<B_n_cols; ++col)
      {
      arrayops::copy( out.colptr(col), tmp.colptr(col), A_n_cols );
      }
    
    return (info == 0);
    }
  #else
    {
    arma_ignore(out);
    arma_ignore(A);
    arma_ignore(X);
    arma_stop("solve(): use of LAPACK needs to be enabled");
    return false;
    }
  #endif
  }



//
// solve_tr

template<typename eT>
inline
bool
auxlib::solve_tr(Mat<eT>& out, const Mat<eT>& A, const Mat<eT>& B, const uword layout)
  {
  arma_extra_debug_sigprint();
  
  #if defined(ARMA_USE_LAPACK)
    {
    if(A.is_empty() || B.is_empty())
      {
      out.zeros(A.n_cols, B.n_cols);
      return true;
      }
    
    out = B;
    
    char     uplo  = (layout == 0) ? 'U' : 'L';
    char     trans = 'N';
    char     diag  = 'N';
    blas_int n     = blas_int(A.n_rows);
    blas_int nrhs  = blas_int(B.n_cols);
    blas_int info  = 0;
    
    lapack::trtrs<eT>(&uplo, &trans, &diag, &n, &nrhs, A.memptr(), &n, out.memptr(), &n, &info);
    
    return (info == 0);
    }
  #else
    {
    arma_ignore(out);
    arma_ignore(A);
    arma_ignore(B);
    arma_ignore(layout);
    arma_stop("solve(): use of LAPACK needs to be enabled");
    return false;
    }
  #endif
  }



//
// Schur decomposition

template<typename eT>
inline
bool
auxlib::schur_dec(Mat<eT>& Z, Mat<eT>& T, const Mat<eT>& A)
  {
  arma_extra_debug_sigprint();
  
  #if defined(ARMA_USE_LAPACK)
    {
    arma_debug_check( (A.is_square() == false), "schur_dec(): given matrix is not square" );
    
    if(A.is_empty())
      {
      Z.reset();
      T.reset();
      return true;
      }
    
    const uword A_n_rows = A.n_rows;
    
    Z.set_size(A_n_rows, A_n_rows);
    T = A;
    
    char    jobvs    = 'V';                // get Schur vectors (Z)
    char     sort    = 'N';                // do not sort eigenvalues/vectors
    blas_int* select = 0;                  // pointer to sorting function
    blas_int    n    = blas_int(A_n_rows);
    blas_int sdim    = 0;                  // output for sorting
    blas_int lwork   = 3 * ( (std::max)(blas_int(1), 3*n) );
    blas_int info    = 0;
    
    podarray<eT>       work( static_cast<uword>(lwork) );
    podarray<blas_int> bwork(A_n_rows);
    
    podarray<eT> wr(A_n_rows);             // output for eigenvalues
    podarray<eT> wi(A_n_rows);             // output for eigenvalues
    
    lapack::gees(&jobvs, &sort, select, &n, T.memptr(), &n, &sdim, wr.memptr(), wi.memptr(), Z.memptr(), &n, work.memptr(), &lwork, bwork.memptr(), &info);
    
    return (info == 0);
    }
  #else
    {
    arma_ignore(Z);
    arma_ignore(T);
    arma_ignore(A);
    
    arma_stop("schur_dec(): use of LAPACK needs to be enabled");
    return false;
    }
  #endif
  }



template<typename cT>
inline
bool
auxlib::schur_dec(Mat<std::complex<cT> >& Z, Mat<std::complex<cT> >& T, const Mat<std::complex<cT> >& A)
  {
  arma_extra_debug_sigprint();
  
  #if defined(ARMA_USE_LAPACK)
    {
    arma_debug_check( (A.is_square() == false), "schur_dec(): matrix A is not square" );
    
    if(A.is_empty())
      {
      Z.reset();
      T.reset();
      return true;
      }
    
    typedef std::complex<cT> eT;
    
    const uword A_n_rows = A.n_rows;
    
    Z.set_size(A_n_rows, A_n_rows);
    T = A;
    
    char        jobvs = 'V';                // get Schur vectors (Z)
    char         sort = 'N';                // do not sort eigenvalues/vectors
    blas_int*  select = 0;                  // pointer to sorting function
    blas_int        n = blas_int(A_n_rows);
    blas_int     sdim = 0;                  // output for sorting
    blas_int lwork    = 3 * ( (std::max)(blas_int(1), 2*n) );
    blas_int info     = 0;
    
    podarray<eT>       work( static_cast<uword>(lwork) );
    podarray<blas_int> bwork(A_n_rows);
    
    podarray<eT>     w(A_n_rows);           // output for eigenvalues
    podarray<cT> rwork(A_n_rows);
    
    lapack::cx_gees(&jobvs, &sort, select, &n, T.memptr(), &n, &sdim, w.memptr(), Z.memptr(), &n, work.memptr(), &lwork, rwork.memptr(), bwork.memptr(), &info);
    
    return (info == 0);
    }
  #else
    {
    arma_ignore(Z);
    arma_ignore(T);
    arma_ignore(A);
    
    arma_stop("schur_dec(): use of LAPACK needs to be enabled");
    return false;
    }
  #endif
  }



//
// syl (solution of the Sylvester equation AX + XB = C)

template<typename eT>
inline
bool
auxlib::syl(Mat<eT>& X, const Mat<eT>& A, const Mat<eT>& B, const Mat<eT>& C)
  {
  arma_extra_debug_sigprint();
  
  arma_debug_check
    (
    (A.is_square() == false) || (B.is_square() == false),
    "syl(): given matrix is not square"
    );
    
  arma_debug_check
    (
    (C.n_rows != A.n_rows) || (C.n_cols != B.n_cols),
    "syl(): matrices are not conformant"
    );
  
  if(A.is_empty() || B.is_empty() || C.is_empty())
    {
    X.reset();
    return true;
    }
  
  #if defined(ARMA_USE_LAPACK)
    {
    Mat<eT> Z1, Z2, T1, T2;
    
    const bool status_sd1 = auxlib::schur_dec(Z1, T1, A);
    const bool status_sd2 = auxlib::schur_dec(Z2, T2, B);
    
    if( (status_sd1 == false) || (status_sd2 == false) )
      {
      return false;
      }
    
    char     trana = 'N';
    char     tranb = 'N';
    blas_int  isgn = +1;
    blas_int     m = blas_int(T1.n_rows);
    blas_int     n = blas_int(T2.n_cols);
    
    eT       scale = eT(0);
    blas_int  info = 0;
    
    Mat<eT> Y = trans(Z1) * C * Z2;
    
    lapack::trsyl<eT>(&trana, &tranb, &isgn, &m, &n, T1.memptr(), &m, T2.memptr(), &n, Y.memptr(), &m, &scale, &info);
    
    //Y /= scale;
    Y /= (-scale);
    
    X = Z1 * Y * trans(Z2);
    
    return (info >= 0);
    }
  #else
    {
    arma_stop("syl(): use of LAPACK needs to be enabled");
    return false;
    }
  #endif
  }
  
  
  
//
// lyap (solution of the continuous Lyapunov equation AX + XA^H + Q = 0)

template<typename eT>
inline
bool
auxlib::lyap(Mat<eT>& X, const Mat<eT>& A, const Mat<eT>& Q)
  {
  arma_extra_debug_sigprint();
  
  arma_debug_check( (A.is_square() == false), "lyap(): matrix A is not square");
  arma_debug_check( (Q.is_square() == false), "lyap(): matrix Q is not square");
  arma_debug_check( (A.n_rows != Q.n_rows),   "lyap(): matrices A and Q have different dimensions");
  
  Mat<eT> htransA;
  op_htrans::apply_noalias(htransA, A);
  
  const Mat<eT> mQ = -Q;
  
  return auxlib::syl(X, A, htransA, mQ);
  }



//
// dlyap (solution of the discrete Lyapunov equation AXA^H - X + Q = 0)

template<typename eT>
inline
bool
auxlib::dlyap(Mat<eT>& X, const Mat<eT>& A, const Mat<eT>& Q)
  {
  arma_extra_debug_sigprint();
  
  arma_debug_check( (A.is_square() == false), "dlyap(): matrix A is not square");
  arma_debug_check( (Q.is_square() == false), "dlyap(): matrix Q is not square");
  arma_debug_check( (A.n_rows != Q.n_rows),   "dlyap(): matrices A and Q have different dimensions");
  
  const Col<eT> vecQ = reshape(Q, Q.n_elem, 1);
  
  const Mat<eT> M = eye< Mat<eT> >(Q.n_elem, Q.n_elem) - kron(conj(A), A);
  
  Col<eT> vecX;
  
  const bool status = solve(vecX, M, vecQ);
  
  if(status == true)
    {
    X = reshape(vecX, Q.n_rows, Q.n_cols);
    return true;
    }
  else
    {
    X.reset();
    return false;
    }
  }



//! @}