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"
51 #include "Epetra_MpiComm.h"
54 #include "Epetra_SerialComm.h"
55 #include "Epetra_Version.h"
60 int N1,
int NRHS1,
double OneNorm1,
61 double *
B1,
int LDB1,
62 double * X1,
int LDX1,
67 bool Residual(
int N,
int NRHS,
double *
A,
int LDA,
bool Transpose,
68 double * X,
int LDX,
double *
B,
int LDB,
double * resid);
78 int main(
int argc,
char *argv[])
82 MPI_Init(&argc,&argv);
93 if (verbose && Comm.
MyPID()==0)
96 int rank = Comm.
MyPID();
98 if (verbose) cout << Comm <<endl;
101 if (verbose && rank!=0) verbose =
false;
105 double * X =
new double[NRHS];
106 double * ed_X =
new double[NRHS];
108 double *
B =
new double[NRHS];
109 double * ed_B =
new double[NRHS];
117 bool Transpose =
false;
128 for(
int i=0;i<
N;++i) {
132 Matrix->
DL()[i]=-1.0;
133 if(i!=2) Matrix->
DU()[i]=-1.0;
137 Matrix->
Print(std::cout);
139 double * ed_a = ed_Matrix->
A();
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;
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 ) {
170 std::cout <<
" not equal for "<<i<<
" delta "<< LHS(i)- ed_LHS(i)<<std::endl;
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);
191 cout <<
"Test `TestRelaxation.exe' failed!" << endl;
200 cout <<
"Test `TestRelaxation.exe' passed!" << endl;
202 return(EXIT_SUCCESS);
206 bool Residual(
int N,
int NRHS,
double *
A,
int LDA,
bool Transpose,
207 double * X,
int LDX,
double *
B,
int LDB,
double * resid) {
211 if (Transpose) Transa =
'T';
212 Blas.
GEMM(Transa,
'N', N, NRHS, N, -1.0, A, LDA,
213 X, LDX, 1.0, B, LDB);
215 for (
int i=0; i<NRHS; i++) {
216 resid[i] = Blas.
NRM2(N, B+i*LDB);
217 if (resid[i]>1.0
E-7) OK =
false;
230 cout <<
"\n==================================================================\n";
231 cout << heading << endl;
232 cout <<
"==================================================================\n";
239 cout <<
"*** " << name <<
" ***" << endl;
247 cout <<
"user array (size " << length <<
"): ";
248 for(
int i = 0; i < length; i++)
249 cout << array[i] <<
" ";
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)
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)
std::string Epetra_Version()
Ifpack_SerialTriDiMatrix * FactoredMatrix() const
Returns pointer to factored matrix (assuming factorization has been performed).
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[])
void SolveToRefinedSolution(bool Flag)
virtual void Print(std::ostream &os) const
Print service methods; defines behavior of ostream << 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)