AbstractLinAlgPack: C++ Interfaces For Vectors, Matrices And Related Linear Algebra Objects  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
AbstractLinAlgPack_MatrixOpSerial.cpp
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
5 // Copyright (2003) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov)
38 //
39 // ***********************************************************************
40 // @HEADER
41 
42 // ToDo: 7/27/99: Give these default implementations and test them.
43 
44 #include <typeinfo>
45 
46 #include "AbstractLinAlgPack_MatrixOpSerial.hpp"
47 #include "AbstractLinAlgPack_VectorDenseEncap.hpp"
48 #include "AbstractLinAlgPack_MatrixOpGetGMSMutable.hpp"
49 #include "AbstractLinAlgPack_MatrixOpGetGMSTri.hpp"
50 #include "AbstractLinAlgPack_MatrixSymOpGetGMSSymMutable.hpp"
51 #include "AbstractLinAlgPack_SpVectorOp.hpp"
52 #include "AbstractLinAlgPack_GenPermMatrixSliceOp.hpp"
53 #include "AbstractLinAlgPack_EtaVector.hpp"
54 #include "AbstractLinAlgPack_GenPermMatrixSlice.hpp"
55 #include "AbstractLinAlgPack_LinAlgOpPackHack.hpp"
56 #include "DenseLinAlgPack_DMatrixClass.hpp"
57 #include "DenseLinAlgPack_DMatrixOut.hpp"
58 #include "DenseLinAlgPack_AssertOp.hpp"
59 #include "Teuchos_Workspace.hpp"
60 #include "Teuchos_dyn_cast.hpp"
61 
62 namespace LinAlgOpPack {
66 }
67 
68 namespace AbstractLinAlgPack {
69 
70 // Level-1 BLAS
71 
72 void MatrixOpSerial::Mp_StM(DMatrixSlice* gms_lhs, value_type alpha
73  , BLAS_Cpp::Transp trans_rhs) const
74 {
75  DenseLinAlgPack::Mp_M_assert_sizes( gms_lhs->rows(), gms_lhs->cols(), BLAS_Cpp::no_trans
76  , rows(), cols(), trans_rhs );
77  const size_type
78  m = gms_lhs->rows(),
79  n = gms_lhs->cols();
80  //
81  // Use sparse matrix-vector multiplication to perform this operation.
82  // C += a * B = a * B * I = [ a*B*e(1), a*B*e(2), ..., a*B*e(m) ]
83  //
84  SpVector rhs;
85  rhs.uninitialized_resize( n, 1, 1 );
86  for( size_type j = 1; j <=n; ++j ) {
87  rhs.begin()->initialize( j, 1.0 ); // e(j)
88  this->Vp_StMtV( &gms_lhs->col(j), alpha, trans_rhs, rhs(), 1.0 );
89  }
90 }
91 
93  , BLAS_Cpp::Transp M_trans
94  , const GenPermMatrixSlice& P, BLAS_Cpp::Transp P_trans
95  ) const
96 {
97  // C += a * op(M) * op(P)
98  TEUCHOS_TEST_FOR_EXCEPT(true); // Implement this!
99 }
100 
102  , const GenPermMatrixSlice& P, BLAS_Cpp::Transp P_trans
103  , BLAS_Cpp::Transp M_trans
104  ) const
105 {
106  // C += a * op(P) * op(M)
107  TEUCHOS_TEST_FOR_EXCEPT(true); // Implement this!
108 }
109 
111  , const GenPermMatrixSlice& P1, BLAS_Cpp::Transp P1_trans
112  , BLAS_Cpp::Transp M_trans
113  , const GenPermMatrixSlice& P2, BLAS_Cpp::Transp P2_trans
114  ) const
115 {
116  // C += a * op(P1) * op(M) * op(P2)
117  TEUCHOS_TEST_FOR_EXCEPT(true); // Implement this!
118 }
119 
120 // Level-2 BLAS
121 
122 void MatrixOpSerial::Vp_StMtV(DVectorSlice* vs_lhs, value_type alpha
123  , BLAS_Cpp::Transp trans_rhs1, const SpVectorSlice& sv_rhs2, value_type beta) const
124 {
125  Vp_MtV_assert_sizes( vs_lhs->dim(), rows(), cols(), trans_rhs1, sv_rhs2.dim() );
126  if( !sv_rhs2.nz() ) {
127  // vs_lhs = beta * vs_lhs
128  if(beta==0.0) *vs_lhs = 0.0;
129  else if(beta!=1.0) DenseLinAlgPack::Vt_S(vs_lhs,beta);
130  }
131  else {
132  // Convert to dense by default.
133  if( sv_rhs2.dim() == sv_rhs2.nz() && sv_rhs2.is_sorted() ) {
134  const DVectorSlice vs_rhs2 = AbstractLinAlgPack::dense_view(sv_rhs2);
135  this->Vp_StMtV( vs_lhs, alpha, trans_rhs1, vs_rhs2, beta );
136  }
137  else {
138  DVector v_rhs2;
139  LinAlgOpPack::assign( &v_rhs2, sv_rhs2 );
140  this->Vp_StMtV( vs_lhs, alpha, trans_rhs1, v_rhs2(), beta );
141  }
142  }
143 }
144 
146  , const GenPermMatrixSlice& P, BLAS_Cpp::Transp P_trans
147  , BLAS_Cpp::Transp M_trans
148  , const DVectorSlice& x, value_type b) const
149 {
150  using Teuchos::Workspace;
152 
153  Workspace<value_type> t_ws(wss,BLAS_Cpp::cols(P.rows(),P.cols(),P_trans));
154  DVectorSlice t(&t_ws[0],t_ws.size());
155  LinAlgOpPack::V_StMtV(&t,a,*this,M_trans,x);
156  LinAlgOpPack::Vp_MtV( y, P, P_trans, t, b );
157 }
158 
160  , const GenPermMatrixSlice& P, BLAS_Cpp::Transp P_trans
161  , BLAS_Cpp::Transp M_trans
162  , const SpVectorSlice& x, value_type b) const
163 {
164  using Teuchos::Workspace;
166 
167  Workspace<value_type> t_ws(wss,BLAS_Cpp::cols(P.rows(),P.cols(),P_trans));
168  DVectorSlice t(&t_ws[0],t_ws.size());
169  LinAlgOpPack::V_StMtV(&t,a,*this,M_trans,x);
170  LinAlgOpPack::Vp_MtV( y, P, P_trans, t, b );
171 }
172 
174  , BLAS_Cpp::Transp M_trans, const DVectorSlice& x2) const
175 {
176  DenseLinAlgPack::Vp_MtV_assert_sizes( x1.dim(), rows(), cols(), M_trans, x2.dim() );
177  DVector tmp(x1.dim());
178  this->Vp_StMtV( &tmp(), 1.0, M_trans, x2, 0.0 );
179  return DenseLinAlgPack::dot( x1, tmp() );
180 }
181 
183  , BLAS_Cpp::Transp M_trans, const SpVectorSlice& x2) const
184 {
185  DenseLinAlgPack::Vp_MtV_assert_sizes( x1.dim(), rows(), cols(), M_trans, x2.dim() );
186  if( !x1.nz() || !x2.nz() ) {
187  return 0.0;
188  }
189  else {
190  if( x1.overlap(x2) == DenseLinAlgPack::SAME_MEM && x1.dim() == x1.nz() && x1.is_sorted() ) {
191  const DVectorSlice x1_d = AbstractLinAlgPack::dense_view(x1);
192  return this->transVtMtV( x1_d, M_trans, x1_d );
193  }
194  DVector tmp(x1.dim());
195  this->Vp_StMtV( &tmp(), 1.0, M_trans, x2, 0.0 );
196  return AbstractLinAlgPack::dot( x1, tmp() );
197  }
198 }
199 
201  BLAS_Cpp::Transp M_trans_in, value_type a
202  , const GenPermMatrixSlice& P1_in, BLAS_Cpp::Transp P1_trans
203  , const GenPermMatrixSlice& P2_in, BLAS_Cpp::Transp P2_trans
204  , value_type b, DMatrixSliceSym* S ) const
205 {
206  using BLAS_Cpp::no_trans;
207  using BLAS_Cpp::trans;
208  using BLAS_Cpp::trans_not;
209  using BLAS_Cpp::rows;
210  using BLAS_Cpp::cols;
211  //
212  // S = b * S
213  //
214  // S += a*op(P1')*op(M)*op(P2) + a*op(P2')*op(M')*op(P1)
215  //
216  // We will start by renaming P1 and P2 such that op(P1).rows() >= op(P2).rows().
217  // This is because we are going to store some temparary vectors and we don't
218  // want them to be too big.
219  //
220  // We will perform the above operation by working with columns of:
221  //
222  // op(P1)(:,j(k)) = e(i(k)) <: R^n
223  //
224  // Then for each column in op(P1) we will perform:
225  //
226  //
227  // for k = 1...P1.nz()
228  //
229  // [ . ]
230  // S += a * [ e(i(k))' ] * op(M)*op(P2) + a * op(P2') * op(M') * [ ... e(i(k)) ... ]
231  // [ . ]
232  // row j(k) col j(k)
233  // =>
234  // [ . ]
235  // S += a * [ y_k' ] + a * [ ... y_k ... ]
236  // [ . ]
237  // row j(k) col j(k)
238  //
239  // where: y_k = a * op(P2') * op(M') * e(i(k)) <: R^m
240  //
241  // Of course above we only need to set the row and column elements for S(j(k),:) and S(:,j(k))
242  // for the portion of the symmetric S that is being stored.
243  //
244  const size_type
245  M_rows = this->rows(),
246  M_cols = this->cols(),
247  P1_cols = cols( P1_in.rows(), P1_in.cols(), P1_trans );
248  DenseLinAlgPack::MtM_assert_sizes(
249  M_rows, M_cols, trans_not(M_trans_in)
250  , P1_in.rows(), P1_in.cols(), P1_trans );
251  DenseLinAlgPack::MtM_assert_sizes(
252  M_rows, M_cols, M_trans_in
253  , P2_in.rows(), P2_in.cols(), P2_trans );
254  DenseLinAlgPack::Mp_M_assert_sizes(
255  S->rows(), S->cols(), no_trans
256  , P1_cols, P1_cols, no_trans );
257  // Rename P1 and P2 so that op(P1).rows() >= op(P2).rows()
258  const bool
259  reorder_P1_P2 = ( rows( P1_in.rows(), P1_in.cols(), P1_trans )
260  < rows( P2_in.rows(), P2_in.cols(), P2_trans ) );
261  const GenPermMatrixSlice
262  &P1 = reorder_P1_P2 ? P2_in : P1_in,
263  &P2 = reorder_P1_P2 ? P1_in : P2_in;
264  const BLAS_Cpp::Transp
265  M_trans = reorder_P1_P2 ? trans_not(M_trans_in) : M_trans_in;
266  // Set rows and columns of S
267  const size_type
268  m = S->rows(),
269  n = rows( P1.rows(), P1.cols(), P1_trans );
270  DVector y_k_store(m); // ToDo: use workspace allocator!
271  DVectorSlice y_k = y_k_store();
272  for( GenPermMatrixSlice::const_iterator P1_itr = P1.begin(); P1_itr != P1.end(); ++P1_itr )
273  {
274  const size_type
275  i_k = P1_trans == no_trans ? P1_itr->row_i() : P1_itr->col_j(),
276  j_k = P1_trans == no_trans ? P1_itr->col_j() : P1_itr->row_i();
277  // e(i(k))
278  EtaVector
279  e_i_k(i_k,n);
280  // y_k = op(P2')*op(M')*e(i(k))
281  AbstractLinAlgPack::Vp_StPtMtV( &y_k, 1.0, P2, trans_not(P2_trans), *this, trans_not(M_trans), e_i_k(), 0.0 );
282  // S(j(k),:) += a*y_k'
283  if( S->uplo() == BLAS_Cpp::upper )
284  DenseLinAlgPack::Vp_StV( &S->gms().row(j_k)(1,j_k), a, y_k(1,j_k) );
285  else
286  DenseLinAlgPack::Vp_StV( &S->gms().row(j_k)(j_k,m), a, y_k(j_k,m) );
287  // S(:,j(k)) += a*y_k
288  if( S->uplo() == BLAS_Cpp::upper )
289  DenseLinAlgPack::Vp_StV( &S->gms().col(j_k)(1,j_k), a, y_k(1,j_k) );
290  else
291  DenseLinAlgPack::Vp_StV( &S->gms().col(j_k)(j_k,m), a, y_k(j_k,m) );
292  }
293 }
294 
295 // Level-3 BLAS
296 
298  , BLAS_Cpp::Transp A_trans, const DMatrixSlice& B
299  , BLAS_Cpp::Transp B_trans, value_type b) const
300 {
301  DenseLinAlgPack::Mp_MtM_assert_sizes( C->rows(), C->cols(), BLAS_Cpp::no_trans
302  , rows(), cols(), A_trans, B.rows(), B.cols(), B_trans );
303  //
304  // C = b*C + a*op(A)*op(B)
305  //
306  // C.col(j) = b*col(j) + a*op(A)*op(B).col(j)
307  //
308 
309  // ToDo: Add the ability to also perform by row if faster
310 
311  for( size_type j = 1; j <= C->cols(); ++j )
312  AbstractLinAlgPack::Vp_StMtV( &C->col(j), a, *this, A_trans, DenseLinAlgPack::col( B, B_trans, j ), b );
313 }
314 
315 void MatrixOpSerial::Mp_StMtM(DMatrixSlice* C, value_type a, const DMatrixSlice& A
316  , BLAS_Cpp::Transp A_trans, BLAS_Cpp::Transp B_trans, value_type b) const
317 {
318  DenseLinAlgPack::Mp_MtM_assert_sizes( C->rows(), C->cols(), BLAS_Cpp::no_trans
319  , A.rows(), A.cols(), A_trans, rows(), cols(), B_trans );
320  //
321  // C = b*C + a*op(A)*op(B)
322  //
323  // C.row(i) = b*row(i) + a*op(B)'*op(A).row(i)
324  //
325 
326  // ToDo: Add the ability to also perform by column if faster
327 
328  for( size_type i = 1; i <= C->rows(); ++i )
329  AbstractLinAlgPack::Vp_StMtV( &C->row(i), a, *this, BLAS_Cpp::trans_not(A_trans)
330  , DenseLinAlgPack::row(A,A_trans,i) , b );
331 }
332 
334  , BLAS_Cpp::Transp A_trans, const MatrixOpSerial& B
335  , BLAS_Cpp::Transp B_trans, value_type b) const
336 {
337  using LinAlgOpPack::assign;
338  // C = b*C + a*op(A)*op(B)
339  DenseLinAlgPack::Mp_MtM_assert_sizes( C->rows(), C->cols(), BLAS_Cpp::no_trans
340  , rows(), cols(), A_trans, B.rows(), B.cols(), B_trans );
341  // Convert one of the matrices to dense, which ever one is the smallest!
342  const size_type
343  size_A = rows() * cols(),
344  size_B = B.rows() * B.cols();
345  if( size_A < size_B ) {
346  DMatrix A_dense;
347  assign( &A_dense, *this, BLAS_Cpp::no_trans );
348  AbstractLinAlgPack::Mp_StMtM( C, a, A_dense(), A_trans, B, B_trans, b );
349  }
350  else {
351  DMatrix B_dense;
352  assign( &B_dense, B, BLAS_Cpp::no_trans );
353  AbstractLinAlgPack::Mp_StMtM( C, a, *this, A_trans, B_dense(), B_trans, b );
354  }
355 }
356 
357 void MatrixOpSerial::Mp_StMtM(DMatrixSlice* gms_lhs, value_type alpha
358  , BLAS_Cpp::Transp trans_rhs1, const DMatrixSliceSym& sym_rhs2
359  , BLAS_Cpp::Transp trans_rhs2, value_type beta) const
360 {
361  TEUCHOS_TEST_FOR_EXCEPT(true); // Todo: Implement!
362 }
363 
364 void MatrixOpSerial::Mp_StMtM(DMatrixSlice* gms_lhs, value_type alpha, const DMatrixSliceSym& sym_rhs1
365  , BLAS_Cpp::Transp trans_rhs1, BLAS_Cpp::Transp trans_rhs2, value_type beta) const
366 {
367  TEUCHOS_TEST_FOR_EXCEPT(true); // Todo: Implement!
368 }
369 
370 void MatrixOpSerial::Mp_StMtM(DMatrixSlice* gms_lhs, value_type alpha
371  , BLAS_Cpp::Transp trans_rhs1, const DMatrixSliceTri& tri_rhs2
372  , BLAS_Cpp::Transp trans_rhs2, value_type beta) const
373 {
374  TEUCHOS_TEST_FOR_EXCEPT(true); // Todo: Implement!
375 }
376 
377 void MatrixOpSerial::Mp_StMtM(DMatrixSlice* gms_lhs, value_type alpha, const DMatrixSliceTri& tri_rhs1
378  , BLAS_Cpp::Transp trans_rhs1, BLAS_Cpp::Transp trans_rhs2, value_type beta) const
379 {
380  TEUCHOS_TEST_FOR_EXCEPT(true); // Todo: Implement!
381 }
382 
384  BLAS_Cpp::Transp M_trans, value_type a
385  , value_type b, DMatrixSliceSym* S ) const
386 {
387  using BLAS_Cpp::no_trans;
388  using BLAS_Cpp::trans;
389  using BLAS_Cpp::trans_not;
390  using BLAS_Cpp::rows;
391  using BLAS_Cpp::cols;
392  using Teuchos::Workspace;
394  //
395  // S = b*S + a*op(M)*op(M')
396  //
397  const size_type
398  M_rows = this->rows(),
399  M_cols = this->cols(),
400  opM_rows = rows( M_rows, M_cols, M_trans ),
401  opM_cols = cols( M_rows, M_cols, M_trans ),
402  m = opM_rows;
403  DenseLinAlgPack::Mp_MtM_assert_sizes(
404  S->rows(), S->cols(), no_trans
405  ,M_rows, M_cols, M_trans
406  ,M_rows, M_cols, trans_not(M_trans) );
407  // S = b*S
408  DenseLinAlgPack::Mt_S( &DenseLinAlgPack::nonconst_tri_ele(S->gms(),S->uplo()), b );
409  //
410  // Form this matrix one column at a time by multiplying by e(j):
411  //
412  // S(:,j) += a*op(M)*(op(M')*e(j))
413  //
414  // j = 1 ... opM_rows
415  //
416  Workspace<value_type> t1_ws(wss,opM_cols),
417  t2_ws(wss,opM_rows);
418  DVectorSlice t1(&t1_ws[0],t1_ws.size()),
419  t2(&t2_ws[0],t2_ws.size());
420  for( size_type j = 1; j <= opM_rows; ++j ) {
421  EtaVector e_j(j,opM_rows);
422  LinAlgOpPack::V_MtV(&t1,*this,trans_not(M_trans),e_j()); // t1 = op(M')*e(j)
423  LinAlgOpPack::V_MtV(&t2,*this,M_trans,t1); // t2 = op(M)*t1
424  // S(j,:) += a*t2'
425  if( S->uplo() == BLAS_Cpp::upper )
426  DenseLinAlgPack::Vp_StV( &S->gms().row(j)(1,j), a, t2(1,j) );
427  else
428  DenseLinAlgPack::Vp_StV( &S->gms().row(j)(j,m), a, t2(j,m) );
429  // S(:,j) += a*t2
430  if( S->uplo() == BLAS_Cpp::upper )
431  DenseLinAlgPack::Vp_StV( &S->gms().col(j)(1,j), a, t2(1,j) );
432  else
433  DenseLinAlgPack::Vp_StV( &S->gms().col(j)(j,m), a, t2(j,m) );
434  }
435 }
436 
437 // Overridden from MatrixOp
438 
440 {
441  const size_type rows = this->rows();
442  if(space_cols_.dim() != rows)
443  space_cols_.initialize(rows);
444  return space_cols_;
445 }
446 
448 {
449  const size_type cols = this->cols();
450  if(space_rows_.dim() != cols)
451  space_rows_.initialize(cols);
452  return space_rows_;
453 }
454 
455 std::ostream& MatrixOpSerial::output(std::ostream& out) const {
456  DMatrix tmp( 0.0, rows(), cols() );
457  this->Mp_StM( &tmp(), 1.0 , BLAS_Cpp::no_trans );
458  return out << tmp();
459 }
460 
462  MatrixOp* mwo_lhs, value_type alpha
463  ,BLAS_Cpp::Transp trans_rhs
464  ) const
465 {
467  *mwo_gms_lhs = dynamic_cast<MatrixOpGetGMSMutable*>(mwo_lhs);
468  if(!mwo_gms_lhs)
469  return MatrixOp::Mp_StM(mwo_lhs,alpha,trans_rhs); // boot it!
470  this->Mp_StM( &MatrixDenseMutableEncap(mwo_gms_lhs)(), alpha, trans_rhs );
471  return true;
472 }
473 
475  MatrixOp* mwo_lhs, value_type alpha
476  , BLAS_Cpp::Transp M_trans
477  , const GenPermMatrixSlice& P_rhs, BLAS_Cpp::Transp P_rhs_trans
478  ) const
479 {
481  *mwo_gms_lhs = dynamic_cast<MatrixOpGetGMSMutable*>(mwo_lhs);
482  if(!mwo_gms_lhs)
483  return MatrixOp::Mp_StMtP(mwo_lhs,alpha,M_trans,P_rhs,P_rhs_trans); // boot it!
484  this->Mp_StMtP(&MatrixDenseMutableEncap(mwo_gms_lhs)(),alpha,M_trans,P_rhs,P_rhs_trans);
485  return true;
486 }
487 
489  MatrixOp* mwo_lhs, value_type alpha
490  , const GenPermMatrixSlice& P_rhs, BLAS_Cpp::Transp P_rhs_trans
491  , BLAS_Cpp::Transp M_trans
492  ) const
493 {
495  *mwo_gms_lhs = dynamic_cast<MatrixOpGetGMSMutable*>(mwo_lhs);
496  if(!mwo_gms_lhs)
497  return MatrixOp::Mp_StPtM(mwo_lhs,alpha,P_rhs,P_rhs_trans,M_trans); // boot it!
498  this->Mp_StPtM(&MatrixDenseMutableEncap(mwo_gms_lhs)(),alpha,P_rhs,P_rhs_trans,M_trans);
499  return true;
500 }
501 
503  MatrixOp* mwo_lhs, value_type alpha
504  ,const GenPermMatrixSlice& P_rhs1, BLAS_Cpp::Transp P_rhs1_trans
505  ,BLAS_Cpp::Transp M_trans
506  ,const GenPermMatrixSlice& P_rhs2, BLAS_Cpp::Transp P_rhs2_trans
507  ) const
508 {
510  *mwo_gms_lhs = dynamic_cast<MatrixOpGetGMSMutable*>(mwo_lhs);
511  if(!mwo_gms_lhs)
512  return MatrixOp::Mp_StPtMtP(mwo_lhs,alpha,P_rhs1,P_rhs1_trans,M_trans,P_rhs2,P_rhs2_trans); // boot it!
513  this->Mp_StPtMtP(&MatrixDenseMutableEncap(mwo_gms_lhs)(),alpha,P_rhs1,P_rhs1_trans,M_trans,P_rhs2,P_rhs2_trans);
514  return true;
515 }
516 
518  VectorMutable* v_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs1
519  , const Vector& v_rhs2, value_type beta) const
520 {
521  VectorDenseMutableEncap vs_lhs(*v_lhs);
522  VectorDenseEncap vs_rhs2(v_rhs2);
523  this->Vp_StMtV( &vs_lhs(), alpha, trans_rhs1, vs_rhs2(), beta );
524 }
525 
527  VectorMutable* v_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs1
528  , const SpVectorSlice& sv_rhs2, value_type beta) const
529 {
530  VectorDenseMutableEncap vs_lhs(*v_lhs);
531  this->Vp_StMtV( &vs_lhs(), alpha, trans_rhs1, sv_rhs2, beta );
532 }
533 
535  VectorMutable* v_lhs, value_type alpha
536  , const GenPermMatrixSlice& P_rhs1, BLAS_Cpp::Transp P_rhs1_trans
537  , BLAS_Cpp::Transp M_rhs2_trans
538  , const Vector& v_rhs3, value_type beta) const
539 {
540  VectorDenseMutableEncap vs_lhs(*v_lhs);
541  VectorDenseEncap vs_rhs3(v_rhs3);
542  this->Vp_StPtMtV( &vs_lhs(), alpha, P_rhs1, P_rhs1_trans, M_rhs2_trans, vs_rhs3(), beta );
543 }
544 
546  VectorMutable* v_lhs, value_type alpha
547  , const GenPermMatrixSlice& P_rhs1, BLAS_Cpp::Transp P_rhs1_trans
548  , BLAS_Cpp::Transp M_rhs2_trans
549  , const SpVectorSlice& sv_rhs3, value_type beta) const
550 {
551  VectorDenseMutableEncap vs_lhs(*v_lhs);
552  this->Vp_StPtMtV( &vs_lhs(), alpha, P_rhs1, P_rhs1_trans, M_rhs2_trans, sv_rhs3, beta );
553 }
554 
556  const Vector& v_rhs1, BLAS_Cpp::Transp trans_rhs2
557  , const Vector& v_rhs3) const
558 {
559  VectorDenseEncap vs_rhs1(v_rhs1);
560  VectorDenseEncap vs_rhs3(v_rhs3);
561  return this->transVtMtV(vs_rhs1(),trans_rhs2,vs_rhs3());
562 }
563 
565  BLAS_Cpp::Transp M_trans, value_type alpha
566  , const GenPermMatrixSlice& P1, BLAS_Cpp::Transp P1_trans
567  , const GenPermMatrixSlice& P2, BLAS_Cpp::Transp P2_trans
568  , value_type beta, MatrixSymOp* symwo_lhs ) const
569 {
571  *symwo_gms_lhs = dynamic_cast<MatrixSymOpGetGMSSymMutable*>(symwo_lhs);
572  if(!symwo_gms_lhs) {
573  MatrixOp::syr2k(M_trans,alpha,P1,P1_trans,P2,P2_trans,beta,symwo_lhs); // Boot it
574  return;
575  }
576  this->syr2k(
577  M_trans,alpha,P1,P1_trans,P2,P2_trans,beta
578  ,&MatrixDenseSymMutableEncap(symwo_gms_lhs)()
579  );
580 }
581 
583  MatrixOp* mwo_lhs, value_type alpha
584  , BLAS_Cpp::Transp trans_rhs1, const MatrixOp& mwo_rhs2
585  , BLAS_Cpp::Transp trans_rhs2, value_type beta ) const
586 {
588  *mwo_gms_lhs = dynamic_cast<MatrixOpGetGMSMutable*>(mwo_lhs);
589  if(mwo_gms_lhs) {
590  // Try to match the rhs arguments to known serial interfaces
591  if(const MatrixOpGetGMS* mwo_gms_rhs2 = dynamic_cast<const MatrixOpGetGMS*>(&mwo_rhs2)) {
592  this->Mp_StMtM(
593  &MatrixDenseMutableEncap(mwo_gms_lhs)(),alpha,trans_rhs1
594  ,MatrixDenseEncap(*mwo_gms_rhs2)(),trans_rhs2,beta );
595  return true;
596  }
597  if(const MatrixSymOpGetGMSSym* mwo_sym_gms_rhs2 = dynamic_cast<const MatrixSymOpGetGMSSym*>(&mwo_rhs2)) {
598  this->Mp_StMtM(
599  &MatrixDenseMutableEncap(mwo_gms_lhs)(),alpha,trans_rhs1
600  ,MatrixDenseEncap(*mwo_sym_gms_rhs2)(),trans_rhs2,beta );
601  return true;
602  }
603  if(const MatrixOpGetGMSTri* mwo_tri_gms_rhs2 = dynamic_cast<const MatrixOpGetGMSTri*>(&mwo_rhs2)) {
604  this->Mp_StMtM(
605  &MatrixDenseMutableEncap(mwo_gms_lhs)(),alpha,trans_rhs1
606  ,MatrixDenseEncap(*mwo_tri_gms_rhs2)(),trans_rhs2,beta );
607  return true;
608  }
609  // If we get here, the matrix arguments did not match up so we have to give up (I think?)
610  }
611  // Let the default implementation try to find matrix arguments that can handle this!
612  return MatrixOp::Mp_StMtM(mwo_lhs,alpha,trans_rhs1,mwo_rhs2,trans_rhs2,beta); // Boot it!
613 }
614 
616  BLAS_Cpp::Transp M_trans, value_type alpha
617  , value_type beta, MatrixSymOp* symwo_lhs ) const
618 {
620  *symwo_gms_lhs = dynamic_cast<MatrixSymOpGetGMSSymMutable*>(symwo_lhs);
621  if(!symwo_gms_lhs) {
622  return MatrixOp::syrk(M_trans,alpha,beta,symwo_lhs); // Boot it
623  }
624  this->syrk(M_trans,alpha,beta,&MatrixDenseSymMutableEncap(symwo_gms_lhs)());
625  return true;
626 }
627 
628 } // end namespace AbstractLinAlgPack
std::ostream & output(std::ostream &out) const
friend void syr2k(const MatrixOp &M, BLAS_Cpp::Transp M_trans, value_type alpha, const GenPermMatrixSlice &P1, BLAS_Cpp::Transp P1_trans, const GenPermMatrixSlice &P2, BLAS_Cpp::Transp P2_trans, value_type beta, MatrixSymOp *symwo_lhs)
Abstract interface for immutable, finite dimensional, coordinate vectors {abstract}.
virtual void Mp_StPtM(DMatrixSlice *gms_lhs, value_type alpha, const GenPermMatrixSlice &P_rhs, BLAS_Cpp::Transp P_rhs_trans, BLAS_Cpp::Transp M_trans) const
gms_lhs += alpha * op(P) * op(M_rhs)
size_type dim() const
Return the number of elements in the full vector.
void Mp_StMtM(MatrixOp *mwo_lhs, value_type alpha, const MatrixOp &mwo_rhs1, BLAS_Cpp::Transp trans_rhs1, const MatrixOp &mwo_rhs2, BLAS_Cpp::Transp trans_rhs2, value_type beta=1.0)
mwo_lhs = alpha * op(mwo_rhs1) * op(mwo_rhs2) + beta * mwo_lhs (right) (xGEMM).
void Vt_S(VectorMutable *v_lhs, const value_type &alpha)
v_lhs *= alpha
Mix-in interface that allows the extraction of a const DenseLinAlgPack::DMatrixSliceTri view of an no...
void Vp_StPtMtV(VectorMutable *v_lhs, value_type alpha, const GenPermMatrixSlice &P_rhs1, BLAS_Cpp::Transp P_rhs1_trans, const MatrixOp &M_rhs2, BLAS_Cpp::Transp M_rhs2_trans, const Vector &v_rhs3, value_type beta=1.0)
v_lhs = alpha * op(P_rhs1) * op(M_rhs2) * v_rhs3 + beta * v_rhs
virtual void Mp_StMtP(DMatrixSlice *gms_lhs, value_type alpha, BLAS_Cpp::Transp M_trans, const GenPermMatrixSlice &P_rhs, BLAS_Cpp::Transp P_rhs_trans) const
gms_lhs += alpha * op(M_rhs) * op(P_rhs)
void Vp_StV(VectorMutable *v_lhs, const value_type &alpha, const Vector &v_rhs)
v_lhs = alpha * v_rhs + v_lhs
size_type rows(size_type rows, size_type cols, BLAS_Cpp::Transp _trans)
void assign(VectorMutable *v_lhs, const V &V_rhs)
v_lhs = V_rhs.
Abstract interface that allows the extraction of a non-const DenseLinAlgPack::DMatrixSliceSym view of...
const_iterator end() const
Return the end of this->const_iterator_begin().
virtual void syr2k(BLAS_Cpp::Transp M_trans, value_type alpha, const GenPermMatrixSlice &P1, BLAS_Cpp::Transp P1_trans, const GenPermMatrixSlice &P2, BLAS_Cpp::Transp P2_trans, value_type beta, DMatrixSliceSym *sym_lhs) const
Perform a specialized rank-2k update of a dense symmetric matrix of the form:
BLAS_Cpp::Uplo uplo() const
bool is_sorted() const
Return true if the sequence is assumed sorted.
friend void Mp_StPtMtP(MatrixOp *mwo_lhs, value_type alpha, const GenPermMatrixSlice &P_rhs1, BLAS_Cpp::Transp P_rhs1_trans, const MatrixOp &M_rhs, BLAS_Cpp::Transp trans_rhs, const GenPermMatrixSlice &P_rhs2, BLAS_Cpp::Transp P_rhs2_trans)
Base class for all matrices implemented in a shared memory address space.
virtual void Mp_StM(DMatrixSlice *gms_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs) const
gms_lhs += alpha * op(M_rhs) (BLAS xAXPY)
Interface adding operations specific for a symmetric matrix {abstract}.
virtual value_type transVtMtV(const DVectorSlice &vs_rhs1, BLAS_Cpp::Transp trans_rhs2, const DVectorSlice &vs_rhs3) const
result = vs_rhs1' * op(M_rhs2) * vs_rhs3
virtual void Mp_StMtM(DMatrixSlice *gms_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs1, const DMatrixSlice &gms_rhs2, BLAS_Cpp::Transp trans_rhs2, value_type beta) const
gms_lhs = alpha * op(M_rhs1) * op(gms_rhs2) + beta * gms_lhs (right) (xGEMM)
virtual size_type cols() const
Return the number of columns in the matrix.
This is a full random access iterator for accessing row and colunmn indices.
Abstract interface for objects that represent a space for mutable coordinate vectors.
const_iterator begin() const
Return a random access iterator for accessing which row and column that each nonzero 1...
friend void Mp_StMtM(MatrixOp *mwo_lhs, value_type alpha, const MatrixOp &mwo_rhs1, BLAS_Cpp::Transp trans_rhs1, const MatrixOp &mwo_rhs2, BLAS_Cpp::Transp trans_rhs2, value_type beta)
Extract a constant DenseLinAlgPack::DVectorSlice view of a Vector object.
Helper class type that simplifies the usage of the MatrixSymOpGetGMSSymMutable interface for clients...
index_type dim() const
Returns 0 if uninitialized.
Create an eta vector (scaled by alpha = default 1).
void Mp_StM(MatrixOp *mwo_lhs, value_type alpha, const MatrixOp &M_rhs, BLAS_Cpp::Transp trans_rhs)
void Vp_StMtV(VectorMutable *v_lhs, value_type alpha, const MatrixOp &M_rhs1, BLAS_Cpp::Transp trans_rhs1, const Vector &v_rhs2, value_type beta=1.0)
v_lhs = alpha * op(M_rhs1) * v_rhs2 + beta * v_lhs (BLAS xGEMV)
void uninitialized_resize(size_type size, size_type nz, size_type max_nz, difference_type offset=0)
Resize to #size# with a max_nz# uninitialized non-zero elements.
friend void Mp_StM(MatrixOp *mwo_lhs, value_type alpha, const MatrixOp &M_rhs, BLAS_Cpp::Transp trans_rhs)
Abstract interface that allows the extraction of a const DMatrixSlice view of an abstract matrix...
Helper class type that simplifies the usage of the MatrixOpGetGMS interface for clients.
EOverLap overlap(const SparseVectorSlice< T_Element > &sv) const
Base class for all matrices that support basic matrix operations.
virtual void syrk(BLAS_Cpp::Transp M_trans, value_type alpha, value_type beta, DMatrixSliceSym *sym_lhs) const
Perform a rank-k update of a dense symmetric matrix of the form:
void Mt_S(MatrixOp *mwo_lhs, value_type alpha)
virtual void Vp_StPtMtV(DVectorSlice *vs_lhs, value_type alpha, const GenPermMatrixSlice &P_rhs1, BLAS_Cpp::Transp P_rhs1_trans, BLAS_Cpp::Transp M_rhs2_trans, const DVectorSlice &vs_rhs3, value_type beta) const
vs_lhs = alpha * op(P_rhs1) * op(M_rhs2) * vs_rhs3 + beta * vs_rhs
value_type dot(const Vector &v_rhs1, const Vector &v_rhs2)
result = v_rhs1' * v_rhs2
Abstract interface that allows the extraction of a const DenseLinAlgPack::DMatrixSliceSym view of an ...
Transp trans_not(Transp _trans)
void initialize(size_type dim)
Initialize given the dimension of the vector space.
virtual void Vp_StMtV(DVectorSlice *vs_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs1, const DVectorSlice &vs_rhs2, value_type beta) const =0
vs_lhs = alpha * op(M_rhs1) * vs_rhs2 + beta * vs_lhs (BLAS xGEMV)
friend void Mp_StMtP(MatrixOp *mwo_lhs, value_type alpha, const MatrixOp &M_rhs, BLAS_Cpp::Transp M_trans, const GenPermMatrixSlice &P_rhs, BLAS_Cpp::Transp P_rhs_trans)
Abstract interface that allows the extraction of a non-const DMatrixSlice view of an abstract matrix...
virtual void Mp_StPtMtP(DMatrixSlice *gms_lhs, value_type alpha, const GenPermMatrixSlice &P_rhs1, BLAS_Cpp::Transp P_rhs1_trans, BLAS_Cpp::Transp M_trans, const GenPermMatrixSlice &P_rhs2, BLAS_Cpp::Transp P_rhs2_trans) const
gms_lhs += alpha * op(P_rhs1) * op(M_rhs) * op(P_rhs2)
Abstract interface for mutable coordinate vectors {abstract}.
void initialize(index_type index, value_type value)
Initialize.
Extract a non-const DenseLinAlgPack::DVectorSlice view of a VectorMutable object. ...
virtual size_type rows() const
Return the number of rows in the matrix.
friend void syrk(const MatrixOp &mwo_rhs, BLAS_Cpp::Transp M_trans, value_type alpha, value_type beta, MatrixSymOp *sym_lhs)
DVectorSlice col(size_type j)
Transp
size_type cols(size_type rows, size_type cols, BLAS_Cpp::Transp _trans)
iterator begin()
Returns iterator that iterates forward through the nonzero elements.
size_type nz() const
Return the number of non-zero elements.
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
TEUCHOSCORE_LIB_DLL_EXPORT Teuchos::RCP< WorkspaceStore > get_default_workspace_store()
Concrete matrix type to represent general permutation (mapping) matrices.
friend void Mp_StPtM(MatrixOp *mwo_lhs, value_type alpha, const GenPermMatrixSlice &P_rhs, BLAS_Cpp::Transp P_rhs_trans, const MatrixOp &M_rhs, BLAS_Cpp::Transp M_trans)
Helper class type that simplifies the usage of the MatrixOpGetGMSMutable interface for clients...
DVectorSlice row(size_type i)