#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"
#include "Epetra_Import.h"
#include "Epetra_IntSerialDenseVector.h"
int find( int key, int vector[], int Length )
{
for( int i=0 ; i<Length ; ++i ) {
if( vector[i] == key ) return i;
}
return -1;
}
class TriDiagonalOperator : public Epetra_Operator
{
public:
TriDiagonalOperator( double diag_minus_one,
double diag,
double diag_plus_one,
const Epetra_Map & Map) :
Map_( Map ),
diag_minus_one_(diag_minus_one),
diag_(diag),
diag_plus_one_(diag_plus_one)
{
NumMyElements_ = Map_.NumMyElements();
NumGlobalElements_ = Map_.NumGlobalElements();
int* MyGlobalElements = new int[NumMyElements_];
Map_.MyGlobalElements(MyGlobalElements);
int count=0;
for( int i=0 ; i<NumMyElements_ ; ++i ) {
int globalIndex = MyGlobalElements[i];
if( globalIndex>0 )
if( Map.LID(globalIndex-1) == -1 ) ++count;
if( globalIndex<NumGlobalElements_-1 )
if( Map.LID(globalIndex+1) == -1 ) ++count;
++count;
}
int Length = count;
int* ListOfNodes = new int[Length];
count=0;
for( int i=0 ; i<NumMyElements_ ; ++i ) {
int globalIndex = MyGlobalElements[i];
if( globalIndex>0 ) {
if( Map.LID(globalIndex-1) == -1 )
if( find( globalIndex-1, ListOfNodes, Length) == -1 ) {
ListOfNodes[count] = globalIndex-1;
++count;
}
}
if( globalIndex<NumGlobalElements_-1 ) {
if( Map.LID(globalIndex+1) == -1 ) {
if( find( globalIndex+1, ListOfNodes, Length) == -1 ) {
ListOfNodes[count] = globalIndex+1;
++count;
}
}
}
ListOfNodes[count] = globalIndex;
++count;
}
ImportMap_ = new Epetra_Map(-1,count,ListOfNodes,0,Map_.Comm());
Importer_ = new Epetra_Import(*ImportMap_,Map_);
delete[] MyGlobalElements;
delete[] ListOfNodes;
return;
}
int Apply( const Epetra_MultiVector & X,
Epetra_MultiVector & Y ) const
{
cout << X;
Epetra_MultiVector Xext((*ImportMap_),X.NumVectors());
Xext.Import(X,*Importer_,Insert);
for( int i=0 ; i<X.MyLength() ; ++i ) {
int globalRow = Map_.GID(i);
int iMinusOne = (*ImportMap_).LID(globalRow-1);
int iPlusOne = (*ImportMap_).LID(globalRow+1);
printf("%d %d %d\n", globalRow, iMinusOne, iPlusOne);
for( int vec=0 ; vec<X.NumVectors() ; ++vec ) {
Y[vec][i] = diag_ * X[vec][i];
if( iMinusOne != -1 )
Y[vec][i] += diag_minus_one_ * Xext[vec][iMinusOne];
if( iPlusOne != -1 )
Y[vec][i] += diag_plus_one_ * Xext[vec][iPlusOne];
}
}
return true;
}
int SetUseTranspose( bool UseTranspose)
{
return(0);
}
int ApplyInverse( const Epetra_MultiVector & X,
Epetra_MultiVector & Y ) const
{
return 0;
}
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_;
int NumMyElements_;
int NumGlobalElements_;
double diag_minus_one_;
double diag_;
double diag_plus_one_;
Epetra_Import *Importer_;
Epetra_Map *ImportMap_;
};
int main(int argc, char *argv[]) {
#ifdef HAVE_MPI
MPI_Init(&argc, &argv);
Epetra_MpiComm Comm(MPI_COMM_WORLD);
#else
Epetra_SerialComm Comm;
#endif
int NumGlobalElements( 5 );
Epetra_Map Map(NumGlobalElements,0,Comm );
Epetra_Vector x(Map);
Epetra_Vector y(Map);
int NumMyElements = Map.NumMyElements();
Epetra_IntSerialDenseVector MyGlobalElements(NumMyElements);
Map.MyGlobalElements( MyGlobalElements.Values() );
for( int i=0 ; i<NumMyElements ; ++i )
x[i] = 1.0*MyGlobalElements[i];
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