op_diagmat_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_diagmat
00017 //! @{
00018 
00019 
00020 
00021 template<typename eT>
00022 inline
00023 void
00024 op_diagmat::zero_offdiag(Mat<eT>& X)
00025   {
00026   for(u32 col=0; col<X.n_cols; ++col)
00027     {
00028     eT* colmem = X.colptr(col);
00029     
00030     // above the diagonal
00031     
00032     for(u32 row=0; row<col; ++row)
00033       {
00034       colmem[row] = eT(0);
00035       }
00036     
00037     // below the diagonal
00038     
00039     for(u32 row=col+1; row<X.n_rows; ++row)
00040       {
00041       colmem[row] = eT(0);
00042       }
00043     }
00044   }
00045 
00046 
00047 
00048 template<typename T1>
00049 inline
00050 void
00051 op_diagmat::apply(Mat<typename T1::elem_type>& out, const Op<T1,op_diagmat>& 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   arma_debug_check( !X.is_square(), "diagmat(): matrix must be square" );
00061   
00062   if(&out != &X)
00063     {
00064     out.zeros(X.n_rows, X.n_rows);
00065     
00066     for(u32 i=0; i<X.n_rows; ++i)
00067       {
00068       out.at(i,i) = X.at(i,i);
00069       }
00070     }
00071   else
00072     {
00073     op_diagmat::zero_offdiag(out);
00074     }
00075   }
00076 
00077 
00078 
00079 template<typename T1, typename T2>
00080 inline
00081 void
00082 op_diagmat::apply(Mat<typename T1::elem_type>& out, const Op< Glue<T1,T2, glue_div>, op_diagmat>& in)
00083   {
00084   arma_extra_debug_sigprint();
00085   
00086   const unwrap<T1> tmp1(in.m.A);
00087   const unwrap<T2> tmp2(in.m.B);
00088   
00089   typedef typename T1::elem_type eT;
00090   
00091   const Mat<eT>& A = tmp1.M;
00092   const Mat<eT>& B = tmp2.M;
00093   
00094   arma_debug_check( (!A.is_square() || !B.is_square()), "diagmat(): matrices must be square" );
00095   arma_debug_assert_same_size(A, B, "matrix element-wise division");
00096     
00097   // not using out.zeros() as 'out' might be an alias of A and/or B.
00098   // the off-diagonal elements are zeroed at the end
00099   out.set_size(A.n_rows, A.n_rows);
00100   
00101   for(u32 i=0; i<A.n_rows; ++i)
00102     {
00103     out.at(i,i) = A.at(i,i) / B.at(i,i);
00104     }
00105   
00106   op_diagmat::zero_offdiag(out);
00107   }
00108 
00109 
00110 
00111 template<typename T1, typename T2>
00112 inline
00113 void
00114 op_diagmat::apply(Mat<typename T1::elem_type>& out, const Op< Glue<T1,T2, glue_minus>, op_diagmat>& in)
00115   {
00116   arma_extra_debug_sigprint();
00117   
00118   const unwrap<T1> tmp1(in.m.A);
00119   const unwrap<T2> tmp2(in.m.B);
00120   
00121   typedef typename T1::elem_type eT;
00122   
00123   const Mat<eT>& A = tmp1.M;
00124   const Mat<eT>& B = tmp2.M;
00125   
00126   arma_debug_check( (!A.is_square() || !B.is_square()), "diagmat(): matrices must be square" );
00127   arma_debug_assert_same_size(A, B, "matrix subtraction");
00128     
00129   // not using out.zeros() as 'out' might be an alias of A and/or B.
00130   // the off-diagonal elements are zeroed at the end
00131   out.set_size(A.n_rows, A.n_rows);
00132   
00133   for(u32 i=0; i<A.n_rows; ++i)
00134     {
00135     out.at(i,i) = A.at(i,i) - B.at(i,i);
00136     }
00137   
00138   op_diagmat::zero_offdiag(out);
00139   }
00140 
00141 
00142 
00143 template<typename T1, typename T2>
00144 inline
00145 void
00146 op_diagmat::apply(Mat<typename T1::elem_type>& out, const Op< Glue<T1,T2, glue_plus>, op_diagmat>& in)
00147   {
00148   arma_extra_debug_sigprint();
00149   
00150   const unwrap<T1> tmp1(in.m.A);
00151   const unwrap<T2> tmp2(in.m.B);
00152   
00153   typedef typename T1::elem_type eT;
00154   
00155   const Mat<eT>& A = tmp1.M;
00156   const Mat<eT>& B = tmp2.M;
00157   
00158   arma_debug_check( (!A.is_square() || !B.is_square()), "diagmat(): matrices must be square" );
00159   arma_debug_assert_same_size(A, B, "matrix addition");
00160     
00161   // not using out.zeros() as 'out' might be an alias of A and/or B.
00162   // the off-diagonal elements are zeroed at the end
00163   out.set_size(A.n_rows, A.n_rows);
00164   
00165   for(u32 i=0; i<A.n_rows; ++i)
00166     {
00167     out.at(i,i) = A.at(i,i) + B.at(i,i);
00168     }
00169   
00170   op_diagmat::zero_offdiag(out);
00171   }
00172 
00173 
00174 
00175 template<typename T1, typename T2>
00176 inline
00177 void
00178 op_diagmat::apply(Mat<typename T1::elem_type>& out, const Op< Glue<T1,T2, glue_schur>, op_diagmat>& in)
00179   {
00180   arma_extra_debug_sigprint();
00181   
00182   const unwrap<T1> tmp1(in.m.A);
00183   const unwrap<T2> tmp2(in.m.B);
00184   
00185   typedef typename T1::elem_type eT;
00186   
00187   const Mat<eT>& A = tmp1.M;
00188   const Mat<eT>& B = tmp2.M;
00189   
00190   arma_debug_check( (!A.is_square() || !B.is_square()), "diagmat(): matrices must be square" );
00191   arma_debug_assert_same_size(A, B, "matrix schur product");
00192     
00193   // not using out.zeros() as 'out' might be an alias of A and/or B.
00194   // the off-diagonal elements are zeroed at the end
00195   out.set_size(A.n_rows, A.n_rows);
00196   
00197   for(u32 i=0; i<A.n_rows; ++i)
00198     {
00199     out.at(i,i) = A.at(i,i) * B.at(i,i);
00200     }
00201   
00202   op_diagmat::zero_offdiag(out);
00203   }
00204 
00205 
00206 
00207 template<typename T1, typename T2>
00208 inline
00209 void
00210 op_diagmat::apply(Mat<typename T1::elem_type>& out, const Op< Glue<T1,T2, glue_times>, op_diagmat>& in)
00211   {
00212   arma_extra_debug_sigprint();
00213   
00214   const unwrap_check<T1> tmp1(in.m.A, out);
00215   const unwrap_check<T2> tmp2(in.m.B, out);
00216   
00217   typedef typename T1::elem_type eT;
00218   
00219   const Mat<eT>& A = tmp1.M;
00220   const Mat<eT>& B = tmp2.M;
00221   
00222   arma_debug_check( (A.n_rows != B.n_cols), "diagmat(): result of multiplication is not square" );
00223   arma_debug_assert_mul_size(A, B, "matrix multiplication");
00224     
00225   // out is cleared here, as we've made sure that A and B are not aliases of 'out'
00226   out.zeros(A.n_rows, A.n_rows);
00227   
00228   for(u32 i=0; i<A.n_rows; ++i)
00229     {
00230     const eT* B_colmem = B.colptr(i);
00231     
00232     eT val = eT(0);
00233     for(u32 j=0; j<A.n_cols; ++j)
00234       {
00235       val += A.at(i,j) * B_colmem[j];
00236       }
00237     
00238     out.at(i,i) = val;
00239     }
00240   }
00241 
00242 
00243 
00244 //
00245 //
00246 //
00247 
00248 
00249 template<typename T1>
00250 inline
00251 void
00252 op_diagmat_vec::apply(Mat<typename T1::elem_type>& out, const Op<T1, op_diagmat_vec>& in)
00253   {
00254   arma_extra_debug_sigprint();
00255   
00256   const unwrap<T1> tmp(in.m);
00257   
00258   typedef typename T1::elem_type eT;
00259   const Mat<eT>& X = tmp.M;
00260   
00261   if(&out != &X)
00262     {
00263     out.zeros(X.n_elem, X.n_elem);
00264     
00265     const eT* X_mem = X.mem;
00266     
00267     for(u32 i=0; i<X.n_elem; ++i)
00268       {
00269       out.at(i,i) = X_mem[i];
00270       }
00271     }
00272   else
00273     {
00274     podarray<eT> tmp_array(X.n_elem);
00275     
00276     for(u32 i=0; i<X.n_elem; ++i)
00277       {
00278       tmp_array[i] = X[i];
00279       }
00280     
00281     out.zeros(tmp_array.n_elem, tmp_array.n_elem);
00282     
00283     for(u32 i=0; i<tmp_array.n_elem; ++i)
00284       {
00285       out.at(i,i) = tmp_array[i];
00286       }
00287     
00288     }
00289   }
00290 
00291 
00292 
00293 //! @}