Stokhos Package Browser (Single Doxygen Collection)  Version of the Day
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 
48 template <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 
60 template <typename ordinal_type, typename value_type>
63 {
64 }
65 
66 template <typename ordinal_type, typename value_type>
69 fact (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 
78 template <typename ordinal_type, typename value_type>
82 {
83  //n is the polynomial order and m is the number of random variables
84  // return (fact(n+m)/(fact(n)*fact(m)));
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 
103 template <typename ordinal_type, typename value_type>
106 ApplyInverse(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 
KOKKOS_INLINE_FUNCTION PCE< Storage > min(const typename PCE< Storage >::value_type &a, const PCE< Storage > &b)
SchurPreconditioner(const Teuchos::SerialDenseMatrix< ordinal_type, value_type > &K, const ordinal_type p, const ordinal_type m, const ordinal_type diag)
Constructor.
expr expr expr expr j
ordinal_type size(ordinal_type n, ordinal_type m) const
ordinal_type fact(ordinal_type n) 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...