examples/matmatmult_bug.cpp

Go to the documentation of this file.

//
// C++ Implementation: matmatmult_bug
//
// Description:
//
//
// Author: Roman Geus <roman.geus@psi.ch>, (C) 2004
//

#include <iostream>
#include <fstream>
#include <cstdlib>
#include <unistd.h>

#include "Epetra_ConfigDefs.h"
#ifdef HAVE_MPI
#include <mpi.h>
#include "Epetra_MpiComm.h"
#else
#include "Epetra_SerialComm.h"
#endif
#include "Epetra_CrsMatrix.h"
#include "Epetra_Map.h"
#include "EpetraExt_MatrixMatrix.h"

#include "balancedmtxreader.h"
#include "mtxreader.h"
#include "export_Epetra_CrsMatrix.h"

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.MyPID() == 0)
        cout << "running on " << comm.NumProc() << " processor(s)..." << endl;

    BaseMtxReader *reader = 0;

    // Y matrix
    string filename("Y.mtx");
    if (access(filename.c_str(), R_OK)) {
        cerr << "Unable to access Y matrix file.\n";
        exit(1);
    }
    reader = new BalancedMtxReader(filename, comm);
    Epetra_CrsMatrix* y = reader->read();
    delete reader;
    const Epetra_Map& short_map = y->DomainMap();
    const Epetra_Map& long_map = y->RangeMap();

    // Y_transp matrix
    filename = "Y_transp.mtx";
    if (access(filename.c_str(), R_OK)) {
        cerr << "Unable to access Y_transp matrix file.\n";
        exit(1);
    }
    reader = new MtxReader(filename, short_map, long_map, short_map);
    Epetra_CrsMatrix* y_transp = reader->read();
    delete reader;

    // C matrix (with compatible maps)
    filename = "C.mtx";
    if (access(filename.c_str(), R_OK)) {
        cerr << "Unable to access C matrix file.\n";
        exit(1);
    }
    reader = new MtxReader(filename, long_map, short_map, long_map);
    Epetra_CrsMatrix* c = reader->read();
    delete reader;

    Epetra_CrsMatrix* ha = new Epetra_CrsMatrix (Copy, short_map, 0);
    EpetraExt::MatrixMatrix::Multiply(*y, true, *c, false, *ha);

    Epetra_CrsMatrix* hb = new Epetra_CrsMatrix (Copy, short_map, 0);
    EpetraExt::MatrixMatrix::Multiply(*y_transp, false, *c, false, *hb);

    export_Epetra_CrsMatrix(*ha, 5, "ha.mtx");
    export_Epetra_CrsMatrix(*hb, 5, "hb.mtx");

    delete y;
    delete y_transp;
    delete c;
    delete ha;
    delete hb;

    return 0;
}

---