00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
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
00031
00032 for(u32 row=0; row<col; ++row)
00033 {
00034 colmem[row] = eT(0);
00035 }
00036
00037
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
00098
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
00130
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
00162
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
00194
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
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