#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_CrsMatrix.h"
class TridiagonalCrsMatrix : public Epetra_CrsMatrix {
public:
TridiagonalCrsMatrix(const Epetra_Map & Map,
double a,
double diag, double c) :
Epetra_CrsMatrix(Copy,Map,3)
{
int NumGlobalElements = Map.NumGlobalElements();
int NumMyElements = Map.NumMyElements();
int * MyGlobalElements = new int [NumMyElements];
Map.MyGlobalElements( MyGlobalElements );
double *Values = new double[2];
Values[0] = a; Values[1] = c;
int *Indices = new int[2];
int NumEntries;
for( int i=0 ; i<NumMyElements; ++i ) {
if (MyGlobalElements[i]==0) {
Indices[0] = 1;
NumEntries = 1;
} else if (MyGlobalElements[i] == NumGlobalElements-1) {
Indices[0] = NumGlobalElements-2;
NumEntries = 1;
} else {
Indices[0] = MyGlobalElements[i]-1;
Indices[1] = MyGlobalElements[i]+1;
NumEntries = 2;
}
InsertGlobalValues(MyGlobalElements[i], NumEntries, Values, Indices);
InsertGlobalValues(MyGlobalElements[i], 1, &diag, MyGlobalElements+i);
}
FillComplete();
}
};
bool CrsMatrix2MATLAB( const Epetra_CrsMatrix & A )
{
int MyPID = A.Comm().MyPID();
int NumProc = A.Comm().NumProc();
if( A.IndicesAreLocal() == false ) {
if( MyPID == 0 ) {
cerr << "ERROR in "<< __FILE__ << ", line " << __LINE__ << endl;
cerr << "Function CrsMatrix2MATLAB accepts\n";
cerr << "transformed matrices ONLY. Please call A.FillComplete()\n";
cerr << "on your matrix A to that purpose.\n";
cerr << "Now returning...\n";
}
return false;
}
int NumMyRows = A.NumMyRows();
int NumNzRow;
int NumEntries;
int NumGlobalRows;
int GlobalRow;
int NumGlobalNonzeros;
NumGlobalRows = A.NumGlobalRows();
NumGlobalNonzeros = A.NumGlobalNonzeros();
int IndexBase = A.IndexBase();
if( IndexBase == 0 ) IndexBase = 1;
if( MyPID==0 ) {
cout << "A = spalloc(";
cout << NumGlobalRows << ',' << NumGlobalRows;
cout << ',' << NumGlobalNonzeros << ");\n";
}
for( int Proc=0 ; Proc<NumProc ; ++Proc ) {
if( MyPID == Proc ) {
cout << "% On proc " << Proc << ": ";
cout << NumMyRows << " rows and ";
cout << A.NumMyNonzeros() << " nonzeros\n";
for( int MyRow=0 ; MyRow<NumMyRows ; ++MyRow ) {
GlobalRow = A.GRID(MyRow);
NumNzRow = A.NumMyEntries(MyRow);
double *Values = new double[NumNzRow];
int *Indices = new int[NumNzRow];
A.ExtractMyRowCopy(MyRow, NumNzRow,
NumEntries, Values, Indices);
for( int j=0 ; j<NumEntries ; ++j ) {
cout << "A(" << GlobalRow + IndexBase
<< "," << A.GCID(Indices[j]) + IndexBase
<< ") = " << Values[j] << ";\n";
}
delete Values;
delete Indices;
}
}
A.Comm().Barrier();
if( MyPID == 0 ) {
cout << " %End of Matrix Output\n";
}
}
return true;
}
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);
TridiagonalCrsMatrix A( Map, -1.0, 2.0, -1.0);
CrsMatrix2MATLAB( A );
#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