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_MatrixNonsing.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: 3/6/00: Provide default implementations for these
43 // operations.
44 
45 #include <assert.h>
46 
47 #include "AbstractLinAlgPack_MatrixNonsing.hpp"
48 #include "AbstractLinAlgPack_MultiVectorMutable.hpp"
49 #include "AbstractLinAlgPack_VectorSpace.hpp"
50 #include "AbstractLinAlgPack_SpVectorClass.hpp"
51 #include "AbstractLinAlgPack_SpVectorView.hpp"
52 #include "AbstractLinAlgPack_EtaVector.hpp"
53 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
54 #include "Teuchos_Assert.hpp"
55 #include "Teuchos_dyn_cast.hpp"
56 
57 namespace AbstractLinAlgPack {
58 
59 // Clone
60 
61 MatrixNonsing::mat_mns_mut_ptr_t
63 {
64  return Teuchos::null;
65 }
66 
67 MatrixNonsing::mat_mns_ptr_t
69 {
70  return const_cast<MatrixNonsing*>(this)->clone_mns(); // Implicit conversion to const
71 }
72 
73 // Level-2 BLAS
74 
76  VectorMutable* y, BLAS_Cpp::Transp M_trans, const SpVectorSlice& sx
77  ) const
78 {
79  if( sx.nz() ) {
81  x = (M_trans == BLAS_Cpp::no_trans
82  ? this->space_cols()
83  : this->space_rows()
84  ).create_member();
85  x->set_sub_vector(sub_vec_view(sx));
86  this->V_InvMtV(y,M_trans,*x);
87  }
88  else {
89  *y = 0.0;
90  }
91 }
92 
94  const Vector& v_rhs1, BLAS_Cpp::Transp trans_rhs2, const Vector& v_rhs3
95  ) const
96 {
98  v = (trans_rhs2 == BLAS_Cpp::no_trans
99  ? this->space_rows()
100  : this->space_cols()
101  ).create_member();
102  this->V_InvMtV( v.get(), trans_rhs2, v_rhs3 );
103  return dot(v_rhs1,*v);
104 }
105 
107  const SpVectorSlice& sv_rhs1, BLAS_Cpp::Transp trans_rhs2, const SpVectorSlice& sv_rhs3
108  ) const
109 {
111  v = (trans_rhs2 == BLAS_Cpp::no_trans
112  ? this->space_rows()
113  : this->space_cols()
114  ).create_member();
115  this->V_InvMtV( v.get(), trans_rhs2, sv_rhs3 );
116  return dot(sv_rhs1,*v);
117 }
118 
119 // Level-3 BLAS
120 
122  MatrixOp* C_lhs, value_type alpha
123  ,BLAS_Cpp::Transp M_trans
124  ,const MatrixOp& B, BLAS_Cpp::Transp B_trans
125  ) const
126 {
127  //
128  // C = a * inv(op(M)) * op(B)
129  //
130  using Teuchos::dyn_cast;
131  using BLAS_Cpp::no_trans;
132  using BLAS_Cpp::trans;
133 #ifdef TEUCHOS_DEBUG
135  C_lhs == NULL, std::invalid_argument
136  ,"MatrixNonsing::M_StInvMtM(...) : Error!" );
137 
138 #endif
139  const size_type
140  C_rows = C_lhs->rows(),
141  C_cols = C_lhs->cols();
142  const size_type
143  op_B_cols = BLAS_Cpp::cols( B.rows(), B.cols(), B_trans );
144 #ifdef TEUCHOS_DEBUG
145  // We can't check vector spaces since *this may not support MatrixOp
146  // However, we could dynamic cast to see if MatrixOp is supported and then
147  // be able to use Mp_MtM_assert_compatibility() but this is okay for now.
148  const size_type
149  M_rows = this->rows(),
150  M_cols = this->cols(),
151  op_B_rows = BLAS_Cpp::rows( B.rows(), B.cols(), B_trans );
153  C_rows != M_rows || M_rows != M_cols || M_cols != op_B_rows || C_cols != op_B_cols
154  , std::invalid_argument
155  ,"MatrixNonsing::M_StInvMtM(...) : Error!" );
156 #endif
157  //
158  // Compute C = a * inv(op(M)) * op(B) one column at a time:
159  //
160  // C(:,j) = inv(op(M)) * a * op(B) * e(j) , for j = 1...C.cols()
161  // \______________/
162  // t_j
163  //
166  t_j = ( B_trans == no_trans ? B.space_cols() : B.space_rows() ).create_member();
167  for( size_type j = 1; j <= C_cols; ++j ) {
168  // t_j = alpha * op(B) * e_j
169  EtaVector e_j( j, op_B_cols );
170  LinAlgOpPack::V_StMtV( t_j.get(), alpha, B, B_trans, e_j() );
171  // C(:,j) = inv(op(M)) * t_j
172  AbstractLinAlgPack::V_InvMtV( C.col(j).get(), *this, M_trans, *t_j );
173  }
174 }
175 
177  MatrixOp* g_lhs, value_type alpha
178  ,const MatrixOp& mwo_rhs1, BLAS_Cpp::Transp trans_rhs1
179  ,BLAS_Cpp::Transp trans_rhs2
180  ) const
181 {
182  TEUCHOS_TEST_FOR_EXCEPT(true); // ToDo: Implement!
183 }
184 
185 } // end namespace AbstractLinAlgPack
virtual mat_mns_mut_ptr_t clone_mns()
Clone the non-const matrix object (if supported).
virtual const VectorSpace & space_rows() const =0
Vector space for vectors that are compatible with the rows of the matrix.
Abstract interface for immutable, finite dimensional, coordinate vectors {abstract}.
friend void V_InvMtV(VectorMutable *v_lhs, const MatrixNonsing &M_rhs1, BLAS_Cpp::Transp trans_rhs1, const Vector &v_rhs2)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
size_type rows(size_type rows, size_type cols, BLAS_Cpp::Transp _trans)
T_To & dyn_cast(T_From &from)
T * get() const
virtual size_type cols() const
Return the number of columns in the matrix.
friend 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)
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
Base class for all matrices that support basic matrix operations.
friend value_type transVtInvMtV(const Vector &v_rhs1, const MatrixNonsing &M_rhs2, BLAS_Cpp::Transp trans_rhs2, const Vector &v_rhs3)
value_type dot(const Vector &v_rhs1, const Vector &v_rhs2)
result = v_rhs1' * v_rhs2
Interface for a collection of mutable vectors (multi-vector, matrix).
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}.
virtual size_type rows() const
Return the number of rows in the matrix.
friend 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)
virtual vec_mut_ptr_t col(index_type j)=0
Get a mutable column vector.
Abstract base class for all nonsingular polymorphic matrices that can solve for linear system with bu...
Transp
size_type cols(size_type rows, size_type cols, BLAS_Cpp::Transp _trans)
size_type nz() const
Return the number of non-zero elements.
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)