44 #include "Teuchos_SerialDenseMatrix.hpp" 45 #include "Teuchos_SerialDenseSolver.hpp" 46 #include "Teuchos_RCP.hpp" 48 template <
typename ordinal_type,
typename value_type>
51 const Teuchos::SerialDenseMatrix<ordinal_type, value_type>& K_,
60 template <
typename ordinal_type,
typename value_type>
66 template <
typename ordinal_type,
typename value_type>
78 template <
typename ordinal_type,
typename value_type>
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,
111 Teuchos::SerialDenseMatrix<ordinal_type, value_type> Resid(Teuchos::Copy, Input);
119 Teuchos::SerialDenseSolver<ordinal_type, value_type> solver;
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);
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);
136 r(i,0)=r(i,0)/D(i,i);
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));
144 RR = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<ordinal_type, value_type> (r));
148 solver.setMatrix(DD);
149 solver.setVectors(RR, RR);
151 if (solver.shouldEquilibrate()) {
152 solver.factorWithEquilibration(
true);
153 solver.equilibrateMatrix();
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);
176 U(0,0)=rm(0,0)/K(0,0);
183 Teuchos::SerialDenseMatrix<ordinal_type, value_type> S(Teuchos::Copy, K, s, s);
184 Teuchos::SerialDenseMatrix<ordinal_type, value_type> BinvD(s,c-s);
187 BinvD(
j,i)=B(
j,i)/D(i,i);
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));
196 Teuchos::SerialDenseSolver<ordinal_type, value_type> solver2;
197 solver2.setMatrix(SS);
198 solver2.setVectors(UU, RR);
200 if (solver2.shouldEquilibrate()) {
201 solver2.factorWithEquilibration(
true);
202 solver2.equilibrateMatrix();
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);
220 Teuchos::SerialDenseMatrix<ordinal_type, value_type> r(Teuchos::Copy, Resid, c-s, 1, s, 0);
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);
228 U(i,0)=r(-s+i,0)/D(-s+i,-s+i);
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));
238 solver.setMatrix(DD);
239 solver.setVectors(RR, RR);
241 if (solver.shouldEquilibrate()) {
242 solver.factorWithEquilibration(
true);
243 solver.equilibrateMatrix();
247 U(i,0)=(*RR)(-s+i,0);
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.
ordinal_type size(ordinal_type n, ordinal_type m) const
ordinal_type fact(ordinal_type n) const
virtual ~SchurPreconditioner()
Destructor.
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...