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_MatrixKKTFullSpaceRelaxed.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 // ToDo: 6/2/00: Implement the relaxed version in the future.
45 
46 #include <assert.h>
47 
48 #include <sstream>
49 #include <typeinfo>
50 
51 #include "ConstrainedOptPack_MatrixKKTFullSpaceRelaxed.hpp"
52 #include "AbstractLinAlgPack/src/DirectSparseFortranCompatibleSolver.h"
53 #include "AbstractLinAlgPack/src/MatrixConvertToSparseFortranCompatible.hpp"
54 #include "AbstractLinAlgPack/test/TestMatrixConvertToSparseFortranCompatible.hpp"
55 #include "DenseLinAlgPack_DVectorClass.hpp"
56 #include "DenseLinAlgPack_AssertOp.hpp"
57 
58 namespace ConstrainedOptPack {
59 
61  const direct_solver_ptr_t& direct_solver )
62  :
63  direct_solver_(direct_solver)
64  , initialized_(false)
65  , n_(0), m_(0)
66  , G_(NULL), convG_(0), G_nz_(0)
67  , A_(NULL), convA_(0), A_nz_(0)
68 {}
69 
71  const MatrixOp& G, const MatrixOp& A
72  , std::ostream* out, EPrintMoreOrLess print_what, ERunTests test_what )
73 {
74  // Set some members first
75  out_ = out;
76  print_what_ = print_what;
77  test_what_ = test_what;
78 
79  // Validate that the matrices check out and get the conversion
80  // interfaces.
81  validate_and_set_matrices(G,A);
82 
83  use_relaxation_ = false;
84 
85  // Factor this matrix.
86  //
87  // If the structure of G and A looks to be the same we will reuse the
88  // previous factorization structure if we have one.
89 
90  // ToDo: Finish this
91 
92  // Now we are initialized and ready to go.
93  initialized_ = true;
94 }
95 
97  const MatrixOp& G, const MatrixOp& A
98  , const DVectorSlice& c, value_type M
99  , std::ostream* out, EPrintMoreOrLess print_what, ERunTests test_what )
100 {
101  // ToDo: implement the relaxation in the future!
103 }
104 
106 {
107  G_ = NULL;
108  A_ = NULL;
109  initialized_ = false;
110 }
111 
113 {
114  direct_solver().release_memory();
115  n_ = m_ = 0;
116  G_nz_ = A_nz_ = 0;
118 }
119 
120 // Overridden from Matrix
121 
123 {
124  assert_initialized();
125  return n_ + m_ + (use_relaxation_ ? 1 : 0 );
126 }
127 
129 {
130  assert_initialized();
131  return n_ + m_ + (use_relaxation_ ? 1 : 0 );
132 }
133 
134 // Overridden from MatrixOp
135 
136 std::ostream& MatrixKKTFullSpaceRelaxed::output(std::ostream& out) const
137 {
138  assert_initialized();
139  // ToDo: Finish me!
141  return out;
142 }
143 
144 MatrixOp& MatrixKKTFullSpaceRelaxed::operator=(const MatrixOp& m)
145 {
146  assert_initialized();
147  // ToDo: Finish me!
149  return *this;
150 }
151 
153  DVectorSlice* vs_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs1
154  , const DVectorSlice& vs_rhs2, value_type beta) const
155 {
157 
158  assert_initialized();
159 
160  DenseLinAlgPack::Vp_MtV_assert_sizes(vs_lhs->size(),rows(),cols(),BLAS_Cpp::no_trans
161  ,vs_rhs2.size());
162 
163  if( use_relaxation_ ) {
164  // ToDo: Implement the relaxation in the future
166  }
167  else {
168  // y = b*y + a*K*x
169  //
170  // [y1] = b * [y1] + a * [ G A ] * [x1]
171  // [y2] [y2] [ A' ] [x2]
172  //
173  // y1 = b*y1 + a*G*x1 + a*A*x2
174  //
175  // y2 = b*y2 + a*A'*x1
176 
177  DVectorSlice
178  y1 = (*vs_lhs)(1,n_),
179  y2 = (*vs_lhs)(n_+1,n_+m_);
180  const DVectorSlice
181  x1 = vs_rhs2(1,n_),
182  x2 = vs_rhs2(n_+1,n_+m_);
183 
184  // y1 = b*y1 + a*G*x1 + a*A*x2
185 
186  // y1 = a*G*x1 + b*y1
187  Vp_StMtV( &y1, alpha, *G_, BLAS_Cpp::no_trans, x1, beta );
188  // y1 += a*A*x2
189  Vp_StMtV( &y1, alpha, *A_, BLAS_Cpp::no_trans, x2 );
190 
191  // y2 = a*A'*x1 + b*y2
192  Vp_StMtV( &y2, alpha, *A_, BLAS_Cpp::trans, x1, beta );
193  }
194 }
195 
196 // Overridden from MatrixFactorized
197 
198 void MatrixKKTFullSpaceRelaxed::V_InvMtV( DVectorSlice* v_lhs, BLAS_Cpp::Transp trans_rhs1
199  , const DVectorSlice& vs_rhs2) const
200 {
201  assert_initialized();
202  // ToDo: Finish me!
204 }
205 
206 // Overridden from MatrixConvertToSparseFortranCompatible
207 
208 FortranTypes::f_int
209 MatrixKKTFullSpaceRelaxed::num_nonzeros( EExtractRegion extract_region ) const
210 {
211  assert_matrices_set();
212 
213  FortranTypes::f_int
214  nz;
215  if( use_relaxation_ ) {
216  // ToDo: Add support for the relaxation.
218  }
219  else {
220  // Return the number of nonzeros in a region of :
221  //
222  // K = [ G A ]
223  // [ A' ]
224  typedef MatrixConvertToSparseFortranCompatible MCTSFC_t;
225  const FortranTypes::f_int
226  A_nz = convA_->num_nonzeros(MCTSFC_t::EXTRACT_FULL_MATRIX);
227  nz = convG_->num_nonzeros( extract_region )
228  + ( extract_region == MCTSFC_t::EXTRACT_FULL_MATRIX ? 2 * A_nz : A_nz );
229  }
230  return nz;
231 }
232 
234  EExtractRegion extract_region
235  , const FortranTypes::f_int len_Aval
236  , FortranTypes::f_dbl_prec Aval[]
237  , const FortranTypes::f_int len_Aij
238  , FortranTypes::f_int Arow[]
239  , FortranTypes::f_int Acol[]
240  , const FortranTypes::f_int row_offset
241  , const FortranTypes::f_int col_offset
242  ) const
243 {
244  assert_matrices_set();
245 
246  if( len_Aval == 0 && len_Aij == 0 ) {
247  if(*out_)
248  *out_
249  << "\n*** MatrixKKTFullSpaceRelaxed::coor_extract_nonzeros(...) : "
250  << "Warning, nothing to compute: len_Aval==0 && len_Aij==0\n";
251  return;
252  }
253 
254  // Validate the conversion interfaces if asked to.
255  if( test_what_ == RUN_TESTS ) {
256 
257  typedef MatrixConvertToSparseFortranCompatible
258  MCTSFC_t;
259  namespace slaptp = AbstractLinAlgPack::TestingPack;
260  using slaptp::TestMatrixConvertToSparseFortranCompatible;
261 
262  const value_type
263  warning_tol = 1e-10, // There should be very little roundoff error
264  error_tol = 1e-6;
265 
266  // Test G
267  {
268  slaptp::ETestSparseFortranFormat
269  sparse_format = slaptp::TEST_COOR_FORMAT;
270 
271  if(out_)
272  *out_
273  << "\n*** Testing conversion interface for submatrix G ...\n";
274 
275  bool result = TestMatrixConvertToSparseFortranCompatible(
276  extract_region,sparse_format,*convG_,*G_
277  ,warning_tol,error_tol,true,out_
278  ,print_what_==PRINT_MORE );
279 
280  if( !result) {
281  char omsg[] = "MatrixKKTFullSpaceRelaxed : Error, the conversion "
282  "interface for G did not check out\n";
283  if(out_)
284  *out_
285  << std::endl << omsg;
286  throw std::logic_error(omsg);
287  }
288  else {
289  if(out_)
290  *out_
291  << "\n*** Conversion interface for G checked out!\n";
292  }
293  }
294 
295  // Test A
296  {
297  MCTSFC_t::EExtractRegion
298  extract_region = MCTSFC_t::EXTRACT_FULL_MATRIX;
299  slaptp::ETestSparseFortranFormat
300  sparse_format = slaptp::TEST_COOR_FORMAT;
301 
302  if(out_)
303  *out_
304  << "\n*** Testing conversion interface for submatrix A ...\n";
305 
306  bool result = TestMatrixConvertToSparseFortranCompatible(
307  extract_region,sparse_format,*convA_,*A_
308  ,warning_tol,error_tol,true,out_
309  ,print_what_==PRINT_MORE );
310 
311  if( !result) {
312  char omsg[] = "MatrixKKTFullSpaceRelaxed : Error, the conversion "
313  "interface for A did not check out\n";
314  if(out_)
315  *out_
316  << std::endl << omsg;
317  throw std::logic_error(omsg);
318  }
319  else {
320  if(out_)
321  *out_
322  << "\n*** Conversion interface for A checked out!\n";
323  }
324  }
325 
326  }
327 
328  // Extract the nonzero entries
329  if( use_relaxation_ ) {
330  // ToDo: Add support for the relaxation.
332  }
333  else {
334  // Extract the nonzero entries in a region of :
335  //
336  // K = [ G A ]
337  // [ A' ]
338 
339  typedef MatrixConvertToSparseFortranCompatible MCTSFC_t;
340 
341  switch(extract_region) {
342  case MCTSFC_t::EXTRACT_FULL_MATRIX:
343  TEUCHOS_TEST_FOR_EXCEPT(true); // We don't support this yet!
344  break;
345  case MCTSFC_t::EXTRACT_UPPER_TRIANGULAR:
346  TEUCHOS_TEST_FOR_EXCEPT(true); // We don't support this yet!
347  break;
348  case MCTSFC_t::EXTRACT_LOWER_TRIANGULAR:
349  {
350  // Set elements for
351  //
352  // K_lower = [ G_lower ] n
353  // [ A' ] m
354  // n m
355 
356  // Get the numbers of nonzeros
357  const FortranTypes::f_int
358  G_lo_nz = convG_->num_nonzeros( MCTSFC_t::EXTRACT_LOWER_TRIANGULAR ),
359  A_nz = convA_->num_nonzeros( MCTSFC_t::EXTRACT_FULL_MATRIX);
360  TEUCHOS_TEST_FOR_EXCEPT( !( (len_Aval == 0 || len_Aval == G_lo_nz + A_nz ) )
361  && (len_Aij == 0 || len_Aij == G_lo_nz + A_nz) );
362  // Set the elements for G
363  convG_->coor_extract_nonzeros(
364  MCTSFC_t::EXTRACT_LOWER_TRIANGULAR
365  , (len_Aval > 0 ? G_lo_nz : 0), Aval
366  , (len_Aij > 0 ? G_lo_nz : 0), Arow, Acol, 0, 0 );
367  // Set the elements for A'
368  // To do this we will reverse the row and column indices
369  // and the row offset will be n (col_offset in argument list)
370  convA_->coor_extract_nonzeros(
371  MCTSFC_t::EXTRACT_FULL_MATRIX
372  , (len_Aval > 0 ? A_nz : 0), Aval+G_lo_nz
373  , (len_Aij > 0 ? A_nz: 0), Acol+G_lo_nz, Arow+G_lo_nz, 0, n_ );
374  break;
375  }
376  default:
378  break;
379  }
380  }
381 }
382 
383 // Private member functions
384 
385 void MatrixKKTFullSpaceRelaxed::assert_initialized() const
386 {
387  if(!initialized_) {
388  throw NotInitializedException("MatrixKKTFullSpaceRelaxed::assert_initialized() : "
389  "Error, called a member function without initialize(...) or "
390  "initialize_relaxed(...) being called first probably" );
391  }
392 }
393 
394 void MatrixKKTFullSpaceRelaxed::assert_matrices_set() const
395 {
396  if( !G_ || !convG_ || !A_ || !convA_ ) {
397  throw NotInitializedException("MatrixKKTFullSpaceRelaxed::assert_matrices_set() : "
398  "Error, the matrices G and A have not been set yet" );
399  }
400 }
401 
402 void MatrixKKTFullSpaceRelaxed::validate_and_set_matrices(
403  const MatrixOp& G, const MatrixOp& A )
404 {
405  const size_type
406  n = G.rows(),
407  m = A.cols();
408 
409  if( G.cols() != n ) {
410  std::ostringstream omsg;
411  omsg
412  << "MatrixKKTFullSpaceRelaxed::validate_and_set_matrices(...) : Error, "
413  << "The matrix G with rows = " << n
414  << " is not square with cols = " << G.cols();
415  throw std::length_error(omsg.str());
416  }
417  if( A.rows() != n ) {
418  std::ostringstream omsg;
419  omsg
420  << "MatrixKKTFullSpaceRelaxed::validate_and_set_matrices(...) : Error, "
421  << "The number of rows in the matrix A with rows = " << A.rows()
422  << ", cols = " << m
423  << " does not match the size of G with rows = cols = " << n;
424  throw std::length_error(omsg.str());
425  }
426  if( !(m < n) || n == 0 || m == 0 ) {
427  std::ostringstream omsg;
428  omsg
429  << "MatrixKKTFullSpaceRelaxed::validate_and_set_matrices(...) : Error, "
430  << "the size of the matrices G nxn and A nxm is not valid with "
431  << "n = " << n << " and m = " << m;
432  throw std::length_error(omsg.str());
433  }
434 
435  const MatrixConvertToSparseFortranCompatible
436  *convG = dynamic_cast<const MatrixConvertToSparseFortranCompatible*>(&G);
437  if( ! convG ) {
438  std::ostringstream omsg;
439  omsg
440  << "MatrixKKTFullSpaceRelaxed::validate_and_set_matrices(...) : Error, "
441  << "The matrix G with concrete type " << typeName(G)
442  << " does not support the MatrixConvertToSparseFortranCompatible "
443  << "interface";
444  throw InvalidMatrixType(omsg.str());
445  }
446  const MatrixConvertToSparseFortranCompatible
447  *convA = dynamic_cast<const MatrixConvertToSparseFortranCompatible*>(&A);
448  if( ! convA ) {
449  std::ostringstream omsg;
450  omsg
451  << "MatrixKKTFullSpaceRelaxed::validate_and_set_matrices(...) : Error, "
452  << "The matrix A with concrete type " << typeName(A)
453  << " does not support the MatrixConvertToSparseFortranCompatible "
454  << "interface";
455  throw InvalidMatrixType(omsg.str());
456  }
457 
458  // The input matrices checkout so set this stuff now
459  G_ = &G;
460  convG_ = convG;
461  G_nz_ = G.nz();
462  A_ = &A;
463  convA_ = convA;
464  A_nz_ = A.nz();
465  n_ = n;
466  m_ = m;
467 }
468 
469 
470 } // end namespace ConstrainedOptPack
471 
472 #endif // 0
void set_uninitialized()
Set the matrix to uninitialized.
void coor_extract_nonzeros(EExtractRegion extract_region, const FortranTypes::f_int len_Aval, FortranTypes::f_dbl_prec Aval[], const FortranTypes::f_int len_Aij, FortranTypes::f_int Arow[], FortranTypes::f_int Acol[], const FortranTypes::f_int row_offset, const FortranTypes::f_int col_offset) const
FortranTypes::f_int num_nonzeros(EExtractRegion extract_region) const
void Vp_StMtV(DVectorSlice *vs_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs1, const DVectorSlice &vs_rhs2, value_type beta) const
(2) vs_lhs = alpha * op(M_rhs1) * vs_rhs2 + beta * vs_lhs (BLAS xGEMV)
void release_memory()
Clear all allocated storage.
MatrixOp & operator=(const MatrixOp &m)
MatrixKKTFullSpaceRelaxed(const direct_solver_ptr_t &direct_solver=0)
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)
void V_InvMtV(DVectorSlice *v_lhs, BLAS_Cpp::Transp trans_rhs1, const DVectorSlice &vs_rhs2) const
(1) v_lhs = inv(op(M_rhs1)) * vs_rhs2
void initialize(const MatrixOp &G, const MatrixOp &A, std::ostream *out=0, EPrintMoreOrLess print_what=PRINT_LESS, ERunTests test_what=NO_TESTS)
Initialize the nonrelaxed matrix.
void initialize_relaxed(const MatrixOp &G, const MatrixOp &A, const DVectorSlice &c, value_type bigM=1e+10, std::ostream *out=0, EPrintMoreOrLess print_what=PRINT_LESS, ERunTests test_what=NO_TESTS)
Initialize the relaxed matrix.
std::ostream & output(std::ostream &out) const
Transp
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
std::string typeName(const T &t)