IFPACK  Development
 All Classes Namespaces Files Functions Variables Enumerations Friends Pages
Ifpack_CrsRiluk.cpp
1 //@HEADER
2 // ***********************************************************************
3 //
4 // Ifpack: Object-Oriented Algebraic Preconditioner Package
5 // Copyright (2002) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38 //
39 // ***********************************************************************
40 //@HEADER
41 
42 #include "Ifpack_CrsRiluk.h"
43 #include "Epetra_ConfigDefs.h"
44 #include "Epetra_Comm.h"
45 #include "Epetra_Map.h"
46 #include "Epetra_CrsGraph.h"
47 #include "Epetra_VbrMatrix.h"
48 #include "Epetra_RowMatrix.h"
49 #include "Epetra_MultiVector.h"
50 
51 #include <Teuchos_ParameterList.hpp>
52 #include <ifp_parameters.h>
53 
54 //==============================================================================
56  : UserMatrixIsVbr_(false),
57  UserMatrixIsCrs_(false),
58  Graph_(Graph_in),
59  Comm_(Graph_in.Comm()),
60  UseTranspose_(false),
61  NumMyDiagonals_(0),
62  Allocated_(false),
63  ValuesInitialized_(false),
64  Factored_(false),
65  RelaxValue_(0.0),
66  Athresh_(0.0),
67  Rthresh_(1.0),
68  Condest_(-1.0),
69  OverlapMode_(Zero)
70 {
71  // Test for non-trivial overlap here so we can use it later.
72  IsOverlapped_ = (Graph_in.LevelOverlap()>0 && Graph_in.DomainMap().DistributedGlobal());
73 }
74 
75 //==============================================================================
77  : UserMatrixIsVbr_(FactoredMatrix.UserMatrixIsVbr_),
78  UserMatrixIsCrs_(FactoredMatrix.UserMatrixIsCrs_),
79  IsOverlapped_(FactoredMatrix.IsOverlapped_),
80  Graph_(FactoredMatrix.Graph_),
81  IlukRowMap_(FactoredMatrix.IlukRowMap_),
82  IlukDomainMap_(FactoredMatrix.IlukDomainMap_),
83  IlukRangeMap_(FactoredMatrix.IlukRangeMap_),
84  Comm_(FactoredMatrix.Comm_),
85  UseTranspose_(FactoredMatrix.UseTranspose_),
86  NumMyDiagonals_(FactoredMatrix.NumMyDiagonals_),
87  Allocated_(FactoredMatrix.Allocated_),
88  ValuesInitialized_(FactoredMatrix.ValuesInitialized_),
89  Factored_(FactoredMatrix.Factored_),
90  RelaxValue_(FactoredMatrix.RelaxValue_),
91  Athresh_(FactoredMatrix.Athresh_),
92  Rthresh_(FactoredMatrix.Rthresh_),
93  Condest_(FactoredMatrix.Condest_),
94  OverlapMode_(FactoredMatrix.OverlapMode_)
95 {
96  L_ = Teuchos::rcp( new Epetra_CrsMatrix(FactoredMatrix.L()) );
97  U_ = Teuchos::rcp( new Epetra_CrsMatrix(FactoredMatrix.U()) );
98  D_ = Teuchos::rcp( new Epetra_Vector(FactoredMatrix.D()) );
99  if (IlukRowMap_!=Teuchos::null) IlukRowMap_ = Teuchos::rcp( new Epetra_Map(*IlukRowMap_) );
100  if (IlukDomainMap_!=Teuchos::null) IlukDomainMap_ = Teuchos::rcp( new Epetra_Map(*IlukDomainMap_) );
101  if (IlukRangeMap_!=Teuchos::null) IlukRangeMap_ = Teuchos::rcp( new Epetra_Map(*IlukRangeMap_) );
102 
103 }
104 //==============================================================================
106 
107  ValuesInitialized_ = false;
108  Factored_ = false;
109  Allocated_ = false;
110 }
111 //==============================================================================
112 int Ifpack_CrsRiluk::AllocateCrs() {
113 
114  // Allocate Epetra_CrsMatrix using ILUK graphs
115  L_ = Teuchos::rcp( new Epetra_CrsMatrix(Copy, Graph_.L_Graph()) );
116  U_ = Teuchos::rcp( new Epetra_CrsMatrix(Copy, Graph_.U_Graph()) );
117  D_ = Teuchos::rcp( new Epetra_Vector(Graph_.L_Graph().RowMap()) );
118  L_Graph_ = Teuchos::null;
119  U_Graph_ = Teuchos::null;
120  SetAllocated(true);
121  return(0);
122 }
123 //==============================================================================
124 int Ifpack_CrsRiluk::AllocateVbr() {
125 
126  // First we need to create a set of Epetra_Maps that has the same number of points as the
127  // Epetra_BlockMaps associated with the Overlap Graph.
128  EPETRA_CHK_ERR(BlockMap2PointMap(Graph_.L_Graph().RowMap(), &IlukRowMap_));
129  EPETRA_CHK_ERR(BlockMap2PointMap(Graph_.U_Graph().DomainMap(), &IlukDomainMap_));
130  EPETRA_CHK_ERR(BlockMap2PointMap(Graph_.L_Graph().RangeMap(), &IlukRangeMap_));
131 
132  // Set L range map and U domain map
133  U_DomainMap_ = IlukDomainMap_;
134  L_RangeMap_ = IlukRangeMap_;
135  // If there is fill, then pre-build the L and U structures from the Block version of L and U.
136  if (Graph().LevelFill()) {
137  L_Graph_ = Teuchos::rcp( new Epetra_CrsGraph(Copy, *IlukRowMap_, *IlukRowMap_, 0) );
138  U_Graph_ = Teuchos::rcp( new Epetra_CrsGraph(Copy, *IlukRowMap_, *IlukRowMap_, 0) );
139  EPETRA_CHK_ERR(BlockGraph2PointGraph(Graph_.L_Graph(), *L_Graph_, false));
140  EPETRA_CHK_ERR(BlockGraph2PointGraph(Graph_.U_Graph(), *U_Graph_, true));
141 
142  L_Graph_->FillComplete(*IlukRowMap_, *IlukRangeMap_);
143  U_Graph_->FillComplete(*IlukDomainMap_, *IlukRowMap_);
144 
145  L_ = Teuchos::rcp( new Epetra_CrsMatrix(Copy, *L_Graph_) );
146  U_ = Teuchos::rcp( new Epetra_CrsMatrix(Copy, *U_Graph_) );
147  D_ = Teuchos::rcp( new Epetra_Vector(*IlukRowMap_) );
148  }
149  else {
150  // Allocate Epetra_CrsMatrix using ILUK graphs
151  L_ = Teuchos::rcp( new Epetra_CrsMatrix(Copy, *IlukRowMap_, *IlukRowMap_, 0) );
152  U_ = Teuchos::rcp( new Epetra_CrsMatrix(Copy, *IlukRowMap_, *IlukRowMap_, 0) );
153  D_ = Teuchos::rcp( new Epetra_Vector(*IlukRowMap_) );
154  L_Graph_ = Teuchos::null;
155  U_Graph_ = Teuchos::null;
156  }
157  SetAllocated(true);
158  return(0);
159 }
160 
161 //==========================================================================
162 int Ifpack_CrsRiluk::SetParameters(const Teuchos::ParameterList& parameterlist,
163  bool cerr_warning_if_unused)
164 {
165  Ifpack::param_struct params;
166  params.double_params[Ifpack::relax_value] = RelaxValue_;
167  params.double_params[Ifpack::absolute_threshold] = Athresh_;
168  params.double_params[Ifpack::relative_threshold] = Rthresh_;
169  params.overlap_mode = OverlapMode_;
170 
171  Ifpack::set_parameters(parameterlist, params, cerr_warning_if_unused);
172 
173  RelaxValue_ = params.double_params[Ifpack::relax_value];
174  Athresh_ = params.double_params[Ifpack::absolute_threshold];
175  Rthresh_ = params.double_params[Ifpack::relative_threshold];
176  OverlapMode_ = params.overlap_mode;
177 
178  return(0);
179 }
180 
181 //==========================================================================
183 
184  UserMatrixIsCrs_ = true;
185 
186  if (!Allocated()) AllocateCrs();
187 
188  Teuchos::RefCountPtr<Epetra_CrsMatrix> OverlapA = Teuchos::rcp( (Epetra_CrsMatrix *) &A, false );
189 
190  if (IsOverlapped_) {
191 
192  OverlapA = Teuchos::rcp( new Epetra_CrsMatrix(Copy, *Graph_.OverlapGraph()) );
193  EPETRA_CHK_ERR(OverlapA->Import(A, *Graph_.OverlapImporter(), Insert));
194  EPETRA_CHK_ERR(OverlapA->FillComplete());
195  }
196 
197  // Get Maximun Row length
198  int MaxNumEntries = OverlapA->MaxNumEntries();
199 
200  // Set L range map and U domain map
201  U_DomainMap_ = Teuchos::rcp( &(A.DomainMap()), false );
202  L_RangeMap_ = Teuchos::rcp( &(A.RangeMap()), false );
203  // Do the rest using generic Epetra_RowMatrix interface
204 
205  EPETRA_CHK_ERR(InitAllValues(*OverlapA, MaxNumEntries));
206 
207  return(0);
208 }
209 
210 //==========================================================================
211 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES // FIXME LONG LONG
213 
214  UserMatrixIsVbr_ = true;
215 
216  if (!Allocated()) AllocateVbr();
217 
218  //cout << "Original Graph " << endl << A.Graph() << endl << flush;
219  //A.Comm().Barrier();
220  //if (A.Comm().MyPID()==0) cout << "*****************************************************" <<endl;
221  //cout << "Original Matrix " << endl << A << endl << flush;
222  //A.Comm().Barrier();
223  //if (A.Comm().MyPID()==0) cout << "*****************************************************" <<endl;
224  //cout << "Overlap Graph " << endl << *Graph_.OverlapGraph() << endl << flush;
225  //A.Comm().Barrier();
226  //if (A.Comm().MyPID()==0) cout << "*****************************************************" <<endl;
227 
228  Teuchos::RefCountPtr<Epetra_VbrMatrix> OverlapA = Teuchos::rcp( (Epetra_VbrMatrix *) &A, false );
229 
230  if (IsOverlapped_) {
231 
232  OverlapA = Teuchos::rcp( new Epetra_VbrMatrix(Copy, *Graph_.OverlapGraph()) );
233  EPETRA_CHK_ERR(OverlapA->Import(A, *Graph_.OverlapImporter(), Insert));
234  EPETRA_CHK_ERR(OverlapA->FillComplete());
235  }
236 
237  //cout << "Overlap Matrix " << endl << *OverlapA << endl << flush;
238 
239  // Get Maximun Row length
240  int MaxNumEntries = OverlapA->MaxNumNonzeros();
241 
242  // Do the rest using generic Epetra_RowMatrix interface
243 
244  EPETRA_CHK_ERR(InitAllValues(*OverlapA, MaxNumEntries));
245 
246  return(0);
247 }
248 #endif
249 //==========================================================================
250 
251 int Ifpack_CrsRiluk::InitAllValues(const Epetra_RowMatrix & OverlapA, int MaxNumEntries) {
252 
253  int ierr = 0;
254  int i, j;
255  int NumIn, NumL, NumU;
256  bool DiagFound;
257  int NumNonzeroDiags = 0;
258 
259 
260  std::vector<int> InI(MaxNumEntries); // Allocate temp space
261  std::vector<int> LI(MaxNumEntries);
262  std::vector<int> UI(MaxNumEntries);
263  std::vector<double> InV(MaxNumEntries);
264  std::vector<double> LV(MaxNumEntries);
265  std::vector<double> UV(MaxNumEntries);
266 
267  bool ReplaceValues = (L_->StaticGraph() || L_->IndicesAreLocal()); // Check if values should be inserted or replaced
268 
269  if (ReplaceValues) {
270  L_->PutScalar(0.0); // Zero out L and U matrices
271  U_->PutScalar(0.0);
272  }
273 
274  D_->PutScalar(0.0); // Set diagonal values to zero
275  double *DV;
276  EPETRA_CHK_ERR(D_->ExtractView(&DV)); // Get view of diagonal
277 
278 
279  // First we copy the user's matrix into L and U, regardless of fill level
280 
281  for (i=0; i< NumMyRows(); i++) {
282 
283  EPETRA_CHK_ERR(OverlapA.ExtractMyRowCopy(i, MaxNumEntries, NumIn, &InV[0], &InI[0])); // Get Values and Indices
284 
285  // Split into L and U (we don't assume that indices are ordered).
286 
287  NumL = 0;
288  NumU = 0;
289  DiagFound = false;
290 
291  for (j=0; j< NumIn; j++) {
292  int k = InI[j];
293 
294  if (k==i) {
295  DiagFound = true;
296  DV[i] += Rthresh_ * InV[j] + EPETRA_SGN(InV[j]) * Athresh_; // Store perturbed diagonal in Epetra_Vector D_
297  }
298 
299  else if (k < 0) {EPETRA_CHK_ERR(-1);} // Out of range
300 
301  else if (k < i) {
302  LI[NumL] = k;
303  LV[NumL] = InV[j];
304  NumL++;
305  }
306  else if (k<NumMyRows()) {
307  UI[NumU] = k;
308  UV[NumU] = InV[j];
309  NumU++;
310  }
311  }
312 
313  // Check in things for this row of L and U
314 
315  if (DiagFound) NumNonzeroDiags++;
316  else DV[i] = Athresh_;
317 
318  if (NumL) {
319  if (ReplaceValues) {
320  EPETRA_CHK_ERR(L_->ReplaceMyValues(i, NumL, &LV[0], &LI[0]));
321  }
322  else {
323  EPETRA_CHK_ERR(L_->InsertMyValues(i, NumL, &LV[0], &LI[0]));
324  }
325  }
326 
327  if (NumU) {
328  if (ReplaceValues) {
329  EPETRA_CHK_ERR(U_->ReplaceMyValues(i, NumU, &UV[0], &UI[0]));
330  }
331  else {
332  EPETRA_CHK_ERR(U_->InsertMyValues(i, NumU, &UV[0], &UI[0]));
333  }
334  }
335 
336  }
337 
338  if (!ReplaceValues) {
339  // The domain of L and the range of U are exactly their own row maps (there is no communication).
340  // The domain of U and the range of L must be the same as those of the original matrix,
341  // However if the original matrix is a VbrMatrix, these two latter maps are translation from
342  // a block map to a point map.
343  EPETRA_CHK_ERR(L_->FillComplete(L_->RowMatrixColMap(), *L_RangeMap_));
344  EPETRA_CHK_ERR(U_->FillComplete(*U_DomainMap_, U_->RowMatrixRowMap()));
345  }
346 
347  // At this point L and U have the values of A in the structure of L and U, and diagonal vector D
348 
349  SetValuesInitialized(true);
350  SetFactored(false);
351 
352  int TotalNonzeroDiags = 0;
353  EPETRA_CHK_ERR(Graph_.L_Graph().RowMap().Comm().SumAll(&NumNonzeroDiags, &TotalNonzeroDiags, 1));
354  NumMyDiagonals_ = NumNonzeroDiags;
355  if (NumNonzeroDiags != NumMyRows()) ierr = 1; // Diagonals are not right, warn user
356 
357  return(ierr);
358 }
359 
360 //==========================================================================
362 
363  // if (!Allocated()) return(-1); // This test is not needed at this time. All constructors allocate.
364  if (!ValuesInitialized()) return(-2); // Must have values initialized.
365  if (Factored()) return(-3); // Can't have already computed factors.
366 
367  SetValuesInitialized(false);
368 
369  // MinMachNum should be officially defined, for now pick something a little
370  // bigger than IEEE underflow value
371 
372  double MinDiagonalValue = Epetra_MinDouble;
373  double MaxDiagonalValue = 1.0/MinDiagonalValue;
374 
375  int ierr = 0;
376  int i, j, k;
377  int * LI=0, * UI = 0;
378  double * LV=0, * UV = 0;
379  int NumIn, NumL, NumU;
380 
381  // Get Maximun Row length
382  int MaxNumEntries = L_->MaxNumEntries() + U_->MaxNumEntries() + 1;
383 
384  std::vector<int> InI(MaxNumEntries); // Allocate temp space
385  std::vector<double> InV(MaxNumEntries);
386  std::vector<int> colflag(NumMyCols());
387 
388  double *DV;
389  ierr = D_->ExtractView(&DV); // Get view of diagonal
390 
391 #ifdef IFPACK_FLOPCOUNTERS
392  int current_madds = 0; // We will count multiply-add as they happen
393 #endif
394 
395  // Now start the factorization.
396 
397  // Need some integer workspace and pointers
398  int NumUU;
399  int * UUI;
400  double * UUV;
401  for (j=0; j<NumMyCols(); j++) colflag[j] = - 1;
402 
403  for(i=0; i<NumMyRows(); i++) {
404 
405  // Fill InV, InI with current row of L, D and U combined
406 
407  NumIn = MaxNumEntries;
408  EPETRA_CHK_ERR(L_->ExtractMyRowCopy(i, NumIn, NumL, &InV[0], &InI[0]));
409  LV = &InV[0];
410  LI = &InI[0];
411 
412  InV[NumL] = DV[i]; // Put in diagonal
413  InI[NumL] = i;
414 
415  EPETRA_CHK_ERR(U_->ExtractMyRowCopy(i, NumIn-NumL-1, NumU, &InV[NumL+1], &InI[NumL+1]));
416  NumIn = NumL+NumU+1;
417  UV = &InV[NumL+1];
418  UI = &InI[NumL+1];
419 
420  // Set column flags
421  for (j=0; j<NumIn; j++) colflag[InI[j]] = j;
422 
423  double diagmod = 0.0; // Off-diagonal accumulator
424 
425  for (int jj=0; jj<NumL; jj++) {
426  j = InI[jj];
427  double multiplier = InV[jj]; // current_mults++;
428 
429  InV[jj] *= DV[j];
430 
431  EPETRA_CHK_ERR(U_->ExtractMyRowView(j, NumUU, UUV, UUI)); // View of row above
432 
433  if (RelaxValue_==0.0) {
434  for (k=0; k<NumUU; k++) {
435  int kk = colflag[UUI[k]];
436  if (kk>-1) {
437  InV[kk] -= multiplier*UUV[k];
438 #ifdef IFPACK_FLOPCOUNTERS
439  current_madds++;
440 #endif
441  }
442  }
443  }
444  else {
445  for (k=0; k<NumUU; k++) {
446  int kk = colflag[UUI[k]];
447  if (kk>-1) InV[kk] -= multiplier*UUV[k];
448  else diagmod -= multiplier*UUV[k];
449 #ifdef IFPACK_FLOPCOUNTERS
450  current_madds++;
451 #endif
452  }
453  }
454  }
455  if (NumL) {
456  EPETRA_CHK_ERR(L_->ReplaceMyValues(i, NumL, LV, LI)); // Replace current row of L
457  }
458 
459  DV[i] = InV[NumL]; // Extract Diagonal value
460 
461  if (RelaxValue_!=0.0) {
462  DV[i] += RelaxValue_*diagmod; // Add off diagonal modifications
463  // current_madds++;
464  }
465 
466  if (fabs(DV[i]) > MaxDiagonalValue) {
467  if (DV[i] < 0) DV[i] = - MinDiagonalValue;
468  else DV[i] = MinDiagonalValue;
469  }
470  else
471  DV[i] = 1.0/DV[i]; // Invert diagonal value
472 
473  for (j=0; j<NumU; j++) UV[j] *= DV[i]; // Scale U by inverse of diagonal
474 
475  if (NumU) {
476  EPETRA_CHK_ERR(U_->ReplaceMyValues(i, NumU, UV, UI)); // Replace current row of L and U
477  }
478 
479  // Reset column flags
480  for (j=0; j<NumIn; j++) colflag[InI[j]] = -1;
481  }
482 
483  // Validate that the L and U factors are actually lower and upper triangular
484 
485  if( !L_->LowerTriangular() )
486  EPETRA_CHK_ERR(-2);
487  if( !U_->UpperTriangular() )
488  EPETRA_CHK_ERR(-3);
489 
490 #ifdef IFPACK_FLOPCOUNTERS
491  // Add up flops
492 
493  double current_flops = 2 * current_madds;
494  double total_flops = 0;
495 
496  EPETRA_CHK_ERR(Graph_.L_Graph().RowMap().Comm().SumAll(&current_flops, &total_flops, 1)); // Get total madds across all PEs
497 
498  // Now count the rest
499  total_flops += (double) L_->NumGlobalNonzeros(); // Accounts for multiplier above
500  total_flops += (double) D_->GlobalLength(); // Accounts for reciprocal of diagonal
501  if (RelaxValue_!=0.0) total_flops += 2 * (double)D_->GlobalLength(); // Accounts for relax update of diag
502 
503  UpdateFlops(total_flops); // Update flop count
504 #endif
505 
506  SetFactored(true);
507 
508  return(ierr);
509 
510 }
511 
512 //=============================================================================
514  Epetra_MultiVector& Y) const {
515 //
516 // This function finds Y such that LDU Y = X or U(trans) D L(trans) Y = X for multiple RHS
517 //
518 
519  // First generate X and Y as needed for this function
520  Teuchos::RefCountPtr<Epetra_MultiVector> X1;
521  Teuchos::RefCountPtr<Epetra_MultiVector> Y1;
522  EPETRA_CHK_ERR(GenerateXY(Trans, X, Y, &X1, &Y1));
523 
524  bool Upper = true;
525  bool Lower = false;
526  bool UnitDiagonal = true;
527 
528 #ifdef IFPACK_FLOPCOUNTERS
529  Epetra_Flops * counter = this->GetFlopCounter();
530  if (counter!=0) {
531  L_->SetFlopCounter(*counter);
532  Y1->SetFlopCounter(*counter);
533  U_->SetFlopCounter(*counter);
534  }
535 #endif
536 
537  if (!Trans) {
538 
539  EPETRA_CHK_ERR(L_->Solve(Lower, Trans, UnitDiagonal, *X1, *Y1));
540  EPETRA_CHK_ERR(Y1->Multiply(1.0, *D_, *Y1, 0.0)); // y = D*y (D_ has inverse of diagonal)
541  EPETRA_CHK_ERR(U_->Solve(Upper, Trans, UnitDiagonal, *Y1, *Y1)); // Solve Uy = y
542  if (IsOverlapped_) {EPETRA_CHK_ERR(Y.Export(*Y1,*L_->Exporter(), OverlapMode_));} // Export computed Y values if needed
543  }
544  else {
545  EPETRA_CHK_ERR(U_->Solve(Upper, Trans, UnitDiagonal, *X1, *Y1)); // Solve Uy = y
546  EPETRA_CHK_ERR(Y1->Multiply(1.0, *D_, *Y1, 0.0)); // y = D*y (D_ has inverse of diagonal)
547  EPETRA_CHK_ERR(L_->Solve(Lower, Trans, UnitDiagonal, *Y1, *Y1));
548  if (IsOverlapped_) {EPETRA_CHK_ERR(Y.Export(*Y1,*U_->Importer(), OverlapMode_));} // Export computed Y values if needed
549  }
550 
551  return(0);
552 }
553 //=============================================================================
555  Epetra_MultiVector& Y) const {
556 //
557 // This function finds X such that LDU Y = X or U(trans) D L(trans) Y = X for multiple RHS
558 //
559 
560  // First generate X and Y as needed for this function
561  Teuchos::RefCountPtr<Epetra_MultiVector> X1;
562  Teuchos::RefCountPtr<Epetra_MultiVector> Y1;
563  EPETRA_CHK_ERR(GenerateXY(Trans, X, Y, &X1, &Y1));
564 
565 #ifdef IFPACK_FLOPCOUNTERS
566  Epetra_Flops * counter = this->GetFlopCounter();
567  if (counter!=0) {
568  L_->SetFlopCounter(*counter);
569  Y1->SetFlopCounter(*counter);
570  U_->SetFlopCounter(*counter);
571  }
572 #endif
573 
574  if (!Trans) {
575  EPETRA_CHK_ERR(U_->Multiply(Trans, *X1, *Y1)); //
576  EPETRA_CHK_ERR(Y1->Update(1.0, *X1, 1.0)); // Y1 = Y1 + X1 (account for implicit unit diagonal)
577  EPETRA_CHK_ERR(Y1->ReciprocalMultiply(1.0, *D_, *Y1, 0.0)); // y = D*y (D_ has inverse of diagonal)
578  Epetra_MultiVector Y1temp(*Y1); // Need a temp copy of Y1
579  EPETRA_CHK_ERR(L_->Multiply(Trans, Y1temp, *Y1));
580  EPETRA_CHK_ERR(Y1->Update(1.0, Y1temp, 1.0)); // (account for implicit unit diagonal)
581  if (IsOverlapped_) {EPETRA_CHK_ERR(Y.Export(*Y1,*L_->Exporter(), OverlapMode_));} // Export computed Y values if needed
582  }
583  else {
584 
585  EPETRA_CHK_ERR(L_->Multiply(Trans, *X1, *Y1));
586  EPETRA_CHK_ERR(Y1->Update(1.0, *X1, 1.0)); // Y1 = Y1 + X1 (account for implicit unit diagonal)
587  EPETRA_CHK_ERR(Y1->ReciprocalMultiply(1.0, *D_, *Y1, 0.0)); // y = D*y (D_ has inverse of diagonal)
588  Epetra_MultiVector Y1temp(*Y1); // Need a temp copy of Y1
589  EPETRA_CHK_ERR(U_->Multiply(Trans, Y1temp, *Y1));
590  EPETRA_CHK_ERR(Y1->Update(1.0, Y1temp, 1.0)); // (account for implicit unit diagonal)
591  if (IsOverlapped_) {EPETRA_CHK_ERR(Y.Export(*Y1,*L_->Exporter(), OverlapMode_));}
592  }
593  return(0);
594 }
595 //=============================================================================
596 int Ifpack_CrsRiluk::Condest(bool Trans, double & ConditionNumberEstimate) const {
597 
598  if (Condest_>=0.0) {
599  ConditionNumberEstimate = Condest_;
600  return(0);
601  }
602  // Create a vector with all values equal to one
603  Epetra_Vector Ones(U_->DomainMap());
604  Epetra_Vector OnesResult(L_->RangeMap());
605  Ones.PutScalar(1.0);
606 
607  EPETRA_CHK_ERR(Solve(Trans, Ones, OnesResult)); // Compute the effect of the solve on the vector of ones
608  EPETRA_CHK_ERR(OnesResult.Abs(OnesResult)); // Make all values non-negative
609  EPETRA_CHK_ERR(OnesResult.MaxValue(&ConditionNumberEstimate)); // Get the maximum value across all processors
610  Condest_ = ConditionNumberEstimate; // Save value for possible later calls
611  return(0);
612 }
613 //==============================================================================
614 int Ifpack_CrsRiluk::BlockGraph2PointGraph(const Epetra_CrsGraph & BG, Epetra_CrsGraph & PG, bool Upper) {
615 
616  if (!BG.IndicesAreLocal()) {EPETRA_CHK_ERR(-1);} // Must have done FillComplete on BG
617 
618  int * ColFirstPointInElementList = BG.RowMap().FirstPointInElementList();
619  int * ColElementSizeList = BG.RowMap().ElementSizeList();
620  if (BG.Importer()!=0) {
621  ColFirstPointInElementList = BG.ImportMap().FirstPointInElementList();
622  ColElementSizeList = BG.ImportMap().ElementSizeList();
623  }
624 
625  int Length = (BG.MaxNumIndices()+1) * BG.ImportMap().MaxMyElementSize();
626  std::vector<int> tmpIndices(Length);
627 
628  int BlockRow, BlockOffset, NumEntries;
629  int NumBlockEntries;
630  int * BlockIndices;
631 
632  int NumMyRows_tmp = PG.NumMyRows();
633 
634  for (int i=0; i<NumMyRows_tmp; i++) {
635  EPETRA_CHK_ERR(BG.RowMap().FindLocalElementID(i, BlockRow, BlockOffset));
636  EPETRA_CHK_ERR(BG.ExtractMyRowView(BlockRow, NumBlockEntries, BlockIndices));
637 
638  int * ptr = &tmpIndices[0]; // Set pointer to beginning of buffer
639 
640  int RowDim = BG.RowMap().ElementSize(BlockRow);
641  NumEntries = 0;
642 
643  // This next line make sure that the off-diagonal entries in the block diagonal of the
644  // original block entry matrix are included in the nonzero pattern of the point graph
645  if (Upper) {
646  int jstart = i+1;
647  int jstop = EPETRA_MIN(NumMyRows_tmp,i+RowDim-BlockOffset);
648  for (int j= jstart; j< jstop; j++) {*ptr++ = j; NumEntries++;}
649  }
650 
651  for (int j=0; j<NumBlockEntries; j++) {
652  int ColDim = ColElementSizeList[BlockIndices[j]];
653  NumEntries += ColDim;
654  assert(NumEntries<=Length); // Sanity test
655  int Index = ColFirstPointInElementList[BlockIndices[j]];
656  for (int k=0; k < ColDim; k++) *ptr++ = Index++;
657  }
658 
659  // This next line make sure that the off-diagonal entries in the block diagonal of the
660  // original block entry matrix are included in the nonzero pattern of the point graph
661  if (!Upper) {
662  int jstart = EPETRA_MAX(0,i-RowDim+1);
663  int jstop = i;
664  for (int j = jstart; j < jstop; j++) {*ptr++ = j; NumEntries++;}
665  }
666 
667  EPETRA_CHK_ERR(PG.InsertMyIndices(i, NumEntries, &tmpIndices[0]));
668  }
669 
670  SetAllocated(true);
671 
672  return(0);
673 }
674 //=========================================================================
675 int Ifpack_CrsRiluk::BlockMap2PointMap(const Epetra_BlockMap & BlockMap, Teuchos::RefCountPtr<Epetra_Map>* PointMap) {
676  // Generate an Epetra_Map that has the same number and distribution of points
677  // as the input Epetra_BlockMap object. The global IDs for the output PointMap
678  // are computed by using the MaxElementSize of the BlockMap. For variable block
679  // sizes this will create gaps in the GID space, but that is OK for Epetra_Maps.
680 
681  int MaxElementSize = BlockMap.MaxElementSize();
682  int PtNumMyElements = BlockMap.NumMyPoints();
683 
684  std::vector<int> PtMyGlobalElements_int;
685  std::vector<long long> PtMyGlobalElements_LL;
686 
687  int NumMyElements = BlockMap.NumMyElements();
688  int curID = 0;
689 
690  if (PtNumMyElements>0) {
691 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
692  if(BlockMap.GlobalIndicesInt()) {
693  PtMyGlobalElements_int.resize(PtNumMyElements);
694  for (int i=0; i<NumMyElements; i++) {
695  int StartID = BlockMap.GID(i)*MaxElementSize;
696  int ElementSize = BlockMap.ElementSize(i);
697  for (int j=0; j<ElementSize; j++) PtMyGlobalElements_int[curID++] = StartID+j;
698  }
699  }
700  else
701 #endif
702 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
703  if(BlockMap.GlobalIndicesLongLong()) {
704  PtMyGlobalElements_LL.resize(PtNumMyElements);
705  for (int i=0; i<NumMyElements; i++) {
706  long long StartID = BlockMap.GID64(i)*MaxElementSize;
707  int ElementSize = BlockMap.ElementSize(i);
708  for (int j=0; j<ElementSize; j++) PtMyGlobalElements_LL[curID++] = StartID+j;
709  }
710  }
711  else
712 #endif
713  throw "Ifpack_CrsRiluk::BlockMap2PointMap: GlobalIndices type unknown";
714  }
715 
716  assert(curID==PtNumMyElements); // Sanity test
717 
718 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
719  if(BlockMap.GlobalIndicesInt())
720  (*PointMap) = Teuchos::rcp( new Epetra_Map(-1, PtNumMyElements, &PtMyGlobalElements_int[0], BlockMap.IndexBase(), BlockMap.Comm()) );
721  else
722 #endif
723 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
724  if(BlockMap.GlobalIndicesLongLong())
725  (*PointMap) = Teuchos::rcp( new Epetra_Map(-1LL, PtNumMyElements, &PtMyGlobalElements_LL[0], BlockMap.IndexBase64(), BlockMap.Comm()) );
726  else
727 #endif
728  throw "Ifpack_CrsRiluk::BlockMap2PointMap: GlobalIndices type unknown";
729 
730  if (!BlockMap.PointSameAs(*(*PointMap))) {EPETRA_CHK_ERR(-1);} // Maps not compatible
731  return(0);
732 }
733 //=========================================================================
734 int Ifpack_CrsRiluk::GenerateXY(bool Trans,
735  const Epetra_MultiVector& Xin, const Epetra_MultiVector& Yin,
736  Teuchos::RefCountPtr<Epetra_MultiVector>* Xout,
737  Teuchos::RefCountPtr<Epetra_MultiVector>* Yout) const {
738 
739  // Generate an X and Y suitable for performing Solve() and Multiply() methods
740 
741  if (Xin.NumVectors()!=Yin.NumVectors()) EPETRA_CHK_ERR(-1); // Return error: X and Y not the same size
742 
743  //cout << "Xin = " << Xin << endl;
744  (*Xout) = Teuchos::rcp( (Epetra_MultiVector *) &Xin, false );
745  (*Yout) = Teuchos::rcp( (Epetra_MultiVector *) &Yin, false );
746  if (!IsOverlapped_ && UserMatrixIsCrs_) return(0); // Nothing more to do
747 
748  if (UserMatrixIsVbr_) {
749  if (VbrX_!=Teuchos::null) {
750  if (VbrX_->NumVectors()!=Xin.NumVectors()) {
751  VbrX_ = Teuchos::null;
752  VbrY_ = Teuchos::null;
753  }
754  }
755  if (VbrX_==Teuchos::null) { // Need to allocate space for overlap X and Y
756  VbrX_ = Teuchos::rcp( new Epetra_MultiVector(View, *U_DomainMap_, (*Xout)->Pointers(), (*Xout)->NumVectors()) );
757  VbrY_ = Teuchos::rcp( new Epetra_MultiVector(View, *L_RangeMap_, (*Yout)->Pointers(), (*Yout)->NumVectors()) );
758  }
759  else {
760  EPETRA_CHK_ERR(VbrX_->ResetView((*Xout)->Pointers()));
761  EPETRA_CHK_ERR(VbrY_->ResetView((*Yout)->Pointers()));
762  }
763  (*Xout) = VbrX_;
764  (*Yout) = VbrY_;
765  }
766 
767  if (IsOverlapped_) {
768  // Make sure the number of vectors in the multivector is the same as before.
769  if (OverlapX_!=Teuchos::null) {
770  if (OverlapX_->NumVectors()!=Xin.NumVectors()) {
771  OverlapX_ = Teuchos::null;
772  OverlapY_ = Teuchos::null;
773  }
774  }
775  if (OverlapX_==Teuchos::null) { // Need to allocate space for overlap X and Y
776  OverlapX_ = Teuchos::rcp( new Epetra_MultiVector(U_->RowMatrixColMap(), (*Xout)->NumVectors()) );
777  OverlapY_ = Teuchos::rcp( new Epetra_MultiVector(L_->RowMatrixRowMap(), (*Yout)->NumVectors()) );
778  }
779  if (!Trans) {
780  EPETRA_CHK_ERR(OverlapX_->Import(*(*Xout),*U_->Importer(), Insert)); // Import X values for solve
781  }
782  else {
783  EPETRA_CHK_ERR(OverlapX_->Import(*(*Xout),*L_->Exporter(), Insert)); // Import X values for solve
784  }
785  (*Xout) = OverlapX_;
786  (*Yout) = OverlapY_; // Set pointers for Xout and Yout to point to overlap space
787  //cout << "OverlapX_ = " << *OverlapX_ << endl;
788  }
789 
790  return(0);
791 }
792 //=============================================================================
793 // Non-member functions
794 
795 std::ostream& operator << (std::ostream& os, const Ifpack_CrsRiluk& A)
796 {
797  using std::endl;
798 
799 /* Epetra_fmtflags olda = os.setf(ios::right,ios::adjustfield);
800  Epetra_fmtflags oldf = os.setf(ios::scientific,ios::floatfield);
801  int oldp = os.precision(12); */
802  int LevelFill = A.Graph().LevelFill();
803  int LevelOverlap = A.Graph().LevelOverlap();
804  Epetra_CrsMatrix & L = (Epetra_CrsMatrix &) A.L();
805  Epetra_CrsMatrix & U = (Epetra_CrsMatrix &) A.U();
806  Epetra_Vector & D = (Epetra_Vector &) A.D();
807 
808  os.width(14);
809  os << endl;
810  os << " Level of Fill = "; os << LevelFill;
811  os << endl;
812  os.width(14);
813  os << " Level of Overlap = "; os << LevelOverlap;
814  os << endl;
815 
816  os.width(14);
817  os << " Lower Triangle = ";
818  os << endl;
819  os << L; // Let Epetra_CrsMatrix handle the rest.
820  os << endl;
821 
822  os.width(14);
823  os << " Inverse of Diagonal = ";
824  os << endl;
825  os << D; // Let Epetra_Vector handle the rest.
826  os << endl;
827 
828  os.width(14);
829  os << " Upper Triangle = ";
830  os << endl;
831  os << U; // Let Epetra_CrsMatrix handle the rest.
832  os << endl;
833 
834  // Reset os flags
835 
836 /* os.setf(olda,ios::adjustfield);
837  os.setf(oldf,ios::floatfield);
838  os.precision(oldp); */
839 
840  return os;
841 }
bool PointSameAs(const Epetra_BlockMap &Map) const
const Epetra_BlockMap & RangeMap() const
int Multiply(bool Trans, const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Returns the result of multiplying U, D and L in that order on an Epetra_MultiVector X in Y...
bool DistributedGlobal() const
const Epetra_Map & RangeMap() const
int ElementSize() const
Ifpack_CrsRiluk: A class for constructing and using an incomplete lower/upper (ILU) factorization of ...
int Condest(bool Trans, double &ConditionNumberEstimate) const
Returns the maximum over all the condition number estimate for each local ILU set of factors...
bool GlobalIndicesLongLong() const
const Epetra_CrsMatrix & L() const
Returns the address of the L factor associated with this factored matrix.
virtual Epetra_CrsGraph & L_Graph()
Returns the graph of lower triangle of the ILU(k) graph as a Epetra_CrsGraph.
virtual Epetra_CrsGraph * OverlapGraph() const
Returns the the overlapped graph.
int InsertMyIndices(int LocalRow, int NumIndices, int *Indices)
const Epetra_Import * Importer() const
int NumMyRows() const
Ifpack_IlukGraph: A class for constructing level filled graphs for use with ILU(k) class precondition...
const Epetra_BlockMap & ImportMap() const
bool GlobalIndicesInt() const
const Epetra_BlockMap & DomainMap() const
int * FirstPointInElementList() const
virtual const Epetra_BlockMap & DomainMap() const
Returns the Epetra_BlockMap object associated with the domain of this matrix operator.
int IndexBase() const
virtual int LevelFill() const
Returns the level of fill used to construct this graph.
virtual int SumAll(double *PartialSums, double *GlobalSums, int Count) const =0
int NumMyElements() const
virtual Epetra_CrsGraph & U_Graph()
Returns the graph of upper triangle of the ILU(k) graph as a Epetra_CrsGraph.
int SetParameters(const Teuchos::ParameterList &parameterlist, bool cerr_warning_if_unused=false)
Set parameters using a Teuchos::ParameterList object.
int GID(int LID) const
int NumMyCols() const
Returns the number of local matrix columns.
int InitValues(const Epetra_CrsMatrix &A)
Initialize L and U with values from user matrix A.
const Epetra_CrsMatrix & U() const
Returns the address of the L factor associated with this factored matrix.
int * ElementSizeList() const
const Epetra_BlockMap & RowMap() const
int FindLocalElementID(int PointID, int &ElementID, int &ElementOffset) const
virtual int LevelOverlap() const
Returns the level of overlap used to construct this graph.
const Epetra_Comm & Comm() const
virtual Epetra_Import * OverlapImporter() const
Returns the importer used to create the overlapped graph.
int ExtractMyRowView(int LocalRow, int &NumIndices, int *&Indices) const
int Solve(bool Trans, const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Returns the result of a Ifpack_CrsRiluk forward/back solve on a Epetra_MultiVector X in Y (works for ...
bool Factored() const
If factor is completed, this query returns true, otherwise it returns false.
int NumMyRows() const
Returns the number of local matrix rows.
bool ValuesInitialized() const
If values have been initialized, this query returns true, otherwise it returns false.
const Epetra_Vector & D() const
Returns the address of the D factor associated with this factored matrix.
int MaxNumIndices() const
const Epetra_Map & DomainMap() const
int MaxElementSize() const
int MaxMyElementSize() const
virtual int ExtractMyRowCopy(int MyRow, int Length, int &NumEntries, double *Values, int *Indices) const =0
int NumMyPoints() const
virtual ~Ifpack_CrsRiluk()
Ifpack_CrsRiluk Destructor.
int Factor()
Compute ILU factors L and U using the specified graph, diagonal perturbation thresholds and relaxatio...
const Ifpack_IlukGraph & Graph() const
returns the address of the Ifpack_IlukGraph associated with this factored matrix. ...
Ifpack_CrsRiluk(const Ifpack_IlukGraph &Graph_in)
Ifpack_CrsRiluk constuctor with variable number of indices per row.
bool IndicesAreLocal() const