Ifpack Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
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
59int 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
67bool Residual( int N, int NRHS, double * A, int LDA, bool Transpose,
68 double * X, int LDX, double * B, int LDB, double * resid);
69
70int matrixCpyCtr(bool verbose, bool debug);
71
72void printHeading(const char* heading);
73void printMat(const char* name, Ifpack_SerialTriDiMatrix& matrix);
74void printArray(double* array, int length);
75
76using namespace std;
77
78int 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
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
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
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
206bool 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
229void printHeading(const char* heading) {
230 cout << "\n==================================================================\n";
231 cout << heading << endl;
232 cout << "==================================================================\n";
233}
234
235//=========================================================================
236// prints SerialTriDiMatrix/Vector with formatting
237void 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
246void 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
Copy
std::string Epetra_Version()
#define RHS(a)
Definition: MatGenFD.c:60
float NRM2(const int N, const float *X, const int INCX=1) const
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
int MyPID() const
double * A() const
virtual void Print(std::ostream &os) const
virtual int Solve(void)
void SolveToRefinedSolution(bool Flag)
int SetMatrix(Epetra_SerialDenseMatrix &A)
void SolveWithTranspose(bool Flag)
Epetra_SerialDenseMatrix * FactoredMatrix() const
int SetVectors(Epetra_SerialDenseMatrix &X, Epetra_SerialDenseMatrix &B)
Ifpack_SerialTriDiMatrix: A class for constructing and using real double precision general TriDi matr...
double * DL()
Returns pointer to the this matrix.
virtual void Print(std::ostream &os) const
Print service methods; defines behavior of ostream << operator.
Ifpack_SerialTriDiSolver: A class for solving TriDi linear problems.
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()....
Ifpack_SerialTriDiMatrix * FactoredMatrix() const
Returns pointer to factored matrix (assuming factorization has been performed).
int SetVectors(Epetra_SerialDenseMatrix &X, Epetra_SerialDenseMatrix &B)
Sets the pointers for left and right hand side vector(s).
void SolveToRefinedSolution(bool Flag)
Causes all solves to compute solution to best ability using iterative refinement.
void SolveWithTranspose(bool Flag)
Causes equilibration to be called just before the matrix factorization as part of the call to Factor.
const int N
int main(int argc, char *argv[])
Definition: tridi_main.cpp:78
int matrixCpyCtr(bool verbose, bool debug)
int check(Ifpack_SerialTriDiSolver &solver, double *A1, int LDA, int N1, int NRHS1, double OneNorm1, double *B1, int LDB1, double *X1, int LDX1, bool Transpose, bool verbose)
void printArray(double *array, int length)
Definition: tridi_main.cpp:246
void printMat(const char *name, Ifpack_SerialTriDiMatrix &matrix)
Definition: tridi_main.cpp:237
void printHeading(const char *heading)
Definition: tridi_main.cpp:229
bool Residual(int N, int NRHS, double *A, int LDA, bool Transpose, double *X, int LDX, double *B, int LDB, double *resid)
Definition: tridi_main.cpp:206