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 terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00009 // license for use of this work by or on behalf of the U.S. Government.
00010 //
00011 // This library is free software; you can redistribute it and/or modify
00012 // it under the terms of the GNU Lesser General Public License as
00013 // published by the Free Software Foundation; either version 2.1 of the
00014 // License, or (at your option) any later version.
00015 //
00016 // This library is distributed in the hope that it will be useful, but
00017 // WITHOUT ANY WARRANTY; without even the implied warranty of
00018 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00019 // Lesser General Public License for more details.
00020 //
00021 // You should have received a copy of the GNU Lesser General Public
00022 // License along with this library; if not, write to the Free Software
00023 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00024 // USA
00025 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00026 //
00027 // ************************************************************************
00028 //@HEADER
00029 
00030 #include <Epetra_matrix_data.h>
00031 #include <Epetra_Map.h>
00032 #include <Epetra_CrsMatrix.h>
00033 #include <Epetra_Util.h>
00034 
00035 namespace epetra_test {
00036 
00037 matrix_data::matrix_data(int num_rows,
00038                          int* rowlengths,
00039                          int blocksize)
00040  : numrows_(num_rows),
00041    numcols_(0),
00042    rows_(0),
00043    rowlengths_(0),
00044    blocksize_(blocksize),
00045    colindices_(0),
00046    coefs_(0)
00047 {
00048   if (numrows_ > 0) {
00049     rows_ = new int[numrows_];
00050     rowlengths_ = new int[numrows_];
00051     colindices_ = new int*[numrows_];
00052     coefs_ = new double*[numrows_];
00053     int dim = blocksize*blocksize;
00054     for(int i=0; i<numrows_; ++i) {
00055       rows_[i] = i;
00056       rowlengths_[i] = rowlengths[i];
00057       colindices_[i] = new int[rowlengths_[i]];
00058       coefs_[i] = new double[rowlengths_[i]*dim];
00059 
00060       for(int j=0; j<rowlengths_[i]; ++j) {
00061         colindices_[i][j] = 0;
00062         for(int k=0; k<dim; ++k) coefs_[i][j*dim+k] = 0.0;
00063       }
00064     }
00065   }
00066 }
00067 
00068 matrix_data::matrix_data(int num_rows,
00069                          int num_cols,
00070                          int num_off_diagonals,
00071                          int blocksize)
00072  : numrows_(num_rows),
00073    numcols_(num_cols),
00074    rows_(0),
00075    rowlengths_(0),
00076    blocksize_(blocksize),
00077    colindices_(0),
00078    coefs_(0)
00079 {
00080   if (numrows_ > 0) {
00081     rows_       = new int[numrows_];
00082     rowlengths_ = new int[numrows_];
00083     colindices_ = new int*[numrows_];
00084     coefs_      = new double*[numrows_];
00085 
00086     int max_row_length = 1+num_off_diagonals*2;
00087 
00088     for(int i=0; i<numrows_; ++i) {
00089       rows_[i] = i;
00090       if (i >= num_off_diagonals && numrows_-i > num_off_diagonals) {
00091         rowlengths_[i] = max_row_length;
00092       }
00093       else {
00094         if (i<num_off_diagonals) {
00095           rowlengths_[i] = 1+max_row_length/2+i;
00096         }
00097         else {
00098           rowlengths_[i] = 1+max_row_length/2+numrows_-i-1;
00099         }
00100       }
00101       colindices_[i] = new int[rowlengths_[i]];
00102       int dim = blocksize*blocksize;
00103       coefs_[i] = new double[rowlengths_[i]*dim];
00104 
00105       int first_col = i - max_row_length/2;
00106       if (first_col < 0) first_col = 0;
00107 
00108       for(int j=0; j<rowlengths_[i]; ++j) {
00109         colindices_[i][j] = first_col+j;
00110         for(int k=0; k<dim; ++k) {
00111           coefs_[i][j*dim+k] = 1.0;
00112         }
00113       }
00114     }
00115   }
00116 }
00117 
00118 static const int nodes_per_elem = 4;
00119 
00120 void get_node_ids(int elem_id, int* node_ids)
00121 {
00122   int first_node = 2*elem_id;
00123   for(int i=0; i<nodes_per_elem; ++i) node_ids[i] = first_node+i;
00124 }
00125 
00126 matrix_data::matrix_data(int num_quad_elements,
00127                          int num_dof_per_node,
00128        bool make_numerically_nonsymmetric)
00129  : numrows_(0),
00130    numcols_(0),
00131    rows_(0),
00132    rowlengths_(0),
00133    blocksize_(num_dof_per_node),
00134    colindices_(0),
00135    coefs_(0)
00136 {
00137   //Set up matrix-data representing a simple finite-element
00138   //mesh containing 2-D quad elements
00139   //
00140   //   *-----*-----*-----*
00141   //  0|    2|    4|    6|
00142   //   | 0   | 1   | ne-1|
00143   //   |     |     |     |
00144   //   *-----*-----*-----*
00145   //  1     3     5     7
00146   //
00147   //In the above drawing, 'ne' means num-elems. node-numbers are to the
00148   //lower-left of each node (*).
00149 
00150   numrows_ = num_quad_elements*2+2;
00151 
00152   if (numrows_ > 0) {
00153     rows_       = new int[numrows_];
00154     rowlengths_ = new int[numrows_];
00155     colindices_ = new int*[numrows_];
00156     coefs_      = new double*[numrows_];
00157 
00158     int i, j, k;
00159     for(i=0; i<numrows_; ++i) {
00160       rows_[i] = i;
00161       rowlengths_[i] = 0;
00162     }
00163 
00164     int* nodes = new int[nodes_per_elem];
00165     for(i=0; i<num_quad_elements; ++i) {
00166       get_node_ids(i, nodes);
00167 
00168       for(j=0; j<nodes_per_elem; ++j) {
00169         int node_j = nodes[j];
00170         for(k=0; k<nodes_per_elem; ++k) {
00171           int insertPoint = -1;
00172           int alloclen = rowlengths_[node_j];
00173           int offset = Epetra_Util_binary_search(nodes[k], colindices_[node_j],
00174                                                  rowlengths_[node_j], insertPoint);
00175           if (offset < 0) {
00176             Epetra_Util_insert(nodes[k], insertPoint,
00177                                colindices_[node_j], rowlengths_[node_j],
00178                                alloclen);
00179           }
00180         }
00181       }
00182     }
00183 
00184     int dim = blocksize_*blocksize_;
00185     for(i=0; i<numrows_; ++i) {
00186       int len = rowlengths_[i]*dim;
00187       coefs_[i] = new double[len];
00188       for(j=0; j<len; ++j) {
00189   if (make_numerically_nonsymmetric) {
00190     coefs_[i][j] = 1.0*(j+1);
00191   }
00192   else {
00193     coefs_[i][j] = 1.0;
00194   }
00195       }
00196     }
00197   }
00198 }
00199 
00200 matrix_data::~matrix_data()
00201 {
00202   for(int i=0; i<numrows_; ++i) {
00203     delete [] colindices_[i];
00204     delete [] coefs_[i];
00205   }
00206 
00207   delete [] colindices_; colindices_ = 0;
00208   delete [] coefs_; coefs_ = 0;
00209   delete [] rowlengths_; rowlengths_ = 0;
00210   delete [] rows_; rows_ = 0;
00211   numrows_ = 0;
00212 }
00213 
00214 double* matrix_data::coefs(int row, int col)
00215 {
00216   int insertPoint = -1;
00217   int row_idx = Epetra_Util_binary_search(row, rows_, numrows_,
00218                                           insertPoint);
00219   if (row_idx < 0) {
00220     cerr << "ERROR, row " << row
00221          << " not found in matrix_data"<<endl;
00222     return 0;
00223   }
00224 
00225   int col_idx = Epetra_Util_binary_search(col, colindices_[row_idx],
00226                                           rowlengths_[row_idx], insertPoint);
00227   if (col_idx < 0) {
00228     cerr << "ERROR, col " << col
00229          << " not found in matrix_data"<<endl;
00230     return 0;
00231   }
00232 
00233   int dim = blocksize_*blocksize_;
00234   return( &(coefs_[row_idx][col_idx*dim]) );
00235 }
00236 
00237 void matrix_data::copy_local_data_to_matrix(Epetra_CrsMatrix& A)
00238 {
00239   const Epetra_Map& rowmap = A.RowMap();
00240 
00241   for(int i=0; i<numrows_; ++i) {
00242     int row = rows_[i];
00243     if (rowmap.MyGID(row)) {
00244       int err = A.ReplaceGlobalValues(row, rowlengths_[i],
00245               coefs_[i], colindices_[i]);
00246       if (err < 0) {
00247   err = A.InsertGlobalValues(row, rowlengths_[i],
00248            coefs_[i], colindices_[i]);
00249       }
00250     }
00251   }
00252 }
00253 
00254 bool matrix_data::compare_local_data(const Epetra_CrsMatrix& A)
00255 {
00256   const Epetra_Map& map = A.RowMap();
00257   int numMyRows = map.NumMyElements();
00258   int* myRows = map.MyGlobalElements();
00259 
00260   Epetra_Util util;
00261 
00262   for(int i=0; i<numMyRows; ++i) {
00263     int row = myRows[i];
00264     int rowLen = A.NumGlobalEntries(row);
00265     if (rowLen != rowlengths_[row]) {
00266       return(false);
00267     }
00268 
00269     int* indices = new int[rowLen];
00270     double* values = new double[rowLen];
00271     A.ExtractGlobalRowCopy(row, rowLen, rowLen, values, indices);
00272 
00273     util.Sort(true, rowLen, indices, 1, &values, 0, 0);
00274 
00275     bool same = true;
00276     int* this_indices = colindices_[row];
00277     double* this_coefs = coefs_[row];
00278 
00279     for(int j=0; j<rowLen; ++j) {
00280       if (indices[j] != this_indices[j]) {
00281         same = false; break;
00282       }
00283       if (values[j] != this_coefs[j]) {
00284         same = false; break;
00285       }
00286     }
00287 
00288     delete [] indices;
00289     delete [] values;
00290 
00291     if (!same) return(false);
00292   }
00293 
00294   return(true);
00295 }
00296 
00297 }//namespace epetra_test
00298 

Generated on Thu Sep 18 12:37:58 2008 for Epetra Package Browser (Single Doxygen Collection) by doxygen 1.3.9.1