IFPACK Development
Loading...
Searching...
No Matches
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//==============================================================================
112int 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//==============================================================================
124int 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//==========================================================================
162int Ifpack_CrsRiluk::SetParameters(const Teuchos::ParameterList& parameterlist,
163 bool cerr_warning_if_unused)
164{
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
251int 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//=============================================================================
596int 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//==============================================================================
614int 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//=========================================================================
675int 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//=========================================================================
734int 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
795std::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}
Insert
Zero
Copy
int MaxElementSize() const
bool PointSameAs(const Epetra_BlockMap &Map) const
int * ElementSizeList() const
bool DistributedGlobal() const
int NumMyPoints() const
int GID(int LID) const
int IndexBase() const
int ElementSize() const
int * FirstPointInElementList() const
int MaxMyElementSize() const
bool GlobalIndicesInt() const
int FindLocalElementID(int PointID, int &ElementID, int &ElementOffset) const
const Epetra_Comm & Comm() const
int NumMyElements() const
bool GlobalIndicesLongLong() const
virtual int SumAll(double *PartialSums, double *GlobalSums, int Count) const=0
Epetra_Flops * GetFlopCounter() const
void UpdateFlops(int Flops_in) const
int ExtractMyRowView(int LocalRow, int &NumIndices, int *&Indices) const
int MaxNumIndices() const
const Epetra_BlockMap & DomainMap() const
bool IndicesAreLocal() const
const Epetra_BlockMap & RowMap() const
int InsertMyIndices(int LocalRow, int NumIndices, int *Indices)
int NumMyRows() const
const Epetra_BlockMap & ImportMap() const
const Epetra_BlockMap & RangeMap() const
const Epetra_Import * Importer() const
const Epetra_Map & RangeMap() const
const Epetra_Map & DomainMap() const
int Export(const Epetra_SrcDistObject &A, const Epetra_Import &Importer, Epetra_CombineMode CombineMode, const Epetra_OffsetIndex *Indexor=0)
int NumVectors() const
int Abs(const Epetra_MultiVector &A)
int MaxValue(double *Result) const
int PutScalar(double ScalarConstant)
virtual int ExtractMyRowCopy(int MyRow, int Length, int &NumEntries, double *Values, int *Indices) const=0
Ifpack_CrsRiluk: A class for constructing and using an incomplete lower/upper (ILU) factorization of ...
const Epetra_Vector & D() const
Returns the address of the D factor associated with this factored matrix.
virtual ~Ifpack_CrsRiluk()
Ifpack_CrsRiluk Destructor.
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.
int Factor()
Compute ILU factors L and U using the specified graph, diagonal perturbation thresholds and relaxatio...
int SetParameters(const Teuchos::ParameterList &parameterlist, bool cerr_warning_if_unused=false)
Set parameters using a Teuchos::ParameterList object.
int NumMyRows() const
Returns the number of local matrix rows.
const Epetra_CrsMatrix & U() const
Returns the address of the L factor associated with this factored matrix.
bool Factored() const
If factor is completed, this query returns true, otherwise it returns false.
int InitValues(const Epetra_CrsMatrix &A)
Initialize L and U with values from user matrix A.
Ifpack_CrsRiluk(const Ifpack_IlukGraph &Graph_in)
Ifpack_CrsRiluk constuctor with variable number of indices per row.
bool ValuesInitialized() const
If values have been initialized, this query returns true, otherwise it returns false.
const Epetra_CrsMatrix & L() const
Returns the address of the L factor associated with this factored matrix.
int NumMyCols() const
Returns the number of local matrix columns.
int Condest(bool Trans, double &ConditionNumberEstimate) const
Returns the maximum over all the condition number estimate for each local ILU set of factors.
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 ...
const Ifpack_IlukGraph & Graph() const
returns the address of the Ifpack_IlukGraph associated with this factored matrix.
Ifpack_IlukGraph: A class for constructing level filled graphs for use with ILU(k) class precondition...
virtual int LevelFill() const
Returns the level of fill used to construct this graph.
virtual const Epetra_BlockMap & DomainMap() const
Returns the Epetra_BlockMap object associated with the domain of this matrix operator.
virtual Epetra_CrsGraph & U_Graph()
Returns the graph of upper triangle of the ILU(k) graph as a Epetra_CrsGraph.
virtual Epetra_CrsGraph * OverlapGraph() const
Returns the the overlapped graph.
virtual int LevelOverlap() const
Returns the level of overlap used to construct this graph.
virtual Epetra_CrsGraph & L_Graph()
Returns the graph of lower triangle of the ILU(k) graph as a Epetra_CrsGraph.
virtual Epetra_Import * OverlapImporter() const
Returns the importer used to create the overlapped graph.