Go to the documentation of this file.00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
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
00037
00038 const u32 N = A.n_elem;
00039 const eT* A_mem = A.memptr();
00040
00041 if(&out != &A)
00042 {
00043
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
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
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