Teko Version of the Day
Loading...
Searching...
No Matches
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
51namespace Teko {
52
53using Teuchos::RCP;
54
65BlockUpperTriInverseOp::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
84void BlockUpperTriInverseOp::implicitApply(const BlockedMultiVector & src, BlockedMultiVector & dst,
85 const double alpha, const double beta) const
86{
87 // call the no tranpose version
88 implicitApply(Thyra::NOTRANS,src,dst,alpha,beta);
89}
90
103void BlockUpperTriInverseOp::implicitApply(const Thyra::EOpTransp M_trans,
104 const BlockedMultiVector & src, BlockedMultiVector & dst,
105 const double alpha, const double beta) const
106{
107 int blocks = blockCount(src);
108
109 TEUCHOS_ASSERT(blocks==blockRowCount(U_));
110 TEUCHOS_ASSERT(blocks==blockCount(dst));
111
112 // build a scrap vector for storing work
113 srcScrap_ = datacopy(src,srcScrap_);
114 BlockedMultiVector dstCopy;
115 if(beta!=0.0) {
116 dstScrap_ = datacopy(dst,dstScrap_);
117 dstCopy = dstScrap_;
118 }
119 else
120 dstCopy = dst; // shallow copy
121
122 // extract the blocks components from
123 // the source and destination vectors
124 std::vector<MultiVector> dstVec;
125 std::vector<MultiVector> scrapVec;
126 for(int b=0;b<blocks;b++) {
127 dstVec.push_back(getBlock(b,dstCopy));
128 scrapVec.push_back(getBlock(b,srcScrap_));
129 }
130
131 // run back-substituion: run over each column
132 // From Heath pg. 66
133 if(M_trans==Thyra::NOTRANS) {
134 for(int b=blocks-1;b>=0;b--) {
135 applyOp(invDiag_[b], scrapVec[b], dstVec[b]);
136
137 // loop over each row
138 for(int i=0;i<b;i++) {
139 LinearOp u_ib = getBlock(i,b,U_);
140 if(u_ib!=Teuchos::null) {
141 applyOp(u_ib,dstVec[b],scrapVec[i],-1.0,1.0);
142 }
143 }
144 }
145 }
146 else if(M_trans==Thyra::TRANS || M_trans==Thyra::CONJTRANS) {
147 for(int b=0;b<blocks;b++) {
148 applyTransposeOp(invDiag_[b], scrapVec[b], dstVec[b]);
149
150 // loop over each row
151 for(int i=b+1;i<blocks;i++) {
152 LinearOp u_bi = getBlock(b,i,U_);
153 if(u_bi!=Teuchos::null) {
154 applyTransposeOp(u_bi,dstVec[b],scrapVec[i],-1.0,1.0);
155 }
156 }
157 }
158 }
159 else {
160 TEUCHOS_TEST_FOR_EXCEPT(true);
161 }
162
163 // scale result by alpha
164 if(beta!=0)
165 update(alpha,dstCopy,beta,dst); // dst = alpha * dstCopy + beta * dst
166 else if(alpha!=1.0)
167 scale(alpha,dst); // dst = alpha * dst
168}
169
170void BlockUpperTriInverseOp::describe(Teuchos::FancyOStream & out_arg,
171 const Teuchos::EVerbosityLevel verbLevel) const
172{
173 using Teuchos::OSTab;
174
175 RCP<Teuchos::FancyOStream> out = rcp(&out_arg,false);
176 OSTab tab(out);
177 switch(verbLevel) {
178 case Teuchos::VERB_DEFAULT:
179 case Teuchos::VERB_LOW:
180 *out << this->description() << std::endl;
181 break;
182 case Teuchos::VERB_MEDIUM:
183 case Teuchos::VERB_HIGH:
184 case Teuchos::VERB_EXTREME:
185 {
186 *out << Teuchos::Describable::description() << "{"
187 << "rangeDim=" << this->range()->dim()
188 << ",domainDim=" << this->domain()->dim()
189 << ",rows=" << blockRowCount(U_)
190 << ",cols=" << blockColCount(U_)
191 << "}\n";
192 {
193 OSTab tab2(out);
194 *out << "[U Operator] = ";
195 *out << Teuchos::describe(*U_,verbLevel);
196 }
197 {
198 OSTab tab2(out);
199 *out << "[invDiag Operators]:\n";
200 tab.incrTab();
201 for(int i=0;i<blockRowCount(U_);i++) {
202 *out << "[invD(" << i << ")] = ";
203 *out << Teuchos::describe(*invDiag_[i],verbLevel);
204 }
205 }
206 break;
207 }
208 default:
209 TEUCHOS_TEST_FOR_EXCEPT(true); // Should never get here!
210 }
211}
212
213} // end namespace Teko
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.
Teuchos::RCP< const Thyra::ProductVectorSpaceBase< double > > productDomain_
Domain vector space.
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
virtual VectorSpace range() const
Range space of this operator.