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,
63 bool Transpose,
bool verbose);
67bool Residual(
int N,
int NRHS,
double * A,
int LDA,
bool Transpose,
68 double * X,
int LDX,
double * B,
int LDB,
double * resid);
78int main(
int argc,
char *argv[])
82 MPI_Init(&argc,&argv);
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];
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);
206bool 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.0E-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] <<
" ";
std::string Epetra_Version()
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
virtual void Print(std::ostream &os) const
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.
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)