@@ -162,6 +162,17 @@ glue_cubemul::apply_noalias(Cube<eT>& out, const Cube<eT>& A, const Cube<eT>& B)
162162 {
163163 arma_debug_sigprint ();
164164
165+ constexpr bool has_xtrans = (do_strans_A || do_htrans_A || do_strans_B || do_htrans_B);
166+
167+ if ( has_xtrans && (&A == &B) )
168+ {
169+ arma_debug_print (" optimisation for symmetric/hermitian result" );
170+
171+ glue_cubemul::apply_noalias_single<eT, do_strans_A, do_htrans_A, do_strans_B, do_htrans_B>(out, A);
172+
173+ return ;
174+ }
175+
165176 const uword A_n_slices = A.n_slices ;
166177
167178 arma_conform_check ( (A_n_slices != B.n_slices ), " cubemul(): given cubes must have the same number of slices" );
@@ -226,6 +237,76 @@ glue_cubemul::apply_noalias(Cube<eT>& out, const Cube<eT>& A, const Cube<eT>& B)
226237
227238
228239
240+ template <typename eT, bool do_strans_A, bool do_htrans_A, bool do_strans_B, bool do_htrans_B>
241+ inline
242+ void
243+ glue_cubemul::apply_noalias_single (Cube<eT>& out, const Cube<eT>& A)
244+ {
245+ arma_debug_sigprint ();
246+
247+ const uword A_n_slices = A.n_slices ;
248+
249+ constexpr bool do_xtrans_A = (do_strans_A || do_htrans_A);
250+ constexpr bool do_xtrans_B = (do_strans_B || do_htrans_B);
251+
252+ const uword final_A_n_rows = (do_xtrans_A == false ) ? A.n_rows : A.n_cols ;
253+ const uword final_A_n_cols = (do_xtrans_A == false ) ? A.n_cols : A.n_rows ;
254+
255+ const uword final_B_n_rows = (do_xtrans_B == false ) ? A.n_rows : A.n_cols ;
256+ const uword final_B_n_cols = (do_xtrans_B == false ) ? A.n_cols : A.n_rows ;
257+
258+ if ( (arma_config::check_conform) && (final_A_n_cols != final_B_n_rows) )
259+ {
260+ std::ostringstream tmp;
261+
262+ tmp << " cubemul(): incompatible dimensions: " << final_A_n_rows << ' x' << final_A_n_cols << ' x' << A_n_slices << " and " << final_B_n_rows << ' x' << final_B_n_cols << ' x' << A_n_slices;
263+
264+ arma_stop_logic_error (tmp.str ());
265+ }
266+
267+ out.set_size (final_A_n_rows, final_B_n_cols, A_n_slices);
268+
269+ if (out.is_empty ()) { return ; }
270+
271+ if (A.is_empty ()) { out.zeros (); return ; }
272+
273+ for (uword s=0 ; s < A_n_slices; ++s)
274+ {
275+ const Mat<eT> A_slice_s (const_cast <eT*>(A.slice_memptr (s)), A.n_rows , A.n_cols , false , true );
276+
277+ Mat<eT> out_slice_s (out.slice_memptr (s), final_A_n_rows, final_B_n_cols, false , true );
278+
279+ // glue_times has optimised handling for A*A' and A'*A
280+
281+ if ( (do_xtrans_A == false ) && (do_xtrans_B == false ) )
282+ {
283+ out_slice_s = A_slice_s * A_slice_s;
284+ }
285+ else
286+ if ( (do_xtrans_A == false ) && (do_xtrans_B == true ) )
287+ {
288+ if (do_strans_B) { out_slice_s = A_slice_s * strans (A_slice_s); }
289+ else if (do_htrans_B) { out_slice_s = A_slice_s * htrans (A_slice_s); }
290+ }
291+ else
292+ if ( (do_xtrans_A == true ) && (do_xtrans_B == false ) )
293+ {
294+ if (do_strans_A) { out_slice_s = strans (A_slice_s) * A_slice_s; }
295+ else if (do_htrans_A) { out_slice_s = htrans (A_slice_s) * A_slice_s; }
296+ }
297+ else
298+ if ( (do_xtrans_A == true ) && (do_xtrans_B == true ) )
299+ {
300+ if ( (do_strans_A) && (do_strans_B) ) { out_slice_s = strans (A_slice_s) * strans (A_slice_s); }
301+ else if ( (do_strans_A) && (do_htrans_B) ) { out_slice_s = strans (A_slice_s) * htrans (A_slice_s); }
302+ else if ( (do_htrans_A) && (do_strans_B) ) { out_slice_s = htrans (A_slice_s) * strans (A_slice_s); }
303+ else if ( (do_htrans_A) && (do_htrans_B) ) { out_slice_s = htrans (A_slice_s) * htrans (A_slice_s); }
304+ }
305+ }
306+ }
307+
308+
309+
229310//
230311
231312
0 commit comments