Teko Version of the Day
Loading...
Searching...
No Matches
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
51namespace Teko {
52
53using Teuchos::RCP;
54
55BlockLowerTriInverseOp::BlockLowerTriInverseOp(BlockedLinearOp & L,const std::vector<LinearOp> & invDiag)
56 : L_(L)
57{
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
86void BlockLowerTriInverseOp::implicitApply(const BlockedMultiVector & src, BlockedMultiVector & dst,
87 const double alpha, const double beta) const
88{
89 int blocks = blockCount(src);
90
91 TEUCHOS_ASSERT(blocks==blockRowCount(L_));
92 TEUCHOS_ASSERT(blocks==blockCount(dst));
93
94 // build a scrap vector for storing work
95 srcScrap_ = datacopy(src,srcScrap_);
96 BlockedMultiVector dstCopy;
97 if(beta!=0.0) {
98 dstScrap_ = datacopy(dst,dstScrap_);
99 dstCopy = dstScrap_;
100 }
101 else
102 dstCopy = dst; // shallow copy
103
104 // extract the blocks componets from
105 // the source and destination vectors
106 std::vector<MultiVector> dstVec;
107 std::vector<MultiVector> scrapVec;
108 for(int b=0;b<blocks;b++) {
109 dstVec.push_back(getBlock(b,dstCopy));
110 scrapVec.push_back(getBlock(b,srcScrap_));
111 }
112
113 // run forward-substituion: run over each column
114 // From Heath pg. 65
115 for(int b=0;b<blocks;b++) {
116 applyOp(invDiag_[b], scrapVec[b], dstVec[b]);
117
118 // loop over each row
119 for(int i=b+1;i<blocks;i++) {
120 LinearOp u_ib = getBlock(i,b,L_);
121 if(u_ib!=Teuchos::null) {
122 applyOp(u_ib,dstVec[b],scrapVec[i],-1.0,1.0);
123 }
124 }
125 }
126
127 // scale result by alpha
128 if(beta!=0)
129 update(alpha,dstCopy,beta,dst); // dst = alpha * dstCopy + beta * dst
130 else if(alpha!=1.0)
131 scale(alpha,dst); // dst = alpha * dst
132}
133
134void BlockLowerTriInverseOp::describe(Teuchos::FancyOStream & out_arg,
135 const Teuchos::EVerbosityLevel verbLevel) const
136{
137 using Teuchos::OSTab;
138
139 RCP<Teuchos::FancyOStream> out = rcp(&out_arg,false);
140 OSTab tab(out);
141 switch(verbLevel) {
142 case Teuchos::VERB_DEFAULT:
143 case Teuchos::VERB_LOW:
144 *out << this->description() << std::endl;
145 break;
146 case Teuchos::VERB_MEDIUM:
147 case Teuchos::VERB_HIGH:
148 case Teuchos::VERB_EXTREME:
149 {
150 *out << Teuchos::Describable::description() << "{"
151 << "rangeDim=" << this->range()->dim()
152 << ",domainDim=" << this->domain()->dim()
153 << ",rows=" << blockRowCount(L_)
154 << ",cols=" << blockColCount(L_)
155 << "}\n";
156 {
157 OSTab tab2(out);
158 *out << "[L Operator] = ";
159 *out << Teuchos::describe(*L_,verbLevel);
160 }
161 {
162 OSTab tab2(out);
163 *out << "[invDiag Operators]:\n";
164 tab.incrTab();
165 for(int i=0;i<blockRowCount(L_);i++) {
166 *out << "[invD(" << i << ")] = ";
167 *out << Teuchos::describe(*invDiag_[i],verbLevel);
168 }
169 }
170 break;
171 }
172 default:
173 TEUCHOS_TEST_FOR_EXCEPT(true); // Should never get here!
174 }
175}
176
177} // end namespace Teko
void scale(const double alpha, MultiVector &x)
Scale a multivector by a constant.
MultiVector datacopy(const MultiVector &src, MultiVector &dst)
Copy the contents of a multivector to a destination vector.
int blockRowCount(const BlockedLinearOp &blo)
Get the row count in a block linear operator.
int blockCount(const BlockedMultiVector &bmv)
Get the column count in a block linear operator.
MultiVector getBlock(int i, const BlockedMultiVector &bmv)
Get the ith block from a BlockedMultiVector object.
int blockColCount(const BlockedLinearOp &blo)
Get the column count in a block linear operator.
virtual VectorSpace range() const
Range space of this operator.
virtual VectorSpace domain() const
Domain space of this operator.
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.
Teuchos::RCP< const Thyra::ProductVectorSpaceBase< double > > productRange_
Range vector space.
std::vector< LinearOp > invDiag_
(Approximate) Inverses of the diagonal operators
Teuchos::RCP< const Thyra::ProductVectorSpaceBase< double > > productDomain_
Domain vector space.