src/export_Epetra_CrsMatrix.cpp

Go to the documentation of this file.

/***************************************************************************
                          export_Epetra_CrsMatrix.cpp  -  description
                             -------------------
    begin                : Thu Feb 5 2004
    copyright            : (C) 2004 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 <algorithm>
#include <iostream>
#include <fstream>
#include <vector>

#include "Epetra_Comm.h"
#include "Epetra_Export.h"
#include "Epetra_Map.h"
#include "Epetra_CrsMatrix.h"

class colval {
public:
    colval() {}
    colval(int col, double val)
        : col_(col),
          val_(val)
    { }
    bool operator<(const colval& rhs) const {
        return col_ < rhs.col_;
    }
    const int& col() const {
        return col_;
    }
    const double& val() const {
        return val_;
    }
private:
    int col_;
    double val_;
};

void export_Epetra_CrsMatrix(const Epetra_CrsMatrix& mat, int precision, std::ostream& ostr) {
    int info;
    const Epetra_Map& RowMap = mat.RowMap();
    const Epetra_Map& ColMap = mat.ColMap();
    const Epetra_Comm& Comm = RowMap.Comm();

    assert(RowMap.LinearMap() || Comm.NumProc() == 1);

    int oldPrecision = ostr.precision();
    std::ostream::fmtflags oldFlags = ostr.setf(ostr.scientific);
    ostr.precision(precision);

    ostr << "%%MatrixMarket matrix coordinate real general\n";
    ostr << mat.NumGlobalRows() << ' ' << mat.NumGlobalCols()
         << ' ' << mat.NumGlobalNonzeros() << std::endl;

    for (int i = 0; i < mat.NumMyRows(); i ++) {
        int nof_entries;
        int* col;
        double* val;
        info = mat.ExtractMyRowView(i, nof_entries, val, col);
        assert(info == 0);

        std::vector<colval> sparse_row(nof_entries);
        for (int k = 0; k < nof_entries; k ++)
            sparse_row[k] = colval(ColMap.GID(col[k]), val[k]);
        std::sort(sparse_row.begin(), sparse_row.end());

        for (int k = 0; k < nof_entries; k ++)
            ostr << RowMap.GID(i)+1 << ' ' << sparse_row[k].col()+1 <<
                ' ' << sparse_row[k].val() << std::endl;
    }

    // restore old format state
    ostr.flags(oldFlags);
    ostr.precision(oldPrecision);
}

void export_Epetra_CrsMatrix(const Epetra_CrsMatrix& mat_in, int precision, const std::string& file_name) {

    int info;
    const Epetra_Comm& Comm = mat_in.Comm();
    const Epetra_CrsMatrix* matp(0);
    Epetra_CrsMatrix* mat_new(0);

    // Redistribute matrix using a linear map if necessary
    if (~mat_in.RowMap().LinearMap() && Comm.NumProc() != 1) {
        Epetra_Map linear_map(mat_in.NumGlobalRows(), 0, Comm);
        Epetra_Export exporter(mat_in.RowMap(), linear_map);
        mat_new = new Epetra_CrsMatrix(Copy, linear_map, 0);
        mat_new->Export(mat_in, exporter, Insert);
        mat_new->FillComplete(mat_in.DomainMap(), linear_map);
        matp = mat_new;
    } else
        matp = &mat_in;    
    const Epetra_CrsMatrix& mat(*matp);
    
    const Epetra_Map& RowMap = mat.RowMap();
    const Epetra_Map& ColMap = mat.ColMap();
    assert(RowMap.LinearMap() || Comm.NumProc() == 1);

    std::ofstream file;
    file.setf(file.scientific);
    file.precision(precision);

    // P0 writes header
    if (Comm.MyPID() == 0) {
        file.open(file_name.c_str(), std::ios_base::out);
        file << "%%MatrixMarket matrix coordinate real general\n";
        file << mat.NumGlobalRows() << ' ' << mat.NumGlobalCols()
             << ' ' << mat.NumGlobalNonzeros() << std::endl;
        file.close();
    }

    // Processors append matrix data one after the other
    for (int p = 0; p < Comm.NumProc(); p ++) {
        if (p == Comm.MyPID()) {
            file.open(file_name.c_str(), std::ios_base::app);
            for (int i = 0; i < mat.NumMyRows(); i ++) {
                int nof_entries;
                int* col;
                double* val;
                info = mat.ExtractMyRowView(i, nof_entries, val, col);
                assert(info == 0);

                std::vector<colval> sparse_row(nof_entries);
                for (int k = 0; k < nof_entries; k ++)
                    sparse_row[k] = colval(ColMap.GID(col[k]), val[k]);
                std::sort(sparse_row.begin(), sparse_row.end());

                for (int k = 0; k < nof_entries; k ++)
                    file << RowMap.GID(i)+1 << ' ' << sparse_row[k].col()+1 <<
                        ' ' << sparse_row[k].val() << std::endl;
            }
            file.close();
        }
        Comm.Barrier();
    }

    delete mat_new;
}

void 
export_Epetra_Map(const Epetra_Map& map,
                  std::ostream& ostr)
{
    const Epetra_Comm& comm(map.Comm());
    const int mpi_size(comm.NumProc());
    const int mpi_rank(comm.MyPID());
       
    for (int i = 0; i < mpi_size; ++ i) {
        if (i == mpi_rank)
            ostr << map;
        comm.Barrier();
    }
}

void 
export_Epetra_Map(const Epetra_Map& map,
                  const std::string& file_name)
{
    const Epetra_Comm& comm(map.Comm());
    const int mpi_size(comm.NumProc());
    const int mpi_rank(comm.MyPID());
    std::ofstream file;

    // P0 creates empty file
    if (mpi_rank == 0) {
        file.open(file_name.c_str(), std::ios_base::out);
        file.close();
    }

    // Processors append map data one after the other
    for (int p = 0; p < mpi_size; p ++) {
        if (p == mpi_rank) {
            file.open(file_name.c_str(), std::ios_base::app);
            file << map;
            file.close();
        }
        comm.Barrier();
    }
}

---