ConstrainedOptPack_MatrixVarReductImplicit.cpp

Go to the documentation of this file.
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 #include "ConstrainedOptPack_MatrixVarReductImplicit.hpp"
00030 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
00031 #include "AbstractLinAlgPack_MatrixOpNonsing.hpp"
00032 #include "AbstractLinAlgPack_MatrixOpOut.hpp"
00033 #include "AbstractLinAlgPack_GenPermMatrixSlice.hpp"
00034 #include "AbstractLinAlgPack_SpVectorClass.hpp"
00035 #include "AbstractLinAlgPack_AssertOp.hpp"
00036 #include "Teuchos_Workspace.hpp"
00037 #include "Teuchos_TestForException.hpp"
00038 
00039 namespace {
00040 
00041 //
00042 // Implicit matrix-vector multiplication:
00043 //
00044 // y = b*y + a*op(inv(C)*N)*x
00045 //
00046 template<class V>
00047 void imp_Vp_StMtV_implicit(
00048   AbstractLinAlgPack::VectorMutable               *y
00049   ,AbstractLinAlgPack::value_type                       a
00050   ,const AbstractLinAlgPack::MatrixOpNonsing    &C
00051   ,const AbstractLinAlgPack::MatrixOp               &N
00052   ,BLAS_Cpp::Transp                                     D_trans
00053   ,const V                                              &x
00054   ,DenseLinAlgPack::value_type                               b
00055   )
00056 {
00057   using BLAS_Cpp::no_trans;
00058   using BLAS_Cpp::trans;
00059   namespace alap = AbstractLinAlgPack;
00060   using Teuchos::Workspace;
00061   Teuchos::WorkspaceStore* wss = Teuchos::get_default_workspace_store().get();
00062 
00063   const DenseLinAlgPack::size_type
00064     r   = C.rows(),
00065     dof = N.cols();
00066 
00067   // ToDo: Pass in workspace vectors to save some allocation time!
00068 
00069   if( D_trans == no_trans ) {
00070     // y = b*y
00071     alap::Vt_S(y,b);
00072     //
00073     // y += -a * inv(C) * ( N * x )
00074     //
00075     alap::VectorSpace::vec_mut_ptr_t
00076       t1 = N.space_cols().create_member(),
00077       t2 = C.space_rows().create_member();
00078     // t1 = N*x
00079     LinAlgOpPack::V_MtV( t1.get(), N, no_trans, x );
00080     // t2 = inv(C) * t1
00081     alap::V_InvMtV( t2.get(), C, no_trans, *t1 );
00082     // y += a*t2
00083     alap::Vp_StV( y, -a, *t2 );
00084   }
00085   else {
00086     //
00087     // y = b*y - a * N' * ( inv(C') * x )
00088     //
00089     alap::VectorSpace::vec_mut_ptr_t
00090       t1 = C.space_cols().create_member();
00091     // t1 = inv(C')*x
00092     alap::V_InvMtV( t1.get(), C, trans, x );
00093     // y = b*y - a*N'*t1
00094       alap::Vp_StMtV( y, -a,  N, trans, *t1, b );
00095   }
00096 }
00097 
00098 /*
00099 
00100 //
00101 // Generate a row of inv(C)*N if not already computed.
00102 //
00103 void imp_compute_InvCtN_row(
00104   DenseLinAlgPack::size_type                                                 r
00105   ,const ConstrainedOptPack::DecompositionSystemVarReduct      &decomp_sys
00106   ,DenseLinAlgPack::size_type                                                j
00107   ,DenseLinAlgPack::DVectorSlice                                              *e_j  // Set to all zeros on input and output!
00108   ,ConstrainedOptPack::MatrixVarReductImplicit::InvCtN_rows_t  *InvCtN_rows
00109   )
00110 {
00111   typedef  DenseLinAlgPack::value_type  value_type;
00112   using DenseLinAlgPack::DVectorSlice;
00113   if( (*InvCtN_rows)[j-1] == NULL ) {
00114     // Generate row j of inv(C)*N
00115     value_type *vec = (*InvCtN_rows)[j-1] = new value_type[r]; // ToDo: We may want to allocate more vectors at once!
00116     DVectorSlice row_j(vec,r);
00117     // row_j = N'*inv(C')*e_j
00118     (*e_j)(j) = 1.0;
00119     imp_Vp_StMtV_implicit( &row_j, 1.0, decomp_sys, BLAS_Cpp::trans, *e_j, 0.0 );
00120     (*e_j)(j) = 0.0;
00121   }
00122 }
00123 
00124 //
00125 // Perform the matrix-vector multiplication:
00126 // 
00127 // y = b*y -a * op(P) * [inv(C) * N] * x
00128 //
00129 // by generating rows [inv(C)*N](j,:) for each nonzero entry op(P)(i,j).
00130 //
00131 template<class V>
00132 void imp_Vp_StPtMtV_by_row(
00133   DenseLinAlgPack::DVectorSlice                                               *y
00134   ,DenseLinAlgPack::value_type                                               a
00135   ,const ConstrainedOptPack::GenPermMatrixSlice                &P
00136   ,BLAS_Cpp::Transp                                                     P_trans
00137   ,const ConstrainedOptPack::DecompositionSystemVarReduct      &decomp_sys
00138   ,const V                                                              &x
00139   ,DenseLinAlgPack::value_type                                               b
00140   ,ConstrainedOptPack::MatrixVarReductImplicit::InvCtN_rows_t *InvCtN_rows
00141   )
00142 {
00143   using BLAS_Cpp::no_trans;
00144   using BLAS_Cpp::trans;
00145   using BLAS_Cpp::trans_not;
00146   using BLAS_Cpp::rows;
00147   using BLAS_Cpp::cols;
00148   using DenseLinAlgPack::dot;
00149   using DenseLinAlgPack::DVectorSlice;
00150   using AbstractLinAlgPack::dot;
00151   using AbstractLinAlgPack::Vp_StMtV;
00152   using AbstractLinAlgPack::GenPermMatrixSlice;
00153   using Teuchos::Workspace;
00154   Teuchos::WorkspaceStore* wss = Teuchos::get_default_workspace_store().get();
00155   
00156   const DenseLinAlgPack::size_type
00157     D_rows = decomp_sys.C().rows(),
00158     D_cols = decomp_sys.N().cols();
00159   // y = b*y
00160   if(b==0.0)       *y = 0.0;
00161   else if(b!=1.0)  DenseLinAlgPack::Vt_S(y,b);
00162   // Compute t = N'*inv(C')*e(j) then y(i) += -a*t'*x where op(P)(i,j) = 1.0
00163   Workspace<DenseLinAlgPack::value_type>   e_j_ws(wss,D_rows);
00164   DVectorSlice                              e_j(&e_j_ws[0],e_j_ws.size());
00165   e_j = 0.0;
00166   for( GenPermMatrixSlice::const_iterator itr = P.begin(); itr != P.end(); ++itr ) {
00167     const DenseLinAlgPack::size_type
00168       i = P_trans == no_trans ? itr->row_i() : itr->col_j(),
00169       j = P_trans == no_trans ? itr->col_j() : itr->row_i();
00170     // t = op(M') * e(j)
00171     imp_compute_InvCtN_row(D_rows,decomp_sys,j,&e_j,InvCtN_rows);
00172     DVectorSlice t((*InvCtN_rows)[j-1],D_cols);
00173     // y(i) += -a*t'*x
00174     (*y)(i) += (-a) * dot(t,x);
00175   }
00176 }
00177 
00178 */
00179 
00180 } // end namespace
00181 
00182 namespace ConstrainedOptPack {
00183 
00184 void MatrixVarReductImplicit::initialize(
00185   const mat_nonsing_ptr_t          &C
00186   ,const mat_ptr_t                 &N
00187   ,const mat_ptr_t                 &D_direct
00188   )
00189 {
00190   namespace rcp = MemMngPack;
00191   // Validate the inputs
00192   TEST_FOR_EXCEPTION(
00193     C.get() == NULL, std::invalid_argument
00194     ,"MatrixVarReductImplicit::initialize(...): Error, "
00195     "C.get() must not be NULL" );
00196   TEST_FOR_EXCEPTION(
00197     N.get() == NULL, std::invalid_argument
00198     ,"MatrixVarReductImplicit::initialize(...): Error, "
00199     "N.get() must not be NULL" );
00200   if( D_direct.get() ) {
00201     const bool is_compatible_cols = D_direct->space_cols().is_compatible(C->space_cols());
00202     TEST_FOR_EXCEPTION(
00203       !is_compatible_cols, VectorSpace::IncompatibleVectorSpaces
00204       ,"MatrixVarReductImplicit::initialize(...): Error, "
00205       "D_direct->space_cols() is not compatible with C->space_cols()" );
00206     const bool is_compatible_rows = D_direct->space_rows().is_compatible(N->space_rows());
00207     TEST_FOR_EXCEPTION(
00208       !is_compatible_rows, VectorSpace::IncompatibleVectorSpaces
00209       ,"MatrixVarReductImplicit::initialize(...): Error, "
00210       "D_direct->space_rows() is not compatible with N->space_rows()" );
00211   }
00212   // Set the members
00213   C_        = C;
00214   N_        = N;
00215   D_direct_ = D_direct;
00216   if(!InvCtN_rows_set_list_.empty()) { // Free previously allocated vectors
00217     for( InvCtN_rows_set_list_t::iterator itr = InvCtN_rows_set_list_.begin();
00218        itr != InvCtN_rows_set_list_.end(); ++itr )
00219         {
00220       InvCtN_rows_[*itr] = Teuchos::null;
00221     }
00222     InvCtN_rows_set_list_.clear();
00223   }
00224 }
00225 
00226 void MatrixVarReductImplicit::set_uninitialized()
00227 {
00228   namespace rcp = MemMngPack;
00229   C_        = Teuchos::null;
00230   N_        = Teuchos::null;
00231   D_direct_ = Teuchos::null;
00232 }
00233 
00234 // Overridden from MatrixBase
00235 
00236 size_type MatrixVarReductImplicit::rows() const
00237 {
00238   return C_.get() ? C_->rows() : 0;
00239 }
00240 
00241 size_type MatrixVarReductImplicit::cols() const
00242 {
00243   return N_.get() ? N_->cols() : 0;
00244 }
00245 
00246 // Overridden from MatrixOp
00247 
00248 const VectorSpace& MatrixVarReductImplicit::space_cols() const
00249 {
00250   assert_initialized();
00251   return C_->space_cols();
00252 }
00253 
00254 const VectorSpace& MatrixVarReductImplicit::space_rows() const
00255 {
00256   assert_initialized();
00257   return N_->space_rows();
00258 }
00259 
00260 MatrixOp& MatrixVarReductImplicit::operator=(const MatrixOp& M)
00261 {
00262   assert_initialized();
00263   TEST_FOR_EXCEPT(true); // ToDo: Finish!
00264   return *this;
00265 }
00266 
00267 std::ostream& MatrixVarReductImplicit::output(std::ostream& o) const
00268 {
00269   o << "This is a " << this->rows() << " x " << this->cols()
00270     << " variable reduction matrix D = -inv(C)*N where C and N are:\n"
00271     << "C =\n" << *C_
00272     << "N =\n" << *N_;
00273   return o;
00274 }
00275 
00276 void MatrixVarReductImplicit::Vp_StMtV(
00277   VectorMutable* y, value_type a
00278   ,BLAS_Cpp::Transp D_trans
00279   ,const Vector& x, value_type b
00280   ) const
00281 {
00282   assert_initialized();
00283   AbstractLinAlgPack::Vp_MtV_assert_compatibility(y,*this,D_trans,x);
00284   imp_Vp_StMtV_implicit( y, a, *C_, *N_, D_trans, x, b );
00285 }
00286 
00287 void MatrixVarReductImplicit::Vp_StMtV(
00288   VectorMutable* y, value_type a
00289   ,BLAS_Cpp::Transp D_trans
00290   ,const SpVectorSlice& x, value_type b
00291   ) const
00292 {
00293   using BLAS_Cpp::rows;
00294   using BLAS_Cpp::cols;
00295   using Teuchos::Workspace;
00296   Teuchos::WorkspaceStore* wss = Teuchos::get_default_workspace_store().get();
00297 
00298   assert_initialized();
00299   AbstractLinAlgPack::Vp_MtV_assert_compatibility(y,*this,D_trans,x);
00300   imp_Vp_StMtV_implicit( y, a, *C_, *N_, D_trans, x, b );
00301 /*
00302   const size_type
00303     D_rows = this->rows(), D_cols = this->cols(),
00304     opD_rows = rows(D_rows,D_cols,D_trans), opD_cols = cols(D_rows,D_cols,D_trans);
00305   DenseLinAlgPack::Vp_MtV_assert_sizes(y->size(),D_rows,D_cols,D_trans,x.size());
00306   if( use_dense_mat_vec_ && D_dense_.rows() > 0 ) {
00307     LinAlgOpPack::Vp_StMtV( y, a, D_dense_, D_trans, x, b );
00308   }
00309   else {
00310     if( x.nz() == x.size() ) {  // This is B.S.  Should use MatrixWithOpFactorized object for C!
00311       DVectorSlice dx = AbstractLinAlgPack::dense_view(x);
00312       imp_Vp_StMtV_implicit( y, -a, *decomp_sys_, D_trans, dx, b );
00313     }
00314     else if( D_trans == BLAS_Cpp::trans && x.nz() < D_cols ) {
00315       //
00316       // y = b*y + (-a)*[N'*inv(C')]*x
00317       //
00318       // We can do something crafty here.  We can generate columns of N'*inv(C')
00319       // and then perform y += -a*[N'*inv(C')](:,j)*x(j) for nonzero x(j)
00320       //
00321       Workspace<DenseLinAlgPack::value_type>   e_j_ws(wss,D_rows);
00322       DVectorSlice                              e_j(&e_j_ws[0],e_j_ws.size());
00323       e_j = 0.0;
00324       // y = b*y
00325       if(b==0.0)       *y = 0.0;
00326       else if(b!=1.0)  Vt_S(y,b);
00327       // y += -a*[N'*inv(C')](:,j)*x(j), for j <: { j | x(j) != 0.0 }
00328       const SpVectorSlice::difference_type o = x.offset();
00329       for( SpVectorSlice::const_iterator itr = x.begin(); itr != x.end(); ++itr ) {
00330         const size_type j = itr->indice() + o;
00331         imp_compute_InvCtN_row(D_rows,*decomp_sys_,j,&e_j,&InvCtN_rows_);
00332         DenseLinAlgPack::Vp_StV( y, -a * itr->value(), DVectorSlice(InvCtN_rows_[j-1],D_cols) );
00333       }
00334     }
00335     else {   // This is B.S.  Should use MatrixWithOpFactorized object for C!
00336       DVector dx;
00337       LinAlgOpPack::assign( &dx, x );
00338       imp_Vp_StMtV_implicit( y, -a, *decomp_sys_, D_trans, dx(), b );
00339     }
00340   }
00341 */
00342   // ToDo: Consider implementing the above specialized implementation!
00343 }
00344 
00345 void MatrixVarReductImplicit::Vp_StPtMtV(
00346   VectorMutable* y, value_type a
00347   ,const GenPermMatrixSlice& P, BLAS_Cpp::Transp P_trans
00348   ,BLAS_Cpp::Transp D_trans
00349   ,const Vector& x, value_type b
00350   ) const
00351 {
00352   using BLAS_Cpp::rows;
00353   using BLAS_Cpp::cols;
00354   assert_initialized();
00355 /*
00356   const size_type
00357     D_rows = this->rows(), D_cols = this->cols(),
00358     opD_rows = rows(D_rows,D_cols,D_trans), opD_cols = cols(D_rows,D_cols,D_trans);
00359   DenseLinAlgPack::Vp_MtV_assert_sizes(y->size(),P.rows(),P.cols(),P_trans,opD_rows);
00360   DenseLinAlgPack::Vp_MtV_assert_sizes(cols(P.rows(),P.cols(),P_trans),D_rows,D_cols,D_trans,x.size());
00361   if( D_dense_.rows() > 0 ) {
00362     AbstractLinAlgPack::dense_Vp_StPtMtV(y,a,P,P_trans,D_dense_,D_trans,x,b);
00363   }
00364   else if( P.nz() > D_cols || D_trans == BLAS_Cpp::trans ) {
00365     // Just use the default implementation
00366     MatrixOp::Vp_StPtMtV(y,a,P,P_trans,D_trans,x,b);
00367   }
00368   else {
00369     imp_Vp_StPtMtV_by_row(y,a,P,P_trans,*decomp_sys_,x,b,&InvCtN_rows_);
00370   }
00371 */
00372   MatrixOp::Vp_StPtMtV(y,a,P,P_trans,D_trans,x,b); // ToDo:Update specialized implementation above!
00373 }
00374 
00375 void MatrixVarReductImplicit::Vp_StPtMtV(
00376   VectorMutable* y, value_type a
00377   ,const GenPermMatrixSlice& P, BLAS_Cpp::Transp P_trans
00378   ,BLAS_Cpp::Transp D_trans
00379   ,const SpVectorSlice& x, value_type b
00380   ) const
00381 {
00382   using BLAS_Cpp::rows;
00383   using BLAS_Cpp::cols;
00384   assert_initialized();
00385 /*
00386   const size_type
00387     D_rows = this->rows(), D_cols = this->cols(),
00388     opD_rows = rows(D_rows,D_cols,D_trans), opD_cols = cols(D_rows,D_cols,D_trans);
00389   DenseLinAlgPack::Vp_MtV_assert_sizes(y->size(),P.rows(),P.cols(),P_trans,opD_rows);
00390   DenseLinAlgPack::Vp_MtV_assert_sizes(cols(P.rows(),P.cols(),P_trans),D_rows,D_cols,D_trans,x.size());
00391   if( D_dense_.rows() > 0 ) {
00392     AbstractLinAlgPack::dense_Vp_StPtMtV(y,a,P,P_trans,D_dense_,D_trans,x,b);
00393   }
00394   else if( P.nz() > D_cols || D_trans == BLAS_Cpp::trans ) {
00395     // Just use the default implementation
00396     MatrixOp::Vp_StPtMtV(y,a,P,P_trans,D_trans,x,b);
00397   }
00398   else {
00399     imp_Vp_StPtMtV_by_row(y,a,P,P_trans,*decomp_sys_,x,b,&InvCtN_rows_);
00400   }
00401 */
00402   MatrixOp::Vp_StPtMtV(y,a,P,P_trans,D_trans,x,b); // ToDo:Update specialized implementation above!
00403 }
00404 
00405 // Private member functions
00406 
00407 void MatrixVarReductImplicit::assert_initialized() const
00408 {
00409   TEST_FOR_EXCEPTION(
00410     C_.get() == NULL, std::logic_error
00411     ,"MatrixVarReductImplicit::assert_initialized(): Error, "
00412     "initialize(...) has not been called yet!" );
00413 }
00414 
00415 } // end namespace ConstrainedOptPack 

Generated on Thu Sep 18 12:35:14 2008 for MOOCHO (Single Doxygen Collection) by doxygen 1.3.9.1