Ifpack Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
tridi_main.cpp
Go to the documentation of this file.
1 //@HEADER
2 // ************************************************************************
3 //
4 // Epetra: Linear Algebra Services Package
5 // Copyright 2011 Sandia Corporation
6 //
7 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8 // the U.S. Government retains certain rights in this software.
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 "Epetra_Map.h"
44 #include "Epetra_Time.h"
46 #include "Epetra_SerialDenseVector.h"
47 #include "Epetra_SerialDenseSolver.h"
48 #include "Epetra_SerialDenseMatrix.h"
50 #ifdef EPETRA_MPI
51 #include "Epetra_MpiComm.h"
52 #include <mpi.h>
53 #endif
54 #include "Epetra_SerialComm.h"
55 #include "Epetra_Version.h"
56 
57 // prototypes
58 
59 int check(Ifpack_SerialTriDiSolver & solver, double * A1, int LDA,
60  int N1, int NRHS1, double OneNorm1,
61  double * B1, int LDB1,
62  double * X1, int LDX1,
63  bool Transpose, bool verbose);
64 
65 
66 
67 bool Residual( int N, int NRHS, double * A, int LDA, bool Transpose,
68  double * X, int LDX, double * B, int LDB, double * resid);
69 
70 int matrixCpyCtr(bool verbose, bool debug);
71 
72 void printHeading(const char* heading);
73 void printMat(const char* name, Ifpack_SerialTriDiMatrix& matrix);
74 void printArray(double* array, int length);
75 
76 using namespace std;
77 
78 int main(int argc, char *argv[])
79 {
80 
81 #ifdef HAVE_MPI
82  MPI_Init(&argc,&argv);
83  Epetra_MpiComm Comm( MPI_COMM_WORLD );
84 #else
85  Epetra_SerialComm Comm;
86 #endif
87 
88  bool verbose = false;
89 
90  // Check if we should print results to standard out
91  verbose = true;
92 
93  if (verbose && Comm.MyPID()==0)
94  cout << Epetra_Version() << endl << endl;
95 
96  int rank = Comm.MyPID();
97 
98  if (verbose) cout << Comm <<endl;
99 
100  // Redefine verbose to only print on PE 0
101  if (verbose && rank!=0) verbose = false;
102 
103  int N = 5;
104  int NRHS = 5;
105  double * X = new double[NRHS];
106  double * ed_X = new double[NRHS];
107 
108  double * B = new double[NRHS];
109  double * ed_B = new double[NRHS];
110 
112  Ifpack_SerialTriDiMatrix * Matrix;
113 
114  Epetra_SerialDenseSolver ed_solver;
115  Epetra_SerialDenseMatrix * ed_Matrix;
116 
117  bool Transpose = false;
118  bool Refine = false;
119  solver.SolveWithTranspose(Transpose);
120  solver.SolveToRefinedSolution(Refine);
121 
122  ed_solver.SolveWithTranspose(Transpose);
123  ed_solver.SolveToRefinedSolution(Refine);
124 
125  Matrix = new Ifpack_SerialTriDiMatrix(5,true);
126  ed_Matrix = new Epetra_SerialDenseMatrix(5,5);
127 
128  for(int i=0;i<N;++i) {
129  B[i] = ed_B[i] =2;
130  Matrix->D()[i]=2.0;
131  if(i<(N-1)) {
132  Matrix->DL()[i]=-1.0;
133  if(i!=2) Matrix->DU()[i]=-1.0;
134  }
135  }
136 
137  Matrix->Print(std::cout);
138 
139  double * ed_a = ed_Matrix->A();
140  for(int i=0;i<N;++i)
141  for(int j=0;j<N;++j) {
142  if(i==j) ed_a[j*N+i] = 2.0;
143  else if(abs(i-j) == 1) ed_a[j*N+i] = -1.0;
144  else ed_a[j*N + i] = 0;
145  if(i==2&&j==3) ed_a[j*N+i] = 0.0;
146  }
147 
148 
149  Epetra_SerialDenseVector LHS(Copy, X, N);
151 
152  Epetra_SerialDenseVector ed_LHS(Copy, ed_X, N);
153  Epetra_SerialDenseVector ed_RHS(Copy, ed_B, N);
154 
155  solver.SetMatrix(*Matrix);
156  solver.SetVectors(LHS, RHS);
157 
158  ed_solver.SetMatrix(*ed_Matrix);
159  ed_solver.SetVectors(ed_LHS, ed_RHS);
160 
161  solver.Solve();
162  ed_solver.Solve();
163 
164  std::cout << " LHS vals are: "<<std::endl;
165  bool TestPassed=true;
166  for(int i=0;i<N;++i) {
167  std::cout << "["<<i<<"] "<< LHS(i)<<" "<<ed_LHS(i)<<" delta "<<LHS(i)-ed_LHS(i)<<std::endl;
168  if( fabs( (LHS(i)- ed_LHS(i))/(LHS(i)+ ed_LHS(i)) ) > 1.0e-12 ) {
169  TestPassed = false;
170  std::cout << " not equal for "<<i<<" delta "<< LHS(i)- ed_LHS(i)<<std::endl;
171  }
172  }
173 
174  Ifpack_SerialTriDiMatrix * tdfac = solver.FactoredMatrix();
175  Epetra_SerialDenseMatrix * sdfac = ed_solver.FactoredMatrix();
176 
177  std::cout << " Tri Di Factored "<<std::endl;
178  tdfac->Print(std::cout);
179  std::cout << " Dense Factored "<<std::endl;
180  sdfac->Print(std::cout);
181 
182  delete Matrix;
183  delete ed_Matrix;
184  delete [] X;
185  delete [] ed_X;
186  delete [] B;
187  delete [] ed_B;
188 
189 
190  if (!TestPassed) {
191  cout << "Test `TestRelaxation.exe' failed!" << endl;
192  exit(EXIT_FAILURE);
193  }
194 
195 #ifdef HAVE_MPI
196  MPI_Finalize();
197 #endif
198 
199  cout << endl;
200  cout << "Test `TestRelaxation.exe' passed!" << endl;
201  cout << endl;
202  return(EXIT_SUCCESS);
203 }
204 
205 
206 bool Residual( int N, int NRHS, double * A, int LDA, bool Transpose,
207  double * X, int LDX, double * B, int LDB, double * resid) {
208 
209  Epetra_BLAS Blas;
210  char Transa = 'N';
211  if (Transpose) Transa = 'T';
212  Blas.GEMM(Transa, 'N', N, NRHS, N, -1.0, A, LDA,
213  X, LDX, 1.0, B, LDB);
214  bool OK = true;
215  for (int i=0; i<NRHS; i++) {
216  resid[i] = Blas.NRM2(N, B+i*LDB);
217  if (resid[i]>1.0E-7) OK = false;
218  }
219 
220  return(OK);
221 }
222 
223 
224 //=========================================================================
225 
226 //=========================================================================
227 //=========================================================================
228 // prints section heading with spacers/formatting
229 void printHeading(const char* heading) {
230  cout << "\n==================================================================\n";
231  cout << heading << endl;
232  cout << "==================================================================\n";
233 }
234 
235 //=========================================================================
236 // prints SerialTriDiMatrix/Vector with formatting
237 void printMat(const char* name, Ifpack_SerialTriDiMatrix& matrix) {
238  //cout << "--------------------" << endl;
239  cout << "*** " << name << " ***" << endl;
240  cout << matrix;
241  //cout << "--------------------" << endl;
242 }
243 
244 //=========================================================================
245 // prints double* array with formatting
246 void printArray(double* array, int length) {
247  cout << "user array (size " << length << "): ";
248  for(int i = 0; i < length; i++)
249  cout << array[i] << " ";
250  cout << endl;
251 }
252 
void GEMM(const char TRANSA, const char TRANSB, const int M, const int N, const int K, const float ALPHA, const float *A, const int LDA, const float *B, const int LDB, const float BETA, float *C, const int LDC) const
virtual int Solve(void)
Computes the solution X to AX = B for the this matrix and the B provided to SetVectors()..
float NRM2(const int N, const float *X, const int INCX=1) const
void SolveWithTranspose(bool Flag)
Causes equilibration to be called just before the matrix factorization as part of the call to Factor...
void SolveToRefinedSolution(bool Flag)
Causes all solves to compute solution to best ability using iterative refinement. ...
int matrixCpyCtr(bool verbose, bool debug)
int SetMatrix(Ifpack_SerialTriDiMatrix &A)
Sets the pointers for coefficient matrix.
void printMat(const char *name, Ifpack_SerialTriDiMatrix &matrix)
Definition: tridi_main.cpp:237
Ifpack_SerialTriDiSolver: A class for solving TriDi linear problems.
void printHeading(const char *heading)
int SetVectors(Epetra_SerialDenseMatrix &X, Epetra_SerialDenseMatrix &B)
void SolveWithTranspose(bool Flag)
int MyPID() const
std::string Epetra_Version()
Ifpack_SerialTriDiMatrix * FactoredMatrix() const
Returns pointer to factored matrix (assuming factorization has been performed).
double * A() const
double * DL()
Returns pointer to the this matrix.
int check(Epetra_CrsGraph &L, Epetra_CrsGraph &U, Ifpack_IlukGraph &LU, int NumGlobalRows1, int NumMyRows1, int LevelFill1, bool verbose)
int main(int argc, char *argv[])
const int N
virtual int Solve(void)
void SolveToRefinedSolution(bool Flag)
virtual void Print(std::ostream &os) const
Print service methods; defines behavior of ostream &lt;&lt; operator.
int SetMatrix(Epetra_SerialDenseMatrix &A)
Epetra_SerialDenseMatrix * FactoredMatrix() const
virtual void Print(std::ostream &os) const
Ifpack_SerialTriDiMatrix: A class for constructing and using real double precision general TriDi matr...
int SetVectors(Epetra_SerialDenseMatrix &X, Epetra_SerialDenseMatrix &B)
Sets the pointers for left and right hand side vector(s).
bool Residual(int N, int NRHS, double *A, int LDA, bool Transpose, double *X, int LDX, double *B, int LDB, double *resid)
void printArray(double *array, int length)
#define RHS(a)
Definition: MatGenFD.c:60