#include "Didasko_ConfigDefs.h"
#if defined(HAVE_DIDASKO_EPETRA) && defined(HAVE_DIDASKO_TRIUTILS)
#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_MultiVector.h"
#include "Epetra_Vector.h"
#include "Epetra_RowMatrix.h"
#include "Trilinos_Util.h"
class MSRMatrix : public Epetra_Operator
{
public:
MSRMatrix(Epetra_Map Map, int * bindx, double * val) :
bindx_(bindx),
val_(val),
Map_(Map)
{
}
~MSRMatrix()
{}
int Apply(const Epetra_MultiVector & X, Epetra_MultiVector & Y ) const
{
int Nrows = bindx_[0]-1;
for( int i=0 ; i<Nrows ; i++ ) {
for( int vec=0 ; vec<X.NumVectors() ; ++vec ) {
Y[vec][i] = val_[i]*X[vec][i];
}
for( int j=bindx_[i] ; j<bindx_[i+1] ; j++ ) {
for( int vec=0 ; vec<X.NumVectors() ; ++vec ) {
Y[vec][bindx_[j]] += val_[j]*X[vec][bindx_[j]];
}
}
}
return 0;
}
int SetUseTranspose( bool UseTranspose)
{
return(-1);
}
int ApplyInverse( const Epetra_MultiVector & X,
Epetra_MultiVector & Y ) const
{
return(-1);
}
double NormInf() const
{
return -1;
}
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:
int * bindx_;
double * val_;
Epetra_Map Map_;
};
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 << "*ERR* can be used only with one process" << endl;
#ifdef HAVE_MPI
MPI_Finalize();
#endif
exit(EXIT_SUCCESS);
}
int N_global;
int N_nonzeros;
double * val = NULL;
int * bindx = NULL;
double * x = NULL, * b = NULL, * xexact = NULL;
FILE* fp = fopen("../HBMatrices/fidap005.rua", "r");
if (fp == 0)
{
cerr << "Matrix file not available" << endl;
#ifdef HAVE_MPI
MPI_Finalize();
#endif
exit(EXIT_SUCCESS);
}
fclose(fp);
Trilinos_Util_read_hb("../HBMatrices/fidap005.rua", 0,
&N_global, &N_nonzeros,
&val, &bindx,
&x, &b, &xexact);
Epetra_Map Map(N_global,0,Comm);
MSRMatrix A(Map,bindx,val);
Epetra_Vector xxx(Map);
Epetra_Vector yyy(Map);
xxx.Random();
A.Apply(xxx,yyy);
cout << yyy;
double norm2;
yyy.Norm2(&norm2);
cout << norm2 << endl;
if (val != NULL) free((void*)val);
if (bindx != NULL) free((void*)bindx);
if (x != NULL) free((void*)x);
if (b != NULL) free((void*)x);
if (xexact != NULL) free((void*)xexact);;
#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\n"
"--enable-triutils");
return(0);
}
#endif