src/export_Epetra_CrsMatrix.cpp

Go to the documentation of this file.
00001 /***************************************************************************
00002                           export_Epetra_CrsMatrix.cpp  -  description
00003                              -------------------
00004     begin                : Thu Feb 5 2004
00005     copyright            : (C) 2004 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 <algorithm>
00019 #include <iostream>
00020 #include <fstream>
00021 #include <vector>
00022 
00023 #include "Epetra_Comm.h"
00024 #include "Epetra_Export.h"
00025 #include "Epetra_Map.h"
00026 #include "Epetra_CrsMatrix.h"
00027 
00030 class colval {
00031 public:
00032     colval() {}
00033     colval(int col, double val)
00034         : col_(col),
00035           val_(val)
00036     { }
00037     bool operator<(const colval& rhs) const {
00038         return col_ < rhs.col_;
00039     }
00040     const int& col() const {
00041         return col_;
00042     }
00043     const double& val() const {
00044         return val_;
00045     }
00046 private:
00047     int col_;
00048     double val_;
00049 };
00050 
00051 void export_Epetra_CrsMatrix(const Epetra_CrsMatrix& mat, int precision, std::ostream& ostr) {
00052     int info;
00053     const Epetra_Map& RowMap = mat.RowMap();
00054     const Epetra_Map& ColMap = mat.ColMap();
00055     const Epetra_Comm& Comm = RowMap.Comm();
00056 
00057     assert(RowMap.LinearMap() || Comm.NumProc() == 1);
00058 
00059     int oldPrecision = ostr.precision();
00060     std::ostream::fmtflags oldFlags = ostr.setf(ostr.scientific);
00061     ostr.precision(precision);
00062 
00063     ostr << "%%MatrixMarket matrix coordinate real general\n";
00064     ostr << mat.NumGlobalRows() << ' ' << mat.NumGlobalCols()
00065          << ' ' << mat.NumGlobalNonzeros() << std::endl;
00066 
00067     for (int i = 0; i < mat.NumMyRows(); i ++) {
00068         int nof_entries;
00069         int* col;
00070         double* val;
00071         info = mat.ExtractMyRowView(i, nof_entries, val, col);
00072         assert(info == 0);
00073 
00074         std::vector<colval> sparse_row(nof_entries);
00075         for (int k = 0; k < nof_entries; k ++)
00076             sparse_row[k] = colval(ColMap.GID(col[k]), val[k]);
00077         std::sort(sparse_row.begin(), sparse_row.end());
00078 
00079         for (int k = 0; k < nof_entries; k ++)
00080             ostr << RowMap.GID(i)+1 << ' ' << sparse_row[k].col()+1 <<
00081                 ' ' << sparse_row[k].val() << std::endl;
00082     }
00083 
00084     // restore old format state
00085     ostr.flags(oldFlags);
00086     ostr.precision(oldPrecision);
00087 }
00088 
00089 void export_Epetra_CrsMatrix(const Epetra_CrsMatrix& mat_in, int precision, const std::string& file_name) {
00090 
00091     int info;
00092     const Epetra_Comm& Comm = mat_in.Comm();
00093     const Epetra_CrsMatrix* matp(0);
00094     Epetra_CrsMatrix* mat_new(0);
00095 
00096     // Redistribute matrix using a linear map if necessary
00097     if (~mat_in.RowMap().LinearMap() && Comm.NumProc() != 1) {
00098         Epetra_Map linear_map(mat_in.NumGlobalRows(), 0, Comm);
00099         Epetra_Export exporter(mat_in.RowMap(), linear_map);
00100         mat_new = new Epetra_CrsMatrix(Copy, linear_map, 0);
00101         mat_new->Export(mat_in, exporter, Insert);
00102         mat_new->FillComplete(mat_in.DomainMap(), linear_map);
00103         matp = mat_new;
00104     } else
00105         matp = &mat_in;    
00106     const Epetra_CrsMatrix& mat(*matp);
00107     
00108     const Epetra_Map& RowMap = mat.RowMap();
00109     const Epetra_Map& ColMap = mat.ColMap();
00110     assert(RowMap.LinearMap() || Comm.NumProc() == 1);
00111 
00112     std::ofstream file;
00113     file.setf(file.scientific);
00114     file.precision(precision);
00115 
00116     // P0 writes header
00117     if (Comm.MyPID() == 0) {
00118         file.open(file_name.c_str(), std::ios_base::out);
00119         file << "%%MatrixMarket matrix coordinate real general\n";
00120         file << mat.NumGlobalRows() << ' ' << mat.NumGlobalCols()
00121              << ' ' << mat.NumGlobalNonzeros() << std::endl;
00122         file.close();
00123     }
00124 
00125     // Processors append matrix data one after the other
00126     for (int p = 0; p < Comm.NumProc(); p ++) {
00127         if (p == Comm.MyPID()) {
00128             file.open(file_name.c_str(), std::ios_base::app);
00129             for (int i = 0; i < mat.NumMyRows(); i ++) {
00130                 int nof_entries;
00131                 int* col;
00132                 double* val;
00133                 info = mat.ExtractMyRowView(i, nof_entries, val, col);
00134                 assert(info == 0);
00135 
00136                 std::vector<colval> sparse_row(nof_entries);
00137                 for (int k = 0; k < nof_entries; k ++)
00138                     sparse_row[k] = colval(ColMap.GID(col[k]), val[k]);
00139                 std::sort(sparse_row.begin(), sparse_row.end());
00140 
00141                 for (int k = 0; k < nof_entries; k ++)
00142                     file << RowMap.GID(i)+1 << ' ' << sparse_row[k].col()+1 <<
00143                         ' ' << sparse_row[k].val() << std::endl;
00144             }
00145             file.close();
00146         }
00147         Comm.Barrier();
00148     }
00149 
00150     delete mat_new;
00151 }
00152 
00153 void 
00154 export_Epetra_Map(const Epetra_Map& map,
00155                   std::ostream& ostr)
00156 {
00157     const Epetra_Comm& comm(map.Comm());
00158     const int mpi_size(comm.NumProc());
00159     const int mpi_rank(comm.MyPID());
00160        
00161     for (int i = 0; i < mpi_size; ++ i) {
00162         if (i == mpi_rank)
00163             ostr << map;
00164         comm.Barrier();
00165     }
00166 }
00167 
00168 void 
00169 export_Epetra_Map(const Epetra_Map& map,
00170                   const std::string& file_name)
00171 {
00172     const Epetra_Comm& comm(map.Comm());
00173     const int mpi_size(comm.NumProc());
00174     const int mpi_rank(comm.MyPID());
00175     std::ofstream file;
00176 
00177     // P0 creates empty file
00178     if (mpi_rank == 0) {
00179         file.open(file_name.c_str(), std::ios_base::out);
00180         file.close();
00181     }
00182 
00183     // Processors append map data one after the other
00184     for (int p = 0; p < mpi_size; p ++) {
00185         if (p == mpi_rank) {
00186             file.open(file_name.c_str(), std::ios_base::app);
00187             file << map;
00188             file.close();
00189         }
00190         comm.Barrier();
00191     }
00192 }

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