ConstrainedOptPack: C++ Tools for Constrained (and Unconstrained) Optimization  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
ConstrainedOptPack_MatrixGenBanded.cpp
1 #if 0
2 
3 // @HEADER
4 // ***********************************************************************
5 //
6 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
7 // Copyright (2003) Sandia Corporation
8 //
9 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
10 // license for use of this work by or on behalf of the U.S. Government.
11 //
12 // Redistribution and use in source and binary forms, with or without
13 // modification, are permitted provided that the following conditions are
14 // met:
15 //
16 // 1. Redistributions of source code must retain the above copyright
17 // notice, this list of conditions and the following disclaimer.
18 //
19 // 2. Redistributions in binary form must reproduce the above copyright
20 // notice, this list of conditions and the following disclaimer in the
21 // documentation and/or other materials provided with the distribution.
22 //
23 // 3. Neither the name of the Corporation nor the names of the
24 // contributors may be used to endorse or promote products derived from
25 // this software without specific prior written permission.
26 //
27 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
28 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
29 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
30 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
31 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
32 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
33 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
34 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
35 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
36 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
37 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
38 //
39 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov)
40 //
41 // ***********************************************************************
42 // @HEADER
43 
44 #include <assert.h>
45 #include <sstream>
46 
47 #include "ConstrainedOptPack_MatrixGenBanded.hpp"
48 #include "DenseLinAlgPack_AssertOp.hpp"
49 #include "DenseLinAlgPack_LinAlgOpPack.hpp"
50 #include "DenseLinAlgPack_BLAS_Cpp.hpp"
51 #include "MiWorkspacePack.h"
52 
53 namespace ConstrainedOptPack {
54 
56  size_type m
57  ,size_type n
58  ,size_type kl
59  ,size_type ku
60  ,DMatrixSlice *MB
61  ,const release_resource_ptr_t& MB_release_resource_ptr
62  )
63 {
64  initialize(m,n,kl,ku,MB,MB_release_resource_ptr);
65 }
66 
68  size_type m
69  ,size_type n
70  ,size_type kl
71  ,size_type ku
72  ,DMatrixSlice *MB
73  ,const release_resource_ptr_t& MB_release_resource_ptr
74  )
75 {
76  // Validate input
77 
78  if( m == 0 ) {
79  if( n != 0 )
80  throw std::invalid_argument(
81  "MatrixGenBanded::initialize(...): Error, "
82  "n must be 0 if m == 0" );
83  if( kl != 0 )
84  throw std::invalid_argument(
85  "MatrixGenBanded::initialize(...): Error, "
86  "kl must be 0 if m == 0" );
87  if( ku != 0 )
88  throw std::invalid_argument(
89  "MatrixGenBanded::initialize(...): Error, "
90  "ku must be 0 if m == 0" );
91  if( MB != NULL )
92  throw std::invalid_argument(
93  "MatrixGenBanded::initialize(...): Error, "
94  "MB must be NULL if m == 0" );
95  if( MB_release_resource_ptr.get() != NULL )
96  throw std::invalid_argument(
97  "MatrixGenBanded::initialize(...): Error, "
98  "MB_release_resource_ptr.get() must be NULL if m == 0" );
99  }
100  else {
101  if( kl + 1 > m )
102  throw std::invalid_argument(
103  "MatrixGenBanded::initialize(...): Error, "
104  "kl + 1 can not be larger than m" );
105  if( ku + 1 > n )
106  throw std::invalid_argument(
107  "MatrixGenBanded::initialize(...): Error, "
108  "ku + 1 can not be larger than n" );
109  if( MB == NULL )
110  throw std::invalid_argument(
111  "MatrixGenBanded::initialize(...): Error, "
112  "MB must not be NULL if n > 0" );
113  }
114 
115  // Set the members
116 
117  if( m == 0 ) {
118  m_ = 0;
119  n_ = 0;
120  kl_ = 0;
121  ku_ = 0;
122  MB_.bind(DMatrixSlice());
123  MB_release_resource_ptr_ = NULL;
124  }
125  else {
126  // Set the members
127  m_ = m;
128  n_ = n;
129  kl_ = kl;
130  ku_ = ku;
131  MB_.bind(*MB);
132  }
133 }
134 
135 // Overridden from MatrixOp
136 
138 {
139  return m_;
140 }
141 
143 {
144  return n_;
145 }
146 
148 {
149  return (ku_ + kl_ + 1) * n_ - ( (ku_+1) * (ku_+1) - (ku_+1) )/2 - ( (kl_+1) * (kl_+1) - (kl_+1) )/2; // Is correct?
150 }
151 
152 std::ostream& MatrixGenBanded::output(std::ostream& out) const
153 {
154  return MatrixOp::output(out); // ToDo: Implement specialized version later!
155 }
156 
158  DVectorSlice* y, value_type a, BLAS_Cpp::Transp M_trans
159  , const DVectorSlice& x, value_type b) const
160 {
161  assert_initialized();
162  DenseLinAlgPack::Vp_MtV_assert_sizes( y->size(), n_, n_, BLAS_Cpp::no_trans, x.size() );
163  BLAS_Cpp::gbmv(M_trans,m_,n_,kl_,ku_,a,MB_.col_ptr(1),MB_.max_rows(),x.raw_ptr(),x.stride()
164  ,b,y->raw_ptr(),y->stride());
165 }
166 
168  DVectorSlice* y, value_type a, BLAS_Cpp::Transp M_trans
169  , const SpVectorSlice& x, value_type b) const
170 {
171  assert_initialized();
172  MatrixOp::Vp_StMtV(y,a,M_trans,x,b); // ToDo: Implement spacialized operation when needed!
173 }
174 
176  DVectorSlice* y, value_type a
177  , const GenPermMatrixSlice& P, BLAS_Cpp::Transp P_trans
178  , BLAS_Cpp::Transp M_trans
179  , const DVectorSlice& x, value_type b) const
180 {
181  assert_initialized();
182  MatrixOp::Vp_StPtMtV(y,a,P,P_trans,M_trans,x,b); // ToDo: Implement spacialized operation when needed!
183 }
184 
186  DVectorSlice* y, value_type a
187  , const GenPermMatrixSlice& P, BLAS_Cpp::Transp P_trans
188  , BLAS_Cpp::Transp M_trans
189  , const SpVectorSlice& x, value_type b) const
190 {
191  assert_initialized();
192  MatrixOp::Vp_StPtMtV(y,a,P,P_trans,M_trans,x,b); // ToDo: Implement spacialized operation when needed!
193 }
194 
195 // Private member functions
196 
197 void MatrixGenBanded::assert_initialized() const
198 {
199  if( m_ == 0 )
200  throw std::logic_error("MatrixGenBanded::assert_initialized(): Error, "
201  "not initialized!" );
202 }
203 
204 } // end namespace ConstrainedOptPack
205 
206 #endif // 0
void Vp_StPtMtV(DVectorSlice *vs_lhs, value_type alpha, const GenPermMatrixSlice &P_rhs1, BLAS_Cpp::Transp P_rhs1_trans, BLAS_Cpp::Transp M_rhs2_trans, const DVectorSlice &vs_rhs3, value_type beta) const
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)
void initialize(size_type m=0, size_type n=0, size_type kl=0, size_type ku=0, DMatrixSlice *MB=NULL, const release_resource_ptr_t &MB_release_resource_ptr=NULL)
Initialize.
size_t size_type
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)
MatrixGenBanded(size_type m=0, size_type n=0, size_type kl=0, size_type ku=0, DMatrixSlice *MB=NULL, const release_resource_ptr_t &MB_release_resource_ptr=NULL)
Construct and Initialize.
void Vp_StMtV(DVectorSlice *vs_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs1, const DVectorSlice &vs_rhs2, value_type beta) const
std::ostream & output(std::ostream &out) const
Transp