MOOCHO (Single Doxygen Collection) Version of the Day
AbstractLinAlgPack_MultiVector.cpp
Go to the documentation of this file.
00001 /*
00002 // @HEADER
00003 // ***********************************************************************
00004 // 
00005 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
00006 //                  Copyright (2003) Sandia Corporation
00007 // 
00008 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00009 // license for use of this work by or on behalf of the U.S. Government.
00010 // 
00011 // Redistribution and use in source and binary forms, with or without
00012 // modification, are permitted provided that the following conditions are
00013 // met:
00014 //
00015 // 1. Redistributions of source code must retain the above copyright
00016 // notice, this list of conditions and the following disclaimer.
00017 //
00018 // 2. Redistributions in binary form must reproduce the above copyright
00019 // notice, this list of conditions and the following disclaimer in the
00020 // documentation and/or other materials provided with the distribution.
00021 //
00022 // 3. Neither the name of the Corporation nor the names of the
00023 // contributors may be used to endorse or promote products derived from
00024 // this software without specific prior written permission.
00025 //
00026 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00027 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00028 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00029 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00030 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00031 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00032 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00033 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00034 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00035 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00036 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00037 //
00038 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
00039 // 
00040 // ***********************************************************************
00041 // @HEADER
00042 */
00043 
00044 // /////////////////////////////////////////////////////////////////////////
00045 // MultiVector.cpp
00046 
00047 #include <assert.h>
00048 
00049 #include "AbstractLinAlgPack_MultiVectorMutable.hpp"
00050 #include "AbstractLinAlgPack_MatrixSymDiag.hpp"
00051 #include "AbstractLinAlgPack_VectorMutable.hpp"
00052 #include "AbstractLinAlgPack_AssertOp.hpp"
00053 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
00054 #include "Teuchos_Workspace.hpp"
00055 #include "Teuchos_Assert.hpp"
00056 
00057 namespace {
00058 
00059 // Map from EApplyBy to Transp
00060 inline
00061 BLAS_Cpp::Transp
00062 to_trans(AbstractLinAlgPack::EApplyBy apply_by)
00063 {
00064   return ( apply_by == AbstractLinAlgPack::APPLY_BY_ROW
00065       ? BLAS_Cpp::no_trans
00066       : BLAS_Cpp::trans
00067       );
00068 }
00069 
00070 // Return a row or a column vector from a multi-vector
00071 
00072 inline 
00073 AbstractLinAlgPack::MultiVector::vec_ptr_t
00074 vec(
00075   const AbstractLinAlgPack::MultiVector&      multi_vec
00076   ,const AbstractLinAlgPack::size_type        k
00077   ,AbstractLinAlgPack::EApplyBy               apply_by
00078   )
00079 {
00080   return ( apply_by == AbstractLinAlgPack::APPLY_BY_ROW
00081       ? multi_vec.row(k)
00082       : multi_vec.col(k)
00083       );
00084 }
00085 
00086 inline 
00087 AbstractLinAlgPack::MultiVectorMutable::vec_mut_ptr_t
00088 vec(
00089   AbstractLinAlgPack::MultiVectorMutable*         multi_vec
00090   ,const AbstractLinAlgPack::size_type            k
00091   ,AbstractLinAlgPack::EApplyBy                   apply_by
00092   )
00093 {
00094   return ( apply_by == AbstractLinAlgPack::APPLY_BY_ROW
00095       ? multi_vec->row(k)
00096       : multi_vec->col(k)
00097       );
00098 }
00099 
00100 // Implement a matrix-matrix multiplication with a diagonal matrix.
00101 //
00102 // op(C) = b*op(C) + a*D*op(B)
00103 //
00104 bool mat_vec(
00105   const AbstractLinAlgPack::value_type        &a
00106   ,const AbstractLinAlgPack::MatrixOp         &D_mwo  // Diagonal matrix?
00107   ,BLAS_Cpp::Transp                           D_trans
00108   ,const AbstractLinAlgPack::MultiVector      &B
00109   ,BLAS_Cpp::Transp                           B_trans
00110   ,const AbstractLinAlgPack::value_type       &b
00111   ,BLAS_Cpp::Transp                           C_trans
00112   ,AbstractLinAlgPack::MatrixOp               *C
00113   )
00114 {
00115   using BLAS_Cpp::no_trans;
00116   using BLAS_Cpp::trans;
00117 
00118   typedef AbstractLinAlgPack::MultiVector          MV;
00119   typedef AbstractLinAlgPack::MultiVectorMutable   MVM;
00120   using AbstractLinAlgPack::size_type;
00121   using AbstractLinAlgPack::Vector;
00122   using AbstractLinAlgPack::MatrixOp;
00123   using AbstractLinAlgPack::MultiVectorMutable;
00124   using AbstractLinAlgPack::MatrixSymDiag;
00125   using AbstractLinAlgPack::ele_wise_prod;
00126   using LinAlgOpPack::Vt_S;
00127   
00128   AbstractLinAlgPack::Mp_MtM_assert_compatibility(C,C_trans,D_mwo,D_trans,B,B_trans);
00129 
00130   MultiVectorMutable
00131     *Cmv = dynamic_cast<MultiVectorMutable*>(C);
00132   const MatrixSymDiag
00133     *D = dynamic_cast<const MatrixSymDiag*>(&D_mwo);
00134   if( !Cmv || !D || !(Cmv->access_by() & ( C_trans == no_trans ? MV::COL_ACCESS : MV::ROW_ACCESS ))
00135     || !(B.access_by() & ( B_trans == no_trans ? MV::COL_ACCESS : MV::ROW_ACCESS ))
00136     )
00137   {
00138     return false;
00139   }
00140   //
00141   // op(C).col(j) = b*op(C).col(j) + a*ele_wise_prod(D_diag,op(B).col(j)), for j = 1...op(C).cols()
00142   //
00143   const Vector  &D_diag = D->diag();
00144   const size_type
00145     opC_cols = BLAS_Cpp::cols( Cmv->rows(), Cmv->cols(), C_trans );
00146   for( size_type j = 1; j <= opC_cols; ++j ) {
00147     MV::vec_ptr_t
00148       opB_col_j = ( B_trans == no_trans ? B.col(j)    : B.row(j) );
00149     MVM::vec_mut_ptr_t
00150       opC_col_j = ( C_trans == no_trans ? Cmv->col(j) : Cmv->row(j) );
00151     Vt_S( opC_col_j.get(), b );
00152     ele_wise_prod( a, D_diag, *opB_col_j, opC_col_j.get() );
00153   } 
00154   return true;
00155 }
00156 
00157 } // end namespace
00158 
00159 namespace AbstractLinAlgPack {
00160 
00161 MultiVector::multi_vec_ptr_t
00162 MultiVector::mv_clone() const
00163 {
00164   return Teuchos::null;
00165 }
00166 
00167 MultiVector::multi_vec_ptr_t
00168 MultiVector::mv_sub_view(const Range1D& row_rng, const Range1D& col_rng) const
00169 {
00170   TEUCHOS_TEST_FOR_EXCEPT(true); // ToDo: return a MultiVectorSubView object.
00171   // Note that the MultiVectorSubView class should derive from MatrixOpSubView
00172   // so that a client can rely on the MatrixOpSubView interface.
00173   return Teuchos::null;
00174 }
00175 
00176 void MultiVector::apply_op(
00177   EApplyBy apply_by, const RTOpPack::RTOp& prim_op
00178   ,const size_t num_multi_vecs,      const MultiVector*   multi_vecs[]
00179   ,const size_t num_targ_multi_vecs, MultiVectorMutable*  targ_multi_vecs[]
00180   ,RTOpPack::ReductTarget* reduct_objs[]
00181   ,const index_type prim_first_ele_in, const index_type prim_sub_dim_in, const index_type prim_global_offset_in
00182   ,const index_type sec_first_ele_in, const index_type sec_sub_dim_in
00183   ) const
00184 {
00185   using Teuchos::Workspace;
00186   Teuchos::WorkspaceStore* wss = Teuchos::get_default_workspace_store().get();
00187 
00188   // ToDo: Validate the input!
00189 
00190   // Get the primary and secondary dimmensions.
00191   const index_type
00192     prim_dim     = ( apply_by == APPLY_BY_ROW ? rows()          : cols()   ),
00193     sec_dim      = ( apply_by == APPLY_BY_ROW ? cols()          : rows()   ),
00194     prim_sub_dim = ( prim_sub_dim_in != 0     ? prim_sub_dim_in : prim_dim ),
00195     sec_sub_dim  = ( sec_sub_dim_in != 0      ? sec_sub_dim_in  : sec_dim  );
00196   TEUCHOS_TEST_FOR_EXCEPT( !( 0 < prim_sub_dim && prim_sub_dim <= prim_dim  ) );
00197   TEUCHOS_TEST_FOR_EXCEPT( !( 0 < sec_sub_dim  && sec_sub_dim  <= sec_dim  ) );
00198 
00199   //
00200   // Apply the reduction/transformation operator and trnasform the target
00201   // vectors and reduce each of the reduction objects.
00202   //
00203 
00204   Workspace<MultiVector::vec_ptr_t>             vecs_s(wss,num_multi_vecs);
00205   Workspace<const Vector*>                      vecs(wss,num_multi_vecs);
00206   Workspace<MultiVectorMutable::vec_mut_ptr_t>  targ_vecs_s(wss,num_targ_multi_vecs);
00207   Workspace<VectorMutable*>                     targ_vecs(wss,num_targ_multi_vecs);
00208 
00209   {for(size_type j = sec_first_ele_in; j <= sec_first_ele_in - 1 + sec_sub_dim; ++j) {
00210     // Fill the arrays of vector arguments 
00211     {for(size_type k = 0; k < num_multi_vecs; ++k) {
00212       vecs_s[k] = vec( *multi_vecs[k], j, apply_by );
00213       vecs[k] = vecs_s[k].get();
00214     }}
00215     {for(size_type k = 0; k < num_targ_multi_vecs; ++k) {
00216       targ_vecs_s[k] = vec( targ_multi_vecs[k], j, apply_by );
00217       targ_vecs[k] = targ_vecs_s[k].get();
00218     }}
00219     // Apply the reduction/transformation operator
00220     AbstractLinAlgPack::apply_op(
00221       prim_op
00222       ,num_multi_vecs,      num_multi_vecs      ? &vecs[0]      : NULL
00223       ,num_targ_multi_vecs, num_targ_multi_vecs ? &targ_vecs[0] : NULL
00224       ,reduct_objs ? reduct_objs[j-1] : NULL
00225       ,prim_first_ele_in, prim_sub_dim_in, prim_global_offset_in
00226       );
00227   }}
00228 
00229   // At this point all of the designated targ vectors in the target multi-vectors have
00230   // been transformed and all the reduction objects in reduct_obj[] have accumulated
00231   // the reductions.
00232 
00233 }
00234 
00235 void MultiVector::apply_op(
00236   EApplyBy apply_by, const RTOpPack::RTOp& prim_op, const RTOpPack::RTOp& sec_op
00237   ,const size_t num_multi_vecs,      const MultiVector*   multi_vecs[]
00238   ,const size_t num_targ_multi_vecs, MultiVectorMutable*  targ_multi_vecs[]
00239   ,RTOpPack::ReductTarget *reduct_obj
00240   ,const index_type prim_first_ele_in, const index_type prim_sub_dim_in, const index_type prim_global_offset_in
00241   ,const index_type sec_first_ele_in, const index_type sec_sub_dim_in
00242   ) const
00243 {
00244   using Teuchos::Workspace;
00245   Teuchos::WorkspaceStore* wss = Teuchos::get_default_workspace_store().get();
00246 
00247   // ToDo: Validate the input!
00248 
00249   // Get the primary and secondary dimmensions.
00250   const index_type
00251     prim_dim    = ( apply_by == APPLY_BY_ROW ? rows()         : cols()  ),
00252     sec_dim     = ( apply_by == APPLY_BY_ROW ? cols()         : rows()  ),
00253     sec_sub_dim = ( sec_sub_dim_in != 0      ? sec_sub_dim_in : sec_dim );
00254   TEUCHOS_TEST_FOR_EXCEPT( !( 0 < sec_sub_dim && sec_sub_dim <= sec_dim  ) );
00255 
00256   // Create a temporary buffer for the reduction objects of the primary reduction
00257   // so that we can call the companion version of this method.
00258   Workspace<Teuchos::RCP<RTOpPack::ReductTarget> >   rcp_reduct_objs(wss,sec_sub_dim);
00259   Workspace<RTOpPack::ReductTarget*>                         reduct_objs(wss,sec_sub_dim);
00260   for(index_type k = 0; k < sec_sub_dim; ++k) {
00261     rcp_reduct_objs[k] = prim_op.reduct_obj_create();
00262     reduct_objs[k] = &*rcp_reduct_objs[k];
00263   }
00264   
00265   // Call the campanion version that accepts an array of reduction objects
00266   this->apply_op(
00267     apply_by, prim_op
00268     ,num_multi_vecs,       multi_vecs
00269     ,num_targ_multi_vecs,  targ_multi_vecs
00270     ,&reduct_objs[0]
00271     ,prim_first_ele_in, prim_sub_dim_in, prim_global_offset_in
00272     ,sec_first_ele_in,  sec_sub_dim_in
00273     );
00274 
00275   // Reduce all the reduction objects using the secondary reduction operator
00276   // into one reduction object and free the intermedate reduction objects.
00277   for(index_type k = 0; k < sec_sub_dim; ++k) {
00278     sec_op.reduce_reduct_objs( *reduct_objs[k], Teuchos::ptr(reduct_obj) );
00279   }
00280 }
00281 
00282 // Overridden form MatrixOp
00283 
00284 MatrixOp::mat_ptr_t
00285 MultiVector::clone() const
00286 {
00287   return this->mv_clone();
00288 }
00289 
00290 MatrixOp::mat_ptr_t
00291 MultiVector::sub_view(const Range1D& row_rng, const Range1D& col_rng) const
00292 {
00293   return mv_sub_view(row_rng,col_rng);
00294 }
00295 
00296 bool MultiVector::Mp_StMtM(
00297   MatrixOp* mwo_lhs, value_type alpha
00298   ,const MatrixOp& mwo_rhs1, BLAS_Cpp::Transp trans_rhs1
00299   ,BLAS_Cpp::Transp trans_rhs2
00300   ,value_type beta
00301   ) const
00302 {
00303   return mat_vec(
00304     alpha
00305     ,mwo_rhs1,trans_rhs1
00306     ,*this,trans_rhs2
00307     ,beta,BLAS_Cpp::no_trans,mwo_lhs
00308     );
00309 }
00310 
00311 bool MultiVector::Mp_StMtM(
00312   MatrixOp* mwo_lhs, value_type alpha
00313   ,BLAS_Cpp::Transp trans_rhs1
00314   ,const MatrixOp& mwo_rhs2, BLAS_Cpp::Transp trans_rhs2
00315   ,value_type beta
00316   ) const
00317 {
00318   return mat_vec(
00319     alpha
00320     ,mwo_rhs2,BLAS_Cpp::trans_not(trans_rhs2)
00321     ,*this,BLAS_Cpp::trans_not(trans_rhs1)
00322     ,beta,BLAS_Cpp::trans,mwo_lhs
00323     );
00324 }
00325 
00326 } // end namespace AbstractLinAlgPack
00327 
00328 // nonmembers
00329 
00330 void AbstractLinAlgPack::apply_op(
00331   EApplyBy                        apply_by
00332   ,const RTOpPack::RTOp           &primary_op
00333   ,const size_t                   num_multi_vecs
00334   ,const MultiVector*             multi_vecs[]
00335   ,const size_t                   num_targ_multi_vecs
00336   ,MultiVectorMutable*            targ_multi_vecs[]
00337   ,RTOpPack::ReductTarget*        reduct_objs[]
00338   ,const index_type               primary_first_ele
00339   ,const index_type               primary_sub_dim
00340   ,const index_type               primary_global_offset
00341   ,const index_type               secondary_first_ele
00342   ,const index_type               secondary_sub_dim
00343   )
00344 {
00345   if(num_multi_vecs)
00346     multi_vecs[0]->apply_op(
00347       apply_by,primary_op
00348       ,num_multi_vecs,multi_vecs,num_targ_multi_vecs,targ_multi_vecs
00349       ,reduct_objs,primary_first_ele,primary_sub_dim,primary_global_offset
00350       ,secondary_first_ele,secondary_sub_dim
00351       );
00352   else if(num_targ_multi_vecs)
00353     targ_multi_vecs[0]->apply_op(
00354       apply_by,primary_op
00355       ,num_multi_vecs,multi_vecs,num_targ_multi_vecs,targ_multi_vecs
00356       ,reduct_objs,primary_first_ele,primary_sub_dim,primary_global_offset
00357       ,secondary_first_ele,secondary_sub_dim
00358       );
00359 }
00360 
00361 void AbstractLinAlgPack::apply_op(
00362   EApplyBy                        apply_by
00363   ,const RTOpPack::RTOp           &primary_op
00364   ,const RTOpPack::RTOp           &secondary_op
00365   ,const size_t                   num_multi_vecs
00366   ,const MultiVector*             multi_vecs[]
00367   ,const size_t                   num_targ_multi_vecs
00368   ,MultiVectorMutable*            targ_multi_vecs[]
00369   ,RTOpPack::ReductTarget         *reduct_obj
00370   ,const index_type               primary_first_ele
00371   ,const index_type               primary_sub_dim
00372   ,const index_type               primary_global_offset
00373   ,const index_type               secondary_first_ele
00374   ,const index_type               secondary_sub_dim
00375   )
00376 {
00377   if(num_multi_vecs)
00378     multi_vecs[0]->apply_op(
00379       apply_by,primary_op,secondary_op
00380       ,num_multi_vecs,multi_vecs,num_targ_multi_vecs,targ_multi_vecs
00381       ,reduct_obj,primary_first_ele,primary_sub_dim,primary_global_offset
00382       ,secondary_first_ele,secondary_sub_dim
00383       );
00384   else if(num_targ_multi_vecs)
00385     targ_multi_vecs[0]->apply_op(
00386       apply_by,primary_op,secondary_op
00387       ,num_multi_vecs,multi_vecs,num_targ_multi_vecs,targ_multi_vecs
00388       ,reduct_obj,primary_first_ele,primary_sub_dim,primary_global_offset
00389       ,secondary_first_ele,secondary_sub_dim
00390       );
00391 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines