47 #include "Epetra_Vector.h" 48 #include "shylu_probing_operator.h" 49 #include "shylu_local_schur_operator.h" 50 #include <EpetraExt_Reindex_CrsMatrix.h> 58 Epetra_CrsMatrix *G, Epetra_CrsMatrix *R,
59 Epetra_LinearProblem *LP, Amesos_BaseSolver *solver,
60 Ifpack_Preconditioner *ifSolver, Epetra_CrsMatrix *C,
61 Epetra_Map *localDRowMap)
63 double relative_thres = config->relative_threshold;
67 localDRowMap, nvectors);
70 Epetra_Map rMap = G->RowMap();
71 int *rows = rMap.MyGlobalElements();
72 int totalElems = rMap.NumGlobalElements();
73 int localElems = rMap.NumMyElements();
83 G->Comm().ScanSum(&localElems, &prefixSum, 1);
86 int *mySGID =
new int[totalElems];
87 int *allSGID =
new int[totalElems];
89 for (i = 0, j = 0; i < totalElems ; i++)
91 if (i >= prefixSum - localElems && i < prefixSum)
103 C->Comm().SumAll(mySGID, allSGID, totalElems);
115 Copy, rMap, localElems));
117 double *values =
new double[localElems];
118 int *indices =
new int[localElems];
121 double *maxvalue =
new double[nvectors];
128 int findex = totalElems / nvectors ;
129 for (i = 0 ; i < findex*nvectors ; i+=nvectors)
131 Epetra_MultiVector probevec(rMap, nvectors);
132 Epetra_MultiVector Scol(rMap, nvectors);
134 probevec.PutScalar(0.0);
136 for (
int k = 0; k < nvectors; k++)
139 if (cindex >= prefixSum - localElems && cindex < prefixSum)
141 probevec.ReplaceGlobalValue(allSGID[cindex], k, 1.0);
148 probeop.Apply(probevec, Scol);
152 Scol.MaxValue(maxvalue);
153 for (
int k = 0; k < nvectors; k++)
158 for (j = 0 ; j < localElems ; j++)
161 if (allSGID[cindex] == rows[j])
163 values[nentries] = vecvalues[j];
164 indices[nentries] = allSGID[cindex];
166 Sbar->InsertGlobalValues(rows[j], nentries, values, indices);
168 else if (abs(vecvalues[j]/maxvalue[k]) > relative_thres)
170 values[nentries] = vecvalues[j];
171 indices[nentries] = allSGID[cindex];
173 Sbar->InsertGlobalValues(rows[j], nentries, values, indices);
177 if (vecvalues[j] != 0.0) dropped++;
183 probeop.ResetTempVectors(1);
185 for ( ; i < totalElems ; i++)
187 Epetra_MultiVector probevec(rMap, 1);
188 Epetra_MultiVector Scol(rMap, 1);
190 probevec.PutScalar(0.0);
191 if (i >= prefixSum - localElems && i < prefixSum)
193 probevec.ReplaceGlobalValue(allSGID[i], 0, 1.0);
199 probeop.Apply(probevec, Scol);
204 Scol.MaxValue(maxvalue);
206 for (j = 0 ; j < localElems ; j++)
209 if (allSGID[i] == rows[j])
211 values[nentries] = vecvalues[j];
212 indices[nentries] = allSGID[i];
214 Sbar->InsertGlobalValues(rows[j], nentries, values, indices);
216 else if (abs(vecvalues[j]/maxvalue[0]) > relative_thres)
218 values[nentries] = vecvalues[j];
219 indices[nentries] = allSGID[i];
221 Sbar->InsertGlobalValues(rows[j], nentries, values, indices);
225 if (vecvalues[j] != 0.0) dropped++;
231 cout <<
"Time in finding and dropping entries" << ftime.
totalElapsedTime() << endl;
237 Sbar->FillComplete();
238 cout <<
"#dropped entries" << dropped << endl;
251 Epetra_CrsMatrix *G, Epetra_CrsMatrix *R,
252 Epetra_LinearProblem *LP, Amesos_BaseSolver *solver,
253 Ifpack_Preconditioner *ifSolver, Epetra_CrsMatrix *C,
254 Epetra_Map *localDRowMap)
257 double relative_thres = config->relative_threshold;
269 Epetra_Map GrMap = G->RowMap();
270 int *g_rows = GrMap.MyGlobalElements();
272 int g_localElems = GrMap.NumMyElements();
281 Epetra_SerialComm LComm;
282 Epetra_Map G_localRMap (-1, g_localElems, g_rows, 0, LComm);
287 int maxentries = max(C->MaxNumEntries(), R->MaxNumEntries());
288 maxentries = max(maxentries, G->MaxNumEntries());
290 double *values1 =
new double[maxentries];
291 double *values2 =
new double[maxentries];
292 double *values3 =
new double[maxentries];
293 int *indices1 =
new int[maxentries];
294 int *indices2 =
new int[maxentries];
295 int *indices3 =
new int[maxentries];
299 Copy, GrMap, g_localElems));
302 Epetra_CrsMatrix localG(
Copy, G_localRMap, G->MaxNumEntries(),
false);
304 for (i = 0; i < g_localElems ; i++)
307 G->ExtractGlobalRowCopy(gid, maxentries, nentries1, values1, indices1);
311 for (
int j = 0 ; j < nentries1 ; j++)
313 if (G->LRID(indices1[j]) != -1)
315 values2[cnt] = values1[j];
316 indices2[cnt++] = indices1[j];
321 values3[scnt] = values1[j];
322 indices3[scnt++] = indices1[j];
326 localG.InsertGlobalValues(gid, cnt, values2, indices2);
327 Sbar->InsertGlobalValues(gid, scnt, values3, indices3);
329 localG.FillComplete();
336 ifSolver, C, localDRowMap, nvectors);
368 double *values =
new double[nvectors];
369 int *indices =
new int[nvectors];
375 #endif // SHYLU_DEBUG 376 double *maxvalue =
new double[nvectors];
380 int findex = g_localElems / nvectors ;
384 Epetra_MultiVector probevec (G_localRMap, nvectors);
385 Epetra_MultiVector Scol (G_localRMap, nvectors);
386 probevec.PutScalar(0.0);
387 for (i = 0 ; i < findex*nvectors ; i+=nvectors)
390 for (
int k = 0; k < nvectors; k++)
395 probevec.ReplaceGlobalValue(g_rows[cindex], k, 1.0);
403 probeop.Apply(probevec, Scol);
409 for (
int k = 0; k < nvectors; k++)
412 probevec.ReplaceGlobalValue(g_rows[cindex], k, 0.0);
415 Scol.MaxValue(maxvalue);
417 for (
int j = 0 ; j < g_localElems ; j++)
419 for (
int k = 0; k < nvectors; k++)
423 if ((g_rows[cindex] == g_rows[j]) ||
424 (abs(vecvalues[j]/maxvalue[k]) > relative_thres))
427 values[nentries] = vecvalues[j];
428 indices[nentries++] = g_rows[cindex];
431 else if (vecvalues[j] != 0.0)
435 #endif // SHYLU_DEBUG 437 Sbar->InsertGlobalValues(g_rows[j], nentries, values,
443 if (i < g_localElems)
445 nvectors = g_localElems - i;
446 probeop.ResetTempVectors(nvectors);
447 Epetra_MultiVector probevec1 (G_localRMap, nvectors);
448 Epetra_MultiVector Scol1 (G_localRMap, nvectors);
450 probevec1.PutScalar(0.0);
451 for (
int k = 0; k < nvectors; k++)
456 probevec1.ReplaceGlobalValue(g_rows[cindex], k, 1.0);
462 probeop.Apply(probevec1, Scol1);
466 Scol1.MaxValue(maxvalue);
468 for (
int j = 0 ; j < g_localElems ; j++)
471 for (
int k = 0; k < nvectors; k++)
474 vecvalues = Scol1[k];
476 if ((g_rows[cindex] == g_rows[j]) ||
477 (abs(vecvalues[j]/maxvalue[k]) > relative_thres))
480 values[nentries] = vecvalues[j];
481 indices[nentries++] = g_rows[cindex];
484 else if (vecvalues[j] != 0.0)
488 #endif // SHYLU_DEBUG 490 Sbar->InsertGlobalValues(g_rows[j], nentries, values,
502 probeop.PrintTimingInfo();
504 Sbar->FillComplete();
507 Epetra_Map defMap2(-1, g_localElems, 0, C->Comm());
508 EpetraExt::ViewTransform<Epetra_CrsMatrix> * ReIdx_MatTrans2 =
509 new EpetraExt::CrsMatrix_Reindex( defMap2 );
510 Epetra_CrsMatrix t2S = (*ReIdx_MatTrans2)( *Sbar );
511 ReIdx_MatTrans2->fwd();
512 EpetraExt::RowMatrixToMatlabFile(
"Schur.mat", t2S);
516 cout <<
"#dropped entries" << dropped << endl;
538 Epetra_Map *localDRowMap
542 double relative_thres = config->relative_threshold;
544 Epetra_CrsMatrix *G = ssym->G.
getRawPtr();
545 Epetra_CrsMatrix *R = ssym->R.
getRawPtr();
546 Epetra_LinearProblem *LP = ssym->LP.
getRawPtr();
547 Amesos_BaseSolver *solver = ssym->Solver.
getRawPtr();
548 Ifpack_Preconditioner *ifSolver = ssym->ifSolver.
getRawPtr();
549 Epetra_CrsMatrix *C = ssym->C.
getRawPtr();
554 Epetra_Map CrMap = C->RowMap();
555 int *c_rows = CrMap.MyGlobalElements();
556 int *c_cols = (C->ColMap()).MyGlobalElements();
558 int c_localElems = CrMap.NumMyElements();
559 int c_localcolElems = (C->ColMap()).NumMyElements();
561 Epetra_Map GrMap = G->RowMap();
562 int *g_rows = GrMap.MyGlobalElements();
564 int g_localElems = GrMap.NumMyElements();
566 Epetra_Map RrMap = R->RowMap();
567 int *r_rows = RrMap.MyGlobalElements();
568 int *r_cols = (R->ColMap()).MyGlobalElements();
570 int r_localElems = RrMap.NumMyElements();
571 int r_localcolElems = (R->ColMap()).NumMyElements();
573 Epetra_SerialComm LComm;
574 Epetra_Map C_localRMap (-1, c_localElems, c_rows, 0, LComm);
575 Epetra_Map C_localCMap (-1, c_localcolElems, c_cols, 0, LComm);
576 Epetra_Map G_localRMap (-1, g_localElems, g_rows, 0, LComm);
577 Epetra_Map R_localRMap (-1, r_localElems, r_rows, 0, LComm);
578 Epetra_Map R_localCMap (-1, r_localcolElems, r_cols, 0, LComm);
583 cout <<
"DEBUG MODE" << endl;
584 int nrows = C->RowMap().NumMyElements();
585 assert(nrows == localDRowMap->NumGlobalElements());
587 int gids[nrows], gids1[nrows];
588 C_localRMap.MyGlobalElements(gids);
589 localDRowMap->MyGlobalElements(gids1);
590 cout <<
"Comparing R's domain map with D's row map" << endl;
592 for (
int i = 0; i < nrows; i++)
594 assert(gids[i] == gids1[i]);
601 int maxentries = max(C->MaxNumEntries(), R->MaxNumEntries());
602 maxentries = max(maxentries, G->MaxNumEntries());
604 double *values1 =
new double[maxentries];
605 double *values2 =
new double[maxentries];
606 double *values3 =
new double[maxentries];
607 int *indices1 =
new int[maxentries];
608 int *indices2 =
new int[maxentries];
609 int *indices3 =
new int[maxentries];
613 Epetra_CrsMatrix localC(
Copy, C_localRMap, C->MaxNumEntries(),
false);
614 for (i = 0; i < c_localElems ; i++)
617 err = C->ExtractGlobalRowCopy(gid, maxentries, nentries1, values1,
622 err = localC.InsertGlobalValues(gid, nentries1, values1, indices1);
626 localC.FillComplete(G_localRMap, C_localRMap);
630 Epetra_CrsMatrix localR(
Copy, R_localRMap, R->MaxNumEntries(),
false);
631 for (i = 0; i < r_localElems ; i++)
634 R->ExtractGlobalRowCopy(gid, maxentries, nentries1, values1, indices1);
635 localR.InsertGlobalValues(gid, nentries1, values1, indices1);
637 localR.FillComplete(*localDRowMap, R_localRMap);
642 Copy, GrMap, g_localElems));
645 Epetra_CrsMatrix localG(
Copy, G_localRMap, G->MaxNumEntries(),
false);
647 for (i = 0; i < g_localElems ; i++)
650 G->ExtractGlobalRowCopy(gid, maxentries, nentries1, values1, indices1);
654 for (
int j = 0 ; j < nentries1 ; j++)
656 if (G->LRID(indices1[j]) != -1)
658 values2[cnt] = values1[j];
659 indices2[cnt++] = indices1[j];
664 values3[scnt] = values1[j];
665 indices3[scnt++] = indices1[j];
669 localG.InsertGlobalValues(gid, cnt, values2, indices2);
670 Sbar->InsertGlobalValues(gid, scnt, values3, indices3);
672 localG.FillComplete();
673 cout <<
"Created local G matrix" << endl;
677 ifSolver, &localC, localDRowMap, nvectors);
705 new Epetra_CrsGraph(
Copy, G_localRMap, maxentries));
707 if (data->num_compute % config->reset_iter == 0)
713 double *values =
new double[g_localElems];
714 int *indices =
new int[g_localElems];
717 double *maxvalue =
new double[nvectors];
721 int findex = g_localElems / nvectors ;
725 Epetra_MultiVector probevec(G_localRMap, nvectors);
726 Epetra_MultiVector Scol(G_localRMap, nvectors);
727 for (i = 0 ; i < findex*nvectors ; i+=nvectors)
729 probevec.PutScalar(0.0);
730 for (
int k = 0; k < nvectors; k++)
736 probevec.ReplaceGlobalValue(g_rows[cindex], k, 1.0);
744 probeop.Apply(probevec, Scol);
749 Scol.MaxValue(maxvalue);
750 for (
int k = 0; k < nvectors; k++)
755 for (
int j = 0 ; j < g_localElems ; j++)
758 if (g_rows[cindex] == g_rows[j])
760 values[nentries] = vecvalues[j];
761 indices[nentries] = g_rows[cindex];
763 err = Sbar->InsertGlobalValues(g_rows[j], nentries, values,
766 err = lSGraph->InsertGlobalIndices(g_rows[j], nentries,
770 else if (abs(vecvalues[j]/maxvalue[k]) > relative_thres)
772 values[nentries] = vecvalues[j];
773 indices[nentries] = g_rows[cindex];
775 err = Sbar->InsertGlobalValues(g_rows[j], nentries, values,
778 err = lSGraph->InsertGlobalIndices(g_rows[j], nentries,
784 if (vecvalues[j] != 0.0)
795 probeop.ResetTempVectors(1);
797 for ( ; i < g_localElems ; i++)
800 Epetra_MultiVector probevec(G_localRMap, 1);
801 Epetra_MultiVector Scol(G_localRMap, 1);
803 probevec.PutScalar(0.0);
806 probevec.ReplaceGlobalValue(g_rows[i], 0, 1.0);
811 probeop.Apply(probevec, Scol);
816 Scol.MaxValue(maxvalue);
818 for (
int j = 0 ; j < g_localElems ; j++)
821 if (g_rows[i] == g_rows[j])
823 values[nentries] = vecvalues[j];
824 indices[nentries] = g_rows[i];
826 err = Sbar->InsertGlobalValues(g_rows[j], nentries, values, indices);
828 err = lSGraph->InsertGlobalIndices(g_rows[j], nentries, indices);
831 else if (abs(vecvalues[j]/maxvalue[0]) > relative_thres)
833 values[nentries] = vecvalues[j];
834 indices[nentries] = g_rows[i];
836 err = Sbar->InsertGlobalValues(g_rows[j], nentries, values, indices);
838 err = lSGraph->InsertGlobalIndices(g_rows[j], nentries, indices);
843 if (vecvalues[j] != 0.0) dropped++;
850 cout <<
"Time in finding and dropping entries" << ftime.
totalElapsedTime() << endl;
856 probeop.PrintTimingInfo();
857 Sbar->FillComplete();
858 lSGraph->FillComplete();
860 data->localSbargraph = lSGraph;
863 Epetra_Map defMap2(-1, g_localElems, 0, C->Comm());
864 EpetraExt::ViewTransform<Epetra_CrsMatrix> * ReIdx_MatTrans2 =
865 new EpetraExt::CrsMatrix_Reindex( defMap2 );
866 Epetra_CrsMatrix t2S = (*ReIdx_MatTrans2)( *Sbar );
867 ReIdx_MatTrans2->fwd();
868 EpetraExt::RowMatrixToMatlabFile(
"Schur.mat", t2S);
871 cout <<
"#dropped entries" << dropped << endl;
878 if (((data->num_compute-1) % config->reset_iter) == 0)
886 Isorropia::Epetra::Prober(
887 data->localSbargraph.
getRawPtr(), pList,
false));
889 data->guided_prober = gprober;
894 int nvectors = data->guided_prober->getNumOrthogonalVectors();
895 cout <<
"Number of Orthogonal Vectors for guided probing" << nvectors
898 probeop.ResetTempVectors(nvectors);
900 data->guided_prober->probe(probeop);
901 int maxentries = blockdiag_Sbar->GlobalMaxNumEntries();
902 int *indices =
new int[maxentries];
903 double *values =
new double[maxentries];
906 for (
int i = 0; i < blockdiag_Sbar->NumGlobalRows() ; i++)
908 int gid = blockdiag_Sbar->GRID(i);
909 blockdiag_Sbar->ExtractGlobalRowCopy(gid, maxentries, numentries,
911 Sbar->InsertGlobalValues(gid, numentries, values, indices);
914 Sbar->FillComplete();
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)
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)
void start(bool reset=false)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
double totalElapsedTime(bool readCurrentTime=false) const
Main data structure holding needed offset and temp variables.
Main header file of ShyLU (Include main user calls)