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

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