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"
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()
const int N
#define RHS(a)
Definition MatGenFD.c:60
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
virtual void Print(std::ostream &os) const
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 main(int argc, char *argv[])
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)
void printMat(const char *name, Ifpack_SerialTriDiMatrix &matrix)
void printHeading(const char *heading)
bool Residual(int N, int NRHS, double *A, int LDA, bool Transpose, double *X, int LDX, double *B, int LDB, double *resid)