Thyra Version of the Day
Thyra_MultiVectorBase_def.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 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Roscoe A. Bartlett (bartlettra@ornl.gov) 
00038 // 
00039 // ***********************************************************************
00040 // @HEADER
00041 
00042 #ifndef THYRA_MULTI_VECTOR_BASE_HPP
00043 #define THYRA_MULTI_VECTOR_BASE_HPP
00044 
00045 #include "Thyra_MultiVectorBase_decl.hpp"
00046 #include "Thyra_LinearOpBase.hpp"
00047 #include "Thyra_VectorSpaceBase.hpp"
00048 
00049 #include "Thyra_VectorBase.hpp"
00050 #include "Thyra_VectorStdOps_decl.hpp"
00051 
00052 #include "RTOpPack_TOpAbs.hpp"
00053 #include "RTOpPack_ROpNorm1.hpp"
00054 
00055 namespace Thyra {
00056 
00057 
00058 // Provide access to the columns as VectorBase objects
00059 
00060 
00061 template<class Scalar>
00062 RCP<const VectorBase<Scalar> >
00063 MultiVectorBase<Scalar>::colImpl(Ordinal j) const
00064 {
00065   return const_cast<MultiVectorBase*>(this)->nonconstColImpl(j);
00066 }
00067 
00068 
00069 // Overridden methods from LinearOpBase
00070 
00071 
00072 template<class Scalar>
00073 RCP<const LinearOpBase<Scalar> >
00074 MultiVectorBase<Scalar>::clone() const
00075 {
00076   return this->clone_mv();
00077 }
00078 
00079 // Overridden methods from RowStatLinearOpBase
00080 
00081 template<class Scalar>
00082 bool MultiVectorBase<Scalar>::
00083 rowStatIsSupportedImpl(const RowStatLinearOpBaseUtils::ERowStat rowStat) const
00084 {
00085   switch (rowStat) {
00086     case RowStatLinearOpBaseUtils::ROW_STAT_INV_ROW_SUM:
00087     case RowStatLinearOpBaseUtils::ROW_STAT_ROW_SUM:
00088     case RowStatLinearOpBaseUtils::ROW_STAT_INV_COL_SUM:
00089     case RowStatLinearOpBaseUtils::ROW_STAT_COL_SUM:
00090       return true;
00091       break;
00092     default:
00093       TEUCHOS_TEST_FOR_EXCEPT(true);
00094   }
00095 
00096   return false; // will never be called
00097 }
00098 
00099 template<class Scalar>
00100 void MultiVectorBase<Scalar>::
00101 getRowStatImpl(const RowStatLinearOpBaseUtils::ERowStat rowStat,
00102                const Ptr<VectorBase<Scalar> > &rowStatVec) const
00103 {
00104   switch (rowStat) {
00105     case RowStatLinearOpBaseUtils::ROW_STAT_INV_ROW_SUM:
00106       absRowSum(rowStatVec);
00107       ::Thyra::reciprocal<Scalar>(*rowStatVec,rowStatVec.ptr());
00108       break;
00109     case RowStatLinearOpBaseUtils::ROW_STAT_ROW_SUM:
00110       // compute absolute row sum
00111       absRowSum(rowStatVec);
00112       break;
00113     case RowStatLinearOpBaseUtils::ROW_STAT_INV_COL_SUM:
00114       absColSum(rowStatVec);
00115       ::Thyra::reciprocal<Scalar>(*rowStatVec,rowStatVec.ptr());
00116       break;
00117     case RowStatLinearOpBaseUtils::ROW_STAT_COL_SUM:
00118       // compute absolute row sum
00119       absColSum(rowStatVec);
00120       break;
00121     default:
00122       TEUCHOS_TEST_FOR_EXCEPT(true);
00123   }
00124 }
00125 
00126 // Overridden methods from ScaledLinearOpBase
00127 
00128 template<class Scalar>
00129 bool MultiVectorBase<Scalar>::
00130 supportsScaleLeftImpl() const
00131 {
00132   return true;
00133 }
00134   
00135 template<class Scalar>
00136 bool MultiVectorBase<Scalar>::
00137 supportsScaleRightImpl() const
00138 {
00139   return true;
00140 }
00141   
00142 template<class Scalar>
00143 void MultiVectorBase<Scalar>::
00144 scaleLeftImpl(const VectorBase< Scalar > &row_scaling)
00145 {
00146   // loop over each column applying the row scaling
00147   for(Ordinal i=0;i<this->domain()->dim();i++)
00148     ::Thyra::ele_wise_scale<Scalar>(row_scaling,this->col(i).ptr());
00149 }
00150   
00151 template<class Scalar>
00152 void MultiVectorBase<Scalar>::
00153 scaleRightImpl(const VectorBase< Scalar > &col_scaling)
00154 {
00155   // this is probably incorrect if the domain is distrbuted
00156   // but if it is on every processor its probably fine...
00157 
00158   RTOpPack::SubVectorView<Scalar> view;
00159   col_scaling.acquireDetachedView(Thyra::Range1D(),&view);
00160 
00161   Teuchos::ArrayRCP<const Scalar> col_scaling_vec = view.values();
00162 
00163   // check to make sure things match up
00164   TEUCHOS_ASSERT(this->domain()->dim()==col_scaling_vec.size());
00165 
00166   for(Ordinal i=0;i<this->domain()->dim();i++)
00167     ::Thyra::scale<Scalar>(col_scaling_vec[i],this->col(i).ptr());
00168 }
00169 
00170 // helper methods
00171 
00172 template<class Scalar>
00173 void MultiVectorBase<Scalar>::
00174 absRowSum(const Teuchos::Ptr<Thyra::VectorBase<Scalar> > & output) const
00175 {
00176   using Teuchos::RCP;
00177   using Teuchos::ptrFromRef;
00178   using Teuchos::tuple;
00179 
00180   // compute absolute value of multi-vector
00181   RTOpPack::TOpAbs<Scalar> abs_op;
00182   RCP<MultiVectorBase<Scalar> > abs_mv = createMembers(this->range(),this->domain());
00183   ::Thyra::applyOp<Scalar>( abs_op, tuple(ptrFromRef(*this)), tuple(abs_mv.ptr()), Teuchos::null );
00184 
00185   // compute sum over all rows
00186   RCP<VectorBase<Scalar> > ones = Thyra::createMember(this->domain());
00187   ::Thyra::put_scalar<Scalar>(Teuchos::ScalarTraits<Scalar>::one(),ones.ptr());
00188   ::Thyra::apply<Scalar>(*abs_mv,Thyra::NOTRANS,*ones,output);
00189 }
00190 
00191 template<class Scalar>
00192 void MultiVectorBase<Scalar>::
00193 absColSum(const Teuchos::Ptr<Thyra::VectorBase<Scalar> > & output) const
00194 { 
00195   using Teuchos::tuple; 
00196   using Teuchos::ptrInArg; 
00197   using Teuchos::null;
00198   using Teuchos::Array;
00199   using Teuchos::ArrayView;
00200 
00201   RTOpPack::SubVectorView<Scalar> view;
00202   output->acquireDetachedView(Thyra::Range1D(),&view);
00203 
00204   // Thyra::norms_1<Scalar>(*this,view.values()());
00205 
00206   ArrayView<Scalar> norms = view.values()();
00207   RTOpPack::ROpNorm1<Scalar> op;
00208 
00209   const int m = this->domain()->dim();
00210   Array<RCP<RTOpPack::ReductTarget> > rcp_op_targs(m);
00211   Array<Ptr<RTOpPack::ReductTarget> > op_targs(m);
00212   for( int kc = 0; kc < m; ++kc ) {
00213     rcp_op_targs[kc] = op.reduct_obj_create();
00214     op_targs[kc] = rcp_op_targs[kc].ptr();
00215   }
00216   ::Thyra::applyOp<Scalar>(op, tuple(ptrInArg(*this)),
00217     ArrayView<Ptr<MultiVectorBase<Scalar> > >(null),
00218     op_targs );
00219   for( int kc = 0; kc < m; ++kc ) {
00220     norms[kc] = op(*op_targs[kc]);
00221   }
00222   
00223   output->commitDetachedView(&view);
00224 }
00225 
00226 
00227 } // end namespace Thyra
00228 
00229 
00230 #endif // THYRA_MULTI_VECTOR_BASE_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines