Ifpack Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
Ifpack_DenseContainer.cpp
Go to the documentation of this file.
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"
45#include "Epetra_RowMatrix.h"
46
47//==============================================================================
49{
50 return(NumRows_);
51}
52
53//==============================================================================
55{
56
57 IsInitialized_ = false;
58
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 // zero out vector elements
70 for (int i = 0 ; i < NumRows_ ; ++i)
71 for (int j = 0 ; j < NumVectors_ ; ++j) {
72 LHS_(i,j) = 0.0;
73 RHS_(i,j) = 0.0;
74 }
75
76 // Set to -1 ID_'s
77 for (int i = 0 ; i < NumRows_ ; ++i)
78 ID_(i) = -1;
79
80 if (NumRows_ != 0) {
83 }
84
85 IsInitialized_ = true;
86 return(0);
87
88}
89
90//==============================================================================
91double& Ifpack_DenseContainer::LHS(const int i, const int Vector)
92{
93 return(LHS_.A()[Vector * NumRows_ + i]);
94}
95
96//==============================================================================
97double& Ifpack_DenseContainer::RHS(const int i, const int Vector)
98{
99 return(RHS_.A()[Vector * NumRows_ + i]);
100}
101
102//==============================================================================
104SetMatrixElement(const int row, const int col, const double value)
105{
106 if (IsInitialized() == false)
108
109 if ((row < 0) || (row >= NumRows())) {
110 IFPACK_CHK_ERR(-2); // not in range
111 }
112
113 if ((col < 0) || (col >= NumRows())) {
114 IFPACK_CHK_ERR(-2); // not in range
115 }
116
117 Matrix_(row, col) = value;
118
119 return(0);
120
121}
122
123//==============================================================================
125{
126
127 if (!IsComputed()) {
128 IFPACK_CHK_ERR(-1);
129 }
130
131 if (NumRows_ != 0)
133
134#ifdef IFPACK_FLOPCOUNTERS
136#endif
137 return(0);
138}
139
140//==============================================================================
142{
143 return(ID_[i]);
144}
145
146//==============================================================================
147// FIXME: optimize performances of this guy...
149{
150
151 for (int j = 0 ; j < NumRows_ ; ++j) {
152 // be sure that the user has set all the ID's
153 if (ID(j) == -1)
154 IFPACK_CHK_ERR(-2);
155 // be sure that all are local indices
156 if (ID(j) > Matrix_in.NumMyRows())
157 IFPACK_CHK_ERR(-2);
158 }
159
160 // allocate storage to extract matrix rows.
161 int Length = Matrix_in.MaxNumEntries();
162 std::vector<double> Values;
163 Values.resize(Length);
164 std::vector<int> Indices;
165 Indices.resize(Length);
166
167 for (int j = 0 ; j < NumRows_ ; ++j) {
168
169 int LRID = ID(j);
170
171 int NumEntries;
172
173 int ierr =
174 Matrix_in.ExtractMyRowCopy(LRID, Length, NumEntries,
175 &Values[0], &Indices[0]);
176 IFPACK_CHK_ERR(ierr);
177
178 for (int k = 0 ; k < NumEntries ; ++k) {
179
180 int LCID = Indices[k];
181
182 // skip off-processor elements
183 if (LCID >= Matrix_in.NumMyRows())
184 continue;
185
186 // for local column IDs, look for each ID in the list
187 // of columns hosted by this object
188 // FIXME: use STL
189 int jj = -1;
190 for (int kk = 0 ; kk < NumRows_ ; ++kk)
191 if (ID(kk) == LCID)
192 jj = kk;
193
194 if (jj != -1)
195 SetMatrixElement(j,jj,Values[k]);
196
197 }
198 }
199
200 return(0);
201}
202
203//==============================================================================
205{
206 IsComputed_ = false;
207 if (IsInitialized() == false) {
209 }
210
213
214 // extract local rows and columns
215 IFPACK_CHK_ERR(Extract(Matrix_in));
216
219
220 // factorize the matrix using LAPACK
221 if (NumRows_ != 0)
223
224 Label_ = "Ifpack_DenseContainer";
225
226 // not sure of count
227#ifdef IFPACK_FLOPCOUNTERS
228 ComputeFlops_ += 4.0 * NumRows_ * NumRows_ * NumRows_ / 3;
229#endif
230 IsComputed_ = true;
231
232 return(0);
233}
234
235//==============================================================================
237{
238 if (IsComputed() == false)
239 IFPACK_CHK_ERR(-3);
240
243 }
244 else
245 IFPACK_CHK_ERR(RHS_.Multiply('N','N', 1.0,Matrix_,LHS_,0.0));
246
247#ifdef IFPACK_FLOPCOUNTERS
249#endif
250 return(0);
251}
252
253//==============================================================================
254std::ostream& Ifpack_DenseContainer::Print(std::ostream & os) const
255{
256 using std::endl;
257
258 os << "================================================================================" << endl;
259 os << "Ifpack_DenseContainer" << endl;
260 os << "Number of rows = " << NumRows() << endl;
261 os << "Number of vectors = " << NumVectors() << endl;
262 os << "IsInitialized() = " << IsInitialized() << endl;
263 os << "IsComputed() = " << IsComputed() << endl;
264#ifdef IFPACK_FLOPCOUNTERS
265 os << "Flops in Compute() = " << ComputeFlops() << endl;
266 os << "Flops in ApplyInverse() = " << ApplyInverseFlops() << endl;
267#endif
268 os << "================================================================================" << endl;
269 os << endl;
270
271 return(os);
272}
#define IFPACK_CHK_ERR(ifpack_err)
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 Multiply(char TransA, char TransB, double ScalarAB, const Epetra_SerialDenseMatrix &A, const Epetra_SerialDenseMatrix &B, double ScalarThis)
int Reshape(int NumRows, int NumCols)
virtual int Solve(void)
int SetMatrix(Epetra_SerialDenseMatrix &A)
virtual int Factor(void)
int SetVectors(Epetra_SerialDenseMatrix &X, Epetra_SerialDenseMatrix &B)
int NumRows_
Number of rows in the container.
std::string Label_
Label for this object.
virtual int NumRows() const
Returns the number of rows of the matrix and LHS/RHS.
virtual int NumVectors() const
Returns the number of vectors in LHS/RHS.
virtual int ApplyInverse()
Apply the inverse of the matrix to RHS, results are stored in LHS.
virtual bool IsInitialized() const
Returns true is the container has been successfully initialized.
virtual int Apply()
Apply the matrix to RHS, results are stored in LHS.
virtual int SetMatrixElement(const int row, const int col, const double value)
Set the matrix element (row,col) to value.
Epetra_SerialDenseMatrix Matrix_
Dense matrix.
bool IsComputed_
If true, the container has been successfully computed.
virtual double ApplyInverseFlops() const
Returns the flops in ApplyInverse().
bool KeepNonFactoredMatrix_
If true, keeps a copy of the non-factored matrix.
virtual int Compute(const Epetra_RowMatrix &Matrix_in)
Finalizes the linear system matrix and prepares for the application of the inverse.
double ApplyFlops_
Flops in Apply().
Epetra_SerialDenseMatrix LHS_
Dense vector representing the LHS.
virtual int Initialize()
Initialize the container.
virtual bool IsComputed() const
Returns true is the container has been successfully computed.
virtual int Extract(const Epetra_RowMatrix &Matrix_in)
Extract the submatrices identified by the ID set int ID().
virtual double ComputeFlops() const
Returns the flops in Compute().
virtual const Epetra_SerialDenseMatrix & LHS() const
Returns the dense vector containing the LHS.
bool IsInitialized_
If true, the container has been successfully initialized.
int NumVectors_
Number of vectors in the container.
virtual std::ostream & Print(std::ostream &os) const
Prints basic information on iostream. This function is used by operator<<.
virtual int & ID(const int i)
Returns the ID associated to local row i.
Epetra_IntSerialDenseVector ID_
Sets of local rows.
Epetra_SerialDenseSolver Solver_
Dense solver (solution will be get using LAPACK).
Epetra_SerialDenseMatrix RHS_
Dense vector representing the RHS.
virtual const Epetra_SerialDenseMatrix & RHS() const
Returns the dense vector containing the RHS.
double ComputeFlops_
Flops in Compute().
Epetra_SerialDenseMatrix NonFactoredMatrix_
Dense matrix, that contains the non-factored matrix.
double ApplyInverseFlops_
Flops in ApplyInverse().
virtual const Epetra_IntSerialDenseVector & ID() const
Returns the integer dense vector of IDs.