MOOCHO (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AbstractLinAlgPack_MatrixSymNonsingSerial.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 #include <assert.h>
43 
53 #include "Teuchos_dyn_cast.hpp"
54 
55 namespace LinAlgOpPack {
58 }
59 
60 namespace AbstractLinAlgPack {
61 
64  , BLAS_Cpp::Transp B_trans, EMatrixDummyArg ) const
65 {
66  using BLAS_Cpp::trans;
67  using BLAS_Cpp::no_trans;
68  using BLAS_Cpp::trans_not;
74  //
75  // S = a * op(B) * inv(M) * op(B')
76  //
77  // We will form S won column at a time:
78  //
79  // S(:,j) = a * op(B) * inv(M) * op(B') * e(j)
80  //
81  // for j = 1 ... op(B').cols()
82  // t1 = op(B')*e(j)
83  // t2 = inv(M)*t1
84  // t3 = a*op(B)*t2
85  // S(:,j) = t3
86  //
87  // Above we only need to set the lower (lower triangle stored)
88  // or upper (upper triangle stored) part of S(:,k)
89  //
91  , B.rows(), B.cols(), trans_not(B_trans) );
93  , B.rows(), B.cols(), B_trans
94  , B.rows(), B.cols(), trans_not(B_trans) );
95 
96  DVector t1, t2, t3; // ToDo: Use temp workspace!
97  const size_type
98  opBT_cols = BLAS_Cpp::cols( B.cols(), B.rows(), B_trans ),
99  m = S->rows();
100  for( size_type j = 1; j <= m; ++j ) {
101  EtaVector e_j(j,opBT_cols); // e(j)
102  LinAlgOpPack::V_MtV( &t1, B, trans_not(B_trans), e_j() ); // t1 = op(B')*e(j)
103  AbstractLinAlgPack::V_InvMtV( &t2, *this, no_trans, t1() ); // t2 = inv(M)*t1
104  LinAlgOpPack::V_StMtV( &t3, a, B, B_trans, t2() ); // t3 = a*op(B)*t2
105  Range1D
106  rng = ( S->uplo() == BLAS_Cpp::upper ? Range1D(1,j) : Range1D(j,m) );
107  S->gms().col(j)(rng) = t3(rng);
108  }
109 }
110 
111 // Overridden from MatrixSymNonsing
112 
114  MatrixSymOp* symwo_lhs, value_type alpha
115  ,const MatrixOp& mwo, BLAS_Cpp::Transp mwo_trans
116  ,EMatrixDummyArg dummy
117  ) const
118 {
119  using Teuchos::dyn_cast;
120  this->M_StMtInvMtM(
121  &MatrixDenseSymMutableEncap(symwo_lhs)(), alpha
122  ,dyn_cast<const MatrixOpSerial>(mwo), mwo_trans
123  ,dummy );
124 }
125 
126 } // end namespace AbstractLinAlgPack
void M_StMtM(MatrixOp *M_lhs, value_type alpha, const MatrixOp &M_rhs1, BLAS_Cpp::Transp trans_rhs1, const MatrixOp &M_rhs2, BLAS_Cpp::Transp trans_rhs2)
M_lhs = alpha * op(M_rhs1) * op(M_rhs2).
void assign(DMatrix *gm_lhs, value_type alpha)
gm_lhs = alpha (elementwise)
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 MtM_assert_sizes(size_type m_rhs1_rows, size_type m_rhs1_cols, BLAS_Cpp::Transp trans_rhs1, size_type m_rhs2_rows, size_type m_rhs2_cols, BLAS_Cpp::Transp trans_rhs2)
op(m_lhs) += op(m_rhs1)
const DMatrixSliceTriEle tri_ele(const DMatrixSlice &gms, BLAS_Cpp::Uplo uplo)
void Mp_MtM_assert_sizes(size_type m_lhs_rows, size_type m_lhs_cols, BLAS_Cpp::Transp trans_lhs, size_type m_rhs1_rows, size_type m_rhs1_cols, BLAS_Cpp::Transp trans_rhs1, size_type m_rhs2_rows, size_type m_rhs2_cols, BLAS_Cpp::Transp trans_rhs2)
op(m_lhs) += op(m_rhs1) * op(m_rhs2)
void V_StMtV(VectorMutable *v_lhs, value_type alpha, const MatrixOp &M_rhs1, BLAS_Cpp::Transp trans_rhs1, const V &V_rhs2)
v_lhs = alpha * op(M_rhs1) * V_rhs2.
Transposed.
Base class for all matrices implemented in a shared memory address space.
Interface adding operations specific for a symmetric matrix {abstract}.
Not transposed.
virtual size_type cols() const
Return the number of columns in the matrix.
. One-based subregion index range class.
T_To & dyn_cast(T_From &from)
Helper class type that simplifies the usage of the MatrixSymOpGetGMSSymMutable interface for clients...
Create an eta vector (scaled by alpha = default 1).
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
f_dbl_prec f_dbl_prec f_dbl_prec * S
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)
DMatrixSliceTriEle nonconst_tri_ele(DMatrixSlice gms, BLAS_Cpp::Uplo uplo)
Return a triangular element-wise matrix.
Base class for all matrices that support basic matrix operations.
Transp trans_not(Transp _trans)
Return the opposite of the transpose argument.
void V_MtV(VectorMutable *v_lhs, const MatrixOp &M_rhs1, BLAS_Cpp::Transp trans_rhs1, const V &V_rhs2)
v_lhs = op(M_rhs1) * V_rhs2.
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)
virtual size_type rows() const
Return the number of rows in the matrix.
virtual void M_StMtInvMtM(DMatrixSliceSym *sym_gms_lhs, value_type alpha, const MatrixOpSerial &mwo, BLAS_Cpp::Transp mwo_trans, EMatrixDummyArg) const
sym_gms_lhs = alpha * op(mwo) * inv(M) * op(mwo)'.
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.
size_type cols(size_type rows, size_type cols, BLAS_Cpp::Transp _trans)
Return columns of a possible transposed matrix.