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 typename T1::elem_type
00025 accu(const Base<typename T1::elem_type,T1>& X)
00026 {
00027 arma_extra_debug_sigprint();
00028
00029 typedef typename T1::elem_type eT;
00030
00031 const unwrap<T1> tmp(X.get_ref());
00032 const Mat<eT>& A = tmp.M;
00033
00034 const u32 A_n_elem = A.n_elem;
00035 const eT* A_mem = A.mem;
00036
00037 eT val = eT(0);
00038
00039 for(u32 i=0; i<A_n_elem; ++i)
00040 {
00041 val += A_mem[i];
00042 }
00043
00044 return val;
00045 }
00046
00047
00048
00049
00050 template<typename T1>
00051 inline
00052 typename T1::elem_type
00053 accu(const Op<T1, op_diagmat>& X)
00054 {
00055 arma_extra_debug_sigprint();
00056
00057 typedef typename T1::elem_type eT;
00058
00059 const unwrap<T1> tmp(X.m);
00060 const Mat<eT>& A = tmp.M;
00061
00062 arma_debug_check( !A.is_square(), "accu(): sum of diagonal values of a non-square matrix requested" );
00063
00064 eT acc = eT(0);
00065
00066 for(u32 i=0; i<A.n_rows; ++i)
00067 {
00068 acc += A.at(i,i);
00069 }
00070
00071 return acc;
00072 }
00073
00074
00075
00076 template<typename eT>
00077 inline
00078 eT
00079 accu(const Op<Mat<eT>, op_diagmat_vec>& X)
00080 {
00081 arma_extra_debug_sigprint();
00082
00083 const Mat<eT>& A = X.m;
00084 arma_debug_check( !A.is_vec(), "accu(): internal error: expected a vector" );
00085
00086 return accu(A);
00087 }
00088
00089
00090
00091
00092 template<typename T1>
00093 inline
00094 typename T1::elem_type
00095 accu(const Op<T1, op_square>& in)
00096 {
00097 arma_extra_debug_sigprint();
00098
00099 typedef typename T1::elem_type eT;
00100
00101 const unwrap<T1> tmp(in.m);
00102 const Mat<eT>& A = tmp.M;
00103
00104 const u32 A_n_elem = A.n_elem;
00105 const eT* A_mem = A.mem;
00106
00107 eT acc = eT(0);
00108
00109 for(u32 i=0; i<A_n_elem; ++i)
00110 {
00111 const eT val = A_mem[i];
00112 acc += val*val;
00113 }
00114
00115 return acc;
00116 }
00117
00118
00119
00120
00121 template<typename T1>
00122 inline
00123 typename T1::elem_type
00124 accu(const Op<T1, op_sqrt>& in)
00125 {
00126 arma_extra_debug_sigprint();
00127
00128 typedef typename T1::elem_type eT;
00129
00130 const unwrap<T1> tmp(in.m);
00131 const Mat<eT>& A = tmp.M;
00132
00133 const u32 A_n_elem = A.n_elem;
00134 const eT* A_mem = A.mem;
00135
00136 eT acc = eT(0);
00137 for(u32 i=0; i<A_n_elem; ++i)
00138 {
00139 acc += std::sqrt(A_mem[i]);
00140 }
00141
00142 return acc;
00143 }
00144
00145
00146
00147
00148 template<typename T1, typename T2>
00149 inline
00150 typename T1::elem_type
00151 accu(const Op< Glue<T1,T2, glue_minus>, op_square>& in)
00152 {
00153 arma_extra_debug_sigprint();
00154
00155 typedef typename T1::elem_type eT;
00156
00157 const unwrap<T1> tmp1(in.m.A);
00158 const unwrap<T2> tmp2(in.m.B);
00159
00160 const Mat<eT>& A = tmp1.M;
00161 const Mat<eT>& B = tmp2.M;
00162
00163 arma_debug_assert_same_size(A,B, "accu()");
00164
00165 const u32 n_elem = A.n_elem;
00166
00167 eT acc = eT(0);
00168 for(u32 i=0; i<n_elem; ++i)
00169 {
00170 const eT val = A.mem[i] - B.mem[i];
00171 acc += val*val;
00172 }
00173
00174 return acc;
00175 }
00176
00177
00178
00179
00180 template<typename eT>
00181 inline
00182 eT
00183 accu(const subview<eT>& X)
00184 {
00185 arma_extra_debug_sigprint();
00186
00187 eT val = eT(0);
00188 for(u32 col=0; col<X.n_cols; ++col)
00189 {
00190 const eT* coldata = X.colptr(col);
00191
00192 for(u32 row=0; row<X.n_rows; ++row)
00193 {
00194 val += coldata[row];
00195 }
00196
00197 }
00198
00199 return val;
00200 }
00201
00202
00203
00204
00205 template<typename eT>
00206 inline
00207 eT
00208 accu(const diagview<eT>& X)
00209 {
00210 arma_extra_debug_sigprint();
00211
00212 const u32 n_elem = X.n_elem;
00213 eT val = eT(0);
00214
00215 for(u32 i=0; i<n_elem; ++i)
00216 {
00217 val += X[i];
00218 }
00219
00220 return val;
00221 }
00222
00223
00224
00225
00226 template<typename eT>
00227 inline
00228 eT
00229 accu_schur(const Mat<eT>& A, const Mat<eT>& B)
00230 {
00231 arma_extra_debug_sigprint();
00232
00233 arma_debug_assert_same_size(A,B, "accu()");
00234
00235 const eT* const A_mem = A.mem;
00236 const eT* const B_mem = B.mem;
00237
00238 const u32 n_elem = A.n_elem;
00239 eT val = eT(0);
00240
00241 for(u32 i=0; i<n_elem; ++i)
00242 {
00243 val += A_mem[i] * B_mem[i];
00244 }
00245
00246 return val;
00247 }
00248
00249
00250
00251
00252 template<typename eT>
00253 inline
00254 eT
00255 accu(const Glue<Mat<eT>,Mat<eT>,glue_schur>& X)
00256 {
00257 return accu_schur(X.A, X.B);
00258 }
00259
00260
00261
00262
00263 template<typename eT>
00264 inline
00265 eT
00266 accu(const Glue<Glue<Mat<eT>,Mat<eT>,glue_schur>,Mat<eT>,glue_schur>& X)
00267 {
00268 arma_extra_debug_sigprint();
00269
00270 const Mat<eT>& A = X.A.A;
00271 const Mat<eT>& B = X.A.B;
00272 const Mat<eT>& C = X.B;
00273
00274 arma_debug_assert_same_size(A,B, "accu()");
00275 arma_debug_assert_same_size(A,C, "accu()");
00276
00277 const eT* const A_mem = A.mem;
00278 const eT* const B_mem = B.mem;
00279 const eT* const C_mem = C.mem;
00280
00281 const u32 n_elem = A.n_elem;
00282 eT val = eT(0);
00283
00284 for(u32 i=0; i<n_elem; ++i)
00285 {
00286 val += A_mem[i] * B_mem[i] * C_mem[i];
00287 }
00288
00289 return val;
00290 }
00291
00292
00293
00294
00295
00296
00297
00298
00299 template<typename T1, typename T2>
00300 inline
00301 typename T1::elem_type
00302 accu(const Glue<T1,T2,glue_schur>& X)
00303 {
00304 arma_extra_debug_sigprint();
00305
00306 isnt_same_type<typename T1::elem_type, typename T2::elem_type>::check();
00307
00308 typedef typename T1::elem_type eT;
00309
00310 const u32 N_mat = 1 + depth_lhs< glue_schur, Glue<T1,T2,glue_schur> >::num;
00311 arma_extra_debug_print( arma_boost::format("N_mat = %d") % N_mat );
00312
00313 if(N_mat == 2)
00314 {
00315 return accu_schur(Mat<eT>(X.A), Mat<eT>(X.B));
00316 }
00317 else
00318 {
00319 const Mat<eT>* ptrs[N_mat];
00320 bool del[N_mat];
00321
00322 mat_ptrs<glue_schur, Glue<T1,T2,glue_schur> >::get_ptrs(ptrs, del, X);
00323
00324 for(u32 i=0; i<N_mat; ++i) arma_extra_debug_print( arma_boost::format("ptrs[%d] = %x") % i % ptrs[i] );
00325 for(u32 i=0; i<N_mat; ++i) arma_extra_debug_print( arma_boost::format(" del[%d] = %d") % i % del[i] );
00326
00327 const Mat<eT>& tmp_mat = *(ptrs[0]);
00328
00329 for(u32 i=1; i<N_mat; ++i)
00330 {
00331 arma_debug_assert_same_size(tmp_mat, *(ptrs[i]), "accu()");
00332 }
00333
00334
00335
00336
00337 eT val = eT(0);
00338
00339 const u32 n_elem = ptrs[0]->n_elem;
00340
00341 for(u32 j=0; j<n_elem; ++j)
00342 {
00343 eT tmp = ptrs[0]->mem[j];
00344
00345 for(u32 i=1; i<N_mat; ++i)
00346 {
00347 tmp *= ptrs[i]->mem[j];
00348 }
00349
00350 val += tmp;
00351 }
00352
00353
00354 for(u32 i=0; i<N_mat; ++i)
00355 {
00356 if(del[i] == true)
00357 {
00358 arma_extra_debug_print( arma_boost::format("delete mat_ptr[%d]") % i );
00359 delete ptrs[i];
00360 }
00361 }
00362
00363 return val;
00364 }
00365 }
00366
00367
00368
00369 template<typename eT>
00370 inline
00371 eT
00372 accu(const Glue<subview<eT>,Mat<eT>,glue_schur>& X)
00373 {
00374 arma_extra_debug_sigprint();
00375
00376 arma_debug_assert_same_size(X.A, X.B, "accu()");
00377
00378 const Mat<eT>& A = X.A.m;
00379 const Mat<eT>& B = X.B;
00380
00381 const u32 A_sub_n_rows = X.A.n_rows;
00382 const u32 A_sub_n_cols = X.A.n_cols;
00383
00384 const u32 A_aux_row1 = X.A.aux_row1;
00385 const u32 A_aux_col1 = X.A.aux_col1;
00386
00387
00388 eT val = eT(0);
00389
00390 for(u32 col = 0; col<A_sub_n_cols; ++col)
00391 {
00392 const u32 col_mod = A_aux_col1 + col;
00393
00394 for(u32 row = 0; row<A_sub_n_rows; ++row)
00395 {
00396 const u32 row_mod = A_aux_row1 + row;
00397
00398 val += A.at(row_mod, col_mod) * B.at(row,col);
00399 }
00400
00401 }
00402
00403 return val;
00404 }
00405
00406
00407
00408
00409 template<typename eT>
00410 inline
00411 eT
00412 accu(const Glue<Mat<eT>,subview<eT>,glue_schur>& X)
00413 {
00414 arma_extra_debug_sigprint();
00415
00416 arma_debug_assert_same_size(X.A, X.B, "accu()");
00417
00418 const Mat<eT>& A = X.A;
00419 const Mat<eT>& B = X.B.m;
00420
00421
00422
00423
00424 const u32 B_aux_row1 = X.B.aux_row1;
00425 const u32 B_aux_col1 = X.B.aux_col1;
00426
00427
00428 eT val = eT(0);
00429
00430 for(u32 col = 0; col<A.n_cols; ++col)
00431 {
00432 const u32 col_mod = B_aux_col1 + col;
00433
00434 for(u32 row = 0; row<A.n_rows; ++row)
00435 {
00436 const u32 row_mod = B_aux_row1 + row;
00437
00438 val += A.at(row, col) * B.at(row_mod, col_mod);
00439 }
00440
00441 }
00442
00443 return val;
00444 }
00445
00446
00447
00448
00449 template<typename eT>
00450 inline
00451 eT
00452 accu(const Glue<subview<eT>,subview<eT>,glue_schur>& X)
00453 {
00454 arma_extra_debug_sigprint();
00455
00456 arma_debug_assert_same_size(X.A, X.B, "accu()");
00457
00458 const Mat<eT>& A = X.A.m;
00459 const Mat<eT>& B = X.B.m;
00460
00461 const u32 A_sub_n_rows = X.A.n_rows;
00462 const u32 A_sub_n_cols = X.A.n_cols;
00463
00464
00465
00466
00467 const u32 A_aux_row1 = X.A.aux_row1;
00468 const u32 A_aux_col1 = X.A.aux_col1;
00469
00470 const u32 B_aux_row1 = X.B.aux_row1;
00471 const u32 B_aux_col1 = X.B.aux_col1;
00472
00473
00474 eT val = eT(0);
00475
00476 for(u32 col = 0; col<A_sub_n_cols; ++col)
00477 {
00478 const u32 A_col_mod = A_aux_col1 + col;
00479 const u32 B_col_mod = B_aux_col1 + col;
00480
00481 for(u32 row = 0; row<A_sub_n_rows; ++row)
00482 {
00483 const u32 A_row_mod = A_aux_row1 + row;
00484 const u32 B_row_mod = B_aux_row1 + row;
00485
00486 val += A.at(A_row_mod, A_col_mod) * B.at(B_row_mod, B_col_mod);
00487 }
00488
00489 }
00490
00491 return val;
00492 }
00493
00494
00495
00496