ConstrainedOptPack_MatrixKKTFullSpaceRelaxed.cpp

00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
00005 //                  Copyright (2003) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 // @HEADER
00028 //
00029 // ToDo: 6/2/00: Implement the relaxed version in the future.
00030 
00031 #include <assert.h>
00032 
00033 #include <sstream>
00034 #include <typeinfo>
00035 
00036 #include "ConstrainedOptPack_MatrixKKTFullSpaceRelaxed.hpp"
00037 #include "AbstractLinAlgPack/src/DirectSparseFortranCompatibleSolver.h"
00038 #include "AbstractLinAlgPack/src/MatrixConvertToSparseFortranCompatible.hpp"
00039 #include "AbstractLinAlgPack/test/TestMatrixConvertToSparseFortranCompatible.hpp"
00040 #include "DenseLinAlgPack_DVectorClass.hpp"
00041 #include "DenseLinAlgPack_AssertOp.hpp"
00042 
00043 namespace ConstrainedOptPack {
00044 
00045 MatrixKKTFullSpaceRelaxed::MatrixKKTFullSpaceRelaxed(
00046     const direct_solver_ptr_t& direct_solver    )
00047   :
00048       direct_solver_(direct_solver)
00049     , initialized_(false)
00050     , n_(0), m_(0)
00051     , G_(NULL), convG_(0), G_nz_(0)
00052     , A_(NULL), convA_(0), A_nz_(0)
00053 {}
00054 
00055 void MatrixKKTFullSpaceRelaxed::initialize(
00056     const MatrixOp& G, const MatrixOp& A
00057   , std::ostream* out, EPrintMoreOrLess print_what, ERunTests test_what )
00058 {
00059   // Set some members first
00060   out_      = out;
00061   print_what_   = print_what;
00062   test_what_    = test_what;
00063 
00064   // Validate that the matrices check out and get the conversion
00065   // interfaces.
00066   validate_and_set_matrices(G,A);
00067 
00068   use_relaxation_ = false;
00069 
00070   // Factor this matrix.
00071   //
00072   // If the structure of G and A looks to be the same we will reuse the
00073   // previous factorization structure if we have one.
00074 
00075   // ToDo: Finish this
00076 
00077   // Now we are initialized and ready to go.
00078   initialized_ = true;
00079 }
00080 
00081 void MatrixKKTFullSpaceRelaxed::initialize_relaxed(
00082     const MatrixOp& G, const MatrixOp& A
00083   , const DVectorSlice& c, value_type M
00084   , std::ostream* out, EPrintMoreOrLess print_what, ERunTests test_what )
00085 {
00086   // ToDo: implement the relaxation in the future!
00087   TEST_FOR_EXCEPT(true);
00088 }
00089 
00090 void MatrixKKTFullSpaceRelaxed::set_uninitialized()
00091 {
00092   G_ = NULL;
00093   A_ = NULL;
00094   initialized_ = false;
00095 }
00096 
00097 void MatrixKKTFullSpaceRelaxed::release_memory()
00098 {
00099   direct_solver().release_memory();
00100   n_ = m_ = 0;
00101   G_nz_ = A_nz_ = 0;
00102   set_uninitialized();
00103 }
00104 
00105 // Overridden from Matrix
00106 
00107 size_type MatrixKKTFullSpaceRelaxed::rows() const
00108 {
00109   assert_initialized();
00110   return n_ + m_ + (use_relaxation_ ? 1 : 0 );
00111 }
00112 
00113 size_type MatrixKKTFullSpaceRelaxed::cols() const
00114 {
00115   assert_initialized();
00116   return n_ + m_ + (use_relaxation_ ? 1 : 0 );
00117 }
00118 
00119 // Overridden from MatrixOp
00120 
00121 std::ostream& MatrixKKTFullSpaceRelaxed::output(std::ostream& out) const
00122 {
00123   assert_initialized();
00124   // ToDo: Finish me!
00125   TEST_FOR_EXCEPT(true);
00126   return out;
00127 }
00128 
00129 MatrixOp& MatrixKKTFullSpaceRelaxed::operator=(const MatrixOp& m)
00130 {
00131   assert_initialized();
00132   // ToDo: Finish me!
00133   TEST_FOR_EXCEPT(true);
00134   return *this;
00135 }
00136 
00137 void MatrixKKTFullSpaceRelaxed::Vp_StMtV(
00138     DVectorSlice* vs_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs1
00139   , const DVectorSlice& vs_rhs2, value_type beta) const
00140 {
00141   using AbstractLinAlgPack::Vp_StMtV;
00142 
00143   assert_initialized();
00144 
00145   DenseLinAlgPack::Vp_MtV_assert_sizes(vs_lhs->size(),rows(),cols(),BLAS_Cpp::no_trans
00146     ,vs_rhs2.size());
00147 
00148   if( use_relaxation_ ) {
00149     // ToDo: Implement the relaxation in the future
00150     TEST_FOR_EXCEPT(true);
00151   }
00152   else {
00153     // y = b*y + a*K*x
00154     //
00155     // [y1] = b * [y1] + a * [ G  A ] * [x1]    
00156     // [y2]       [y2]       [ A'   ]   [x2]
00157     //
00158     // y1 = b*y1 + a*G*x1 + a*A*x2
00159     //
00160     // y2 = b*y2 + a*A'*x1
00161 
00162     DVectorSlice
00163       y1 = (*vs_lhs)(1,n_),
00164       y2 = (*vs_lhs)(n_+1,n_+m_);
00165     const DVectorSlice
00166       x1 = vs_rhs2(1,n_),
00167       x2 = vs_rhs2(n_+1,n_+m_);
00168 
00169     // y1 = b*y1 + a*G*x1 + a*A*x2
00170 
00171     // y1 = a*G*x1 + b*y1
00172     Vp_StMtV( &y1, alpha, *G_, BLAS_Cpp::no_trans, x1, beta );
00173     // y1 += a*A*x2
00174     Vp_StMtV( &y1, alpha, *A_, BLAS_Cpp::no_trans, x2 );
00175 
00176     // y2 = a*A'*x1 + b*y2
00177     Vp_StMtV( &y2, alpha, *A_, BLAS_Cpp::trans, x1, beta );
00178   }
00179 }
00180 
00181 // Overridden from MatrixFactorized
00182 
00183 void MatrixKKTFullSpaceRelaxed::V_InvMtV( DVectorSlice* v_lhs, BLAS_Cpp::Transp trans_rhs1
00184   , const DVectorSlice& vs_rhs2) const
00185 {
00186   assert_initialized();
00187   // ToDo: Finish me!
00188   TEST_FOR_EXCEPT(true);
00189 }
00190 
00191 // Overridden from MatrixConvertToSparseFortranCompatible
00192 
00193 FortranTypes::f_int
00194 MatrixKKTFullSpaceRelaxed::num_nonzeros( EExtractRegion extract_region ) const
00195 {
00196   assert_matrices_set();
00197 
00198   FortranTypes::f_int
00199     nz;
00200   if( use_relaxation_ ) {
00201     // ToDo: Add support for the relaxation.
00202     TEST_FOR_EXCEPT(true);
00203   }
00204   else {
00205     // Return the number of nonzeros in a region of :
00206     //
00207     // K = [ G  A ]
00208     //     [ A'   ]
00209     typedef MatrixConvertToSparseFortranCompatible MCTSFC_t;
00210     const FortranTypes::f_int
00211       A_nz = convA_->num_nonzeros(MCTSFC_t::EXTRACT_FULL_MATRIX);
00212     nz = convG_->num_nonzeros( extract_region )
00213       + ( extract_region == MCTSFC_t::EXTRACT_FULL_MATRIX ? 2 * A_nz : A_nz );
00214   }
00215   return nz;
00216 }
00217 
00218 void MatrixKKTFullSpaceRelaxed::coor_extract_nonzeros(
00219     EExtractRegion extract_region
00220   , const FortranTypes::f_int len_Aval
00221     , FortranTypes::f_dbl_prec Aval[]
00222   , const FortranTypes::f_int len_Aij
00223     , FortranTypes::f_int Arow[]
00224     , FortranTypes::f_int Acol[]
00225     , const FortranTypes::f_int row_offset
00226     , const FortranTypes::f_int col_offset
00227    ) const
00228 {
00229   assert_matrices_set();
00230 
00231   if( len_Aval == 0 && len_Aij == 0 ) {
00232     if(*out_)
00233       *out_
00234         << "\n*** MatrixKKTFullSpaceRelaxed::coor_extract_nonzeros(...) : "
00235         << "Warning, nothing to compute: len_Aval==0 && len_Aij==0\n";
00236     return;
00237   }
00238 
00239   // Validate the conversion interfaces if asked to.
00240   if( test_what_ == RUN_TESTS ) {
00241 
00242     typedef MatrixConvertToSparseFortranCompatible
00243       MCTSFC_t;
00244     namespace slaptp = AbstractLinAlgPack::TestingPack;
00245     using slaptp::TestMatrixConvertToSparseFortranCompatible;
00246 
00247     const value_type
00248       warning_tol   = 1e-10,  // There should be very little roundoff error
00249       error_tol   = 1e-6;
00250 
00251     // Test G
00252     {
00253       slaptp::ETestSparseFortranFormat
00254         sparse_format = slaptp::TEST_COOR_FORMAT;
00255 
00256       if(out_)
00257         *out_
00258           << "\n*** Testing conversion interface for submatrix G ...\n";
00259 
00260       bool result = TestMatrixConvertToSparseFortranCompatible(
00261          extract_region,sparse_format,*convG_,*G_
00262         ,warning_tol,error_tol,true,out_
00263         ,print_what_==PRINT_MORE );
00264 
00265       if( !result) {
00266         char omsg[] = "MatrixKKTFullSpaceRelaxed : Error, the conversion "
00267           "interface for G did not check out\n";
00268         if(out_)
00269           *out_
00270             << std::endl << omsg;
00271         throw std::logic_error(omsg);
00272       }
00273       else {
00274         if(out_)
00275           *out_
00276             << "\n*** Conversion interface for G checked out!\n";
00277       }
00278     }
00279 
00280     // Test A
00281     {
00282       MCTSFC_t::EExtractRegion
00283         extract_region = MCTSFC_t::EXTRACT_FULL_MATRIX;
00284       slaptp::ETestSparseFortranFormat
00285         sparse_format = slaptp::TEST_COOR_FORMAT;
00286 
00287       if(out_)
00288         *out_
00289           << "\n*** Testing conversion interface for submatrix A ...\n";
00290 
00291       bool result = TestMatrixConvertToSparseFortranCompatible(
00292          extract_region,sparse_format,*convA_,*A_
00293         ,warning_tol,error_tol,true,out_
00294         ,print_what_==PRINT_MORE );
00295 
00296       if( !result) {
00297         char omsg[] = "MatrixKKTFullSpaceRelaxed : Error, the conversion "
00298           "interface for A did not check out\n";
00299         if(out_)
00300           *out_
00301             << std::endl << omsg;
00302         throw std::logic_error(omsg);
00303       }
00304       else {
00305         if(out_)
00306           *out_
00307             << "\n*** Conversion interface for A checked out!\n";
00308       }
00309     }
00310 
00311   }
00312 
00313   // Extract the nonzero entries
00314   if( use_relaxation_ ) {
00315     // ToDo: Add support for the relaxation.
00316     TEST_FOR_EXCEPT(true);
00317   }
00318   else {
00319     // Extract the nonzero entries in a region of :
00320     //
00321     // K = [ G  A ]
00322     //     [ A'   ]
00323 
00324     typedef MatrixConvertToSparseFortranCompatible MCTSFC_t;
00325 
00326     switch(extract_region) {
00327       case MCTSFC_t::EXTRACT_FULL_MATRIX:
00328         TEST_FOR_EXCEPT(true);  // We don't support this yet!
00329         break;
00330       case MCTSFC_t::EXTRACT_UPPER_TRIANGULAR:
00331         TEST_FOR_EXCEPT(true);  // We don't support this yet!
00332         break;
00333       case MCTSFC_t::EXTRACT_LOWER_TRIANGULAR:
00334       {
00335         // Set elements for
00336         //
00337         // K_lower = [ G_lower     ] n
00338         //           [ A'          ] m
00339         //             n        m
00340 
00341         // Get the numbers of nonzeros
00342         const FortranTypes::f_int
00343           G_lo_nz = convG_->num_nonzeros( MCTSFC_t::EXTRACT_LOWER_TRIANGULAR ),
00344           A_nz = convA_->num_nonzeros( MCTSFC_t::EXTRACT_FULL_MATRIX);
00345         assert( (len_Aval == 0 || len_Aval == G_lo_nz + A_nz)
00346             && (len_Aij == 0 || len_Aij == G_lo_nz + A_nz) );
00347         // Set the elements for G
00348         convG_->coor_extract_nonzeros(
00349            MCTSFC_t::EXTRACT_LOWER_TRIANGULAR
00350           , (len_Aval > 0 ? G_lo_nz : 0), Aval
00351           , (len_Aij > 0 ? G_lo_nz : 0), Arow, Acol, 0, 0 );
00352         // Set the elements for A'
00353         // To do this we will reverse the row and column indices
00354         // and the row offset will be n (col_offset in argument list)
00355         convA_->coor_extract_nonzeros(
00356            MCTSFC_t::EXTRACT_FULL_MATRIX
00357           , (len_Aval > 0 ? A_nz : 0), Aval+G_lo_nz
00358           , (len_Aij > 0 ? A_nz: 0), Acol+G_lo_nz, Arow+G_lo_nz, 0, n_ );
00359         break;
00360       }
00361       default:
00362         TEST_FOR_EXCEPT(true);
00363         break;
00364     }
00365   }
00366 }
00367 
00368 // Private member functions
00369 
00370 void MatrixKKTFullSpaceRelaxed::assert_initialized() const
00371 {
00372   if(!initialized_) {
00373     throw NotInitializedException("MatrixKKTFullSpaceRelaxed::assert_initialized() : "
00374       "Error, called a member function without initialize(...) or "
00375       "initialize_relaxed(...) being called first probably" );
00376   }
00377 }
00378 
00379 void MatrixKKTFullSpaceRelaxed::assert_matrices_set() const
00380 {
00381   if( !G_ || !convG_ || !A_ || !convA_ ) {
00382     throw NotInitializedException("MatrixKKTFullSpaceRelaxed::assert_matrices_set() : "
00383       "Error, the matrices G and A have not been set yet" );
00384   }
00385 }
00386 
00387 void MatrixKKTFullSpaceRelaxed::validate_and_set_matrices(
00388     const MatrixOp& G, const MatrixOp& A )
00389 {
00390   const size_type
00391     n = G.rows(),
00392     m = A.cols();
00393 
00394   if( G.cols() != n ) {
00395     std::ostringstream omsg;
00396     omsg
00397       << "MatrixKKTFullSpaceRelaxed::validate_and_set_matrices(...) : Error, "
00398       << "The matrix G with rows = " << n
00399       << " is not square with cols = " << G.cols();
00400     throw std::length_error(omsg.str());
00401   }
00402   if( A.rows() != n ) {
00403     std::ostringstream omsg;
00404     omsg
00405       << "MatrixKKTFullSpaceRelaxed::validate_and_set_matrices(...) : Error, "
00406       << "The number of rows in the matrix A with rows = " << A.rows()
00407       << ", cols = " << m
00408       << " does not match the size of G with rows = cols = " << n;
00409     throw std::length_error(omsg.str());
00410   }
00411   if( !(m < n) ||  n == 0 || m == 0 ) {
00412     std::ostringstream omsg;
00413     omsg
00414       << "MatrixKKTFullSpaceRelaxed::validate_and_set_matrices(...) : Error, "
00415       << "the size of the matrices G nxn and A nxm is not valid with "
00416       << "n = " << n << " and m = " << m;
00417     throw std::length_error(omsg.str());
00418   }
00419 
00420   const MatrixConvertToSparseFortranCompatible
00421     *convG = dynamic_cast<const MatrixConvertToSparseFortranCompatible*>(&G);
00422   if( ! convG ) {
00423     std::ostringstream omsg;
00424     omsg
00425       << "MatrixKKTFullSpaceRelaxed::validate_and_set_matrices(...) : Error, "
00426       << "The matrix G with concrete type " << typeid(G).name()
00427       << " does not support the MatrixConvertToSparseFortranCompatible "
00428       << "interface";
00429     throw InvalidMatrixType(omsg.str());
00430   }
00431   const MatrixConvertToSparseFortranCompatible
00432     *convA = dynamic_cast<const MatrixConvertToSparseFortranCompatible*>(&A);
00433   if( ! convA ) {
00434     std::ostringstream omsg;
00435     omsg
00436       << "MatrixKKTFullSpaceRelaxed::validate_and_set_matrices(...) : Error, "
00437       << "The matrix A with concrete type " << typeid(A).name()
00438       << " does not support the MatrixConvertToSparseFortranCompatible "
00439       << "interface";
00440     throw InvalidMatrixType(omsg.str());
00441   }
00442 
00443   // The input matrices checkout so set this stuff now
00444   G_ = &G;
00445   convG_ = convG;
00446   G_nz_ = G.nz();
00447   A_ = &A;
00448   convA_ = convA;
00449   A_nz_ = A.nz();
00450   n_ = n;
00451   m_ = m;
00452 }
00453 
00454 
00455 } // end namespace ConstrainedOptPack

Generated on Thu Sep 18 12:34:15 2008 for ConstrainedOptPack: C++ Tools for Constrained (and Unconstrained) Optimization by doxygen 1.3.9.1