src/matrixmarket/balancedmtxreader.cpp

Go to the documentation of this file.
00001 /***************************************************************************
00002                           balancedmtxreader.cpp  -  description
00003                              -------------------
00004     begin                : Thu Dec 4 2003
00005     copyright            : (C) 2003 by Roman Geus
00006     email                : roman.geus@psi.ch
00007 ***************************************************************************/
00008 
00009 /***************************************************************************
00010  *                                                                         *
00011  *   This program is free software; you can redistribute it and/or modify  *
00012  *   it under the terms of the GNU General Public License as published by  *
00013  *   the Free Software Foundation; either version 2 of the License, or     *
00014  *   (at your option) any later version.                                   *
00015  *                                                                         *
00016  ***************************************************************************/
00017 
00018 #include <cassert>
00019 #include <vector>
00020 #include "Epetra_config.h"
00021 #ifdef HAVE_MPI
00022 #include <mpi.h>
00023 #endif
00024 #include "Epetra_Map.h"
00025 #include "Epetra_Comm.h"
00026 
00027 #include "mtxdataiterator.h"
00028 #include "balancedmtxreader.h"
00029 
00030 using namespace std;
00031 
00032 Epetra_CrsMatrix* BalancedMtxReader::read() {
00033 
00034     vector<int> start_row(_comm.NumProc() + 1);
00035 
00036     //
00037     // P0 analyses the entire matrix and determines the distribution
00038     //
00039     if (_comm.MyPID() == 0) {
00040         MtxDataIterator it(*_istr, true);
00041         vector<int> row_nnz;
00042         row_nnz.assign(it.get_rows(), 0);
00043 
00044         while (!it.eof()) {
00045             row_nnz[(*it).row] ++;
00046             ++it;
00047         }
00048 
00049         // compute prefix sum
00050         for(int i = 1; i < it.get_rows(); i ++)
00051             row_nnz[i] += row_nnz[i-1];
00052 
00053         start_row[0] = 0;
00054         int i = 0;
00055         for (int p = 1; p < _comm.NumProc(); ++ p) {
00056             while (_comm.NumProc() * row_nnz[i] < p * row_nnz[it.get_rows()-1])
00057                 ++ i;
00058             start_row[p] = i;
00059         }
00060         start_row[_comm.NumProc()] = it.get_rows();
00061 
00062         for (int p = 0; p < _comm.NumProc(); p ++) {
00063             int nnz_local;
00064             if (p == 0)
00065                 nnz_local = row_nnz[start_row[p+1] - 1];
00066             else
00067                 nnz_local = row_nnz[start_row[p+1] - 1] - row_nnz[start_row[p] - 1];
00068             //cout << "processor " << p << " gets " << nnz_local << " non-zeros." << endl;
00069         }
00070 
00071         for (int p = 0; p < _comm.NumProc(); p ++)
00072             start_row[p] = start_row[p+1] - start_row[p];
00073 
00074         for (int p = 0; p < _comm.NumProc(); p ++) {
00075             //cout << "processor " << p << " gets " << start_row[p] << " lines." << endl;
00076         }
00077 
00078         // reset input stream, so the matrix can be read again on P0
00079         _istr->seekg(0); // does this always work?
00080     }
00081 
00082     // each processor receives the number of rows it will hold
00083     int my_num_rows;
00084 #ifdef HAVE_MPI
00085     MPI_Scatter(&start_row[0], 1, MPI_INT, &my_num_rows, 1, MPI_INT, 0, MPI_COMM_WORLD);
00086 #endif
00087 
00088     MtxDataIterator it(*_istr, true);
00089 
00090     // Allocate maps
00091     Epetra_Map range_row_map(it.get_rows(), my_num_rows, 0, _comm);
00092 
00093     // Allocate matrix and read matrix data
00094     Epetra_CrsMatrix* mat = new Epetra_CrsMatrix (Copy, range_row_map, 0);
00095 
00096     // read matrix entries
00097     while (!it.eof()) {
00098         if (range_row_map.MyGID((*it).row))
00099             mat->InsertGlobalValues((*it).row, 1, &(*it).val, &(*it).col);
00100         ++ it;
00101     }
00102 
00103     // finalise Crs matrix
00104     if (it.get_rows() == it.get_cols()) {
00105         // if square, make sure that row, domain and range map are identical
00106         mat->FillComplete(range_row_map, range_row_map);
00107     } else {
00108         // if not square, only row and range map are identical
00109         Epetra_Map domain_map(it.get_cols(), 0, _comm);
00110         mat->FillComplete(domain_map, range_row_map);
00111     }
00112 
00113     return mat;
00114 }

Generated on Fri Oct 26 13:35:12 2007 for FEMAXX (Finite Element Maxwell Eigensolver) by  doxygen 1.4.7