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::RefCountPtr<const ProductVectorSpaceBase<Scalar> >  &productRange
00068   ,const Teuchos::RefCountPtr<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::RefCountPtr<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::RefCountPtr<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="<<j<<" block column was added"
00134         " and we can not complete the block fill!"
00135         );
00136     }
00137 #endif
00138     productRange_
00139       = rcp(
00140         new DefaultProductVectorSpace<Scalar>(
00141           numRowBlocks_,&rangeBlocks_[0]
00142           )
00143         );
00144     productDomain_
00145       = rcp(
00146         new DefaultProductVectorSpace<Scalar>(
00147           numColBlocks_,&domainBlocks_[0]
00148           )
00149         );
00150   }
00151   numRowBlocks_ = productRange_->numBlocks();
00152   numColBlocks_ = productDomain_->numBlocks();
00153   rangeBlocks_.resize(0);
00154   domainBlocks_.resize(0);
00155   // Insert the block LOB objects if doing a flexible fill.
00156   if(Ops_stack_.size()) {
00157     Ops_.resize(numRowBlocks_*numColBlocks_);
00158     for( int k = 0; k < static_cast<int>(Ops_stack_.size()); ++k ) {
00159       const BlockEntry<Scalar> &block_i_j = Ops_stack_[k];
00160       Ops_[numRowBlocks_*block_i_j.j+block_i_j.i] = block_i_j.block;
00161     }
00162     Ops_stack_.resize(0);
00163   }
00164   blockFillIsActive_ = false;
00165 }
00166 
00167 template<class Scalar>
00168 void DefaultBlockedLinearOp<Scalar>::uninitialize()
00169 {
00170   productRange_ = Teuchos::null;
00171   productDomain_ = Teuchos::null;
00172   numRowBlocks_ = 0;
00173   numColBlocks_ = 0;
00174   Ops_.resize(0);
00175   Ops_stack_.resize(0);
00176   rangeBlocks_.resize(0);
00177   domainBlocks_.resize(0);
00178   blockFillIsActive_ = false;
00179 }
00180 
00181 // Overridden from BlockedLinearOpBase
00182 
00183 template<class Scalar>
00184 Teuchos::RefCountPtr<const ProductVectorSpaceBase<Scalar> >
00185 DefaultBlockedLinearOp<Scalar>::productRange() const
00186 {
00187   return productRange_;
00188 }
00189 
00190 template<class Scalar>
00191 Teuchos::RefCountPtr<const ProductVectorSpaceBase<Scalar> >
00192 DefaultBlockedLinearOp<Scalar>::productDomain() const
00193 {
00194   return productDomain_;
00195 }
00196 
00197 template<class Scalar>
00198 bool DefaultBlockedLinearOp<Scalar>::blockExists(
00199   const int i, const int j
00200   ) const
00201 {
00202   assertBlockFillIsActive(false);
00203   assertBlockRowCol(i,j);
00204   return true;
00205 } 
00206 
00207 template<class Scalar>
00208 bool DefaultBlockedLinearOp<Scalar>::blockIsConst(
00209   const int i, const int j
00210   ) const
00211 {
00212 #ifdef TEUCHOS_DEBUG
00213   TEST_FOR_EXCEPT(!blockExists(i,j));
00214 #endif
00215   assertBlockFillIsActive(false);
00216   assertBlockRowCol(i,j);
00217   return Ops_[numRowBlocks_*j+i].isConst();
00218 }
00219 
00220 template<class Scalar>
00221 Teuchos::RefCountPtr<LinearOpBase<Scalar> >
00222 DefaultBlockedLinearOp<Scalar>::getNonconstBlock(const int i, const int j)
00223 {
00224 #ifdef TEUCHOS_DEBUG
00225   TEST_FOR_EXCEPT(!blockExists(i,j));
00226 #endif
00227   assertBlockFillIsActive(false);
00228   assertBlockRowCol(i,j);
00229   return Ops_[numRowBlocks_*j+i].getNonconstObj();
00230 } 
00231 
00232 template<class Scalar>
00233 Teuchos::RefCountPtr<const LinearOpBase<Scalar> >
00234 DefaultBlockedLinearOp<Scalar>::getBlock(const int i, const int j) const
00235 {
00236 #ifdef TEUCHOS_DEBUG
00237   TEST_FOR_EXCEPT(!blockExists(i,j));
00238 #endif
00239   assertBlockFillIsActive(false);
00240   assertBlockRowCol(i,j);
00241   return Ops_[numRowBlocks_*j+i].getConstObj();
00242 } 
00243 
00244 // Overridden from LinearOpBase
00245 
00246 template<class Scalar>
00247 Teuchos::RefCountPtr< const VectorSpaceBase<Scalar> >
00248 DefaultBlockedLinearOp<Scalar>::range() const
00249 {
00250   return productRange_;
00251 }
00252 
00253 template<class Scalar>
00254 Teuchos::RefCountPtr< const VectorSpaceBase<Scalar> >
00255 DefaultBlockedLinearOp<Scalar>::domain() const
00256 {
00257   return productDomain_;
00258 }
00259 
00260 template<class Scalar>
00261 Teuchos::RefCountPtr<const LinearOpBase<Scalar> >
00262 DefaultBlockedLinearOp<Scalar>::clone() const
00263 {
00264   return Teuchos::null; // ToDo: Implement this when needed!
00265 }
00266 
00267 // Overridden from Teuchos::Describable
00268 
00269 template<class Scalar>
00270 std::string DefaultBlockedLinearOp<Scalar>::description() const
00271 {
00272   typedef Teuchos::ScalarTraits<Scalar>  ST;
00273   assertBlockFillIsActive(false);
00274   std::ostringstream oss;
00275   oss
00276     << "Thyra::DefaultBlockedLinearOp<" << ST::name() << ">"
00277     << "{"
00278     << "numRowBlocks="<<numRowBlocks_
00279     << ",numColBlocks="<<numColBlocks_
00280     << "}";
00281   return oss.str();
00282 }
00283 
00284 template<class Scalar>
00285 void DefaultBlockedLinearOp<Scalar>::describe(
00286   Teuchos::FancyOStream                &out_arg
00287   ,const Teuchos::EVerbosityLevel      verbLevel
00288   ) const
00289 {
00290   typedef Teuchos::ScalarTraits<Scalar>  ST;
00291   using Teuchos::RefCountPtr;
00292   using Teuchos::FancyOStream;
00293   using Teuchos::OSTab;
00294   assertBlockFillIsActive(false);
00295   RefCountPtr<FancyOStream> out = rcp(&out_arg,false);
00296   OSTab tab(out);
00297   switch(verbLevel) {
00298     case Teuchos::VERB_DEFAULT:
00299     case Teuchos::VERB_LOW:
00300       *out << this->description() << std::endl;
00301       break;
00302     case Teuchos::VERB_MEDIUM:
00303     case Teuchos::VERB_HIGH:
00304     case Teuchos::VERB_EXTREME:
00305     {
00306       *out
00307         << "Thyra::DefaultBlockedLinearOp<" << ST::name() << ">{"
00308         << "rangeDim=" << this->range()->dim()
00309         << ",domainDim=" << this->domain()->dim()
00310         << ",numRowBlocks=" << numRowBlocks_
00311         << ",numColBlocks=" << numColBlocks_
00312         << "}\n";
00313       OSTab tab(out);
00314       *out
00315         <<  "Constituent LinearOpBase objects for M = [ Op[0,0] ... ; ... ; ... 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::RefCountPtr<const LinearOpBase<Scalar> >
00321             block_i_j = getBlock(i,j);
00322           if(block_i_j.get())
00323             *out << "\n" << 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::RefCountPtr<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::RefCountPtr;
00366   using Teuchos::dyn_cast;
00367   typedef Teuchos::ScalarTraits<Scalar>                ST;
00368   typedef RefCountPtr<MultiVectorBase<Scalar> >        MultiVectorPtr;
00369   typedef RefCountPtr<const MultiVectorBase<Scalar> >  ConstMultiVectorPtr;
00370   typedef RefCountPtr<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     Op_i_j = ( !struct_transp ? getBlock(i,j) : getBlock(j,i) );
00398       ConstMultiVectorPtr  X_j    = X.getMultiVectorBlock(j);
00399       if(j==0) {
00400         if(Op_i_j.get())  Thyra::apply(*Op_i_j,M_trans,*X_j,&*Y_i,alpha,beta);
00401         else              scale(beta,&*Y_i);
00402       }
00403       else {
00404         if(Op_i_j.get())  Thyra::apply(*Op_i_j,M_trans,*X_j,&*Y_i,alpha,ST::one());
00405       }
00406     }
00407   }
00408 }
00409 
00410 // private
00411 
00412 template<class Scalar>
00413 void DefaultBlockedLinearOp<Scalar>::resetStorage(
00414   const int numRowBlocks, const int numColBlocks
00415   )
00416 {
00417   uninitialize();
00418   numRowBlocks_ = numRowBlocks;
00419   numColBlocks_ = numColBlocks;
00420   Ops_.resize(numRowBlocks_*numColBlocks_);
00421   if(!productRange_.get()) {
00422     rangeBlocks_.resize(numRowBlocks);
00423     domainBlocks_.resize(numColBlocks);
00424   }
00425   blockFillIsActive_ = true;
00426 }
00427 
00428 template<class Scalar>
00429 void DefaultBlockedLinearOp<Scalar>::assertBlockFillIsActive(
00430   bool blockFillIsActive
00431   ) const
00432 {
00433 #ifdef TEUCHOS_DEBUG
00434   TEST_FOR_EXCEPT(!(blockFillIsActive_==blockFillIsActive));
00435 #endif
00436 }
00437 
00438 template<class Scalar>
00439 void DefaultBlockedLinearOp<Scalar>::assertBlockRowCol(
00440   const int i, const int j
00441   ) const
00442 {
00443 #ifdef TEUCHOS_DEBUG
00444   TEST_FOR_EXCEPTION(
00445     !( 0 <= i ), std::logic_error
00446     ,"Error, i="<<i<<" is invalid!"
00447     );
00448   TEST_FOR_EXCEPTION(
00449     !( 0 <= j ), std::logic_error
00450     ,"Error, j="<<j<<" is invalid!"
00451     );
00452   // Only validate upper range if the number of row and column blocks is
00453   // fixed!
00454   if(Ops_.size()) {
00455     TEST_FOR_EXCEPTION(
00456       !( 0 <= i && i < numRowBlocks_ ), std::logic_error
00457       ,"Error, i="<<i<<" does not fall in the range [0,"<<numRowBlocks_-1<<"]!"
00458       );
00459     TEST_FOR_EXCEPTION(
00460       !( 0 <= j && j < numColBlocks_ ), std::logic_error
00461       ,"Error, j="<<j<<" does not fall in the range [0,"<<numColBlocks_-1<<"]!"
00462       );
00463   }
00464 #endif
00465 }
00466 
00467 template<class Scalar>
00468 void DefaultBlockedLinearOp<Scalar>::setBlockSpaces(
00469   const int i, const int j, const LinearOpBase<Scalar> &block
00470   )
00471 {
00472   typedef std::string s;
00473   using Teuchos::toString;
00474   assertBlockFillIsActive(true);
00475   assertBlockRowCol(i,j);
00476   // Validate that if the vector space block is already set that it is
00477   // compatible with the block that is being set.
00478   if( i < numRowBlocks_ && j < numColBlocks_ ) {
00479 #ifdef TEUCHOS_DEBUG
00480     Teuchos::RefCountPtr<const VectorSpaceBase<Scalar> >
00481       rangeBlock = (
00482         productRange_.get()
00483         ? productRange_->getBlock(i)
00484         : rangeBlocks_[i]
00485         ),
00486       domainBlock = (
00487         productDomain_.get()
00488         ? productDomain_->getBlock(j)
00489         : domainBlocks_[j]
00490         );
00491     if(rangeBlock.get()) {
00492       THYRA_ASSERT_VEC_SPACES_NAMES(
00493         "DefaultBlockedLinearOp<Scalar>::setBlockSpaces(i,j,block)"
00494         ,*rangeBlock,("(*productRange->getBlock("+toString(i)+"))")
00495         ,*block.range(),("(*block["+toString(i)+","+toString(j)+"].range())")
00496         );
00497     }
00498     if(domainBlock.get()) {
00499       THYRA_ASSERT_VEC_SPACES_NAMES(
00500         "DefaultBlockedLinearOp<Scalar>::setBlockSpaces(i,j,block)"
00501         ,*domainBlock,("(*productDomain->getBlock("+toString(j)+"))")
00502         ,*block.domain(),("(*block["+toString(i)+","+toString(j)+"].domain())")
00503         );
00504     }
00505 #endif // TEUCHOS_DEBUG
00506   }
00507   // Add spaces missing range and domain space blocks if we are doing a
00508   // flexible fill (otherwise these loops will not be executed)
00509   for( int k = numRowBlocks_; k <= i; ++k )
00510     rangeBlocks_.push_back(Teuchos::null);
00511   for( int k = numColBlocks_; k <= j; ++k )
00512     domainBlocks_.push_back(Teuchos::null);
00513   // Set the incoming range and domain blocks if not already set
00514   if(!productRange_.get()) {
00515     if(!rangeBlocks_[i].get())
00516       rangeBlocks_[i] = block.range();
00517     if(!domainBlocks_[j].get()) {
00518       domainBlocks_[j] = block.domain();
00519     }
00520   }
00521   // Update the current number of row and columns blocks if doing a flexible
00522   // fill.
00523   if(!Ops_.size()) {
00524     numRowBlocks_ = rangeBlocks_.size();
00525     numColBlocks_ = domainBlocks_.size();
00526   }
00527 }
00528 
00529 template<class Scalar>
00530 template<class LinearOpType>
00531 void DefaultBlockedLinearOp<Scalar>::setBlockImpl(
00532   const int i, const int j
00533   ,const Teuchos::RefCountPtr<LinearOpType> &block
00534   )
00535 {
00536   setBlockSpaces(i,j,*block);
00537   if(Ops_.size()) {
00538     // We are doing a fill with a fixed number of row and column blocks so we
00539     // can just set this.
00540     Ops_[numRowBlocks_*j+i] = block;
00541   }
00542   else {
00543     // We are doing a flexible fill so add the block to the stack of blocks or
00544     // replace a block that already exists.
00545     bool foundBlock = false;
00546     for( unsigned int k = 0; k < Ops_stack_.size(); ++k ) {
00547       BlockEntry<Scalar> &block_i_j = Ops_stack_[k];
00548       if( block_i_j.i == i && block_i_j.j == j ) {
00549         block_i_j.block = block;
00550         foundBlock = true;
00551         break;
00552       }
00553     }
00554     if(!foundBlock)
00555       Ops_stack_.push_back(BlockEntry<Scalar>(i,j,block));
00556   }
00557 }
00558 
00559 } // namespace Thyra
00560 
00561 template<class Scalar>
00562 Teuchos::RefCountPtr<const Thyra::LinearOpBase<Scalar> >
00563 Thyra::block2x2(
00564   const Teuchos::RefCountPtr<const LinearOpBase<Scalar> >    &A00
00565   ,const Teuchos::RefCountPtr<const LinearOpBase<Scalar> >   &A01
00566   ,const Teuchos::RefCountPtr<const LinearOpBase<Scalar> >   &A10
00567   ,const Teuchos::RefCountPtr<const LinearOpBase<Scalar> >   &A11
00568   )
00569 {
00570   Teuchos::RefCountPtr<PhysicallyBlockedLinearOpBase<Scalar> >
00571     M = Teuchos::rcp(new DefaultBlockedLinearOp<Scalar>());
00572   M->beginBlockFill(2,2);
00573   if(A00.get()) M->setBlock(0,0,A00);
00574   if(A01.get()) M->setBlock(0,1,A01);
00575   if(A10.get()) M->setBlock(1,0,A10);
00576   if(A11.get()) M->setBlock(1,1,A11);
00577   M->endBlockFill();
00578   return M;
00579 }
00580 
00581 template<class Scalar>
00582 Teuchos::RefCountPtr<const Thyra::LinearOpBase<Scalar> >
00583 Thyra::block2x1(
00584   const Teuchos::RefCountPtr<const LinearOpBase<Scalar> >    &A00
00585   ,const Teuchos::RefCountPtr<const LinearOpBase<Scalar> >   &A10
00586   )
00587 {
00588   Teuchos::RefCountPtr<PhysicallyBlockedLinearOpBase<Scalar> >
00589     M = Teuchos::rcp(new DefaultBlockedLinearOp<Scalar>());
00590   M->beginBlockFill(2,1);
00591   if(A00.get()) M->setBlock(0,0,A00);
00592   if(A10.get()) M->setBlock(1,0,A10);
00593   M->endBlockFill();
00594   return M;
00595 }
00596 
00597 template<class Scalar>
00598 Teuchos::RefCountPtr<const Thyra::LinearOpBase<Scalar> >
00599 Thyra::block1x2(
00600   const Teuchos::RefCountPtr<const LinearOpBase<Scalar> >    &A00
00601   ,const Teuchos::RefCountPtr<const LinearOpBase<Scalar> >   &A01
00602   )
00603 {
00604   Teuchos::RefCountPtr<PhysicallyBlockedLinearOpBase<Scalar> >
00605     M = Teuchos::rcp(new DefaultBlockedLinearOp<Scalar>());
00606   M->beginBlockFill(1,2);
00607   if(A00.get()) M->setBlock(0,0,A00);
00608   if(A01.get()) M->setBlock(0,1,A01);
00609   M->endBlockFill();
00610   return M;
00611 }
00612 
00613 template<class Scalar>
00614 Teuchos::RefCountPtr<Thyra::LinearOpBase<Scalar> >
00615 Thyra::nonconstBlock2x2(
00616   const Teuchos::RefCountPtr<LinearOpBase<Scalar> >    &A00
00617   ,const Teuchos::RefCountPtr<LinearOpBase<Scalar> >   &A01
00618   ,const Teuchos::RefCountPtr<LinearOpBase<Scalar> >   &A10
00619   ,const Teuchos::RefCountPtr<LinearOpBase<Scalar> >   &A11
00620   )
00621 {
00622   Teuchos::RefCountPtr<PhysicallyBlockedLinearOpBase<Scalar> >
00623     M = Teuchos::rcp(new DefaultBlockedLinearOp<Scalar>());
00624   M->beginBlockFill(2,2);
00625   if(A00.get()) M->setNonconstBlock(0,0,A00);
00626   if(A01.get()) M->setNonconstBlock(0,1,A01);
00627   if(A10.get()) M->setNonconstBlock(1,0,A10);
00628   if(A11.get()) M->setNonconstBlock(1,1,A11);
00629   M->endBlockFill();
00630   return M;
00631 }
00632 
00633 template<class Scalar>
00634 Teuchos::RefCountPtr<Thyra::LinearOpBase<Scalar> >
00635 Thyra::nonconstBlock2x1(
00636   const Teuchos::RefCountPtr<LinearOpBase<Scalar> >    &A00
00637   ,const Teuchos::RefCountPtr<LinearOpBase<Scalar> >   &A10
00638   )
00639 {
00640   Teuchos::RefCountPtr<PhysicallyBlockedLinearOpBase<Scalar> >
00641     M = Teuchos::rcp(new DefaultBlockedLinearOp<Scalar>());
00642   M->beginBlockFill(2,1);
00643   if(A00.get()) M->setNonconstBlock(0,0,A00);
00644   if(A10.get()) M->setNonconstBlock(1,0,A10);
00645   M->endBlockFill();
00646   return M;
00647 }
00648 
00649 template<class Scalar>
00650 Teuchos::RefCountPtr<Thyra::LinearOpBase<Scalar> >
00651 Thyra::nonconstBlock1x2(
00652   const Teuchos::RefCountPtr<LinearOpBase<Scalar> >    &A00
00653   ,const Teuchos::RefCountPtr<LinearOpBase<Scalar> >   &A01
00654   )
00655 {
00656   Teuchos::RefCountPtr<PhysicallyBlockedLinearOpBase<Scalar> >
00657     M = Teuchos::rcp(new DefaultBlockedLinearOp<Scalar>());
00658   M->beginBlockFill(1,2);
00659   if(A00.get()) M->setNonconstBlock(0,0,A00);
00660   if(A01.get()) M->setNonconstBlock(0,1,A01);
00661   M->endBlockFill();
00662   return M;
00663 }
00664 
00665 #endif  // THYRA_DEFAULT_BLOCKED_LINEAR_OP_HPP

Generated on Thu Sep 18 12:33:02 2008 for Thyra Package Browser (Single Doxygen Collection) by doxygen 1.3.9.1