src/matrixmarket/balancedmtxreader.cpp

Go to the documentation of this file.

/***************************************************************************
                          balancedmtxreader.cpp  -  description
                             -------------------
    begin                : Thu Dec 4 2003
    copyright            : (C) 2003 by Roman Geus
    email                : roman.geus@psi.ch
***************************************************************************/

/***************************************************************************
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 ***************************************************************************/

#include <cassert>
#include <vector>
#include "Epetra_config.h"
#ifdef HAVE_MPI
#include <mpi.h>
#endif
#include "Epetra_Map.h"
#include "Epetra_Comm.h"

#include "mtxdataiterator.h"
#include "balancedmtxreader.h"

using namespace std;

Epetra_CrsMatrix* BalancedMtxReader::read() {

    vector<int> start_row(_comm.NumProc() + 1);

    //
    // P0 analyses the entire matrix and determines the distribution
    //
    if (_comm.MyPID() == 0) {
        MtxDataIterator it(*_istr, true);
        vector<int> row_nnz;
        row_nnz.assign(it.get_rows(), 0);

        while (!it.eof()) {
            row_nnz[(*it).row] ++;
            ++it;
        }

        // compute prefix sum
        for(int i = 1; i < it.get_rows(); i ++)
            row_nnz[i] += row_nnz[i-1];

        start_row[0] = 0;
        int i = 0;
        for (int p = 1; p < _comm.NumProc(); ++ p) {
            while (_comm.NumProc() * row_nnz[i] < p * row_nnz[it.get_rows()-1])
                ++ i;
            start_row[p] = i;
        }
        start_row[_comm.NumProc()] = it.get_rows();

        for (int p = 0; p < _comm.NumProc(); p ++) {
            int nnz_local;
            if (p == 0)
                nnz_local = row_nnz[start_row[p+1] - 1];
            else
                nnz_local = row_nnz[start_row[p+1] - 1] - row_nnz[start_row[p] - 1];
            //cout << "processor " << p << " gets " << nnz_local << " non-zeros." << endl;
        }

        for (int p = 0; p < _comm.NumProc(); p ++)
            start_row[p] = start_row[p+1] - start_row[p];

        for (int p = 0; p < _comm.NumProc(); p ++) {
            //cout << "processor " << p << " gets " << start_row[p] << " lines." << endl;
        }

        // reset input stream, so the matrix can be read again on P0
        _istr->seekg(0); // does this always work?
    }

    // each processor receives the number of rows it will hold
    int my_num_rows;
#ifdef HAVE_MPI
    MPI_Scatter(&start_row[0], 1, MPI_INT, &my_num_rows, 1, MPI_INT, 0, MPI_COMM_WORLD);
#endif

    MtxDataIterator it(*_istr, true);

    // Allocate maps
    Epetra_Map range_row_map(it.get_rows(), my_num_rows, 0, _comm);

    // Allocate matrix and read matrix data
    Epetra_CrsMatrix* mat = new Epetra_CrsMatrix (Copy, range_row_map, 0);

    // read matrix entries
    while (!it.eof()) {
        if (range_row_map.MyGID((*it).row))
            mat->InsertGlobalValues((*it).row, 1, &(*it).val, &(*it).col);
        ++ it;
    }

    // finalise Crs matrix
    if (it.get_rows() == it.get_cols()) {
        // if square, make sure that row, domain and range map are identical
        mat->FillComplete(range_row_map, range_row_map);
    } else {
        // if not square, only row and range map are identical
        Epetra_Map domain_map(it.get_cols(), 0, _comm);
        mat->FillComplete(domain_map, range_row_map);
    }

    return mat;
}

---