Teko  Version of the Day
 All Classes Files Functions Variables Pages
Teko_BlockLowerTriInverseOp.cpp
1 /*
2 // @HEADER
3 //
4 // ***********************************************************************
5 //
6 // Teko: A package for block and physics based preconditioning
7 // Copyright 2010 Sandia Corporation
8 //
9 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
10 // the U.S. Government retains certain rights in this software.
11 //
12 // Redistribution and use in source and binary forms, with or without
13 // modification, are permitted provided that the following conditions are
14 // met:
15 //
16 // 1. Redistributions of source code must retain the above copyright
17 // notice, this list of conditions and the following disclaimer.
18 //
19 // 2. Redistributions in binary form must reproduce the above copyright
20 // notice, this list of conditions and the following disclaimer in the
21 // documentation and/or other materials provided with the distribution.
22 //
23 // 3. Neither the name of the Corporation nor the names of the
24 // contributors may be used to endorse or promote products derived from
25 // this software without specific prior written permission.
26 //
27 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
28 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
29 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
30 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
31 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
32 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
33 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
34 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
35 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
36 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
37 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
38 //
39 // Questions? Contact Eric C. Cyr (eccyr@sandia.gov)
40 //
41 // ***********************************************************************
42 //
43 // @HEADER
44 
45 */
46 
47 #include "Teko_BlockLowerTriInverseOp.hpp"
48 
49 #include "Teuchos_Utils.hpp"
50 
51 namespace Teko {
52 
53 using Teuchos::RCP;
54 
55 BlockLowerTriInverseOp::BlockLowerTriInverseOp(BlockedLinearOp& L,
56  const std::vector<LinearOp>& invDiag)
57  : L_(L) {
58  invDiag_ = invDiag;
59 
60  // sanity check
61  int blocks = blockRowCount(L_);
62  TEUCHOS_ASSERT(blocks > 0);
63  TEUCHOS_ASSERT(blocks == blockColCount(L_));
64  TEUCHOS_ASSERT(blocks == (int)invDiag_.size());
65 
66  // create the range and product space
68 
69  // just flip flop them!
70  productRange_ = L_->productDomain();
71  productDomain_ = L_->productRange();
72 }
73 
86 void BlockLowerTriInverseOp::implicitApply(const BlockedMultiVector& src, BlockedMultiVector& dst,
87  const double alpha, const double beta) const {
88  int blocks = blockCount(src);
89 
90  TEUCHOS_ASSERT(blocks == blockRowCount(L_));
91  TEUCHOS_ASSERT(blocks == blockCount(dst));
92 
93  // build a scrap vector for storing work
94  srcScrap_ = datacopy(src, srcScrap_);
95  BlockedMultiVector dstCopy;
96  if (beta != 0.0) {
97  dstScrap_ = datacopy(dst, dstScrap_);
98  dstCopy = dstScrap_;
99  } else
100  dstCopy = dst; // shallow copy
101 
102  // extract the blocks componets from
103  // the source and destination vectors
104  std::vector<MultiVector> dstVec;
105  std::vector<MultiVector> scrapVec;
106  for (int b = 0; b < blocks; b++) {
107  dstVec.push_back(getBlock(b, dstCopy));
108  scrapVec.push_back(getBlock(b, srcScrap_));
109  }
110 
111  // run forward-substituion: run over each column
112  // From Heath pg. 65
113  for (int b = 0; b < blocks; b++) {
114  applyOp(invDiag_[b], scrapVec[b], dstVec[b]);
115 
116  // loop over each row
117  for (int i = b + 1; i < blocks; i++) {
118  LinearOp u_ib = getBlock(i, b, L_);
119  if (u_ib != Teuchos::null) {
120  applyOp(u_ib, dstVec[b], scrapVec[i], -1.0, 1.0);
121  }
122  }
123  }
124 
125  // scale result by alpha
126  if (beta != 0)
127  update(alpha, dstCopy, beta, dst); // dst = alpha * dstCopy + beta * dst
128  else if (alpha != 1.0)
129  scale(alpha, dst); // dst = alpha * dst
130 }
131 
132 void BlockLowerTriInverseOp::describe(Teuchos::FancyOStream& out_arg,
133  const Teuchos::EVerbosityLevel verbLevel) const {
134  using Teuchos::OSTab;
135 
136  RCP<Teuchos::FancyOStream> out = rcp(&out_arg, false);
137  OSTab tab(out);
138  switch (verbLevel) {
139  case Teuchos::VERB_DEFAULT:
140  case Teuchos::VERB_LOW: *out << this->description() << std::endl; break;
141  case Teuchos::VERB_MEDIUM:
142  case Teuchos::VERB_HIGH:
143  case Teuchos::VERB_EXTREME: {
144  *out << Teuchos::Describable::description() << "{"
145  << "rangeDim=" << this->range()->dim() << ",domainDim=" << this->domain()->dim()
146  << ",rows=" << blockRowCount(L_) << ",cols=" << blockColCount(L_) << "}\n";
147  {
148  OSTab tab2(out);
149  *out << "[L Operator] = ";
150  *out << Teuchos::describe(*L_, verbLevel);
151  }
152  {
153  OSTab tab2(out);
154  *out << "[invDiag Operators]:\n";
155  tab.incrTab();
156  for (int i = 0; i < blockRowCount(L_); i++) {
157  *out << "[invD(" << i << ")] = ";
158  *out << Teuchos::describe(*invDiag_[i], verbLevel);
159  }
160  }
161  break;
162  }
163  default: TEUCHOS_TEST_FOR_EXCEPT(true); // Should never get here!
164  }
165 }
166 
167 } // end namespace Teko
virtual VectorSpace domain() const
Domain space of this operator.
const BlockedLinearOp L_
operator
Teuchos::RCP< const Thyra::ProductVectorSpaceBase< double > > productDomain_
Domain vector space.
virtual void implicitApply(const BlockedMultiVector &x, BlockedMultiVector &y, const double alpha=1.0, const double beta=0.0) const
Perform a matrix vector multiply with this operator.
std::vector< LinearOp > invDiag_
(Approximate) Inverses of the diagonal operators
Teuchos::RCP< const Thyra::ProductVectorSpaceBase< double > > productRange_
Range vector space.
virtual VectorSpace range() const
Range space of this operator.