IFPACK Development
Loading...
Searching...
No Matches
Ifpack_CrsRick.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
43#include "Ifpack_CrsRick.h"
44#include "Epetra_Comm.h"
45#include "Epetra_Map.h"
46#include "Epetra_CrsGraph.h"
47#include "Epetra_CrsMatrix.h"
48#include "Epetra_Vector.h"
49#include "Epetra_MultiVector.h"
50
51#include <Teuchos_ParameterList.hpp>
52#include <ifp_parameters.h>
53
54//==============================================================================
56 : A_(A),
57 Graph_(Graph),
58 UseTranspose_(false),
59 Allocated_(false),
60 ValuesInitialized_(false),
61 Factored_(false),
62 RelaxValue_(0.0),
63 Condest_(-1.0),
64 Athresh_(0.0),
65 Rthresh_(1.0),
66 OverlapX_(0),
67 OverlapY_(0),
68 OverlapMode_(Zero)
69{
70 int ierr = Allocate();
71}
72
73//==============================================================================
75 : A_(FactoredMatrix.A_),
76 Graph_(FactoredMatrix.Graph_),
77 UseTranspose_(FactoredMatrix.UseTranspose_),
78 Allocated_(FactoredMatrix.Allocated_),
79 ValuesInitialized_(FactoredMatrix.ValuesInitialized_),
80 Factored_(FactoredMatrix.Factored_),
81 RelaxValue_(FactoredMatrix.RelaxValue_),
82 Condest_(FactoredMatrix.Condest_),
83 Athresh_(FactoredMatrix.Athresh_),
84 Rthresh_(FactoredMatrix.Rthresh_),
85 OverlapX_(0),
86 OverlapY_(0),
87 OverlapMode_(FactoredMatrix.OverlapMode_)
88{
89 U_ = new Epetra_CrsMatrix(FactoredMatrix.U());
90 D_ = new Epetra_Vector(Graph_.L_Graph().RowMap());
91
92}
93
94//==============================================================================
95int Ifpack_CrsRick::Allocate() {
96
97 // Allocate Epetra_CrsMatrix using ILUK graphs
98 U_ = new Epetra_CrsMatrix(Copy, Graph_.U_Graph());
99 D_ = new Epetra_Vector(Graph_.L_Graph().RowMap());
100
101
102 SetAllocated(true);
103 return(0);
104}
105//==============================================================================
107
108
109 delete U_;
110 delete D_; // Diagonal is stored separately. We store the inverse.
111
112 if (OverlapX_!=0) delete OverlapX_;
113 if (OverlapY_!=0) delete OverlapY_;
114
115 ValuesInitialized_ = false;
116 Factored_ = false;
117 Allocated_ = false;
118}
119
120//==========================================================================
121int Ifpack_CrsRick::SetParameters(const Teuchos::ParameterList& parameterlist,
122 bool cerr_warning_if_unused)
123{
125 params.double_params[Ifpack::relax_value] = RelaxValue_;
126 params.double_params[Ifpack::absolute_threshold] = Athresh_;
127 params.double_params[Ifpack::relative_threshold] = Rthresh_;
128 params.overlap_mode = OverlapMode_;
129
130 Ifpack::set_parameters(parameterlist, params, cerr_warning_if_unused);
131
132 RelaxValue_ = params.double_params[Ifpack::relax_value];
133 Athresh_ = params.double_params[Ifpack::absolute_threshold];
134 Rthresh_ = params.double_params[Ifpack::relative_threshold];
135 OverlapMode_ = params.overlap_mode;
136
137 return(0);
138}
139
140//==========================================================================
142
143 // if (!Allocated()) return(-1); // This test is not needed at this time. All constructors allocate.
144
145 int ierr = 0;
146 int i, j;
147 int * InI=0, * LI=0, * UI = 0;
148 double * InV=0, * LV=0, * UV = 0;
149 int NumIn, NumL, NumU;
150 bool DiagFound;
151 int NumNonzeroDiags = 0;
152
153 Epetra_CrsMatrix * OverlapA = (Epetra_CrsMatrix *) &A_;
154
155 if (Graph_.LevelOverlap()>0 && Graph_.L_Graph().DomainMap().DistributedGlobal()) {
156
157 OverlapA = new Epetra_CrsMatrix(Copy, *Graph_.OverlapGraph());
158 OverlapA->Import(A_, *Graph_.OverlapImporter(), Insert);
159 }
160
161 // Get Maximun Row length
162 int MaxNumEntries = OverlapA->MaxNumEntries();
163
164 InI = new int[MaxNumEntries]; // Allocate temp space
165 UI = new int[MaxNumEntries];
166 InV = new double[MaxNumEntries];
167 UV = new double[MaxNumEntries];
168
169 double *DV;
170 ierr = D_->ExtractView(&DV); // Get view of diagonal
171
172
173 // First we copy the user's matrix into diagonal vector and U, regardless of fill level
174
175 for (i=0; i< NumMyRows(); i++) {
176
177 OverlapA->ExtractMyRowCopy(i, MaxNumEntries, NumIn, InV, InI); // Get Values and Indices
178
179 // Split into L and U (we don't assume that indices are ordered).
180
181 NumL = 0;
182 NumU = 0;
183 DiagFound = false;
184
185 for (j=0; j< NumIn; j++) {
186 int k = InI[j];
187
188 if (k==i) {
189 DiagFound = true;
190 DV[i] += Rthresh_ * InV[j] + EPETRA_SGN(InV[j]) * Athresh_; // Store perturbed diagonal in Epetra_Vector D_
191 }
192
193 else if (k < 0) return(-1); // Out of range
194 else if (k<NumMyRows()) {
195 UI[NumU] = k;
196 UV[NumU] = InV[j];
197 NumU++;
198 }
199 }
200
201 // Check in things for this row of L and U
202
203 if (DiagFound) NumNonzeroDiags++;
204 if (NumU) U_->ReplaceMyValues(i, NumU, UV, UI);
205
206 }
207
208 delete [] UI;
209 delete [] UV;
210 delete [] InI;
211 delete [] InV;
212
213 if (Graph_.LevelOverlap()>0 && Graph_.L_Graph().DomainMap().DistributedGlobal()) delete OverlapA;
214
215
216 U_->FillComplete();
217
218 // At this point L and U have the values of A in the structure of L and U, and diagonal vector D
219
220 SetValuesInitialized(true);
221 SetFactored(false);
222
223 int TotalNonzeroDiags = 0;
224 Graph_.U_Graph().RowMap().Comm().SumAll(&NumNonzeroDiags, &TotalNonzeroDiags, 1);
225 if (Graph_.LevelOverlap()==0 &&
226 ((TotalNonzeroDiags!=NumGlobalRows()) ||
227 (TotalNonzeroDiags!=NumGlobalDiagonals()))) ierr = 1;
228 if (NumNonzeroDiags != NumMyDiagonals()) ierr = 1; // Diagonals are not right, warn user
229
230 return(ierr);
231}
232
233//==========================================================================
235
236 // if (!Allocated()) return(-1); // This test is not needed at this time. All constructors allocate.
237 if (!ValuesInitialized()) return(-2); // Must have values initialized.
238 if (Factored()) return(-3); // Can't have already computed factors.
239
240 SetValuesInitialized(false);
241
242 // MinMachNum should be officially defined, for now pick something a little
243 // bigger than IEEE underflow value
244
245 double MinDiagonalValue = Epetra_MinDouble;
246 double MaxDiagonalValue = 1.0/MinDiagonalValue;
247
248 int ierr = 0;
249 int i, j, k;
250 int * UI = 0;
251 double * UV = 0;
252 int NumIn, NumU;
253
254 // Get Maximun Row length
255 int MaxNumEntries = U_->MaxNumEntries() + 1;
256
257 int * InI = new int[MaxNumEntries]; // Allocate temp space
258 double * InV = new double[MaxNumEntries];
259 int * colflag = new int[NumMyCols()];
260
261 double *DV;
262 ierr = D_->ExtractView(&DV); // Get view of diagonal
263
264#ifdef IFPACK_FLOPCOUNTERS
265 int current_madds = 0; // We will count multiply-add as they happen
266#endif
267
268 // Now start the factorization.
269
270 // Need some integer workspace and pointers
271 int NumUU;
272 int * UUI;
273 double * UUV;
274 for (j=0; j<NumMyCols(); j++) colflag[j] = - 1;
275
276 for(i=0; i<NumMyRows(); i++) {
277
278 // Fill InV, InI with current row of L, D and U combined
279
280 NumIn = MaxNumEntries;
281 IFPACK_CHK_ERR(L_->ExtractMyRowCopy(i, NumIn, NumL, InV, InI)==0);
282 LV = InV;
283 LI = InI;
284
285 InV[NumL] = DV[i]; // Put in diagonal
286 InI[NumL] = i;
287
288 IFPACK_CHK_ERR(U_->ExtractMyRowCopy(i, NumIn-NumL-1, NumU, InV+NumL+1, InI+NumL+1));
289 NumIn = NumL+NumU+1;
290 UV = InV+NumL+1;
291 UI = InI+NumL+1;
292
293 // Set column flags
294 for (j=0; j<NumIn; j++) colflag[InI[j]] = j;
295
296 double diagmod = 0.0; // Off-diagonal accumulator
297
298 for (int jj=0; jj<NumL; jj++) {
299 j = InI[jj];
300 double multiplier = InV[jj]; // current_mults++;
301
302 InV[jj] *= DV[j];
303
304 IFPACK_CHK_ERR(U_->ExtractMyRowView(j, NumUU, UUV, UUI)); // View of row above
305
306 if (RelaxValue_==0.0) {
307 for (k=0; k<NumUU; k++) {
308 int kk = colflag[UUI[k]];
309 if (kk>-1) {
310 InV[kk] -= multiplier*UUV[k];
311#ifdef IFPACK_FLOPCOUNTERS
312 current_madds++;
313#endif
314#endif
315 }
316 }
317 }
318 else {
319 for (k=0; k<NumUU; k++) {
320 int kk = colflag[UUI[k]];
321 if (kk>-1) InV[kk] -= multiplier*UUV[k];
322 else diagmod -= multiplier*UUV[k];
323#ifdef IFPACK_FLOPCOUNTERS
324 current_madds++;
325#endif
326 }
327 }
328 }
329 if (NumL)
330 IFPACK_CHK_ERR(L_->ReplaceMyValues(i, NumL, LV, LI)); // Replace current row of L
331
332 DV[i] = InV[NumL]; // Extract Diagonal value
333
334 if (RelaxValue_!=0.0) {
335 DV[i] += RelaxValue_*diagmod; // Add off diagonal modifications
336 // current_madds++;
337 }
338
339 if (fabs(DV[i]) > MaxDiagonalValue) {
340 if (DV[i] < 0) DV[i] = - MinDiagonalValue;
341 else DV[i] = MinDiagonalValue;
342 }
343 else
344 DV[i] = 1.0/DV[i]; // Invert diagonal value
345
346 for (j=0; j<NumU; j++) UV[j] *= DV[i]; // Scale U by inverse of diagonal
347
348 if (NumU)
349 IFPACK_CHK_ERR(U_->ReplaceMyValues(i, NumU, UV, UI)); // Replace current row of L and U
350
351
352 // Reset column flags
353 for (j=0; j<NumIn; j++) colflag[InI[j]] = -1;
354 }
355
356
357#ifdef IFPACK_FLOPCOUNTERS
358 // Add up flops
359
360 double current_flops = 2 * current_madds;
361 double total_flops = 0;
362
363 Graph_.L_Graph().RowMap().Comm().SumAll(&current_flops, &total_flops, 1); // Get total madds across all PEs
364
365 // Now count the rest
366 total_flops += (double) L_->NumGlobalNonzeros(); // Accounts for multiplier above
367 total_flops += (double) D_->GlobalLength(); // Accounts for reciprocal of diagonal
368 if (RelaxValue_!=0.0) total_flops += 2 * (double)D_->GlobalLength(); // Accounts for relax update of diag
369
370 UpdateFlops(total_flops); // Update flop count
371#endif
372
373 delete [] InI;
374 delete [] InV;
375 delete [] colflag;
376
377 SetFactored(true);
378
379 return(ierr);
380
381}
382
383//=============================================================================
384int Ifpack_CrsRick::Solve(bool Trans, const Epetra_Vector& X,
385 Epetra_Vector& Y) const {
386//
387// This function finds Y such that LDU Y = X or U(trans) D L(trans) Y = X for a single RHS
388//
389
390 bool Upper = true;
391 bool Lower = false;
392 bool UnitDiagonal = true;
393
394 Epetra_Vector * X1 = (Epetra_Vector *) &X;
395 Epetra_Vector * Y1 = (Epetra_Vector *) &Y;
396
397 if (Graph_.LevelOverlap()>0 && Graph_.L_Graph().DomainMap().DistributedGlobal()) {
398 if (OverlapX_==0) { // Need to allocate space for overlap X and Y
399 OverlapX_ = new Epetra_Vector(Graph_.OverlapGraph()->RowMap());
400 OverlapY_ = new Epetra_Vector(Graph_.OverlapGraph()->RowMap());
401 }
402 OverlapX_->Import(X,*Graph_.OverlapImporter(), Insert); // Import X values for solve
403 X1 = (Epetra_Vector *) OverlapX_;
404 Y1 = (Epetra_Vector *) OverlapY_; // Set pointers for X1 and Y1 to point to overlap space
405 }
406
407#ifdef IFPACK_FLOPCOUNTERS
408 Epetra_Flops * counter = this->GetFlopCounter();
409 if (counter!=0) {
410 L_->SetFlopCounter(*counter);
411 Y1->SetFlopCounter(*counter);
412 U_->SetFlopCounter(*counter);
413 }
414#endif
415
416 if (!Trans) {
417
418 L_->Solve(Lower, Trans, UnitDiagonal, *X1, *Y1);
419 Y1->Multiply(1.0, *D_, *Y1, 0.0); // y = D*y (D_ has inverse of diagonal)
420 U_->Solve(Upper, Trans, UnitDiagonal, *Y1, *Y1); // Solve Uy = y
421 }
422 else
423 {
424 U_->Solve(Upper, Trans, UnitDiagonal, *X1, *Y1); // Solve Uy = y
425 Y1->Multiply(1.0, *D_, *Y1, 0.0); // y = D*y (D_ has inverse of diagonal)
426 L_->Solve(Lower, Trans, UnitDiagonal, *Y1, *Y1);
427
428 }
429
430 // Export computed Y values as directed
431 if (Graph_.LevelOverlap()>0 && Graph_.L_Graph().DomainMap().DistributedGlobal())
432 Y.Export(*OverlapY_,*Graph_.OverlapImporter(), OverlapMode_);
433 return(0);
434}
435
436
437//=============================================================================
439 Epetra_MultiVector& Y) const {
440//
441// This function finds Y such that LDU Y = X or U(trans) D L(trans) Y = X for multiple RHS
442//
443
444 if (X.NumVectors()!=Y.NumVectors()) EPETRA_CHK_ERR(-1); // Return error: X and Y not the same size
445
446 bool Upper = true;
447 bool Lower = false;
448 bool UnitDiagonal = true;
449
452
453 if (Graph_.LevelOverlap()>0 && Graph_.L_Graph().DomainMap().DistributedGlobal()) {
454 // Make sure the number of vectors in the multivector is the same as before.
455 if (OverlapX_!=0) {
456 if (OverlapX_->NumVectors()!=X.NumVectors()) {
457 delete OverlapX_; OverlapX_ = 0;
458 delete OverlapY_; OverlapY_ = 0;
459 }
460 }
461 if (OverlapX_==0) { // Need to allocate space for overlap X and Y
462 OverlapX_ = new Epetra_MultiVector(Graph_.OverlapGraph()->RowMap(), X.NumVectors());
463 OverlapY_ = new Epetra_MultiVector(Graph_.OverlapGraph()->RowMap(), Y.NumVectors());
464 }
465 OverlapX_->Import(X,*Graph_.OverlapImporter(), Insert); // Import X values for solve
466 X1 = OverlapX_;
467 Y1 = OverlapY_; // Set pointers for X1 and Y1 to point to overlap space
468 }
469
470#ifdef IFPACK_FLOPCOUNTERS
471 Epetra_Flops * counter = this->GetFlopCounter();
472 if (counter!=0) {
473 L_->SetFlopCounter(*counter);
474 Y1->SetFlopCounter(*counter);
475 U_->SetFlopCounter(*counter);
476 }
477#endif
478
479 if (!Trans) {
480
481 L_->Solve(Lower, Trans, UnitDiagonal, *X1, *Y1);
482 Y1->Multiply(1.0, *D_, *Y1, 0.0); // y = D*y (D_ has inverse of diagonal)
483 U_->Solve(Upper, Trans, UnitDiagonal, *Y1, *Y1); // Solve Uy = y
484 }
485 else
486 {
487 U_->Solve(Upper, Trans, UnitDiagonal, *X1, *Y1); // Solve Uy = y
488 Y1->Multiply(1.0, *D_, *Y1, 0.0); // y = D*y (D_ has inverse of diagonal)
489 L_->Solve(Lower, Trans, UnitDiagonal, *Y1, *Y1);
490
491 }
492
493 // Export computed Y values as directed
494 if (Graph_.LevelOverlap()>0 && Graph_.L_Graph().DomainMap().DistributedGlobal())
495 Y.Export(*OverlapY_,*Graph_.OverlapImporter(), OverlapMode_);
496 return(0);
497}
498//=============================================================================
500 Epetra_MultiVector& Y) const {
501//
502// This function finds X such that LDU Y = X or U(trans) D L(trans) Y = X for multiple RHS
503//
504
505 if (X.NumVectors()!=Y.NumVectors()) EPETRA_CHK_ERR(-1); // Return error: X and Y not the same size
506
507 bool Upper = true;
508 bool Lower = false;
509 bool UnitDiagonal = true;
510
513
514 if (Graph_.LevelOverlap()>0 && Graph_.L_Graph().DomainMap().DistributedGlobal()) {
515 // Make sure the number of vectors in the multivector is the same as before.
516 if (OverlapX_!=0) {
517 if (OverlapX_->NumVectors()!=X.NumVectors()) {
518 delete OverlapX_; OverlapX_ = 0;
519 delete OverlapY_; OverlapY_ = 0;
520 }
521 }
522 if (OverlapX_==0) { // Need to allocate space for overlap X and Y
523 OverlapX_ = new Epetra_MultiVector(Graph_.OverlapGraph()->RowMap(), X.NumVectors());
524 OverlapY_ = new Epetra_MultiVector(Graph_.OverlapGraph()->RowMap(), Y.NumVectors());
525 }
526 OverlapX_->Import(X,*Graph_.OverlapImporter(), Insert); // Import X values for solve
527 X1 = OverlapX_;
528 Y1 = OverlapY_; // Set pointers for X1 and Y1 to point to overlap space
529 }
530
531#ifdef IFPACK_FLOPCOUNTERS
532 Epetra_Flops * counter = this->GetFlopCounter();
533 if (counter!=0) {
534 L_->SetFlopCounter(*counter);
535 Y1->SetFlopCounter(*counter);
536 U_->SetFlopCounter(*counter);
537 }
538#endif
539
540 if (Trans) {
541
542 L_->Multiply(Trans, *X1, *Y1);
543 Y1->Update(1.0, *X1, 1.0); // Y1 = Y1 + X1 (account for implicit unit diagonal)
544 Y1->ReciprocalMultiply(1.0, *D_, *Y1, 0.0); // y = D*y (D_ has inverse of diagonal)
545 Epetra_MultiVector Y1temp(*Y1); // Need a temp copy of Y1
546 U_->Multiply(Trans, Y1temp, *Y1);
547 Y1->Update(1.0, Y1temp, 1.0); // (account for implicit unit diagonal)
548 }
549 else {
550 U_->Multiply(Trans, *X1, *Y1); //
551 Y1->Update(1.0, *X1, 1.0); // Y1 = Y1 + X1 (account for implicit unit diagonal)
552 Y1->ReciprocalMultiply(1.0, *D_, *Y1, 0.0); // y = D*y (D_ has inverse of diagonal)
553 Epetra_MultiVector Y1temp(*Y1); // Need a temp copy of Y1
554 L_->Multiply(Trans, Y1temp, *Y1);
555 Y1->Update(1.0, Y1temp, 1.0); // (account for implicit unit diagonal)
556 }
557
558 // Export computed Y values as directed
559 if (Graph_.LevelOverlap()>0 && Graph_.L_Graph().DomainMap().DistributedGlobal())
560 Y.Export(*OverlapY_,*Graph_.OverlapImporter(), OverlapMode_);
561 return(0);
562}
563//=============================================================================
564int Ifpack_CrsRick::Condest(bool Trans, double & ConditionNumberEstimate) const {
565
566 if (Condest_>=0.0) {
567 ConditionNumberEstimate = Condest_;
568 return(0);
569 }
570 // Create a vector with all values equal to one
571 Epetra_Vector Ones(A_.RowMap());
572 Epetra_Vector OnesResult(Ones);
573 Ones.PutScalar(1.0);
574
575 EPETRA_CHK_ERR(Solve(Trans, Ones, OnesResult)); // Compute the effect of the solve on the vector of ones
576 EPETRA_CHK_ERR(OnesResult.Abs(OnesResult)); // Make all values non-negative
577 EPETRA_CHK_ERR(OnesResult.MaxValue(&ConditionNumberEstimate)); // Get the maximum value across all processors
578 Condest_ = ConditionNumberEstimate; // Save value for possible later calls
579 return(0);
580}
581//=============================================================================
582// Non-member functions
583
584std::ostream& operator << (std::ostream& os, const Ifpack_CrsRick& A)
585{
586/* Epetra_fmtflags olda = os.setf(ios::right,ios::adjustfield);
587 Epetra_fmtflags oldf = os.setf(ios::scientific,ios::floatfield);
588 int oldp = os.precision(12); */
589 int LevelFill = A.Graph().LevelFill();
590 int LevelOverlap = A.Graph().LevelOverlap();
591 Epetra_CrsMatrix & L = (Epetra_CrsMatrix &) A.L();
592 Epetra_CrsMatrix & U = (Epetra_CrsMatrix &) A.U();
593 Epetra_Vector & D = (Epetra_Vector &) A.D();
594
595 os.width(14);
596 os << endl;
597 os << " Level of Fill = "; os << LevelFill;
598 os << endl;
599 os.width(14);
600 os << " Level of Overlap = "; os << LevelOverlap;
601 os << endl;
602
603 os.width(14);
604 os << " Lower Triangle = ";
605 os << endl;
606 os << L; // Let Epetra_CrsMatrix handle the rest.
607 os << endl;
608
609 os.width(14);
610 os << " Inverse of Diagonal = ";
611 os << endl;
612 os << D; // Let Epetra_Vector handle the rest.
613 os << endl;
614
615 os.width(14);
616 os << " Upper Triangle = ";
617 os << endl;
618 os << U; // Let Epetra_CrsMatrix handle the rest.
619 os << endl;
620
621 // Reset os flags
622
623/* os.setf(olda,ios::adjustfield);
624 os.setf(oldf,ios::floatfield);
625 os.precision(oldp); */
626
627 return os;
628}
Insert
Zero
Copy
bool DistributedGlobal() const
const Epetra_Comm & Comm() const
virtual int SumAll(double *PartialSums, double *GlobalSums, int Count) const=0
Epetra_Flops * GetFlopCounter() const
void SetFlopCounter(const Epetra_Flops &FlopCounter_in)
void UpdateFlops(int Flops_in) const
const Epetra_BlockMap & DomainMap() const
const Epetra_BlockMap & RowMap() const
const Epetra_Map & RowMap() const
int ExtractMyRowView(int MyRow, int &NumEntries, double *&Values, int *&Indices) const
int FillComplete(bool OptimizeDataStorage=true)
int Solve(bool Upper, bool Trans, bool UnitDiagonal, const Epetra_Vector &x, Epetra_Vector &y) const
int ExtractMyRowCopy(int MyRow, int Length, int &NumEntries, double *Values, int *Indices) const
int MaxNumEntries() const
int ReplaceMyValues(int MyRow, int NumEntries, const double *Values, const int *Indices)
int Multiply(bool TransA, const Epetra_Vector &x, Epetra_Vector &y) const
int Import(const Epetra_SrcDistObject &A, const Epetra_Import &Importer, Epetra_CombineMode CombineMode, const Epetra_OffsetIndex *Indexor=0)
int Export(const Epetra_SrcDistObject &A, const Epetra_Import &Importer, Epetra_CombineMode CombineMode, const Epetra_OffsetIndex *Indexor=0)
int GlobalLength() const
int NumVectors() const
int Abs(const Epetra_MultiVector &A)
int Multiply(char TransA, char TransB, double ScalarAB, const Epetra_MultiVector &A, const Epetra_MultiVector &B, double ScalarThis)
int Update(double ScalarA, const Epetra_MultiVector &A, double ScalarThis)
int MaxValue(double *Result) const
int PutScalar(double ScalarConstant)
int ReciprocalMultiply(double ScalarAB, const Epetra_MultiVector &A, const Epetra_MultiVector &B, double ScalarThis)
int ExtractView(double **V) const
Ifpack_CrsRick: A class for constructing and using an incomplete Cholesky (IC) factorization of a giv...
int NumMyCols() const
Returns the number of local matrix columns.
Ifpack_CrsRick(const Epetra_CrsMatrix &A, const Ifpack_IlukGraph &Graph)
Ifpack_CrsRick constuctor with variable number of indices per row.
virtual int NumGlobalDiagonals() const
Returns the number of diagonal entries found in the global input graph.
virtual int NumMyDiagonals() const
Returns the number of diagonal entries found in the local input graph.
int SetParameters(const Teuchos::ParameterList &parameterlist, bool cerr_warning_if_unused=false)
Set parameters using a Teuchos::ParameterList object.
int Factor()
Compute IC factor L using the specified graph, diagonal perturbation thresholds and relaxation parame...
const Epetra_Vector & D() const
Returns the address of the D factor associated with this factored matrix.
int NumMyRows() const
Returns the number of local matrix rows.
int InitValues()
Initialize L and U with values from user matrix A.
int NumGlobalRows() const
Returns the number of global matrix rows.
const Ifpack_IlukGraph & Graph() const
Returns the address of the Ifpack_IlukGraph associated with this factored matrix.
virtual ~Ifpack_CrsRick()
Ifpack_CrsRick Destructor.
bool Factored() const
If factor is completed, this query returns true, otherwise it returns false.
const Epetra_CrsMatrix & U() const
Returns the address of the U factor associated with this factored matrix.
int Solve(bool Trans, const Epetra_Vector &x, Epetra_Vector &y) const
Returns the result of a Ifpack_CrsRick forward/back solve on a Epetra_Vector x in y.
int Condest(bool Trans, double &ConditionNumberEstimate) const
Returns the maximum over all the condition number estimate for each local IC set of factors.
int Multiply(bool Trans, const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Returns the result of multiplying U, D and U^T in that order on an Epetra_MultiVector X in Y.
bool ValuesInitialized() const
If values have been initialized, this query returns true, otherwise it returns false.
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 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.