ConstrainedOptPack: C++ Tools for Constrained (and Unconstrained) Optimization Version of the Day
ConstrainedOptPack_MatrixVarReductImplicit.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 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
00038 // 
00039 // ***********************************************************************
00040 // @HEADER
00041 
00042 #include "ConstrainedOptPack_MatrixVarReductImplicit.hpp"
00043 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
00044 #include "AbstractLinAlgPack_MatrixOpNonsing.hpp"
00045 #include "AbstractLinAlgPack_MatrixOpOut.hpp"
00046 #include "AbstractLinAlgPack_GenPermMatrixSlice.hpp"
00047 #include "AbstractLinAlgPack_SpVectorClass.hpp"
00048 #include "AbstractLinAlgPack_AssertOp.hpp"
00049 #include "Teuchos_Workspace.hpp"
00050 #include "Teuchos_Assert.hpp"
00051 
00052 namespace {
00053 
00054 //
00055 // Implicit matrix-vector multiplication:
00056 //
00057 // y = b*y + a*op(inv(C)*N)*x
00058 //
00059 template<class V>
00060 void imp_Vp_StMtV_implicit(
00061   AbstractLinAlgPack::VectorMutable               *y
00062   ,AbstractLinAlgPack::value_type                       a
00063   ,const AbstractLinAlgPack::MatrixOpNonsing    &C
00064   ,const AbstractLinAlgPack::MatrixOp               &N
00065   ,BLAS_Cpp::Transp                                     D_trans
00066   ,const V                                              &x
00067   ,DenseLinAlgPack::value_type                               b
00068   )
00069 {
00070   using BLAS_Cpp::no_trans;
00071   using BLAS_Cpp::trans;
00072   namespace alap = AbstractLinAlgPack;
00073   using Teuchos::Workspace;
00074   Teuchos::WorkspaceStore* wss = Teuchos::get_default_workspace_store().get();
00075 
00076   const DenseLinAlgPack::size_type
00077     r   = C.rows(),
00078     dof = N.cols();
00079 
00080   // ToDo: Pass in workspace vectors to save some allocation time!
00081 
00082   if( D_trans == no_trans ) {
00083     // y = b*y
00084     alap::Vt_S(y,b);
00085     //
00086     // y += -a * inv(C) * ( N * x )
00087     //
00088     alap::VectorSpace::vec_mut_ptr_t
00089       t1 = N.space_cols().create_member(),
00090       t2 = C.space_rows().create_member();
00091     // t1 = N*x
00092     LinAlgOpPack::V_MtV( t1.get(), N, no_trans, x );
00093     // t2 = inv(C) * t1
00094     alap::V_InvMtV( t2.get(), C, no_trans, *t1 );
00095     // y += a*t2
00096     alap::Vp_StV( y, -a, *t2 );
00097   }
00098   else {
00099     //
00100     // y = b*y - a * N' * ( inv(C') * x )
00101     //
00102     alap::VectorSpace::vec_mut_ptr_t
00103       t1 = C.space_cols().create_member();
00104     // t1 = inv(C')*x
00105     alap::V_InvMtV( t1.get(), C, trans, x );
00106     // y = b*y - a*N'*t1
00107       alap::Vp_StMtV( y, -a,  N, trans, *t1, b );
00108   }
00109 }
00110 
00111 /*
00112 
00113 //
00114 // Generate a row of inv(C)*N if not already computed.
00115 //
00116 void imp_compute_InvCtN_row(
00117   DenseLinAlgPack::size_type                                                 r
00118   ,const ConstrainedOptPack::DecompositionSystemVarReduct      &decomp_sys
00119   ,DenseLinAlgPack::size_type                                                j
00120   ,DenseLinAlgPack::DVectorSlice                                              *e_j  // Set to all zeros on input and output!
00121   ,ConstrainedOptPack::MatrixVarReductImplicit::InvCtN_rows_t  *InvCtN_rows
00122   )
00123 {
00124   typedef  DenseLinAlgPack::value_type  value_type;
00125   using DenseLinAlgPack::DVectorSlice;
00126   if( (*InvCtN_rows)[j-1] == NULL ) {
00127     // Generate row j of inv(C)*N
00128     value_type *vec = (*InvCtN_rows)[j-1] = new value_type[r]; // ToDo: We may want to allocate more vectors at once!
00129     DVectorSlice row_j(vec,r);
00130     // row_j = N'*inv(C')*e_j
00131     (*e_j)(j) = 1.0;
00132     imp_Vp_StMtV_implicit( &row_j, 1.0, decomp_sys, BLAS_Cpp::trans, *e_j, 0.0 );
00133     (*e_j)(j) = 0.0;
00134   }
00135 }
00136 
00137 //
00138 // Perform the matrix-vector multiplication:
00139 // 
00140 // y = b*y -a * op(P) * [inv(C) * N] * x
00141 //
00142 // by generating rows [inv(C)*N](j,:) for each nonzero entry op(P)(i,j).
00143 //
00144 template<class V>
00145 void imp_Vp_StPtMtV_by_row(
00146   DenseLinAlgPack::DVectorSlice                                               *y
00147   ,DenseLinAlgPack::value_type                                               a
00148   ,const ConstrainedOptPack::GenPermMatrixSlice                &P
00149   ,BLAS_Cpp::Transp                                                     P_trans
00150   ,const ConstrainedOptPack::DecompositionSystemVarReduct      &decomp_sys
00151   ,const V                                                              &x
00152   ,DenseLinAlgPack::value_type                                               b
00153   ,ConstrainedOptPack::MatrixVarReductImplicit::InvCtN_rows_t *InvCtN_rows
00154   )
00155 {
00156   using BLAS_Cpp::no_trans;
00157   using BLAS_Cpp::trans;
00158   using BLAS_Cpp::trans_not;
00159   using BLAS_Cpp::rows;
00160   using BLAS_Cpp::cols;
00161   using DenseLinAlgPack::dot;
00162   using DenseLinAlgPack::DVectorSlice;
00163   using AbstractLinAlgPack::dot;
00164   using AbstractLinAlgPack::Vp_StMtV;
00165   using AbstractLinAlgPack::GenPermMatrixSlice;
00166   using Teuchos::Workspace;
00167   Teuchos::WorkspaceStore* wss = Teuchos::get_default_workspace_store().get();
00168   
00169   const DenseLinAlgPack::size_type
00170     D_rows = decomp_sys.C().rows(),
00171     D_cols = decomp_sys.N().cols();
00172   // y = b*y
00173   if(b==0.0)       *y = 0.0;
00174   else if(b!=1.0)  DenseLinAlgPack::Vt_S(y,b);
00175   // Compute t = N'*inv(C')*e(j) then y(i) += -a*t'*x where op(P)(i,j) = 1.0
00176   Workspace<DenseLinAlgPack::value_type>   e_j_ws(wss,D_rows);
00177   DVectorSlice                              e_j(&e_j_ws[0],e_j_ws.size());
00178   e_j = 0.0;
00179   for( GenPermMatrixSlice::const_iterator itr = P.begin(); itr != P.end(); ++itr ) {
00180     const DenseLinAlgPack::size_type
00181       i = P_trans == no_trans ? itr->row_i() : itr->col_j(),
00182       j = P_trans == no_trans ? itr->col_j() : itr->row_i();
00183     // t = op(M') * e(j)
00184     imp_compute_InvCtN_row(D_rows,decomp_sys,j,&e_j,InvCtN_rows);
00185     DVectorSlice t((*InvCtN_rows)[j-1],D_cols);
00186     // y(i) += -a*t'*x
00187     (*y)(i) += (-a) * dot(t,x);
00188   }
00189 }
00190 
00191 */
00192 
00193 } // end namespace
00194 
00195 namespace ConstrainedOptPack {
00196 
00197 void MatrixVarReductImplicit::initialize(
00198   const mat_nonsing_ptr_t          &C
00199   ,const mat_ptr_t                 &N
00200   ,const mat_ptr_t                 &D_direct
00201   )
00202 {
00203   namespace rcp = MemMngPack;
00204   // Validate the inputs
00205   TEUCHOS_TEST_FOR_EXCEPTION(
00206     C.get() == NULL, std::invalid_argument
00207     ,"MatrixVarReductImplicit::initialize(...): Error, "
00208     "C.get() must not be NULL" );
00209   TEUCHOS_TEST_FOR_EXCEPTION(
00210     N.get() == NULL, std::invalid_argument
00211     ,"MatrixVarReductImplicit::initialize(...): Error, "
00212     "N.get() must not be NULL" );
00213   if( D_direct.get() ) {
00214     const bool is_compatible_cols = D_direct->space_cols().is_compatible(C->space_cols());
00215     TEUCHOS_TEST_FOR_EXCEPTION(
00216       !is_compatible_cols, VectorSpace::IncompatibleVectorSpaces
00217       ,"MatrixVarReductImplicit::initialize(...): Error, "
00218       "D_direct->space_cols() is not compatible with C->space_cols()" );
00219     const bool is_compatible_rows = D_direct->space_rows().is_compatible(N->space_rows());
00220     TEUCHOS_TEST_FOR_EXCEPTION(
00221       !is_compatible_rows, VectorSpace::IncompatibleVectorSpaces
00222       ,"MatrixVarReductImplicit::initialize(...): Error, "
00223       "D_direct->space_rows() is not compatible with N->space_rows()" );
00224   }
00225   // Set the members
00226   C_        = C;
00227   N_        = N;
00228   D_direct_ = D_direct;
00229   if(!InvCtN_rows_set_list_.empty()) { // Free previously allocated vectors
00230     for( InvCtN_rows_set_list_t::iterator itr = InvCtN_rows_set_list_.begin();
00231        itr != InvCtN_rows_set_list_.end(); ++itr )
00232         {
00233       InvCtN_rows_[*itr] = Teuchos::null;
00234     }
00235     InvCtN_rows_set_list_.clear();
00236   }
00237 }
00238 
00239 void MatrixVarReductImplicit::set_uninitialized()
00240 {
00241   namespace rcp = MemMngPack;
00242   C_        = Teuchos::null;
00243   N_        = Teuchos::null;
00244   D_direct_ = Teuchos::null;
00245 }
00246 
00247 // Overridden from MatrixBase
00248 
00249 size_type MatrixVarReductImplicit::rows() const
00250 {
00251   return C_.get() ? C_->rows() : 0;
00252 }
00253 
00254 size_type MatrixVarReductImplicit::cols() const
00255 {
00256   return N_.get() ? N_->cols() : 0;
00257 }
00258 
00259 // Overridden from MatrixOp
00260 
00261 const VectorSpace& MatrixVarReductImplicit::space_cols() const
00262 {
00263   assert_initialized();
00264   return C_->space_cols();
00265 }
00266 
00267 const VectorSpace& MatrixVarReductImplicit::space_rows() const
00268 {
00269   assert_initialized();
00270   return N_->space_rows();
00271 }
00272 
00273 MatrixOp& MatrixVarReductImplicit::operator=(const MatrixOp& M)
00274 {
00275   assert_initialized();
00276   TEUCHOS_TEST_FOR_EXCEPT(true); // ToDo: Finish!
00277   return *this;
00278 }
00279 
00280 std::ostream& MatrixVarReductImplicit::output(std::ostream& o) const
00281 {
00282   o << "This is a " << this->rows() << " x " << this->cols()
00283     << " variable reduction matrix D = -inv(C)*N where C and N are:\n"
00284     << "C =\n" << *C_
00285     << "N =\n" << *N_;
00286   return o;
00287 }
00288 
00289 void MatrixVarReductImplicit::Vp_StMtV(
00290   VectorMutable* y, value_type a
00291   ,BLAS_Cpp::Transp D_trans
00292   ,const Vector& x, value_type b
00293   ) const
00294 {
00295   assert_initialized();
00296   AbstractLinAlgPack::Vp_MtV_assert_compatibility(y,*this,D_trans,x);
00297   imp_Vp_StMtV_implicit( y, a, *C_, *N_, D_trans, x, b );
00298 }
00299 
00300 void MatrixVarReductImplicit::Vp_StMtV(
00301   VectorMutable* y, value_type a
00302   ,BLAS_Cpp::Transp D_trans
00303   ,const SpVectorSlice& x, value_type b
00304   ) const
00305 {
00306   using BLAS_Cpp::rows;
00307   using BLAS_Cpp::cols;
00308   using Teuchos::Workspace;
00309   Teuchos::WorkspaceStore* wss = Teuchos::get_default_workspace_store().get();
00310 
00311   assert_initialized();
00312   AbstractLinAlgPack::Vp_MtV_assert_compatibility(y,*this,D_trans,x);
00313   imp_Vp_StMtV_implicit( y, a, *C_, *N_, D_trans, x, b );
00314 /*
00315   const size_type
00316     D_rows = this->rows(), D_cols = this->cols(),
00317     opD_rows = rows(D_rows,D_cols,D_trans), opD_cols = cols(D_rows,D_cols,D_trans);
00318   DenseLinAlgPack::Vp_MtV_assert_sizes(y->size(),D_rows,D_cols,D_trans,x.size());
00319   if( use_dense_mat_vec_ && D_dense_.rows() > 0 ) {
00320     LinAlgOpPack::Vp_StMtV( y, a, D_dense_, D_trans, x, b );
00321   }
00322   else {
00323     if( x.nz() == x.size() ) {  // This is B.S.  Should use MatrixWithOpFactorized object for C!
00324       DVectorSlice dx = AbstractLinAlgPack::dense_view(x);
00325       imp_Vp_StMtV_implicit( y, -a, *decomp_sys_, D_trans, dx, b );
00326     }
00327     else if( D_trans == BLAS_Cpp::trans && x.nz() < D_cols ) {
00328       //
00329       // y = b*y + (-a)*[N'*inv(C')]*x
00330       //
00331       // We can do something crafty here.  We can generate columns of N'*inv(C')
00332       // and then perform y += -a*[N'*inv(C')](:,j)*x(j) for nonzero x(j)
00333       //
00334       Workspace<DenseLinAlgPack::value_type>   e_j_ws(wss,D_rows);
00335       DVectorSlice                              e_j(&e_j_ws[0],e_j_ws.size());
00336       e_j = 0.0;
00337       // y = b*y
00338       if(b==0.0)       *y = 0.0;
00339       else if(b!=1.0)  Vt_S(y,b);
00340       // y += -a*[N'*inv(C')](:,j)*x(j), for j <: { j | x(j) != 0.0 }
00341       const SpVectorSlice::difference_type o = x.offset();
00342       for( SpVectorSlice::const_iterator itr = x.begin(); itr != x.end(); ++itr ) {
00343         const size_type j = itr->indice() + o;
00344         imp_compute_InvCtN_row(D_rows,*decomp_sys_,j,&e_j,&InvCtN_rows_);
00345         DenseLinAlgPack::Vp_StV( y, -a * itr->value(), DVectorSlice(InvCtN_rows_[j-1],D_cols) );
00346       }
00347     }
00348     else {   // This is B.S.  Should use MatrixWithOpFactorized object for C!
00349       DVector dx;
00350       LinAlgOpPack::assign( &dx, x );
00351       imp_Vp_StMtV_implicit( y, -a, *decomp_sys_, D_trans, dx(), b );
00352     }
00353   }
00354 */
00355   // ToDo: Consider implementing the above specialized implementation!
00356 }
00357 
00358 void MatrixVarReductImplicit::Vp_StPtMtV(
00359   VectorMutable* y, value_type a
00360   ,const GenPermMatrixSlice& P, BLAS_Cpp::Transp P_trans
00361   ,BLAS_Cpp::Transp D_trans
00362   ,const Vector& x, value_type b
00363   ) const
00364 {
00365   using BLAS_Cpp::rows;
00366   using BLAS_Cpp::cols;
00367   assert_initialized();
00368 /*
00369   const size_type
00370     D_rows = this->rows(), D_cols = this->cols(),
00371     opD_rows = rows(D_rows,D_cols,D_trans), opD_cols = cols(D_rows,D_cols,D_trans);
00372   DenseLinAlgPack::Vp_MtV_assert_sizes(y->size(),P.rows(),P.cols(),P_trans,opD_rows);
00373   DenseLinAlgPack::Vp_MtV_assert_sizes(cols(P.rows(),P.cols(),P_trans),D_rows,D_cols,D_trans,x.size());
00374   if( D_dense_.rows() > 0 ) {
00375     AbstractLinAlgPack::dense_Vp_StPtMtV(y,a,P,P_trans,D_dense_,D_trans,x,b);
00376   }
00377   else if( P.nz() > D_cols || D_trans == BLAS_Cpp::trans ) {
00378     // Just use the default implementation
00379     MatrixOp::Vp_StPtMtV(y,a,P,P_trans,D_trans,x,b);
00380   }
00381   else {
00382     imp_Vp_StPtMtV_by_row(y,a,P,P_trans,*decomp_sys_,x,b,&InvCtN_rows_);
00383   }
00384 */
00385   MatrixOp::Vp_StPtMtV(y,a,P,P_trans,D_trans,x,b); // ToDo:Update specialized implementation above!
00386 }
00387 
00388 void MatrixVarReductImplicit::Vp_StPtMtV(
00389   VectorMutable* y, value_type a
00390   ,const GenPermMatrixSlice& P, BLAS_Cpp::Transp P_trans
00391   ,BLAS_Cpp::Transp D_trans
00392   ,const SpVectorSlice& x, value_type b
00393   ) const
00394 {
00395   using BLAS_Cpp::rows;
00396   using BLAS_Cpp::cols;
00397   assert_initialized();
00398 /*
00399   const size_type
00400     D_rows = this->rows(), D_cols = this->cols(),
00401     opD_rows = rows(D_rows,D_cols,D_trans), opD_cols = cols(D_rows,D_cols,D_trans);
00402   DenseLinAlgPack::Vp_MtV_assert_sizes(y->size(),P.rows(),P.cols(),P_trans,opD_rows);
00403   DenseLinAlgPack::Vp_MtV_assert_sizes(cols(P.rows(),P.cols(),P_trans),D_rows,D_cols,D_trans,x.size());
00404   if( D_dense_.rows() > 0 ) {
00405     AbstractLinAlgPack::dense_Vp_StPtMtV(y,a,P,P_trans,D_dense_,D_trans,x,b);
00406   }
00407   else if( P.nz() > D_cols || D_trans == BLAS_Cpp::trans ) {
00408     // Just use the default implementation
00409     MatrixOp::Vp_StPtMtV(y,a,P,P_trans,D_trans,x,b);
00410   }
00411   else {
00412     imp_Vp_StPtMtV_by_row(y,a,P,P_trans,*decomp_sys_,x,b,&InvCtN_rows_);
00413   }
00414 */
00415   MatrixOp::Vp_StPtMtV(y,a,P,P_trans,D_trans,x,b); // ToDo:Update specialized implementation above!
00416 }
00417 
00418 // Private member functions
00419 
00420 void MatrixVarReductImplicit::assert_initialized() const
00421 {
00422   TEUCHOS_TEST_FOR_EXCEPTION(
00423     C_.get() == NULL, std::logic_error
00424     ,"MatrixVarReductImplicit::assert_initialized(): Error, "
00425     "initialize(...) has not been called yet!" );
00426 }
00427 
00428 } // end namespace ConstrainedOptPack 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends