IFPACK  Development
 All Classes Namespaces Files Functions Variables Enumerations Friends Pages
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 //=========================================================================
151 int 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 //=============================================================================
165 int 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 }
const Epetra_Map & RangeMap() const
int Multiply(bool TransA, const Epetra_Vector &x, Epetra_Vector &y) const
Ifpack_OverlapSolveObject: Provides Overlapped Forward/back solve services for Ifpack.
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 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...
virtual ~Ifpack_OverlapSolveObject()
Ifpack_OverlapSolveObject Destructor.
int Solve(bool Upper, bool Trans, bool UnitDiagonal, const Epetra_Vector &x, Epetra_Vector &y) const
const Epetra_Import * Importer() const
Ifpack_OverlapSolveObject(char *Label, const Epetra_Comm &Comm)
Constructor.
const Epetra_Export * Exporter() const
int Condest(bool Trans, double &ConditionNumberEstimate) const
Returns the maximum over all the condition number estimate for each local ILU set of factors...
const Epetra_Map & DomainMap() const