ShyLU  Version of the Day
shylu_factor.cpp
Go to the documentation of this file.
1 
2 //@HEADER
3 // ************************************************************************
4 //
5 // ShyLU: Hybrid preconditioner package
6 // Copyright 2012 Sandia Corporation
7 //
8 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9 // the U.S. Government retains certain rights in this software.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
39 //
40 // ************************************************************************
41 //@HEADER
42 
59 /* This define will make S block diagonal, To make this happen even some
60  zero columns/rows of C and R will be stored.
61  */
62 #define BLOCK_DIAGONAL_Si
63 
64 #include "shylu.h"
65 #include "shylu_util.h"
66 #ifdef HAVE_SHYLUCORE_MPI
67 #include <Epetra_MpiComm.h>
68 #else
69 #include <Epetra_SerialComm.h>
70 #endif
71 #include <EpetraExt_Reindex_LinearProblem2.h>
72 
73 #include "Ifpack_config.h"
74 
75 #ifdef HAVE_IFPACK_DYNAMIC_FACTORY
76 #include "Ifpack_DynamicFactory.h"
77 #else
78 #include "Ifpack.h"
79 #endif
80 
81 #include "shylu_internal_gmres.h"
82 #include "shylu_internal_gmres_tools.h"
83 
84 int create_matrices
85 (
86  Epetra_CrsMatrix *A, // i/p: A matrix
87  shylu_symbolic *ssym, // symbolic structure
88  shylu_data *data, // numeric structure, TODO: Required ?
89  shylu_config *config // i/p: library configuration
90 )
91 {
92  int Dnr = data->Dnr;
93  int Snr = data->Snr;
94  int Dnc = data->Dnc;
95  int *DRowElems = data->DRowElems;
96  int *SRowElems = data->SRowElems;
97  int *DColElems = data->DColElems;
98  int *gvals = data->gvals;
99 
100  double Sdiagfactor = config->Sdiagfactor;
101  int sym = config->sym;
102 
103  /* --------------- Create the maps for the DBBD form ------------------- */
104  // Create the local and distributed row map
105 #ifdef HAVE_SHYLUCORE_MPI
106  Epetra_MpiComm LComm(MPI_COMM_SELF);
107 #else
108  Epetra_SerialComm LComm;
109 #endif // HAVE_SHYLUCORE_MPI
110 
111  // Use Serial Comm for the local blocks.
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);
116 
117  // Create the local column map
118  Epetra_Map LocalDColMap(-1, Dnc, DColElems, 0, LComm);
119 
120  /*----------- Compute nnz per row for all five matrices ----------------- */
121  // Compute nnz per row.
122  int *Ai;
123  double *Ax;
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; // Max nnz in any row
130 
131  // Find the required no of diagonals
132  /*int Sdiag = (int) SNumGlobalCols * Sdiagfactor;
133  //cout << "No of diagonals in Sbar =" << Sdiag << endl;
134  Sdiag = MIN(Sdiag, SNumGlobalCols-1);*/
135  int Sdiag = (int) Snr * Sdiagfactor;
136  Sdiag = MIN(Sdiag, Snr-1);
137  Sdiag = MAX(Sdiag, 0);
138  //cout << "No of diagonals in Sbar =" << Sdiag << endl;
139  //assert (Sdiag <= SNumGlobalCols-1);
140  if (Snr != 0) assert (Sdiag <= Snr-1);
141 
142  int dcol, ccol, rcol, scol;
143  int gcid;
144  int gid;
145 
146  int nrows = A->RowMap().NumMyElements();
147  int *rows = A->RowMap().MyGlobalElements();
148  if (sym)
149  {
150 
151  int NumEntries;
152  int dcnt, scnt, ccnt, rcnt;
153  dcol = 0 ; ccol = 0; rcol = 0; scol = 0;
154  // Three things to worry about here, Num entries in each row of S,
155  // size of C and R and the corresponding cols/rows.
156  for (int i = 0; i < nrows; i++)
157  {
158  dcnt = 0; scnt = 0; ccnt = 0; rcnt = 0;
159  // Need to pass local id i to this function
160  int err = A->ExtractMyRowView(i, NumEntries, Ax, Ai);
161  if (err != 0)
162  {
163  assert (err == 0);
164  //config->dm.error("create_matrices: extract_my_row failed");
165  }
166 
167  gid = rows[i];
168  if (gvals[gid] == 1)
169  { // D or C row
170  for (int j = 0 ; j < NumEntries ; j++)
171  { // O(nnz) ! Careful what you do inside
172  gcid = A->GCID(Ai[j]);
173  // Only cols corresponding to local rows belong to D.
174  if (A->LRID(gcid) != -1 && gvals[gcid] == 1)
175  {
176  dcnt++;
177  }
178  else
179  {
180  ccnt++;
181  }
182  }
183  // There should not be a null row in D.
184  assert(dcnt != 0);
185  DNumEntriesPerRow[dcol++] = dcnt;
186  CNumEntriesPerRow[ccol++] = ccnt;
187  }
188  else
189  { // R or S row
190  //cout << "proc/row " << myPID << "/" << gid;
191  //cout << " has Cols = " ;
192  for (int j = 0 ; j < NumEntries ; j++)
193  { // O(nnz) ! Careful what you do inside
194  gcid = A->GCID(Ai[j]);
195  //if (A->LRID(gcid) != -1 && gvals[gcid] == 1) // TBD : Need to change here
196  // No need to check for local rows as above.
197  // All cols with gvals[cols] == 1 belong to R.
198  if (gvals[gcid] == 1)
199  {
200  rcnt++;
201  }
202  else
203  {
204  scnt++;
205  }
206  }
207  // There should not be a null row in S.
208  assert(scnt != 0);
209  assert(scol < Snr);
210  assert(rcol < Snr);
211  GNumEntriesPerRow[scol++] = scnt;
212  RNumEntriesPerRow[rcol++] = rcnt;
213  }
214  Dmaxnnz = max(Dmaxnnz, dcnt);
215  Smaxnnz = max(Smaxnnz, scnt);
216  Rmaxnnz = max(Rmaxnnz, rcnt);
217  Cmaxnnz = max(Cmaxnnz, ccnt);
218  }
219  assert( dcol == Dnr);
220  assert( scol == Snr);
221  assert( ccol == Dnr);
222  assert( rcol == Snr);
223  int sbarcol = 0;
224  for (int i = 0; i < nrows; i++)
225  {
226  gid = rows[i];
227  if (gvals[gid] != 1)
228  { // S row
229  // Add the number of required diagonals in approximate Schur
230  // complement , +1 below for the main diagonal
231  assert(sbarcol < Snr);
232  SBarNumEntriesPerRow[sbarcol] = GNumEntriesPerRow[sbarcol]
233  + Sdiag*2 + 1;
234  sbarcol++;
235  //SBarNumEntriesPerRow[i] = GNumEntriesPerRow[i] ;
236 
237  }
238  }
239  assert( sbarcol == Snr);
240  Smaxnnz = Smaxnnz + Sdiag * 2 + 1;
241  }
242  else
243  {
244  assert (0 == 1);
245  /*assert(Dnr == nrows);
246  // TODO : Need to compute the entries for D and C in one array of size
247  // nrow*/
248  }
249 
250 
251  //Create the local matrices
252  ssym->D = Teuchos::RCP<Epetra_CrsMatrix> (new Epetra_CrsMatrix(Copy,
253  LocalDRowMap, LocalDColMap, DNumEntriesPerRow, true));
254  //config->dm.print(5, "Created D matrix");
255 
256  // Leave the column map out, Let Epetra do the work in the rest of the cases
257  ssym->G = Teuchos::RCP<Epetra_CrsMatrix> (new Epetra_CrsMatrix(Copy,
258  SRowMap, GNumEntriesPerRow, true));
259  //config->dm.print(5, "Created S matrix");
260  if (config->sep_type != 1)
261  {
262 
263  ssym->C = Teuchos::RCP<Epetra_CrsMatrix> (new Epetra_CrsMatrix(Copy,
264  DRowMap, CNumEntriesPerRow, true));
265  //config->dm.print(5, "Created C matrix");
266 
267  ssym->R = Teuchos::RCP<Epetra_CrsMatrix> (new Epetra_CrsMatrix(Copy,
268  SRowMap, RNumEntriesPerRow, true));
269  //config->dm.print(5, "Created R matrix");
270  }
271  else
272  {
273  ssym->C = Teuchos::RCP<Epetra_CrsMatrix> (new Epetra_CrsMatrix(Copy,
274  LocalDRowMap, CNumEntriesPerRow, true));
275  //config->dm.print(5, "Created C matrix");
276 
277  ssym->R = Teuchos::RCP<Epetra_CrsMatrix> (new Epetra_CrsMatrix(Copy,
278  LocalSRowMap, RNumEntriesPerRow, true));
279  //config->dm.print(5, "Created R matrix");
280  }
281 
282  if (config->schurApproxMethod == 1)
283  {
284  // TODO: Check for the local case.
285  ssym->Sg = Teuchos::RCP<Epetra_CrsGraph> (new Epetra_CrsGraph(Copy,
286  SRowMap, SBarNumEntriesPerRow, false));
287  //config->dm.print(5, "Created Sg graph");
288  }
289 
290  data->lmax = max(Dmaxnnz, Rmaxnnz);
291  data->rmax = max(Cmaxnnz, Smaxnnz);
292 
293  delete[] DNumEntriesPerRow;
294  delete[] GNumEntriesPerRow;
295  delete[] CNumEntriesPerRow;
296  delete[] RNumEntriesPerRow;
297  delete[] SBarNumEntriesPerRow;
298  return 0;
299 }
300 
301 int extract_matrices
302 (
303  Epetra_CrsMatrix *A, // i/p: A matrix
304  shylu_symbolic *ssym, // symbolic structure
305  shylu_data *data, // numeric structure, TODO: Required ?
306  shylu_config *config, // i/p: library configuration
307  bool insertValues // true implies values will be inserted and fill
308  // complete will be called. false implies values
309  // will be replaced.
310 )
311 {
312  Teuchos::RCP<Epetra_CrsMatrix> D = ssym->D;
313  Teuchos::RCP<Epetra_CrsMatrix> C = ssym->C;
314  Teuchos::RCP<Epetra_CrsMatrix> R = ssym->R;
315  Teuchos::RCP<Epetra_CrsMatrix> G = ssym->G;
316  Teuchos::RCP<Epetra_CrsGraph> Sg = ssym->Sg;
317  int *DColElems = data->DColElems;
318  int *gvals = data->gvals;
319  double Sdiagfactor = config->Sdiagfactor;
320 
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];
325  int err;
326  int lcnt, rcnt ;
327  int gcid;
328  int gid;
329  int *Ai;
330  double *Ax;
331 
332  int nrows = A->RowMap().NumMyElements();
333  int *rows = A->RowMap().MyGlobalElements();
334 
335  for (int i = 0; i < nrows ; i++)
336  {
337  int NumEntries;
338  err = A->ExtractMyRowView(i, NumEntries, Ax, Ai);
339 
340  lcnt = 0; rcnt = 0;
341  // Place the entry in the correct sub matrix, Works only for sym
342  gid = rows[i];
343  int lcid;
344  for (int j = 0 ; j < NumEntries ; j++)
345  { // O(nnz) ! Careful what you do inside
346  // Row permutation does not matter here
347  gcid = A->GCID(Ai[j]);
348  assert(gcid != -1);
349  //Either in D or R
350  if ((gvals[gid] != 1 && gvals[gcid] == 1)
351  || (gvals[gid] == 1 && A->LRID(gcid) != -1 && gvals[gcid] == 1))
352  {
353  assert(lcnt < data->lmax);
354  if (insertValues)
355  LeftIndex[lcnt] = gcid;
356  else
357  {
358  //local column id
359  lcid = (gvals[gid] == 1 ? D->LCID(gcid) : R->LCID(gcid));
360  assert(lcid != -1);
361  LeftIndex[lcnt] = lcid;
362  }
363  LeftValues[lcnt++] = Ax[j];
364  }
365  else
366  {
367  assert(rcnt < data->rmax);
368  if (insertValues)
369  RightIndex[rcnt] = gcid;
370  else
371  {
372  //local column id
373  lcid = (gvals[gid] == 1 ? C->LCID(gcid) : G->LCID(gcid));
374  assert(lcid != -1);
375  RightIndex[rcnt] = lcid;
376  }
377  RightValues[rcnt++] = Ax[j];
378  }
379  }
380 
381  if (gvals[gid] == 1)
382  { // D or C row
383  if (insertValues)
384  {
385  err = D->InsertGlobalValues(gid, lcnt, LeftValues, LeftIndex);
386  assert(err == 0);
387  err = C->InsertGlobalValues(gid, rcnt, RightValues, RightIndex);
388  assert(err == 0);
389  }
390  else
391  {
392  err = D->ReplaceMyValues(D->LRID(gid), lcnt, LeftValues,
393  LeftIndex);
394  assert(err == 0);
395  err = C->ReplaceMyValues(C->LRID(gid), rcnt, RightValues,
396  RightIndex);
397  assert(err == 0);
398  }
399  }
400  else
401  { // R or S row
402  //assert(lcnt > 0); // TODO: Enable this once using narrow sep.
403  if (insertValues)
404  {
405  assert(rcnt > 0);
406  err = R->InsertGlobalValues(gid, lcnt, LeftValues, LeftIndex);
407  assert(err == 0);
408  err = G->InsertGlobalValues(gid, rcnt, RightValues, RightIndex);
409  assert(err == 0);
410  if (config->schurApproxMethod == 1)
411  {
412  err = Sg->InsertGlobalIndices(gid, rcnt, RightIndex);
413  assert(err == 0);
414  }
415  }
416  else
417  {
418  assert(rcnt > 0);
419  err = R->ReplaceMyValues(R->LRID(gid), lcnt, LeftValues,
420  LeftIndex);
421  assert(err == 0);
422  err = G->ReplaceMyValues(G->LRID(gid), rcnt, RightValues,
423  RightIndex);
424  assert(err == 0);
425  }
426  }
427  }
428 
429  if (insertValues)
430  {
431  /* ------------- Create the maps for the DBBD form ------------------ */
432  Epetra_Map *DRowMap, *SRowMap, *DColMap;
433  Epetra_SerialComm LComm;
434  if (config->sep_type != 1)
435  {
436  DRowMap = new Epetra_Map(-1, data->Dnr, data->DRowElems, 0,
437  A->Comm());
438  SRowMap = new Epetra_Map(-1, data->Snr, data->SRowElems, 0,
439  A->Comm());
440  DColMap = new Epetra_Map(-1, data->Dnc, DColElems, 0,
441  A->Comm());
442  }
443  else
444  {
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);
448  }
449 
450  D->FillComplete();
451  //config->dm.print(5, "Done D fillcomplete");
452 
453  G->FillComplete();
454  //config->dm.print(5, "Done G fillcomplete");
455 
456  C->FillComplete(*SRowMap, *DRowMap); //TODO:Won't work if permutation is
457  // unsymmetric SRowMap
458  //config->dm.print(5, "Done C fillcomplete");
459 
460  R->FillComplete(*DColMap, *SRowMap);
461  //config->dm.print(5, "Done R fillcomplete");
462 
463  int Sdiag = (int) data->Snr * Sdiagfactor;
464  Sdiag = MIN(Sdiag, data->Snr-1);
465  Sdiag = MAX(Sdiag, 0);
466 
467  // Add the diagonals to Sg
468  for (int i = 0; config->schurApproxMethod == 1 && i < nrows ; i++)
469  {
470  gid = rows[i];
471  if (gvals[gid] == 1) continue; // not a row in S
472  if (data->Snr == 0) assert(0 == 1);
473 
474  rcnt = 0;
475  //TODO Will be trouble if SNumGlobalCols != Snc
476  //assert(SNumGlobalCols == Snc);
477  //for (int j = MAX(i-Sdiag,0) ; j<MIN(SNumGlobalCols, i+Sdiag); j++)
478  for (int j = MAX(i-Sdiag, 0) ; j < MIN(data->Snr, i+Sdiag); j++)
479  {
480  // find the adjacent columns from the row map of S
481  //assert (j >= 0 && j < Snr);
482  RightIndex[rcnt++] = data->SRowElems[j];
483  }
484  err = Sg->InsertGlobalIndices(gid, rcnt, RightIndex);
485  assert(err == 0);
486  // Always insert the diagonals, if it is added twice that is fine.
487  err = Sg->InsertGlobalIndices(gid, 1, &gid);
488  assert(err == 0);
489  }
490 
491  if (config->schurApproxMethod == 1)
492  Sg->FillComplete();
493 
494  delete DRowMap;
495  delete SRowMap;
496  delete DColMap;
497  }
498 
499  // A is no longer needed
500  delete[] LeftIndex;
501  delete[] LeftValues;
502  delete[] RightIndex;
503  delete[] RightValues;
504 
505  //cout << msg << "S rows=" << S.NumGlobalRows() << " S cols=" <<
506  //S.NumGlobalCols() << "#cols in column map="<<
507  //S.ColMap().NumMyElements() << endl;
508  //cout << msg << "C rows=" << Cptr->NumGlobalRows() << " C cols=" <<
509  //Cptr->NumGlobalCols() << "#cols in column map="<<
510  //Cptr->ColMap().NumMyElements() << endl;
511  //cout << msg << "D rows=" << D.NumGlobalRows() << " D cols=" <<
512  //D.NumGlobalCols() << "#cols in column map="<<
513  //D.ColMap().NumMyElements() << endl;
514  //cout << msg << "R rows=" << Rptr->NumGlobalRows() << " R cols=" <<
515  //Rptr->NumGlobalCols() << "#cols in column map="<<
516  //Rptr->ColMap().NumMyElements() << endl;
517  // ]
518 
519  return 0;
520 
521 }
522 
523 /* Find the DBBD form */
525 (
526  Epetra_CrsMatrix *A, // i/p: A matrix
527  shylu_symbolic *ssym, // symbolic structure
528  shylu_data *data, // numeric structure, TODO: Required ?
529  shylu_config *config // i/p: library configuration
530 )
531 {
532 #ifdef TIMING_OUTPUT
533  Teuchos::Time symtime("symbolic time");
534  symtime.start();
535 #endif
536  int myPID = A->Comm().MyPID();
537  int n = A->NumGlobalRows();
538 
539  int Dnr;
540  int Snr;
541  int *DRowElems;
542  int *SRowElems;
543  int sym = config->sym;
544 
545  checkMaps(A);
546 
547  // Get column map
548  Epetra_Map AColMap = A->ColMap();
549  int ncols = AColMap.NumMyElements();
550  int *cols = AColMap.MyGlobalElements();
551 
552  // Get row map
553  Epetra_Map ARowMap = A->RowMap();
554  int nrows = ARowMap.NumMyElements();
555  int *rows = ARowMap.MyGlobalElements();
556 
557  // Find all columns in this proc
558  int *gvals = new int[n]; // vector of size n, not ncols !
559  // gvals[local cols] = 1, gvals[shared cols] > 1.
560  int SNumGlobalCols;
561  findLocalColumns(A, gvals, SNumGlobalCols);
562 
563  // See if you can shrink the separator by assigning more rows/columns to
564  // the block diagonals
565  // TODO: This is because of a bug in coloring remove the if once that is
566  // fixed
567  //if (config->schurApproxMethod == 2)
568  if (config->sep_type == 2)
569  findNarrowSeparator(A, gvals);
570 
571  // 3. Assemble diagonal block and the border in convenient form [
572  /* In each processor, we have (in a permuted form)
573  * | D_i C_i |
574  * | R_i S_i |
575  * D_i - diagonal block, C_i - Column Separator, R_i - Row separator
576  * S_i - A22 block corresponding to Schur complement part of A
577  * Assemble all four blocks in local matrices. */
578 
579  ostringstream ssmsg1;
580  ssmsg1 << "PID =" << myPID << " ";
581  string msg = ssmsg1.str();
582  ssmsg1.clear(); ssmsg1.str("");
583 
584  // Find #cols in each block
585  int Dnc = 0; // #cols in diagonal block
586  int Snc = 0; // #cols in the col. separator
587  /* Looping on cols will work only for wide separator
588  * as for narrow sep there will be some sep cols with gvals[col] ==1
589  * */
590  /*for (int i = 0; i < ncols ; i++)
591  {
592  if (gvals[cols[i]] == 1)
593  Dnc++;
594  else
595  Snc++;
596  }
597  // Find #rows in each block
598  Dnr = Dnc; // #rows in square diagonal block
599  Snr = nrows - Dnr; // #rows in the row separator*/
600 
601  // Find #rows in each block
602  Dnr = 0;
603  Snr = 0;
604  for (int i = 0; i < nrows ; i++)
605  {
606  if (gvals[rows[i]] == 1)
607  Dnr++;
608  else
609  Snr++;
610  }
611  Dnc = Dnr;
612  // TODO: Snc is no longer useful, should remove it
613  for (int i = 0; i < ncols ; i++)
614  {
615  if (gvals[cols[i]] != 1)
616  Snc++;
617  }
618 
619  assert(Snc >= 0);
620 
621  // TODO : The above assignment may not be correct in the unsymetric case
622 
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;
628 
629  ostringstream pidstr;
630  pidstr << myPID ;
631  // Create a row map for the D and S blocks [
632  DRowElems = new int[Dnr];
633  SRowElems = new int[Snr];
634  // int gid; // unused
635  // Assemble row ids in two arrays (for D and R blocks)
636  if (sym)
637  {
638  findBlockElems(A, nrows, rows, gvals, Dnr, DRowElems, Snr, SRowElems,
639  "D"+pidstr.str()+"Rows", "S"+pidstr.str()+"Rows", false) ;
640  }
641  else
642  {
643  // SRowElems are not known until factorization, TODO
644  assert(0 == 1);
645  }
646 
647  data->Dnr = Dnr;
648  data->Snr = Snr;
649  data->Dnc = Dnc;
650  data->DRowElems = DRowElems;
651  data->SRowElems = SRowElems;
652 
653  // Create a column map for the D and S blocks [
654  int *DColElems = new int[Dnc]; // Elems in column map of D
655  int *SColElems = new int[Snc]; // Elems in column map of C TODO: Unused
656  // Assemble column ids in two arrays (for D and C blocks)
657  findBlockElems(A, ncols, cols, gvals, Dnc, DColElems, Snc, SColElems,
658  "D"+pidstr.str()+"Cols", "S"+pidstr.str()+"Cols", true) ;
659 
660  data->DColElems = DColElems;
661  data->gvals = gvals;
662 
663  for (int i = 0; i < Snr; i++)
664  {
665  // Epetra guarentees columns corresponding to local rows will be first
666  // in the column map.
667  assert(SRowElems[i] == SColElems[i]);
668  }
669  // ]
670 
671  /*--Create the Epetra Matrices with the maps (does not insert values) --- */
672  create_matrices(A, ssym, data, config);
673 
674  /*--Extract the Epetra Matrices and call fillComplete --- */
675  extract_matrices(A, ssym, data, config, true);
676 
677  delete[] SColElems;
678 
679  Teuchos::RCP<Epetra_LinearProblem> LP = Teuchos::RCP<Epetra_LinearProblem>
680  (new Epetra_LinearProblem());
681  LP->SetOperator((ssym->D).getRawPtr());
682  //LP->SetOperator((ssym->DT).getRawPtr()); // for transpose
683 
684  // Create temp vectors
685  ssym->Dlhs = Teuchos::RCP<Epetra_MultiVector>
686  (new Epetra_MultiVector(ssym->D->RowMap(), 16));
687  ssym->Drhs = Teuchos::RCP<Epetra_MultiVector>
688  (new Epetra_MultiVector(ssym->D->RowMap(), 16));
689  ssym->Gvec = Teuchos::RCP<Epetra_MultiVector>
690  (new Epetra_MultiVector(ssym->G->RowMap(), 16));
691 
692  // Are these even needed, plan to remove them ?
693  data->temp3 = Teuchos::RCP<Epetra_MultiVector>
694  (new Epetra_MultiVector(View, *(ssym->Drhs), 0, 1));
695  data->locallhs = Teuchos::RCP<Epetra_MultiVector>
696  (new Epetra_MultiVector(View, *(ssym->Dlhs), 0, 1));
697 
698  LP->SetRHS(ssym->Drhs.getRawPtr());
699  LP->SetLHS(ssym->Dlhs.getRawPtr());
700 
701  ssym->ReIdx_LP = Teuchos::RCP<
702  EpetraExt::ViewTransform<Epetra_LinearProblem> >
703  (new EpetraExt::LinearProblem_Reindex2(0));
704  ssym->LP = Teuchos::RCP<Epetra_LinearProblem>(&((*(ssym->ReIdx_LP))(*LP)),
705  false);
706 
707  Teuchos::RCP<Amesos_BaseSolver> Solver;
708  Teuchos::RCP<Ifpack_Preconditioner> ifpackSolver;
709  std::size_t found = (config->diagonalBlockSolver).find("Amesos");
710  if (found == 0)
711  {
712  Amesos Factory;
713  const char* SolverType = config->diagonalBlockSolver.c_str();
714  // mfh 25 May 2015: assert() gets defined to nothing in a
715  // release build. Thus, the return value IsAvailable won't
716  // get used, which causes a build warning.
717 #ifdef NDEBUG
718  (void) Factory.Query(SolverType);
719 #else
720  bool IsAvailable = Factory.Query(SolverType);
721  assert(IsAvailable == true);
722 #endif // NDEBUG
723  config->amesosForDiagonal = true;
724  Solver = Teuchos::RCP<Amesos_BaseSolver>
725  (Factory.Create(SolverType, *(ssym->LP)));
726  }
727  else
728  {
729  config->amesosForDiagonal = false;
730 #ifdef HAVE_IFPACK_DYNAMIC_FACTORY
731  Ifpack_DynamicFactory IfpackFactory;
732 #else
733  Ifpack IfpackFactory;
734 #endif
735  ifpackSolver = Teuchos::rcp<Ifpack_Preconditioner>(IfpackFactory.Create
736  (config->diagonalBlockSolver, (ssym->D).getRawPtr(), 0, false));
737  }
738  //config->dm.print(5, "Created the diagonal solver");
739 
740 #ifdef TIMING_OUTPUT
741  Teuchos::Time ftime("setup time");
742  ftime.start();
743 #endif
744  //Solver->SetUseTranspose(true); // for transpose
745  if (config->amesosForDiagonal)
746  {
747  Teuchos::ParameterList aList;
748  aList.set("TrustMe", true);
749  Solver->SetParameters(aList);
750  Solver->SymbolicFactorization();
751  }
752  else
753  {
754  ifpackSolver->Initialize();
755  }
756 
757  //config->dm.print(3, "Symbolic Factorization done");
758 
759 #ifdef TIMING_OUTPUT
760  ftime.stop();
761  cout << "Symbolic Factorization Time" << ftime.totalElapsedTime() << endl;
762  ftime.reset();
763 #endif
764 
765  ssym->OrigLP = LP;
766  //ssym->LP = LP;
767  ssym->Solver = Solver;
768  ssym->ifSolver = ifpackSolver;
769 
770  if (config->schurApproxMethod == 1)
771  {
772  Teuchos::ParameterList pList;
773  Teuchos::RCP<Isorropia::Epetra::Prober> prober =
774  Teuchos::RCP<Isorropia::Epetra::Prober> (new
775  Isorropia::Epetra::Prober((ssym->Sg).getRawPtr(),
776  pList, false));
777  //config->dm.print(3, "Doing Coloring");
778 #ifdef TIMING_OUTPUT
779  ftime.start();
780 #endif
781  prober->color();
782 #ifdef TIMING_OUTPUT
783  ftime.stop();
784  cout << "Time to color" << ftime.totalElapsedTime() << endl;
785  ftime.reset();
786  ftime.start();
787 #endif
788  ssym->prober = prober;
789  }
790 
791  // Set the maps, importers and multivectors to be used in the solve once.
792 #ifdef HAVE_SHYLUCORE_MPI
793  Epetra_MpiComm LComm(MPI_COMM_SELF);
794 #else
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));
801  data->GMap = Teuchos::rcp(new Epetra_Map(-1, data->Snr,
802  data->SRowElems, 0, A->Comm()));
803 
804  // Assuming X and A will have the same rowmap. Should it be domain map ?
805  data->BdImporter = Teuchos::rcp(new Epetra_Import(*(data->LDRowMap),
806  A->RowMap()));
807  data->DistImporter = Teuchos::rcp(new Epetra_Import(*(data->GMap),
808  *(data->LGRowMap)));;
809  // Assuming X and A will have the same rowmap. Should it be domain map ?
810  data->BsImporter = Teuchos::rcp(new Epetra_Import(*(data->GMap),
811  A->RowMap()));
812  data->XsImporter = Teuchos::rcp(new Epetra_Import(*(data->LGRowMap),
813  *(data->GMap)));
814  // Assuming Y and A will have the same rowmap. Should it be range map ?
815  data->XdExporter = Teuchos::rcp(new Epetra_Export(*(data->LDRowMap),
816  A->RowMap()));
817  data->XsExporter = Teuchos::rcp(new Epetra_Export(*(data->LGRowMap),
818  A->RowMap()));
819 
820  // Create multivectors for solve, TODO: Can we do with fewer
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);
828 
829  //data->importExportTime = Teuchos::rcp(new Teuchos::Time("import export time"));
830  //data->innerIterTime = Teuchos::rcp(new Teuchos::Time("innertIter time"));
831  //data->fwdTime = Teuchos::rcp(new Teuchos::Time("reindex fwd time"));
832  //data->amesosSchurTime = Teuchos::rcp(new Teuchos::Time("amesos schur time"));
833 #ifdef TIMING_OUTPUT
834  symtime.stop();
835  cout << "Symbolic Time" << symtime.totalElapsedTime() << endl;
836  symtime.reset();
837 #endif
838 
839  // FIXME (mfh 25 May 2015) Uh, shouldn't this function return an int???
840 }
841 
842 int shylu_factor(Epetra_CrsMatrix *A, shylu_symbolic *ssym, shylu_data *data,
843  shylu_config *config)
844 {
845 #ifdef TIMING_OUTPUT
846  Teuchos::Time fact_time("factor time");
847  fact_time.start();
848 #endif
849 
850  Teuchos::RCP<Epetra_LinearProblem> LP = ssym->LP;
851  Teuchos::RCP<Amesos_BaseSolver> Solver = ssym->Solver;
852  Teuchos::RCP<Ifpack_Preconditioner> ifpackSolver = ssym->ifSolver;
853  Teuchos::RCP<Isorropia::Epetra::Prober> prober = ssym->prober;
854 
855  /*--Extract the Epetra Matrices into already existing matrices --- */
856  extract_matrices(A, ssym, data, config, false);
857 
858  // ======================= Numeric factorization =========================
859 #ifdef TIMING_OUTPUT
860  Teuchos::Time ftime("setup time");
861  ftime.start();
862 #endif
863 
864  //config->dm.print(3, "In Numeric Factorization");
865  if (config->amesosForDiagonal)
866  Solver->NumericFactorization();
867  else
868  ifpackSolver->Compute();
869 
870  //config->dm.print(3, "Numeric Factorization done");
871 
872 #ifdef SHYLU_DEBUG
873  if (config->amesosForDiagonal)
874  Solver->PrintStatus();
875 #endif
876 
877 #ifdef TIMING_OUTPUT
878  ftime.stop();
879  cout << "Time to factor" << ftime.totalElapsedTime() << endl;
880  ftime.reset();
881 #endif
882 
883  // Create the local and distributed row map
884 #ifdef HAVE_SHYLUCORE_MPI
885  Epetra_MpiComm LComm(MPI_COMM_SELF);
886 #else
887  Epetra_SerialComm LComm;
888 #endif // HAVE_SHYLUCORE_MPI
889  Epetra_Map LocalDRowMap(-1, data->Dnr, data->DRowElems, 0, LComm);
890  Teuchos::RCP<Epetra_CrsMatrix> Sbar;
891 
892  data->schur_op = Teuchos::RCP<ShyLU_Probing_Operator> (new
893  ShyLU_Probing_Operator(config, ssym, (ssym->G).getRawPtr(),
894  (ssym->R).getRawPtr(), (ssym->LP).getRawPtr(),
895  (ssym->Solver).getRawPtr(), (ssym->ifSolver).getRawPtr(),
896  (ssym->C).getRawPtr(), &LocalDRowMap, 1));
897 
898  if (config->schurApproxMethod == 1)
899  {
900  int nvectors = prober->getNumOrthogonalVectors();
901  //Set up the probing operator
902  // TODO: Change to RCPs. Call Set vectors on schur_op and remove
903  // probeop
904  ShyLU_Probing_Operator probeop(config, ssym, (ssym->G).getRawPtr(),
905  (ssym->R).getRawPtr(), (ssym->LP).getRawPtr(),
906  (ssym->Solver).getRawPtr(), (ssym->ifSolver).getRawPtr(),
907  (ssym->C).getRawPtr(), &LocalDRowMap, nvectors);
908 
909  //cout << "Doing probing" << endl;
910  Sbar = prober->probe(probeop);
911  //cout << "SIZE of SBAR = " << (*Sbar).NumGlobalRows() << endl;
912 #ifdef TIMING_OUTPUT
913  ftime.stop();
914  cout << "Time to probe" << ftime.totalElapsedTime() << endl;
915  ftime.reset();
916 #endif
917  }
918  else if (config->schurApproxMethod == 2)
919  {
920  // Compute and drop the entries in the Schur complement
921  // Ignore the structure of the Schur complement
922  //cout << "Computing the Approx Schur complement" << endl;
923 #ifdef TIMING_OUTPUT
924  ftime.start();
925 #endif
926 
927  // TODO: Pass the operator rather than all the matrices
928  if (config->sep_type == 2)
929  Sbar = computeApproxSchur(config, ssym, (ssym->G).getRawPtr(),
930  (ssym->R).getRawPtr(), (ssym->LP).getRawPtr(),
931  (ssym->Solver).getRawPtr(), (ssym->ifSolver).getRawPtr(),
932  (ssym->C).getRawPtr(), &LocalDRowMap);
933  else
934  Sbar = computeApproxWideSchur(config, ssym, (ssym->G).getRawPtr(),
935  (ssym->R).getRawPtr(), (ssym->LP).getRawPtr(),
936  (ssym->Solver).getRawPtr(), (ssym->ifSolver).getRawPtr(),
937  (ssym->C).getRawPtr(),
938  &LocalDRowMap);
939 
940  //cout << *Sbar << endl;
941 #ifdef TIMING_OUTPUT
942  ftime.stop();
943  cout << "Time to Compute Approx Schur Complement" << ftime.totalElapsedTime() << endl;
944  ftime.reset();
945 #endif
946  //cout << "Computed Approx Schur complement" << endl;
947  }
948  else if (config->schurApproxMethod == 3)
949  {
950  Sbar = computeSchur_GuidedProbing(config, ssym, data, &LocalDRowMap);
951 
952  // TODO: The above call does not work with narrow separator. We need
953  // distribute probing here, save the Sbar from the iteration in
954  // data in a different variable other than Sbar (as Sbar gets
955  // overwritten after every iteration, but we need the original Sbar
956  // graph for the probing in later iterations) Set the prober, when
957  // new dropped is computed in previous iteration, otherwise just do
958  // probing, same idea from above function this time it is distributed.
959  // This was slower when I implemented it for wide sep (and I deleted
960  // it !).
961  }
962 
963  data->Sbar = Sbar;
964 
965  if (config->schurSolver == "Amesos")
966  {
967  Teuchos::ParameterList aList;
968  //aList.set("Reindex", true);
969 
970  data->Sbarlhs = Teuchos::RCP<Epetra_MultiVector>
971  (new Epetra_MultiVector(data->Sbar->RowMap(), 16));
972  data->Sbarrhs = Teuchos::RCP<Epetra_MultiVector>
973  (new Epetra_MultiVector(data->Sbar->RowMap(), 16));
974 
975  Teuchos::RCP<Epetra_LinearProblem> LP2 =
976  Teuchos::rcp(new Epetra_LinearProblem);
977 
978  LP2->SetOperator(data->Sbar.get());
979  LP2->SetRHS(data->Sbarrhs.getRawPtr());
980  LP2->SetLHS(data->Sbarlhs.getRawPtr());
981 
982  data->ReIdx_LP2 = Teuchos::RCP
983  <EpetraExt::ViewTransform<Epetra_LinearProblem> >
984  (new EpetraExt::LinearProblem_Reindex2(0));
985  data->LP2 =
986  Teuchos::RCP<Epetra_LinearProblem>
987  (&((*(data->ReIdx_LP2))(*LP2)), false);
988 
989  data->OrigLP2 = LP2;
990 
991  Amesos Factory;
992  // mfh 25 May 2015: assert() gets defined to nothing in a
993  // release build. This causes a build warning that
994  // IsAvailable is set but not used. That's why I protect its
995  // declaration with the NDEBUG macro.
996 #ifdef NDEBUG
997  (void) Factory.Query(config->schurAmesosSolver);
998 #else
999  bool IsAvailable = Factory.Query(config->schurAmesosSolver);
1000  assert(IsAvailable == true);
1001 #endif // NDEBUG
1002  data->dsolver = Factory.Create(
1003  config->schurAmesosSolver, *(data->LP2));
1004 
1005  data->dsolver->SetParameters(aList);
1006  //cout << "Created the direct Schur Solver" << endl;
1007 
1008  data->dsolver->SymbolicFactorization();
1009  //cout << "Symbolic Factorization done for schur complement" << endl;
1010 
1011  //cout << "In Numeric Factorization of Schur complement" << endl;
1012  data->dsolver->NumericFactorization();
1013  //cout << "Numeric Factorization done for schur complement" << endl;
1014 
1015  }
1016  else if (config->schurSolver == "AztecOO-Exact")
1017  {
1018  if (config->libName == "AztecOO") {
1019  data->innersolver = new AztecOO();
1020  }
1021  std::string schurPrec = config->schurPreconditioner;
1022 
1023  int err = data->innersolver->SetUserOperator
1024  (data->schur_op.get());
1025  assert (err == 0);
1026 
1027 #ifdef HAVE_IFPACK_DYNAMIC_FACTORY
1028  Ifpack_DynamicFactory IfpackFactory;
1029 #else
1030  Ifpack IfpackFactory;
1031 #endif
1032  data->schur_prec = Teuchos::rcp<Ifpack_Preconditioner>
1033  (IfpackFactory.Create(schurPrec,
1034  Sbar.getRawPtr(), config->overlap, false));
1035 
1036  data->schur_prec->Initialize();
1037  data->schur_prec->Compute();
1038 
1039  err = data->innersolver->SetPrecOperator
1040  (data->schur_prec.get());
1041  assert (err == 0);
1042 
1043  data->innersolver->SetAztecOption(AZ_solver, AZ_gmres);
1044 
1045  if (config->silent_subiter) {
1046  data->innersolver->SetAztecOption(AZ_output, AZ_none);
1047  }
1048 
1049  data->innersolver->SetMatrixName(999);
1050 
1051  }
1052  else if (config->schurSolver == "AztecOO-Inexact")
1053  {
1054  // Doing an inexact Schur complement
1055  if (config->libName == "Belos")
1056  {
1057  int err = data->innersolver->SetUserMatrix
1058  (data->Sbar.get());
1059  assert (err == 0);
1060 
1061  data->innersolver->SetAztecOption(AZ_solver, AZ_gmres);
1062  data->innersolver->SetAztecOption(AZ_precond, AZ_dom_decomp);
1063  data->innersolver->SetAztecOption(AZ_keep_info, 1);
1064 #ifndef DEBUG
1065  data->innersolver->SetAztecOption(AZ_output, AZ_none);
1066 #endif
1067  //solver->SetAztecOption(AZ_overlap, 3);
1068  //solver->SetAztecOption(AZ_subdomain_solve, AZ_ilu);
1069  data->innersolver->SetMatrixName(999);
1070 
1071  double condest;
1072  err = data->innersolver->ConstructPreconditioner(condest);
1073  assert (err == 0);
1074  //cout << "Condition number of inner Sbar" << condest << endl;
1075  }
1076  else
1077  {
1078  // The outer solver is Aztec
1079  // I suspect there is a bug in AztecOO. Doing what we do in the if
1080  // case
1081  // here will cause an error when we use the solver in ApplyInverse
1082  // The error will not happen when we call the dummy JustTryIt()
1083  // below
1084  data->innersolver = NULL;
1085  }
1086  }
1087  else if ((config->schurSolver == "IQR") || (config->schurSolver == "G"))
1088  {
1089  // DO NOTHING
1090  } else
1091  {
1092  assert (0 == 1);
1093  }
1094 
1095  //cout << " Out of factor" << endl ;
1096 #ifdef TIMING_OUTPUT
1097  fact_time.stop();
1098  cout << "Factor Time" << fact_time.totalElapsedTime() << endl;
1099  fact_time.reset();
1100 #endif
1101  return 0;
1102 }
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)
Definition: shylu_schur.cpp:56
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.
int shylu_symbolic_factor(Epetra_CrsMatrix *A, shylu_symbolic *ssym, shylu_data *data, shylu_config *config)
Call symbolic factorization on matrix.
Utilities for ShyLU.
Main data structure holding needed offset and temp variables.
Definition: shylu.h:108
Main header file of ShyLU (Include main user calls)