fei_AztecDMSR_Matrix.cpp

00001 /*--------------------------------------------------------------------*/
00002 /*    Copyright 2005 Sandia Corporation.                              */
00003 /*    Under the terms of Contract DE-AC04-94AL85000, there is a       */
00004 /*    non-exclusive license for use of this work by or on behalf      */
00005 /*    of the U.S. Government.  Export of this program may require     */
00006 /*    a license from the United States Government.                    */
00007 /*--------------------------------------------------------------------*/
00008 
00009 #include <fei_trilinos_macros.hpp>
00010 #include <fei_iostream.hpp>
00011 
00012 #ifdef HAVE_FEI_AZTECOO
00013 
00014 #include <assert.h>
00015 #include <stdlib.h>
00016 #include <math.h>
00017 
00018 #include <fei_mpi.h>
00019 
00020 #ifndef FEI_SER
00021 
00022 #define AZTEC_MPI AZTEC_MPI
00023 #define AZ_MPI AZ_MPI
00024 #ifndef MPI
00025 #define MPI MPI
00026 #endif
00027 
00028 #endif
00029 
00030 #include <az_aztec.h>
00031 
00032 #include <fei_ArrayUtils.hpp>
00033 
00034 #include <fei_Aztec_Map.hpp>
00035 #include <fei_Aztec_Vector.hpp>
00036 #include <fei_AztecDMSR_Matrix.hpp>
00037 
00038 #define ADMSR_LOCAL_ROW_ALLOC_LEN(localRow) 1+bindx[localRow+1]-bindx[localRow]
00039 
00040 #define ADMSR_LOCAL_ROW_LEN(localRow) 1+rowLengths_[localRow]
00041 
00042 #define ADMSR_GET_GLOBAL_COL(col, index) \
00043            col = index<N_update_ ? update_[orderingUpdate_[index]] : \
00044                                    external_[index-N_update_];
00045 
00046 namespace fei_trilinos {
00047 
00048 //==============================================================================
00049 AztecDMSR_Matrix::AztecDMSR_Matrix(Aztec_Map& map, int* update,
00050            bool linearDistribution)
00051   : isFilled_(false),
00052     isAllocated_(false),
00053     linearDistribution_(linearDistribution),
00054     localOffset_(map.localOffset()),
00055     localSize_(map.localSize()),
00056     amap_(map),
00057     Amat_(NULL),
00058     arraysAllocated_(false),
00059     val(NULL),
00060     bindx(NULL),
00061     rowLengths_(NULL),
00062     nnzeros_(0),
00063     N_update_(map.localSize()),
00064     update_(NULL),
00065     update_index_(NULL),
00066     orderingUpdate_(NULL),
00067     external_(NULL),
00068     extern_index_(NULL),
00069     data_org_(NULL),
00070     tmp_array_(0),
00071     tmp_array_len_(0),
00072     dtmp_array_(0),
00073     dtmp_array_len_(0),
00074     azTransformed_(false)
00075 {
00076 //   if (N_update_ <= 0)
00077 //     messageAbort("N_update_ (map.localSize()) <= 0.");
00078 
00079    if (N_update_ > 0) {
00080       update_ = new int[N_update_];
00081       rowLengths_ = new int[N_update_];
00082       for(int i=0; i<N_update_; i++) {
00083          update_[i] = update[i];
00084    rowLengths_[i] = 0;
00085       }
00086    }
00087 
00088    Amat_ = AZ_matrix_create(N_update_);
00089 
00090    update_index_ = NULL;
00091    orderingUpdate_ = NULL;
00092    external_ = NULL;
00093    extern_index_ = NULL;
00094    data_org_ = NULL;
00095 }
00096 
00097 //==============================================================================
00098 AztecDMSR_Matrix::AztecDMSR_Matrix(const AztecDMSR_Matrix& src)
00099   : isFilled_(src.isFilled_),
00100     isAllocated_(src.isAllocated_),
00101     linearDistribution_(src.linearDistribution_),
00102     localOffset_(src.localOffset_),
00103     localSize_(src.localSize_),
00104     amap_(src.amap_),
00105     Amat_(NULL),
00106     arraysAllocated_(src.arraysAllocated_),
00107     val(NULL),
00108     bindx(NULL),
00109     rowLengths_(NULL),
00110     nnzeros_(src.nnzeros_),
00111     N_update_(src.N_update_),
00112     update_(NULL),
00113     update_index_(NULL),
00114     orderingUpdate_(NULL),
00115     external_(NULL),
00116     extern_index_(NULL),
00117     data_org_(NULL),
00118     tmp_array_(0),
00119     tmp_array_len_(0),
00120     dtmp_array_(0),
00121     dtmp_array_len_(0),
00122     azTransformed_(src.azTransformed_)
00123 {
00124   expand_array(tmp_array_, tmp_array_len_, src.tmp_array_len_);
00125   expand_array(dtmp_array_, dtmp_array_len_, src.dtmp_array_len_);
00126 
00127   if (N_update_ > 0) {
00128     update_ = new int[N_update_];
00129     rowLengths_ = new int[N_update_];
00130     for(int i=0; i<N_update_; i++) {
00131       update_[i] = src.update_[i];
00132       rowLengths_[i] = src.rowLengths_[i];
00133     }
00134   }
00135 
00136   Amat_ = AZ_matrix_create(N_update_);
00137 
00138   if (isAllocated_ && nnzeros_ > 0) {
00139     val = new double[nnzeros_+1];
00140     bindx = new int[nnzeros_+1];
00141 
00142     for(int i=0; i<nnzeros_+1; ++i) {
00143       val[i] = src.val[i];
00144       bindx[i] = src.bindx[i];
00145     }
00146 
00147     if (isFilled_) {
00148       int j;
00149       int data_org_len = src.data_org_[AZ_total_send]+AZ_send_list;
00150       data_org_ = (int*)AZ_allocate(data_org_len*sizeof(int));
00151       for(j=0; j<data_org_len; ++j) data_org_[j] = src.data_org_[j];
00152 
00153       update_index_ = (int*)AZ_allocate((N_update_+1)*sizeof(int));
00154       orderingUpdate_ = new int[N_update_];
00155 
00156       for(j=0; j<N_update_+1; ++j) {
00157   update_index_[j] = src.update_index_[j];
00158       }
00159       for(j=0; j<N_update_; ++j) {
00160   orderingUpdate_[j] = src.orderingUpdate_[j];
00161       }
00162 
00163       int ext_len = data_org_[AZ_N_external]+1;
00164       external_ = (int*)AZ_allocate(ext_len*sizeof(int));
00165       extern_index_ = (int*)AZ_allocate(ext_len*sizeof(int));
00166 
00167       for(j=0; j<ext_len; ++j) {
00168   external_[j] = src.external_[j];
00169   extern_index_[j] = src.extern_index_[j];
00170       }
00171 
00172       AZ_set_MSR(Amat_, bindx, val, data_org_, 0, NULL, AZ_LOCAL);
00173     }
00174   }
00175 }
00176 
00177 //==============================================================================
00178 AztecDMSR_Matrix::~AztecDMSR_Matrix()
00179 {
00180   //
00181   //Destructor.
00182   //
00183   //The four arrays
00184   //   update_index_
00185   //   external_
00186   //   extern_index_
00187   //   data_org_
00188   //are allocated by an Aztec C function, and so will be
00189   //freed using 'free' instead of delete[].
00190   //
00191 
00192   if (arraysAllocated_) {
00193     delete [] val;
00194     val = NULL;
00195     delete [] bindx;
00196     bindx = NULL;
00197     arraysAllocated_ = false;
00198   }
00199 
00200   if (N_update_ > 0) {
00201     delete [] update_;
00202     update_ = NULL;
00203     delete [] rowLengths_;
00204     rowLengths_ = NULL;
00205     N_update_ = 0;
00206   }
00207 
00208   if (azTransformed_) {
00209     free(update_index_);
00210     update_index_ = NULL;
00211     delete [] orderingUpdate_;
00212     orderingUpdate_ = NULL;
00213     free(external_);
00214     external_ = NULL;
00215     free(extern_index_);
00216     extern_index_ = NULL;
00217     free(data_org_);
00218     data_org_ = NULL;
00219 
00220     azTransformed_ = false;
00221   }
00222 
00223   AZ_matrix_destroy(&Amat_);
00224   Amat_ = NULL;
00225 
00226   isFilled_ = false;
00227   isAllocated_ = false;
00228   arraysAllocated_ = false;
00229 
00230   if (tmp_array_len_ > 0) {
00231     delete [] tmp_array_;
00232     tmp_array_ = 0;
00233     tmp_array_len_ = 0;
00234   }
00235 }
00236 
00237 //==============================================================================
00238 void AztecDMSR_Matrix::expand_array(int*& array, int& arraylen, int newlen)
00239 {
00240   if (arraylen < newlen) {
00241     delete [] array;
00242     array = new int[newlen];
00243     for(int i=0; i<newlen; ++i) array[i] = -999;
00244     arraylen = newlen;
00245   }
00246 }
00247 
00248 //==============================================================================
00249 void AztecDMSR_Matrix::expand_array(double*& array, int& arraylen, int newlen)
00250 {
00251   if (arraylen < newlen) {
00252     delete [] array;
00253     array = new double[newlen];
00254     for(int i=0; i<newlen; ++i) array[i] = -999.9;
00255     arraylen = newlen;
00256   }
00257 }
00258 
00259 //==============================================================================
00260 void AztecDMSR_Matrix::scale(double scalar)
00261 {
00262   if (scalar == 1.0) return;
00263 
00264   if (val != NULL) {
00265     for(int i=0; i<nnzeros_+1; ++i) {
00266       val[i] *= scalar;
00267     }
00268   }
00269 }
00270 
00271 //==============================================================================
00272 void AztecDMSR_Matrix::matvec(const Aztec_Vector& x,
00273                               Aztec_Vector& y) const
00274 {
00275   //
00276   // This function forms the product y = Ax
00277   //
00278 
00279   assert(isFilled());
00280 
00281   int *proc_config = amap_.getProcConfig();
00282   //    int *idummy = 0, one=1;
00283   double *b = (double*)x.startPointer();
00284   double *c = (double*)y.startPointer();
00285 
00286   AZ_MSR_matvec_mult(b, c, Amat_, proc_config);
00287   //val,idummy,bindx,idummy,idummy,idummy,b,c,one, data_org_);
00288  
00289   return;
00290 }
00291 
00292 //==============================================================================
00293 void AztecDMSR_Matrix::put(double s)
00294 {
00295   if(isAllocated()){
00296     for(int i=0; i<bindx[N_update_]; i++) val[i] = s;
00297   }
00298   else {
00299     FEI_CERR << "AztecDMSR_Matrix::put - ERROR, can't do put until allocated"
00300    << FEI_ENDL;
00301   }
00302   return;
00303 }
00304 
00305 //==============================================================================
00306 int AztecDMSR_Matrix::rowLength(int row) const
00307 {
00308   int localRow;
00309 
00310   int thisRow = row;
00311 
00312   if(inUpdate(thisRow,localRow)){
00313     return(ADMSR_LOCAL_ROW_ALLOC_LEN(localRow));
00314   }
00315   else {
00316     FEI_CERR << "AztecDMSR_Matrix::rowLength: ERROR row " << row 
00317    << " not in local update set." << FEI_ENDL;
00318     abort();
00319     return(-1);
00320   }
00321 }
00322 
00324 int AztecDMSR_Matrix::setDiagEntry(int row, double value)
00325 {
00326   int thisRow = row;
00327   int localRow = -1;
00328   if(!inUpdate(thisRow,localRow)){
00329     FEI_CERR << "AztecDMSR_Matrix::setDiagEntry: ERROR - row " << row 
00330    << " not in local update set." << FEI_ENDL;
00331     abort(); return(-1);
00332   }
00333 
00334   val[localRow] = value;
00335   return(0);
00336 }
00337 
00339 double AztecDMSR_Matrix::getDiagEntry(int row) const
00340 {
00341   int thisRow = row;
00342   int localRow = -1;
00343   if(!inUpdate(thisRow,localRow)){
00344     FEI_CERR << "AztecDMSR_Matrix::getDiagEntry: ERROR - row " << row 
00345    << " not in local update set." << FEI_ENDL;
00346     abort(); return val[0];
00347   }
00348 
00349   return(val[localRow]);
00350 }
00351 
00353 int AztecDMSR_Matrix::getOffDiagRowPointers(int row, int*& colIndices,
00354               double*& coefs,
00355               int& offDiagRowLength)
00356 {
00357   int thisRow = row;
00358   int localRow = -1;
00359   if(!inUpdate(thisRow,localRow)){
00360     FEI_CERR << "AztecDMSR_Matrix::getOffDiagRowPointers: ERROR - row " << row 
00361    << " not in local update set." << FEI_ENDL;
00362     abort(); return(-1);
00363   }
00364 
00365   offDiagRowLength = rowLengths_[localRow];
00366   if (isFilled_) offDiagRowLength = ADMSR_LOCAL_ROW_ALLOC_LEN(localRow)-1;
00367   int start = bindx[localRow];
00368   colIndices = &(bindx[start]);
00369   coefs = &(val[start]);
00370 
00371   return(0);
00372 }
00373 
00375 void AztecDMSR_Matrix::getRow(int row,
00376                               int &length,
00377                               double *coefs,
00378                               int *colInd) const
00379 {
00380   //
00381   //Note to myself:
00382   //getRow, putRow and sumIntoRow are incredibly ugly, and need re-written.
00383   //
00384 
00385   bool foundDiagonal = false;
00386   double dtmp;
00387   int j, localRow, itmp;
00388 
00389   int thisRow = row;
00390 
00391   if(!inUpdate(thisRow,localRow)){
00392     FEI_CERR << "AztecDMSR_Matrix::getRow: ERROR - row " << row 
00393    << " not in local update set." << FEI_ENDL;
00394     length = 0;
00395     return;
00396   }
00397 
00398   int start = bindx[localRow];
00399   int end = bindx[localRow+1]-1;
00400 
00401   j = 0;
00402   for(int i=start; i<=end; i++){
00403     coefs[j] = val[i];
00404     if (isFilled()) {
00405       ADMSR_GET_GLOBAL_COL(colInd[j], bindx[i])
00406     }
00407     else {
00408       colInd[j] = bindx[i];
00409     }
00410 
00411     if (colInd[j]==row) {
00412       //we're at the diagonal element, so put it in.
00413       dtmp = coefs[j];
00414       itmp = colInd[j];
00415       coefs[j] = val[localRow];
00416       colInd[j] = row;
00417       j++;
00418       coefs[j] = dtmp;
00419       colInd[j] = itmp;
00420       foundDiagonal = true;
00421     }
00422     j++;
00423   }
00424 
00425   if(!foundDiagonal){ // still need to put in the diagonal
00426     coefs[j] = val[localRow];
00427     colInd[j] = row;
00428   }
00429 
00430   length = j+1;
00431   return;
00432 }
00433 
00435 int AztecDMSR_Matrix::putRow(int row, int len, const double *coefs,
00436                              const int *colInd)
00437 {
00438   //
00439   //Note to myself:
00440   //getRow, putRow and sumIntoRow are incredibly ugly, and need re-written.
00441   //
00442 
00443   int j, localRow, globalColIndex;
00444 
00445   int thisRow = row;
00446 
00447   if (!inUpdate(thisRow,localRow)){
00448     FEI_CERR << "AztecDMSR_Matrix::putRow: ERROR row " << row
00449    << " not in local update set." << FEI_ENDL;
00450     return(-1);
00451   }
00452 
00453   int offDiagRowLen = rowLengths_[localRow];
00454   int rowAllocLen = ADMSR_LOCAL_ROW_ALLOC_LEN(localRow);
00455   if (isFilled_) offDiagRowLen = rowAllocLen - 1;
00456   int offDiagRowAllocLen = rowAllocLen - 1;
00457 
00458   if (len > rowAllocLen) {
00459     FEI_CERR << "AztecDMSR_Matrix::putRow. too many coefs, row " << row << FEI_ENDL;
00460     return(-1);
00461   }
00462 
00463   int jLimit = bindx[localRow+1] - 1;
00464   int jStart = bindx[localRow];
00465   int* colInds = &(bindx[jStart]);
00466   int colStart = *colInds;
00467   double* rowCoefs = &(val[jStart]);
00468 
00469   for(int i=0; i<len; i++){
00470     if (colInd[i] == row){  //it's on the diagonal
00471       val[localRow] = coefs[i];
00472       continue;
00473     }
00474 
00475     //look along the row until we find the right col index
00476     //or an empty spot
00477     j=jStart;
00478     if (isFilled()){
00479       int colIndex = colStart;
00480       int col = colInd[i];
00481       while (j <= jLimit) {
00482   ADMSR_GET_GLOBAL_COL(globalColIndex, colIndex)
00483 
00484   if (globalColIndex == col) break;
00485 
00486   colIndex = bindx[++j];
00487       }
00488 
00489       //now put the coefficient in if we haven't gone too far
00490       //along the row
00491       if (j <= jLimit) {
00492   val[j] = coefs[i];
00493       }
00494       else {
00495   FEI_CERR << "AztecDMSR_Matrix::putRow: ERROR didn't "
00496        << " find col. index " << colInd[i] << " in "
00497        << "row " << row << FEI_ENDL;
00498   return(-1);
00499       }
00500     }
00501     else { // !isFilled()
00502 
00503       //first, look for colInd[i] in the row
00504       int col = colInd[i];
00505       int insertPoint = -1;
00506       int index = fei::binarySearch<int>(col, colInds, offDiagRowLen, insertPoint);
00507 
00508       if (index >= offDiagRowAllocLen){ //bad index
00509   FEI_CERR << "AztecDMSR_Matrix::putRow, ERROR: " 
00510        << "row " << row << ", colInd["<<i<<"] " << colInd[i]
00511        << ", index = " << index << FEI_ENDL;
00512   return(-1);
00513       }
00514 
00515       if (index >= 0) {
00516   rowCoefs[index] = coefs[i];
00517       }
00518       else {
00519   int tmp = offDiagRowLen;
00520   int err = insert(col, insertPoint, colInds,
00521        tmp, offDiagRowAllocLen);
00522   err += insert(coefs[i], insertPoint, rowCoefs,
00523           offDiagRowLen, offDiagRowAllocLen);
00524   if (err != 0) {
00525     FEI_CERR << "AztecDMSR_Matrix::putRow ERROR: failed to add "
00526          << "value for index " << col << " to row " << row << FEI_ENDL;
00527     return(-1);
00528   }
00529   rowLengths_[localRow]++;
00530       }
00531     }
00532   }
00533 
00534   return(0);
00535 }
00536 
00538 int AztecDMSR_Matrix::sumIntoRow(int numRows, const int* rows,
00539                                  int numCols, const int *colInd,
00540          const double* const* coefs)
00541 {
00542   if (numRows == 0 || numCols == 0) return(0);
00543 
00544   if (!isFilled_) {
00545     int err = 0;
00546     for(int i=0; i<numRows; ++i) {
00547       err = sumIntoRow(rows[i], numCols, coefs[i], colInd);
00548       if (err != 0) return(err);
00549     }
00550 
00551     return(0);
00552   }
00553 
00554   //Now for the harder (but more important) case where isFilled_ == true.
00555 
00556   //first compute max-row-length:
00557   int maxRowLen = 0;
00558   for(int i=0; i<numRows; ++i) {
00559     int row = rows[i];
00560     int localRow;
00561     if (!inUpdate(row, localRow)) {
00562       FEI_CERR << "AztecDMSR_Matrix::sumIntoRow: ERROR row " << row
00563          << " not in local update set [" << update_[0] << " ... "
00564          << update_[N_update_-1] << "]." << FEI_ENDL;
00565       return(-1);
00566     }
00567 
00568     int rowlen = bindx[localRow+1]-bindx[localRow];
00569     if (maxRowLen < rowlen) maxRowLen = rowlen;
00570   }
00571 
00572   if (maxRowLen+2*numCols > tmp_array_len_) {
00573     expand_array(tmp_array_, tmp_array_len_, maxRowLen+2*numCols);
00574   }
00575 
00576   int* incols = &tmp_array_[maxRowLen];
00577   int* indirect = incols+numCols;
00578 
00579   for(int jj=0; jj<numCols; ++jj) {
00580     incols[jj] = colInd[jj];
00581     indirect[jj] = jj;
00582   }
00583 
00584   fei::insertion_sort_with_companions<int>(numCols, incols, indirect);
00585 
00586   int row, localRow;
00587 
00588   for(int i=0; i<numRows; ++i) {
00589     row = rows[i];
00590     if (!inUpdate(row, localRow)) {
00591       FEI_CERR << "AztecDMSR_Matrix::sumIntoRow: ERROR row " << row
00592          << " not in local update set [" << update_[0] << " ... "
00593          << update_[N_update_-1] << "]." << FEI_ENDL;
00594       return(-1);
00595     }
00596 
00597     int jStart = bindx[localRow];
00598     double* rowCoefs = val+jStart;
00599     int* rowColInds = bindx+jStart;
00600     int rowLen= bindx[localRow+1]-jStart;
00601 
00602     for(int jj=0; jj<rowLen; ++jj) {
00603       ADMSR_GET_GLOBAL_COL(tmp_array_[jj], rowColInds[jj]);
00604     }
00605 
00606     const double* coefs_i = coefs[i];
00607 
00608     int inoffset = 0;
00609     int incol = incols[inoffset];
00610     while (incol == row) {
00611       val[localRow] += coefs_i[indirect[inoffset++]];
00612       if (inoffset >= numCols) break;
00613       incol = incols[inoffset];
00614     }
00615 
00616     if (inoffset >= numCols) continue;
00617 
00618     //rowOffset is the offset into the row at which incol appears.
00619 
00620     int rowOffset = fei::binarySearch<int>(incol, tmp_array_, rowLen);
00621     if (rowOffset < 0) {
00622        FEI_CERR << "AztecDMSR_Matrix::sumIntoRow, ERROR: "
00623              << "row " << row << ", col not found: "
00624              << incol << FEI_ENDL;
00625       return(-1);
00626     }
00627 
00628     rowCoefs[rowOffset++] += coefs_i[indirect[inoffset++]];
00629 
00630     //check whether incols has a repeated column-index
00631     if (inoffset>0 && incols[inoffset] == incols[inoffset-1]) --rowOffset;
00632 
00633     while(inoffset < numCols) {
00634       incol = incols[inoffset];
00635 
00636       if (incol == row) {
00637         val[localRow] += coefs_i[indirect[inoffset++]];
00638         continue;
00639       }
00640 
00641       while(tmp_array_[rowOffset] != incol) {
00642         ++rowOffset;
00643         if (rowOffset >= rowLen) {
00644           FEI_CERR << "AztecDMSR_Matrix::sumIntoRow, ERROR, col "
00645              << incol << " not found in row " << row << FEI_ENDL;
00646           return(-1); 
00647         }
00648       }
00649 
00650       rowCoefs[rowOffset++] += coefs_i[indirect[inoffset++]];
00651       if (inoffset>0 && incols[inoffset] == incols[inoffset-1]) --rowOffset;
00652     }
00653   }
00654 
00655   return(0);
00656 }
00657 
00659 int AztecDMSR_Matrix::sumIntoRow(int row, int len, const double *coefs,
00660                                  const int *colInd)
00661 {
00662   //
00663   //Note to myself:
00664   //getRow, putRow and sumIntoRow are incredibly ugly, and need re-written.
00665   //
00666 
00667   int localRow, thisRow = row ;
00668 
00669   if (!inUpdate(thisRow,localRow)) {
00670     FEI_CERR << "AztecDMSR_Matrix::sumIntoRow: ERROR row " << row
00671    << " not in local update set." << FEI_ENDL;
00672     return(-1);
00673   }
00674 
00675   int jLimit = bindx[localRow+1] - 1;
00676   int jStart = bindx[localRow];
00677   int jLen = jLimit-jStart+1;
00678   int* colInds = &(bindx[jStart]);
00679   double* rowCoefs = &(val[jStart]);
00680 
00681   if (isFilled_) {
00682     if (jLen+len > tmp_array_len_) {
00683       expand_array(tmp_array_, tmp_array_len_, jLen+len);
00684       expand_array(dtmp_array_, dtmp_array_len_, len);
00685     }
00686     for(int jj=0; jj<jLen; ++jj) {
00687       ADMSR_GET_GLOBAL_COL(tmp_array_[jj], colInds[jj])
00688     }
00689 
00690     int* incols = &tmp_array_[jLen];
00691     int doffs = 0;
00692     for(int jj=0; jj<len; ++jj) {
00693       int col = colInd[jj];
00694       if (col == row) {
00695         val[localRow] += coefs[jj];
00696       }
00697       else {
00698         incols[doffs] = col;
00699         dtmp_array_[doffs++] = coefs[jj];
00700       }
00701     }
00702     fei::insertion_sort_with_companions<double>(doffs, incols, dtmp_array_);
00703 
00704     int ioffset = 0;
00705     int offset = fei::binarySearch<int>(incols[ioffset], tmp_array_, jLen);
00706     if (offset < 0) {
00707       FEI_CERR << "AztecDMSR_Matrix::sumIntoRow, ERROR: "
00708              << "row " << row << ", col not found: "
00709              << colInd[ioffset] << FEI_ENDL;
00710       return(-1);
00711     }
00712 
00713     rowCoefs[offset++] += dtmp_array_[ioffset++];
00714     if (incols[ioffset] == tmp_array_[offset-1]) --offset;
00715 
00716     while(ioffset < doffs) {
00717       int incol = incols[ioffset];
00718 
00719       while(tmp_array_[offset] != incol) {
00720         ++offset;
00721         if (offset >= jLen) {
00722           FEI_CERR << "AztecDMSR_Matrix::sumIntoRow, ERROR, col "
00723              << incols[ioffset] << " not found in row " << row << FEI_ENDL;
00724           return(-1);
00725         }
00726       }
00727       rowCoefs[offset++] += dtmp_array_[ioffset++];
00728       if (incols[ioffset] == tmp_array_[offset-1]) --offset;
00729     }
00730 
00731     return(0);
00732   }
00733 
00734   //if we get to here, then we know that isFilled_ is false...
00735 
00736   int rowAllocLen = ADMSR_LOCAL_ROW_ALLOC_LEN(localRow);
00737   int offDiagRowLen = isFilled_ ? rowAllocLen - 1 : rowLengths_[localRow];
00738   int offDiagRowAllocLen = rowAllocLen - 1;
00739 
00740   for(int i=0; i<len; i++){
00741     if (colInd[i] == row){  //it's on the diagonal
00742       val[localRow] += coefs[i];
00743       continue;
00744     }
00745 
00746     //find the right col index in the row, or an empty spot
00747     int col = colInd[i];
00748     int insertPoint = -1;
00749     int index = fei::binarySearch<int>(col, colInds, offDiagRowLen, insertPoint);
00750 
00751     if (index >= 0) {
00752       rowCoefs[index] += coefs[i];
00753     }
00754     else {
00755       int tmp = offDiagRowLen;
00756       int err = insert(col, insertPoint, colInds,
00757                        tmp, offDiagRowAllocLen);
00758       err += insert(coefs[i], insertPoint, rowCoefs,
00759                     offDiagRowLen, offDiagRowAllocLen);
00760       if (err != 0) {
00761         FEI_CERR << "AztecDMSR_Matrix::sumIntoRow ERROR: failed to add "
00762                  << "value for index " << col << " to row " << row << FEI_ENDL;
00763         return(-1);
00764       }
00765       rowLengths_[localRow]++;
00766     }
00767   }
00768 
00769   return(0);
00770 }
00771 
00772 //==============================================================================
00773 int AztecDMSR_Matrix::addScaledMatrix(double scalar,
00774               const AztecDMSR_Matrix& source)
00775 {
00776   if (N_update_ != source.N_update_ ||
00777       nnzeros_ != source.nnzeros_ ||
00778       isFilled_ != source.isFilled_) {
00779     FEI_CERR << "AztecDMSR_Matrix::addScaledMatrix ERROR, not compatible"
00780    << FEI_ENDL;
00781     return(-1);
00782   }
00783 
00784   const double* src_val = source.val;
00785   int i;
00786   for(i=0; i<N_update_; ++i) {
00787     val[i] += scalar*src_val[i];
00788   }
00789 
00790   //val[N_update_] is not used.
00791 
00792   if (scalar == 1.0) {
00793     for(i=N_update_+1; i<nnzeros_+1; ++i) {
00794       val[i] += src_val[i];
00795     }
00796   }
00797   else {
00798     for(i=N_update_+1; i<nnzeros_+1; ++i) {
00799       val[i] += scalar*src_val[i];
00800     }
00801   }
00802 
00803   return(0);
00804 }
00805 
00806 //==============================================================================
00807 int AztecDMSR_Matrix::insert(int item, int offset, int* list,
00808            int& len, int allocLen)
00809 {
00810   if (len >= allocLen) return(-1);
00811 
00812   for(int i=len; i>offset; i--) list[i] = list[i-1];
00813 
00814   list[offset] = item;
00815   len++;
00816 
00817   return(0);
00818 }
00819 
00820 //==============================================================================
00821 int AztecDMSR_Matrix::insert(double item, int offset, double* list,
00822            int& len, int allocLen)
00823 {
00824   if (len >= allocLen) return(-1);
00825 
00826   for(int i=len; i>offset; i--) list[i] = list[i-1];
00827 
00828   list[offset] = item;
00829   len++;
00830 
00831   return(0);
00832 }
00833 
00834 //==============================================================================
00835 void AztecDMSR_Matrix::getDiagonal(Aztec_Vector& diagVector) const {
00836   
00840 // have each processor form its piece of the diagonal
00841     double *pdv = (double*)diagVector.startPointer();
00842 
00843     for(int i=0; i<N_update_; i++){
00844         pdv[i] = val[i];
00845     }
00846 }
00847 
00849 void AztecDMSR_Matrix::allocate(int *rowLengths)
00850 {
00851   //
00852   //We assume that 'rowLengths' is of length N_update_.
00853   //
00854   //rowLengths contains the length of each local row, *NOT* including the
00855   //coefficient on the diagonal.
00856   //
00857   int i;
00858 
00859   //first, count how many non-zeros there are in the local submatrix
00860 
00861   nnzeros_ = 0;
00862   for(i=0; i<N_update_; i++){
00863     if (rowLengths[i] < 0) {
00864       messageAbort("allocate: negative row length");
00865     }
00866     nnzeros_ += rowLengths[i] + 1;
00867   }
00868 
00869   if (bindx != NULL) {
00870     delete [] bindx;
00871   }
00872   bindx = new int[nnzeros_+1];
00873 
00874   if (val != NULL) {
00875     delete [] val;
00876   }
00877   val = new double[nnzeros_+1];
00878 
00879   arraysAllocated_ = true;
00880 
00881   for(i=0; i<nnzeros_+1; i++){
00882     val[i] = 0.0;
00883     bindx[i] = -1;
00884   }
00885 
00886   bindx[0] = N_update_+1;
00887 
00888   for(i=0; i<N_update_; i++){
00889     bindx[i+1] = bindx[i] + rowLengths[i];
00890     if (bindx[i+1] < 0) {
00891       messageAbort("allocate: bindx row length negative.");
00892     }
00893   }
00894 
00895   //val[N_update_] not used by aztec but we'll initialize it anyway...
00896   val[N_update_] = 0.0;
00897 
00898   AZ_set_MSR(Amat_, bindx, val, data_org_, 0, NULL, AZ_LOCAL);
00899 
00900   setAllocated(true);
00901   return;
00902 }
00903 
00905 void AztecDMSR_Matrix::allocate(int *rowLengths,
00906         const int* const* colIndices)
00907 {
00908   allocate(rowLengths);
00909   int col;
00910   
00911   int offset = N_update_+1;
00912   for(int i=0; i<N_update_; ++i) {
00913     const int* row_cols = colIndices[i];
00914     int row_len = rowLengths[i];
00915     rowLengths_[i] = row_len-1;
00916 
00917     int prev_col = -999;
00918     int coffset = 0;
00919     for(int j=0; j<row_len; ++j) {
00920       col = row_cols[coffset++];
00921 
00922       if (col == localOffset_+i) {
00923   col = row_cols[coffset++];
00924       }
00925 
00926       if (col <= prev_col) {
00927   messageAbort("allocate: column-indices not sorted.");
00928       }
00929 
00930       prev_col = col;
00931 
00932       bindx[offset++] = col;
00933     }
00934   }
00935 }
00936 
00938 double AztecDMSR_Matrix::rowMax(int row) const {
00939     int localRow;
00940     double max = 0.0;
00941 
00942     if(!inUpdate(row,localRow)){
00943         FEI_CERR << "AztecDMSR_Matrix::rowMax: ERROR row " << row 
00944              << " not in local update set." << FEI_ENDL;
00945         return(-1.0);
00946     }
00947 
00948     max = fabs(val[localRow]);
00949 
00950     for(int i=bindx[localRow]; i<bindx[localRow+1]; i++)
00951         if(fabs(val[i])>max)max = fabs(val[i]);
00952 
00953     return(max);
00954 }
00955 
00957 void AztecDMSR_Matrix::fillComplete() {
00958 /*
00959    This is where we call the Aztec function AZ_transform, which calculates
00960    communication parameters and re-orders the equations for use as a
00961    global distributed matrix.
00962 */
00963   if (isFilled_ || azTransformed_) {
00964     isFilled_ = true;
00965     azTransformed_ = true;
00966     return;
00967   }
00968 
00969     int *proc_config = amap_.getProcConfig();
00970     int *dummy = 0;
00971 
00972    //before we turn Aztec loose on the matrix, lets do a quick check on the
00973    //indices to try to make sure none of them are garbage...
00974    int globalSize = amap_.globalSize();
00975    for(int i=N_update_+1; i<nnzeros_+1; i++) {
00976       if (bindx[i] < 0 || bindx[i] >= globalSize) {
00977          FEI_CERR << "AztecDMSR_Matrix: ERROR, bindx["<<i<<"]: " << bindx[i]
00978               << ", globalSize: " << globalSize << FEI_ENDL;
00979 #ifndef FEI_SER
00980          MPI_Comm thisComm = amap_.getCommunicator();
00981          MPI_Abort(thisComm, -1);
00982 #endif
00983       }
00984    }
00985 
00986     AZ_transform(proc_config, &external_, bindx, val, update_, &update_index_,
00987                  &extern_index_, &data_org_, N_update_,
00988                  dummy, dummy, dummy, &dummy, AZ_MSR_MATRIX);
00989 
00990 //AZ_transform allocates these arrays:
00991 //  external_
00992 //  update_index_
00993 //  extern_index_
00994 //  data_org_
00995 //
00996 //On return from AZ_transform, the array update_index contains a mapping
00997 //to the local re-ordering of the indices of the update array. Now we will fill
00998 //the orderingUpdate array with the reverse of that mapping. i.e., a record
00999 //of how to get back to the original ordering of the update indices.
01000 
01001     AZ_set_MSR(Amat_, bindx, val, data_org_, 0, NULL, AZ_LOCAL);
01002 
01003     orderingUpdate_ = new int[N_update_];
01004     for(int ii=0; ii<N_update_; ii++) {
01005       orderingUpdate_[update_index_[ii]] = ii;
01006     }
01007 
01008     azTransformed_ = true;
01009 
01010     setFilled(true);
01011     return;
01012 }
01013 
01014 //==============================================================================
01015 void AztecDMSR_Matrix::copyStructure(AztecDMSR_Matrix& source)
01016 {
01017   //
01018   //This function copies the structure (essentially just the bindx and
01019   //rowLengths_ arrays) and other relevant variables from the 'source' matrix.
01020   //The result is that 'this' matrix is laid out the same as 'source'.
01021   //
01022   //This is not the greatest C++, it would be better to use a copy-
01023   //constructor, but there are also reasons why a copy-constructor isn't
01024   //always a good thing for high-performance computing. (e.g., I don't want
01025   //people passing one of these matrices by value.) For now, this
01026   //function is a suitable band-aid solution.
01027   //
01028   nnzeros_ = source.nnzeros_;
01029 
01030   if (arraysAllocated_) {
01031     delete [] val;
01032     delete [] bindx;
01033     arraysAllocated_ = false;
01034   }
01035 
01036   val = new double[nnzeros_+1];
01037   bindx = new int[nnzeros_+1];
01038 
01039   int i;
01040   for(i=0; i<nnzeros_+1; i++) {
01041     val[i] = 0.0;
01042     bindx[i] = source.bindx[i];
01043   }
01044 
01045   for(i=0; i<N_update_; ++i) rowLengths_[i] = source.rowLengths_[i];
01046 
01047   AZ_set_MSR(Amat_, bindx, val, data_org_, 0, NULL, AZ_LOCAL);
01048 
01049   if (source.isFilled_) {
01050     if (!isFilled_) {
01051 
01052       update_index_ = (int*)AZ_allocate(N_update_*sizeof(int));
01053       orderingUpdate_ = new int[N_update_];
01054       for(i=0; i<N_update_; ++i) {
01055   update_index_[i] = source.update_index_[i];
01056   orderingUpdate_[i] = source.orderingUpdate_[i];
01057       }
01058 
01059       int dlen = AZ_COMM_SIZE + source.data_org_[AZ_total_send];
01060       data_org_ = (int*)AZ_allocate(dlen*sizeof(int));
01061       for(i=0; i<dlen; ++i) data_org_[i] = source.data_org_[i];
01062 
01063       external_ = (int*)AZ_allocate(data_org_[AZ_N_external]*sizeof(int));
01064       extern_index_ = (int*)AZ_allocate(data_org_[AZ_N_external]*sizeof(int));
01065       for(i=0; i<data_org_[AZ_N_external]; ++i) {
01066   external_[i] = source.external_[i];
01067   extern_index_[i] = source.extern_index_[i];
01068       }
01069       isFilled_ = source.isFilled_;
01070       azTransformed_ = true;
01071     }
01072   }
01073 
01074   arraysAllocated_ = true;
01075   setAllocated(true);
01076 }
01077 
01078 //=============================================================================
01079 bool AztecDMSR_Matrix::readFromFile(const char *filename)
01080 {
01081   /*
01082     This function reads the matrix data from a matrix-market format data file,
01083     which looks like:
01084 
01085     n
01086     i j val
01087     .
01088     .
01089 
01090     Important note: we are going to assume that the index-base of the indices in
01091     the file are 0.
01092   */
01093   int i, j, dummy;
01094   double value;
01095   char line[128];
01096 
01097   FILE *mfp = fopen(filename,"r");
01098 
01099   if(!mfp){
01100     FEI_CERR << "AztecDMSR_Matrix::readFromFile - couldn't open matrix file."
01101    << FEI_ENDL;
01102     return(false);
01103   }
01104 
01105   if (strstr(filename, ".mtx") == NULL) {
01106     FEI_CERR << "AztecDMSR_Matrix::readFromFile: filename doesn't contain "
01107    << "'.mtx'. File should be a MatrixMarket file." << FEI_ENDL;
01108     return(false);
01109   }
01110 
01111   do {
01112     fgets(line,128,mfp);
01113   } while(strchr(line,'%'));
01114   
01115   while(!feof(mfp)){
01116     do {
01117       fgets(line,128,mfp);
01118     } while(strchr(line,'%'));
01119     if(feof(mfp)){
01120       fclose(mfp);
01121       return(true);
01122     }
01123     sscanf(line,"%d %d %le",&i,&j,&value);
01124 
01125     if(inUpdate(i, dummy)) {
01126       if (putRow(i, 1, &value, &j) != 0) return(false);
01127     }
01128   }
01129 
01130   fclose(mfp);
01131   return(true);
01132 }
01133 
01135 bool AztecDMSR_Matrix::writeToFile(const char *fileName) const
01136 {
01137   /* Write the matrix into the file "fileName", using the format:
01138      n
01139      i j val
01140      .
01141      .
01142   */
01143 
01144   int numProcs = amap_.getProcConfig()[AZ_N_procs];
01145   int thisProc = amap_.getProcConfig()[AZ_node];
01146   int masterRank = 0;
01147 
01148 
01149   int localNNZ = nnzeros_;
01150   int globalNNZ = localNNZ;
01151 #ifndef FEI_SER
01152   MPI_Comm thisComm = amap_.getCommunicator();
01153   MPI_Allreduce(&localNNZ, &globalNNZ, 1, MPI_INT, MPI_SUM, thisComm);
01154 #endif
01155 
01156   for(int p=0; p<numProcs; p++){
01157 
01158     //A barrier inside the loop so each processor waits its turn.
01159 #ifndef FEI_SER
01160     MPI_Barrier(thisComm);
01161 #endif
01162  
01163     if (p == thisProc) {
01164       FILE *file = NULL;
01165 
01166       if (masterRank == thisProc) {
01167   //This is the master processor, open a new file.
01168   file = fopen(fileName,"w");
01169 
01170   //Write the matrix dimensions n and n (rows==cols) into the file,
01171   //along with the global number-of-nonzeros globalNNZ.
01172 
01173   int n = amap_.globalSize();
01174   fprintf(file,"%d %d %d\n",n, n, globalNNZ);
01175       }
01176       else {
01177   //This is a non-master node, open the file for appending to.
01178   file = fopen(fileName,"a");
01179       }
01180 
01181       //Now loop over the local portion of the matrix.
01182       for(int i=0; i<N_update_; i++){
01183   int row = localOffset_+i;
01184 
01185   int localRow = -1;
01186   if (!inUpdate(row, localRow)) return(false);
01187 
01188   int offDiagRowLen = ADMSR_LOCAL_ROW_LEN(localRow) - 1;
01189   if (isFilled_) offDiagRowLen = ADMSR_LOCAL_ROW_ALLOC_LEN(localRow) - 1;
01190   int* colInds = &(bindx[bindx[localRow]]);
01191   double* coefs = &(val[bindx[localRow]]);
01192 
01193   bool wroteDiagonal = false;
01194   for(int j=0; j<offDiagRowLen; j++) {
01195     int col = colInds[j];
01196     int globalCol = col;
01197     if (isFilled()) {
01198       ADMSR_GET_GLOBAL_COL(globalCol, col)
01199     }
01200 
01201     if (globalCol >= row && !wroteDiagonal) {
01202       fprintf(file,"%d %d %20.13e\n", row, row,
01203         val[localRow]);
01204       wroteDiagonal = true;
01205     }
01206 
01207     fprintf(file,"%d %d %20.13e\n", row, globalCol, coefs[j]);
01208   }
01209   if (!wroteDiagonal) {
01210     fprintf(file,"%d %d %20.13e\n", row, row,
01211       val[localRow]);
01212     wroteDiagonal = true;
01213   }
01214       }
01215       fclose(file);
01216     }
01217   }
01218 
01219   return(true);
01220 }
01221 
01222 //==============================================================================
01223 int AztecDMSR_Matrix::inUpdate(int globalIndex, int& localIndex) const
01224 {
01225   //
01226   // This function determines whether globalIndex is in the local update set,
01227   // and if it is, returns in localIndex the local index for it. If update_index_
01228   // has already been allocated and set (by AZ_transform) then localIndex is
01229   // taken from there.
01230   // If globalIndex is not in the update set, inUpdate returns 0.
01231   //
01232   if (linearDistribution_) {
01233     localIndex = globalIndex - localOffset_;
01234   }
01235   else {
01236     localIndex = AZ_find_index(globalIndex, update_, N_update_);
01237   }
01238 
01239   if(localIndex<0 || localIndex>=localSize_) {localIndex = -1; return(0); }
01240 
01241   if(azTransformed_){
01242     localIndex = update_index_[localIndex];
01243   }
01244 
01245   return(1);
01246 }
01247 
01248 //==============================================================================
01249 void AztecDMSR_Matrix::messageAbort(const char* mesg) {
01250     FEI_CERR << "AztecDMSR_Matrix: ERROR: " << mesg << " Aborting." << FEI_ENDL;
01251     abort();
01252 }
01253 
01254 }//namespace fei_trilinos
01255 
01256 #endif
01257 //HAVE_FEI_AZTECOO
01258 

Generated on Tue Jul 13 09:27:44 2010 for FEI by  doxygen 1.4.7