44 #include "Teuchos_SerialDenseMatrix.hpp" 45 #include "Teuchos_SerialDenseSolver.hpp" 49 template <
typename ordinal_type,
typename value_type>
59 template <
typename ordinal_type,
typename value_type>
65 template <
typename ordinal_type,
typename value_type>
71 return (n * facto(n-1));
76 template <
typename ordinal_type,
typename value_type>
82 return (facto(n+m)/(facto(n)*facto(m)));
85 template <
typename ordinal_type,
typename value_type>
88 ApplyInverse(
const Teuchos::SerialDenseMatrix<ordinal_type, value_type>& Input,
89 Teuchos::SerialDenseMatrix<ordinal_type, value_type>& Result,
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);
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);
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);
108 Teuchos::SerialDenseMatrix<ordinal_type, value_type> Dr(c-s,1);
111 Dr(i,0)=r2(i,0)/D(i,i);
113 ordinal_type ret = r1.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS, -1.0, B, Dr, 1.0);
114 TEUCHOS_ASSERT(ret == 0);
117 Teuchos::SerialDenseMatrix<ordinal_type, value_type> S(Teuchos::Copy, K, s, s);
119 Teuchos::SerialDenseMatrix<ordinal_type, value_type> BinvD(s,c-s);
122 BinvD(
j,i)=B(
j,i)/D(i,i);
124 S.multiply(Teuchos::NO_TRANS,Teuchos::TRANS, -1.0, BinvD, B, 1.0);
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));
133 Teuchos::SerialDenseSolver<ordinal_type, value_type> solver;
134 solver.setMatrix(SS);
135 solver.setVectors(w, rr);
137 if (solver.shouldEquilibrate()) {
138 solver.factorWithEquilibration(
true);
139 solver.equilibrateMatrix();
144 Result(i,0)=(*w)(i,0);
146 ret = r2.multiply(Teuchos::TRANS,Teuchos::NO_TRANS, -1.0, B, *w, 1.0);
147 TEUCHOS_ASSERT(ret == 0);
150 Result(i,0)=r2(-s+i,0)/D(-s+i, -s+i);
virtual ~BlockPreconditioner()
Destructor.
ordinal_type facto(ordinal_type n) const
BlockPreconditioner(const Teuchos::SerialDenseMatrix< ordinal_type, value_type > &K, const ordinal_type p, const ordinal_type m)
Constructor.
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...