Teko  Version of the Day
Teko_BlockUpperTriInverseOp.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_BlockUpperTriInverseOp.hpp"
48 
49 #include "Teuchos_Utils.hpp"
50 
51 namespace Teko {
52 
53 using Teuchos::RCP;
54 
65 BlockUpperTriInverseOp::BlockUpperTriInverseOp(BlockedLinearOp & U,const std::vector<LinearOp> & invDiag)
66  : U_(U)
67 {
68  invDiag_ = invDiag;
69 
70  // sanity check
71  int blocks = blockRowCount(U_);
72  TEUCHOS_ASSERT(blocks>0);
73  TEUCHOS_ASSERT(blocks==blockColCount(U_));
74  TEUCHOS_ASSERT(blocks==(int) invDiag_.size());
75 
76  // create the range and product space
78 
79  // just flip flop them!
80  productRange_ = U->productDomain();
81  productDomain_ = U->productRange();
82 }
83 
96 void BlockUpperTriInverseOp::implicitApply(const BlockedMultiVector & src, BlockedMultiVector & dst,
97  const double alpha, const double beta) const
98 {
99  int blocks = blockCount(src);
100 
101  TEUCHOS_ASSERT(blocks==blockRowCount(U_));
102  TEUCHOS_ASSERT(blocks==blockCount(dst));
103 
104  // build a scrap vector for storing work
105  srcScrap_ = datacopy(src,srcScrap_);
106  BlockedMultiVector dstCopy;
107  if(beta!=0.0) {
108  dstScrap_ = datacopy(dst,dstScrap_);
109  dstCopy = dstScrap_;
110  }
111  else
112  dstCopy = dst; // shallow copy
113 
114  // extract the blocks components from
115  // the source and destination vectors
116  std::vector<MultiVector> dstVec;
117  std::vector<MultiVector> scrapVec;
118  for(int b=0;b<blocks;b++) {
119  dstVec.push_back(getBlock(b,dstCopy));
120  scrapVec.push_back(getBlock(b,srcScrap_));
121  }
122 
123  // run back-substituion: run over each column
124  // From Heath pg. 66
125  for(int b=blocks-1;b>=0;b--) {
126  applyOp(invDiag_[b], scrapVec[b], dstVec[b]);
127 
128  // loop over each row
129  for(int i=0;i<b;i++) {
130  LinearOp u_ib = getBlock(i,b,U_);
131  if(u_ib!=Teuchos::null) {
132  applyOp(u_ib,dstVec[b],scrapVec[i],-1.0,1.0);
133  }
134  }
135  }
136 
137  // scale result by alpha
138  if(beta!=0)
139  update(alpha,dstCopy,beta,dst); // dst = alpha * dstCopy + beta * dst
140  else if(alpha!=1.0)
141  scale(alpha,dst); // dst = alpha * dst
142 }
143 
144 void BlockUpperTriInverseOp::describe(Teuchos::FancyOStream & out_arg,
145  const Teuchos::EVerbosityLevel verbLevel) const
146 {
147  using Teuchos::OSTab;
148 
149  RCP<Teuchos::FancyOStream> out = rcp(&out_arg,false);
150  OSTab tab(out);
151  switch(verbLevel) {
152  case Teuchos::VERB_DEFAULT:
153  case Teuchos::VERB_LOW:
154  *out << this->description() << std::endl;
155  break;
156  case Teuchos::VERB_MEDIUM:
157  case Teuchos::VERB_HIGH:
158  case Teuchos::VERB_EXTREME:
159  {
160  *out << Teuchos::Describable::description() << "{"
161  << "rangeDim=" << this->range()->dim()
162  << ",domainDim=" << this->domain()->dim()
163  << ",rows=" << blockRowCount(U_)
164  << ",cols=" << blockColCount(U_)
165  << "}\n";
166  {
167  OSTab tab(out);
168  *out << "[U Operator] = ";
169  *out << Teuchos::describe(*U_,verbLevel);
170  }
171  {
172  OSTab tab(out);
173  *out << "[invDiag Operators]:\n";
174  tab.incrTab();
175  for(int i=0;i<blockRowCount(U_);i++) {
176  *out << "[invD(" << i << ")] = ";
177  *out << Teuchos::describe(*invDiag_[i],verbLevel);
178  }
179  }
180  break;
181  }
182  default:
183  TEUCHOS_TEST_FOR_EXCEPT(true); // Should never get here!
184  }
185 }
186 
187 } // end namespace Teko
int blockCount(const BlockedMultiVector &bmv)
Get the column count in a block linear operator.
virtual VectorSpace domain() const
Domain space of this operator.
int blockRowCount(const BlockedLinearOp &blo)
Get the row count in a block linear operator.
void scale(const double alpha, MultiVector &x)
Scale a multivector by a constant.
std::vector< LinearOp > invDiag_
(Approximate) Inverses of the diagonal operators
virtual VectorSpace range() const
Range space of this operator.
Teuchos::RCP< const Thyra::ProductVectorSpaceBase< double > > productRange_
Range 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.
const BlockedLinearOp U_
operator
MultiVector datacopy(const MultiVector &src, MultiVector &dst)
Copy the contents of a multivector to a destination vector.
MultiVector getBlock(int i, const BlockedMultiVector &bmv)
Get the ith block from a BlockedMultiVector object.
Teuchos::RCP< const Thyra::ProductVectorSpaceBase< double > > productDomain_
Domain vector space.
int blockColCount(const BlockedLinearOp &blo)
Get the column count in a block linear operator.