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);
83 Epetra_MpiComm Comm( MPI_COMM_WORLD );
85 Epetra_SerialComm Comm;
94 cout << Epetra_Version() << endl << endl;
96 int rank = Comm.MyPID();
98 if (
verbose) cout << Comm <<endl;
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];
114 Epetra_SerialDenseSolver ed_solver;
115 Epetra_SerialDenseMatrix * ed_Matrix;
117 bool Transpose =
false;
122 ed_solver.SolveWithTranspose(Transpose);
123 ed_solver.SolveToRefinedSolution(Refine);
126 ed_Matrix =
new Epetra_SerialDenseMatrix(5,5);
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;
149 Epetra_SerialDenseVector LHS(
Copy, X,
N);
150 Epetra_SerialDenseVector
RHS(
Copy,
B,
N);
152 Epetra_SerialDenseVector ed_LHS(
Copy, ed_X,
N);
153 Epetra_SerialDenseVector ed_RHS(
Copy, ed_B,
N);
158 ed_solver.SetMatrix(*ed_Matrix);
159 ed_solver.SetVectors(ed_LHS, ed_RHS);
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;
175 Epetra_SerialDenseMatrix * sdfac = ed_solver.FactoredMatrix();
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] <<
" ";
virtual int Solve(void)
Computes the solution X to AX = B for the this matrix and the B provided to SetVectors()..
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. ...
Ifpack_SerialTriDiMatrix * FactoredMatrix() const
Returns pointer to factored matrix (assuming factorization has been performed).
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)
int main(int argc, char *argv[])
Ifpack_SerialTriDiSolver: A class for solving TriDi linear problems.
void printHeading(const char *heading)
virtual void Print(std::ostream &os) const
Print service methods; defines behavior of ostream << operator.
double * DL()
Returns pointer to the this matrix.
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)
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)