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