Ifpack Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Ifpack_TriDiContainer.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"
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 
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) {
87  }
88 
89  IsInitialized_ = true;
90  return(0);
91 
92 }
93 
94 //==============================================================================
95 double& Ifpack_TriDiContainer::LHS(const int i, const int Vector)
96 {
97  return(LHS_.A()[Vector * NumRows_ + i]);
98 }
99 
100 //==============================================================================
101 double& Ifpack_TriDiContainer::RHS(const int i, const int Vector)
102 {
103  return(RHS_.A()[Vector * NumRows_ + i]);
104 }
105 
106 //==============================================================================
108 SetMatrixElement(const int row, const int col, const double value)
109 {
110  if (IsInitialized() == false)
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)
137 
138 #ifdef IFPACK_FLOPCOUNTERS
140 #endif
141  return(0);
142 }
143 
144 //==============================================================================
145 int& Ifpack_TriDiContainer::ID(const int i)
146 {
147  return(ID_[i]);
148 }
149 
150 //==============================================================================
151 // FIXME: optimize performances of this guy...
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) {
213  }
214 
217 
218  // extract local rows and columns
219  IFPACK_CHK_ERR(Extract(Matrix_in));
220 
223 
224  // factorize the matrix using LAPACK
225  if (NumRows_ != 0)
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 //==============================================================================
259 std::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 }
virtual int NumRows() const
Returns the number of rows of the matrix and LHS/RHS.
virtual std::ostream & Print(std::ostream &os) const
Prints basic information on iostream. This function is used by operator&lt;&lt;.
std::string Label_
Label for this object.
int Shape(int NumRowCol)
Set dimensions of a Ifpack_SerialTriDiMatrix object; init values to zero.
double ApplyInverseFlops_
Flops in ApplyInverse().
virtual int Solve(void)
Computes the solution X to AX = B for the this matrix and the B provided to SetVectors()..
virtual int NumVectors() const
Returns the number of vectors in LHS/RHS.
virtual const Epetra_IntSerialDenseVector & ID() const
Returns the integer dense vector of IDs.
virtual int Factor(void)
Computes the in-place LU factorization of the matrix using the LAPACK routine DGETRF.
virtual int SetMatrixElement(const int row, const int col, const double value)
Set the matrix element (row,col) to value.
bool IsInitialized_
If true, the container has been successfully initialized.
bool KeepNonFactoredMatrix_
If true, keeps a copy of the non-factored matrix.
virtual int Extract(const Epetra_RowMatrix &Matrix_in)
Extract the submatrices identified by the ID set int ID().
virtual const Epetra_SerialDenseMatrix & RHS() const
Returns the dense vector containing the RHS.
int SetMatrix(Ifpack_SerialTriDiMatrix &A)
Sets the pointers for coefficient matrix.
Epetra_SerialDenseMatrix RHS_
SerialDense vector representing the RHS.
virtual const Epetra_SerialDenseMatrix & LHS() const
Returns the dense vector containing the LHS.
int NumVectors_
Number of vectors in the container.
Ifpack_SerialTriDiSolver Solver_
TriDi solver (solution will be get using LAPACK).
int NumRows_
Number of rows in the container.
double * A() const
virtual int MaxNumEntries() const =0
virtual double ComputeFlops() const
Returns the flops in Compute().
Epetra_IntSerialDenseVector ID_
Sets of local rows.
double * A() const
Returns pointer to the this matrix.
virtual int NumMyRows() const =0
int Reshape(int NumRows, int NumCols)
virtual int Apply()
Apply the matrix to RHS, results are stored in LHS.
virtual double ApplyInverseFlops() const
Returns the flops in ApplyInverse().
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.
Ifpack_SerialTriDiMatrix Matrix_
TriDi matrix.
virtual int Initialize()
Initialize the container.
int Reshape(int NumRows, int NumCols)
Ifpack_SerialTriDiMatrix NonFactoredMatrix_
TriDi matrix, that contains the non-factored matrix.
int SetVectors(Epetra_SerialDenseMatrix &X, Epetra_SerialDenseMatrix &B)
Sets the pointers for left and right hand side vector(s).
virtual int ExtractMyRowCopy(int MyRow, int Length, int &NumEntries, double *Values, int *Indices) const =0
bool IsComputed_
If true, the container has been successfully computed.
#define IFPACK_CHK_ERR(ifpack_err)
double ComputeFlops_
Flops in Compute().
virtual bool IsInitialized() const
Returns true is the container has been successfully initialized.
virtual bool IsComputed() const
Returns true is the container has been successfully computed.
Epetra_SerialDenseMatrix LHS_
SerialDense vector representing the LHS.