Thyra Package Browser (Single Doxygen Collection) Version of the Day
Thyra_DefaultBlockedLinearOp_def.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 
00030 #ifndef THYRA_DEFAULT_BLOCKED_LINEAR_OP_DEF_HPP
00031 #define THYRA_DEFAULT_BLOCKED_LINEAR_OP_DEF_HPP
00032 
00033 
00034 #include "Thyra_DefaultBlockedLinearOp_decl.hpp"
00035 #include "Thyra_DefaultProductVectorSpace.hpp"
00036 #include "Thyra_DefaultProductVector.hpp"
00037 #include "Thyra_DefaultProductMultiVector.hpp"
00038 #include "Thyra_MultiVectorStdOps.hpp"
00039 #include "Thyra_AssertOp.hpp"
00040 
00041 
00042 namespace Thyra {
00043 
00044 
00045 // Constructors
00046 
00047 
00048 template<class Scalar>
00049 DefaultBlockedLinearOp<Scalar>::DefaultBlockedLinearOp()
00050   :numRowBlocks_(0), numColBlocks_(0), blockFillIsActive_(false)
00051 {}
00052 
00053 
00054 // Overridden from PhysicallyBlockedLinearOpBase
00055 
00056 
00057 template<class Scalar>
00058 void DefaultBlockedLinearOp<Scalar>::beginBlockFill()
00059 {
00060   assertBlockFillIsActive(false);
00061   uninitialize();
00062   resetStorage(0,0);
00063 }
00064 
00065 
00066 template<class Scalar>
00067 void DefaultBlockedLinearOp<Scalar>::beginBlockFill(
00068   const int numRowBlocks, const int numColBlocks
00069   )
00070 {
00071   assertBlockFillIsActive(false);
00072   uninitialize();
00073   resetStorage(numRowBlocks,numColBlocks);
00074 }
00075 
00076 
00077 template<class Scalar>
00078 void DefaultBlockedLinearOp<Scalar>::beginBlockFill(
00079   const RCP<const ProductVectorSpaceBase<Scalar> > &new_productRange
00080   ,const RCP<const ProductVectorSpaceBase<Scalar> > &new_productDomain
00081   )
00082 {
00083   using Teuchos::rcp_dynamic_cast;
00084   assertBlockFillIsActive(false);
00085   uninitialize();
00086   productRange_ = new_productRange.assert_not_null();
00087   productDomain_ = new_productDomain.assert_not_null();
00088   defaultProductRange_ =
00089     rcp_dynamic_cast<const DefaultProductVectorSpace<Scalar> >(productRange_);
00090   defaultProductDomain_ =
00091     rcp_dynamic_cast<const DefaultProductVectorSpace<Scalar> >(productDomain_);
00092   // Note: the above spaces must be set before calling the next function!
00093   resetStorage(productRange_->numBlocks(), productDomain_->numBlocks());
00094 }
00095 
00096 
00097 template<class Scalar>
00098 bool DefaultBlockedLinearOp<Scalar>::blockFillIsActive() const
00099 {
00100   return blockFillIsActive_;
00101 }
00102 
00103 
00104 template<class Scalar>
00105 bool DefaultBlockedLinearOp<Scalar>::acceptsBlock(
00106   const int i, const int j
00107   ) const
00108 {
00109   assertBlockFillIsActive(true);
00110   assertBlockRowCol(i,j);
00111   return true;
00112 }
00113 
00114 
00115 template<class Scalar>
00116 void DefaultBlockedLinearOp<Scalar>::setNonconstBlock(
00117   const int i, const int j
00118   ,const RCP<LinearOpBase<Scalar> > &block
00119   )
00120 {
00121   setBlockImpl(i, j, block);
00122 }
00123 
00124 
00125 template<class Scalar>
00126 void DefaultBlockedLinearOp<Scalar>::setBlock(
00127   const int i, const int j
00128   ,const RCP<const LinearOpBase<Scalar> > &block
00129   )
00130 {
00131   setBlockImpl(i, j, block);
00132 }
00133 
00134 
00135 template<class Scalar>
00136 void DefaultBlockedLinearOp<Scalar>::endBlockFill()
00137 {
00138 
00139   using Teuchos::as;
00140 
00141   assertBlockFillIsActive(true);
00142 
00143   // 2009/05/06: rabartl: ToDo: When doing a flexible block fill
00144   // (Ops_stack_.size() > 0), we need to assert that all of the block rows and
00145   // columns have been filled in.  I don't think we do that here.
00146 
00147   // Get the number of block rows and columns
00148   if (nonnull(productRange_)) {
00149     numRowBlocks_ = productRange_->numBlocks();
00150     numColBlocks_ = productDomain_->numBlocks();
00151   }
00152   else {
00153     numRowBlocks_ = rangeBlocks_.size();
00154     numColBlocks_ = domainBlocks_.size();
00155     // NOTE: Above, whether doing a flexible fill or not, all of the blocks
00156     // must be set in order to have a valid filled operator so this
00157     // calculation should be correct.
00158   }
00159 
00160   // Assert that all of the block rows and columns have at least one entry if
00161   // the spaces were not given up front.
00162 #ifdef TEUCHOS_DEBUG
00163   if (is_null(productRange_)) {
00164     for (int i = 0; i < numRowBlocks_; ++i) {
00165       TEST_FOR_EXCEPTION(
00166         !rangeBlocks_[i].get(), std::logic_error
00167         ,"DefaultBlockedLinearOp<Scalar>::endBlockFill():"
00168         " Error, no linear operator block for the i="<<i<<" block row was added"
00169         " and we can not complete the block fill!"
00170         );
00171     }
00172     for(int j = 0; j < numColBlocks_; ++j) {
00173       TEST_FOR_EXCEPTION(
00174         !domainBlocks_[j].get(), std::logic_error
00175         ,"DefaultBlockedLinearOp<Scalar>::endBlockFill():"
00176         " Error, no linear operator block for the j="
00177         <<j<<" block column was added"
00178         " and we can not complete the block fill!"
00179         );
00180     }
00181   }
00182 #endif
00183   
00184   // Insert the block LOB objects if doing a flexible fill.
00185   if (Ops_stack_.size()) {
00186     Ops_.resize(numRowBlocks_*numColBlocks_);
00187     for ( int k = 0; k < as<int>(Ops_stack_.size()); ++k ) {
00188       const BlockEntry<Scalar> &block_i_j = Ops_stack_[k];
00189       Ops_[numRowBlocks_*block_i_j.j + block_i_j.i] = block_i_j.block;
00190     }
00191     Ops_stack_.resize(0);
00192   }
00193 
00194   // Set the product range and domain spaces if not already set
00195   if (is_null(productRange_)) {
00196     adjustBlockSpaces();
00197     defaultProductRange_ = productVectorSpace<Scalar>(rangeBlocks_());
00198     defaultProductDomain_ = productVectorSpace<Scalar>(domainBlocks_());
00199     productRange_ = defaultProductRange_;
00200     productDomain_ = defaultProductDomain_;
00201   }
00202 
00203   rangeBlocks_.resize(0);
00204   domainBlocks_.resize(0);
00205 
00206   blockFillIsActive_ = false;
00207 
00208 }
00209 
00210 
00211 template<class Scalar>
00212 void DefaultBlockedLinearOp<Scalar>::uninitialize()
00213 {
00214   productRange_ = Teuchos::null;
00215   productDomain_ = Teuchos::null;
00216   numRowBlocks_ = 0;
00217   numColBlocks_ = 0;
00218   Ops_.resize(0);
00219   Ops_stack_.resize(0);
00220   rangeBlocks_.resize(0);
00221   domainBlocks_.resize(0);
00222   blockFillIsActive_ = false;
00223 }
00224 
00225 
00226 // Overridden from BlockedLinearOpBase
00227 
00228 
00229 template<class Scalar>
00230 RCP<const ProductVectorSpaceBase<Scalar> >
00231 DefaultBlockedLinearOp<Scalar>::productRange() const
00232 {
00233   return productRange_;
00234 }
00235 
00236 
00237 template<class Scalar>
00238 RCP<const ProductVectorSpaceBase<Scalar> >
00239 DefaultBlockedLinearOp<Scalar>::productDomain() const
00240 {
00241   return productDomain_;
00242 }
00243 
00244 
00245 template<class Scalar>
00246 bool DefaultBlockedLinearOp<Scalar>::blockExists(
00247   const int i, const int j
00248   ) const
00249 {
00250   assertBlockFillIsActive(false);
00251   assertBlockRowCol(i,j);
00252   return true;
00253 } 
00254 
00255 
00256 template<class Scalar>
00257 bool DefaultBlockedLinearOp<Scalar>::blockIsConst(
00258   const int i, const int j
00259   ) const
00260 {
00261 #ifdef TEUCHOS_DEBUG
00262   TEST_FOR_EXCEPT(!blockExists(i,j));
00263 #endif
00264   assertBlockFillIsActive(false);
00265   assertBlockRowCol(i,j);
00266   return Ops_[numRowBlocks_*j+i].isConst();
00267 }
00268 
00269 
00270 template<class Scalar>
00271 RCP<LinearOpBase<Scalar> >
00272 DefaultBlockedLinearOp<Scalar>::getNonconstBlock(const int i, const int j)
00273 {
00274 #ifdef TEUCHOS_DEBUG
00275   TEST_FOR_EXCEPT(!blockExists(i,j));
00276 #endif
00277   assertBlockFillIsActive(false);
00278   assertBlockRowCol(i,j);
00279   return Ops_[numRowBlocks_*j+i].getNonconstObj();
00280 } 
00281 
00282 
00283 template<class Scalar>
00284 RCP<const LinearOpBase<Scalar> >
00285 DefaultBlockedLinearOp<Scalar>::getBlock(const int i, const int j) const
00286 {
00287 #ifdef TEUCHOS_DEBUG
00288   TEST_FOR_EXCEPT(!blockExists(i,j));
00289 #endif
00290   assertBlockFillIsActive(false);
00291   assertBlockRowCol(i,j);
00292   return Ops_[numRowBlocks_*j+i];
00293 } 
00294 
00295 
00296 // Overridden from LinearOpBase
00297 
00298 
00299 template<class Scalar>
00300 RCP< const VectorSpaceBase<Scalar> >
00301 DefaultBlockedLinearOp<Scalar>::range() const
00302 {
00303   return productRange_;
00304 }
00305 
00306 
00307 template<class Scalar>
00308 RCP< const VectorSpaceBase<Scalar> >
00309 DefaultBlockedLinearOp<Scalar>::domain() const
00310 {
00311   return productDomain_;
00312 }
00313 
00314 
00315 template<class Scalar>
00316 RCP<const LinearOpBase<Scalar> >
00317 DefaultBlockedLinearOp<Scalar>::clone() const
00318 {
00319   return Teuchos::null; // ToDo: Implement this when needed!
00320 }
00321 
00322 
00323 // Overridden from Teuchos::Describable
00324 
00325 
00326 template<class Scalar>
00327 std::string DefaultBlockedLinearOp<Scalar>::description() const
00328 {
00329   assertBlockFillIsActive(false);
00330   std::ostringstream oss;
00331   oss
00332     << Teuchos::Describable::description() << "{"
00333     << "numRowBlocks="<<numRowBlocks_
00334     << ",numColBlocks="<<numColBlocks_
00335     << "}";
00336   return oss.str();
00337 }
00338 
00339 
00340 template<class Scalar>
00341 void DefaultBlockedLinearOp<Scalar>::describe(
00342   Teuchos::FancyOStream &out_arg
00343   ,const Teuchos::EVerbosityLevel verbLevel
00344   ) const
00345 {
00346   typedef Teuchos::ScalarTraits<Scalar> ST;
00347   using Teuchos::rcpFromRef;
00348   using Teuchos::FancyOStream;
00349   using Teuchos::OSTab;
00350   assertBlockFillIsActive(false);
00351   RCP<FancyOStream> out = rcpFromRef(out_arg);
00352   OSTab tab1(out);
00353   switch(verbLevel) {
00354     case Teuchos::VERB_DEFAULT:
00355     case Teuchos::VERB_LOW:
00356       *out << this->description() << std::endl;
00357       break;
00358     case Teuchos::VERB_MEDIUM:
00359     case Teuchos::VERB_HIGH:
00360     case Teuchos::VERB_EXTREME:
00361     {
00362       *out
00363         << Teuchos::Describable::description() << "{"
00364         << "rangeDim=" << this->range()->dim()
00365         << ",domainDim=" << this->domain()->dim()
00366         << ",numRowBlocks=" << numRowBlocks_
00367         << ",numColBlocks=" << numColBlocks_
00368         << "}\n";
00369       OSTab tab2(out);
00370       *out
00371         << "Constituent LinearOpBase objects for M = [ Op[0,0] ..."
00372         << " ; ... ; ... Op[numRowBlocks-1,numColBlocks-1] ]:\n";
00373       tab2.incrTab();
00374       for( int i = 0; i < numRowBlocks_; ++i ) {
00375         for( int j = 0; j < numColBlocks_; ++j ) {
00376           *out << "Op["<<i<<","<<j<<"] = ";
00377           RCP<const LinearOpBase<Scalar> >
00378             block_i_j = getBlock(i,j);
00379           if(block_i_j.get())
00380             *out << Teuchos::describe(*getBlock(i,j),verbLevel);
00381           else
00382             *out << "NULL\n";
00383         }
00384       }
00385       break;
00386     }
00387     default:
00388       TEST_FOR_EXCEPT(true); // Should never get here!
00389   }
00390 }
00391 
00392 
00393 // protected
00394 
00395 
00396 // Overridden from LinearOpBase
00397 
00398 
00399 template<class Scalar>
00400 bool DefaultBlockedLinearOp<Scalar>::opSupportedImpl(EOpTransp M_trans) const
00401 {
00402   bool supported = true;
00403   for( int i = 0; i < numRowBlocks_; ++i ) {
00404     for( int j = 0; j < numColBlocks_; ++j ) {
00405       RCP<const LinearOpBase<Scalar> >
00406         block_i_j = getBlock(i,j);
00407       if( block_i_j.get() && !Thyra::opSupported(*block_i_j,M_trans) )
00408         supported = false;
00409     }
00410   }
00411   return supported;
00412 }
00413 
00414 
00415 template<class Scalar>
00416 void DefaultBlockedLinearOp<Scalar>::applyImpl(
00417   const EOpTransp M_trans,
00418   const MultiVectorBase<Scalar> &X_in,
00419   const Ptr<MultiVectorBase<Scalar> > &Y_inout,
00420   const Scalar alpha,
00421   const Scalar beta
00422   ) const
00423 {
00424 
00425   using Teuchos::rcpFromRef;
00426   typedef Teuchos::ScalarTraits<Scalar> ST;
00427   typedef RCP<MultiVectorBase<Scalar> > MultiVectorPtr;
00428   typedef RCP<const MultiVectorBase<Scalar> > ConstMultiVectorPtr;
00429   typedef RCP<const LinearOpBase<Scalar> > ConstLinearOpPtr;
00430 
00431 #ifdef TEUCHOS_DEBUG
00432   THYRA_ASSERT_LINEAR_OP_MULTIVEC_APPLY_SPACES(
00433     "DefaultBlockedLinearOp<Scalar>::apply(...)", *this, M_trans, X_in, &*Y_inout
00434     );
00435 #endif // TEUCHOS_DEBUG 
00436   
00437   const bool
00438     struct_transp = (real_trans(M_trans)!=NOTRANS); // Structural transpose?
00439   const int
00440     opNumRowBlocks = ( !struct_transp ? numRowBlocks_ : numColBlocks_ ), 
00441     opNumColBlocks = ( !struct_transp ? numColBlocks_ : numRowBlocks_ ); 
00442   
00443   //
00444   // Y = alpha * op(M) * X + beta*Y
00445   //
00446   // =>
00447   //
00448   // Y[i] = beta+Y[i] + sum(alpha*op(Op)[i,j]*X[j],j=0...opNumColBlocks-1)
00449   //
00450   // , for i=0...opNumRowBlocks-1
00451   //
00452 
00453   const RCP<const DefaultProductVectorSpace<Scalar> >
00454     defaultProductRange_op = ( real_trans(M_trans)==NOTRANS
00455       ? defaultProductRange_ : defaultProductDomain_ ),
00456     defaultProductDomain_op = ( real_trans(M_trans)==NOTRANS
00457       ? defaultProductDomain_ : defaultProductRange_ );
00458 
00459   const RCP<const ProductMultiVectorBase<Scalar> >
00460     X = castOrCreateSingleBlockProductMultiVector<Scalar>(
00461       defaultProductDomain_op, rcpFromRef(X_in));
00462   const RCP<ProductMultiVectorBase<Scalar> >
00463     Y = nonconstCastOrCreateSingleBlockProductMultiVector<Scalar>(
00464       defaultProductRange_op, rcpFromPtr(Y_inout));
00465 
00466   for( int i = 0; i < opNumRowBlocks; ++i ) {
00467     MultiVectorPtr Y_i = Y->getNonconstMultiVectorBlock(i);
00468     for( int j = 0; j < opNumColBlocks; ++j ) {
00469       ConstLinearOpPtr
00470         Op_i_j = ( !struct_transp ? getBlock(i,j) : getBlock(j,i) );
00471       ConstMultiVectorPtr
00472         X_j = X->getMultiVectorBlock(j);
00473       if(j==0) {
00474         if (nonnull(Op_i_j))
00475           Thyra::apply(*Op_i_j, M_trans,* X_j, Y_i.ptr(), alpha, beta);
00476         else
00477           scale(beta, Y_i.ptr());
00478       }
00479       else {
00480         if (nonnull(Op_i_j))
00481           Thyra::apply(*Op_i_j, M_trans, *X_j, Y_i.ptr(), alpha, ST::one());
00482       }
00483     }
00484   }
00485 }
00486 
00487 
00488 // private
00489 
00490 
00491 template<class Scalar>
00492 void DefaultBlockedLinearOp<Scalar>::resetStorage(
00493   const int numRowBlocks, const int numColBlocks
00494   )
00495 {
00496   numRowBlocks_ = numRowBlocks;
00497   numColBlocks_ = numColBlocks;
00498   Ops_.resize(numRowBlocks_*numColBlocks_);
00499   if (is_null(productRange_)) {
00500     rangeBlocks_.resize(numRowBlocks);
00501     domainBlocks_.resize(numColBlocks);
00502   }
00503   blockFillIsActive_ = true;
00504 }
00505 
00506 
00507 template<class Scalar>
00508 void DefaultBlockedLinearOp<Scalar>::assertBlockFillIsActive(
00509   bool wantedValue
00510   ) const
00511 {
00512 #ifdef TEUCHOS_DEBUG
00513   TEST_FOR_EXCEPT(!(blockFillIsActive_==wantedValue));
00514 #endif
00515 }
00516 
00517 
00518 template<class Scalar>
00519 void DefaultBlockedLinearOp<Scalar>::assertBlockRowCol(
00520   const int i, const int j
00521   ) const
00522 {
00523 #ifdef TEUCHOS_DEBUG
00524   TEST_FOR_EXCEPTION(
00525     !( 0 <= i ), std::logic_error
00526     ,"Error, i="<<i<<" is invalid!"
00527     );
00528   TEST_FOR_EXCEPTION(
00529     !( 0 <= j ), std::logic_error
00530     ,"Error, j="<<j<<" is invalid!"
00531     );
00532   // Only validate upper range if the number of row and column blocks is
00533   // fixed!
00534   if(Ops_.size()) {
00535     TEST_FOR_EXCEPTION(
00536       !( 0 <= i && i < numRowBlocks_ ), std::logic_error
00537       ,"Error, i="<<i<<" does not fall in the range [0,"<<numRowBlocks_-1<<"]!"
00538       );
00539     TEST_FOR_EXCEPTION(
00540       !( 0 <= j && j < numColBlocks_ ), std::logic_error
00541       ,"Error, j="<<j<<" does not fall in the range [0,"<<numColBlocks_-1<<"]!"
00542       );
00543   }
00544 #endif
00545 }
00546 
00547 
00548 template<class Scalar>
00549 void DefaultBlockedLinearOp<Scalar>::setBlockSpaces(
00550   const int i, const int j, const LinearOpBase<Scalar> &block
00551   )
00552 {
00553 
00554   typedef std::string s;
00555   using Teuchos::toString;
00556   assertBlockFillIsActive(true);
00557   assertBlockRowCol(i,j);
00558 
00559   // Validate that if the vector space block is already set that it is
00560   // compatible with the block that is being set.
00561   if( i < numRowBlocks_ && j < numColBlocks_ ) {
00562 #ifdef TEUCHOS_DEBUG
00563     RCP<const VectorSpaceBase<Scalar> >
00564       rangeBlock = (
00565         productRange_.get()
00566         ? productRange_->getBlock(i)
00567         : rangeBlocks_[i]
00568         ),
00569       domainBlock = (
00570         productDomain_.get()
00571         ? productDomain_->getBlock(j)
00572         : domainBlocks_[j]
00573         );
00574     if(rangeBlock.get()) {
00575       THYRA_ASSERT_VEC_SPACES_NAMES(
00576         "DefaultBlockedLinearOp<Scalar>::setBlockSpaces(i,j,block):\n\n"
00577         "Adding block: " + block.description(),
00578         *rangeBlock,("(*productRange->getBlock("+toString(i)+"))"),
00579         *block.range(),("(*block["+toString(i)+","+toString(j)+"].range())")
00580         );
00581     }
00582     if(domainBlock.get()) {
00583       THYRA_ASSERT_VEC_SPACES_NAMES(
00584         "DefaultBlockedLinearOp<Scalar>::setBlockSpaces(i,j,block):\n\n"
00585         "Adding block: " + block.description(),
00586         *domainBlock,("(*productDomain->getBlock("+toString(j)+"))"),
00587         *block.domain(),("(*block["+toString(i)+","+toString(j)+"].domain())")
00588         );
00589     }
00590 #endif // TEUCHOS_DEBUG
00591   }
00592 
00593   // Add spaces missing range and domain space blocks if we are doing a
00594   // flexible fill (otherwise these loops will not be executed)
00595   for( int k = numRowBlocks_; k <= i; ++k )
00596     rangeBlocks_.push_back(Teuchos::null);
00597   for( int k = numColBlocks_; k <= j; ++k )
00598     domainBlocks_.push_back(Teuchos::null);
00599 
00600   // Set the incoming range and domain blocks if not already set
00601   if(!productRange_.get()) {
00602     if(!rangeBlocks_[i].get())
00603       rangeBlocks_[i] = block.range().assert_not_null();
00604     if(!domainBlocks_[j].get()) {
00605       domainBlocks_[j] = block.domain().assert_not_null();
00606     }
00607   }
00608 
00609   // Update the current number of row and columns blocks if doing a flexible
00610   // fill.
00611   if(!Ops_.size()) {
00612     numRowBlocks_ = rangeBlocks_.size();
00613     numColBlocks_ = domainBlocks_.size();
00614   }
00615 
00616 }
00617 
00618 
00619 template<class Scalar>
00620 template<class LinearOpType>
00621 void DefaultBlockedLinearOp<Scalar>::setBlockImpl(
00622   const int i, const int j,
00623   const RCP<LinearOpType> &block
00624   )
00625 {
00626   setBlockSpaces(i, j, *block);
00627   if (Ops_.size()) {
00628     // We are doing a fill with a fixed number of row and column blocks so we
00629     // can just set this.
00630     Ops_[numRowBlocks_*j+i] = block;
00631   }
00632   else {
00633     // We are doing a flexible fill so add the block to the stack of blocks or
00634     // replace a block that already exists.
00635     bool foundBlock = false;
00636     for( unsigned int k = 0; k < Ops_stack_.size(); ++k ) {
00637       BlockEntry<Scalar> &block_i_j = Ops_stack_[k];
00638       if( block_i_j.i == i && block_i_j.j == j ) {
00639         block_i_j.block = block;
00640         foundBlock = true;
00641         break;
00642       }
00643     }
00644     if(!foundBlock)
00645       Ops_stack_.push_back(BlockEntry<Scalar>(i,j,block));
00646   }
00647 }
00648 
00649 
00650 template<class Scalar>
00651 void DefaultBlockedLinearOp<Scalar>::adjustBlockSpaces()
00652 {
00653 
00654 #ifdef TEUCHOS_DEBUG
00655   TEUCHOS_ASSERT_INEQUALITY(Ops_.size(), !=, 0);
00656 #endif
00657 
00658   //
00659   // Loop through the rows and columns looking for rows with mixed
00660   // single-space range and/or domain spaces on operators and set the single
00661   // spaces as the block space if it exists.
00662   //
00663   // NOTE: Once we get here, we can safely assume that all of the operators
00664   // are compatible w.r.t. their spaces so if there are rows and/or columns
00665   // with mixed product and single vector spaces that we can just pick the
00666   // single vector space for the whole row and/or column.
00667   //
00668 
00669   // Adjust blocks in the range space
00670   for (int i = 0; i < numRowBlocks_; ++i) {
00671     for (int j = 0; j < numColBlocks_; ++j) {
00672       const RCP<const LinearOpBase<Scalar> >
00673         op_i_j = Ops_[numRowBlocks_*j+i];
00674       if (is_null(op_i_j))
00675         continue;
00676       const RCP<const VectorSpaceBase<Scalar> > range_i_j = op_i_j->range();
00677       if (is_null(productVectorSpaceBase<Scalar>(range_i_j, false))) {
00678         rangeBlocks_[i] = range_i_j;
00679         break;
00680       }
00681     }
00682   }
00683 
00684   // Adjust blocks in the domain space
00685   for (int j = 0; j < numColBlocks_; ++j) {
00686     for (int i = 0; i < numRowBlocks_; ++i) {
00687       const RCP<const LinearOpBase<Scalar> >
00688         op_i_j = Ops_[numRowBlocks_*j+i];
00689       if (is_null(op_i_j))
00690         continue;
00691       const RCP<const VectorSpaceBase<Scalar> >
00692         domain_i_j = op_i_j->domain();
00693       if (is_null(productVectorSpaceBase<Scalar>(domain_i_j, false))) {
00694         domainBlocks_[j] = domain_i_j;
00695         break;
00696       }
00697     }
00698   }
00699 
00700 }
00701 
00702 
00703 } // namespace Thyra
00704 
00705 
00706 template<class Scalar>
00707 Teuchos::RCP<Thyra::DefaultBlockedLinearOp<Scalar> >
00708 Thyra::defaultBlockedLinearOp()
00709 {
00710   return Teuchos::rcp(new DefaultBlockedLinearOp<Scalar>());
00711 }
00712 
00713 
00714 template<class Scalar>
00715 Teuchos::RCP<const Thyra::LinearOpBase<Scalar> >
00716 Thyra::block1x1(
00717   const RCP<const LinearOpBase<Scalar> > &A00,
00718   const std::string &label
00719   )
00720 {
00721   RCP<PhysicallyBlockedLinearOpBase<Scalar> >
00722     M = defaultBlockedLinearOp<Scalar>();
00723   M->beginBlockFill(1,1);
00724   M->setBlock(0, 0, A00);
00725   M->endBlockFill();
00726   if (label.length())
00727     M->setObjectLabel(label);
00728   return M;
00729 }
00730 
00731 
00732 template<class Scalar>
00733 Teuchos::RCP<const Thyra::LinearOpBase<Scalar> >
00734 Thyra::block1x2(
00735   const RCP<const LinearOpBase<Scalar> > &A00,
00736   const RCP<const LinearOpBase<Scalar> > &A01,
00737   const std::string &label
00738   )
00739 {
00740   RCP<PhysicallyBlockedLinearOpBase<Scalar> >
00741     M = defaultBlockedLinearOp<Scalar>();
00742   M->beginBlockFill(1,2);
00743   if(A00.get()) M->setBlock(0,0,A00);
00744   if(A01.get()) M->setBlock(0,1,A01);
00745   M->endBlockFill();
00746   if (label.length())
00747     M->setObjectLabel(label);
00748   return M;
00749 }
00750 
00751 
00752 template<class Scalar>
00753 Teuchos::RCP<const Thyra::LinearOpBase<Scalar> >
00754 Thyra::block2x1(
00755   const RCP<const LinearOpBase<Scalar> > &A00,
00756   const RCP<const LinearOpBase<Scalar> > &A10,
00757   const std::string &label
00758   )
00759 {
00760   RCP<PhysicallyBlockedLinearOpBase<Scalar> >
00761     M = defaultBlockedLinearOp<Scalar>();
00762   M->beginBlockFill(2,1);
00763   if(A00.get()) M->setBlock(0,0,A00);
00764   if(A10.get()) M->setBlock(1,0,A10);
00765   M->endBlockFill();
00766   if (label.length())
00767     M->setObjectLabel(label);
00768   return M;
00769 }
00770 
00771 
00772 template<class Scalar>
00773 Teuchos::RCP<const Thyra::LinearOpBase<Scalar> >
00774 Thyra::block2x2(
00775   const RCP<const LinearOpBase<Scalar> > &A00,
00776   const RCP<const LinearOpBase<Scalar> > &A01,
00777   const RCP<const LinearOpBase<Scalar> > &A10,
00778   const RCP<const LinearOpBase<Scalar> > &A11,
00779   const std::string &label
00780   )
00781 {
00782   RCP<PhysicallyBlockedLinearOpBase<Scalar> >
00783     M = defaultBlockedLinearOp<Scalar>();
00784   M->beginBlockFill(2,2);
00785   if(A00.get()) M->setBlock(0,0,A00);
00786   if(A01.get()) M->setBlock(0,1,A01);
00787   if(A10.get()) M->setBlock(1,0,A10);
00788   if(A11.get()) M->setBlock(1,1,A11);
00789   M->endBlockFill();
00790   if (label.length())
00791     M->setObjectLabel(label);
00792   return M;
00793 }
00794 
00795 
00796 template<class Scalar>
00797 Teuchos::RCP<Thyra::LinearOpBase<Scalar> >
00798 Thyra::nonconstBlock1x1(
00799   const RCP<LinearOpBase<Scalar> > &A00,
00800   const std::string &label
00801   )
00802 {
00803   RCP<PhysicallyBlockedLinearOpBase<Scalar> >
00804     M = defaultBlockedLinearOp<Scalar>();
00805   M->beginBlockFill(1, 1);
00806   M->setNonconstBlock(0, 0, A00);
00807   M->endBlockFill();
00808   if (label.length())
00809     M->setObjectLabel(label);
00810   return M;
00811 }
00812 
00813 
00814 template<class Scalar>
00815 Teuchos::RCP<Thyra::LinearOpBase<Scalar> >
00816 Thyra::nonconstBlock1x2(
00817   const RCP<LinearOpBase<Scalar> > &A00,
00818   const RCP<LinearOpBase<Scalar> > &A01,
00819   const std::string &label
00820   )
00821 {
00822   RCP<PhysicallyBlockedLinearOpBase<Scalar> >
00823     M = defaultBlockedLinearOp<Scalar>();
00824   M->beginBlockFill(1,2);
00825   if(A00.get()) M->setNonconstBlock(0,0,A00);
00826   if(A01.get()) M->setNonconstBlock(0,1,A01);
00827   M->endBlockFill();
00828   if (label.length())
00829     M->setObjectLabel(label);
00830   return M;
00831 }
00832 
00833 
00834 template<class Scalar>
00835 Teuchos::RCP<Thyra::LinearOpBase<Scalar> >
00836 Thyra::nonconstBlock2x1(
00837   const RCP<LinearOpBase<Scalar> > &A00,
00838   const RCP<LinearOpBase<Scalar> > &A10,
00839   const std::string &label
00840   )
00841 {
00842   RCP<PhysicallyBlockedLinearOpBase<Scalar> >
00843     M = defaultBlockedLinearOp<Scalar>();
00844   M->beginBlockFill(2,1);
00845   if(A00.get()) M->setNonconstBlock(0,0,A00);
00846   if(A10.get()) M->setNonconstBlock(1,0,A10);
00847   M->endBlockFill();
00848   if (label.length())
00849     M->setObjectLabel(label);
00850   return M;
00851 }
00852 
00853 
00854 template<class Scalar>
00855 Teuchos::RCP<Thyra::LinearOpBase<Scalar> >
00856 Thyra::nonconstBlock2x2(
00857   const RCP<LinearOpBase<Scalar> > &A00,
00858   const RCP<LinearOpBase<Scalar> > &A01,
00859   const RCP<LinearOpBase<Scalar> > &A10,
00860   const RCP<LinearOpBase<Scalar> > &A11,
00861   const std::string &label
00862   )
00863 {
00864   RCP<PhysicallyBlockedLinearOpBase<Scalar> >
00865     M = defaultBlockedLinearOp<Scalar>();
00866   M->beginBlockFill(2,2);
00867   if(A00.get()) M->setNonconstBlock(0,0,A00);
00868   if(A01.get()) M->setNonconstBlock(0,1,A01);
00869   if(A10.get()) M->setNonconstBlock(1,0,A10);
00870   if(A11.get()) M->setNonconstBlock(1,1,A11);
00871   M->endBlockFill();
00872   if (label.length())
00873     M->setObjectLabel(label);
00874   return M;
00875 }
00876 
00877 
00878 //
00879 // Explicit instantiation macro
00880 //
00881 // Must be expanded from within the Thyra namespace!
00882 //
00883 
00884 
00885 #define THYRA_DEFAULT_BLOCKED_LINEAR_OP_INSTANT(SCALAR) \
00886   \
00887   template class DefaultBlockedLinearOp<SCALAR >; \
00888    \
00889   template RCP<DefaultBlockedLinearOp<SCALAR > > \
00890   defaultBlockedLinearOp<SCALAR >(); \
00891    \
00892   template RCP<const LinearOpBase<SCALAR > > \
00893   block1x1( \
00894     const RCP<const LinearOpBase<SCALAR > > &A00, \
00895     const std::string &label \
00896     ); \
00897    \
00898   template RCP<const LinearOpBase<SCALAR > > \
00899   block1x2( \
00900     const RCP<const LinearOpBase<SCALAR > > &A00, \
00901     const RCP<const LinearOpBase<SCALAR > > &A01, \
00902     const std::string &label \
00903     ); \
00904    \
00905   template RCP<const LinearOpBase<SCALAR > > \
00906   block2x1( \
00907     const RCP<const LinearOpBase<SCALAR > > &A00, \
00908     const RCP<const LinearOpBase<SCALAR > > &A10, \
00909     const std::string &label \
00910     ); \
00911    \
00912   template RCP<const LinearOpBase<SCALAR > > \
00913   block2x2( \
00914     const RCP<const LinearOpBase<SCALAR > > &A00, \
00915     const RCP<const LinearOpBase<SCALAR > > &A01, \
00916     const RCP<const LinearOpBase<SCALAR > > &A10, \
00917     const RCP<const LinearOpBase<SCALAR > > &A11, \
00918     const std::string &label \
00919     ); \
00920    \
00921   template RCP<LinearOpBase<SCALAR > > \
00922   nonconstBlock1x1( \
00923     const RCP<LinearOpBase<SCALAR > > &A00, \
00924     const std::string &label \
00925     ); \
00926    \
00927   template RCP<LinearOpBase<SCALAR > > \
00928   nonconstBlock1x2( \
00929     const RCP<LinearOpBase<SCALAR > > &A00, \
00930     const RCP<LinearOpBase<SCALAR > > &A01, \
00931     const std::string &label \
00932     ); \
00933    \
00934   template RCP<LinearOpBase<SCALAR > > \
00935   nonconstBlock2x1( \
00936     const RCP<LinearOpBase<SCALAR > > &A00, \
00937     const RCP<LinearOpBase<SCALAR > > &A10, \
00938     const std::string &label \
00939     ); \
00940    \
00941   template RCP<LinearOpBase<SCALAR > > \
00942   nonconstBlock2x2( \
00943     const RCP<LinearOpBase<SCALAR > > &A00, \
00944     const RCP<LinearOpBase<SCALAR > > &A01, \
00945     const RCP<LinearOpBase<SCALAR > > &A10, \
00946     const RCP<LinearOpBase<SCALAR > > &A11, \
00947     const std::string &label \
00948     );
00949 
00950 
00951 #endif  // THYRA_DEFAULT_BLOCKED_LINEAR_OP_DEF_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines