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_MatrixOpNonsingAggr.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 #include "AbstractLinAlgPack_MatrixOpNonsingAggr.hpp"
43 #include "AbstractLinAlgPack_MatrixOpOut.hpp"
44 #include "AbstractLinAlgPack_VectorSpace.hpp"
45 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
46 #include "Teuchos_Assert.hpp"
47 #include "Teuchos_dyn_cast.hpp"
48 
49 namespace AbstractLinAlgPack {
50 
51 // Constructors / initializers
52 
54 {} // Nothing to explicitly initialize
55 
57  const mwo_ptr_t &mwo
58  ,BLAS_Cpp::Transp mwo_trans
59  ,const mns_ptr_t &mns
60  ,BLAS_Cpp::Transp mns_trans
61  )
62 {
63  this->initialize(mwo,mwo_trans,mns,mns_trans);
64 }
65 
67  const mwo_ptr_t &mwo
68  ,BLAS_Cpp::Transp mwo_trans
69  ,const mns_ptr_t &mns
70  ,BLAS_Cpp::Transp mns_trans
71  )
72 {
73 #ifdef TEUCHOS_DEBUG
75  mwo.get() == NULL, std::invalid_argument
76  ,"MatrixOpNonsingAggr::initialize(...): Error!" );
78  mns.get() == NULL, std::invalid_argument
79  ,"MatrixOpNonsingAggr::initialize(...): Error!" );
80  const size_type
81  mwo_rows = mwo->rows(),
82  mwo_cols = mwo->cols(),
83  mns_rows = mns->rows(),
84  mns_cols = mns->cols();
86  mwo_rows != mwo_cols, std::invalid_argument
87  ,"MatrixOpNonsingAggr::initialize(...): Error!" );
89  mns_rows != mns_cols, std::invalid_argument
90  ,"MatrixOpNonsingAggr::initialize(...): Error!" );
92  mwo_rows != mns_rows, std::invalid_argument
93  ,"MatrixOpNonsingAggr::initialize(...): Error!" );
94 #endif
95  mwo_ = mwo;
96  mwo_trans_ = mwo_trans;
97  mns_ = mns;
98  mns_trans_ = mns_trans;
99 }
100 
102 {
103  namespace rcp = MemMngPack;
104  mwo_ = Teuchos::null;
105  mwo_trans_ = BLAS_Cpp::no_trans;
106  mns_ = Teuchos::null;
107  mns_trans_ = BLAS_Cpp::no_trans;
108 }
109 
110 // Overridden from MatrixBase
111 
112 size_type MatrixOpNonsingAggr::rows() const
113 {
114  return mwo_.get() ? mwo_->rows() : 0; // square matrix!
115 }
116 
117 size_type MatrixOpNonsingAggr::cols() const
118 {
119  return mwo_.get() ? mwo_->rows() : 0; // square matrix!
120 }
121 
122 size_type MatrixOpNonsingAggr::nz() const
123 {
124  return mwo_.get() ? mwo_->nz() : 0;
125 }
126 
127 // Overridden from MatrixOp
128 
130 {
131  return mwo_trans_ == BLAS_Cpp::no_trans ? mwo_->space_cols() : mwo_->space_rows();
132 }
133 
135 {
136  return mwo_trans_ == BLAS_Cpp::no_trans ? mwo_->space_rows() : mwo_->space_cols();
137 }
138 
139 MatrixOp::mat_ptr_t
140 MatrixOpNonsingAggr::sub_view(const Range1D& row_rng, const Range1D& col_rng) const
141 {
142  return MatrixOp::sub_view(row_rng,col_rng); // ToDo: Speicalize!
143 }
144 
146 {
147  using Teuchos::dyn_cast;
148  const MatrixOpNonsingAggr
149  Mp = dyn_cast<const MatrixOpNonsingAggr>(M);
150  if( this == &Mp )
151  return *this; // Assignment to self
152  // Shallow copy is okay as long as client is careful!
153  mwo_ = Mp.mwo_;
154  mwo_trans_ = Mp.mwo_trans_;
155  mns_ = Mp.mns_;
156  mns_trans_ = Mp.mns_trans_;
157  return *this;
158 }
159 
160 std::ostream& MatrixOpNonsingAggr::output(std::ostream& out) const
161 {
162  out << "Aggregate nonsingular matrix:\n";
163  out << (mwo_trans_ == BLAS_Cpp::no_trans ? "mwo" : "mwo\'") << " =\n" << *mwo_;
164  out << (mns_trans_ == BLAS_Cpp::no_trans ? "mns" : "mns\'") << " = ???\n";
165  return out;
166 }
167 
169  MatrixOp* mwo_lhs, value_type alpha
170  , BLAS_Cpp::Transp trans_rhs) const
171 {
172  AbstractLinAlgPack::Mp_StM(mwo_lhs,alpha,*mwo_,BLAS_Cpp::trans_trans(mwo_trans_,trans_rhs));
173  return true;
174 }
175 
177  MatrixOp* mwo_lhs, value_type alpha
178  , BLAS_Cpp::Transp M_trans
179  , const GenPermMatrixSlice& P_rhs, BLAS_Cpp::Transp P_rhs_trans
180  ) const
181 {
183  mwo_lhs,alpha,*mwo_,BLAS_Cpp::trans_trans(mwo_trans_,M_trans)
184  ,P_rhs,P_rhs_trans);
185  return true;
186 }
187 
189  MatrixOp* mwo_lhs, value_type alpha
190  , const GenPermMatrixSlice& P_rhs, BLAS_Cpp::Transp P_rhs_trans
191  , BLAS_Cpp::Transp M_trans
192  ) const
193 {
195  mwo_lhs,alpha,P_rhs,P_rhs_trans
196  ,*mwo_,BLAS_Cpp::trans_trans(mwo_trans_,M_trans));
197  return true;
198 }
199 
201  MatrixOp* mwo_lhs, value_type alpha
202  ,const GenPermMatrixSlice& P_rhs1, BLAS_Cpp::Transp P_rhs1_trans
203  ,BLAS_Cpp::Transp M_trans
204  ,const GenPermMatrixSlice& P_rhs2, BLAS_Cpp::Transp P_rhs2_trans
205  ) const
206 {
208  mwo_lhs,alpha,P_rhs1,P_rhs1_trans
209  ,*mwo_,BLAS_Cpp::trans_trans(mwo_trans_,M_trans)
210  ,P_rhs2,P_rhs2_trans);
211  return true;
212 }
213 
215  VectorMutable* y, value_type a, BLAS_Cpp::Transp M_trans
216  , const Vector& x, value_type b) const
217 {
218  AbstractLinAlgPack::Vp_StMtV(y,a,*mwo_,BLAS_Cpp::trans_trans(mwo_trans_,M_trans),x,b);
219 }
220 
222  VectorMutable* y, value_type a, BLAS_Cpp::Transp M_trans
223  , const SpVectorSlice& x, value_type b) const
224 {
225  AbstractLinAlgPack::Vp_StMtV(y,a,*mwo_,BLAS_Cpp::trans_trans(mwo_trans_,M_trans),x,b);
226 }
227 
229  VectorMutable* vs_lhs, value_type alpha
230  , const GenPermMatrixSlice& P_rhs1, BLAS_Cpp::Transp P_rhs1_trans
231  , BLAS_Cpp::Transp M_rhs2_trans
232  , const Vector& v_rhs3, value_type beta) const
233 {
235  vs_lhs,alpha,P_rhs1,P_rhs1_trans
236  ,*mwo_,BLAS_Cpp::trans_trans(mwo_trans_,M_rhs2_trans)
237  ,v_rhs3,beta);
238 }
239 
241  VectorMutable* vs_lhs, value_type alpha
242  , const GenPermMatrixSlice& P_rhs1, BLAS_Cpp::Transp P_rhs1_trans
243  , BLAS_Cpp::Transp M_rhs2_trans
244  , const SpVectorSlice& sv_rhs3, value_type beta) const
245 {
247  vs_lhs,alpha,P_rhs1,P_rhs1_trans
248  ,*mwo_,BLAS_Cpp::trans_trans(mwo_trans_,M_rhs2_trans)
249  ,sv_rhs3,beta);
250 }
251 
253  const Vector& v_rhs1, BLAS_Cpp::Transp trans_rhs2
254  , const Vector& v_rhs3) const
255 {
256  return AbstractLinAlgPack::transVtMtV(v_rhs1,*mwo_,BLAS_Cpp::trans_trans(mwo_trans_,trans_rhs2),v_rhs3);
257 }
258 
260  const SpVectorSlice& sv_rhs1, BLAS_Cpp::Transp trans_rhs2
261  ,const SpVectorSlice& sv_rhs3
262  ) const
263 {
264  return AbstractLinAlgPack::transVtMtV(sv_rhs1,*mwo_,BLAS_Cpp::trans_trans(mwo_trans_,trans_rhs2),sv_rhs3);
265 }
266 
268  BLAS_Cpp::Transp M_trans, value_type alpha
269  ,const GenPermMatrixSlice& P1, BLAS_Cpp::Transp P1_trans
270  ,const GenPermMatrixSlice& P2, BLAS_Cpp::Transp P2_trans
271  ,value_type beta, MatrixSymOp* symwo_lhs
272  ) const
273 {
275  *mwo_,BLAS_Cpp::trans_trans(mwo_trans_,M_trans)
276  ,alpha,P1,P1_trans,P2,P2_trans,beta,symwo_lhs);
277 }
278 
280  MatrixOp* mwo_lhs, value_type alpha
281  ,BLAS_Cpp::Transp trans_rhs1, const MatrixOp& mwo_rhs2
282  ,BLAS_Cpp::Transp trans_rhs2, value_type beta
283  ) const
284 {
286  mwo_lhs,alpha,*mwo_,trans_rhs1
287  ,mwo_rhs2,BLAS_Cpp::trans_trans(mwo_trans_,trans_rhs2),beta);
288  return true;
289 }
290 
292  MatrixOp* mwo_lhs, value_type alpha
293  ,const MatrixOp& mwo_rhs1, BLAS_Cpp::Transp trans_rhs1
294  ,BLAS_Cpp::Transp trans_rhs2, value_type beta
295  ) const
296 {
298  mwo_lhs,alpha,mwo_rhs1,trans_rhs1
299  ,*mwo_,BLAS_Cpp::trans_trans(mwo_trans_,trans_rhs2),beta);
300  return true;
301 }
302 
304  BLAS_Cpp::Transp M_trans, value_type alpha
305  ,value_type beta, MatrixSymOp* sym_lhs
306  ) const
307 {
308  AbstractLinAlgPack::syrk(*mwo_,BLAS_Cpp::trans_trans(mwo_trans_,M_trans),alpha,beta,sym_lhs);
309  return true;
310 }
311 
312 // Overridden from MatrixNonsing */
313 
315  VectorMutable* v_lhs, BLAS_Cpp::Transp trans_rhs1
316  ,const Vector& v_rhs2
317  ) const
318 {
319  AbstractLinAlgPack::V_InvMtV(v_lhs,*mns_,BLAS_Cpp::trans_trans(mns_trans_,trans_rhs1),v_rhs2);
320 }
321 
323  VectorMutable* v_lhs, BLAS_Cpp::Transp trans_rhs1
324  ,const SpVectorSlice& sv_rhs2
325  ) const
326 {
327  AbstractLinAlgPack::V_InvMtV(v_lhs,*mns_,BLAS_Cpp::trans_trans(mns_trans_,trans_rhs1),sv_rhs2);
328 }
329 
331  const Vector& v_rhs1
332  ,BLAS_Cpp::Transp trans_rhs2, const Vector& v_rhs3
333  ) const
334 {
335  return AbstractLinAlgPack::transVtInvMtV(v_rhs1,*mns_,BLAS_Cpp::trans_trans(mns_trans_,trans_rhs2),v_rhs3);
336 }
337 
339  const SpVectorSlice& sv_rhs1
340  ,BLAS_Cpp::Transp trans_rhs2, const SpVectorSlice& sv_rhs3
341  ) const
342 {
343  return AbstractLinAlgPack::transVtInvMtV(sv_rhs1,*mns_,BLAS_Cpp::trans_trans(mns_trans_,trans_rhs2),sv_rhs3);
344 }
345 
347  MatrixOp* m_lhs, value_type alpha
348  ,BLAS_Cpp::Transp trans_rhs1
349  ,const MatrixOp& mwo_rhs2, BLAS_Cpp::Transp trans_rhs2
350  ) const
351 {
352  AbstractLinAlgPack::M_StInvMtM(m_lhs,alpha,*mns_,BLAS_Cpp::trans_trans(mns_trans_,trans_rhs1),mwo_rhs2,trans_rhs2);
353 }
354 
356  MatrixOp* m_lhs, value_type alpha
357  ,const MatrixOp& mwo_rhs1, BLAS_Cpp::Transp trans_rhs1
358  ,BLAS_Cpp::Transp trans_rhs2
359  ) const
360 {
361  AbstractLinAlgPack::M_StMtInvM(m_lhs,alpha,mwo_rhs1,trans_rhs1,*mns_,BLAS_Cpp::trans_trans(mns_trans_,trans_rhs1));
362 }
363 
364 } // end namespace AbstractLinAlgPack
bool Mp_StPtMtP(MatrixOp *mwo_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
value_type transVtInvMtV(const Vector &v_rhs1, const MatrixNonsing &M_rhs2, BLAS_Cpp::Transp trans_rhs2, const Vector &v_rhs3)
result = v_rhs1' * inv(op(M_rhs2)) * v_rhs3
Abstract interface for immutable, finite dimensional, coordinate vectors {abstract}.
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).
bool Mp_StPtM(MatrixOp *mwo_lhs, value_type alpha, const GenPermMatrixSlice &P_rhs, BLAS_Cpp::Transp P_rhs_trans, BLAS_Cpp::Transp M_trans) const
void V_InvMtV(VectorMutable *v_lhs, BLAS_Cpp::Transp trans_rhs1, const Vector &v_rhs2) const
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 mat_ptr_t sub_view(const Range1D &row_rng, const Range1D &col_rng) const
Create a transient constant sub-matrix view of this matrix (if supported).
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
T_To & dyn_cast(T_From &from)
void M_StMtInvM(MatrixOp *m_lhs, value_type alpha, const MatrixOp &mwo_rhs1, BLAS_Cpp::Transp trans_rhs1, const MatrixNonsing &M_rhs2, BLAS_Cpp::Transp trans_rhs2)
m_lhs = alpha * op(mwo_rhs1) * inv(op(M_rhs2)) (left)
void Vp_StPtMtV(VectorMutable *vs_lhs, value_type alpha, const GenPermMatrixSlice &P_rhs1, BLAS_Cpp::Transp P_rhs1_trans, BLAS_Cpp::Transp M_rhs2_trans, const Vector &v_rhs3, value_type beta) const
const mns_ptr_t & mns() const
T * get() const
void initialize(const mwo_ptr_t &mwo, BLAS_Cpp::Transp mwo_trans, const mns_ptr_t &mns, BLAS_Cpp::Transp mns_trans)
Initialize.
MatrixOp::mat_ptr_t sub_view(const Range1D &row_rng, const Range1D &col_rng) const
bool syrk(BLAS_Cpp::Transp M_trans, value_type alpha, value_type beta, MatrixSymOp *sym_lhs) const
void Vp_StMtV(VectorMutable *vs_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs1, const Vector &v_rhs2, value_type beta) const
value_type transVtMtV(const Vector &v_rhs1, const MatrixOp &M_rhs2, BLAS_Cpp::Transp trans_rhs2, const Vector &v_rhs3)
result = v_rhs1' * op(M_rhs2) * v_rhs3
Interface adding operations specific for a symmetric matrix {abstract}.
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)
symwo_lhs += alpha*op(P1')*op(M)*op(P2) + alpha*op(P2')*op(M')*op(P1) + beta*symwo_lhs ...
value_type transVtMtV(const Vector &v_rhs1, BLAS_Cpp::Transp trans_rhs2, const Vector &v_rhs3) const
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)
mwo_lhs += alpha * op(P_rhs1) * op(M_rhs) * op(P_rhs2).
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
value_type transVtInvMtV(const Vector &v_rhs1, BLAS_Cpp::Transp trans_rhs2, const Vector &v_rhs3) const
bool Mp_StMtP(MatrixOp *mwo_lhs, value_type alpha, BLAS_Cpp::Transp M_trans, const GenPermMatrixSlice &P_rhs, BLAS_Cpp::Transp P_rhs_trans) const
Abstract interface for objects that represent a space for mutable coordinate vectors.
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, MatrixSymOp *symwo_lhs) const
void syrk(const MatrixOp &mwo_rhs, BLAS_Cpp::Transp M_trans, value_type alpha, value_type beta, MatrixSymOp *sym_lhs)
Perform a rank-k update of a symmetric matrix of the form:
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)
mwo_lhs += alpha * op(M_rhs) * op(P_rhs).
bool Mp_StM(MatrixOp *mwo_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs) const
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)
mwo_lhs += alpha * op(P) * op(M_rhs).
void V_InvMtV(VectorMutable *v_lhs, const MatrixNonsing &M_rhs1, BLAS_Cpp::Transp trans_rhs1, const Vector &v_rhs2)
v_lhs = inv(op(M_rhs1)) * v_rhs2
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)
Transp trans_trans(Transp _trans1, Transp _trans2)
void M_StMtInvM(MatrixOp *m_lhs, value_type alpha, const MatrixOp &mwo_rhs1, BLAS_Cpp::Transp trans_rhs1, BLAS_Cpp::Transp trans_rhs2) const
Base class for all matrices that support basic matrix operations.
void M_StInvMtM(MatrixOp *m_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs1, const MatrixOp &mwo_rhs2, BLAS_Cpp::Transp trans_rhs2) const
const mwo_ptr_t & mwo() const
Abstract interface for mutable coordinate vectors {abstract}.
void M_StInvMtM(MatrixOp *m_lhs, value_type alpha, const MatrixNonsing &M_rhs1, BLAS_Cpp::Transp trans_rhs1, const MatrixOp &mwo_rhs2, BLAS_Cpp::Transp trans_rhs2)
m_lhs = alpha * inv(op(mwo_rhs1)) * op(mwo_rhs2) (right)
Transp
Aggregate matrix class pulling together a MatrixOp object and a MatrixNonsing object into a unified m...
bool Mp_StMtM(MatrixOp *mwo_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs1, const MatrixOp &mwo_rhs2, BLAS_Cpp::Transp trans_rhs2, value_type beta) const
Concrete matrix type to represent general permutation (mapping) matrices.