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