62 #define BLOCK_DIAGONAL_Si 66 #ifdef HAVE_SHYLUCORE_MPI 67 #include <Epetra_MpiComm.h> 69 #include <Epetra_SerialComm.h> 71 #include <EpetraExt_Reindex_LinearProblem2.h> 73 #include "Ifpack_config.h" 75 #ifdef HAVE_IFPACK_DYNAMIC_FACTORY 76 #include "Ifpack_DynamicFactory.h" 81 #include "shylu_internal_gmres.h" 82 #include "shylu_internal_gmres_tools.h" 95 int *DRowElems = data->DRowElems;
96 int *SRowElems = data->SRowElems;
97 int *DColElems = data->DColElems;
98 int *gvals = data->gvals;
100 double Sdiagfactor = config->Sdiagfactor;
101 int sym = config->sym;
105 #ifdef HAVE_SHYLUCORE_MPI 106 Epetra_MpiComm LComm(MPI_COMM_SELF);
108 Epetra_SerialComm LComm;
109 #endif // HAVE_SHYLUCORE_MPI 112 Epetra_Map LocalDRowMap(-1, Dnr, DRowElems, 0, LComm);
113 Epetra_Map DRowMap(-1, Dnr, DRowElems, 0, A->Comm());
114 Epetra_Map SRowMap(-1, Snr, SRowElems, 0, A->Comm());
115 Epetra_Map LocalSRowMap(-1, Snr, SRowElems, 0, LComm);
118 Epetra_Map LocalDColMap(-1, Dnc, DColElems, 0, LComm);
124 int *DNumEntriesPerRow =
new int[Dnr];
125 int *GNumEntriesPerRow =
new int[Snr];
126 int *CNumEntriesPerRow =
new int[Dnr];
127 int *RNumEntriesPerRow =
new int[Snr];
128 int *SBarNumEntriesPerRow =
new int[Snr];
129 int Dmaxnnz=0, Cmaxnnz=0, Rmaxnnz=0, Smaxnnz=0;
135 int Sdiag = (int) Snr * Sdiagfactor;
136 Sdiag = MIN(Sdiag, Snr-1);
137 Sdiag = MAX(Sdiag, 0);
140 if (Snr != 0) assert (Sdiag <= Snr-1);
142 int dcol, ccol, rcol, scol;
146 int nrows = A->RowMap().NumMyElements();
147 int *rows = A->RowMap().MyGlobalElements();
152 int dcnt, scnt, ccnt, rcnt;
153 dcol = 0 ; ccol = 0; rcol = 0; scol = 0;
156 for (
int i = 0; i < nrows; i++)
158 dcnt = 0; scnt = 0; ccnt = 0; rcnt = 0;
160 int err = A->ExtractMyRowView(i, NumEntries, Ax, Ai);
170 for (
int j = 0 ; j < NumEntries ; j++)
172 gcid = A->GCID(Ai[j]);
174 if (A->LRID(gcid) != -1 && gvals[gcid] == 1)
185 DNumEntriesPerRow[dcol++] = dcnt;
186 CNumEntriesPerRow[ccol++] = ccnt;
192 for (
int j = 0 ; j < NumEntries ; j++)
194 gcid = A->GCID(Ai[j]);
198 if (gvals[gcid] == 1)
211 GNumEntriesPerRow[scol++] = scnt;
212 RNumEntriesPerRow[rcol++] = rcnt;
214 Dmaxnnz = max(Dmaxnnz, dcnt);
215 Smaxnnz = max(Smaxnnz, scnt);
216 Rmaxnnz = max(Rmaxnnz, rcnt);
217 Cmaxnnz = max(Cmaxnnz, ccnt);
219 assert( dcol == Dnr);
220 assert( scol == Snr);
221 assert( ccol == Dnr);
222 assert( rcol == Snr);
224 for (
int i = 0; i < nrows; i++)
231 assert(sbarcol < Snr);
232 SBarNumEntriesPerRow[sbarcol] = GNumEntriesPerRow[sbarcol]
239 assert( sbarcol == Snr);
240 Smaxnnz = Smaxnnz + Sdiag * 2 + 1;
253 LocalDRowMap, LocalDColMap, DNumEntriesPerRow,
true));
258 SRowMap, GNumEntriesPerRow,
true));
260 if (config->sep_type != 1)
264 DRowMap, CNumEntriesPerRow,
true));
268 SRowMap, RNumEntriesPerRow,
true));
274 LocalDRowMap, CNumEntriesPerRow,
true));
278 LocalSRowMap, RNumEntriesPerRow,
true));
282 if (config->schurApproxMethod == 1)
286 SRowMap, SBarNumEntriesPerRow,
false));
290 data->lmax = max(Dmaxnnz, Rmaxnnz);
291 data->rmax = max(Cmaxnnz, Smaxnnz);
293 delete[] DNumEntriesPerRow;
294 delete[] GNumEntriesPerRow;
295 delete[] CNumEntriesPerRow;
296 delete[] RNumEntriesPerRow;
297 delete[] SBarNumEntriesPerRow;
317 int *DColElems = data->DColElems;
318 int *gvals = data->gvals;
319 double Sdiagfactor = config->Sdiagfactor;
321 int *LeftIndex =
new int[data->lmax];
322 double *LeftValues =
new double[data->lmax];
323 int *RightIndex =
new int[data->rmax];
324 double *RightValues =
new double[data->rmax];
332 int nrows = A->RowMap().NumMyElements();
333 int *rows = A->RowMap().MyGlobalElements();
335 for (
int i = 0; i < nrows ; i++)
338 err = A->ExtractMyRowView(i, NumEntries, Ax, Ai);
344 for (
int j = 0 ; j < NumEntries ; j++)
347 gcid = A->GCID(Ai[j]);
350 if ((gvals[gid] != 1 && gvals[gcid] == 1)
351 || (gvals[gid] == 1 && A->LRID(gcid) != -1 && gvals[gcid] == 1))
353 assert(lcnt < data->lmax);
355 LeftIndex[lcnt] = gcid;
359 lcid = (gvals[gid] == 1 ? D->LCID(gcid) : R->LCID(gcid));
361 LeftIndex[lcnt] = lcid;
363 LeftValues[lcnt++] = Ax[j];
367 assert(rcnt < data->rmax);
369 RightIndex[rcnt] = gcid;
373 lcid = (gvals[gid] == 1 ? C->LCID(gcid) : G->LCID(gcid));
375 RightIndex[rcnt] = lcid;
377 RightValues[rcnt++] = Ax[j];
385 err = D->InsertGlobalValues(gid, lcnt, LeftValues, LeftIndex);
387 err = C->InsertGlobalValues(gid, rcnt, RightValues, RightIndex);
392 err = D->ReplaceMyValues(D->LRID(gid), lcnt, LeftValues,
395 err = C->ReplaceMyValues(C->LRID(gid), rcnt, RightValues,
406 err = R->InsertGlobalValues(gid, lcnt, LeftValues, LeftIndex);
408 err = G->InsertGlobalValues(gid, rcnt, RightValues, RightIndex);
410 if (config->schurApproxMethod == 1)
412 err = Sg->InsertGlobalIndices(gid, rcnt, RightIndex);
419 err = R->ReplaceMyValues(R->LRID(gid), lcnt, LeftValues,
422 err = G->ReplaceMyValues(G->LRID(gid), rcnt, RightValues,
432 Epetra_Map *DRowMap, *SRowMap, *DColMap;
433 Epetra_SerialComm LComm;
434 if (config->sep_type != 1)
436 DRowMap =
new Epetra_Map(-1, data->Dnr, data->DRowElems, 0,
438 SRowMap =
new Epetra_Map(-1, data->Snr, data->SRowElems, 0,
440 DColMap =
new Epetra_Map(-1, data->Dnc, DColElems, 0,
445 DRowMap =
new Epetra_Map(-1, data->Dnr, data->DRowElems, 0, LComm);
446 SRowMap =
new Epetra_Map(-1, data->Snr, data->SRowElems, 0, LComm);
447 DColMap =
new Epetra_Map(-1, data->Dnc, DColElems, 0, LComm);
456 C->FillComplete(*SRowMap, *DRowMap);
460 R->FillComplete(*DColMap, *SRowMap);
463 int Sdiag = (int) data->Snr * Sdiagfactor;
464 Sdiag = MIN(Sdiag, data->Snr-1);
465 Sdiag = MAX(Sdiag, 0);
468 for (
int i = 0; config->schurApproxMethod == 1 && i < nrows ; i++)
471 if (gvals[gid] == 1)
continue;
472 if (data->Snr == 0) assert(0 == 1);
478 for (
int j = MAX(i-Sdiag, 0) ; j < MIN(data->Snr, i+Sdiag); j++)
482 RightIndex[rcnt++] = data->SRowElems[j];
484 err = Sg->InsertGlobalIndices(gid, rcnt, RightIndex);
487 err = Sg->InsertGlobalIndices(gid, 1, &gid);
491 if (config->schurApproxMethod == 1)
503 delete[] RightValues;
536 int myPID = A->Comm().MyPID();
537 int n = A->NumGlobalRows();
543 int sym = config->sym;
548 Epetra_Map AColMap = A->ColMap();
549 int ncols = AColMap.NumMyElements();
550 int *cols = AColMap.MyGlobalElements();
553 Epetra_Map ARowMap = A->RowMap();
554 int nrows = ARowMap.NumMyElements();
555 int *rows = ARowMap.MyGlobalElements();
558 int *gvals =
new int[n];
561 findLocalColumns(A, gvals, SNumGlobalCols);
568 if (config->sep_type == 2)
569 findNarrowSeparator(A, gvals);
579 ostringstream ssmsg1;
580 ssmsg1 <<
"PID =" << myPID <<
" ";
581 string msg = ssmsg1.str();
582 ssmsg1.clear(); ssmsg1.str(
"");
604 for (
int i = 0; i < nrows ; i++)
606 if (gvals[rows[i]] == 1)
613 for (
int i = 0; i < ncols ; i++)
615 if (gvals[cols[i]] != 1)
624 cout << msg <<
" Mycols="<< ncols <<
"Myrows ="<< nrows << endl;
625 cout << msg <<
" #rows and #cols in diagonal blk ="<< Dnr << endl;
626 cout << msg <<
" #columns in S ="<< Snc << endl;
627 cout << msg <<
" #rows in S ="<< Snr << endl;
629 ostringstream pidstr;
632 DRowElems =
new int[Dnr];
633 SRowElems =
new int[Snr];
638 findBlockElems(A, nrows, rows, gvals, Dnr, DRowElems, Snr, SRowElems,
639 "D"+pidstr.str()+
"Rows",
"S"+pidstr.str()+
"Rows",
false) ;
650 data->DRowElems = DRowElems;
651 data->SRowElems = SRowElems;
654 int *DColElems =
new int[Dnc];
655 int *SColElems =
new int[Snc];
657 findBlockElems(A, ncols, cols, gvals, Dnc, DColElems, Snc, SColElems,
658 "D"+pidstr.str()+
"Cols",
"S"+pidstr.str()+
"Cols",
true) ;
660 data->DColElems = DColElems;
663 for (
int i = 0; i < Snr; i++)
667 assert(SRowElems[i] == SColElems[i]);
672 create_matrices(A, ssym, data, config);
675 extract_matrices(A, ssym, data, config,
true);
680 (
new Epetra_LinearProblem());
681 LP->SetOperator((ssym->D).getRawPtr());
686 (
new Epetra_MultiVector(ssym->D->RowMap(), 16));
688 (
new Epetra_MultiVector(ssym->D->RowMap(), 16));
690 (
new Epetra_MultiVector(ssym->G->RowMap(), 16));
694 (
new Epetra_MultiVector(
View, *(ssym->Drhs), 0, 1));
696 (
new Epetra_MultiVector(
View, *(ssym->Dlhs), 0, 1));
702 EpetraExt::ViewTransform<Epetra_LinearProblem> >
703 (
new EpetraExt::LinearProblem_Reindex2(0));
709 std::size_t found = (config->diagonalBlockSolver).find(
"Amesos");
713 const char* SolverType = config->diagonalBlockSolver.c_str();
718 (void) Factory.Query(SolverType);
720 bool IsAvailable = Factory.Query(SolverType);
721 assert(IsAvailable ==
true);
723 config->amesosForDiagonal =
true;
725 (Factory.Create(SolverType, *(ssym->LP)));
729 config->amesosForDiagonal =
false;
730 #ifdef HAVE_IFPACK_DYNAMIC_FACTORY 731 Ifpack_DynamicFactory IfpackFactory;
733 Ifpack IfpackFactory;
735 ifpackSolver = Teuchos::rcp<Ifpack_Preconditioner>(IfpackFactory.Create
736 (config->diagonalBlockSolver, (ssym->D).getRawPtr(), 0,
false));
745 if (config->amesosForDiagonal)
748 aList.
set(
"TrustMe",
true);
749 Solver->SetParameters(aList);
750 Solver->SymbolicFactorization();
754 ifpackSolver->Initialize();
767 ssym->Solver = Solver;
768 ssym->ifSolver = ifpackSolver;
770 if (config->schurApproxMethod == 1)
775 Isorropia::Epetra::Prober((ssym->Sg).getRawPtr(),
788 ssym->prober = prober;
792 #ifdef HAVE_SHYLUCORE_MPI 793 Epetra_MpiComm LComm(MPI_COMM_SELF);
795 Epetra_SerialComm LComm;
796 #endif // HAVE_SHYLUCORE_MPI 797 data->LDRowMap =
Teuchos::rcp(
new Epetra_Map(-1, data->Dnr,
798 data->DRowElems, 0, LComm));
799 data->LGRowMap =
Teuchos::rcp(
new Epetra_Map(-1, data->Snr,
800 data->SRowElems, 0, LComm));
802 data->SRowElems, 0, A->Comm()));
805 data->BdImporter =
Teuchos::rcp(
new Epetra_Import(*(data->LDRowMap),
807 data->DistImporter =
Teuchos::rcp(
new Epetra_Import(*(data->GMap),
808 *(data->LGRowMap)));;
810 data->BsImporter =
Teuchos::rcp(
new Epetra_Import(*(data->GMap),
812 data->XsImporter =
Teuchos::rcp(
new Epetra_Import(*(data->LGRowMap),
815 data->XdExporter =
Teuchos::rcp(
new Epetra_Export(*(data->LDRowMap),
817 data->XsExporter =
Teuchos::rcp(
new Epetra_Export(*(data->LGRowMap),
821 data->localrhs =
Teuchos::rcp(
new Epetra_MultiVector(*(data->LDRowMap), 1));
822 data->temp1 =
Teuchos::rcp(
new Epetra_MultiVector(*(data->LGRowMap), 1));
823 data->temp2 =
Teuchos::rcp(
new Epetra_MultiVector(*(data->GMap), 1));
824 data->Bs =
Teuchos::rcp(
new Epetra_MultiVector(*(data->GMap), 1));
825 data->Xs =
Teuchos::rcp(
new Epetra_MultiVector(*(data->GMap), 1));
826 data->LocalXs =
Teuchos::rcp(
new Epetra_MultiVector(*(data->LGRowMap), 1));
827 data->Xs->PutScalar(0.0);
857 extract_matrices(A, ssym, data, config,
false);
866 if (config->amesosForDiagonal)
867 Solver->NumericFactorization();
869 ifpackSolver->Compute();
874 if (config->amesosForDiagonal)
875 Solver->PrintStatus();
885 #ifdef HAVE_SHYLUCORE_MPI 886 Epetra_MpiComm LComm(MPI_COMM_SELF);
888 Epetra_SerialComm LComm;
889 #endif // HAVE_SHYLUCORE_MPI 890 Epetra_Map LocalDRowMap(-1, data->Dnr, data->DRowElems, 0, LComm);
895 (ssym->R).getRawPtr(), (ssym->LP).getRawPtr(),
896 (ssym->Solver).getRawPtr(), (ssym->ifSolver).getRawPtr(),
897 (ssym->C).getRawPtr(), &LocalDRowMap, 1));
899 if (config->schurApproxMethod == 1)
901 int nvectors = prober->getNumOrthogonalVectors();
906 (ssym->R).getRawPtr(), (ssym->LP).getRawPtr(),
907 (ssym->Solver).getRawPtr(), (ssym->ifSolver).getRawPtr(),
908 (ssym->C).getRawPtr(), &LocalDRowMap, nvectors);
911 Sbar = prober->probe(probeop);
919 else if (config->schurApproxMethod == 2)
929 if (config->sep_type == 2)
931 (ssym->R).getRawPtr(), (ssym->LP).getRawPtr(),
932 (ssym->Solver).getRawPtr(), (ssym->ifSolver).getRawPtr(),
933 (ssym->C).getRawPtr(), &LocalDRowMap);
936 (ssym->R).getRawPtr(), (ssym->LP).getRawPtr(),
937 (ssym->Solver).getRawPtr(), (ssym->ifSolver).getRawPtr(),
938 (ssym->C).getRawPtr(),
944 cout <<
"Time to Compute Approx Schur Complement" << ftime.
totalElapsedTime() << endl;
949 else if (config->schurApproxMethod == 3)
966 if (config->schurSolver ==
"Amesos")
972 (
new Epetra_MultiVector(data->Sbar->RowMap(), 16));
974 (
new Epetra_MultiVector(data->Sbar->RowMap(), 16));
979 LP2->SetOperator(data->Sbar.
get());
984 <EpetraExt::ViewTransform<Epetra_LinearProblem> >
985 (
new EpetraExt::LinearProblem_Reindex2(0));
988 (&((*(data->ReIdx_LP2))(*LP2)),
false);
998 (void) Factory.Query(config->schurAmesosSolver);
1000 bool IsAvailable = Factory.Query(config->schurAmesosSolver);
1001 assert(IsAvailable ==
true);
1003 data->dsolver = Factory.Create(
1004 config->schurAmesosSolver, *(data->LP2));
1006 data->dsolver->SetParameters(aList);
1009 data->dsolver->SymbolicFactorization();
1013 data->dsolver->NumericFactorization();
1017 else if (config->schurSolver ==
"AztecOO-Exact")
1019 if (config->libName ==
"AztecOO") {
1020 data->innersolver =
new AztecOO();
1022 std::string schurPrec = config->schurPreconditioner;
1024 int err = data->innersolver->SetUserOperator
1025 (data->schur_op.
get());
1028 #ifdef HAVE_IFPACK_DYNAMIC_FACTORY 1029 Ifpack_DynamicFactory IfpackFactory;
1031 Ifpack IfpackFactory;
1033 data->schur_prec = Teuchos::rcp<Ifpack_Preconditioner>
1034 (IfpackFactory.Create(schurPrec,
1035 Sbar.
getRawPtr(), config->overlap,
false));
1037 data->schur_prec->Initialize();
1038 data->schur_prec->Compute();
1040 err = data->innersolver->SetPrecOperator
1041 (data->schur_prec.
get());
1044 data->innersolver->SetAztecOption(AZ_solver, AZ_gmres);
1046 if (config->silent_subiter) {
1047 data->innersolver->SetAztecOption(AZ_output, AZ_none);
1050 data->innersolver->SetMatrixName(999);
1053 else if (config->schurSolver ==
"AztecOO-Inexact")
1056 if (config->libName ==
"Belos")
1058 int err = data->innersolver->SetUserMatrix
1062 data->innersolver->SetAztecOption(AZ_solver, AZ_gmres);
1063 data->innersolver->SetAztecOption(AZ_precond, AZ_dom_decomp);
1064 data->innersolver->SetAztecOption(AZ_keep_info, 1);
1066 data->innersolver->SetAztecOption(AZ_output, AZ_none);
1070 data->innersolver->SetMatrixName(999);
1073 err = data->innersolver->ConstructPreconditioner(condest);
1085 data->innersolver = NULL;
1088 else if ((config->schurSolver ==
"IQR") || (config->schurSolver ==
"G"))
1097 #ifdef TIMING_OUTPUT
Teuchos::RCP< Epetra_CrsMatrix > computeSchur_GuidedProbing(shylu_config *config, shylu_symbolic *ssym, shylu_data *data, Epetra_Map *localDRowMap)
Compute an approximate Schur Complement using the option of Guided Probing.
Teuchos::RCP< Epetra_CrsMatrix > computeApproxSchur(shylu_config *config, shylu_symbolic *ssym, Epetra_CrsMatrix *G, Epetra_CrsMatrix *R, Epetra_LinearProblem *LP, Amesos_BaseSolver *solver, Ifpack_Preconditioner *ifSolver, Epetra_CrsMatrix *C, Epetra_Map *localDRowMap)
Compute an approximate Schur Complement (Narrow Sep)
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
Teuchos::RCP< Epetra_CrsMatrix > computeApproxWideSchur(shylu_config *config, shylu_symbolic *ssym, Epetra_CrsMatrix *G, Epetra_CrsMatrix *R, Epetra_LinearProblem *LP, Amesos_BaseSolver *solver, Ifpack_Preconditioner *ifSolver, Epetra_CrsMatrix *C, Epetra_Map *localDRowMap)
Compute an approximate Shur Complete (Wide Sep)
int shylu_factor(Epetra_CrsMatrix *A, shylu_symbolic *ssym, shylu_data *data, shylu_config *config)
Main function call into ShylU.
void start(bool reset=false)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
double totalElapsedTime(bool readCurrentTime=false) const
int shylu_symbolic_factor(Epetra_CrsMatrix *A, shylu_symbolic *ssym, shylu_data *data, shylu_config *config)
Call symbolic factorization on matrix.
Main data structure holding needed offset and temp variables.
Main header file of ShyLU (Include main user calls)