IFPACK Development
Loading...
Searching...
No Matches
Ifpack_TriDiContainer.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_ConfigDefs.h"
44#include "Ifpack_TriDiContainer.h"
45#include "Epetra_RowMatrix.h"
46
47//==============================================================================
49{
50 return(NumRows_);
51}
52
53//==============================================================================
55{
56
57 IsInitialized_ = false;
58
59 IFPACK_CHK_ERR(LHS_.Reshape(NumRows_,NumVectors_));
60 IFPACK_CHK_ERR(RHS_.Reshape(NumRows_,NumVectors_));
61 IFPACK_CHK_ERR(ID_.Reshape(NumRows_,NumVectors_));
62 IFPACK_CHK_ERR(Matrix_.Reshape(NumRows_,NumRows_));
63
64 // zero out matrix elements
65 // for (int i = 0 ; i < NumRows_ ; ++i)
66 // for (int j = 0 ; j < NumRows_ ; ++j)
67 // Matrix_(i,j) = 0.0;
68
69 if(Matrix_.A() == NULL ) Matrix_.Shape(NumRows_);
70 int size = (NumRows_ == 1) ? 1 : 4*(NumRows_ -1);
71 memset(Matrix_.A(),0,sizeof(double)*size);
72
73 // zero out vector elements
74 for (int i = 0 ; i < NumRows_ ; ++i)
75 for (int j = 0 ; j < NumVectors_ ; ++j) {
76 LHS_(i,j) = 0.0;
77 RHS_(i,j) = 0.0;
78 }
79
80 // Set to -1 ID_'s
81 for (int i = 0 ; i < NumRows_ ; ++i)
82 ID_(i) = -1;
83
84 if (NumRows_ != 0) {
85 IFPACK_CHK_ERR(Solver_.SetMatrix(Matrix_));
86 IFPACK_CHK_ERR(Solver_.SetVectors(LHS_,RHS_));
87 }
88
89 IsInitialized_ = true;
90 return(0);
91
92}
93
94//==============================================================================
95double& Ifpack_TriDiContainer::LHS(const int i, const int Vector)
96{
97 return(LHS_.A()[Vector * NumRows_ + i]);
98}
99
100//==============================================================================
101double& Ifpack_TriDiContainer::RHS(const int i, const int Vector)
102{
103 return(RHS_.A()[Vector * NumRows_ + i]);
104}
105
106//==============================================================================
108SetMatrixElement(const int row, const int col, const double value)
109{
110 if (IsInitialized() == false)
111 IFPACK_CHK_ERR(Initialize());
112
113 if ((row < 0) || (row >= NumRows())) {
114 IFPACK_CHK_ERR(-2); // not in range
115 }
116
117 if ((col < 0) || (col >= NumRows())) {
118 IFPACK_CHK_ERR(-2); // not in range
119 }
120
121 Matrix_(row, col) = value;
122
123 return(0);
124
125}
126
127//==============================================================================
129{
130
131 if (!IsComputed()) {
132 IFPACK_CHK_ERR(-1);
133 }
134
135 if (NumRows_ != 0)
136 IFPACK_CHK_ERR(Solver_.Solve());
137
138#ifdef IFPACK_FLOPCOUNTERS
139 ApplyInverseFlops_ += 2.0 * NumVectors_ * NumRows_ * NumRows_;
140#endif
141 return(0);
142}
143
144//==============================================================================
146{
147 return(ID_[i]);
148}
149
150//==============================================================================
151// FIXME: optimize performances of this guy...
152int Ifpack_TriDiContainer::Extract(const Epetra_RowMatrix& Matrix_in)
153{
154
155 for (int j = 0 ; j < NumRows_ ; ++j) {
156 // be sure that the user has set all the ID's
157 if (ID(j) == -1)
158 IFPACK_CHK_ERR(-2);
159 // be sure that all are local indices
160 if (ID(j) > Matrix_in.NumMyRows())
161 IFPACK_CHK_ERR(-2);
162 }
163
164 // allocate storage to extract matrix rows.
165 int Length = Matrix_in.MaxNumEntries();
166 std::vector<double> Values;
167 Values.resize(Length);
168 std::vector<int> Indices;
169 Indices.resize(Length);
170
171 for (int j = 0 ; j < NumRows_ ; ++j) {
172
173 int LRID = ID(j);
174
175 int NumEntries;
176
177 int ierr =
178 Matrix_in.ExtractMyRowCopy(LRID, Length, NumEntries,
179 &Values[0], &Indices[0]);
180 IFPACK_CHK_ERR(ierr);
181
182 for (int k = 0 ; k < NumEntries ; ++k) {
183
184 int LCID = Indices[k];
185
186 // skip off-processor elements
187 if (LCID >= Matrix_in.NumMyRows())
188 continue;
189
190 // for local column IDs, look for each ID in the list
191 // of columns hosted by this object
192 // FIXME: use STL
193 int jj = -1;
194 for (int kk = 0 ; kk < NumRows_ ; ++kk)
195 if (ID(kk) == LCID)
196 jj = kk;
197
198 if (jj != -1)
199 SetMatrixElement(j,jj,Values[k]);
200
201 }
202 }
203
204 return(0);
205}
206
207//==============================================================================
209{
210 IsComputed_ = false;
211 if (IsInitialized() == false) {
212 IFPACK_CHK_ERR(Initialize());
213 }
214
215 if (KeepNonFactoredMatrix_)
216 NonFactoredMatrix_ = Matrix_;
217
218 // extract local rows and columns
219 IFPACK_CHK_ERR(Extract(Matrix_in));
220
221 if (KeepNonFactoredMatrix_)
222 NonFactoredMatrix_ = Matrix_;
223
224 // factorize the matrix using LAPACK
225 if (NumRows_ != 0)
226 IFPACK_CHK_ERR(Solver_.Factor());
227
228 Label_ = "Ifpack_TriDiContainer";
229
230 // not sure of count
231#ifdef IFPACK_FLOPCOUNTERS
232 ComputeFlops_ += 4.0 * NumRows_ * NumRows_ * NumRows_ / 3;
233#endif
234 IsComputed_ = true;
235
236 return(0);
237}
238
239//==============================================================================
241 {
242 IFPACK_CHK_ERR(-300);
243 // if (IsComputed() == false)
244 // IFPACK_CHK_ERR(-3);
245
246 // if (KeepNonFactoredMatrix_) {
247 // IFPACK_CHK_ERR(RHS_.Multiply('N','N', 1.0,NonFactoredMatrix_,LHS_,0.0));
248 // }
249 // else
250 // IFPACK_CHK_ERR(RHS_.Multiply('N','N', 1.0,Matrix_,LHS_,0.0));
251
252 // #ifdef IFPACK_FLOPCOUNTERS
253 // ApplyFlops_ += 2 * NumRows_ * NumRows_;
254 // #endif
255 //return(0); // unreachable
256 }
257
258//==============================================================================
259std::ostream& Ifpack_TriDiContainer::Print(std::ostream & os) const
260{
261 using std::endl;
262
263 os << "================================================================================" << endl;
264 os << "Ifpack_TriDiContainer" << endl;
265 os << "Number of rows = " << NumRows() << endl;
266 os << "Number of vectors = " << NumVectors() << endl;
267 os << "IsInitialized() = " << IsInitialized() << endl;
268 os << "IsComputed() = " << IsComputed() << endl;
269#ifdef IFPACK_FLOPCOUNTERS
270 os << "Flops in Compute() = " << ComputeFlops() << endl;
271 os << "Flops in ApplyInverse() = " << ApplyInverseFlops() << endl;
272#endif
273 os << "================================================================================" << endl;
274 os << endl;
275
276 return(os);
277}
int Reshape(int NumRows, int NumCols)
virtual int NumMyRows() const=0
virtual int MaxNumEntries() const=0
virtual int ExtractMyRowCopy(int MyRow, int Length, int &NumEntries, double *Values, int *Indices) const=0
double * A() const
int Reshape(int NumRows, int NumCols)
int Shape(int NumRowCol)
Set dimensions of a Ifpack_SerialTriDiMatrix object; init values to zero.
double * A() const
Returns pointer to the this matrix.
int SetMatrix(Ifpack_SerialTriDiMatrix &A)
Sets the pointers for coefficient matrix.
virtual int Solve(void)
Computes the solution X to AX = B for the this matrix and the B provided to SetVectors()....
virtual int Factor(void)
Computes the in-place LU factorization of the matrix using the LAPACK routine DGETRF.
int SetVectors(Epetra_SerialDenseMatrix &X, Epetra_SerialDenseMatrix &B)
Sets the pointers for left and right hand side vector(s).
virtual double ComputeFlops() const
Returns the flops in Compute().
virtual bool IsInitialized() const
Returns true is the container has been successfully initialized.
virtual int Compute(const Epetra_RowMatrix &Matrix_in)
Finalizes the linear system matrix and prepares for the application of the inverse.
virtual int ApplyInverse()
Apply the inverse of the matrix to RHS, results are stored in LHS.
virtual const Epetra_SerialDenseMatrix & RHS() const
Returns the dense vector containing the RHS.
virtual int Initialize()
Initialize the container.
virtual int Apply()
Apply the matrix to RHS, results are stored in LHS.
virtual double ApplyInverseFlops() const
Returns the flops in ApplyInverse().
virtual int SetMatrixElement(const int row, const int col, const double value)
Set the matrix element (row,col) to value.
virtual const Epetra_IntSerialDenseVector & ID() const
Returns the integer dense vector of IDs.
virtual int NumRows() const
Returns the number of rows of the matrix and LHS/RHS.
virtual bool IsComputed() const
Returns true is the container has been successfully computed.
virtual const Epetra_SerialDenseMatrix & LHS() const
Returns the dense vector containing the LHS.
virtual std::ostream & Print(std::ostream &os) const
Prints basic information on iostream. This function is used by operator<<.
virtual int NumVectors() const
Returns the number of vectors in LHS/RHS.