MOOCHO (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AbstractLinAlgPack_LinAlgOpPackHack.cpp
Go to the documentation of this file.
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 
54 
55 
57  DMatrixSlice* gms_lhs, const MatrixOp& M_rhs,
58  BLAS_Cpp::Transp trans_rhs
59  )
60 {
62  gms_lhs->rows(), gms_lhs->cols(), BLAS_Cpp::no_trans,
63  M_rhs.rows(), M_rhs.cols(), trans_rhs );
64  (*gms_lhs) = 0.0;
65  Mp_StM(gms_lhs,1.0,M_rhs,trans_rhs);
66 }
67 
68 
71  ,const MatrixOp& B, BLAS_Cpp::Transp B_trans
72  )
73 {
78  const MatrixOpGetGMS
79  *B_get_gms = dynamic_cast<const MatrixOpGetGMS*>(&B);
80  if(B_get_gms) {
81  DenseLinAlgPack::Mp_StM( C, a, MatrixDenseEncap(*B_get_gms)(), B_trans );
82  }
83  else {
84  const size_type num_cols = C->cols();
86  B_mv = ( B_trans == BLAS_Cpp::no_trans
87  ? B.space_cols() : B.space_rows()
88  ).create_members(num_cols);
89  assign(B_mv.get(),B,B_trans);
90  for( size_type j = 1; j <= num_cols; ++j ) {
91  DenseLinAlgPack::Vp_StV(&C->col(j),a,VectorDenseEncap(*B_mv->col(j))());
92  }
93  }
94 }
95 
97  DVectorSlice* y, value_type a, const MatrixOp& M
98  ,BLAS_Cpp::Transp M_trans, const DVectorSlice& x, value_type b
99  )
100 {
101  using BLAS_Cpp::no_trans;
104  ay = ( M_trans == no_trans ? M.space_cols() : M.space_rows() ).create_member(),
105  ax = ( M_trans == no_trans ? M.space_rows() : M.space_cols() ).create_member();
106  (VectorDenseMutableEncap(*ay))() = *y;
107  (VectorDenseMutableEncap(*ax))() = x;
108  AbstractLinAlgPack::Vp_StMtV( ay.get(), a, M, M_trans, *ax, b );
109  *y = VectorDenseMutableEncap(*ay)();
110 }
111 
113  DVectorSlice* y, value_type a, const MatrixOp& M
114  ,BLAS_Cpp::Transp M_trans, const SpVectorSlice& x, value_type b
115  )
116 {
117  using BLAS_Cpp::no_trans;
120  ay = ( M_trans == no_trans ? M.space_cols() : M.space_rows() ).create_member();
121  (VectorDenseMutableEncap(*ay))() = *y;
122  AbstractLinAlgPack::Vp_StMtV( ay.get(), a, M, M_trans, x, b );
123  *y = VectorDenseMutableEncap(*ay)();
124 }
125 
127  DVectorSlice* y, const MatrixOpNonsing& M
128  ,BLAS_Cpp::Transp M_trans, const DVectorSlice& x
129  )
130 {
131  using BLAS_Cpp::trans;
134  ay = ( M_trans == trans ? M.space_cols() : M.space_rows() ).create_member(),
135  ax = ( M_trans == trans ? M.space_rows() : M.space_cols() ).create_member();
136  (VectorDenseMutableEncap(*ay))() = *y;
137  (VectorDenseMutableEncap(*ax))() = x;
138  AbstractLinAlgPack::V_InvMtV( ay.get(), M, M_trans, *ax );
139  *y = VectorDenseMutableEncap(*ay)();
140 }
141 
143  DVector* y, const MatrixOpNonsing& M
144  ,BLAS_Cpp::Transp M_trans, const DVectorSlice& x
145  )
146 {
147  using BLAS_Cpp::trans;
150  ay = ( M_trans == trans ? M.space_cols() : M.space_rows() ).create_member(),
151  ax = ( M_trans == trans ? M.space_rows() : M.space_cols() ).create_member();
152  (VectorDenseMutableEncap(*ax))() = x;
153  AbstractLinAlgPack::V_InvMtV( ay.get(), M, M_trans, *ax );
154  *y = VectorDenseMutableEncap(*ay)();
155 }
156 
158  DVectorSlice* y, const MatrixOpNonsing& M
159  ,BLAS_Cpp::Transp M_trans, const SpVectorSlice& x
160  )
161 {
162  using BLAS_Cpp::trans;
165  ay = ( M_trans == trans ? M.space_cols() : M.space_rows() ).create_member();
166  AbstractLinAlgPack::V_InvMtV( ay.get(), M, M_trans, x );
167  *y = VectorDenseMutableEncap(*ay)();
168 }
169 
171  DVector* y, const MatrixOpNonsing& M
172  ,BLAS_Cpp::Transp M_trans, const SpVectorSlice& x
173  )
174 {
175  y->resize(M.rows());
176  LinAlgOpPack::V_InvMtV( &(*y)(), M, M_trans, x );
177 }
178 
179 // These methods below are a real problem to implement in general.
180 //
181 // If the column space of op(M) could not return the above vector space
182 // for op(M).space_cols().space(P,P_trans) then we will try, as a last
183 // resort, to a dense serial vector and see what happens.
184 
187  ,const GenPermMatrixSlice& P, BLAS_Cpp::Transp P_trans
188  ,const MatrixOp& M, BLAS_Cpp::Transp M_trans
189  ,const DVectorSlice& x, value_type b
190  )
191 {
192  namespace mmp = MemMngPack;
193  using BLAS_Cpp::no_trans;
198  ay_space = ( M_trans == no_trans ? M.space_cols() : M.space_rows() ).space(P,P_trans);
200  ay = ( ay_space.get()
201  ? ay_space->create_member()
202  : Teuchos::rcp_implicit_cast<VectorMutable>(
203  Teuchos::rcp(new VectorMutableDense(BLAS_Cpp::rows(P.rows(),P.cols(),P_trans)))
204  ) ),
205  ax = ( M_trans == no_trans ? M.space_rows() : M.space_cols() ).create_member();
206  (VectorDenseMutableEncap(*ay))() = *y;
207  (VectorDenseMutableEncap(*ax))() = x;
208  Vp_StPtMtV( ay.get(), a, P, P_trans, M, M_trans, *ax, b );
209  *y = VectorDenseMutableEncap(*ay)();
210 }
211 
214  ,const GenPermMatrixSlice& P, BLAS_Cpp::Transp P_trans
215  ,const MatrixOp& M, BLAS_Cpp::Transp M_trans
216  ,const SpVectorSlice& x, value_type b
217  )
218 {
219  using BLAS_Cpp::no_trans;
224  ay = ( M_trans == no_trans ? M.space_cols() : M.space_rows() ).space(P,P_trans)->create_member();
225  (VectorDenseMutableEncap(*ay))() = *y;
226  Vp_StPtMtV( ay.get(), a, P, P_trans, M, M_trans, x, b );
227  *y = VectorDenseMutableEncap(*ay)();
228 }
AbstractLinAlgPack::size_type size_type
DVector "Adaptor" subclass for DenseLinAlgPack::DVectorSlice or DenseLinAlgPack::DVector objects...
void Vp_StMtV(DVectorSlice *vs_lhs, value_type alpha, const MatrixOp &mwo_rhs1, BLAS_Cpp::Transp trans_rhs1, const DVectorSlice &vs_rhs2, value_type beta=1.0)
vs_lhs = alpha * op(mwo_rhs1) * vs_rhs2 + beta * vs_lhs.
virtual const VectorSpace & space_rows() const =0
Vector space for vectors that are compatible with the rows of the matrix.
void Vp_StV(DVectorSlice *vs_lhs, value_type alpha, const DVectorSlice &vs_rhs)
vs_lhs += alpha * vs_rhs (BLAS xAXPY)
size_type cols() const
Return the number of columns.
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
size_type rows(size_type rows, size_type cols, BLAS_Cpp::Transp _trans)
Return rows of a possible transposed matrix.
void assign(VectorMutable *v_lhs, const V &V_rhs)
v_lhs = V_rhs.
T * get() const
void Mp_StM(DMatrixSlice *vs_lhs, value_type alpha, const MatrixOp &mwo_rhs1, BLAS_Cpp::Transp trans_rhs1)
m_lhs += alpha * op(mwo_rhs1).
void Mp_M_assert_sizes(size_type m_lhs_rows, size_type m_lhs_cols, BLAS_Cpp::Transp trans_lhs, size_type m_rhs_rows, size_type m_rhs_cols, BLAS_Cpp::Transp trans_rhs)
op(m_lhs) += op op(m_rhs)
Transposed.
void V_InvMtV(DVectorSlice *vs_lhs, const MatrixOpNonsing &mwo_rhs1, BLAS_Cpp::Transp trans_rhs1, const DVectorSlice &vs_rhs2)
vs_lhs = inv(op(mwo_rhs1)) * vs_rhs2.
Not transposed.
void Vp_StPtMtV(DVectorSlice *vs_lhs, value_type alpha, const GenPermMatrixSlice &gpms_rhs1, BLAS_Cpp::Transp trans_rhs1, const MatrixOp &mwo_rhs2, BLAS_Cpp::Transp trans_rhs2, const DVectorSlice &vs_rhs3, value_type beta=1.0)
vs_lhs = alpha * op(gpms_rhs1) * op(mwo_rhs2) * vs_rhs3 + beta * vs_lhs.
virtual size_type cols() const
Return the number of columns in the matrix.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Abstract interface for objects that represent a space for mutable coordinate vectors.
Extract a constant DenseLinAlgPack::DVectorSlice view of a Vector object.
const LAPACK_C_Decl::f_int const LAPACK_C_Decl::f_int const LAPACK_C_Decl::f_dbl_prec const LAPACK_C_Decl::f_int const LAPACK_C_Decl::f_int LAPACK_C_Decl::f_dbl_prec B[]
const LAPACK_C_Decl::f_int & M
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)
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.
Base class for all matrices that support basic matrix operations.
SparseVectorSlice< SparseElement< index_type, value_type > > SpVectorSlice
DenseLinAlgPack::VectorSliceTmpl< value_type > DVectorSlice
size_type rows() const
Return the number of rows.
DenseLinAlgPack::VectorTmpl< value_type > DVector
AbstractLinAlgPack::value_type value_type
virtual const VectorSpace & space_cols() const =0
Vector space for vectors that are compatible with the columns of the matrix.
Abstract interface for mutable coordinate vectors {abstract}.
void Mp_StM(DMatrixSliceTriEle *tri_lhs, value_type alpha, const DMatrixSliceTriEle &tri_rhs)
tri_lhs += alpha * tri_rhs (BLAS xAXPY)
Abstract base class for all nonsingular polymorphic matrices that can be used to compute matrix-vecto...
Extract a non-const DenseLinAlgPack::DVectorSlice view of a VectorMutable object. ...
virtual size_type rows() const
Return the number of rows in the matrix.
DVectorSlice col(size_type j)
Return DVectorSlice object representing the jth column (1-based; 1,2,..,#this->cols()#, or throw std::out_of_range)
Transp
TRANS.
Concrete matrix type to represent general permutation (mapping) matrices.