Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Stokhos_SchurPreconditionerImp.hpp
Go to the documentation of this file.
1/*
2// @HEADER
3// ***********************************************************************
4//
5// Stokhos Package
6// Copyright (2009) Sandia Corporation
7//
8// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
9// license for use of this work by or on behalf of the U.S. Government.
10//
11// Redistribution and use in source and binary forms, with or without
12// modification, are permitted provided that the following conditions are
13// met:
14//
15// 1. Redistributions of source code must retain the above copyright
16// notice, this list of conditions and the following disclaimer.
17//
18// 2. Redistributions in binary form must reproduce the above copyright
19// notice, this list of conditions and the following disclaimer in the
20// documentation and/or other materials provided with the distribution.
21//
22// 3. Neither the name of the Corporation nor the names of the
23// contributors may be used to endorse or promote products derived from
24// this software without specific prior written permission.
25//
26// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37//
38// Questions? Contact Eric T. Phipps (etphipp@sandia.gov).
39//
40// ***********************************************************************
41// @HEADER
42*/
43
44#include "Teuchos_SerialDenseMatrix.hpp"
45#include "Teuchos_SerialDenseSolver.hpp"
46#include "Teuchos_RCP.hpp"
47
48template <typename ordinal_type, typename value_type>
51 const Teuchos::SerialDenseMatrix<ordinal_type, value_type>& K_,
52 const ordinal_type p_, const ordinal_type m_, const ordinal_type diag_) :
53 K(K_),
54 p(p_),
55 m(m_),
56 diag(diag_)
57{
58}
59
60template <typename ordinal_type, typename value_type>
65
66template <typename ordinal_type, typename value_type>
67ordinal_type
69fact (ordinal_type n) const
70{
71 if (n > 1)
72 n = n*fact (n-1);
73 else
74 n = 1;
75 return n;
76}
77
78template <typename ordinal_type, typename value_type>
79ordinal_type
81size (ordinal_type n, ordinal_type m) const
82{
83 //n is the polynomial order and m is the number of random variables
84 // return (fact(n+m)/(fact(n)*fact(m)));
85 ordinal_type min;
86 if (n == 0 )
87 return 1;
88 else {
89 if (n<=m){
90 min = n;
91 }
92 else {
93 min = m;
94 }
95
96 ordinal_type num = n+m;
97 for (ordinal_type i=1; i<=min-1; i++)
98 num = num*(n+m-i);
99 return num/fact(min);
100 }
101}
102
103template <typename ordinal_type, typename value_type>
104ordinal_type
106ApplyInverse(const Teuchos::SerialDenseMatrix<ordinal_type, value_type>& Input,
107 Teuchos::SerialDenseMatrix<ordinal_type, value_type>& U,
108 const ordinal_type n) const
109{
110 //p: polynomial order; m: number of random variables; n: level to stop and form actual Schur Complement; diag=0: Use only diagonal of block D
111 Teuchos::SerialDenseMatrix<ordinal_type, value_type> Resid(Teuchos::Copy, Input);
112 ordinal_type ret;
113 ordinal_type lmin;
114 if (n<=1)
115 lmin=1;
116 else
117 lmin=n;
118
119 Teuchos::SerialDenseSolver<ordinal_type, value_type> solver;
120 for (ordinal_type l=p; l>=lmin; l--){
121 ordinal_type c=size(l,m);
122 ordinal_type s = size(l-1,m);
123
124 //split residual
125 Teuchos::SerialDenseMatrix<ordinal_type, value_type> rminus(Teuchos::View, Resid, s, 1);
126 Teuchos::SerialDenseMatrix<ordinal_type, value_type> r(Teuchos::Copy, Resid, c-s, 1, s, 0);
127
128 //Compute pre-correction
129 Teuchos::SerialDenseMatrix<ordinal_type, value_type> B(Teuchos::View, K, s, c-s, 0, s);
130 Teuchos::SerialDenseMatrix<ordinal_type, value_type> D(Teuchos::View, K, c-s, c-s, s,s);
131
132 //Computing inv(D)r
133 if (diag == 0){
134 //For D diagonal
135 for (ordinal_type i=0; i<c-s; i++)
136 r(i,0)=r(i,0)/D(i,i);
137
138 }
139 else{
140 //For full D block
141 Teuchos::RCP< Teuchos::SerialDenseMatrix<ordinal_type, value_type> > DD, RR;
142 DD = Teuchos::rcp(new Teuchos::SerialDenseMatrix<ordinal_type, value_type> (Teuchos::Copy, D));
143
144 RR = Teuchos::rcp(new Teuchos::SerialDenseMatrix<ordinal_type, value_type> (r));
145
146
147 // Setup solver
148 solver.setMatrix(DD);
149 solver.setVectors(RR, RR);
150 //Solve D*r=r
151 if (solver.shouldEquilibrate()) {
152 solver.factorWithEquilibration(true);
153 solver.equilibrateMatrix();
154 }
155 solver.solve();
156
157
158 for (ordinal_type i=0; i<c-s; i++){
159 r(i,0)=(*RR)(i,0);
160 }
161
162 }
163
164 Teuchos::SerialDenseMatrix<ordinal_type, value_type> rm(Teuchos::Copy, rminus);
165 ret = rm.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS, -1.0, B, r, 1.0);
166 TEUCHOS_ASSERT(ret == 0);
167
168 //set r(l-1)=g(l-1)
169 if (l>lmin){
170 for (ordinal_type i=0; i<s; i++)
171 Resid(i,0)=rm(i,0);
172 }
173 else {
174 //For l=1, solve A0u0=g0
175 if (n<1){
176 U(0,0)=rm(0,0)/K(0,0);
177 }
178 //Compute the Schur completement
179 else {
180 for (ordinal_type i=0; i<s; i++)
181 Resid(i,0)=rm(i,0);
182
183 Teuchos::SerialDenseMatrix<ordinal_type, value_type> S(Teuchos::Copy, K, s, s);
184 Teuchos::SerialDenseMatrix<ordinal_type, value_type> BinvD(s,c-s);
185 for (ordinal_type i=0; i<c-s; i++)
186 for (ordinal_type j=0; j<s; j++)
187 BinvD(j,i)=B(j,i)/D(i,i);
188
189 S.multiply(Teuchos::NO_TRANS,Teuchos::TRANS, -1.0, BinvD, B, 1.0);
190 Teuchos::SerialDenseMatrix<ordinal_type, value_type> Ul(Teuchos::View, U, s, 1);
191 Teuchos::RCP< Teuchos::SerialDenseMatrix<ordinal_type, value_type> > SS, UU, RR;
192 SS = Teuchos::rcp(new Teuchos::SerialDenseMatrix<ordinal_type, value_type> (S));
193 UU = Teuchos::rcp(new Teuchos::SerialDenseMatrix<ordinal_type, value_type> (Teuchos::Copy, Ul));
194 RR = Teuchos::rcp(new Teuchos::SerialDenseMatrix<ordinal_type, value_type> (rm));
195 //Setup solver
196 Teuchos::SerialDenseSolver<ordinal_type, value_type> solver2;
197 solver2.setMatrix(SS);
198 solver2.setVectors(UU, RR);
199 //Solve S*u=rm
200 if (solver2.shouldEquilibrate()) {
201 solver2.factorWithEquilibration(true);
202 solver2.equilibrateMatrix();
203 }
204 solver2.solve();
205
206 for (ordinal_type i=0; i<s; i++)
207 U(i,0)=(*UU)(i,0);
208 }
209 }
210 }
211
212 for (ordinal_type l=lmin; l<=p; l++){
213 //compute post-correction
214 ordinal_type c = size(l,m);
215 ordinal_type s = size(l-1,m);
216
217 Teuchos::SerialDenseMatrix<ordinal_type, value_type> B(Teuchos::View, K, s, c-s, 0, s);
218 Teuchos::SerialDenseMatrix<ordinal_type, value_type> D(Teuchos::Copy, K, c-s, c-s, s,s);
219 //split residual
220 Teuchos::SerialDenseMatrix<ordinal_type, value_type> r(Teuchos::Copy, Resid, c-s, 1, s, 0);
221 //Get first s values of U
222 Teuchos::SerialDenseMatrix<ordinal_type, value_type> Ul(Teuchos::View, U, s, 1);
223 ret = r.multiply(Teuchos::TRANS,Teuchos::NO_TRANS, -1.0, B, Ul, 1.0);
224 if (diag == 0) {
225 //For diag D
226 //Concatenate U
227 for (ordinal_type i=s; i<c; i++)
228 U(i,0)=r(-s+i,0)/D(-s+i,-s+i);
229
230 }
231 else {
232 //For full block D
233
234 Teuchos::RCP< Teuchos::SerialDenseMatrix<ordinal_type, value_type> > DD, RR;
235 DD = Teuchos::rcp(new Teuchos::SerialDenseMatrix<ordinal_type, value_type> (D));
236 RR = Teuchos::rcp(new Teuchos::SerialDenseMatrix<ordinal_type, value_type> (r));
237 // Setup solver
238 solver.setMatrix(DD);
239 solver.setVectors(RR, RR);
240 //Solve D*r=r
241 if (solver.shouldEquilibrate()) {
242 solver.factorWithEquilibration(true);
243 solver.equilibrateMatrix();
244 }
245 solver.solve();
246 for (ordinal_type i=s; i<c; i++)
247 U(i,0)=(*RR)(-s+i,0);
248 }
249
250
251 }
252
253 return 0;
254}
255
256
ordinal_type fact(ordinal_type n) const
SchurPreconditioner(const Teuchos::SerialDenseMatrix< ordinal_type, value_type > &K, const ordinal_type p, const ordinal_type m, const ordinal_type diag)
Constructor.
ordinal_type size(ordinal_type n, ordinal_type m) const
virtual ordinal_type ApplyInverse(const Teuchos::SerialDenseMatrix< ordinal_type, value_type > &Input, Teuchos::SerialDenseMatrix< ordinal_type, value_type > &Result, ordinal_type prec_iters) const
Returns the result of a Operator inverse applied to a Teuchos::SerialDenseMatrix Input in Result.