MOOCHO (Single Doxygen Collection) Version of the Day
AbstractLinAlgPack_MultiVectorMutableCols.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 // 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 "AbstractLinAlgPack_MultiVectorMutableCols.hpp"
00043 #include "AbstractLinAlgPack_MatrixSymOpGetGMSSymMutable.hpp"
00044 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
00045 #include "AbstractLinAlgPack_SpVectorClass.hpp"
00046 #include "Teuchos_dyn_cast.hpp"
00047 #include "Teuchos_Assert.hpp"
00048 
00049 namespace AbstractLinAlgPack {
00050 
00051 // Constructors/Initializers
00052 
00053 MultiVectorMutableCols::MultiVectorMutableCols()
00054 {}
00055 
00056 MultiVectorMutableCols::MultiVectorMutableCols(
00057   const  Teuchos::RCP<const VectorSpace>   &space_cols
00058   ,const  Teuchos::RCP<const VectorSpace>  &space_rows
00059   ,Teuchos::RCP<VectorMutable>       col_vecs[]
00060   )
00061 {
00062   this->initialize(space_cols,space_rows,col_vecs);
00063 }
00064   
00065 void MultiVectorMutableCols::initialize(
00066   const  Teuchos::RCP<const VectorSpace>   &space_cols
00067   ,const  Teuchos::RCP<const VectorSpace>  &space_rows
00068   ,Teuchos::RCP<VectorMutable>       col_vecs[]
00069   )
00070 {
00071   space_cols_ = space_cols;
00072   space_rows_ = space_rows;
00073   const size_type num_cols = space_rows->dim();
00074   col_vecs_.resize(num_cols);
00075   if(col_vecs) {
00076     for( size_type j = 1; j <= num_cols; ++j )
00077       col_vecs_[j-1] = col_vecs[j-1];
00078   }
00079   else {
00080     for( size_type j = 1; j <= num_cols; ++j )
00081       col_vecs_[j-1] = space_cols->create_member();
00082   }
00083 }
00084 
00085 void MultiVectorMutableCols::set_uninitialized()
00086 {
00087   col_vecs_.resize(0);
00088   space_cols_ = Teuchos::null;
00089   space_rows_ = Teuchos::null;
00090 }
00091 
00092 // Overridden from MatrixBase
00093 
00094 size_type MultiVectorMutableCols::rows() const
00095 {
00096   return space_cols_.get() ? space_cols_->dim() : 0;
00097 }
00098 
00099 size_type MultiVectorMutableCols::cols() const
00100 {
00101   return space_rows_.get() ? space_rows_->dim() : 0;
00102 }
00103 
00104 // Overridden from MatrixOp
00105 
00106 const VectorSpace& MultiVectorMutableCols::space_cols() const
00107 {
00108   return *space_cols_;
00109 }
00110 
00111 const VectorSpace& MultiVectorMutableCols::space_rows() const
00112 {
00113   return *space_rows_;
00114 }
00115 
00116 void MultiVectorMutableCols::zero_out()
00117 {
00118   for( size_type j = 1; j <= col_vecs_.size(); ++j )
00119     col_vecs_[j-1]->zero();
00120 }
00121 
00122 void MultiVectorMutableCols::Mt_S( value_type alpha )
00123 {
00124   for( size_type j = 1; j <= col_vecs_.size(); ++j )
00125     LinAlgOpPack::Vt_S(col_vecs_[j-1].get(),alpha);
00126 }
00127 
00128 MatrixOp&
00129 MultiVectorMutableCols::operator=(const MatrixOp& mwo_rhs)
00130 {
00131   const MultiVector
00132     *mv_rhs = dynamic_cast<const MultiVector*>(&mwo_rhs);
00133   if(!mv_rhs)
00134     return MatrixOp::operator=(mwo_rhs);
00135   for( size_type j = 1; j <= col_vecs_.size(); ++j )
00136     *col_vecs_[j-1] = *mv_rhs->col(j);
00137   return *this;
00138 }
00139 
00140 MatrixOp::mat_mut_ptr_t
00141 MultiVectorMutableCols::clone()
00142 {
00143   TEUCHOS_TEST_FOR_EXCEPT(true); // ToDo: Implement!
00144   return Teuchos::null;
00145 }
00146 
00147 MatrixOp::mat_ptr_t
00148 MultiVectorMutableCols::clone() const
00149 {
00150   TEUCHOS_TEST_FOR_EXCEPT(true); // ToDo: Implement!
00151   return Teuchos::null;
00152 }
00153 
00154 void MultiVectorMutableCols::Vp_StMtV(
00155   VectorMutable* y, value_type a, BLAS_Cpp::Transp M_trans
00156   ,const Vector& x, value_type b
00157   ) const
00158 {
00159   using AbstractLinAlgPack::dot;
00160 
00161   // y = b*y
00162   LinAlgOpPack::Vt_S(y,b);
00163 
00164   if( M_trans == BLAS_Cpp::no_trans ) {
00165     //
00166     // y += a*M*x
00167     //
00168     // =>
00169     //
00170     // y += a * x(1) * M(:,1) + a * x(2) * M(:,2) + ...
00171     //
00172     for( size_type j = 1; j <= col_vecs_.size(); ++j )
00173       LinAlgOpPack::Vp_StV( y, a * x.get_ele(j), *col_vecs_[j-1] );
00174   }
00175   else {
00176     //
00177     // y += a*M'*x
00178     //
00179     // =>
00180     //
00181     // y(1) += a M(:,1)*x
00182     // y(2) += a M(:,2)*x
00183     // ...
00184     //
00185     for( size_type j = 1; j <= col_vecs_.size(); ++j )
00186       y->set_ele(
00187         j
00188         ,y->get_ele(j) + a * dot(*col_vecs_[j-1],x)
00189         );
00190   }
00191 }
00192 
00193 void MultiVectorMutableCols::Vp_StMtV(
00194   VectorMutable* y, value_type a, BLAS_Cpp::Transp M_trans
00195   ,const SpVectorSlice& x, value_type b
00196   ) const
00197 {
00198   using AbstractLinAlgPack::dot;
00199 
00200   // y = b*y
00201   LinAlgOpPack::Vt_S(y,b);
00202 
00203   if( M_trans == BLAS_Cpp::no_trans ) {
00204     //
00205     // y += a*M*x
00206     //
00207     // =>
00208     //
00209     // y += a * x(1) * M(:,1) + a * x(2) * M(:,2) + ...
00210     //
00211     SpVectorSlice::difference_type o = x.offset();
00212     for( SpVectorSlice::const_iterator itr = x.begin(); itr != x.end(); ++itr ) {
00213       const size_type j = itr->index() + o;
00214       LinAlgOpPack::Vp_StV( y, a * itr->value(), *col_vecs_[j-1] );
00215     }
00216   }
00217   else {
00218     //
00219     // y += a*M'*x
00220     //
00221     // =>
00222     //
00223     // y(1) += a M(:,1)*x
00224     // y(2) += a M(:,2)*x
00225     // ...
00226     //
00227     for( size_type j = 1; j <= col_vecs_.size(); ++j )
00228       y->set_ele(
00229         j
00230         ,y->get_ele(j) + a * dot(*col_vecs_[j-1],x)
00231         );
00232   }
00233 }
00234 
00235 bool MultiVectorMutableCols::syrk(
00236   BLAS_Cpp::Transp M_trans, value_type alpha
00237   , value_type beta, MatrixSymOp* sym_lhs ) const
00238 {
00239   using LinAlgOpPack::dot;
00240   MatrixSymOpGetGMSSymMutable
00241     *symwo_gms_lhs = dynamic_cast<MatrixSymOpGetGMSSymMutable*>(sym_lhs);
00242   if(!symwo_gms_lhs) {
00243     return MatrixOp::syrk(M_trans,alpha,beta,sym_lhs); // Boot it
00244   }
00245   MatrixDenseSymMutableEncap  DMatrixSliceSym(symwo_gms_lhs);
00246   const int num_vecs = this->col_vecs_.size();
00247   TEUCHOS_TEST_FOR_EXCEPTION(
00248     num_vecs != DMatrixSliceSym().rows(), std::logic_error
00249     ,"MultiVectorMutableCols::syrk(...) : Error, sizes do not match up!" );
00250   // Fill the upper or lower triangular region.
00251   if( DMatrixSliceSym().uplo() == BLAS_Cpp::upper ) {
00252     for( int i = 1; i <= num_vecs; ++i ) {
00253       for( int j = i; j <= num_vecs; ++j ) { // Upper triangle!
00254         DMatrixSliceSym().gms()(i,j) = beta * DMatrixSliceSym().gms()(i,j) + alpha * dot(*col_vecs_[i-1],*col_vecs_[j-1]);
00255       }
00256     }
00257   }
00258   else {
00259     for( int i = 1; i <= num_vecs; ++i ) {
00260       for( int j = 1; j <= i; ++j ) { // Lower triangle!
00261         DMatrixSliceSym().gms()(i,j) = beta * DMatrixSliceSym().gms()(i,j) + alpha * dot(*col_vecs_[i-1],*col_vecs_[j-1]);
00262       }
00263     }
00264   }
00265   return true;
00266 }
00267 
00268 // Overridden from MultiVector
00269 
00270 MultiVector::access_by_t
00271 MultiVectorMutableCols::access_by() const
00272 {
00273   return COL_ACCESS;
00274 }
00275 
00276 // Overridden from MultiVectorMutable
00277 
00278 MultiVectorMutable::vec_mut_ptr_t
00279 MultiVectorMutableCols::col(index_type j)
00280 {
00281   TEUCHOS_TEST_FOR_EXCEPTION( !(  1 <= j  && j <= col_vecs_.size() ), std::logic_error, "Error!" );
00282   return col_vecs_[j-1];
00283 }
00284 
00285 MultiVectorMutable::vec_mut_ptr_t
00286 MultiVectorMutableCols::row(index_type i)
00287 {
00288   return Teuchos::null;
00289 }
00290 
00291 MultiVectorMutable::vec_mut_ptr_t
00292 MultiVectorMutableCols::diag(int k)
00293 {
00294   return Teuchos::null;
00295 }
00296 
00297 MultiVectorMutable::multi_vec_mut_ptr_t
00298 MultiVectorMutableCols::mv_sub_view(const Range1D& row_rng, const Range1D& col_rng)
00299 {
00300 #ifdef TEUCHOS_DEBUG
00301   const size_type rows = this->rows();
00302   TEUCHOS_TEST_FOR_EXCEPTION(
00303     !( row_rng.full_range() || (row_rng.lbound() == 1 && row_rng.ubound() == rows) )
00304     ,std::logic_error, "Error, can't handle subrange on vectors yet!" );
00305 #endif
00306   return Teuchos::rcp(
00307     new MultiVectorMutableCols(
00308       space_cols_,space_rows_->sub_space(col_rng),&col_vecs_[col_rng.lbound()-1]
00309       ) );
00310 }
00311   
00312 } // end namespace AbstractLinAlgPack
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines