106ApplyInverse(
const Teuchos::SerialDenseMatrix<ordinal_type, value_type>& Input,
107 Teuchos::SerialDenseMatrix<ordinal_type, value_type>& U,
108 const ordinal_type n)
const
111 Teuchos::SerialDenseMatrix<ordinal_type, value_type> Resid(Teuchos::Copy, Input);
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);
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);
135 for (ordinal_type i=0; i<c-s; i++)
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();
158 for (ordinal_type i=0; i<c-s; i++){
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);
170 for (ordinal_type i=0; i<s; i++)
176 U(0,0)=rm(0,0)/K(0,0);
180 for (ordinal_type i=0; i<s; i++)
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);
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();
206 for (ordinal_type i=0; i<s; i++)
212 for (ordinal_type l=lmin; l<=p; l++){
214 ordinal_type c = size(l,m);
215 ordinal_type s = size(l-1,m);
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);
227 for (ordinal_type i=s; i<c; i++)
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();
246 for (ordinal_type i=s; i<c; i++)
247 U(i,0)=(*RR)(-s+i,0);