#include "Didasko_ConfigDefs.h"
#if defined(HAVE_DIDASKO_EPETRA)
#include "Epetra_ConfigDefs.h"
#ifdef HAVE_MPI
#include "mpi.h"
#include "Epetra_MpiComm.h"
#else
#include "Epetra_SerialComm.h"
#endif
#include "Epetra_Map.h"
#include "Epetra_Vector.h"
#include "Epetra_MultiVector.h"
#include "Epetra_Operator.h"
class TriDiagonalOperator : public Epetra_Operator
{
public:
TriDiagonalOperator( double diag_minus_one,
double diag,
double diag_plus_one,
Epetra_Map & Map) :
Map_(Map),
diag_minus_one_(diag_minus_one),
diag_(diag),
diag_plus_one_(diag_plus_one)
{}
int Apply( const Epetra_MultiVector & X,
Epetra_MultiVector & Y ) const
{
int Length = X.MyLength();
for( int vec=0 ; vec<X.NumVectors() ; ++vec ) {
if( Length == 1 ) {
Y[vec][0] = diag_ * X[vec][0];
break;
}
Y[vec][0] = diag_ * X[vec][0] + diag_plus_one_ * X[vec][1];
for( int i=1 ; i<Length-1 ; ++i ) {
Y[vec][i] = diag_ * X[vec][i] + diag_plus_one_ * X[vec][i+1]
+ diag_minus_one_ * X[vec][i-1];
}
Y[vec][Length-1] = diag_ * X[vec][Length-1]
+ diag_minus_one_ * X[vec][Length-2];
}
return true;
}
int SetUseTranspose( bool UseTranspose)
{
return(-1);
}
int ApplyInverse( const Epetra_MultiVector & X,
Epetra_MultiVector & Y ) const
{
return(-1);
}
double NormInf() const
{
return(abs(diag_) + abs(diag_minus_one_) + abs(diag_plus_one_));
}
const char * Label () const
{
return("TriDiagonalOperator");
}
bool UseTranspose() const
{
return(false);
}
bool HasNormInf () const
{
return(true);
}
const Epetra_Comm & Comm() const
{
return(Map_.Comm());
}
const Epetra_Map & OperatorDomainMap() const
{
return(Map_);
}
const Epetra_Map & OperatorRangeMap() const
{
return(Map_);
}
private:
Epetra_Map Map_;
double diag_minus_one_;
double diag_;
double diag_plus_one_;
};
int main(int argc, char *argv[]) {
#ifdef HAVE_MPI
MPI_Init(&argc, &argv);
Epetra_MpiComm Comm(MPI_COMM_WORLD);
#else
Epetra_SerialComm Comm;
#endif
if( Comm.NumProc() != 1 ) {
if( Comm.MyPID() == 0 ) {
cerr << "This is mono-process example\n"
<< "Please run with one processo only\n";
}
#ifdef HAVE_MPI
MPI_Finalize();
#endif
exit(EXIT_SUCCESS);
}
int NumGlobalElements( 5 );
Epetra_Map Map(NumGlobalElements,0,Comm );
Epetra_Vector x(Map);
Epetra_Vector y(Map);
x.PutScalar(1.0);
TriDiagonalOperator TriDiagOp(-1.0,2.0,-1.0,Map);
TriDiagOp.Apply(x,y);
cout << x;
cout << y;
#ifdef HAVE_MPI
MPI_Finalize();
#endif
return( EXIT_SUCCESS );
}
#else
#include <stdlib.h>
#include <stdio.h>
int main(int argc, char *argv[])
{
puts("Please configure Didasko with:\n"
"--enable-epetra");
return 0;
}
#endif