Thyra_MultiVectorStdOps.hpp

00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 //    Thyra: Interfaces and Support for Abstract Numerical Algorithms
00005 //                 Copyright (2004) 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 Michael A. Heroux (maherou@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 #ifndef THYRA_MULTI_VECTOR_STD_OPS_HPP
00030 #define THYRA_MULTI_VECTOR_STD_OPS_HPP
00031 
00032 #include "Thyra_MultiVectorStdOpsDecl.hpp"
00033 #include "Thyra_VectorStdOps.hpp"
00034 #include "Thyra_VectorSpaceBase.hpp"
00035 #include "Thyra_VectorStdOps.hpp"
00036 #include "Thyra_MultiVectorBase.hpp"
00037 #include "RTOpPack_ROpSum.hpp"
00038 #include "RTOpPack_ROpDotProd.hpp"
00039 #include "RTOpPack_ROpNorm1.hpp"
00040 #include "RTOpPack_ROpNormInf.hpp"
00041 #include "RTOpPack_TOpAssignScalar.hpp"
00042 #include "RTOpPack_TOpAssignVectors.hpp"
00043 #include "RTOpPack_TOpAXPY.hpp"
00044 #include "RTOpPack_TOpLinearCombination.hpp"
00045 #include "RTOpPack_TOpScaleVector.hpp"
00046 #include "Teuchos_TestForException.hpp"
00047 #include "Teuchos_Assert.hpp"
00048 #include "Teuchos_as.hpp"
00049 
00050 
00051 template<class Scalar>
00052 void Thyra::norms( const MultiVectorBase<Scalar>& V,
00053   const ArrayView<typename ScalarTraits<Scalar>::magnitudeType> &norms )
00054 {
00055   typedef Teuchos::ScalarTraits<Scalar> ST;
00056   const int m = V.domain()->dim();
00057   Array<Scalar> prods(m);
00058   V.range()->scalarProds(V,V,&prods[0]);
00059   for ( int j = 0; j < m; ++j )
00060     norms[j] = ST::magnitude(ST::squareroot(prods[j]));
00061 }
00062 
00063 
00064 template<class Scalar, class NormOp>
00065 void Thyra::reductions( const MultiVectorBase<Scalar>& V, const NormOp &op,
00066   const ArrayView<typename ScalarTraits<Scalar>::magnitudeType> &norms )
00067 {
00068   using Teuchos::tuple; using Teuchos::ptrInArg; using Teuchos::null;
00069   const int m = V.domain()->dim();
00070   Array<RCP<RTOpPack::ReductTarget> > rcp_op_targs(m);
00071   Array<Ptr<RTOpPack::ReductTarget> > op_targs(m);
00072   for( int kc = 0; kc < m; ++kc ) {
00073     rcp_op_targs[kc] = op.reduct_obj_create();
00074     op_targs[kc] = rcp_op_targs[kc].ptr();
00075   }
00076   applyOp<Scalar>(op, tuple(ptrInArg(V)),
00077     ArrayView<Ptr<MultiVectorBase<Scalar> > >(null),
00078     op_targs );
00079   for( int kc = 0; kc < m; ++kc ) {
00080     norms[kc] = op(*op_targs[kc]);
00081   }
00082 }
00083 
00084 
00085 template<class Scalar>
00086 void Thyra::dots( const MultiVectorBase<Scalar>& V1, const MultiVectorBase<Scalar>& V2,
00087   const ArrayView<Scalar> &dots )
00088 {
00089   using Teuchos::tuple; using Teuchos::ptrInArg; using Teuchos::null;
00090   const int m = V1.domain()->dim();
00091   RTOpPack::ROpDotProd<Scalar> dot_op;
00092   Array<RCP<RTOpPack::ReductTarget> > rcp_dot_targs(m);
00093   Array<Ptr<RTOpPack::ReductTarget> > dot_targs(m);
00094   for( int kc = 0; kc < m; ++kc ) {
00095     rcp_dot_targs[kc] = dot_op.reduct_obj_create();
00096     dot_targs[kc] = rcp_dot_targs[kc].ptr();
00097   }
00098   applyOp<Scalar>( dot_op, tuple(ptrInArg(V1), ptrInArg(V2)),
00099     ArrayView<Ptr<MultiVectorBase<Scalar> > >(null),
00100     dot_targs );
00101   for( int kc = 0; kc < m; ++kc ) {
00102     dots[kc] = dot_op(*dot_targs[kc]);
00103   }
00104 }
00105 
00106 
00107 template<class Scalar>
00108 void Thyra::sums( const MultiVectorBase<Scalar>& V, const ArrayView<Scalar> &sums )
00109 {
00110   using Teuchos::tuple; using Teuchos::ptrInArg; using Teuchos::null;
00111   const int m = V.domain()->dim();
00112   RTOpPack::ROpSum<Scalar> sum_op;
00113   Array<RCP<RTOpPack::ReductTarget> > rcp_op_targs(m);
00114   Array<Ptr<RTOpPack::ReductTarget> > op_targs(m);
00115   for( int kc = 0; kc < m; ++kc ) {
00116     rcp_op_targs[kc] = sum_op.reduct_obj_create();
00117     op_targs[kc] = rcp_op_targs[kc].ptr();
00118   }
00119   applyOp<Scalar>(sum_op, tuple(ptrInArg(V)), null, op_targs );
00120   for( int kc = 0; kc < m; ++kc ) {
00121     sums[kc] = sum_op(*op_targs[kc]);
00122   }
00123 }
00124 
00125 
00126 template<class Scalar>
00127 typename Teuchos::ScalarTraits<Scalar>::magnitudeType
00128 Thyra::norm_1( const MultiVectorBase<Scalar>& V )
00129 {
00130   using Teuchos::tuple; using Teuchos::ptrInArg; using Teuchos::null;
00131   // Primary column-wise reduction (sum of absolute values)
00132   RTOpPack::ROpNorm1<Scalar> sum_abs_op;
00133   // Secondary reduction (max over all columns = induced norm_1 matrix norm)
00134   RTOpPack::ROpNormInf<Scalar> max_op;
00135   // Reduction object (must be same for both sum_abs and max_targ objects)
00136   RCP<RTOpPack::ReductTarget>
00137     max_targ = max_op.reduct_obj_create();
00138   // Perform the reductions
00139   applyOp<Scalar>(sum_abs_op, max_op, tuple(ptrInArg(V)), null, max_targ.ptr() );
00140   // Return the final value
00141   return max_op(*max_targ);
00142 }
00143 
00144 
00145 template<class Scalar>
00146 void Thyra::scale( Scalar alpha, const Ptr<MultiVectorBase<Scalar> > &V )
00147 {
00148   using Teuchos::tuple; using Teuchos::null;
00149   typedef ScalarTraits<Scalar> ST;
00150   if (alpha==ST::zero()) {
00151     assign( V, ST::zero() );
00152     return;
00153   }
00154   if (alpha==ST::one()) {
00155     return;
00156   }
00157   RTOpPack::TOpScaleVector<Scalar> scale_vector_op(alpha);
00158   applyOp<Scalar>(scale_vector_op,
00159     ArrayView<Ptr<const MultiVectorBase<Scalar> > >(null),
00160     tuple(V), null );
00161 }
00162 
00163 
00164 template<class Scalar>
00165 void Thyra::scaleUpdate( const VectorBase<Scalar>& a,
00166   const MultiVectorBase<Scalar>& U, const Ptr<MultiVectorBase<Scalar> > &V )
00167 {
00168 #ifdef TEUCHOS_DEBUG
00169   bool is_compatible = U.range()->isCompatible(*a.space());
00170   TEST_FOR_EXCEPTION(
00171     !is_compatible, Exceptions::IncompatibleVectorSpaces,
00172     "update(...), Error, U.range()->isCompatible(*a.space())==false" );
00173   is_compatible = U.range()->isCompatible(*V->range());
00174   TEST_FOR_EXCEPTION(
00175     !is_compatible, Exceptions::IncompatibleVectorSpaces,
00176     "update(...), Error, U.range()->isCompatible((V->range())==false" );
00177   is_compatible = U.domain()->isCompatible(*V->domain());
00178   TEST_FOR_EXCEPTION(
00179     !is_compatible, Exceptions::IncompatibleVectorSpaces,
00180     "update(...), Error, U.domain().isCompatible(V->domain())==false" );
00181 #endif
00182   const int m = U.domain()->dim();
00183   for( int j = 0; j < m; ++j ) {
00184     ele_wise_prod<Scalar>( 1.0, a, *U.col(j), V->col(j).ptr() ); 
00185   }
00186 }
00187 
00188 
00189 template<class Scalar>
00190 void Thyra::assign( const Ptr<MultiVectorBase<Scalar> > &V, Scalar alpha )
00191 {
00192   using Teuchos::tuple; using Teuchos::null;
00193   RTOpPack::TOpAssignScalar<Scalar> assign_scalar_op(alpha);
00194   applyOp<Scalar>(assign_scalar_op,
00195     ArrayView<Ptr<const MultiVectorBase<Scalar> > >(null),
00196     tuple(V), null);
00197 }
00198 
00199 
00200 template<class Scalar>
00201 void Thyra::assign( const Ptr<MultiVectorBase<Scalar> > &V,
00202   const MultiVectorBase<Scalar>& U )
00203 {
00204   using Teuchos::tuple; using Teuchos::ptrInArg; using Teuchos::null;
00205   RTOpPack::TOpAssignVectors<Scalar> assign_vectors_op;
00206   applyOp<Scalar>( assign_vectors_op, tuple(ptrInArg(U)), tuple(V), null );
00207 }
00208 
00209 
00210 template<class Scalar>
00211 void Thyra::update( Scalar alpha, const MultiVectorBase<Scalar>& U,
00212   const Ptr<MultiVectorBase<Scalar> > &V )
00213 {
00214   using Teuchos::tuple; using Teuchos::null;
00215   RTOpPack::TOpAXPY<Scalar> axpy_op(alpha);
00216   applyOp<Scalar>( axpy_op, tuple(ptrInArg(U)), tuple(V), null );
00217 }
00218 
00219 
00220 template<class Scalar>
00221 void Thyra::update( const ArrayView<const Scalar> &alpha, Scalar beta,
00222   const MultiVectorBase<Scalar>& U, const Ptr<MultiVectorBase<Scalar> > &V )
00223 {
00224 #ifdef TEUCHOS_DEBUG
00225   bool is_compatible = U.range()->isCompatible(*V->range());
00226   TEST_FOR_EXCEPTION(
00227     !is_compatible, Exceptions::IncompatibleVectorSpaces,
00228     "update(...), Error, U.range()->isCompatible((V->range())==false");
00229   is_compatible = U.domain()->isCompatible(*V->domain());
00230   TEST_FOR_EXCEPTION(
00231     !is_compatible, Exceptions::IncompatibleVectorSpaces,
00232     "update(...), Error, U.domain().isCompatible(V->domain())==false");
00233 #endif
00234   const int m = U.domain()->dim();
00235   for( int j = 0; j < m; ++j )
00236     Vp_StV<Scalar>( V->col(j).ptr(), alpha[j]*beta, *U.col(j) );
00237 }
00238 
00239 
00240 template<class Scalar>
00241 void Thyra::update( const MultiVectorBase<Scalar>& U,
00242   const ArrayView<const Scalar> &alpha, Scalar beta,
00243   const Ptr<MultiVectorBase<Scalar> > &V )
00244 {
00245 #ifdef TEUCHOS_DEBUG
00246   bool is_compatible = U.range()->isCompatible(*V->range());
00247     TEST_FOR_EXCEPTION(
00248       !is_compatible, Exceptions::IncompatibleVectorSpaces,
00249       "update(...), Error, U.range()->isCompatible((V->range())==false");
00250     is_compatible = U.domain()->isCompatible(*V->domain());
00251     TEST_FOR_EXCEPTION(
00252       !is_compatible, Exceptions::IncompatibleVectorSpaces,
00253       "update(...), Error, U.domain().isCompatible(V->domain())==false");
00254 #endif
00255   const int m = U.domain()->dim();
00256   for( int j = 0; j < m; ++j ) {
00257     Vt_S<Scalar>( V->col(j).ptr(), alpha[j]*beta );
00258     Vp_StV<Scalar>( V->col(j).ptr(), 1.0, *U.col(j) );
00259   }
00260 }
00261 
00262 
00263 template<class Scalar>
00264 void Thyra::linear_combination(
00265   const ArrayView<const Scalar> &alpha,
00266   const ArrayView<const Ptr<const MultiVectorBase<Scalar> > > &X,
00267   const Scalar &beta,
00268   const Ptr<MultiVectorBase<Scalar> > &Y
00269   )
00270 {
00271   using Teuchos::tuple; using Teuchos::null;
00272 #ifdef TEUCHOS_DEBUG
00273   TEUCHOS_ASSERT_EQUALITY(alpha.size(),X.size());
00274 #endif
00275   const int m = alpha.size();
00276   if ( beta == ScalarTraits<Scalar>::one() && m == 1 ) {
00277     update( alpha[0], *X[0], Y );
00278     return;
00279   }
00280   else if (m == 0) {
00281     scale( beta, Y );
00282     return;
00283   }
00284   RTOpPack::TOpLinearCombination<Scalar> lin_comb_op(alpha, beta);
00285   Thyra::applyOp<Scalar>( lin_comb_op, X, tuple(Y), null );
00286 }
00287 
00288 
00289 template<class Scalar>
00290 void Thyra::randomize( Scalar l, Scalar u,
00291   const Ptr<MultiVectorBase<Scalar> > &V )
00292 {
00293   const int m = V->domain()->dim();
00294   for( int j = 0; j < m; ++j )
00295     randomize( l, u, V->col(j).ptr() );
00296   // Todo: call applyOp(...) directly!
00297 }
00298 
00299 
00300 template<class Scalar>
00301 void Thyra::Vt_S( const Ptr<MultiVectorBase<Scalar> > &Z,
00302   const Scalar& alpha )
00303 {
00304   const int m = Z->domain()->dim();
00305   for( int j = 0; j < m; ++j )
00306     Vt_S( Z->col(j).ptr(), alpha );
00307   // Todo: call applyOp(...) directly!
00308 }
00309 
00310 
00311 template<class Scalar>
00312 void Thyra::Vp_S( const Ptr<MultiVectorBase<Scalar> > &Z,
00313   const Scalar& alpha )
00314 {
00315   const int m = Z->domain()->dim();
00316   for( int j = 0; j < m; ++j )
00317     Vp_S( Z->col(j).ptr(), alpha );
00318   // Todo: call applyOp(...) directly!
00319 }
00320 
00321 
00322 template<class Scalar>
00323 void Thyra::Vp_V( const Ptr<MultiVectorBase<Scalar> > &Z,
00324   const MultiVectorBase<Scalar>& X )
00325 {
00326   using Teuchos::tuple; using Teuchos::ptrInArg;
00327   typedef Teuchos::ScalarTraits<Scalar> ST;
00328   linear_combination<Scalar>( tuple(ST::one()), tuple(ptrInArg(X)),
00329     ST::one(), Z );
00330 }
00331 
00332 
00333 template<class Scalar>
00334 void Thyra::V_VpV( const Ptr<MultiVectorBase<Scalar> > &Z,
00335   const MultiVectorBase<Scalar>& X, const MultiVectorBase<Scalar>& Y )
00336 {
00337   using Teuchos::tuple; using Teuchos::ptrInArg;
00338   typedef Teuchos::ScalarTraits<Scalar> ST;
00339   linear_combination<Scalar>(
00340     tuple(ST::one(), ST::one()), tuple(ptrInArg(X), ptrInArg(Y)),
00341     ST::zero(), Z
00342     );
00343 }
00344 
00345 
00346 template<class Scalar>
00347 void Thyra::V_VmV( const Ptr<MultiVectorBase<Scalar> > &Z,
00348   const MultiVectorBase<Scalar>& X, const MultiVectorBase<Scalar>& Y )
00349 {
00350   using Teuchos::tuple; using Teuchos::ptrInArg; using Teuchos::as;
00351   typedef Teuchos::ScalarTraits<Scalar> ST;
00352   linear_combination<Scalar>(
00353     tuple(ST::one(), as<Scalar>(-ST::one())), tuple(ptrInArg(X), ptrInArg(Y)),
00354     ST::zero(), Z
00355     );
00356 }
00357 
00358 
00359 template<class Scalar>
00360 void Thyra::V_StVpV(
00361   const Ptr<MultiVectorBase<Scalar> > &Z, const Scalar &alpha,
00362   const MultiVectorBase<Scalar>& X, const MultiVectorBase<Scalar>& Y 
00363   )
00364 {
00365   using Teuchos::tuple; using Teuchos::ptrInArg;
00366   typedef Teuchos::ScalarTraits<Scalar> ST;
00367   linear_combination<Scalar>(
00368     tuple(alpha, ST::one()),  tuple(ptrInArg(X), ptrInArg(Y)),
00369     ST::zero(), Z
00370     );
00371 }
00372 
00373 
00374 #endif // THYRA_MULTI_VECTOR_STD_OPS_HPP

Generated on Wed May 12 21:42:28 2010 for Thyra Operator/Vector Support by  doxygen 1.4.7