Epetra Package Browser (Single Doxygen Collection) Development
Epetra_matrix_data.cpp
Go to the documentation of this file.
00001 
00002 //@HEADER
00003 // ************************************************************************
00004 // 
00005 //               Epetra: Linear Algebra Services Package 
00006 //                 Copyright 2001 Sandia Corporation
00007 // 
00008 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00009 // the U.S. Government retains certain rights in this software.
00010 //
00011 // Redistribution and use in source and binary forms, with or without
00012 // modification, are permitted provided that the following conditions are
00013 // met:
00014 //
00015 // 1. Redistributions of source code must retain the above copyright
00016 // notice, this list of conditions and the following disclaimer.
00017 //
00018 // 2. Redistributions in binary form must reproduce the above copyright
00019 // notice, this list of conditions and the following disclaimer in the
00020 // documentation and/or other materials provided with the distribution.
00021 //
00022 // 3. Neither the name of the Corporation nor the names of the
00023 // contributors may be used to endorse or promote products derived from
00024 // this software without specific prior written permission.
00025 //
00026 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00027 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00028 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00029 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00030 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00031 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00032 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00033 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00034 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00035 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00036 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00037 //
00038 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00039 // 
00040 // ************************************************************************
00041 //@HEADER
00042 
00043 #include <Epetra_matrix_data.h>
00044 #include <Epetra_Map.h>
00045 #include <Epetra_CrsMatrix.h>
00046 #include <Epetra_Util.h>
00047 
00048 namespace epetra_test {
00049 
00050 matrix_data::matrix_data(int num_rows,
00051                          int* rowlengths,
00052                          int blocksize)
00053  : numrows_(num_rows),
00054    numcols_(0),
00055    rows_(0),
00056    rowlengths_(0),
00057    blocksize_(blocksize),
00058    colindices_(0),
00059    coefs_(0)
00060 {
00061   if (numrows_ > 0) {
00062     rows_ = new int[numrows_];
00063     rowlengths_ = new int[numrows_];
00064     colindices_ = new int*[numrows_];
00065     coefs_ = new double*[numrows_];
00066     int dim = blocksize*blocksize;
00067     for(int i=0; i<numrows_; ++i) {
00068       rows_[i] = i;
00069       rowlengths_[i] = rowlengths[i];
00070       colindices_[i] = new int[rowlengths_[i]];
00071       coefs_[i] = new double[rowlengths_[i]*dim];
00072 
00073       for(int j=0; j<rowlengths_[i]; ++j) {
00074         colindices_[i][j] = 0;
00075         for(int k=0; k<dim; ++k) coefs_[i][j*dim+k] = 0.0;
00076       }
00077     }
00078   }
00079 }
00080 
00081 matrix_data::matrix_data(int num_rows,
00082                          int num_cols,
00083                          int num_off_diagonals,
00084                          int blocksize)
00085  : numrows_(num_rows),
00086    numcols_(num_cols),
00087    rows_(0),
00088    rowlengths_(0),
00089    blocksize_(blocksize),
00090    colindices_(0),
00091    coefs_(0)
00092 {
00093   if (numrows_ > 0) {
00094     rows_       = new int[numrows_];
00095     rowlengths_ = new int[numrows_];
00096     colindices_ = new int*[numrows_];
00097     coefs_      = new double*[numrows_];
00098 
00099     int max_row_length = 1+num_off_diagonals*2;
00100 
00101     for(int i=0; i<numrows_; ++i) {
00102       rows_[i] = i;
00103       if (i >= num_off_diagonals && numrows_-i > num_off_diagonals) {
00104         rowlengths_[i] = max_row_length;
00105       }
00106       else {
00107         if (i<num_off_diagonals) {
00108           rowlengths_[i] = 1+max_row_length/2+i;
00109         }
00110         else {
00111           rowlengths_[i] = 1+max_row_length/2+numrows_-i-1;
00112         }
00113       }
00114       colindices_[i] = new int[rowlengths_[i]];
00115       int dim = blocksize*blocksize;
00116       coefs_[i] = new double[rowlengths_[i]*dim];
00117 
00118       int first_col = i - max_row_length/2;
00119       if (first_col < 0) first_col = 0;
00120 
00121       for(int j=0; j<rowlengths_[i]; ++j) {
00122         colindices_[i][j] = first_col+j;
00123         for(int k=0; k<dim; ++k) {
00124           coefs_[i][j*dim+k] = 1.0;
00125         }
00126       }
00127     }
00128   }
00129 }
00130 
00131 static const int nodes_per_elem = 4;
00132 
00133 void get_node_ids(int elem_id, int* node_ids)
00134 {
00135   int first_node = 2*elem_id;
00136   for(int i=0; i<nodes_per_elem; ++i) node_ids[i] = first_node+i;
00137 }
00138 
00139 matrix_data::matrix_data(int num_quad_elements,
00140                          int num_dof_per_node,
00141        bool make_numerically_nonsymmetric)
00142  : numrows_(0),
00143    numcols_(0),
00144    rows_(0),
00145    rowlengths_(0),
00146    blocksize_(num_dof_per_node),
00147    colindices_(0),
00148    coefs_(0)
00149 {
00150   //Set up matrix-data representing a simple finite-element
00151   //mesh containing 2-D quad elements
00152   //
00153   //   *-----*-----*-----*
00154   //  0|    2|    4|    6|
00155   //   | 0   | 1   | ne-1|
00156   //   |     |     |     |
00157   //   *-----*-----*-----*
00158   //  1     3     5     7
00159   //
00160   //In the above drawing, 'ne' means num-elems. node-numbers are to the
00161   //lower-left of each node (*).
00162 
00163   numrows_ = num_quad_elements*2+2;
00164 
00165   if (numrows_ > 0) {
00166     rows_       = new int[numrows_];
00167     rowlengths_ = new int[numrows_];
00168     colindices_ = new int*[numrows_];
00169     coefs_      = new double*[numrows_];
00170 
00171     int i, j, k;
00172     for(i=0; i<numrows_; ++i) {
00173       rows_[i] = i;
00174       rowlengths_[i] = 0;
00175     }
00176 
00177     int* nodes = new int[nodes_per_elem];
00178     for(i=0; i<num_quad_elements; ++i) {
00179       get_node_ids(i, nodes);
00180 
00181       for(j=0; j<nodes_per_elem; ++j) {
00182         int node_j = nodes[j];
00183         for(k=0; k<nodes_per_elem; ++k) {
00184           int insertPoint = -1;
00185           int alloclen = rowlengths_[node_j];
00186           int offset = Epetra_Util_binary_search(nodes[k], colindices_[node_j],
00187                                                  rowlengths_[node_j], insertPoint);
00188           if (offset < 0) {
00189             Epetra_Util_insert(nodes[k], insertPoint,
00190                                colindices_[node_j], rowlengths_[node_j],
00191                                alloclen);
00192           }
00193         }
00194       }
00195     }
00196 
00197     int dim = blocksize_*blocksize_;
00198     for(i=0; i<numrows_; ++i) {
00199       int len = rowlengths_[i]*dim;
00200       coefs_[i] = new double[len];
00201       for(j=0; j<len; ++j) {
00202   if (make_numerically_nonsymmetric) {
00203     coefs_[i][j] = 1.0*(j+1);
00204   }
00205   else {
00206     coefs_[i][j] = 1.0;
00207   }
00208       }
00209     }
00210   }
00211 }
00212 
00213 matrix_data::~matrix_data()
00214 {
00215   for(int i=0; i<numrows_; ++i) {
00216     delete [] colindices_[i];
00217     delete [] coefs_[i];
00218   }
00219 
00220   delete [] colindices_; colindices_ = 0;
00221   delete [] coefs_; coefs_ = 0;
00222   delete [] rowlengths_; rowlengths_ = 0;
00223   delete [] rows_; rows_ = 0;
00224   numrows_ = 0;
00225 }
00226 
00227 double* matrix_data::coefs(int row, int col)
00228 {
00229   int insertPoint = -1;
00230   int row_idx = Epetra_Util_binary_search(row, rows_, numrows_,
00231                                           insertPoint);
00232   if (row_idx < 0) {
00233     cerr << "ERROR, row " << row
00234          << " not found in matrix_data"<<endl;
00235     return 0;
00236   }
00237 
00238   int col_idx = Epetra_Util_binary_search(col, colindices_[row_idx],
00239                                           rowlengths_[row_idx], insertPoint);
00240   if (col_idx < 0) {
00241     cerr << "ERROR, col " << col
00242          << " not found in matrix_data"<<endl;
00243     return 0;
00244   }
00245 
00246   int dim = blocksize_*blocksize_;
00247   return( &(coefs_[row_idx][col_idx*dim]) );
00248 }
00249 
00250 void matrix_data::copy_local_data_to_matrix(Epetra_CrsMatrix& A)
00251 {
00252   const Epetra_Map& rowmap = A.RowMap();
00253 
00254   for(int i=0; i<numrows_; ++i) {
00255     int row = rows_[i];
00256     if (rowmap.MyGID(row)) {
00257       int err = A.ReplaceGlobalValues(row, rowlengths_[i],
00258               coefs_[i], colindices_[i]);
00259       if (err < 0) {
00260   err = A.InsertGlobalValues(row, rowlengths_[i],
00261            coefs_[i], colindices_[i]);
00262       }
00263     }
00264   }
00265 }
00266 
00267 bool matrix_data::compare_local_data(const Epetra_CrsMatrix& A)
00268 {
00269   const Epetra_Map& map = A.RowMap();
00270   int numMyRows = map.NumMyElements();
00271   int* myRows = map.MyGlobalElements();
00272 
00273   Epetra_Util util;
00274 
00275   for(int i=0; i<numMyRows; ++i) {
00276     int row = myRows[i];
00277     int rowLen = A.NumGlobalEntries(row);
00278     if (rowLen != rowlengths_[row]) {
00279       return(false);
00280     }
00281 
00282     int* indices = new int[rowLen];
00283     double* values = new double[rowLen];
00284     A.ExtractGlobalRowCopy(row, rowLen, rowLen, values, indices);
00285 
00286     util.Sort(true, rowLen, indices, 1, &values, 0, 0);
00287 
00288     bool same = true;
00289     int* this_indices = colindices_[row];
00290     double* this_coefs = coefs_[row];
00291 
00292     for(int j=0; j<rowLen; ++j) {
00293       if (indices[j] != this_indices[j]) {
00294         same = false; break;
00295       }
00296       if (values[j] != this_coefs[j]) {
00297         same = false; break;
00298       }
00299     }
00300 
00301     delete [] indices;
00302     delete [] values;
00303 
00304     if (!same) return(false);
00305   }
00306 
00307   return(true);
00308 }
00309 
00310 }//namespace epetra_test
00311 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines