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