IFPACK  Development
 All Classes Namespaces Files Functions Variables Enumerations Friends Pages
Ifpack_DenseContainer.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_DenseContainer.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  // 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) {
81  IFPACK_CHK_ERR(Solver_.SetMatrix(Matrix_));
82  IFPACK_CHK_ERR(Solver_.SetVectors(LHS_,RHS_));
83  }
84 
85  IsInitialized_ = true;
86  return(0);
87 
88 }
89 
90 //==============================================================================
91 double& Ifpack_DenseContainer::LHS(const int i, const int Vector)
92 {
93  return(LHS_.A()[Vector * NumRows_ + i]);
94 }
95 
96 //==============================================================================
97 double& Ifpack_DenseContainer::RHS(const int i, const int Vector)
98 {
99  return(RHS_.A()[Vector * NumRows_ + i]);
100 }
101 
102 //==============================================================================
104 SetMatrixElement(const int row, const int col, const double value)
105 {
106  if (IsInitialized() == false)
107  IFPACK_CHK_ERR(Initialize());
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)
132  IFPACK_CHK_ERR(Solver_.Solve());
133 
134 #ifdef IFPACK_FLOPCOUNTERS
135  ApplyInverseFlops_ += 2.0 * NumVectors_ * NumRows_ * NumRows_;
136 #endif
137  return(0);
138 }
139 
140 //==============================================================================
141 int& Ifpack_DenseContainer::ID(const int i)
142 {
143  return(ID_[i]);
144 }
145 
146 //==============================================================================
147 // FIXME: optimize performances of this guy...
148 int Ifpack_DenseContainer::Extract(const Epetra_RowMatrix& Matrix_in)
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) {
208  IFPACK_CHK_ERR(Initialize());
209  }
210 
211  if (KeepNonFactoredMatrix_)
212  NonFactoredMatrix_ = Matrix_;
213 
214  // extract local rows and columns
215  IFPACK_CHK_ERR(Extract(Matrix_in));
216 
217  if (KeepNonFactoredMatrix_)
218  NonFactoredMatrix_ = Matrix_;
219 
220  // factorize the matrix using LAPACK
221  if (NumRows_ != 0)
222  IFPACK_CHK_ERR(Solver_.Factor());
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 
241  if (KeepNonFactoredMatrix_) {
242  IFPACK_CHK_ERR(RHS_.Multiply('N','N', 1.0,NonFactoredMatrix_,LHS_,0.0));
243  }
244  else
245  IFPACK_CHK_ERR(RHS_.Multiply('N','N', 1.0,Matrix_,LHS_,0.0));
246 
247 #ifdef IFPACK_FLOPCOUNTERS
248  ApplyFlops_ += 2 * NumRows_ * NumRows_;
249 #endif
250  return(0);
251 }
252 
253 //==============================================================================
254 std::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 }
virtual double ComputeFlops() const
Returns the flops in Compute().
virtual int Factor(void)
virtual const Epetra_IntSerialDenseVector & ID() const
Returns the integer dense vector of IDs.
virtual int Compute(const Epetra_RowMatrix &Matrix_in)
Finalizes the linear system matrix and prepares for the application of the inverse.
virtual std::ostream & Print(std::ostream &os) const
Prints basic information on iostream. This function is used by operator&lt;&lt;.
virtual int SetMatrixElement(const int row, const int col, const double value)
Set the matrix element (row,col) to value.
virtual int Solve(void)
int SetVectors(Epetra_SerialDenseMatrix &X, Epetra_SerialDenseMatrix &B)
virtual int NumVectors() const
Returns the number of vectors in LHS/RHS.
virtual double ApplyInverseFlops() const
Returns the flops in ApplyInverse().
double * A() const
virtual bool IsInitialized() const
Returns true is the container has been successfully initialized.
virtual int MaxNumEntries() const =0
virtual int Apply()
Apply the matrix to RHS, results are stored in LHS.
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 int NumMyRows() const =0
virtual int Initialize()
Initialize the container.
int Reshape(int NumRows, int NumCols)
virtual int ApplyInverse()
Apply the inverse of the matrix to RHS, results are stored in LHS.
virtual int NumRows() const
Returns the number of rows of the matrix and LHS/RHS.
int SetMatrix(Epetra_SerialDenseMatrix &A)
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 ExtractMyRowCopy(int MyRow, int Length, int &NumEntries, double *Values, int *Indices) const =0
virtual const Epetra_SerialDenseMatrix & RHS() const
Returns the dense vector containing the RHS.