44 #include <Epetra_Operator.h> 45 #include <Epetra_CrsMatrix.h> 46 #include <Epetra_LinearProblem.h> 47 #include <Amesos_BaseSolver.h> 48 #include <Ifpack_Preconditioner.h> 49 #include <Epetra_MultiVector.h> 51 #include "shylu_probing_operator.h" 55 ShyLU_Probing_Operator::ShyLU_Probing_Operator(
60 Epetra_LinearProblem *LP, Amesos_BaseSolver *solver,
61 Ifpack_Preconditioner *ifSolver, Epetra_CrsMatrix *C,
62 Epetra_Map *localDRowMap,
int nvectors)
72 localDRowMap_ = localDRowMap;
73 ResetTempVectors(nvectors);
81 matvec_time_ = RCP<Time>(
new Time(
"First 2 Matvec time in prober: "));
82 localize_time_ = RCP<Time>(
new Time(
"Localize time in prober: "));
83 trisolve_time_ = RCP<Time>(
new Time(
"Triangular solve time in prober: "));
84 dist_time_ = RCP<Time>(
new Time(
"Distribute time in prober: "));
85 matvec2_time_ = RCP<Time>(
new Time(
"Second 1 matvec time in prober: "));
86 apply_time_ = RCP<Time>(
new Time(
"Apply time in prober: "));
87 update_time_ = RCP<Time>(
new Time(
"Update time in prober: "));
91 int ShyLU_Probing_Operator::SetUseTranspose(
bool useTranspose)
96 int ShyLU_Probing_Operator::Apply(
const Epetra_MultiVector &X,
97 Epetra_MultiVector &Y)
const 100 apply_time_->start();
103 int nvectors = X.NumVectors();
104 bool local = (C_->Comm().NumProc() == 1);
109 matvec_time_->start();
112 err = G_->Multiply(
false, X, *temp2);
115 err = C_->Multiply(
false, X, *temp);
121 X.ExtractView(&values, &mylda);
123 Epetra_SerialComm LComm;
124 Epetra_Map SerialMap(X.Map().NumMyElements(), X.Map().NumMyElements(),
125 X.Map().MyGlobalElements(), 0, LComm);
126 Epetra_MultiVector Xl(View, SerialMap, values, mylda, X.NumVectors());
127 err = C_->Multiply(
false, Xl, *temp);
132 matvec_time_->stop();
135 int nrows = C_->RowMap().NumMyElements();
138 cout <<
"DEBUG MODE" << endl;
139 assert(nrows == localDRowMap_->NumGlobalElements());
141 int gids[nrows], gids1[nrows];
142 C_->RowMap().MyGlobalElements(gids);
143 localDRowMap_->MyGlobalElements(gids1);
145 for (
int i = 0; i < nrows; i++)
147 assert(gids[i] == gids1[i]);
152 localize_time_->start();
160 err = temp->ExtractView(&values, &lda);
164 assert(lda == nrows);
167 for (
int v = 0; v < nvectors; v++)
169 for (
int i = 0; i < nrows; i++)
171 err = ltemp->ReplaceMyValue(i, v, values[i+v*lda]);
178 localize_time_->stop();
179 trisolve_time_->start();
193 if (config_->amesosForDiagonal)
195 ssym_->OrigLP->SetLHS(localX.
getRawPtr());
197 ssym_->ReIdx_LP->fwd();
201 ifSolver_->ApplyInverse(*temp, *localX);
204 trisolve_time_->stop();
210 err = localX->ExtractView(&values, &lda);
215 for (
int v = 0; v < nvectors; v++)
217 for (
int i = 0; i < nrows; i++)
219 err = temp->ReplaceMyValue(i, v, values[i+v*lda]);
227 matvec2_time_->start();
232 R_->Multiply(
false, *temp, Y);
240 Y.ExtractView(&values, &mylda);
242 Epetra_SerialComm LComm;
243 Epetra_Map SerialMap(Y.Map().NumMyElements(), Y.Map().NumMyElements(),
244 Y.Map().MyGlobalElements(), 0, LComm);
245 Epetra_MultiVector Yl(View, SerialMap, values, mylda, Y.NumVectors());
246 R_->Multiply(
false, *localX, Yl);
250 matvec2_time_->stop();
251 update_time_->start();
254 err = Y.Update(1.0, *temp2, -1.0);
259 update_time_->stop();
266 void ShyLU_Probing_Operator::PrintTimingInfo()
269 cout << matvec_time_->name()<< matvec_time_->totalElapsedTime() << endl;
270 cout << localize_time_->name()<< localize_time_->totalElapsedTime() << endl;
271 cout << trisolve_time_->name()<< trisolve_time_->totalElapsedTime() << endl;
272 cout << dist_time_->name() <<dist_time_->totalElapsedTime() << endl;
273 cout << matvec2_time_->name() <<matvec2_time_->totalElapsedTime() << endl;
274 cout << update_time_->name() <<update_time_->totalElapsedTime() << endl;
275 cout << apply_time_->name() <<apply_time_->totalElapsedTime() << endl;
279 void ShyLU_Probing_Operator::ResetTempVectors(
int nvectors)
282 nvectors_ = nvectors;
283 if (nvectors_ == 0)
return;
285 if (nvectors <= ssym_->Drhs->NumVectors())
289 (
new Epetra_MultiVector(View, *(ssym_->Drhs), 0, nvectors));
291 (
new Epetra_MultiVector(View, *(ssym_->Dlhs), 0, nvectors));
293 (
new Epetra_MultiVector(View, *(ssym_->Gvec), 0, nvectors));
298 (
new Epetra_MultiVector((*(ssym_->Drhs)).Map(), nvectors));
300 (
new Epetra_MultiVector((*(ssym_->Dlhs)).Map(), nvectors));
302 (
new Epetra_MultiVector((*(ssym_->Gvec)).Map(), nvectors));
304 temp = RCP<Epetra_MultiVector>(
new Epetra_MultiVector(C_->RowMap(),
306 ssym_->OrigLP->SetLHS(localX.
getRawPtr());
308 ssym_->ReIdx_LP->fwd();
311 int ShyLU_Probing_Operator::ApplyInverse(
const Epetra_MultiVector &X,
312 Epetra_MultiVector &Y)
const 317 double ShyLU_Probing_Operator::NormInf()
const 322 const char* ShyLU_Probing_Operator::Label()
const 324 return "Shylu probing";
327 bool ShyLU_Probing_Operator::UseTranspose()
const 332 bool ShyLU_Probing_Operator::HasNormInf()
const 337 const Epetra_Comm& ShyLU_Probing_Operator::Comm()
const 342 const Epetra_Map& ShyLU_Probing_Operator::OperatorDomainMap()
const 345 return G_->DomainMap();
348 const Epetra_Map& ShyLU_Probing_Operator::OperatorRangeMap()
const 351 return G_->RangeMap();