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