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_MatrixOpNonsing.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 <math.h>
43 
44 #include "AbstractLinAlgPack_MatrixOpNonsing.hpp"
45 #include "AbstractLinAlgPack_VectorSpace.hpp"
46 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
47 #include "Teuchos_Assert.hpp"
48 
49 namespace AbstractLinAlgPack {
50 
51 MatrixOpNonsing::mat_mwons_mut_ptr_t
53 {
54  return Teuchos::null;
55 }
56 
57 MatrixOpNonsing::mat_mwons_ptr_t
59 {
60  return Teuchos::null;
61 }
62 
65  EMatNormType requested_norm_type
66  ,bool allow_replacement
67  ) const
68 {
69  using BLAS_Cpp::no_trans;
70  using BLAS_Cpp::trans;
72  const VectorSpace
73  &space_cols = this->space_cols(),
74  &space_rows = this->space_rows();
75  const index_type
76  num_cols = space_rows.dim();
78  !(requested_norm_type == MAT_NORM_1 || requested_norm_type == MAT_NORM_INF), MethodNotImplemented
79  ,"MatrixOp::calc_norm(...): Error, This default implemenation can only "
80  "compute the one norm or the infinity norm!"
81  );
82  //
83  // Here we implement Algorithm 2.5 in "Applied Numerical Linear Algebra", Demmel (1997)
84  // using the momenclature in the text. This is applied to the inverse matrix.
85  //
86  const MatrixOpNonsing
87  &B = *this;
88  bool
89  do_trans = requested_norm_type == MAT_NORM_INF;
91  x = (do_trans ? space_rows : space_cols).create_member(1.0/num_cols),
92  w = (do_trans ? space_cols : space_rows).create_member(),
93  zeta = (do_trans ? space_cols : space_rows).create_member(),
94  z = (do_trans ? space_rows : space_cols).create_member();
95  const index_type max_iter = 5; // Recommended by Highman 1988, (see Demmel's reference)
96  value_type w_nrm = 0.0;
97  for( index_type k = 0; k <= max_iter; ++k ) {
98  V_InvMtV( w.get(), B, !do_trans ? no_trans : trans, *x ); // w = B*x
99  sign( *w, zeta.get() ); // zeta = sign(w)
100  V_InvMtV( z.get(), B, !do_trans ? trans : no_trans, *zeta ); // z = B'*zeta
101  value_type z_j = 0.0; // max |z(j)| = ||z||inf
102  index_type j = 0;
103  max_abs_ele( *z, &z_j, &j );
104  const value_type zTx = dot(*z,*x); // z'*x
105  w_nrm = w->norm_1(); // ||w||1
106  if( ::fabs(z_j) <= zTx ) { // Update
107  break;
108  }
109  else {
110  *x = 0.0;
111  x->set_ele(j,1.0);
112  }
113  }
114  const MatNorm M_nrm = this->calc_norm(requested_norm_type);
115  return MatNorm( w_nrm * M_nrm.value ,requested_norm_type );
116 }
117 
118 // Overridden from MatrixOp
119 
120 MatrixOpNonsing::mat_mut_ptr_t
122 {
123  return clone_mwons();
124 }
125 
126 MatrixOpNonsing::mat_ptr_t
128 {
129  return clone_mwons();
130 }
131 
132 // Overridden from MatrixNonsing
133 
134 MatrixOpNonsing::mat_mns_mut_ptr_t
136 {
137  return clone_mwons();
138 }
139 
140 MatrixOpNonsing::mat_mns_ptr_t
142 {
143  return clone_mwons();
144 }
145 
146 } // end namespace AbstractLinAlgPack
virtual const VectorSpace & space_rows() const =0
Vector space for vectors that are compatible with the rows of the matrix.
mat_mns_mut_ptr_t clone_mns()
Returns this->clone_mwons().
void sign(const Vector &v, VectorMutable *z)
Compute the sign of each element in an input vector.
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)
The induced infinity norm ||M||inf, i.e. max abs row sum.
Abstract interface for objects that represent a space for mutable coordinate vectors.
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
const MatNorm calc_norm(EMatNormType requested_norm_type=MAT_NORM_1, bool allow_replacement=false) const
Compute a norm of this matrix.
value_type dot(const Vector &v_rhs1, const Vector &v_rhs2)
result = v_rhs1' * v_rhs2
virtual index_type dim() const =0
Return the dimmension of the vector space.
The induced one norm ||M||1, i.e. max abs col sum.
virtual mat_mwons_mut_ptr_t clone_mwons()
Clone the non-const matrix object (if supported).
const MatNorm calc_cond_num(EMatNormType requested_norm_type=MAT_NORM_1, bool allow_replacement=false) const
Compute an estimate of the condition number of this matrix.
virtual const VectorSpace & space_cols() const =0
Vector space for vectors that are compatible with the columns of the matrix.
void max_abs_ele(const Vector &v, value_type *max_v_j, index_type *max_j)
Compute the maximum element in a vector.
Abstract base class for all nonsingular polymorphic matrices that can be used to compute matrix-vecto...
mat_mut_ptr_t clone()
Returns this->clone_mwons().