00001
00002
00003
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
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