EpetraExt Package Browser (Single Doxygen Collection) Development
EpetraExt_MMHelpers.cpp
Go to the documentation of this file.
00001 //@HEADER
00002 // ***********************************************************************
00003 //
00004 //     EpetraExt: Epetra Extended - Linear Algebra Services Package
00005 //                 Copyright (2011) Sandia Corporation
00006 //
00007 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00008 // the U.S. Government retains certain rights in this software.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00038 //
00039 // ***********************************************************************
00040 //@HEADER
00041 
00042 #include <EpetraExt_ConfigDefs.h>
00043 #include <EpetraExt_MMHelpers.h>
00044 #include <Epetra_Comm.h>
00045 #include <Epetra_Map.h>
00046 #include <Epetra_CrsMatrix.h>
00047 
00048 namespace EpetraExt {
00049 
00050 CrsMatrixStruct::CrsMatrixStruct()
00051  : numRows(0), numEntriesPerRow(NULL), indices(NULL), values(NULL),
00052       remote(NULL), numRemote(0), rowMap(NULL), colMap(NULL),
00053       domainMap(NULL), importColMap(NULL), importMatrix(NULL)
00054 {
00055 }
00056 
00057 CrsMatrixStruct::~CrsMatrixStruct()
00058 {
00059   deleteContents();
00060 }
00061 
00062 void CrsMatrixStruct::deleteContents()
00063 {
00064   numRows = 0;
00065   delete [] numEntriesPerRow; numEntriesPerRow = NULL;
00066   delete [] indices; indices = NULL;
00067   delete [] values; values = NULL;
00068   delete [] remote; remote = NULL;
00069   numRemote = 0;
00070   delete importMatrix;
00071 }
00072 
00073 int dumpCrsMatrixStruct(const CrsMatrixStruct& M)
00074 {
00075   cout << "proc " << M.rowMap->Comm().MyPID()<<endl;
00076   cout << "numRows: " << M.numRows<<endl;
00077   for(int i=0; i<M.numRows; ++i) {
00078     for(int j=0; j<M.numEntriesPerRow[i]; ++j) {
00079       if (M.remote[i]) {
00080         cout << "  *"<<M.rowMap->GID(i)<<"   "
00081              <<M.importColMap->GID(M.indices[i][j])<<"   "<<M.values[i][j]<<endl;
00082       }
00083       else {
00084         cout << "   "<<M.rowMap->GID(i)<<"   "
00085              <<M.colMap->GID(M.indices[i][j])<<"   "<<M.values[i][j]<<endl;
00086       }
00087     }
00088   }
00089   return(0);
00090 }
00091 
00092 CrsWrapper_Epetra_CrsMatrix::CrsWrapper_Epetra_CrsMatrix(Epetra_CrsMatrix& epetracrsmatrix)
00093  : ecrsmat_(epetracrsmatrix)
00094 {
00095 }
00096 
00097 CrsWrapper_Epetra_CrsMatrix::~CrsWrapper_Epetra_CrsMatrix()
00098 {
00099 }
00100 
00101 const Epetra_Map&
00102 CrsWrapper_Epetra_CrsMatrix::RowMap() const
00103 {
00104   return ecrsmat_.RowMap();
00105 }
00106 
00107 bool CrsWrapper_Epetra_CrsMatrix::Filled()
00108 {
00109   return ecrsmat_.Filled();
00110 }
00111 
00112 int
00113 CrsWrapper_Epetra_CrsMatrix::InsertGlobalValues(int GlobalRow, int NumEntries, double* Values, int* Indices)
00114 {
00115   return ecrsmat_.InsertGlobalValues(GlobalRow, NumEntries, Values, Indices);
00116 }
00117 
00118 int
00119 CrsWrapper_Epetra_CrsMatrix::SumIntoGlobalValues(int GlobalRow, int NumEntries, double* Values, int* Indices)
00120 {
00121   return ecrsmat_.SumIntoGlobalValues(GlobalRow, NumEntries, Values, Indices);
00122 }
00123 
00124 
00125 //------------------------------------
00126 
00127 CrsWrapper_GraphBuilder::CrsWrapper_GraphBuilder(const Epetra_Map& emap)
00128  : graph_(),
00129    rowmap_(emap),
00130    max_row_length_(0)
00131 {
00132   int num_rows = emap.NumMyElements();
00133   int* rows = emap.MyGlobalElements();
00134 
00135   for(int i=0; i<num_rows; ++i) {
00136     graph_[rows[i]] = new std::set<int>;
00137   }
00138 }
00139 
00140 CrsWrapper_GraphBuilder::~CrsWrapper_GraphBuilder()
00141 {
00142   std::map<int,std::set<int>*>::iterator
00143     iter = graph_.begin(), iter_end = graph_.end();
00144   for(; iter!=iter_end; ++iter) {
00145     delete iter->second;
00146   }
00147 
00148   graph_.clear();
00149 }
00150 
00151 bool CrsWrapper_GraphBuilder::Filled()
00152 {
00153   return false;
00154 }
00155 
00156 int
00157 CrsWrapper_GraphBuilder::InsertGlobalValues(int GlobalRow, int NumEntries, double* Values, int* Indices)
00158 {
00159   std::map<int,std::set<int>*>::iterator
00160     iter = graph_.find(GlobalRow);
00161 
00162   if (iter == graph_.end()) return(-1);
00163 
00164   std::set<int>& cols = *(iter->second);
00165 
00166   for(int i=0; i<NumEntries; ++i) {
00167     cols.insert(Indices[i]);
00168   }
00169 
00170   int row_length = cols.size();
00171   if (row_length > max_row_length_) max_row_length_ = row_length;
00172 
00173   return(0);
00174 }
00175 
00176 int
00177 CrsWrapper_GraphBuilder::SumIntoGlobalValues(int GlobalRow, int NumEntries, double* Values, int* Indices)
00178 {
00179   return InsertGlobalValues(GlobalRow, NumEntries, Values, Indices);
00180 }
00181 
00182 std::map<int,std::set<int>*>&
00183 CrsWrapper_GraphBuilder::get_graph()
00184 {
00185   return graph_;
00186 }
00187 
00188 void insert_matrix_locations(CrsWrapper_GraphBuilder& graphbuilder,
00189                               Epetra_CrsMatrix& C)
00190 {
00191   int max_row_length = graphbuilder.get_max_row_length();
00192   if (max_row_length < 1) return;
00193 
00194   std::vector<int> indices(max_row_length);
00195   int* indices_ptr = &indices[0];
00196   std::vector<double> zeros(max_row_length, 0.0);
00197   double* zeros_ptr = &zeros[0];
00198 
00199   std::map<int,std::set<int>*>& graph = graphbuilder.get_graph();
00200 
00201   std::map<int,std::set<int>*>::iterator
00202     iter = graph.begin(), iter_end = graph.end();
00203 
00204   for(; iter!=iter_end; ++iter) {
00205     int row = iter->first;
00206     std::set<int>& cols = *(iter->second);
00207     int num_entries = cols.size();
00208 
00209     std::set<int>::iterator
00210       col_iter = cols.begin(), col_end = cols.end();
00211     for(int j=0; col_iter!=col_end; ++col_iter, ++j) {
00212       indices_ptr[j] = *col_iter;
00213     }
00214 
00215     C.InsertGlobalValues(row, num_entries, zeros_ptr, indices_ptr);
00216   }
00217 }
00218 
00219 void pack_outgoing_rows(const Epetra_CrsMatrix& mtx,
00220                         const std::vector<int>& proc_col_ranges,
00221                         std::vector<int>& send_rows,
00222                         std::vector<int>& rows_per_send_proc)
00223 {
00224   const Epetra_Map& rowmap = mtx.RowMap();
00225   int numrows = mtx.NumMyRows();
00226   const Epetra_CrsGraph& graph = mtx.Graph();
00227   int rowlen = 0;
00228   int* col_indices = NULL;
00229   int num_col_ranges = proc_col_ranges.size()/2;
00230   rows_per_send_proc.resize(num_col_ranges);
00231   send_rows.clear();
00232   for(int nc=0; nc<num_col_ranges; ++nc) {
00233     int first_col = proc_col_ranges[nc*2];
00234     int last_col = proc_col_ranges[nc*2+1];
00235     int num_send_rows = 0;
00236     for(int i=0; i<numrows; ++i) {
00237       int grow = rowmap.GID(i);
00238       if (mtx.Filled()) {
00239         const Epetra_Map& colmap = mtx.ColMap();
00240         graph.ExtractMyRowView(i, rowlen, col_indices);
00241         for(int j=0; j<rowlen; ++j) {
00242           int col = colmap.GID(col_indices[j]);
00243           if (first_col <= col && last_col >= col) {
00244             ++num_send_rows;
00245             send_rows.push_back(grow);
00246             break;
00247           }
00248         }
00249       }
00250       else {
00251         graph.ExtractGlobalRowView(grow, rowlen, col_indices);
00252         for(int j=0; j<rowlen; ++j) {
00253           if (first_col <= col_indices[j] && last_col >= col_indices[j]) {
00254             ++num_send_rows;
00255             send_rows.push_back(grow);
00256             break;
00257           }
00258         }
00259       }
00260     }
00261     rows_per_send_proc[nc] = num_send_rows;
00262   }
00263 }
00264 
00265 std::pair<int,int> get_col_range(const Epetra_CrsMatrix& mtx)
00266 {
00267   std::pair<int,int> col_range;
00268   if (mtx.Filled()) {
00269     col_range = get_col_range(mtx.ColMap());
00270   }
00271   else {
00272     const Epetra_Map& row_map = mtx.RowMap();
00273     col_range.first = row_map.MaxMyGID();
00274     col_range.second = row_map.MinMyGID();
00275     int rowlen = 0;
00276     int* col_indices = NULL;
00277     const Epetra_CrsGraph& graph = mtx.Graph();
00278     for(int i=0; i<row_map.NumMyElements(); ++i) {
00279       graph.ExtractGlobalRowView(row_map.GID(i), rowlen, col_indices);
00280       for(int j=0; j<rowlen; ++j) {
00281         if (col_indices[j] < col_range.first) col_range.first = col_indices[j];
00282         if (col_indices[j] > col_range.second) col_range.second = col_indices[j];
00283       }
00284     }
00285   }
00286 
00287   return col_range;
00288 }
00289 
00290 }//namespace EpetraExt
00291 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines