ConstrainedOptPack_MatrixKKTFullSpaceRelaxed.cpp

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

Generated on Tue Jul 13 09:30:51 2010 for MOOCHO (Single Doxygen Collection) by  doxygen 1.4.7