Thyra_DefaultBlockedLinearOp.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_DEFAULT_BLOCKED_LINEAR_OP_HPP
00030 #define THYRA_DEFAULT_BLOCKED_LINEAR_OP_HPP
00031 
00032 #include "Thyra_DefaultBlockedLinearOpDecl.hpp"
00033 #include "Thyra_DefaultProductVectorSpace.hpp"
00034 #include "Thyra_MultiVectorStdOps.hpp"
00035 #include "Thyra_AssertOp.hpp"
00036 #include "Teuchos_Utils.hpp"
00037 
00038 namespace Thyra {
00039 
00040 // Constructors
00041 
00042 template<class Scalar>
00043 DefaultBlockedLinearOp<Scalar>::DefaultBlockedLinearOp()
00044   :blockFillIsActive_(false)
00045 {}
00046 
00047 // Overridden from PhysicallyBlockedLinearOpBase
00048 
00049 template<class Scalar>
00050 void DefaultBlockedLinearOp<Scalar>::beginBlockFill()
00051 {
00052   assertBlockFillIsActive(false);
00053   resetStorage(0,0);
00054 }
00055 
00056 template<class Scalar>
00057 void DefaultBlockedLinearOp<Scalar>::beginBlockFill(
00058   const int numRowBlocks, const int numColBlocks
00059   )
00060 {
00061   assertBlockFillIsActive(false);
00062   resetStorage(numRowBlocks,numColBlocks);
00063 }
00064 
00065 template<class Scalar>
00066 void DefaultBlockedLinearOp<Scalar>::beginBlockFill(
00067   const Teuchos::RCP<const ProductVectorSpaceBase<Scalar> >  &productRange
00068   ,const Teuchos::RCP<const ProductVectorSpaceBase<Scalar> > &productDomain
00069   )
00070 {
00071   assertBlockFillIsActive(false);
00072   productRange_ = productRange.assert_not_null();
00073   productDomain_ = productDomain.assert_not_null();
00074   // Note: the above spaces must be set before calling the next function!
00075   resetStorage(productRange_->numBlocks(),productDomain_->numBlocks());
00076 }
00077 
00078 template<class Scalar>
00079 bool DefaultBlockedLinearOp<Scalar>::blockFillIsActive() const
00080 {
00081   return blockFillIsActive_;
00082 }
00083 
00084 template<class Scalar>
00085 bool DefaultBlockedLinearOp<Scalar>::acceptsBlock(
00086   const int i, const int j
00087   ) const
00088 {
00089   assertBlockFillIsActive(true);
00090   assertBlockRowCol(i,j);
00091   return true;
00092 }
00093 
00094 template<class Scalar>
00095 void DefaultBlockedLinearOp<Scalar>::setNonconstBlock(
00096   const int i, const int j
00097   ,const Teuchos::RCP<LinearOpBase<Scalar> > &block
00098   )
00099 {
00100   setBlockImpl(i,j,block);
00101 }
00102 
00103 template<class Scalar>
00104 void DefaultBlockedLinearOp<Scalar>::setBlock(
00105   const int i, const int j
00106   ,const Teuchos::RCP<const LinearOpBase<Scalar> > &block
00107   )
00108 {
00109   setBlockImpl(i,j,block);
00110 }
00111 
00112 template<class Scalar>
00113 void DefaultBlockedLinearOp<Scalar>::endBlockFill()
00114 {
00115   using Teuchos::rcp;
00116   using Teuchos::arrayArg;
00117   assertBlockFillIsActive(true);
00118   // Create the product range and domain spaces if these are not already set.
00119   if(!productRange_.get()) {
00120 #ifdef TEUCHOS_DEBUG
00121     for(int i = 0; i < numRowBlocks_; ++i) {
00122       TEST_FOR_EXCEPTION(
00123         !rangeBlocks_[i].get(), std::logic_error
00124         ,"DefaultBlockedLinearOp<Scalar>::endBlockFill():"
00125         " Error, no linear operator block for the i="<<i<<" block row was added"
00126         " and we can not complete the block fill!"
00127         );
00128     }
00129     for(int j = 0; j < numColBlocks_; ++j) {
00130       TEST_FOR_EXCEPTION(
00131         !domainBlocks_[j].get(), std::logic_error
00132         ,"DefaultBlockedLinearOp<Scalar>::endBlockFill():"
00133         " Error, no linear operator block for the j="
00134         <<j<<" block column was added"
00135         " and we can not complete the block fill!"
00136         );
00137     }
00138 #endif
00139     productRange_
00140       = rcp(
00141         new DefaultProductVectorSpace<Scalar>(
00142           numRowBlocks_,&rangeBlocks_[0]
00143           )
00144         );
00145     productDomain_
00146       = rcp(
00147         new DefaultProductVectorSpace<Scalar>(
00148           numColBlocks_,&domainBlocks_[0]
00149           )
00150         );
00151   }
00152   numRowBlocks_ = productRange_->numBlocks();
00153   numColBlocks_ = productDomain_->numBlocks();
00154   rangeBlocks_.resize(0);
00155   domainBlocks_.resize(0);
00156   // Insert the block LOB objects if doing a flexible fill.
00157   if(Ops_stack_.size()) {
00158     Ops_.resize(numRowBlocks_*numColBlocks_);
00159     for( int k = 0; k < static_cast<int>(Ops_stack_.size()); ++k ) {
00160       const BlockEntry<Scalar> &block_i_j = Ops_stack_[k];
00161       Ops_[numRowBlocks_*block_i_j.j+block_i_j.i] = block_i_j.block;
00162     }
00163     Ops_stack_.resize(0);
00164   }
00165   blockFillIsActive_ = false;
00166 }
00167 
00168 template<class Scalar>
00169 void DefaultBlockedLinearOp<Scalar>::uninitialize()
00170 {
00171   productRange_ = Teuchos::null;
00172   productDomain_ = Teuchos::null;
00173   numRowBlocks_ = 0;
00174   numColBlocks_ = 0;
00175   Ops_.resize(0);
00176   Ops_stack_.resize(0);
00177   rangeBlocks_.resize(0);
00178   domainBlocks_.resize(0);
00179   blockFillIsActive_ = false;
00180 }
00181 
00182 // Overridden from BlockedLinearOpBase
00183 
00184 template<class Scalar>
00185 Teuchos::RCP<const ProductVectorSpaceBase<Scalar> >
00186 DefaultBlockedLinearOp<Scalar>::productRange() const
00187 {
00188   return productRange_;
00189 }
00190 
00191 template<class Scalar>
00192 Teuchos::RCP<const ProductVectorSpaceBase<Scalar> >
00193 DefaultBlockedLinearOp<Scalar>::productDomain() const
00194 {
00195   return productDomain_;
00196 }
00197 
00198 template<class Scalar>
00199 bool DefaultBlockedLinearOp<Scalar>::blockExists(
00200   const int i, const int j
00201   ) const
00202 {
00203   assertBlockFillIsActive(false);
00204   assertBlockRowCol(i,j);
00205   return true;
00206 } 
00207 
00208 template<class Scalar>
00209 bool DefaultBlockedLinearOp<Scalar>::blockIsConst(
00210   const int i, const int j
00211   ) const
00212 {
00213 #ifdef TEUCHOS_DEBUG
00214   TEST_FOR_EXCEPT(!blockExists(i,j));
00215 #endif
00216   assertBlockFillIsActive(false);
00217   assertBlockRowCol(i,j);
00218   return Ops_[numRowBlocks_*j+i].isConst();
00219 }
00220 
00221 template<class Scalar>
00222 Teuchos::RCP<LinearOpBase<Scalar> >
00223 DefaultBlockedLinearOp<Scalar>::getNonconstBlock(const int i, const int j)
00224 {
00225 #ifdef TEUCHOS_DEBUG
00226   TEST_FOR_EXCEPT(!blockExists(i,j));
00227 #endif
00228   assertBlockFillIsActive(false);
00229   assertBlockRowCol(i,j);
00230   return Ops_[numRowBlocks_*j+i].getNonconstObj();
00231 } 
00232 
00233 template<class Scalar>
00234 Teuchos::RCP<const LinearOpBase<Scalar> >
00235 DefaultBlockedLinearOp<Scalar>::getBlock(const int i, const int j) const
00236 {
00237 #ifdef TEUCHOS_DEBUG
00238   TEST_FOR_EXCEPT(!blockExists(i,j));
00239 #endif
00240   assertBlockFillIsActive(false);
00241   assertBlockRowCol(i,j);
00242   return Ops_[numRowBlocks_*j+i].getConstObj();
00243 } 
00244 
00245 // Overridden from LinearOpBase
00246 
00247 template<class Scalar>
00248 Teuchos::RCP< const VectorSpaceBase<Scalar> >
00249 DefaultBlockedLinearOp<Scalar>::range() const
00250 {
00251   return productRange_;
00252 }
00253 
00254 template<class Scalar>
00255 Teuchos::RCP< const VectorSpaceBase<Scalar> >
00256 DefaultBlockedLinearOp<Scalar>::domain() const
00257 {
00258   return productDomain_;
00259 }
00260 
00261 template<class Scalar>
00262 Teuchos::RCP<const LinearOpBase<Scalar> >
00263 DefaultBlockedLinearOp<Scalar>::clone() const
00264 {
00265   return Teuchos::null; // ToDo: Implement this when needed!
00266 }
00267 
00268 // Overridden from Teuchos::Describable
00269 
00270 template<class Scalar>
00271 std::string DefaultBlockedLinearOp<Scalar>::description() const
00272 {
00273   assertBlockFillIsActive(false);
00274   std::ostringstream oss;
00275   oss
00276     << Teuchos::Describable::description() << "{"
00277     << "numRowBlocks="<<numRowBlocks_
00278     << ",numColBlocks="<<numColBlocks_
00279     << "}";
00280   return oss.str();
00281 }
00282 
00283 template<class Scalar>
00284 void DefaultBlockedLinearOp<Scalar>::describe(
00285   Teuchos::FancyOStream                &out_arg
00286   ,const Teuchos::EVerbosityLevel      verbLevel
00287   ) const
00288 {
00289   typedef Teuchos::ScalarTraits<Scalar>  ST;
00290   using Teuchos::RCP;
00291   using Teuchos::FancyOStream;
00292   using Teuchos::OSTab;
00293   assertBlockFillIsActive(false);
00294   RCP<FancyOStream> out = rcp(&out_arg,false);
00295   OSTab tab(out);
00296   switch(verbLevel) {
00297     case Teuchos::VERB_DEFAULT:
00298     case Teuchos::VERB_LOW:
00299       *out << this->description() << std::endl;
00300       break;
00301     case Teuchos::VERB_MEDIUM:
00302     case Teuchos::VERB_HIGH:
00303     case Teuchos::VERB_EXTREME:
00304     {
00305       *out
00306         << Teuchos::Describable::description() << "{"
00307         << "rangeDim=" << this->range()->dim()
00308         << ",domainDim=" << this->domain()->dim()
00309         << ",numRowBlocks=" << numRowBlocks_
00310         << ",numColBlocks=" << numColBlocks_
00311         << "}\n";
00312       OSTab tab(out);
00313       *out
00314         <<  "Constituent LinearOpBase objects for M = [ Op[0,0] ..."
00315   << " ; ... ; ... Op[numRowBlocks-1,numColBlocks-1]:\n";
00316       tab.incrTab();
00317       for( int i = 0; i < numRowBlocks_; ++i ) {
00318         for( int j = 0; j < numColBlocks_; ++j ) {
00319           *out << "Op["<<i<<","<<j<<"] = ";
00320           Teuchos::RCP<const LinearOpBase<Scalar> >
00321             block_i_j = getBlock(i,j);
00322           if(block_i_j.get())
00323             *out << Teuchos::describe(*getBlock(i,j),verbLevel);
00324           else
00325             *out << "NULL\n";
00326         }
00327       }
00328       break;
00329     }
00330     default:
00331       TEST_FOR_EXCEPT(true); // Should never get here!
00332   }
00333 }
00334 
00335 // protected
00336 
00337 // Overridden from SingleScalarLinearOpBase
00338 
00339 template<class Scalar>
00340 bool DefaultBlockedLinearOp<Scalar>::opSupported(
00341   ETransp M_trans
00342   ) const
00343 {
00344   bool opSupported = true;
00345   for( int i = 0; i < numRowBlocks_; ++i ) {
00346     for( int j = 0; j < numColBlocks_; ++j ) {
00347       Teuchos::RCP<const LinearOpBase<Scalar> >
00348         block_i_j = getBlock(i,j);
00349       if( block_i_j.get() && !Thyra::opSupported(*block_i_j,M_trans) )
00350         opSupported = false;
00351     }
00352   }
00353   return opSupported;
00354 }
00355 
00356 template<class Scalar>
00357 void DefaultBlockedLinearOp<Scalar>::apply(
00358   const ETransp                     M_trans
00359   ,const MultiVectorBase<Scalar>    &X_in
00360   ,MultiVectorBase<Scalar>          *Y_inout
00361   ,const Scalar                     alpha
00362   ,const Scalar                     beta
00363   ) const
00364 {
00365   using Teuchos::RCP;
00366   using Teuchos::dyn_cast;
00367   typedef Teuchos::ScalarTraits<Scalar>                ST;
00368   typedef RCP<MultiVectorBase<Scalar> >        MultiVectorPtr;
00369   typedef RCP<const MultiVectorBase<Scalar> >  ConstMultiVectorPtr;
00370   typedef RCP<const LinearOpBase<Scalar> >     ConstLinearOpPtr;
00371 #ifdef TEUCHOS_DEBUG
00372   THYRA_ASSERT_LINEAR_OP_MULTIVEC_APPLY_SPACES(
00373     "DefaultBlockedLinearOp<Scalar>::apply(...)",*this,M_trans,X_in,Y_inout
00374     );
00375 #endif // TEUCHOS_DEBUG  
00376   const bool
00377     struct_transp = (real_trans(M_trans)!=NOTRANS); // Structural transpose?
00378   const int
00379     opNumRowBlocks = ( !struct_transp ? numRowBlocks_ : numColBlocks_ ), 
00380     opNumColBlocks = ( !struct_transp ? numColBlocks_ : numRowBlocks_ ); 
00381   //
00382   // Y = alpha * op(M) * X + beta*Y
00383   //
00384   // =>
00385   //
00386   // Y[i] = beta+Y[i] + sum(alpha*op(Op)[i,j]*X[j],j=0...opNumColBlocks-1)
00387   //
00388   //   , for i=0...opNumRowBlocks-1
00389   //
00390   const ProductMultiVectorBase<Scalar>
00391     &X = dyn_cast<const ProductMultiVectorBase<Scalar> >(X_in);
00392   ProductMultiVectorBase<Scalar>
00393     &Y = dyn_cast<ProductMultiVectorBase<Scalar> >(*Y_inout);
00394   for( int i = 0; i < opNumRowBlocks; ++i ) {
00395     MultiVectorPtr Y_i = Y.getNonconstMultiVectorBlock(i);
00396     for( int j = 0; j < opNumColBlocks; ++j ) {
00397       ConstLinearOpPtr
00398         Op_i_j = ( !struct_transp ? getBlock(i,j) : getBlock(j,i) );
00399       ConstMultiVectorPtr
00400         X_j    = X.getMultiVectorBlock(j);
00401       if(j==0) {
00402         if(Op_i_j.get())
00403           Thyra::apply(*Op_i_j,M_trans,*X_j,&*Y_i,alpha,beta);
00404         else
00405           scale(beta,&*Y_i);
00406       }
00407       else {
00408         if(Op_i_j.get())
00409           Thyra::apply(*Op_i_j,M_trans,*X_j,&*Y_i,alpha,ST::one());
00410       }
00411     }
00412   }
00413 }
00414 
00415 // private
00416 
00417 template<class Scalar>
00418 void DefaultBlockedLinearOp<Scalar>::resetStorage(
00419   const int numRowBlocks, const int numColBlocks
00420   )
00421 {
00422   uninitialize();
00423   numRowBlocks_ = numRowBlocks;
00424   numColBlocks_ = numColBlocks;
00425   Ops_.resize(numRowBlocks_*numColBlocks_);
00426   if(!productRange_.get()) {
00427     rangeBlocks_.resize(numRowBlocks);
00428     domainBlocks_.resize(numColBlocks);
00429   }
00430   blockFillIsActive_ = true;
00431 }
00432 
00433 template<class Scalar>
00434 void DefaultBlockedLinearOp<Scalar>::assertBlockFillIsActive(
00435   bool blockFillIsActive
00436   ) const
00437 {
00438 #ifdef TEUCHOS_DEBUG
00439   TEST_FOR_EXCEPT(!(blockFillIsActive_==blockFillIsActive));
00440 #endif
00441 }
00442 
00443 template<class Scalar>
00444 void DefaultBlockedLinearOp<Scalar>::assertBlockRowCol(
00445   const int i, const int j
00446   ) const
00447 {
00448 #ifdef TEUCHOS_DEBUG
00449   TEST_FOR_EXCEPTION(
00450     !( 0 <= i ), std::logic_error
00451     ,"Error, i="<<i<<" is invalid!"
00452     );
00453   TEST_FOR_EXCEPTION(
00454     !( 0 <= j ), std::logic_error
00455     ,"Error, j="<<j<<" is invalid!"
00456     );
00457   // Only validate upper range if the number of row and column blocks is
00458   // fixed!
00459   if(Ops_.size()) {
00460     TEST_FOR_EXCEPTION(
00461       !( 0 <= i && i < numRowBlocks_ ), std::logic_error
00462       ,"Error, i="<<i<<" does not fall in the range [0,"<<numRowBlocks_-1<<"]!"
00463       );
00464     TEST_FOR_EXCEPTION(
00465       !( 0 <= j && j < numColBlocks_ ), std::logic_error
00466       ,"Error, j="<<j<<" does not fall in the range [0,"<<numColBlocks_-1<<"]!"
00467       );
00468   }
00469 #endif
00470 }
00471 
00472 template<class Scalar>
00473 void DefaultBlockedLinearOp<Scalar>::setBlockSpaces(
00474   const int i, const int j, const LinearOpBase<Scalar> &block
00475   )
00476 {
00477   typedef std::string s;
00478   using Teuchos::toString;
00479   assertBlockFillIsActive(true);
00480   assertBlockRowCol(i,j);
00481   // Validate that if the vector space block is already set that it is
00482   // compatible with the block that is being set.
00483   if( i < numRowBlocks_ && j < numColBlocks_ ) {
00484 #ifdef TEUCHOS_DEBUG
00485     Teuchos::RCP<const VectorSpaceBase<Scalar> >
00486       rangeBlock = (
00487         productRange_.get()
00488         ? productRange_->getBlock(i)
00489         : rangeBlocks_[i]
00490         ),
00491       domainBlock = (
00492         productDomain_.get()
00493         ? productDomain_->getBlock(j)
00494         : domainBlocks_[j]
00495         );
00496     if(rangeBlock.get()) {
00497       THYRA_ASSERT_VEC_SPACES_NAMES(
00498         "DefaultBlockedLinearOp<Scalar>::setBlockSpaces(i,j,block):\n\n"
00499         "Adding block: " + block.description(),
00500         *rangeBlock,("(*productRange->getBlock("+toString(i)+"))"),
00501         *block.range(),("(*block["+toString(i)+","+toString(j)+"].range())")
00502         );
00503     }
00504     if(domainBlock.get()) {
00505       THYRA_ASSERT_VEC_SPACES_NAMES(
00506         "DefaultBlockedLinearOp<Scalar>::setBlockSpaces(i,j,block):\n\n"
00507         "Adding block: " + block.description(),
00508         *domainBlock,("(*productDomain->getBlock("+toString(j)+"))"),
00509         *block.domain(),("(*block["+toString(i)+","+toString(j)+"].domain())")
00510         );
00511     }
00512 #endif // TEUCHOS_DEBUG
00513   }
00514   // Add spaces missing range and domain space blocks if we are doing a
00515   // flexible fill (otherwise these loops will not be executed)
00516   for( int k = numRowBlocks_; k <= i; ++k )
00517     rangeBlocks_.push_back(Teuchos::null);
00518   for( int k = numColBlocks_; k <= j; ++k )
00519     domainBlocks_.push_back(Teuchos::null);
00520   // Set the incoming range and domain blocks if not already set
00521   if(!productRange_.get()) {
00522     if(!rangeBlocks_[i].get())
00523       rangeBlocks_[i] = block.range();
00524     if(!domainBlocks_[j].get()) {
00525       domainBlocks_[j] = block.domain();
00526     }
00527   }
00528   // Update the current number of row and columns blocks if doing a flexible
00529   // fill.
00530   if(!Ops_.size()) {
00531     numRowBlocks_ = rangeBlocks_.size();
00532     numColBlocks_ = domainBlocks_.size();
00533   }
00534 }
00535 
00536 template<class Scalar>
00537 template<class LinearOpType>
00538 void DefaultBlockedLinearOp<Scalar>::setBlockImpl(
00539   const int i, const int j
00540   ,const Teuchos::RCP<LinearOpType> &block
00541   )
00542 {
00543   setBlockSpaces(i,j,*block);
00544   if(Ops_.size()) {
00545     // We are doing a fill with a fixed number of row and column blocks so we
00546     // can just set this.
00547     Ops_[numRowBlocks_*j+i] = block;
00548   }
00549   else {
00550     // We are doing a flexible fill so add the block to the stack of blocks or
00551     // replace a block that already exists.
00552     bool foundBlock = false;
00553     for( unsigned int k = 0; k < Ops_stack_.size(); ++k ) {
00554       BlockEntry<Scalar> &block_i_j = Ops_stack_[k];
00555       if( block_i_j.i == i && block_i_j.j == j ) {
00556         block_i_j.block = block;
00557         foundBlock = true;
00558         break;
00559       }
00560     }
00561     if(!foundBlock)
00562       Ops_stack_.push_back(BlockEntry<Scalar>(i,j,block));
00563   }
00564 }
00565 
00566 } // namespace Thyra
00567 
00568 template<class Scalar>
00569 Teuchos::RCP<const Thyra::LinearOpBase<Scalar> >
00570 Thyra::block2x2(
00571   const Teuchos::RCP<const LinearOpBase<Scalar> >    &A00
00572   ,const Teuchos::RCP<const LinearOpBase<Scalar> >   &A01
00573   ,const Teuchos::RCP<const LinearOpBase<Scalar> >   &A10
00574   ,const Teuchos::RCP<const LinearOpBase<Scalar> >   &A11
00575   )
00576 {
00577   Teuchos::RCP<PhysicallyBlockedLinearOpBase<Scalar> >
00578     M = Teuchos::rcp(new DefaultBlockedLinearOp<Scalar>());
00579   M->beginBlockFill(2,2);
00580   if(A00.get()) M->setBlock(0,0,A00);
00581   if(A01.get()) M->setBlock(0,1,A01);
00582   if(A10.get()) M->setBlock(1,0,A10);
00583   if(A11.get()) M->setBlock(1,1,A11);
00584   M->endBlockFill();
00585   return M;
00586 }
00587 
00588 template<class Scalar>
00589 Teuchos::RCP<const Thyra::LinearOpBase<Scalar> >
00590 Thyra::block2x1(
00591   const Teuchos::RCP<const LinearOpBase<Scalar> >    &A00
00592   ,const Teuchos::RCP<const LinearOpBase<Scalar> >   &A10
00593   )
00594 {
00595   Teuchos::RCP<PhysicallyBlockedLinearOpBase<Scalar> >
00596     M = Teuchos::rcp(new DefaultBlockedLinearOp<Scalar>());
00597   M->beginBlockFill(2,1);
00598   if(A00.get()) M->setBlock(0,0,A00);
00599   if(A10.get()) M->setBlock(1,0,A10);
00600   M->endBlockFill();
00601   return M;
00602 }
00603 
00604 template<class Scalar>
00605 Teuchos::RCP<const Thyra::LinearOpBase<Scalar> >
00606 Thyra::block1x2(
00607   const Teuchos::RCP<const LinearOpBase<Scalar> >    &A00
00608   ,const Teuchos::RCP<const LinearOpBase<Scalar> >   &A01
00609   )
00610 {
00611   Teuchos::RCP<PhysicallyBlockedLinearOpBase<Scalar> >
00612     M = Teuchos::rcp(new DefaultBlockedLinearOp<Scalar>());
00613   M->beginBlockFill(1,2);
00614   if(A00.get()) M->setBlock(0,0,A00);
00615   if(A01.get()) M->setBlock(0,1,A01);
00616   M->endBlockFill();
00617   return M;
00618 }
00619 
00620 template<class Scalar>
00621 Teuchos::RCP<Thyra::LinearOpBase<Scalar> >
00622 Thyra::nonconstBlock2x2(
00623   const Teuchos::RCP<LinearOpBase<Scalar> >    &A00
00624   ,const Teuchos::RCP<LinearOpBase<Scalar> >   &A01
00625   ,const Teuchos::RCP<LinearOpBase<Scalar> >   &A10
00626   ,const Teuchos::RCP<LinearOpBase<Scalar> >   &A11
00627   )
00628 {
00629   Teuchos::RCP<PhysicallyBlockedLinearOpBase<Scalar> >
00630     M = Teuchos::rcp(new DefaultBlockedLinearOp<Scalar>());
00631   M->beginBlockFill(2,2);
00632   if(A00.get()) M->setNonconstBlock(0,0,A00);
00633   if(A01.get()) M->setNonconstBlock(0,1,A01);
00634   if(A10.get()) M->setNonconstBlock(1,0,A10);
00635   if(A11.get()) M->setNonconstBlock(1,1,A11);
00636   M->endBlockFill();
00637   return M;
00638 }
00639 
00640 template<class Scalar>
00641 Teuchos::RCP<Thyra::LinearOpBase<Scalar> >
00642 Thyra::nonconstBlock2x1(
00643   const Teuchos::RCP<LinearOpBase<Scalar> >    &A00
00644   ,const Teuchos::RCP<LinearOpBase<Scalar> >   &A10
00645   )
00646 {
00647   Teuchos::RCP<PhysicallyBlockedLinearOpBase<Scalar> >
00648     M = Teuchos::rcp(new DefaultBlockedLinearOp<Scalar>());
00649   M->beginBlockFill(2,1);
00650   if(A00.get()) M->setNonconstBlock(0,0,A00);
00651   if(A10.get()) M->setNonconstBlock(1,0,A10);
00652   M->endBlockFill();
00653   return M;
00654 }
00655 
00656 template<class Scalar>
00657 Teuchos::RCP<Thyra::LinearOpBase<Scalar> >
00658 Thyra::nonconstBlock1x2(
00659   const Teuchos::RCP<LinearOpBase<Scalar> >    &A00
00660   ,const Teuchos::RCP<LinearOpBase<Scalar> >   &A01
00661   )
00662 {
00663   Teuchos::RCP<PhysicallyBlockedLinearOpBase<Scalar> >
00664     M = Teuchos::rcp(new DefaultBlockedLinearOp<Scalar>());
00665   M->beginBlockFill(1,2);
00666   if(A00.get()) M->setNonconstBlock(0,0,A00);
00667   if(A01.get()) M->setNonconstBlock(0,1,A01);
00668   M->endBlockFill();
00669   return M;
00670 }
00671 
00672 #endif  // THYRA_DEFAULT_BLOCKED_LINEAR_OP_HPP

Generated on Tue Oct 20 12:47:25 2009 for Thyra Package Browser (Single Doxygen Collection) by doxygen 1.4.7