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