PlayaSimpleBlockOpImpl.hpp

00001 /* @HEADER@ */
00002 //   
00003  /* @HEADER@ */
00004 
00005 #ifndef PLAYA_SIMPLEBLOCKOP_IMPL_HPP
00006 #define PLAYA_SIMPLEBLOCKOP_IMPL_HPP
00007 
00008 #include "PlayaTabs.hpp"
00009 #include "PlayaExceptions.hpp"
00010 #include "PlayaSimpleBlockOpDecl.hpp"
00011 #include "PlayaSimpleZeroOpDecl.hpp"
00012 
00013 #ifndef HAVE_TEUCHOS_EXPLICIT_INSTANTIATION
00014 #include "PlayaLinearOperatorImpl.hpp"
00015 #include "PlayaSimpleZeroOpImpl.hpp"
00016 #endif
00017 
00018 
00019 
00020 namespace Playa
00021 {
00022 using namespace Teuchos;
00023 
00024 
00025 
00026 /* ---- Simplified linear op with spaces ------- */
00027 
00028 template <class Scalar> inline
00029 SimpleBlockOp<Scalar>::SimpleBlockOp(
00030   const VectorSpace<Scalar>& domain,
00031   const VectorSpace<Scalar>& range)
00032   : LinearOpWithSpaces<Scalar>(domain, range), blocks_(range.numBlocks())
00033 {
00034   for (int i=0; i<blocks_.size(); i++) 
00035   {
00036     blocks_[i] = Array<LinearOperator<Scalar> >(domain.numBlocks());
00037     for (int j=0; j<blocks_[i].size(); j++)
00038     {
00039       blocks_[i][j] = zeroOperator(domain.getBlock(j), range.getBlock(i));
00040     }
00041   }
00042 }
00043 
00044 template <class Scalar> inline
00045 int SimpleBlockOp<Scalar>::numBlockRows() const
00046 {
00047   return blocks_.size();
00048 }
00049 
00050 template <class Scalar> inline
00051 int SimpleBlockOp<Scalar>::numBlockCols() const
00052 {
00053   return blocks_[0].size();
00054 }
00055 
00056 template <class Scalar> inline
00057 const LinearOperator<Scalar>& SimpleBlockOp<Scalar>::getBlock(int i, int j) const 
00058 {
00059   return blocks_[i][j];
00060 }
00061 
00062 template <class Scalar> inline
00063 LinearOperator<Scalar> SimpleBlockOp<Scalar>::getNonconstBlock(int i, int j) 
00064 {
00065   return blocks_[i][j];
00066 }
00067 
00068 
00069 template <class Scalar> inline
00070 void SimpleBlockOp<Scalar>::setBlock(int i, int j, 
00071   const LinearOperator<Scalar>& Aij) 
00072 {
00073   blocks_[i][j] = Aij;
00074 }
00075 
00076 template <class Scalar> inline
00077 void SimpleBlockOp<Scalar>::apply(Teuchos::ETransp transApplyType,
00078   const Vector<Scalar>& in,
00079   Vector<Scalar> out) const
00080 {
00081   Tabs tab(0);
00082   PLAYA_MSG2(this->verb(), tab << "SimpleBlockOp::apply()");
00083   if (transApplyType == Teuchos::NO_TRANS)
00084   {
00085     out.zero();
00086     for (int i=0; i<this->numBlockRows(); i++)
00087     {
00088       for (int j=0; j<this->numBlockCols(); j++)
00089       {
00090         Vector<Scalar> tmp; 
00091         blocks_[i][j].apply(in.getBlock(j), tmp);
00092         out.getNonConstBlock(i).update(1.0, tmp);
00093       }
00094     }
00095   }
00096 
00097   else if (transApplyType == Teuchos::TRANS)
00098   {
00099     for (int i=0; i<this->numBlockCols(); i++)
00100     {
00101       out.zero();
00102       for (int j=0; j<this->numBlockRows(); j++)
00103       {
00104         Vector<Scalar> tmp;
00105         blocks_[j][i].applyTranspose(in.getBlock(j),tmp);
00106         out.getNonConstBlock(i).update(1.0, tmp);
00107       }
00108     }
00109   }
00110   else
00111   {
00112     TEUCHOS_TEST_FOR_EXCEPT(transApplyType != Teuchos::TRANS && transApplyType != Teuchos::NO_TRANS);
00113   }
00114   PLAYA_MSG2(this->verb(), tab << "done SimpleBlockOp::apply()");
00115 }
00116 
00117 
00118 template <class Scalar> inline
00119 LinearOperator<Scalar> makeBlockOperator(
00120   const VectorSpace<Scalar>& domain,
00121   const VectorSpace<Scalar>& range
00122   )
00123 {
00124   RCP<SimpleBlockOp<Scalar> > b = 
00125     rcp(new SimpleBlockOp<Scalar>(domain, range));
00126   RCP<LinearOperatorBase<Scalar> > p = b;
00127   return p;
00128 }
00129 
00130 
00131 
00132 }
00133 
00134 #endif
doxygen