Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Stokhos_BlockPreconditionerImp.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
47//Computes the exact Schur complement block LU decomposition
48
49template <typename ordinal_type, typename value_type>
52 const Teuchos::SerialDenseMatrix<ordinal_type, value_type>& K_,const ordinal_type p_, const ordinal_type m_) :
53 K(K_),
54 p(p_),
55 m(m_)
56{
57}
58
59template <typename ordinal_type, typename value_type>
62{
63}
64
65template <typename ordinal_type, typename value_type>
66ordinal_type
68facto(ordinal_type n) const
69{
70 if (n > 1)
71 return (n * facto(n-1));
72 else
73 return (1);
74}
75
76template <typename ordinal_type, typename value_type>
77ordinal_type
79siz (ordinal_type n, ordinal_type m) const
80{
81 //n is the polynomial order and m is the number of random variables
82 return (facto(n+m)/(facto(n)*facto(m)));
83 }
84
85template <typename ordinal_type, typename value_type>
86ordinal_type
88ApplyInverse(const Teuchos::SerialDenseMatrix<ordinal_type, value_type>& Input,
89 Teuchos::SerialDenseMatrix<ordinal_type, value_type>& Result,
90 ordinal_type n) const
91{ //Solve M*Result=Input
92 ordinal_type c=siz(p,m);
93 ordinal_type s = siz(p-1,m);
94
95 //Split residual
96 Teuchos::SerialDenseMatrix<ordinal_type, value_type> r1(Teuchos::Copy, Input, s, 1);
97 Teuchos::SerialDenseMatrix<ordinal_type, value_type> r2(Teuchos::Copy, Input, c-s, 1, s, 0);
98
99 //Split Result
100 Teuchos::SerialDenseMatrix<ordinal_type, value_type> u1(Teuchos::Copy, Result, s, 1);
101 Teuchos::SerialDenseMatrix<ordinal_type, value_type> u2(Teuchos::Copy, Result, c-s, 1, s, 0);
102
103 Teuchos::SerialDenseMatrix<ordinal_type, value_type> B(Teuchos::View, K, s, c-s, 0, s);
104 Teuchos::SerialDenseMatrix<ordinal_type, value_type> D(Teuchos::View, K, c-s, c-s, s,s);
105
106 //rD=inv(D)r2
107
108 Teuchos::SerialDenseMatrix<ordinal_type, value_type> Dr(c-s,1);
109
110 for (ordinal_type i=0; i<c-s; i++)
111 Dr(i,0)=r2(i,0)/D(i,i);
112
113 ordinal_type ret = r1.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS, -1.0, B, Dr, 1.0);
114 TEUCHOS_ASSERT(ret == 0);
115
116 //Compute S=A-B*inv(D)*Bt
117 Teuchos::SerialDenseMatrix<ordinal_type, value_type> S(Teuchos::Copy, K, s, s);
118 //Compute B*inv(D)
119 Teuchos::SerialDenseMatrix<ordinal_type, value_type> BinvD(s,c-s);
120 for (ordinal_type i=0; i<c-s; i++) //col
121 for (ordinal_type j=0; j<s; j++) //row
122 BinvD(j,i)=B(j,i)/D(i,i);
123
124 S.multiply(Teuchos::NO_TRANS,Teuchos::TRANS, -1.0, BinvD, B, 1.0);
125
126 Teuchos::RCP< Teuchos::SerialDenseMatrix<ordinal_type, value_type> > SS, w, rr;
127 SS = Teuchos::rcp(new Teuchos::SerialDenseMatrix<ordinal_type, value_type> (S));
128 w = Teuchos::rcp(new Teuchos::SerialDenseMatrix<ordinal_type, value_type> (s,1));
129 rr = Teuchos::rcp(new Teuchos::SerialDenseMatrix<ordinal_type, value_type> (r1));
130
131
132 // Setup solver
133 Teuchos::SerialDenseSolver<ordinal_type, value_type> solver;
134 solver.setMatrix(SS);
135 solver.setVectors(w, rr);
136 //Solve S*w=r1
137 if (solver.shouldEquilibrate()) {
138 solver.factorWithEquilibration(true);
139 solver.equilibrateMatrix();
140 }
141 solver.solve();
142
143 for (ordinal_type i=0; i<s; i++)
144 Result(i,0)=(*w)(i,0);
145
146 ret = r2.multiply(Teuchos::TRANS,Teuchos::NO_TRANS, -1.0, B, *w, 1.0);
147 TEUCHOS_ASSERT(ret == 0);
148
149 for (ordinal_type i=s; i<c; i++)
150 Result(i,0)=r2(-s+i,0)/D(-s+i, -s+i);
151
152 return 0;
153}
ordinal_type facto(ordinal_type n) const
ordinal_type siz(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 m) const
Returns the result of a Operator inverse applied to a Teuchos::SerialDenseMatrix Input in Result.
BlockPreconditioner(const Teuchos::SerialDenseMatrix< ordinal_type, value_type > &K, const ordinal_type p, const ordinal_type m)
Constructor.