Thyra Version of the Day
Thyra_DefaultBlockedTriangularLinearOpWithSolve_def.hpp
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 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Roscoe A. Bartlett (bartlettra@ornl.gov) 
00038 // 
00039 // ***********************************************************************
00040 // @HEADER
00041 
00042 #ifndef THYRA_DEFAULT_BLOCKED_TRIANGULAR_LINEAR_OP_WITH_SOLVE_DEF_HPP
00043 #define THYRA_DEFAULT_BLOCKED_TRIANGULAR_LINEAR_OP_WITH_SOLVE_DEF_HPP
00044 
00045 
00046 #include "Thyra_DefaultBlockedTriangularLinearOpWithSolve_decl.hpp"
00047 #include "Thyra_ProductMultiVectorBase.hpp"
00048 #include "Thyra_DefaultProductVectorSpace.hpp"
00049 #include "Thyra_AssertOp.hpp"
00050 
00051 
00052 namespace Thyra {
00053 
00054 
00055 // public
00056 
00057 
00058 // Constructors/Initializers/Accessors
00059 
00060 
00061 template<class Scalar>
00062 DefaultBlockedTriangularLinearOpWithSolve<Scalar>::DefaultBlockedTriangularLinearOpWithSolve()
00063   : blockFillIsActive_(false), numDiagBlocks_(0)
00064 {}
00065 
00066 
00067 template<class Scalar>
00068 void DefaultBlockedTriangularLinearOpWithSolve<Scalar>::setNonconstBlocks(
00069   const RCP<PhysicallyBlockedLinearOpBase<Scalar> > &blocks
00070   )
00071 {
00072   assertAndSetBlockStructure(*blocks);
00073   blocks_.initialize(blocks);
00074 }
00075 
00076 
00077 template<class Scalar>
00078 void DefaultBlockedTriangularLinearOpWithSolve<Scalar>::setBlocks(
00079   const RCP<const PhysicallyBlockedLinearOpBase<Scalar> > &blocks
00080   )
00081 {
00082   assertAndSetBlockStructure(*blocks);
00083   blocks_.initialize(blocks);
00084 }
00085 
00086 
00087 template<class Scalar>
00088 RCP<PhysicallyBlockedLinearOpBase<Scalar> >
00089 DefaultBlockedTriangularLinearOpWithSolve<Scalar>::getNonconstBlocks()
00090 {
00091   return blocks_.getNonconstObj();
00092 }
00093 
00094 
00095 template<class Scalar>
00096 RCP<const PhysicallyBlockedLinearOpBase<Scalar> >
00097 DefaultBlockedTriangularLinearOpWithSolve<Scalar>::getBlocks()
00098 {
00099   return blocks_.getConstObj();
00100 }
00101 
00102 
00103 // Overridden from PhysicallyBlockedLinearOpWithSolveBase
00104 
00105 
00106 template<class Scalar>
00107 bool DefaultBlockedTriangularLinearOpWithSolve<Scalar>::acceptsLOWSBlock(
00108   const int i, const int j
00109   ) const
00110 {
00111   assertBlockFillIsActive(true);
00112   assertBlockRowCol(i,j);
00113   return i==j; // Only accept LOWS blocks along the diagonal!
00114 }
00115 
00116 template<class Scalar>
00117 void DefaultBlockedTriangularLinearOpWithSolve<Scalar>::setNonconstLOWSBlock(
00118   const int i, const int j,
00119   const Teuchos::RCP<LinearOpWithSolveBase<Scalar> > &block
00120   )
00121 {
00122   setLOWSBlockImpl(i,j,block);
00123 }
00124 
00125 
00126 template<class Scalar>
00127 void DefaultBlockedTriangularLinearOpWithSolve<Scalar>::setLOWSBlock(
00128   const int i, const int j,
00129   const Teuchos::RCP<const LinearOpWithSolveBase<Scalar> > &block
00130   )
00131 {
00132   setLOWSBlockImpl(i,j,block);
00133 }
00134 
00135 
00136 // Overridden from PhysicallyBlockedLinearOpBase
00137 
00138 
00139 template<class Scalar>
00140 void DefaultBlockedTriangularLinearOpWithSolve<Scalar>::beginBlockFill()
00141 {
00142  assertBlockFillIsActive(false);
00143  TEUCHOS_TEST_FOR_EXCEPT("ToDo: Have not implemented flexible block fill yet!");
00144 }
00145 
00146 
00147 template<class Scalar>
00148 void DefaultBlockedTriangularLinearOpWithSolve<Scalar>::beginBlockFill(
00149   const int numRowBlocks, const int numColBlocks
00150   )
00151 {
00152   using Teuchos::null;
00153 #ifdef THYRA_DEBUG
00154   TEUCHOS_ASSERT_EQUALITY(numRowBlocks, numColBlocks);
00155 #endif
00156   assertBlockFillIsActive(false);
00157   numDiagBlocks_ = numRowBlocks;
00158   diagonalBlocks_.resize(numDiagBlocks_);
00159   productRange_ = null;
00160   productDomain_ = null;
00161   blockFillIsActive_ = true;
00162 }
00163 
00164 
00165 template<class Scalar>
00166 void DefaultBlockedTriangularLinearOpWithSolve<Scalar>::beginBlockFill(
00167   const Teuchos::RCP<const ProductVectorSpaceBase<Scalar> > &productRange_in,
00168   const Teuchos::RCP<const ProductVectorSpaceBase<Scalar> > &productDomain_in
00169   )
00170 {
00171 #ifdef THYRA_DEBUG
00172   TEUCHOS_TEST_FOR_EXCEPT( is_null(productRange_in) );
00173   TEUCHOS_TEST_FOR_EXCEPT( is_null(productDomain_in) );
00174   TEUCHOS_TEST_FOR_EXCEPT( productRange_in->numBlocks() != productDomain_in->numBlocks() );
00175 #endif
00176   assertBlockFillIsActive(false);
00177   productRange_ = productRange_in;
00178   productDomain_ = productDomain_in;
00179   numDiagBlocks_ = productRange_in->numBlocks();
00180   diagonalBlocks_.resize(numDiagBlocks_);
00181   blockFillIsActive_ = true;
00182 }
00183 
00184 
00185 template<class Scalar>
00186 bool DefaultBlockedTriangularLinearOpWithSolve<Scalar>::blockFillIsActive() const
00187 {
00188   return blockFillIsActive_;
00189 }
00190 
00191 
00192 template<class Scalar>
00193 bool DefaultBlockedTriangularLinearOpWithSolve<Scalar>::acceptsBlock(
00194   const int i, const int j
00195   ) const
00196 {
00197   assertBlockFillIsActive(true);
00198   assertBlockRowCol(i,j);
00199   return false; // ToDo: Change this once we accept off-diagonal blocks
00200 }
00201 
00202 
00203 template<class Scalar>
00204 void DefaultBlockedTriangularLinearOpWithSolve<Scalar>::setNonconstBlock(
00205   const int i, const int j,
00206   const Teuchos::RCP<LinearOpBase<Scalar> > &block
00207   )
00208 {
00209   assertBlockFillIsActive(true);
00210   TEUCHOS_TEST_FOR_EXCEPT("Error, we don't support off-diagonal LOB objects yet!");
00211 }
00212 
00213 
00214 template<class Scalar>
00215 void DefaultBlockedTriangularLinearOpWithSolve<Scalar>::setBlock(
00216   const int i, const int j,
00217   const Teuchos::RCP<const LinearOpBase<Scalar> > &block
00218   )
00219 {
00220   assertBlockFillIsActive(true);
00221   TEUCHOS_TEST_FOR_EXCEPT("Error, we don't support off-diagonal LOB objects yet!");
00222 }
00223 
00224 
00225 template<class Scalar>
00226 void DefaultBlockedTriangularLinearOpWithSolve<Scalar>::endBlockFill()
00227 {
00228   assertBlockFillIsActive(true);
00229   Array<RCP<const VectorSpaceBase<Scalar> > > rangeSpaces;
00230   Array<RCP<const VectorSpaceBase<Scalar> > > domainSpaces;
00231   for ( int k = 0; k < numDiagBlocks_; ++k ) {
00232     const RCP<const LinearOpWithSolveBase<Scalar> > lows_k = 
00233       diagonalBlocks_[k].getConstObj();
00234     TEUCHOS_TEST_FOR_EXCEPTION(is_null(lows_k), std::logic_error,
00235       "Error, the block diagonal k="<<k<<" can not be null when ending block fill!"
00236       );
00237     if (is_null(productRange_)) {
00238       rangeSpaces.push_back(lows_k->range());
00239       domainSpaces.push_back(lows_k->domain());
00240     }
00241   }
00242   if (is_null(productRange_)) {
00243     productRange_ = productVectorSpace<Scalar>(rangeSpaces());
00244     productDomain_ = productVectorSpace<Scalar>(domainSpaces());
00245   }
00246   blockFillIsActive_ = false;
00247 }
00248 
00249 
00250 template<class Scalar>
00251 void DefaultBlockedTriangularLinearOpWithSolve<Scalar>::uninitialize()
00252 {
00253   assertBlockFillIsActive(false);
00254   productRange_ = Teuchos::null;
00255   productDomain_ = Teuchos::null;
00256   numDiagBlocks_ = 0;
00257   diagonalBlocks_.resize(0);
00258 }
00259 
00260 
00261 // Overridden from BlockedLinearOpWithSolveBase
00262 
00263 
00264 template<class Scalar>
00265 Teuchos::RCP<LinearOpWithSolveBase<Scalar> >
00266 DefaultBlockedTriangularLinearOpWithSolve<Scalar>::getNonconstLOWSBlock(
00267   const int i, const int j
00268   )
00269 {
00270   assertBlockFillIsActive(false);
00271   assertBlockRowCol(i,j);
00272   if (i!=j)
00273     return Teuchos::null;
00274   return diagonalBlocks_[i].getNonconstObj();
00275 } 
00276 
00277 
00278 template<class Scalar>
00279 Teuchos::RCP<const LinearOpWithSolveBase<Scalar> >
00280 DefaultBlockedTriangularLinearOpWithSolve<Scalar>::getLOWSBlock(
00281   const int i, const int j
00282   ) const
00283 {
00284   assertBlockFillIsActive(false);
00285   assertBlockRowCol(i, j);
00286   if (i != j)
00287     return Teuchos::null;
00288   return diagonalBlocks_[i].getConstObj();
00289 }
00290 
00291 
00292 // Overridden from BlockedLinearOpBase
00293 
00294 
00295 template<class Scalar>
00296 Teuchos::RCP<const ProductVectorSpaceBase<Scalar> >
00297 DefaultBlockedTriangularLinearOpWithSolve<Scalar>::productRange() const
00298 {
00299   return productRange_;
00300 }
00301 
00302 
00303 template<class Scalar>
00304 Teuchos::RCP<const ProductVectorSpaceBase<Scalar> >
00305 DefaultBlockedTriangularLinearOpWithSolve<Scalar>::productDomain() const
00306 {
00307   return productDomain_;
00308 }
00309 
00310 
00311 template<class Scalar>
00312 bool DefaultBlockedTriangularLinearOpWithSolve<Scalar>::blockExists(
00313   const int i, const int j
00314   ) const
00315 {
00316   assertBlockFillIsActive(false);
00317   assertBlockRowCol(i,j);
00318   if (i!=j)
00319     return false; // ToDo: Update this when off-diagonals are supported!
00320   return !is_null(diagonalBlocks_[i].getConstObj());
00321 } 
00322 
00323 
00324 template<class Scalar>
00325 bool DefaultBlockedTriangularLinearOpWithSolve<Scalar>::blockIsConst(
00326   const int i, const int j
00327   ) const
00328 {
00329   assertBlockFillIsActive(true);
00330   assertBlockRowCol(i,j);
00331   return diagonalBlocks_[i].isConst();
00332 } 
00333 
00334 
00335 template<class Scalar>
00336 Teuchos::RCP<LinearOpBase<Scalar> >
00337 DefaultBlockedTriangularLinearOpWithSolve<Scalar>::getNonconstBlock(
00338   const int i, const int j
00339   )
00340 {
00341   assertBlockFillIsActive(true);
00342   assertBlockRowCol(i,j);
00343   if (i!=j)
00344     return Teuchos::null; // ToDo: Update when off-diagonals are supported!
00345   return this->getNonconstLOWSBlock(i,j);
00346 } 
00347 
00348 
00349 template<class Scalar>
00350 Teuchos::RCP<const LinearOpBase<Scalar> >
00351 DefaultBlockedTriangularLinearOpWithSolve<Scalar>::getBlock(
00352   const int i, const int j
00353   ) const
00354 {
00355   assertBlockFillIsActive(true);
00356   assertBlockRowCol(i,j);
00357   if (i!=j)
00358     return Teuchos::null; // ToDo: Update when off-diagonals are supported!
00359   return this->getLOWSBlock(i,j);
00360 } 
00361 
00362 
00363 // Overridden from LinearOpBase
00364 
00365 
00366 template<class Scalar>
00367 Teuchos::RCP< const VectorSpaceBase<Scalar> >
00368 DefaultBlockedTriangularLinearOpWithSolve<Scalar>::range() const
00369 {
00370   return this->productRange();
00371 }
00372 
00373 
00374 template<class Scalar>
00375 Teuchos::RCP< const VectorSpaceBase<Scalar> >
00376 DefaultBlockedTriangularLinearOpWithSolve<Scalar>::domain() const
00377 {
00378   return this->productDomain();
00379 }
00380 
00381 
00382 template<class Scalar>
00383 Teuchos::RCP<const LinearOpBase<Scalar> >
00384 DefaultBlockedTriangularLinearOpWithSolve<Scalar>::clone() const
00385 {
00386   return Teuchos::null;  // ToDo: Implement clone when needed!
00387 }
00388 
00389 
00390 // Overridden from Teuchos::Describable
00391 
00392                                                 
00393 template<class Scalar>
00394 std::string
00395 DefaultBlockedTriangularLinearOpWithSolve<Scalar>::description() const
00396 {
00397   typedef Teuchos::ScalarTraits<Scalar>  ST;
00398   assertBlockFillIsActive(false);
00399   std::ostringstream oss;
00400   oss
00401     << Teuchos::Describable::description() << "{"
00402     << "numDiagBlocks="<<numDiagBlocks_
00403     << "}";
00404   return oss.str();
00405 }
00406 
00407 
00408 template<class Scalar>
00409 void DefaultBlockedTriangularLinearOpWithSolve<Scalar>::describe(
00410   Teuchos::FancyOStream &out,
00411   const Teuchos::EVerbosityLevel verbLevel
00412   ) const
00413 {
00414   assertBlockFillIsActive(false);
00415   Teuchos::Describable::describe(out, verbLevel);
00416   // ToDo: Fill in a better version of this!
00417 }
00418 
00419 
00420 // protected
00421 
00422 
00423 // Overridden from LinearOpBase
00424 
00425 
00426 template<class Scalar>
00427 bool DefaultBlockedTriangularLinearOpWithSolve<Scalar>::opSupportedImpl(
00428   EOpTransp M_trans
00429   ) const
00430 {
00431   using Thyra::opSupported;
00432   assertBlockFillIsActive(false);
00433   for ( int k = 0; k < numDiagBlocks_; ++k ) {
00434     if ( !opSupported(*diagonalBlocks_[k].getConstObj(),M_trans) )
00435       return false;
00436   }
00437   return true;
00438   // ToDo: To be safe we really should do a collective reduction of all
00439   // clusters of processes.  However, for the typical use case, every block
00440   // will return the same info and we should be safe!
00441 }
00442 
00443 
00444 template<class Scalar>
00445 void DefaultBlockedTriangularLinearOpWithSolve<Scalar>::applyImpl(
00446   const EOpTransp M_trans,
00447   const MultiVectorBase<Scalar> &X_in,
00448   const Ptr<MultiVectorBase<Scalar> > &Y_inout,
00449   const Scalar alpha,
00450   const Scalar beta
00451   ) const
00452 {
00453 
00454   using Teuchos::RCP;
00455   using Teuchos::dyn_cast;
00456   typedef Teuchos::ScalarTraits<Scalar> ST;
00457   using Thyra::apply;
00458 
00459 #ifdef THYRA_DEBUG
00460   THYRA_ASSERT_LINEAR_OP_MULTIVEC_APPLY_SPACES(
00461     "DefaultBlockedTriangularLinearOpWithSolve<Scalar>::apply(...)",
00462     *this, M_trans, X_in, &*Y_inout
00463     );
00464 #endif // THYRA_DEBUG  
00465 
00466   //
00467   // Y = alpha * op(M) * X + beta*Y
00468   //
00469   // =>
00470   //
00471   // Y[i] = beta+Y[i] + alpha*op(Op)[i]*X[i], for i=0...numDiagBlocks-1
00472   //
00473   // ToDo: Update to handle off diagonal blocks when needed!
00474   //
00475 
00476   const ProductMultiVectorBase<Scalar>
00477     &X = dyn_cast<const ProductMultiVectorBase<Scalar> >(X_in);
00478   ProductMultiVectorBase<Scalar>
00479     &Y = dyn_cast<ProductMultiVectorBase<Scalar> >(*Y_inout);
00480   
00481   for ( int i = 0; i < numDiagBlocks_; ++ i ) {
00482     Thyra::apply( *diagonalBlocks_[i].getConstObj(), M_trans,
00483       *X.getMultiVectorBlock(i), Y.getNonconstMultiVectorBlock(i).ptr(),
00484       alpha, beta
00485       );
00486   }
00487 
00488 }
00489 
00490 
00491 // Overridden from LinearOpWithSolveBase
00492 
00493 
00494 template<class Scalar>
00495 bool
00496 DefaultBlockedTriangularLinearOpWithSolve<Scalar>::solveSupportsImpl(
00497   EOpTransp M_trans
00498   ) const
00499 {
00500   assertBlockFillIsActive(false);
00501   for ( int k = 0; k < numDiagBlocks_; ++k ) {
00502     if ( !Thyra::solveSupports( *diagonalBlocks_[k].getConstObj(), M_trans ) )
00503       return false;
00504   }
00505   return true;
00506   // ToDo: To be safe we really should do a collective reduction of all
00507   // clusters of processes.  However, for the typical use case, every block
00508   // will return the same info and we should be safe!
00509 }
00510 
00511 
00512 template<class Scalar>
00513 bool
00514 DefaultBlockedTriangularLinearOpWithSolve<Scalar>::solveSupportsSolveMeasureTypeImpl(
00515   EOpTransp M_trans, const SolveMeasureType& solveMeasureType
00516   ) const
00517 {
00518   using Thyra::solveSupportsSolveMeasureType;
00519   assertBlockFillIsActive(false);
00520   for ( int k = 0; k < numDiagBlocks_; ++k ) {
00521     if (
00522       !solveSupportsSolveMeasureType(
00523         *diagonalBlocks_[k].getConstObj(),
00524         M_trans, solveMeasureType
00525         )
00526       )
00527     {
00528       return false;
00529     }
00530   }
00531   return true;
00532 }
00533 
00534 
00535 template<class Scalar>
00536 SolveStatus<Scalar>
00537 DefaultBlockedTriangularLinearOpWithSolve<Scalar>::solveImpl(
00538   const EOpTransp M_trans,
00539   const MultiVectorBase<Scalar> &B_in,
00540   const Ptr<MultiVectorBase<Scalar> > &X_inout,
00541   const Ptr<const SolveCriteria<Scalar> > solveCriteria
00542   ) const
00543 {
00544 
00545   using Teuchos::RCP;
00546   using Teuchos::dyn_cast;
00547   typedef Teuchos::ScalarTraits<Scalar> ST;
00548   using Thyra::solve;
00549 
00550 #ifdef THYRA_DEBUG
00551   THYRA_ASSERT_LINEAR_OP_MULTIVEC_APPLY_SPACES(
00552     "DefaultBlockedTriangularLinearOpWithSolve<Scalar>::apply(...)",
00553     *this, M_trans, *X_inout, &B_in
00554     );
00555   TEUCHOS_TEST_FOR_EXCEPT(!this->solveSupportsImpl(M_trans));
00556   TEUCHOS_TEST_FOR_EXCEPTION(
00557     nonnull(solveCriteria) && !solveCriteria->solveMeasureType.useDefault(),
00558     std::logic_error,
00559     "Error, we can't handle any non-default solve criteria yet!"
00560     );
00561   // ToDo: If solve criteria is to be handled, then we will have to be very
00562   // carefull how it it interpreted in terms of the individual period solves!
00563 #endif // THYRA_DEBUG  
00564 
00565   //
00566   // Y = alpha * inv(op(M)) * X + beta*Y
00567   //
00568   // =>
00569   //
00570   // X[i] = inv(op(Op[i]))*B[i], for i=0...numDiagBlocks-1
00571   //
00572   // ToDo: Update to handle off diagonal blocks when needed!
00573   //
00574 
00575   const ProductMultiVectorBase<Scalar>
00576     &B = dyn_cast<const ProductMultiVectorBase<Scalar> >(B_in);
00577   ProductMultiVectorBase<Scalar>
00578     &X = dyn_cast<ProductMultiVectorBase<Scalar> >(*X_inout);
00579   
00580   for ( int i = 0; i < numDiagBlocks_; ++ i ) {
00581     const RCP<const LinearOpWithSolveBase<Scalar> >
00582       Op_k = diagonalBlocks_[i].getConstObj();
00583     Op_k->setOStream(this->getOStream());
00584     Op_k->setVerbLevel(this->getVerbLevel());
00585     Thyra::solve( *Op_k, M_trans, *B.getMultiVectorBlock(i),
00586       X.getNonconstMultiVectorBlock(i).ptr() );
00587     // ToDo: Pass in solve criteria when needed!
00588   }
00589 
00590   return SolveStatus<Scalar>();
00591 
00592 }
00593 
00594 
00595 
00596 // private
00597 
00598 
00599 template<class Scalar>
00600 void DefaultBlockedTriangularLinearOpWithSolve<Scalar>::assertBlockFillIsActive(
00601   bool blockFillIsActive_in
00602   ) const
00603 {
00604 #ifdef THYRA_DEBUG
00605   TEUCHOS_TEST_FOR_EXCEPT(!(blockFillIsActive_==blockFillIsActive_in));
00606 #endif
00607 }
00608 
00609 
00610 template<class Scalar>
00611 void DefaultBlockedTriangularLinearOpWithSolve<Scalar>::assertBlockRowCol(
00612   const int i, const int j
00613   ) const
00614 {
00615 #ifdef THYRA_DEBUG
00616   TEUCHOS_TEST_FOR_EXCEPTION(
00617     !( 0 <= i && i < numDiagBlocks_ ), std::logic_error,
00618     "Error, i="<<i<<" does not fall in the range [0,"<<numDiagBlocks_-1<<"]!"
00619       );
00620   TEUCHOS_TEST_FOR_EXCEPTION(
00621     !( 0 <= j && j < numDiagBlocks_ ), std::logic_error,
00622     "Error, j="<<j<<" does not fall in the range [0,"<<numDiagBlocks_-1<<"]!"
00623       );
00624 #endif
00625 }
00626 
00627 
00628 template<class Scalar>
00629 template<class LinearOpWithSolveType>
00630 void DefaultBlockedTriangularLinearOpWithSolve<Scalar>::setLOWSBlockImpl(
00631   const int i, const int j,
00632   const Teuchos::RCP<LinearOpWithSolveType> &block
00633   )
00634 {
00635   assertBlockFillIsActive(true);
00636 #ifdef THYRA_DEBUG
00637   TEUCHOS_ASSERT_INEQUALITY( i, >=, 0 );
00638   TEUCHOS_ASSERT_INEQUALITY( j, >=, 0 );
00639   TEUCHOS_ASSERT_INEQUALITY( i, <, numDiagBlocks_ );
00640   TEUCHOS_ASSERT_INEQUALITY( j, <, numDiagBlocks_ );
00641   TEUCHOS_TEST_FOR_EXCEPTION(
00642     !this->acceptsLOWSBlock(i,j), std::logic_error,
00643     "Error, this DefaultBlockedTriangularLinearOpWithSolve does not accept\n"
00644     "LOWSB objects for the block i="<<i<<", j="<<j<<"!"
00645     );
00646 #endif
00647   diagonalBlocks_[i] = block;
00648 }
00649 
00650 
00651 template<class Scalar>
00652 void DefaultBlockedTriangularLinearOpWithSolve<Scalar>::assertAndSetBlockStructure(
00653   const PhysicallyBlockedLinearOpBase<Scalar>& blocks
00654   )
00655 {
00656 #ifdef THYRA_DEBUG
00657   THYRA_ASSERT_VEC_SPACES(
00658     "DefaultBlockedTriangularLinearOpWithSolve<Scalar>::assertAndSetBlockStructure(blocks)",
00659     *blocks.range(), *this->range()
00660     );
00661   THYRA_ASSERT_VEC_SPACES(
00662     "DefaultBlockedTriangularLinearOpWithSolve<Scalar>::assertAndSetBlockStructure(blocks)",
00663     *blocks.domain(), *this->domain()
00664     );
00665   // ToDo: Make sure that all of the blocks are above or below the diagonal
00666   // but not both!
00667 #endif
00668   // ToDo: Set if this is an upper or lower triangular block operator.
00669 }
00670 
00671 
00672 } // namespace Thyra
00673 
00674 
00675 #endif  // THYRA_DEFAULT_BLOCKED_TRIANGULAR_LINEAR_OP_WITH_SOLVE_DEF_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines