56 static int shylu_dist_solve(
60 const Epetra_MultiVector& X,
66 assert(X.Map().SameAs(Y.Map()));
68 const Epetra_MultiVector *newX;
72 int nvectors = newX->NumVectors();
75 Epetra_Map BsMap(-1, data->Snr, data->SRowElems, 0, X.Comm());
76 Epetra_Map BdMap(-1, data->Dnr, data->DRowElems, 0, X.Comm());
78 Epetra_MultiVector Bs(BsMap, nvectors);
79 Epetra_Import BsImporter(BsMap, newX->Map());
81 assert(BsImporter.SourceMap().SameAs(newX->Map()));
82 assert((newX->Map()).SameAs(BsImporter.SourceMap()));
84 Bs.Import(*newX, BsImporter, Insert);
85 Epetra_MultiVector Xs(BsMap, nvectors);
87 Epetra_SerialComm LComm;
88 Epetra_Map LocalBdMap(-1, data->Dnr, data->DRowElems, 0, LComm);
89 Epetra_MultiVector localrhs(LocalBdMap, nvectors);
90 Epetra_MultiVector locallhs(LocalBdMap, nvectors);
92 Epetra_MultiVector Z(BdMap, nvectors);
94 Epetra_MultiVector Bd(BdMap, nvectors);
95 Epetra_Import BdImporter(BdMap, newX->Map());
96 assert(BdImporter.SourceMap().SameAs(newX->Map()));
97 assert((newX->Map()).SameAs(BdImporter.SourceMap()));
98 Bd.Import(*newX, BdImporter, Insert);
102 err = Bd.ExtractView(&values, &lda);
104 int nrows = ssym->C->RowMap().NumMyElements();
107 assert(lda == nrows);
108 for (
int v = 0; v < nvectors; v++)
110 for (
int i = 0; i < nrows; i++)
112 err = localrhs.ReplaceMyValue(i, v, values[i+v*lda]);
118 if (config->amesosForDiagonal)
120 ssym->LP->SetRHS(&localrhs);
121 ssym->LP->SetLHS(&locallhs);
122 ssym->Solver->Solve();
126 ssym->ifSolver->ApplyInverse(localrhs, locallhs);
129 err = locallhs.ExtractView(&values, &lda);
133 assert(lda == nrows);
134 for (
int v = 0; v < nvectors; v++)
136 for (
int i = 0; i < nrows; i++)
138 err = Z.ReplaceMyValue(i, v, values[i+v*lda]);
143 Epetra_MultiVector temp1(BsMap, nvectors);
144 ssym->R->Multiply(
false, Z, temp1);
145 Bs.Update(-1.0, temp1, 1.0);
149 Epetra_LinearProblem Problem(data->Sbar.get(), &Xs, &Bs);
150 if (config->schurSolver ==
"Amesos")
152 Amesos_BaseSolver *solver2 = data->dsolver;
153 data->LP2->SetLHS(&Xs);
154 data->LP2->SetRHS(&Bs);
161 if (config->libName ==
"Belos")
163 solver = data->innersolver;
171 solver =
new AztecOO();
173 solver->SetAztecOption(AZ_solver, AZ_gmres);
175 solver->SetAztecOption(AZ_precond, AZ_dom_decomp);
183 solver->SetProblem(Problem);
187 solver->Iterate(config->inner_maxiters, config->inner_tolerance);
190 Epetra_MultiVector temp(BdMap, nvectors);
191 ssym->C->Multiply(
false, Xs, temp);
192 temp.Update(1.0, Bd, -1.0);
201 err = temp.ExtractView(&values, &lda);
206 assert(lda == nrows);
207 for (
int v = 0; v < nvectors; v++)
209 for (
int i = 0; i < nrows; i++)
211 err = localrhs.ReplaceMyValue(i, v, values[i+v*lda]);
216 if (config->amesosForDiagonal)
218 ssym->LP->SetRHS(&localrhs);
219 ssym->LP->SetLHS(&locallhs);
220 ssym->Solver->Solve();
224 ssym->ifSolver->ApplyInverse(localrhs, locallhs);
227 err = locallhs.ExtractView(&values, &lda);
231 assert(lda == nrows);
232 for (
int v = 0; v < nvectors; v++)
234 for (
int i = 0; i < nrows; i++)
236 err = temp.ReplaceMyValue(i, v, values[i+v*lda]);
244 Epetra_Export XdExporter(BdMap, Y.Map());
245 Y.Export(temp, XdExporter, Insert);
247 Epetra_Export XsExporter(BsMap, Y.Map());
248 Y.Export(Xs, XsExporter, Insert);
250 if (config->libName ==
"Belos" || config->schurSolver ==
"Amesos")
261 static int shylu_local_solve(
265 const Epetra_MultiVector& X,
266 Epetra_MultiVector& Y
271 int nvectors = X.NumVectors();
272 assert (nvectors == data->localrhs->NumVectors());
276 data->Xs->PutScalar(0.0);
279 data->localrhs->Import(X, *(data->BdImporter), Insert);
282 if (config->amesosForDiagonal) {
283 ssym->OrigLP->SetRHS((data->localrhs).getRawPtr());
284 ssym->OrigLP->SetLHS((data->locallhs).getRawPtr());
285 ssym->ReIdx_LP->fwd();
286 ssym->Solver->Solve();
289 ssym->ifSolver->ApplyInverse(*(data->localrhs), *(data->locallhs));
292 err = ssym->R->Multiply(
false, *(data->locallhs), *(data->temp1));
296 data->temp2->Import(*(data->temp1), *(data->DistImporter), Insert);
299 data->Bs->Import(X, *(data->BsImporter), Insert);
300 data->Bs->Update(-1.0, *(data->temp2), 1.0);
303 Epetra_LinearProblem Problem(data->Sbar.get(),
304 (data->Xs).getRawPtr(), (data->Bs).getRawPtr());
305 if ((config->schurSolver ==
"G") || (config->schurSolver ==
"IQR"))
307 IFPACK_CHK_ERR(data->iqrSolver->Solve(*(data->schur_op),
308 *(data->Bs), *(data->Xs)));
310 else if (config->schurSolver ==
"Amesos")
312 Amesos_BaseSolver *solver2 = data->dsolver;
313 data->OrigLP2->SetLHS((data->Xs).getRawPtr());
314 data->OrigLP2->SetRHS((data->Bs).getRawPtr());
315 data->ReIdx_LP2->fwd();
322 if (config->libName ==
"Belos")
324 solver = data->innersolver;
325 solver->SetLHS((data->Xs).getRawPtr());
326 solver->SetRHS((data->Bs).getRawPtr());
332 solver =
new AztecOO();
334 solver->SetAztecOption(AZ_solver, AZ_gmres);
336 solver->SetAztecOption(AZ_precond, AZ_dom_decomp);
344 solver->SetProblem(Problem);
348 solver->Iterate(config->inner_maxiters, config->inner_tolerance);
352 data->LocalXs->Import(*(data->Xs), *(data->XsImporter), Insert);
354 err = ssym->C->Multiply(
false, *(data->LocalXs), *(data->temp3));
356 data->temp3->Update(1.0, *(data->localrhs), -1.0);
358 if (config->amesosForDiagonal) {
359 ssym->OrigLP->SetRHS((data->temp3).getRawPtr());
360 ssym->OrigLP->SetLHS((data->locallhs).getRawPtr());
361 ssym->ReIdx_LP->fwd();
362 ssym->Solver->Solve();
365 ssym->ifSolver->ApplyInverse(*(data->temp3), *(data->locallhs));
368 Y.Export(*(data->locallhs), *(data->XdExporter), Insert);
369 Y.Export(*(data->LocalXs), *(data->XsExporter), Insert);
371 if (config->libName ==
"Belos" || config->schurSolver ==
"Amesos")
386 const Epetra_MultiVector& X,
387 Epetra_MultiVector& Y
390 if (config->sep_type != 1)
391 shylu_dist_solve(ssym, data, config, X, Y);
393 shylu_local_solve(ssym, data, config, X, Y);
int shylu_solve(shylu_symbolic *ssym, shylu_data *data, shylu_config *config, const Epetra_MultiVector &X, Epetra_MultiVector &Y)
Call solve on multiple RHS.
Main data structure holding needed offset and temp variables.
Main header file of ShyLU (Include main user calls)