Ifpack Package Browser (Single Doxygen Collection)  Development
Ifpack_CrsRiluk.cpp
Go to the documentation of this file.
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 //==============================================================================
113 
114  // Allocate Epetra_CrsMatrix using ILUK graphs
118  L_Graph_ = Teuchos::null;
119  U_Graph_ = Teuchos::null;
120  SetAllocated(true);
121  return(0);
122 }
123 //==============================================================================
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.
131 
132  // Set L range map and U domain map
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()) {
141 
142  L_Graph_->FillComplete(*IlukRowMap_, *IlukRangeMap_);
143  U_Graph_->FillComplete(*IlukDomainMap_, *IlukRowMap_);
144 
148  }
149  else {
150  // Allocate Epetra_CrsMatrix using ILUK graphs
154  L_Graph_ = Teuchos::null;
155  U_Graph_ = Teuchos::null;
156  }
157  SetAllocated(true);
158  return(0);
159 }
160 
161 //==========================================================================
163  bool cerr_warning_if_unused)
164 {
165  Ifpack::param_struct params;
169  params.overlap_mode = OverlapMode_;
170 
171  Ifpack::set_parameters(parameterlist, params, cerr_warning_if_unused);
172 
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 //==============================================================================
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 //=========================================================================
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 }
const Epetra_Import * Importer() const
Teuchos::RefCountPtr< const Epetra_Map > L_RangeMap_
void UpdateFlops(int Flops_in) const
int NumMyRows() const
Returns the number of local matrix rows.
int BlockMap2PointMap(const Epetra_BlockMap &BlockMap, Teuchos::RefCountPtr< Epetra_Map > *PointMap)
Ifpack_CrsRiluk: A class for constructing and using an incomplete lower/upper (ILU) factorization of ...
void SetFactored(bool Flag)
double double_params[FIRST_INT_PARAM]
virtual Epetra_CrsGraph * OverlapGraph() const
Returns the the overlapped graph.
bool DistributedGlobal() const
Teuchos::RefCountPtr< Epetra_CrsGraph > U_Graph_
bool Allocated() const
Teuchos::RefCountPtr< const Epetra_Map > U_DomainMap_
const Epetra_CrsMatrix & L() const
Returns the address of the L factor associated with this factored matrix.
Teuchos::RefCountPtr< Epetra_MultiVector > OverlapY_
virtual Epetra_CrsGraph & L_Graph()
Returns the graph of lower triangle of the ILU(k) graph as a Epetra_CrsGraph.
Epetra_Flops * GetFlopCounter() const
int PutScalar(double ScalarConstant)
int InsertMyIndices(int LocalRow, int NumIndices, int *Indices)
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 ...
#define EPETRA_CHK_ERR(a)
int IndexBase() const
Zero
virtual int LevelOverlap() const
Returns the level of overlap used to construct this graph.
Ifpack_IlukGraph: A class for constructing level filled graphs for use with ILU(k) class precondition...
Insert
int * FirstPointInElementList() const
#define EPETRA_MIN(x, y)
bool ValuesInitialized() const
If values have been initialized, this query returns true, otherwise it returns false.
int BlockGraph2PointGraph(const Epetra_CrsGraph &BG, Epetra_CrsGraph &PG, bool Upper)
const Epetra_BlockMap & RowMap() const
int SetAllocated(bool Flag)
int * ElementSizeList() const
int Condest(bool Trans, double &ConditionNumberEstimate) const
Returns the maximum over all the condition number estimate for each local ILU set of factors...
int NumMyRows() const
void SetValuesInitialized(bool Flag)
virtual int SumAll(double *PartialSums, double *GlobalSums, int Count) const=0
std::ostream & operator<<(std::ostream &os, const Ifpack_CrsRiluk &A)
<< operator will work for Ifpack_CrsRiluk.
virtual Epetra_CrsGraph & U_Graph()
Returns the graph of upper triangle of the ILU(k) graph as a Epetra_CrsGraph.
int Export(const Epetra_SrcDistObject &A, const Epetra_Import &Importer, Epetra_CombineMode CombineMode, const Epetra_OffsetIndex *Indexor=0)
Teuchos::RefCountPtr< Epetra_CrsGraph > L_Graph_
Teuchos::RefCountPtr< Epetra_Map > IlukDomainMap_
Epetra_CombineMode overlap_mode
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
long long GID64(int LID) const
const Epetra_BlockMap & RangeMap() const
const Ifpack_IlukGraph & Graph() const
returns the address of the Ifpack_IlukGraph associated with this factored matrix. ...
Teuchos::RefCountPtr< Epetra_CrsMatrix > U_
#define EPETRA_SGN(x)
int FindLocalElementID(int PointID, int &ElementID, int &ElementOffset) const
int InitAllValues(const Epetra_RowMatrix &A, int MaxNumEntries)
View
int NumVectors() const
int SetParameters(const Teuchos::ParameterList &parameterlist, bool cerr_warning_if_unused=false)
Set parameters using a Teuchos::ParameterList object.
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 GlobalIndicesInt() const
int NumMyElements() const
const Ifpack_IlukGraph & Graph_
int NumMyCols() const
Returns the number of local matrix columns.
Teuchos::RefCountPtr< Epetra_Map > IlukRowMap_
virtual const Epetra_BlockMap & DomainMap() const
Returns the Epetra_BlockMap object associated with the domain of this matrix operator.
int ElementSize() const
Teuchos::RefCountPtr< Epetra_MultiVector > OverlapX_
int InitValues(const Epetra_CrsMatrix &A)
Initialize L and U with values from user matrix A.
const Epetra_Comm & Comm() const
bool IndicesAreLocal() const
#define false
Teuchos::RefCountPtr< Epetra_CrsMatrix > L_
int MaxMyElementSize() const
const double Epetra_MinDouble
const Epetra_BlockMap & ImportMap() const
#define LL
long long IndexBase64() const
const Epetra_BlockMap & DomainMap() const
Teuchos::RefCountPtr< Epetra_MultiVector > VbrY_
int MaxNumIndices() const
int MaxElementSize() const
Copy
const Epetra_CrsMatrix & U() const
Returns the address of the L factor associated with this factored matrix.
virtual int ExtractMyRowCopy(int MyRow, int Length, int &NumEntries, double *Values, int *Indices) const=0
int ExtractMyRowView(int LocalRow, int &NumIndices, int *&Indices) const
bool Factored() const
If factor is completed, this query returns true, otherwise it returns false.
int GenerateXY(bool Trans, const Epetra_MultiVector &Xin, const Epetra_MultiVector &Yin, Teuchos::RefCountPtr< Epetra_MultiVector > *Xout, Teuchos::RefCountPtr< Epetra_MultiVector > *Yout) const
void set_parameters(const Teuchos::ParameterList &parameterlist, param_struct &params, bool cerr_warning_if_unused)
int NumMyPoints() const
Teuchos::RefCountPtr< Epetra_Vector > D_
Teuchos::RefCountPtr< Epetra_MultiVector > VbrX_
virtual ~Ifpack_CrsRiluk()
Ifpack_CrsRiluk Destructor.
bool GlobalIndicesLongLong() const
int GID(int LID) const
const Epetra_Vector & D() const
Returns the address of the D factor associated with this factored matrix.
bool PointSameAs(const Epetra_BlockMap &Map) const
Epetra_CombineMode OverlapMode_
#define EPETRA_MAX(x, y)
Teuchos::RefCountPtr< Epetra_Map > IlukRangeMap_
int Factor()
Compute ILU factors L and U using the specified graph, diagonal perturbation thresholds and relaxatio...
Ifpack_CrsRiluk(const Ifpack_IlukGraph &Graph_in)
Ifpack_CrsRiluk constuctor with variable number of indices per row.
virtual Epetra_Import * OverlapImporter() const
Returns the importer used to create the overlapped graph.