IFPACK Development
Loading...
Searching...
No Matches
Ifpack_OverlapSolveObject.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_OverlapSolveObject.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#include "Epetra_Flops.h"
51
52//==============================================================================
54 : Label_(Label),
55 L_(0),
56 UseLTrans_(false),
57 D_(0),
58 U_(0),
59 UseUTrans_(false),
60 UseTranspose_(false),
61 Comm_(Comm),
62 Condest_(-1.0),
63 OverlapMode_(Zero)
64{
65}
66
67//==============================================================================
69 : Label_(Source.Label_),
70 L_(Source.L_),
71 UseLTrans_(Source.UseLTrans_),
72 D_(Source.D_),
73 U_(Source.U_),
74 UseUTrans_(Source.UseUTrans_),
75 UseTranspose_(Source.UseTranspose_),
76 Comm_(Source.Comm_),
77 Condest_(Source.Condest_),
78 OverlapMode_(Source.OverlapMode_)
79{
80}
81//==============================================================================
83}
84//==============================================================================
85//=============================================================================
87 Epetra_MultiVector& Y) const {
88//
89// This function finds Y such that LDU Y = X or U(trans) D L(trans) Y = X for multiple RHS
90//
91
92 // First generate X and Y as needed for this function
93 Epetra_MultiVector * X1 = 0;
94 Epetra_MultiVector * Y1 = 0;
95 EPETRA_CHK_ERR(SetupXY(Trans, X, Y, X1, Y1));
96
97 bool Upper = true;
98 bool Lower = false;
99 bool UnitDiagonal = true;
100
101 if (!Trans) {
102
103 EPETRA_CHK_ERR(L_->Solve(Lower, Trans, UnitDiagonal, *X1, *Y1));
104 EPETRA_CHK_ERR(Y1->Multiply(1.0, *D_, *Y1, 0.0)); // y = D*y (D_ has inverse of diagonal)
105 EPETRA_CHK_ERR(U_->Solve(Upper, Trans, UnitDiagonal, *Y1, *Y1)); // Solve Uy = y
106 if (L_->Exporter()!=0) {EPETRA_CHK_ERR(Y.Export(*Y1,*L_->Exporter(), OverlapMode_));} // Export computed Y values if needed
107 }
108 else {
109 EPETRA_CHK_ERR(U_->Solve(Upper, Trans, UnitDiagonal, *X1, *Y1)); // Solve Uy = y
110 EPETRA_CHK_ERR(Y1->Multiply(1.0, *D_, *Y1, 0.0)); // y = D*y (D_ has inverse of diagonal)
111 EPETRA_CHK_ERR(L_->Solve(Lower, Trans, UnitDiagonal, *Y1, *Y1));
112 if (U_->Importer()!=0) {EPETRA_CHK_ERR(Y.Export(*Y1,*U_->Importer(), OverlapMode_));} // Export computed Y values if needed
113 }
114
115 return(0);
116}
117//=============================================================================
119 Epetra_MultiVector& Y) const {
120//
121// This function finds X such that LDU Y = X or U(trans) D L(trans) Y = X for multiple RHS
122//
123
124 // First generate X and Y as needed for this function
125 Epetra_MultiVector * X1 = 0;
126 Epetra_MultiVector * Y1 = 0;
127 EPETRA_CHK_ERR(SetupXY(Trans, X, Y, X1, Y1));
128
129 if (!Trans) {
130 EPETRA_CHK_ERR(U_->Multiply(Trans, *X1, *Y1)); //
131 EPETRA_CHK_ERR(Y1->Update(1.0, *X1, 1.0)); // Y1 = Y1 + X1 (account for implicit unit diagonal)
132 EPETRA_CHK_ERR(Y1->ReciprocalMultiply(1.0, *D_, *Y1, 0.0)); // y = D*y (D_ has inverse of diagonal)
133 Epetra_MultiVector Y1temp(*Y1); // Need a temp copy of Y1
134 EPETRA_CHK_ERR(L_->Multiply(Trans, Y1temp, *Y1));
135 EPETRA_CHK_ERR(Y1->Update(1.0, Y1temp, 1.0)); // (account for implicit unit diagonal)
136 if (L_->Exporter()!=0) {EPETRA_CHK_ERR(Y.Export(*Y1,*L_->Exporter(), OverlapMode_));} // Export computed Y values if needed
137 }
138 else {
139
140 EPETRA_CHK_ERR(L_->Multiply(Trans, *X1, *Y1));
141 EPETRA_CHK_ERR(Y1->Update(1.0, *X1, 1.0)); // Y1 = Y1 + X1 (account for implicit unit diagonal)
142 EPETRA_CHK_ERR(Y1->ReciprocalMultiply(1.0, *D_, *Y1, 0.0)); // y = D*y (D_ has inverse of diagonal)
143 Epetra_MultiVector Y1temp(*Y1); // Need a temp copy of Y1
144 EPETRA_CHK_ERR(U_->Multiply(Trans, Y1temp, *Y1));
145 EPETRA_CHK_ERR(Y1->Update(1.0, Y1temp, 1.0)); // (account for implicit unit diagonal)
146 if (L_->Exporter()!=0) {EPETRA_CHK_ERR(Y.Export(*Y1,*L_->Exporter(), OverlapMode_));}
147 }
148 return(0);
149}
150//=========================================================================
151int Ifpack_OverlapSolveObject::SetupXY(bool Trans,
152 const Epetra_MultiVector& Xin, const Epetra_MultiVector& Yin,
153 Epetra_MultiVector * & Xout, Epetra_MultiVector * & Yout) const {
154
155 // Generate an X and Y suitable for performing Solve() and Multiply() methods
156
157 if (Xin.NumVectors()!=Yin.NumVectors()) EPETRA_CHK_ERR(-1); // Return error: X and Y not the same size
158
159 //cout << "Xin = " << Xin << endl;
160 Xout = (Epetra_MultiVector *) &Xin;
161 Yout = (Epetra_MultiVector *) &Yin;
162 return(0);
163}
164//=============================================================================
165int Ifpack_OverlapSolveObject::Condest(bool Trans, double & ConditionNumberEstimate) const {
166
167 if (Condest_>=0.0) {
168 ConditionNumberEstimate = Condest_;
169 return(0);
170 }
171 // Create a vector with all values equal to one
172 Epetra_Vector Ones(U_->DomainMap());
173 Epetra_Vector OnesResult(L_->RangeMap());
174 Ones.PutScalar(1.0);
175
176 EPETRA_CHK_ERR(Solve(Trans, Ones, OnesResult)); // Compute the effect of the solve on the vector of ones
177 EPETRA_CHK_ERR(OnesResult.Abs(OnesResult)); // Make all values non-negative
178 EPETRA_CHK_ERR(OnesResult.MaxValue(&ConditionNumberEstimate)); // Get the maximum value across all processors
179 Condest_ = ConditionNumberEstimate; // Save value for possible later calls
180 return(0);
181}
Zero
const Epetra_Export * Exporter() const
int Solve(bool Upper, bool Trans, bool UnitDiagonal, const Epetra_Vector &x, Epetra_Vector &y) const
const Epetra_Map & RangeMap() const
const Epetra_Map & DomainMap() const
const Epetra_Import * Importer() const
int Multiply(bool TransA, const Epetra_Vector &x, Epetra_Vector &y) 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 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)
Ifpack_OverlapSolveObject: Provides Overlapped Forward/back solve services for Ifpack.
Ifpack_OverlapSolveObject(char *Label, const Epetra_Comm &Comm)
Constructor.
int Solve(bool Trans, const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Returns the result of a Ifpack_CrsIlut forward/back solve on a Epetra_MultiVector X in Y (works for E...
int Condest(bool Trans, double &ConditionNumberEstimate) const
Returns the maximum over all the condition number estimate for each local ILU set of factors.
virtual ~Ifpack_OverlapSolveObject()
Ifpack_OverlapSolveObject 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.