op_inv_meat.hpp

Go to the documentation of this file.
00001 // Copyright (C) 2009 NICTA
00002 // 
00003 // Authors:
00004 // - Conrad Sanderson (conradsand at ieee dot org)
00005 // 
00006 // This file is part of the Armadillo C++ library.
00007 // It is provided without any warranty of fitness
00008 // for any purpose. You can redistribute this file
00009 // and/or modify it under the terms of the GNU
00010 // Lesser General Public License (LGPL) as published
00011 // by the Free Software Foundation, either version 3
00012 // of the License or (at your option) any later version.
00013 // (see http://www.opensource.org/licenses for more info)
00014 
00015 
00016 //! \addtogroup op_inv
00017 //! @{
00018 
00019 
00020 //! immediate inverse of a matrix, storing the result in a dense matrix
00021 template<typename eT>
00022 inline
00023 void
00024 op_inv::apply(Mat<eT>& out, const Mat<eT>& A)
00025   {
00026   arma_extra_debug_sigprint();
00027   
00028   // no need to check for aliasing, due to:
00029   // - auxlib::inv() copies A to out before inversion
00030   // - for 2x2 and 3x3 matrices the code is alias safe
00031   
00032   arma_debug_check( !A.is_square(), "op_inv::apply(): matrix must be square" );
00033   
00034   if(&out != &A)
00035     {
00036     auxlib::inv_noalias(out, A);
00037     }
00038   else
00039     {
00040     auxlib::inv_inplace(out);
00041     }
00042   
00043   }
00044 
00045 
00046 
00047 //! immediate inverse of T1, storing the result in a dense matrix
00048 template<typename T1>
00049 inline
00050 void
00051 op_inv::apply(Mat<typename T1::elem_type>& out, const Op<T1,op_inv>& in)
00052   {
00053   arma_extra_debug_sigprint();
00054   
00055   const unwrap<T1> tmp(in.m);
00056   
00057   typedef typename T1::elem_type eT;
00058   const Mat<eT>& X = tmp.M;
00059   
00060   op_inv::apply(out, X);
00061   }
00062 
00063 
00064 
00065 //! inverse of diagmat(mat)
00066 template<typename T1>
00067 inline
00068 void
00069 op_inv::apply(Mat<typename T1::elem_type>& out, const Op< Op<T1,op_diagmat>, op_inv>& in)
00070   {
00071   arma_extra_debug_sigprint();
00072   
00073   const unwrap<T1> X_tmp(in.m.m);
00074   
00075   typedef typename T1::elem_type eT;
00076   const Mat<eT>& X = X_tmp.M;
00077   
00078   arma_debug_check( !X.is_square(), "op_inv::apply(): matrix must be square" );
00079   
00080   if(&out != &X)
00081     {
00082     out.zeros(X.n_rows, X.n_rows);
00083     
00084     for(u32 i=0; i<X.n_rows; ++i)
00085       {
00086       out.at(i,i) = 1.0 / X.at(i,i);
00087       }
00088     }
00089   else
00090     {
00091     podarray<eT> tmp(X.n_rows);
00092     
00093     for(u32 i=0; i<X.n_rows; ++i)
00094       {
00095       tmp[i] = X.at(i,i);
00096       }
00097       
00098     out.zeros(X.n_rows, X.n_rows);
00099     
00100     for(u32 i=0; i<X.n_rows; ++i)
00101       {
00102       out.at(i,i) = eT(1) / tmp.mem[i];
00103       }
00104     
00105     }
00106   
00107   }
00108 
00109 
00110 
00111 template<typename eT>
00112 inline
00113 void
00114 op_inv::apply_diagvec(Mat<eT>& out, const Mat<eT>& X)
00115   {
00116   arma_extra_debug_sigprint();
00117 
00118   arma_debug_check( !X.is_vec(), "op_inv::apply_diagvec(): internal error: can't interpret as a vector");
00119   
00120   if(&out != &X)
00121     {
00122     out.zeros(X.n_elem, X.n_elem);
00123     
00124     for(u32 i=0; i<X.n_elem; ++i)
00125       {
00126       out.at(i,i) = eT(1) / X.mem[i];
00127       }
00128     }
00129   else
00130     {
00131     podarray<eT> tmp(X.n_elem);
00132     
00133     for(u32 i=0; i<X.n_elem; ++i)
00134       {
00135       tmp[i] = X.mem[i];
00136       }
00137       
00138     out.zeros(X.n_elem, X.n_elem);
00139     
00140     for(u32 i=0; i<X.n_elem; ++i)
00141       {
00142       out.at(i,i) = eT(1) / tmp.mem[i];
00143       }
00144     
00145     }
00146 
00147   }
00148 
00149 
00150 
00151 //! inverse of diagmat(colvec or rowvec)
00152 template<typename eT>
00153 inline
00154 void
00155 op_inv::apply(Mat<eT>& out, const Op< Op<Mat<eT>,op_diagmat_vec>, op_inv>& in)
00156   {
00157   arma_extra_debug_sigprint();
00158   
00159   const Mat<eT>& X = in.m.m;
00160   op_inv::apply_diagvec(out, X);
00161   }
00162 
00163 
00164 
00165 //! @}