46 #ifdef HAVE_STOKHOS_ML 51 Stokhos::MLPrecOp::MLPrecOp(
const Epetra_CrsMatrix& mean_op,
const Teuchos::Array<double>& norms,
const Epetra_Comm& Comm,
const Epetra_Map& DMap,
const Epetra_Map& RMap)
59 Label_ =
"Stochos::MLPrecOp";
63 Teuchos::ParameterList MLList;
70 ML_Epetra::SetDefaults(
"SA",MLList);
78 MLList.set(
"ML output", 10);
80 MLList.set(
"max levels",5);
82 MLList.set(
"increasing or decreasing",
"increasing");
85 MLList.set(
"aggregation: type",
"Uncoupled");
91 MLList.set(
"smoother: type",
"Chebyshev");
92 MLList.set(
"smoother: sweeps",3);
95 MLList.set(
"smoother: pre or post",
"both");
100 MLList.set(
"coarse: type",
"Amesos-KLU");
105 MLList.set(
"coarse: type",
"Jacobi");
113 MLPrec =
new ML_Epetra::MultiLevelPreconditioner(mean_op, MLList);
117 MLPrec->PrintUnused(0);
126 int Stokhos::MLPrecOp::ApplyInverse(
const Epetra_MultiVector& X, Epetra_MultiVector& Y)
const{
128 if (!X.Map().SameAs(OperatorDomainMap()))
129 std::cout <<
"!X.Map().SameAs(OperatorDomainMap())\n";
130 if (!Y.Map().SameAs(OperatorRangeMap()))
131 std::cout<<
"!Y.Map().SameAs(OperatorRangeMap())\n";
132 if (Y.NumVectors()!=X.NumVectors())
133 std::cout<<
"Y.NumVectors()!=X.NumVectors()\n";
137 for(
int mm = 0; mm< X.NumVectors(); mm++){
140 int N_xi = norms_.size();
141 int N_x = X.MyLength()/N_xi;
145 Epetra_Map Map(MLPrec->OperatorDomainMap());
148 Epetra_MultiVector xBlock(Map,N_xi);
149 Epetra_MultiVector yBlock(MLPrec->OperatorRangeMap(),N_xi);
152 int MyLength = xBlock.MyLength();
153 for(
int c=0; c<N_xi ; c++){
154 for(
int i=0; i<MyLength; i++){
155 xBlock[c][i] = (X)[mm][c*N_x + i];
160 Epetra_MultiVector blockProducts(MLPrec->OperatorRangeMap(),N_xi);
161 MLPrec->ApplyInverse(xBlock, blockProducts);
164 for(
int j = 0;
j<N_xi;
j++){
165 (*yBlock(
j)).Update(1/norms_[
j],*blockProducts(
j), 1.0);
168 for(
int c=0; c<N_xi ; c++){
169 for(
int i=0; i<MyLength; i++){
170 (Y)[mm][c*N_x + i] = yBlock[c][i];