44#include "Teuchos_SerialDenseMatrix.hpp"
45#include "Teuchos_SerialDenseSolver.hpp"
49template <
typename ordinal_type,
typename value_type>
52 const Teuchos::SerialDenseMatrix<ordinal_type, value_type>& K_,
const ordinal_type p_,
const ordinal_type m_) :
59template <
typename ordinal_type,
typename value_type>
65template <
typename ordinal_type,
typename value_type>
68facto(ordinal_type n)
const
71 return (n * facto(n-1));
76template <
typename ordinal_type,
typename value_type>
79siz (ordinal_type n, ordinal_type m)
const
82 return (facto(n+m)/(facto(n)*facto(m)));
85template <
typename ordinal_type,
typename value_type>
88ApplyInverse(
const Teuchos::SerialDenseMatrix<ordinal_type, value_type>& Input,
89 Teuchos::SerialDenseMatrix<ordinal_type, value_type>& Result,
92 ordinal_type c=siz(p,m);
93 ordinal_type s = siz(p-1,m);
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);
110 for (ordinal_type i=0; i<c-s; i++)
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);
120 for (ordinal_type i=0; i<c-s; i++)
121 for (ordinal_type
j=0;
j<s;
j++)
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();
143 for (ordinal_type i=0; i<s; i++)
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);
149 for (ordinal_type i=s; i<c; i++)
150 Result(i,0)=r2(-s+i,0)/D(-s+i, -s+i);
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.
virtual ~BlockPreconditioner()
Destructor.