AbstractLinAlgPack_MultiVector.cpp

Go to the documentation of this file.
00001 // /////////////////////////////////////////////////////////////////////////
00002 // MultiVector.cpp
00003 
00004 #include <assert.h>
00005 
00006 #include "AbstractLinAlgPack_MultiVectorMutable.hpp"
00007 #include "AbstractLinAlgPack_MatrixSymDiag.hpp"
00008 #include "AbstractLinAlgPack_VectorMutable.hpp"
00009 #include "AbstractLinAlgPack_AssertOp.hpp"
00010 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
00011 #include "Teuchos_Workspace.hpp"
00012 #include "Teuchos_TestForException.hpp"
00013 
00014 namespace {
00015 
00016 // Map from EApplyBy to Transp
00017 inline
00018 BLAS_Cpp::Transp
00019 to_trans(AbstractLinAlgPack::EApplyBy apply_by)
00020 {
00021   return ( apply_by == AbstractLinAlgPack::APPLY_BY_ROW
00022       ? BLAS_Cpp::no_trans
00023       : BLAS_Cpp::trans
00024       );
00025 }
00026 
00027 // Return a row or a column vector from a multi-vector
00028 
00029 inline 
00030 AbstractLinAlgPack::MultiVector::vec_ptr_t
00031 vec(
00032   const AbstractLinAlgPack::MultiVector&      multi_vec
00033   ,const AbstractLinAlgPack::size_type        k
00034   ,AbstractLinAlgPack::EApplyBy               apply_by
00035   )
00036 {
00037   return ( apply_by == AbstractLinAlgPack::APPLY_BY_ROW
00038       ? multi_vec.row(k)
00039       : multi_vec.col(k)
00040       );
00041 }
00042 
00043 inline 
00044 AbstractLinAlgPack::MultiVectorMutable::vec_mut_ptr_t
00045 vec(
00046   AbstractLinAlgPack::MultiVectorMutable*         multi_vec
00047   ,const AbstractLinAlgPack::size_type            k
00048   ,AbstractLinAlgPack::EApplyBy                   apply_by
00049   )
00050 {
00051   return ( apply_by == AbstractLinAlgPack::APPLY_BY_ROW
00052       ? multi_vec->row(k)
00053       : multi_vec->col(k)
00054       );
00055 }
00056 
00057 // Implement a matrix-matrix multiplication with a diagonal matrix.
00058 //
00059 // op(C) = b*op(C) + a*D*op(B)
00060 //
00061 bool mat_vec(
00062   const AbstractLinAlgPack::value_type        &a
00063   ,const AbstractLinAlgPack::MatrixOp         &D_mwo  // Diagonal matrix?
00064   ,BLAS_Cpp::Transp                           D_trans
00065   ,const AbstractLinAlgPack::MultiVector      &B
00066   ,BLAS_Cpp::Transp                           B_trans
00067   ,const AbstractLinAlgPack::value_type       &b
00068   ,BLAS_Cpp::Transp                           C_trans
00069   ,AbstractLinAlgPack::MatrixOp               *C
00070   )
00071 {
00072   using BLAS_Cpp::no_trans;
00073   using BLAS_Cpp::trans;
00074 
00075   typedef AbstractLinAlgPack::MultiVector          MV;
00076   typedef AbstractLinAlgPack::MultiVectorMutable   MVM;
00077   using AbstractLinAlgPack::size_type;
00078   using AbstractLinAlgPack::Vector;
00079   using AbstractLinAlgPack::MatrixOp;
00080   using AbstractLinAlgPack::MultiVectorMutable;
00081   using AbstractLinAlgPack::MatrixSymDiag;
00082   using AbstractLinAlgPack::ele_wise_prod;
00083   using LinAlgOpPack::Vt_S;
00084   
00085   AbstractLinAlgPack::Mp_MtM_assert_compatibility(C,C_trans,D_mwo,D_trans,B,B_trans);
00086 
00087   MultiVectorMutable
00088     *Cmv = dynamic_cast<MultiVectorMutable*>(C);
00089   const MatrixSymDiag
00090     *D = dynamic_cast<const MatrixSymDiag*>(&D_mwo);
00091   if( !Cmv || !D || !(Cmv->access_by() & ( C_trans == no_trans ? MV::COL_ACCESS : MV::ROW_ACCESS ))
00092     || !(B.access_by() & ( B_trans == no_trans ? MV::COL_ACCESS : MV::ROW_ACCESS ))
00093     )
00094   {
00095     return false;
00096   }
00097   //
00098   // 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()
00099   //
00100   const Vector  &D_diag = D->diag();
00101   const size_type
00102     opC_cols = BLAS_Cpp::cols( Cmv->rows(), Cmv->cols(), C_trans );
00103   for( size_type j = 1; j <= opC_cols; ++j ) {
00104     MV::vec_ptr_t
00105       opB_col_j = ( B_trans == no_trans ? B.col(j)    : B.row(j) );
00106     MVM::vec_mut_ptr_t
00107       opC_col_j = ( C_trans == no_trans ? Cmv->col(j) : Cmv->row(j) );
00108     Vt_S( opC_col_j.get(), b );
00109     ele_wise_prod( a, D_diag, *opB_col_j, opC_col_j.get() );
00110   } 
00111   return true;
00112 }
00113 
00114 } // end namespace
00115 
00116 namespace AbstractLinAlgPack {
00117 
00118 MultiVector::multi_vec_ptr_t
00119 MultiVector::mv_clone() const
00120 {
00121   return Teuchos::null;
00122 }
00123 
00124 MultiVector::multi_vec_ptr_t
00125 MultiVector::mv_sub_view(const Range1D& row_rng, const Range1D& col_rng) const
00126 {
00127   TEST_FOR_EXCEPT(true); // ToDo: return a MultiVectorSubView object.
00128   // Note that the MultiVectorSubView class should derive from MatrixOpSubView
00129   // so that a client can rely on the MatrixOpSubView interface.
00130   return Teuchos::null;
00131 }
00132 
00133 void MultiVector::apply_op(
00134   EApplyBy apply_by, const RTOpPack::RTOp& prim_op
00135   ,const size_t num_multi_vecs,      const MultiVector*   multi_vecs[]
00136   ,const size_t num_targ_multi_vecs, MultiVectorMutable*  targ_multi_vecs[]
00137   ,RTOpPack::ReductTarget* reduct_objs[]
00138   ,const index_type prim_first_ele_in, const index_type prim_sub_dim_in, const index_type prim_global_offset_in
00139   ,const index_type sec_first_ele_in, const index_type sec_sub_dim_in
00140   ) const
00141 {
00142   using Teuchos::Workspace;
00143   Teuchos::WorkspaceStore* wss = Teuchos::get_default_workspace_store().get();
00144 
00145   // ToDo: Validate the input!
00146 
00147   // Get the primary and secondary dimmensions.
00148   const index_type
00149     prim_dim     = ( apply_by == APPLY_BY_ROW ? rows()          : cols()   ),
00150     sec_dim      = ( apply_by == APPLY_BY_ROW ? cols()          : rows()   ),
00151     prim_sub_dim = ( prim_sub_dim_in != 0     ? prim_sub_dim_in : prim_dim ),
00152     sec_sub_dim  = ( sec_sub_dim_in != 0      ? sec_sub_dim_in  : sec_dim  );
00153   TEST_FOR_EXCEPT( !( 0 < prim_sub_dim && prim_sub_dim <= prim_dim  ) );
00154   TEST_FOR_EXCEPT( !( 0 < sec_sub_dim  && sec_sub_dim  <= sec_dim  ) );
00155 
00156   //
00157   // Apply the reduction/transformation operator and trnasform the target
00158   // vectors and reduce each of the reduction objects.
00159   //
00160 
00161   Workspace<MultiVector::vec_ptr_t>             vecs_s(wss,num_multi_vecs);
00162   Workspace<const Vector*>                      vecs(wss,num_multi_vecs);
00163   Workspace<MultiVectorMutable::vec_mut_ptr_t>  targ_vecs_s(wss,num_targ_multi_vecs);
00164   Workspace<VectorMutable*>                     targ_vecs(wss,num_targ_multi_vecs);
00165 
00166   {for(size_type j = sec_first_ele_in; j <= sec_first_ele_in - 1 + sec_sub_dim; ++j) {
00167     // Fill the arrays of vector arguments 
00168     {for(size_type k = 0; k < num_multi_vecs; ++k) {
00169       vecs_s[k] = vec( *multi_vecs[k], j, apply_by );
00170       vecs[k] = vecs_s[k].get();
00171     }}
00172     {for(size_type k = 0; k < num_targ_multi_vecs; ++k) {
00173       targ_vecs_s[k] = vec( targ_multi_vecs[k], j, apply_by );
00174       targ_vecs[k] = targ_vecs_s[k].get();
00175     }}
00176     // Apply the reduction/transformation operator
00177     AbstractLinAlgPack::apply_op(
00178       prim_op
00179       ,num_multi_vecs,      num_multi_vecs      ? &vecs[0]      : NULL
00180       ,num_targ_multi_vecs, num_targ_multi_vecs ? &targ_vecs[0] : NULL
00181       ,reduct_objs ? reduct_objs[j-1] : NULL
00182       ,prim_first_ele_in, prim_sub_dim_in, prim_global_offset_in
00183       );
00184   }}
00185 
00186   // At this point all of the designated targ vectors in the target multi-vectors have
00187   // been transformed and all the reduction objects in reduct_obj[] have accumulated
00188   // the reductions.
00189 
00190 }
00191 
00192 void MultiVector::apply_op(
00193   EApplyBy apply_by, const RTOpPack::RTOp& prim_op, const RTOpPack::RTOp& sec_op
00194   ,const size_t num_multi_vecs,      const MultiVector*   multi_vecs[]
00195   ,const size_t num_targ_multi_vecs, MultiVectorMutable*  targ_multi_vecs[]
00196   ,RTOpPack::ReductTarget *reduct_obj
00197   ,const index_type prim_first_ele_in, const index_type prim_sub_dim_in, const index_type prim_global_offset_in
00198   ,const index_type sec_first_ele_in, const index_type sec_sub_dim_in
00199   ) const
00200 {
00201   using Teuchos::Workspace;
00202   Teuchos::WorkspaceStore* wss = Teuchos::get_default_workspace_store().get();
00203 
00204   // ToDo: Validate the input!
00205 
00206   // Get the primary and secondary dimmensions.
00207   const index_type
00208     prim_dim    = ( apply_by == APPLY_BY_ROW ? rows()         : cols()  ),
00209     sec_dim     = ( apply_by == APPLY_BY_ROW ? cols()         : rows()  ),
00210     sec_sub_dim = ( sec_sub_dim_in != 0      ? sec_sub_dim_in : sec_dim );
00211   TEST_FOR_EXCEPT( !( 0 < sec_sub_dim && sec_sub_dim <= sec_dim  ) );
00212 
00213   // Create a temporary buffer for the reduction objects of the primary reduction
00214   // so that we can call the companion version of this method.
00215   Workspace<Teuchos::RCP<RTOpPack::ReductTarget> >   rcp_reduct_objs(wss,sec_sub_dim);
00216   Workspace<RTOpPack::ReductTarget*>                         reduct_objs(wss,sec_sub_dim);
00217   for(index_type k = 0; k < sec_sub_dim; ++k) {
00218     rcp_reduct_objs[k] = prim_op.reduct_obj_create();
00219     reduct_objs[k] = &*rcp_reduct_objs[k];
00220   }
00221   
00222   // Call the campanion version that accepts an array of reduction objects
00223   this->apply_op(
00224     apply_by, prim_op
00225     ,num_multi_vecs,       multi_vecs
00226     ,num_targ_multi_vecs,  targ_multi_vecs
00227     ,&reduct_objs[0]
00228     ,prim_first_ele_in, prim_sub_dim_in, prim_global_offset_in
00229     ,sec_first_ele_in,  sec_sub_dim_in
00230     );
00231 
00232   // Reduce all the reduction objects using the secondary reduction operator
00233   // into one reduction object and free the intermedate reduction objects.
00234   for(index_type k = 0; k < sec_sub_dim; ++k) {
00235     sec_op.reduce_reduct_objs( *reduct_objs[k], reduct_obj );
00236   }
00237 }
00238 
00239 // Overridden form MatrixOp
00240 
00241 MatrixOp::mat_ptr_t
00242 MultiVector::clone() const
00243 {
00244   return this->mv_clone();
00245 }
00246 
00247 MatrixOp::mat_ptr_t
00248 MultiVector::sub_view(const Range1D& row_rng, const Range1D& col_rng) const
00249 {
00250   return mv_sub_view(row_rng,col_rng);
00251 }
00252 
00253 bool MultiVector::Mp_StMtM(
00254   MatrixOp* mwo_lhs, value_type alpha
00255   ,const MatrixOp& mwo_rhs1, BLAS_Cpp::Transp trans_rhs1
00256   ,BLAS_Cpp::Transp trans_rhs2
00257   ,value_type beta
00258   ) const
00259 {
00260   return mat_vec(
00261     alpha
00262     ,mwo_rhs1,trans_rhs1
00263     ,*this,trans_rhs2
00264     ,beta,BLAS_Cpp::no_trans,mwo_lhs
00265     );
00266 }
00267 
00268 bool MultiVector::Mp_StMtM(
00269   MatrixOp* mwo_lhs, value_type alpha
00270   ,BLAS_Cpp::Transp trans_rhs1
00271   ,const MatrixOp& mwo_rhs2, BLAS_Cpp::Transp trans_rhs2
00272   ,value_type beta
00273   ) const
00274 {
00275   return mat_vec(
00276     alpha
00277     ,mwo_rhs2,BLAS_Cpp::trans_not(trans_rhs2)
00278     ,*this,BLAS_Cpp::trans_not(trans_rhs1)
00279     ,beta,BLAS_Cpp::trans,mwo_lhs
00280     );
00281 }
00282 
00283 } // end namespace AbstractLinAlgPack
00284 
00285 // nonmembers
00286 
00287 void AbstractLinAlgPack::apply_op(
00288   EApplyBy                        apply_by
00289   ,const RTOpPack::RTOp           &primary_op
00290   ,const size_t                   num_multi_vecs
00291   ,const MultiVector*             multi_vecs[]
00292   ,const size_t                   num_targ_multi_vecs
00293   ,MultiVectorMutable*            targ_multi_vecs[]
00294   ,RTOpPack::ReductTarget*        reduct_objs[]
00295   ,const index_type               primary_first_ele
00296   ,const index_type               primary_sub_dim
00297   ,const index_type               primary_global_offset
00298   ,const index_type               secondary_first_ele
00299   ,const index_type               secondary_sub_dim
00300   )
00301 {
00302   if(num_multi_vecs)
00303     multi_vecs[0]->apply_op(
00304       apply_by,primary_op
00305       ,num_multi_vecs,multi_vecs,num_targ_multi_vecs,targ_multi_vecs
00306       ,reduct_objs,primary_first_ele,primary_sub_dim,primary_global_offset
00307       ,secondary_first_ele,secondary_sub_dim
00308       );
00309   else if(num_targ_multi_vecs)
00310     targ_multi_vecs[0]->apply_op(
00311       apply_by,primary_op
00312       ,num_multi_vecs,multi_vecs,num_targ_multi_vecs,targ_multi_vecs
00313       ,reduct_objs,primary_first_ele,primary_sub_dim,primary_global_offset
00314       ,secondary_first_ele,secondary_sub_dim
00315       );
00316 }
00317 
00318 void AbstractLinAlgPack::apply_op(
00319   EApplyBy                        apply_by
00320   ,const RTOpPack::RTOp           &primary_op
00321   ,const RTOpPack::RTOp           &secondary_op
00322   ,const size_t                   num_multi_vecs
00323   ,const MultiVector*             multi_vecs[]
00324   ,const size_t                   num_targ_multi_vecs
00325   ,MultiVectorMutable*            targ_multi_vecs[]
00326   ,RTOpPack::ReductTarget         *reduct_obj
00327   ,const index_type               primary_first_ele
00328   ,const index_type               primary_sub_dim
00329   ,const index_type               primary_global_offset
00330   ,const index_type               secondary_first_ele
00331   ,const index_type               secondary_sub_dim
00332   )
00333 {
00334   if(num_multi_vecs)
00335     multi_vecs[0]->apply_op(
00336       apply_by,primary_op,secondary_op
00337       ,num_multi_vecs,multi_vecs,num_targ_multi_vecs,targ_multi_vecs
00338       ,reduct_obj,primary_first_ele,primary_sub_dim,primary_global_offset
00339       ,secondary_first_ele,secondary_sub_dim
00340       );
00341   else if(num_targ_multi_vecs)
00342     targ_multi_vecs[0]->apply_op(
00343       apply_by,primary_op,secondary_op
00344       ,num_multi_vecs,multi_vecs,num_targ_multi_vecs,targ_multi_vecs
00345       ,reduct_obj,primary_first_ele,primary_sub_dim,primary_global_offset
00346       ,secondary_first_ele,secondary_sub_dim
00347       );
00348 }

Generated on Tue Jul 13 09:30:50 2010 for MOOCHO (Single Doxygen Collection) by  doxygen 1.4.7