op_diagmat_meat.hpp

Go to the documentation of this file.
00001 // Copyright (C) 2010 NICTA and the authors listed below
00002 // http://nicta.com.au
00003 // 
00004 // Authors:
00005 // - Conrad Sanderson (conradsand at ieee dot org)
00006 // 
00007 // This file is part of the Armadillo C++ library.
00008 // It is provided without any warranty of fitness
00009 // for any purpose. You can redistribute this file
00010 // and/or modify it under the terms of the GNU
00011 // Lesser General Public License (LGPL) as published
00012 // by the Free Software Foundation, either version 3
00013 // of the License or (at your option) any later version.
00014 // (see http://www.opensource.org/licenses for more info)
00015 
00016 
00017 //! \addtogroup op_diagmat
00018 //! @{
00019 
00020 
00021 
00022 template<typename T1>
00023 inline
00024 void
00025 op_diagmat::apply(Mat<typename T1::elem_type>& out, const Op<T1, op_diagmat>& X)
00026   {
00027   arma_extra_debug_sigprint();
00028   
00029   typedef typename T1::elem_type eT;
00030   
00031   const unwrap<T1>   tmp(X.m);
00032   const Mat<eT>& A = tmp.M;
00033   
00034   if(A.is_vec() == true)
00035     {
00036     // generate a diagonal matrix out of a vector
00037     
00038     const u32 N     = A.n_elem;
00039     const eT* A_mem = A.memptr();
00040     
00041     if(&out != &A)
00042       {
00043       // no aliasing
00044       out.zeros(N,N);
00045       
00046       for(u32 i=0; i<N; ++i)
00047         {
00048         out.at(i,i) = A_mem[i];
00049         }
00050       }
00051     else
00052       {
00053       // aliasing
00054       
00055       const podarray<eT> tmp(A_mem, N);
00056       
00057       const eT* tmp_mem = tmp.memptr();
00058       
00059       out.zeros(N,N);
00060       
00061       for(u32 i=0; i<N; ++i)
00062         {
00063         out.at(i,i) = tmp_mem[i];
00064         }
00065       }
00066     }
00067   else
00068     {
00069     // generate a diagonal matrix out of a matrix
00070     
00071     arma_debug_check( (A.is_square() == false), "diagmat(): given matrix is not square" );
00072     
00073     const u32 N = A.n_rows;
00074     
00075     out.set_size(N,N);
00076     
00077     for(u32 col=0; col<N; ++col)
00078       {
00079       for(u32 row=0;     row<col; ++row) { out.at(row,col) = eT(0); }
00080       
00081       out.at(col,col) = A.at(col,col);
00082       
00083       for(u32 row=col+1; row<N;   ++row) { out.at(row,col) = eT(0); }
00084       }
00085     }
00086   }
00087 
00088 
00089 
00090 //! @}