DenseLinAlgPack: Concreate C++ Classes for Dense Blas-Compatible Linear Algebra  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
DenseLinAlgPack_MatVecCompare.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 <limits>
43 #include <algorithm>
44 
45 #include "DenseLinAlgPack_MatVecCompare.hpp"
46 #include "DenseLinAlgPack_DMatrixClass.hpp"
47 #include "DenseLinAlgPack_DMatrixAsTriSym.hpp"
48 
49 namespace {
50 //
51 using DenseLinAlgPack::value_type;
52 //
53 using DenseLinAlgPack::sqrt_eps;
54 //
55 template< class T >
56 inline
57 T my_max( const T& v1, const T& v2 ) { return v1 > v2 ? v1 : v2; }
58 //
59 template< class T >
60 inline
61 T my_min( const T& v1, const T& v2 ) { return v1 < v2 ? v1 : v2; }
62 //
63 inline
64 bool _comp(value_type val1, value_type val2)
65 {
66  const value_type denom = my_max( ::fabs(val1), 1.0 ); // compare relative errors.
67  return ::fabs(val1 - val2) / denom < sqrt_eps;
68 }
69 
70 } // end namespace
71 
72 bool DenseLinAlgPack::comp(const DVectorSlice& vs1, const DVectorSlice& vs2) {
73  DVectorSlice::const_iterator
74  vs1_itr = vs1.begin(),
75  vs2_itr = vs2.begin();
76  for(; vs1_itr != vs1.end() && vs2_itr != vs2.end(); ++vs1_itr, ++vs2_itr)
77  if( !_comp(*vs1_itr,*vs2_itr) ) return false;
78  return true;
79 }
80 
81 bool DenseLinAlgPack::comp(const DVectorSlice& vs, value_type alpha) {
82  DVectorSlice::const_iterator vs_itr = vs.begin();
83  for(; vs_itr != vs.end(); ++vs_itr)
84  if( !_comp(*vs_itr,alpha) ) return false;
85  return true;
86 }
87 
88 bool DenseLinAlgPack::comp(const DMatrixSlice& gms1, BLAS_Cpp::Transp trans1
89  , const DMatrixSlice& gms2, BLAS_Cpp::Transp trans2)
90 {
91  for(size_type i = 1; i < my_min(gms1.cols(),gms2.cols()); ++i)
92  if( !comp( col(gms1,trans1,i) , col( gms2, trans2, i ) ) ) return false;
93  return true;
94 }
95 
96 bool DenseLinAlgPack::comp(const DMatrixSlice& gms, value_type alpha)
97 {
98  for(size_type i = 1; i < gms.cols(); ++i)
99  if( !comp( gms.col(i) , alpha ) ) return false;
100  return true;
101 }
102 
103 bool DenseLinAlgPack::comp(const DMatrixSliceTriEle& tri_gms1, const DMatrixSliceTriEle& tri_gms2)
104 {
107  trans1 = bool_to_trans(tri_gms1.uplo() != BLAS_Cpp::lower),
108  trans2 = bool_to_trans(tri_gms2.uplo() != BLAS_Cpp::lower);
109 
110  const size_type n = tri_gms1.rows(); // same as cols()
111  for(size_type i = 1; i < n; ++i) {
112  if( !comp( col(tri_gms1.gms(),trans1,i)(i,n), col(tri_gms2.gms(),trans2,i)(i,n) ) )
113  return false;
114  }
115  return true;
116 }
117 
118 bool DenseLinAlgPack::comp(const DMatrixSliceTriEle& tri_gms1, value_type alpha)
119 {
122  trans1 = bool_to_trans(tri_gms1.uplo() != BLAS_Cpp::lower);
123 
124  const size_type n = tri_gms1.rows(); // same as cols()
125  for(size_type i = 1; i < n; ++i) {
126  if( !comp( col(tri_gms1.gms(),trans1,i)(i,n), alpha ) )
127  return false;
128  }
129  return true;
130 }
131 
132 bool DenseLinAlgPack::comp_less(const DVectorSlice& vs, value_type alpha)
133 {
134  DVectorSlice::const_iterator vs_itr = vs.begin();
135  const value_type denom = my_max( ::fabs(alpha), 1.0 );
136  for(; vs_itr != vs.end(); ++vs_itr)
137  if( *vs_itr > alpha ) return false;
138  return true;
139 }
size_t size_type
Transp bool_to_trans(bool return_trans)
Transp