171 const Teuchos::ParameterList& pL = GetParameterList();
173 bool M_is_diagonal =
false;
174 if(pL.isParameter(
"rap: shift diagonal M"))
175 M_is_diagonal = pL.get<
bool>(
"rap: shift diagonal M");
178 bool use_low_storage =
false;
179 if(pL.isParameter(
"rap: shift low storage")) {
180 use_low_storage = pL.get<
bool>(
"rap: shift low storage");
181 M_is_diagonal = use_low_storage ? true : M_is_diagonal;
184 Teuchos::ArrayView<const double> doubleShifts;
185 Teuchos::ArrayRCP<double> myshifts;
186 if(pL.isParameter(
"rap: shift array") && pL.get<Teuchos::Array<double> >(
"rap: shift array").size() > 0 ) {
188 doubleShifts = pL.get<Teuchos::Array<double> >(
"rap: shift array")();
190 if(pL.isParameter(
"rap: cfl array") && pL.get<Teuchos::Array<double> >(
"rap: cfl array").size() > 0) {
192 Teuchos::ArrayView<const double> CFLs = pL.get<Teuchos::Array<double> >(
"rap: cfl array")();
194 double dt = Get<double>(fineLevel,
"deltaT");
195 double cfl = Get<double>(fineLevel,
"cfl");
196 double ts_at_cfl1 = dt / cfl;
197 myshifts.resize(CFLs.size());
198 Teuchos::Array<double> myCFLs(CFLs.size());
202 for(
int i=1; i<(int)CFLs.size(); i++)
203 myCFLs[i] = (CFLs[i]> cfl) ? cfl : CFLs[i];
206 std::ostringstream ofs;
207 ofs<<
"RAPShiftFactory: CFL schedule = ";
208 for(
int i=0; i<(int)CFLs.size(); i++)
212 GetOStream(
Statistics0)<<
"RAPShiftFactory: Timestep at CFL=1 is "<< ts_at_cfl1 <<
" " <<std::endl;
215 for(
int i=0; i<(int)myshifts.size(); i++)
216 myshifts[i] = 1.0 / (ts_at_cfl1*myCFLs[i]);
217 doubleShifts = myshifts();
220 std::ostringstream ofs;
221 ofs<<
"RAPShiftFactory: shift schedule = ";
222 for(
int i=0; i<(int)doubleShifts.size(); i++)
223 ofs<<
" "<<doubleShifts[i];
226 Set(coarseLevel,
"cfl-based shift array",myshifts);
229 myshifts = Get<Teuchos::ArrayRCP<double> > (fineLevel,
"cfl-based shift array");
230 doubleShifts = myshifts();
231 Set(coarseLevel,
"cfl-based shift array",myshifts);
242 if(use_low_storage) K = Get< RCP<Matrix> >(fineLevel,
"A");
243 else K = Get< RCP<Matrix> >(fineLevel,
"K");
244 if(!M_is_diagonal) M = Get< RCP<Matrix> >(fineLevel,
"M");
245 else Mdiag = Get< RCP<Vector> >(fineLevel,
"Mdiag");
247 RCP<Matrix> P = Get< RCP<Matrix> >(coarseLevel,
"P");
261 KP = Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*K,
false, *P,
false, KP, GetOStream(
Statistics2));
263 MP = Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*M,
false, *P,
false, MP, GetOStream(
Statistics2));
266 MP = Xpetra::MatrixFactory2<Scalar, LocalOrdinal, GlobalOrdinal, Node>::BuildCopy(P);
267 MP->leftScale(*Mdiag);
270 Set(coarseLevel,
"AP Pattern", KP);
273 bool doOptimizedStorage =
true;
275 RCP<Matrix> Ac, Kc, Mc;
281 bool doFillComplete=
true;
282 if (implicitTranspose_) {
284 Kc = Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*P,
true, *KP,
false, Kc, GetOStream(
Statistics2), doFillComplete, doOptimizedStorage);
285 Mc = Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*P,
true, *MP,
false, Mc, GetOStream(
Statistics2), doFillComplete, doOptimizedStorage);
288 RCP<Matrix> R = Get< RCP<Matrix> >(coarseLevel,
"R");
290 Kc = Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*R,
false, *KP,
false, Kc, GetOStream(
Statistics2), doFillComplete, doOptimizedStorage);
291 Mc = Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*R,
false, *MP,
false, Mc, GetOStream(
Statistics2), doFillComplete, doOptimizedStorage);
298 int level = coarseLevel.GetLevelID();
299 Scalar shift = Teuchos::ScalarTraits<Scalar>::zero();
300 if(!use_low_storage) {
302 if(level < (
int)shifts_.size()) shift = shifts_[level];
303 else shift = Teuchos::as<Scalar>(pL.get<
double>(
"rap: shift"));
307 if(level < (
int)shifts_.size()) {
308 if(level==1) shift = shifts_[level];
310 Scalar prod1 = Teuchos::ScalarTraits<Scalar>::one();
311 for(
int i=1; i < level-1; i++) {
314 shift = (prod1 * shifts_[level] - prod1);
317 else if(doubleShifts.size() != 0) {
318 double d_shift = 0.0;
319 if(level < doubleShifts.size())
320 d_shift = doubleShifts[level] - doubleShifts[level-1];
323 GetOStream(
Warnings1) <<
"WARNING: RAPShiftFactory has detected a negative shift... This implies a less stable coarse grid."<<std::endl;
324 shift = Teuchos::as<Scalar>(d_shift);
327 double base_shift = pL.get<
double>(
"rap: shift");
328 if(level == 1) shift = Teuchos::as<Scalar>(base_shift);
329 else shift = Teuchos::as<Scalar>(pow(base_shift,level) - pow(base_shift,level-1));
332 GetOStream(
Runtime0) <<
"RAPShiftFactory: Using shift " << shift << std::endl;
338 Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::TwoMatrixAdd(*Kc,
false, Teuchos::ScalarTraits<Scalar>::one(), *Mc,
false, shift, Ac, GetOStream(
Statistics2));
342 Teuchos::ArrayView<const double> relativeFloor = pL.get<Teuchos::Array<double> >(
"rap: relative diagonal floor")();
343 if(relativeFloor.size() > 0)
344 Xpetra::MatrixUtils<SC,LO,GO,NO>::RelativeDiagonalBoost(Ac, relativeFloor,GetOStream(
Statistics2));
347 bool repairZeroDiagonals = pL.get<
bool>(
"RepairMainDiagonal") || pL.get<
bool>(
"rap: fix zero diagonals");
348 bool checkAc = pL.get<
bool>(
"CheckMainDiagonal")|| pL.get<
bool>(
"rap: fix zero diagonals"); ;
349 if (checkAc || repairZeroDiagonals)
350 Xpetra::MatrixUtils<SC,LO,GO,NO>::CheckRepairMainDiagonal(Ac, repairZeroDiagonals, GetOStream(
Warnings1));
352 RCP<ParameterList> params = rcp(
new ParameterList());;
353 params->set(
"printLoadBalancingInfo",
true);
356 Set(coarseLevel,
"A", Ac);
359 Set(coarseLevel,
"K", Kc);
362 Set(coarseLevel,
"M", Mc);
367 RCP<Vector> Mcv = Xpetra::VectorFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(Mc->getRowMap(),
false);
368 Mc->getLocalDiagCopy(*Mcv);
369 Set(coarseLevel,
"Mdiag", Mcv);
375 if (transferFacts_.begin() != transferFacts_.end()) {
379 for (std::vector<RCP<const FactoryBase> >::const_iterator it = transferFacts_.begin(); it != transferFacts_.end(); ++it) {
380 RCP<const FactoryBase> fac = *it;
381 GetOStream(
Runtime0) <<
"RAPShiftFactory: call transfer factory: " << fac->description() << std::endl;
382 fac->CallBuild(coarseLevel);