fei_AztecDVBR_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 <stdio.h>
00017 #include <cstring>
00018 
00019 #include <fei_mpi.h>
00020 
00021 #ifndef FEI_SER
00022 
00023 #define AZTEC_MPI AZTEC_MPI
00024 #define AZ_MPI AZ_MPI
00025 #ifndef MPI
00026 #define MPI MPI
00027 #endif
00028 
00029 #endif
00030 
00031 #include <az_aztec.h>
00032 #include <fei_Aztec_Map.hpp>
00033 #include <fei_Aztec_BlockMap.hpp>
00034 #include <fei_Aztec_LSVector.hpp>
00035 #include <fei_AztecDVBR_Matrix.hpp>
00036 
00037 namespace fei_trilinos {
00038 
00039 //==============================================================================
00040 AztecDVBR_Matrix::AztecDVBR_Matrix(fei::SharedPtr<Aztec_BlockMap> map)
00041   : amap_(map),
00042     Amat_(NULL),
00043     N_update_(map->getNumLocalBlocks()),
00044     external_(NULL),
00045     extern_index_(NULL),
00046     update_index_(NULL),
00047     data_org_(NULL),
00048     orderingUpdate_(NULL),
00049     isLoaded_(false),
00050     isAllocated_(false),
00051     localNNZ_(0),
00052     nnzPerRow_(NULL),
00053     numRemoteBlocks_(0),
00054     remoteInds_(NULL),
00055     remoteBlockSizes_(NULL)
00056 {
00057     nnzPerRow_ = new int[N_update_];
00058 
00059     for(int i=0; i<N_update_; i++) {
00060        nnzPerRow_[i] = 0;
00061     }
00062 
00063     Amat_ = AZ_matrix_create(N_update_);
00064     Amat_->matrix_type = AZ_VBR_MATRIX;
00065     Amat_->matvec =
00066             (void(*)(double*,double*,AZ_MATRIX_STRUCT*,int*))AZ_VBR_matvec_mult;
00067 
00068     //now we can allocate and fill the rpntr array.
00069     Amat_->rpntr = NULL;
00070     calcRpntr();
00071 }
00072 
00073 //==============================================================================
00074 AztecDVBR_Matrix::AztecDVBR_Matrix(const AztecDVBR_Matrix& src)
00075  : amap_(src.amap_),
00076    Amat_(NULL),
00077    N_update_(src.N_update_),
00078    external_(NULL),
00079    extern_index_(NULL),
00080    update_index_(NULL),
00081    data_org_(NULL),
00082    orderingUpdate_(NULL),
00083    isLoaded_(src.isLoaded_),
00084    isAllocated_(src.isAllocated_),
00085    localNNZ_(src.localNNZ_),
00086    nnzPerRow_(NULL),
00087    numRemoteBlocks_(0),
00088    remoteInds_(NULL),
00089    remoteBlockSizes_(NULL)
00090 {
00091 //
00092 //This copy constructor just takes a reference to src's amap_ (above), but
00093 //'deep'-copies everything else. i.e., all arrays are allocated here and copied
00094 //from src, etc.
00095 //
00096 //When this constructor completes, this matrix should be in the same state as
00097 //the 'src' matrix. i.e., if src is already allocated, then this matrix will
00098 //have its structure allocated. If src is already loaded (AZ_transform'd) then
00099 //this matrix will have all arrays resulting from AZ_transform allocated too.
00100 //The only thing this matrix won't get is the coefficient data from src.
00101 //
00102    nnzPerRow_ = new int[N_update_];
00103 
00104    int i;
00105    for(i=0; i<N_update_; i++) {
00106       nnzPerRow_[i] = src.nnzPerRow_[i];
00107    }
00108 
00109    Amat_ = AZ_matrix_create(N_update_);
00110    Amat_->matrix_type = AZ_VBR_MATRIX;
00111    Amat_->matvec =
00112             (void(*)(double*,double*,AZ_MATRIX_STRUCT*,int*))AZ_VBR_matvec_mult;
00113 
00114    Amat_->rpntr = NULL;
00115    calcRpntr();
00116 
00117    if (isAllocated_) {
00118       Amat_->bpntr = new int[N_update_+1];
00119       for(i=0; i<N_update_+1; i++) Amat_->bpntr[i] = src.Amat_->bpntr[i];
00120 
00121       int totalNNZBlks = Amat_->bpntr[N_update_];
00122       Amat_->bindx = new int[totalNNZBlks];
00123       for(i=0; i<totalNNZBlks; i++) Amat_->bindx[i] = src.Amat_->bindx[i];
00124 
00125       Amat_->indx = new int[totalNNZBlks+1];
00126       for(i=0; i<totalNNZBlks+1; i++) Amat_->indx[i] = src.Amat_->indx[i];
00127 
00128       Amat_->val = new double[localNNZ_];
00129       for(i=0; i<localNNZ_; i++) Amat_->val[i] = 0.0;
00130    }
00131 
00132    if (isLoaded_) {
00133       int dataOrgLength = src.data_org_[AZ_total_send] + AZ_send_list;
00134       data_org_ = (int*) AZ_allocate(dataOrgLength * sizeof(int));
00135       for(i=0; i<dataOrgLength; i++) data_org_[i] = src.data_org_[i];
00136 
00137       Amat_->data_org = data_org_;
00138 
00139       int extLength = src.data_org_[AZ_N_ext_blk];
00140       external_ = (int*) AZ_allocate(extLength * sizeof(int));
00141       extern_index_ = (int*) AZ_allocate(extLength * sizeof(int));
00142       for(i=0; i<extLength; i++) {
00143          external_[i] = src.external_[i];
00144          extern_index_[i] = src.extern_index_[i];
00145       }
00146 
00147       update_index_ = (int*) AZ_allocate(N_update_ * sizeof(int));
00148       orderingUpdate_ = new int[N_update_];
00149       for(i=0; i<N_update_; i++) {
00150          update_index_[i] = src.update_index_[i];
00151          orderingUpdate_[i] = src.orderingUpdate_[i];
00152       }
00153 
00154       int cpntrLength = N_update_ + src.data_org_[AZ_N_ext_blk] + 1;
00155       Amat_->cpntr = (int*) AZ_allocate(cpntrLength * sizeof(int));
00156       for(i=0; i<cpntrLength; i++) Amat_->cpntr[i] = src.Amat_->cpntr[i];
00157    }
00158 }
00159 
00160 //==============================================================================
00161 AztecDVBR_Matrix::~AztecDVBR_Matrix(){
00162 
00163    if (isAllocated()) {
00164       delete [] Amat_->val;
00165       delete [] Amat_->bindx;
00166       delete [] Amat_->rpntr;
00167       delete [] Amat_->bpntr;
00168       delete [] Amat_->indx;
00169 
00170       delete [] remoteInds_;
00171       delete [] remoteBlockSizes_;
00172 
00173       setAllocated(false);
00174    }
00175 
00176    if (isLoaded()) {
00177       free(Amat_->cpntr);
00178       free(external_);
00179       free(extern_index_);
00180       free(update_index_);
00181       free(data_org_);
00182       delete [] orderingUpdate_;
00183 
00184       setLoaded(false);
00185    }
00186 
00187    delete [] nnzPerRow_;
00188    localNNZ_ = 0;
00189 
00190    AZ_matrix_destroy(&Amat_);
00191    Amat_ = NULL;
00192 }
00193 
00194 //==============================================================================
00195 int AztecDVBR_Matrix::getNumBlocksPerRow(int blkRow, int& nnzBlksPerRow) const {
00196 //
00197 //On return, nnzBlksPerRow will be the number of nonzero blocks in row blkRow.
00198 //
00199    if (!isAllocated()) return(1);
00200 
00201    int index;
00202 
00203    if (!inUpdate(blkRow, index)) {
00204       FEI_CERR << "AztecDVBR_Matrix::getNumBlocksPerRow: ERROR: blkRow "
00205            << blkRow << " not in local update list." << FEI_ENDL;
00206       return(1);
00207    }
00208 
00209    nnzBlksPerRow = Amat_->bpntr[index+1] - Amat_->bpntr[index];
00210 
00211    return(0);
00212 }
00213 
00214 //==============================================================================
00215 int AztecDVBR_Matrix::getNumBlocksPerRow(int* nnzBlksPerRow) const {
00216 //
00217 //nnzBlksPerRow must be allocated by the calling code.
00218 //
00219 //nnzBlksPerRow is a list of length number-of-local-block-rows, and
00220 //nnzBlksPerRow[i] gives the number of nonzeros blocks in row i.
00221 //
00222    if (!isAllocated()) return(1);
00223 
00224    for(int i=0; i<amap_->getNumLocalBlocks(); i++) {
00225       nnzBlksPerRow[i] = Amat_->bpntr[i+1] - Amat_->bpntr[i];
00226    }
00227 
00228    return(0);
00229 }
00230 
00231 //==============================================================================
00232 int AztecDVBR_Matrix::getNumNonzerosPerRow(int blkRow, int& nnzPerRow) const {
00233 //
00234 //This function finds nnzPerRow, the number of nonzero *point* entries for
00235 //row 'blkRow'.
00236 //
00237    if (!isAllocated()) return(1);
00238 
00239    int index;
00240 
00241    if (!inUpdate(blkRow, index)) {
00242       FEI_CERR << "AztecDVBR_Matrix::getNumNonzerosPerRow: ERROR: blkRow "
00243            << blkRow << " not in local update list." << FEI_ENDL;
00244       return(1);
00245    }
00246 
00247    nnzPerRow = nnzPerRow_[index];
00248 
00249    return(0);
00250 }
00251 
00252 //==============================================================================
00253 int AztecDVBR_Matrix::getNumNonzerosPerRow(int* nnzPerRow) const {
00254 //
00255 //nnzPerRow must be allocated by the calling code,
00256 //length number-of-local-block-rows.
00257 //
00258 //This function fills nnzPerRow so that nnzPerRow[i] gives the
00259 //number of nonzero *point* entries for row i.
00260 //
00261    if (!isAllocated()) return(1);
00262 
00263    for(int i=0; i<amap_->getNumLocalBlocks(); i++) {
00264       nnzPerRow[i] = nnzPerRow_[i];
00265    }
00266 
00267    return(0);
00268 }
00269 
00270 //==============================================================================
00271 int AztecDVBR_Matrix::getBlockSize(int blkRow, int blkCol,
00272                                    int& ptRows, int& ptCols) {
00273    int index;
00274 
00275    ptRows = 0;
00276    ptCols = 0;
00277 
00278    if (!inUpdate(blkRow, index)) {
00279       FEI_CERR << "AztecDVBR_Matrix::getBlockSize: ERROR: blkRow "
00280            << blkRow << " not in local update list." << FEI_ENDL;
00281       return(1);
00282    }
00283 
00284    ptRows = Amat_->rpntr[index+1] - Amat_->rpntr[index];
00285 
00286    int local = inUpdate(blkCol, index);
00287 
00288    if (local) {
00289       ptCols = Amat_->rpntr[index+1] - Amat_->rpntr[index];
00290    }
00291    else {
00292       index = AZ_find_index(blkCol, remoteInds_, numRemoteBlocks_);
00293 
00294       if (index < 0) return(1);
00295       else ptCols = remoteBlockSizes_[index];
00296    }
00297 
00298    return(0);
00299 }
00300 
00301 //==============================================================================
00302 void AztecDVBR_Matrix::matvec(const Aztec_LSVector& x, Aztec_LSVector& y) const {
00303       
00304 // AztecDVBR_Matrix::matvec --- form y = Ax
00305 
00306     assert(isLoaded());
00307 
00308     int *proc_config = amap_->getProcConfig();
00309     double *b = (double*)x.startPointer();
00310     double *c = (double*)y.startPointer();
00311 
00312     AZ_VBR_matvec_mult(b, c, Amat_, proc_config);
00313 
00314     return;
00315 }
00316 
00317 //==============================================================================
00318 void AztecDVBR_Matrix::put(double s){
00319 
00320    if (!isAllocated()) return;
00321 
00322    for(int i=0; i<localNNZ_; i++) {
00323       Amat_->val[i] = s;
00324    }
00325 
00326    return;
00327 }
00328 
00329 //==============================================================================
00330 int AztecDVBR_Matrix::getBlockRow(int blkRow,
00331                                   double* val,
00332                                   int* blkColInds,
00333                                   int numNzBlks) const {
00334 
00335    if (!isAllocated()) return(1);
00336 
00337    int index;
00338 
00339    if (!inUpdate(blkRow, index)) {
00340       FEI_CERR << "AztecDVBR_Matrix::getBlockRow: ERROR: blkRow "
00341            << blkRow << " not in local update list." << FEI_ENDL;
00342       return(1);
00343    }
00344 
00345    //for each block, we need to find its block column index
00346    //in the bindx array, then go to that same position in the indx
00347    //array to find out how many point-entries are in that block.
00348    //We can then use the indx entry to go to the val array and get
00349    //the data.
00350 
00351    int nnzBlks = 0, nnzPts = 0;
00352    int err = getNumBlocksPerRow(blkRow, nnzBlks);
00353    if (err) return(err);
00354    err = getNumNonzerosPerRow(blkRow, nnzPts);
00355    if (err) return(err);
00356    
00357    if (numNzBlks != nnzBlks) return(1);
00358 
00359    int offset = 0;
00360    int blkCounter = 0;
00361    const int* blkUpdate = amap_->getBlockUpdate();
00362    for(int indb = Amat_->bpntr[index]; indb<Amat_->bpntr[index+1]; indb++) {
00363 
00364       int numEntries = Amat_->indx[indb+1] - Amat_->indx[indb];
00365       int valOffset = Amat_->indx[indb];
00366 
00367       if (isLoaded()) {
00368          int ind = Amat_->bindx[indb];
00369          if (ind < N_update_) {
00370             blkColInds[blkCounter++] = blkUpdate[orderingUpdate_[ind]];
00371          }
00372          else {
00373             blkColInds[blkCounter++] = external_[ind-N_update_];
00374          }
00375       }
00376       else {
00377          blkColInds[blkCounter++] = Amat_->bindx[indb];
00378       }
00379 
00380       //ok, now we're ready to get the stuff.
00381       for(int i=0; i<numEntries; i++) {
00382          val[offset + i] = Amat_->val[valOffset + i];
00383       }
00384 
00385       offset += numEntries;
00386    }
00387 
00388    return(0);
00389 }
00390 
00391 //==============================================================================
00392 int AztecDVBR_Matrix::putBlockRow(int blkRow,
00393                                    double* val,
00394                                    int* blkColInds,
00395                                    int numNzBlks) const {
00396 
00397    if (!isAllocated()) return(1);
00398 
00399    int index;
00400 
00401    if (!inUpdate(blkRow, index)) {
00402       FEI_CERR << "AztecDVBR_Matrix::putBlockRow: ERROR: blkRow "
00403            << blkRow << " not in local update list." << FEI_ENDL;
00404       return(1);
00405    }
00406 
00407    //for each incoming block, we need to find its block column index
00408    //in the bindx array, then go to that same position in the indx
00409    //array to find out how many (point) entries are in that block.
00410    //We can then use the indx entry to go to the val array and store
00411    //the data.
00412 
00413    int offset = 0;
00414    for(int blk = 0; blk<numNzBlks; blk++) {
00415       int indb = getBindxOffset(blkColInds[blk],
00416                                 Amat_->bpntr[index], Amat_->bpntr[index+1]-1);
00417 
00418       if (indb < 0) messageAbort("putBlockRow: blk col not found in row.");
00419 
00420       int numEntries = Amat_->indx[indb+1] - Amat_->indx[indb];
00421       int valOffset = Amat_->indx[indb];
00422 
00423       //ok, now we're ready to store the stuff.
00424       for(int i=0; i<numEntries; i++) {
00425          Amat_->val[valOffset + i] = val[offset + i];
00426       }
00427 
00428       offset += numEntries;
00429    }
00430 
00431    return(0);
00432 }
00433 
00434 //==============================================================================
00435 int AztecDVBR_Matrix::sumIntoBlockRow(int blkRow,
00436                                        double* val,
00437                                        int* blkColInds,
00438                                        int numNzBlks) const
00439 {
00440   //
00441   //This function is the same as putBlockRow, except the values
00442   //are summed into any existing values rather than overwriting
00443   //them.
00444   //
00445   if (!isAllocated()) return(1);
00446 
00447   int index;
00448 
00449   if (!inUpdate(blkRow, index)) {
00450     FEI_CERR << "AztecDVBR_Matrix::sumIntoBlockRow: ERROR: blkRow "
00451    << blkRow << " not in local update list." << FEI_ENDL;
00452     return(1);
00453   }
00454 
00455   //for each incoming block, we need to find its block column index
00456   //in the bindx array, then go to that same position in the indx
00457   //array to find out how many (point) entries are in that block.
00458   //We can then use the indx entry to go to the val array and store
00459   //the data.
00460 
00461   int offset = 0;
00462 
00463   for(int blk = 0; blk<numNzBlks; blk++) {
00464     int indb = getBindxOffset(blkColInds[blk],
00465             Amat_->bpntr[index], Amat_->bpntr[index+1]-1);
00466 
00467     if (indb < 0) {
00468       FEI_CERR << "AztecDVBR_Matrix::sumIntoBlockRow: blk col "
00469      << blkColInds[blk] << " not found in row " << blkRow << FEI_ENDL;
00470       abort();
00471     }
00472 
00473     int numEntries = Amat_->indx[indb+1] - Amat_->indx[indb];
00474     int valOffset = Amat_->indx[indb];
00475 
00476     //ok, now we're ready to store the stuff.
00477     for(int i=0; i<numEntries; i++) {
00478       Amat_->val[valOffset + i] += val[offset + i];
00479     }
00480 
00481     offset += numEntries;
00482   }
00483 
00484   return(0);
00485 }
00486  
00487 //==============================================================================
00488 void AztecDVBR_Matrix::allocate(int* numNzBlks, int* blkColInds) {
00489 //
00490 // This function builds the structure of the matrix. i.e., does the
00491 // memory allocation, etc.
00492 //
00493 
00494    //calculate the bpntr array, which holds info about the number of
00495    //nonzero blocks per row.
00496    calcBpntr(numNzBlks);
00497 
00498    //we can now get the total number of nonzero blocks from the last
00499    //entry in bpntr.
00500    int totalNumNzBlks = Amat_->bpntr[N_update_];
00501 
00502    //now we can set the bindx array, which holds block column indices.
00503    setBindx(totalNumNzBlks, blkColInds);
00504 
00505    //and now we're ready to allocate and fill the indx array, which
00506    //holds info on the number of point entries in each nonzero block.
00507    calcIndx(totalNumNzBlks);
00508 
00509    //the last thing we need to do is allocate and initialize the val array.
00510    Amat_->val = new double[localNNZ_];
00511 
00512    for(int i=0; i<localNNZ_; i++) {
00513       Amat_->val[i] = 0.0;
00514    }
00515 
00516    setAllocated(true);
00517    return;
00518 }
00519 
00520 //==============================================================================
00521 void AztecDVBR_Matrix::loadComplete() {
00522 //
00523 // This is where we call the Aztec function AZ_transform, which calculates
00524 // communication parameters and re-orders the equations for use as a
00525 // global distributed matrix.
00526 //
00527    MPI_Comm thisComm = amap_->getCommunicator();
00528 
00529 // Sync processors.
00530 
00531    MPI_Barrier(thisComm);
00532 
00533 #ifndef FEI_SER
00534    int thisProc = 0;
00535    MPI_Comm_rank(thisComm, &thisProc);
00536 #endif
00537 
00538    AZ_transform(amap_->getProcConfig(), &external_, Amat_->bindx, Amat_->val,
00539                 amap_->getBlockUpdate(), &update_index_, &extern_index_, &data_org_,
00540                 N_update_, Amat_->indx, Amat_->bpntr, Amat_->rpntr,
00541                 &(Amat_->cpntr), AZ_VBR_MATRIX);
00542 
00543    data_org_[AZ_internal_use] = 1;
00544 
00545    Amat_->data_org = data_org_;
00546 
00547 //On return from AZ_transform, the array update_index_ contains a mapping
00548 //to the local re-ordering of the indices of the update_ array. Now we will
00549 //fill the orderingUpdate array with the reverse of that mapping. i.e., a
00550 //record of how to get back to the original ordering of the update indices.
00551 
00552    orderingUpdate_ = new int[N_update_];
00553    for(int ii=0; ii<N_update_; ii++)
00554       orderingUpdate_[update_index_[ii]] = ii;
00555 
00556 // Sync processors.
00557 #ifndef FEI_SER
00558    MPI_Barrier(thisComm);
00559 #endif
00560 
00561    setLoaded(true);
00562 
00563    return;
00564 }
00565 
00566 //==============================================================================
00567 bool AztecDVBR_Matrix::readFromFile(const char *filename){
00568 //
00569 //readFromFile should be able to be called after the matrix is constructed,
00570 //and before allocate has been called. i.e., example usage should include:
00571 //
00572 //  AztecDVBR_Matrix A(map, update);
00573 //  A.readFromFile(fileName);
00574 //  A.matvec(b, c);
00575 //
00576 //i.e., readFromFile can take the place of the allocate and loadComplete
00577 //calls.
00578 //
00579    FILE *infile = NULL;
00580    MPI_Comm thisComm = amap_->getCommunicator();
00581 
00582    MPI_Barrier(thisComm);
00583 
00584    infile = fopen(filename, "r");
00585    if (!infile) messageAbort("readFromFile: couldn't open file.");
00586 
00587    int* num_nz_blocks = NULL;
00588    int* blk_col_inds = NULL;
00589 
00590    readAllocateInfo(infile, num_nz_blocks, blk_col_inds);
00591 
00592    allocate(num_nz_blocks, blk_col_inds);
00593 
00594    delete [] num_nz_blocks;
00595    delete [] blk_col_inds;
00596 
00597    fclose(infile);
00598    infile = fopen(filename, "r");
00599     
00600    readMatrixData(infile);
00601 
00602    fclose(infile);
00603 
00604    loadComplete();
00605 
00606    return(true);
00607 }
00608 
00609 //==============================================================================
00610 void AztecDVBR_Matrix::readAllocateInfo(FILE* infile,
00611                                      int*& num_nz_blocks,
00612                                      int*& blk_col_inds) {
00613 //
00614 //This function will read through infile and construct the lists
00615 //num_nz_blocks (which is the number of nonzero blocks per row) and
00616 //blk_col_inds (which is the block-column indices of those blocks).
00617 //
00618 //It is assumed that these two lists are empty when this function is
00619 //called.
00620 
00621    int i;
00622 
00623    if (num_nz_blocks) delete [] num_nz_blocks;
00624    if (blk_col_inds) delete [] blk_col_inds;
00625 
00626    num_nz_blocks = new int[N_update_];
00627 
00628    //we'll use a 2-D array for constructing the set of block column indices,
00629    //because we need to keep them grouped by rows, and we aren't guaranteed
00630    //that they'll be grouped by rows in the file.
00631    int totalNumBlks = 0;
00632    int** blkColInds = new int*[N_update_];
00633 
00634    for(i=0; i<N_update_; i++) {
00635       num_nz_blocks[i] = 0;
00636       blkColInds[i] = NULL;
00637    }
00638 
00639    int blkRows, blkCols, rows, cols;
00640    char line[256];
00641 
00642    do {
00643       fgets(line,256,infile);
00644    } while(strchr(line,'%'));
00645    sscanf(line,"%d %d %d %d",&blkRows, &blkCols, &rows, &cols);
00646 
00647    if ((blkRows != blkCols) || (rows != cols))
00648       messageAbort("readAllocateInfo: non-square matrix not allowed.");
00649 
00650    int br, bc, pr, pc, index;
00651 
00652    while (!feof(infile)) {
00653       do {
00654          fgets(line,256,infile);
00655       } while(strchr(line,'%'));
00656 
00657       if(feof(infile))break;
00658 
00659       sscanf(line, "%d %d %d %d", &br, &bc, &pr, &pc);
00660 
00661       if (inUpdate(br, index)) {
00662          if ((bc < 0) || bc >= blkCols) {
00663             char mesg[80];
00664             sprintf(mesg,"readAllocateInfo: blkCols %d, 0-based col ind %d",
00665                     blkCols, bc);
00666             fclose(infile);
00667             messageAbort(mesg);
00668          }
00669          insertList(bc, blkColInds[index], num_nz_blocks[index]);
00670          totalNumBlks++;
00671       }
00672    }
00673 
00674    //so we've read the whole file, now flatten the 2-D list blkColInds
00675    //into the required 1-D list blk_col_inds.
00676    blk_col_inds = new int[totalNumBlks];
00677 
00678    int offset = 0;
00679    for(i=0; i<N_update_; i++) {
00680       for(int j=0; j<num_nz_blocks[i]; j++) {
00681          blk_col_inds[offset++] = blkColInds[i][j];
00682       }
00683 
00684       delete [] blkColInds[i];
00685    }
00686 
00687    delete [] blkColInds;
00688 }
00689 
00690 //==============================================================================
00691 void AztecDVBR_Matrix::readMatrixData(FILE* infile) {
00692 
00693    int blkRows, blkCols, rows, cols;
00694    int br, bc, pr, pc, nnz, index;
00695    double* blockValues = NULL;
00696    char line[256];
00697 
00698    do {
00699       fgets(line,256,infile);
00700    } while(strchr(line,'%'));
00701    sscanf(line,"%d %d %d %d",&blkRows, &blkCols, &rows, &cols);
00702 
00703    while (!feof(infile))  {
00704       do {
00705          fgets(line,256,infile);
00706       } while(strchr(line,'%'));
00707 
00708       if(feof(infile))break;
00709 
00710       sscanf(line, "%d %d %d %d %d", &br, &bc, &pr, &pc, &nnz);
00711 
00712       if (inUpdate(br, index)){
00713          //br (block-row) is in the local row-space, so let's
00714          //plug this block into our matrix data structure.
00715 
00716          blockValues = new double[nnz];
00717          getValuesFromString(line, std::strlen(line)+1, blockValues, nnz);
00718 
00719          putBlockRow(br, blockValues, &bc, 1);
00720 
00721          delete [] blockValues;
00722       }
00723    }
00724 }
00725 
00727 void AztecDVBR_Matrix::getValuesFromString(char *line, int len, double *values,
00728                                            int lenValues){
00729     (void)len;
00730     int i;
00731 
00732 //first, we know that 'line' contains 5 integers at the beginning, separated
00733 //by spaces, which we want to jump over. So we need the offset of the 5th
00734 //space in 'line'.
00735 
00736     char *offset = &line[0];
00737     for(i=0; i<5; i++){
00738         offset = strchr(offset, ' ');
00739         offset++;
00740     }
00741 
00742 //now we're ready to pick out the numbers to put into the values array.
00743     for(i=0; i<lenValues; i++){
00744         sscanf(offset,"%le",&(values[i]));
00745         offset = strchr(offset, ' ');
00746         offset++;
00747     }
00748 
00749     return;
00750 }
00751 
00753 bool AztecDVBR_Matrix::writeToFile(const char *fileName) const {
00754 
00755    int thisProc = amap_->getProcConfig()[AZ_node];
00756    int numProcs = amap_->getProcConfig()[AZ_N_procs];
00757    MPI_Comm thisComm = amap_->getCommunicator();
00758 
00759    int numGlobalBlocks = amap_->getNumGlobalBlocks();
00760    int numGlobalPtEqns = amap_->globalSize();
00761 
00762    int numLocalBlocks = N_update_;
00763 
00764    FILE* file = NULL;
00765 
00766    for(int p=0; p<numProcs; p++) {
00767       MPI_Barrier(thisComm);
00768 
00769       if (p == thisProc) {
00770          if (thisProc == 0) {
00771             //open the file for writing and write the first line, which is
00772             //num-blk-rows num-block-cols num-pt-rows num-pt-cols
00773 
00774             file = fopen(fileName, "w");
00775             fprintf(file, "%d %d %d %d\n", numGlobalBlocks, numGlobalBlocks,
00776                                            numGlobalPtEqns, numGlobalPtEqns);
00777          }
00778          else {
00779             //open the file for appending
00780             file = fopen(fileName, "a");
00781          }
00782 
00783          //now loop over the local portion of the matrix, writing it to file,
00784          //one nonzer-block per line. Each line of the file will contain these
00785          //numbers separated by spaces:
00786          //blk-row-index
00787          //blk-col-index
00788          //blk-row-size (number of pt-rows in this block-row)
00789          //nnz          (number of nonzero point-entries in this block-entry)
00790          //nonzero1 nonzero2 ... nonzero<nnz> (the nonzeros for this block)
00791 
00792          for(int brow=0; brow<numLocalBlocks; brow++) {
00793             int bcolind1 = Amat_->bpntr[brow];
00794             int bcolind2 = Amat_->bpntr[brow+1];
00795 
00796             for(int ind=bcolind1; ind<bcolind2; ind++) {
00797                int nnzPts = Amat_->indx[ind+1] - Amat_->indx[ind];
00798 
00799                if (nnzPts <= 0) continue;
00800 
00801                int blkRowSize = Amat_->rpntr[brow+1]-Amat_->rpntr[brow];
00802 
00803                int globCol = -1;
00804                int lookup = Amat_->bindx[ind];
00805                if (isLoaded()) {
00806                   if (lookup < N_update_) {
00807                      globCol = amap_->getBlockUpdate()[orderingUpdate_[lookup]];
00808                   }
00809                   else {
00810                      globCol = external_[lookup-N_update_];
00811                   }
00812                }
00813                else {
00814                   globCol = lookup;
00815                }
00816 
00817                int globalRow = amap_->getBlockUpdate()[brow];
00818                if (isLoaded()) globalRow = amap_->getBlockUpdate()[orderingUpdate_[brow]];
00819 
00820                fprintf(file, "%d %d %d %d ", globalRow, globCol,
00821                                           blkRowSize, nnzPts);
00822 
00823                int offset = Amat_->indx[ind];
00824                for(int i=0; i<nnzPts; i++) {
00825                   fprintf(file, "%20.13e ", Amat_->val[offset+i]);
00826                }
00827                fprintf(file, "\n");
00828             }
00829          }
00830 
00831          fclose(file);
00832       }
00833    }
00834 
00835    return(true);
00836 }
00837 
00838 //==============================================================================
00839 int AztecDVBR_Matrix::inUpdate(int globalIndex, int& localIndex) const {
00840 //
00841 // This function determines whether globalIndex is in the local update set,
00842 // and if it is, returns in localIndex the local index for it. If update_index_
00843 // has already been allocated and set (by AZ_transform) then localIndex is
00844 // taken from there.
00845 // If globalIndex is not in the update set, inUpdate returns 0.
00846 //
00847     localIndex = AZ_find_index(globalIndex, amap_->getBlockUpdate(), N_update_);
00848 
00849     if(localIndex==-1)return(0);
00850 
00851     if(isLoaded_){
00852         localIndex = update_index_[localIndex];
00853     }
00854 
00855     return(1);
00856 }
00857 
00858 //==============================================================================
00859 void AztecDVBR_Matrix::calcRpntr() {
00860 //
00861 //This function will use information from the Aztec_BlockMap 'amap_'
00862 //to set the Amat_->rpntr array.
00863 //
00864 //rpntr[0..M] (where M = number-of-blocks)
00865 //rpntr[0] = 0
00866 //rpntr[k+1] - rpntr[k] = size of block k
00867 //
00868    const int* blkSizes = amap_->getBlockSizes();
00869 
00870    Amat_->rpntr = new int[N_update_+1];
00871 
00872    Amat_->rpntr[0] = 0;
00873 
00874    for(int i=0; i<N_update_; i++) {
00875       Amat_->rpntr[i+1] = Amat_->rpntr[i] + blkSizes[i];
00876       if (blkSizes[i] < 0)
00877           messageAbort("allocate: negative block size.");
00878    }
00879 }
00880 
00881 //==============================================================================
00882 void AztecDVBR_Matrix::calcBpntr(int* numNzBlks) {
00883 //
00884 //This function will set the Amat_->bpntr array.
00885 //
00886 //bpntr[0..M] (where M = number-of-blocks)
00887 //bpntr[0] = 0
00888 //bpntr[k+1]-bpntr[k] = number of nonzero blocks in row k
00889 //
00890    Amat_->bpntr = new int[N_update_+1];
00891 
00892    Amat_->bpntr[0] = 0;
00893 
00894    for(int i=0; i<N_update_; i++) {
00895       Amat_->bpntr[i+1] = Amat_->bpntr[i] + numNzBlks[i];
00896    }
00897 }
00898 
00899 //==============================================================================
00900 void AztecDVBR_Matrix::setBindx(int nnzBlks, int* blkColInds) {
00901 //
00902 //This function simply allocates and fills the Amat_->bindx array.
00903 //
00904    Amat_->bindx = new int[nnzBlks];
00905 
00906    for(int i=0; i<nnzBlks; i++) {
00907       Amat_->bindx[i] = blkColInds[i];
00908       if (blkColInds[i] < 0)
00909          messageAbort("setBindx: negative block col index.");
00910    }
00911 }
00912 
00913 //==============================================================================
00914 void AztecDVBR_Matrix::messageAbort(const char* mesg) const {
00915    FEI_CERR << "AztecDVBR_Matrix: ERROR: " << mesg << " Aborting." << FEI_ENDL;
00916    abort();
00917 }
00918 
00919 //==============================================================================
00920 void AztecDVBR_Matrix::calcIndx(int nnzBlks) {
00921 //
00922 //This function allocates and fills the Amat_->indx array, which holds info
00923 //on the number of entries in each nonzero block.
00924 //
00925 //indx[0..bpntr[M]], (where M = number of local block rows)
00926 //indx[0] = 0
00927 //indx[k+1]-indx[k] = number of entries in nonzero block k
00928 //
00929 
00930    Amat_->indx = new int[nnzBlks+1];
00931 
00932    //we need to obtain block sizes for all local nonzero blocks. rpntr
00933    //gives us the sizes for the blocks with column indices in the local
00934    //update set, but we'll have to do some message passing to obtain the
00935    //sizes of blocks with column indices in other procs' update sets.
00936 
00937    int numProcs = amap_->getProcConfig()[AZ_N_procs];
00938 
00939    if (numProcs > 1) {
00940       //form a list of the column indices that are not local.
00941       calcRemoteInds(remoteInds_, numRemoteBlocks_);
00942 
00943       //now get sizes of blocks that correspond to remote rows.
00944       remoteBlockSizes_ = new int[numRemoteBlocks_];
00945       getRemoteBlkSizes(remoteBlockSizes_, remoteInds_, numRemoteBlocks_);
00946    }
00947 
00948    //now we're ready to set the block sizes in Amat_->indx.
00949    int index;
00950 
00951    Amat_->indx[0] = 0;
00952 
00953    for(int i=0; i<amap_->getNumLocalBlocks(); i++) {
00954       int rowBlkSize = Amat_->rpntr[i+1] - Amat_->rpntr[i];
00955 
00956       int colStart = Amat_->bpntr[i];
00957       int colEnd = Amat_->bpntr[i+1] - 1;
00958 
00959       for(int j=colStart; j<=colEnd; j++) {
00960          if (inUpdate(Amat_->bindx[j], index)) {
00961             int colBlkSize = Amat_->rpntr[index+1] - Amat_->rpntr[index];
00962 
00963             Amat_->indx[j+1] = Amat_->indx[j] + rowBlkSize*colBlkSize;
00964          }
00965          else { //it's a remoteIndex
00966             if (numProcs == 1) {
00967                char mesg[80];
00968                sprintf(mesg,"calcIndx: blk col index %d not in update set.",
00969                        Amat_->bindx[j]);
00970                messageAbort(mesg);
00971             }
00972 
00973             index = AZ_find_index(Amat_->bindx[j], remoteInds_,
00974                                   numRemoteBlocks_);
00975             if (index >= 0) {
00976                Amat_->indx[j+1] = Amat_->indx[j] +  
00977                                 rowBlkSize*remoteBlockSizes_[index];
00978             }
00979             else { //if it wasn't in update or remoteInds, then panic!
00980                messageAbort("calcIndx: block column index not found.");
00981             }
00982          }
00983       } // end for j loop
00984 
00985       nnzPerRow_[i] = Amat_->indx[colEnd+1] - Amat_->indx[colStart];
00986    } // end for i loop
00987 
00988    localNNZ_ = Amat_->indx[nnzBlks];
00989 }
00990 
00991 //==============================================================================
00992 int AztecDVBR_Matrix::getBindxOffset(int blkInd,
00993                                      int bpntrStart, int bpntrEnd) const {
00994 //
00995 //This function returns the index of blkInd in the bindx array,
00996 //searching positions bindx[bpntrStart..bpntrEnd].
00997 //
00998 //If blkInd is not found, -1 is returned.
00999 //
01000    for(int i=bpntrStart; i<=bpntrEnd; i++) {
01001       int ind = Amat_->bindx[i];
01002       int globalCol = -1;
01003       if (isLoaded()) {
01004          if (ind < N_update_) {
01005             globalCol = amap_->getBlockUpdate()[orderingUpdate_[ind]];
01006          }
01007          else {
01008             globalCol = external_[ind-N_update_];
01009          }
01010       }
01011       else globalCol = ind;
01012 
01013       if (globalCol == blkInd) return(i);
01014    }
01015 
01016    return(-1);
01017 }
01018 
01019 //==============================================================================
01020 void AztecDVBR_Matrix::calcRemoteInds(int*& remoteInds, int& len) {
01021 //
01022 //Form a list of the block column indices that are not in the local
01023 //update set.
01024 //
01025    int nnzBlks = Amat_->bpntr[amap_->getNumLocalBlocks()];
01026    int local;
01027 
01028    for(int i=0; i<nnzBlks; i++) {
01029       if (!inUpdate(Amat_->bindx[i], local)) {
01030          insertList(Amat_->bindx[i], remoteInds, len);
01031       }
01032    }
01033 }
01034 
01035 //==============================================================================
01036 void AztecDVBR_Matrix::getRemoteBlkSizes(int* remoteBlkSizes,
01037            int* remoteInds,
01038                                          int len)
01039 {
01040   //
01041   //remoteInds is a sorted list of indices that correspond to rows
01042   //in remote processors' update lists. This function will spread the
01043   //indices to all processors so that they can provide the blk sizes,
01044   //then spread that information back to all processors.
01045   //
01046 #ifdef FEI_SER
01047   return;
01048 #else
01049   int numProcs = amap_->getProcConfig()[AZ_N_procs];
01050   int thisProc = amap_->getProcConfig()[AZ_node];
01051   MPI_Comm comm = amap_->getCommunicator();
01052 
01053    int* lengths = new int[numProcs];
01054    lengths[0] = 0;
01055 
01056    //gather up the lengths of the lists that each proc will be sending.
01057    MPI_Allgather(&len, 1, MPI_INT, lengths, 1, MPI_INT, comm);
01058 
01059    //now form a list of the offset at which each proc's contribution will
01060    //be placed in the all-gathered list.
01061    int* offsets = new int[numProcs];
01062 
01063    offsets[0] = 0;
01064    int totalLength = lengths[0];
01065    for(int i=1; i<numProcs; i++) {
01066       offsets[i] = offsets[i-1] + lengths[i-1];
01067       totalLength += lengths[i];
01068    }
01069 
01070    //now we can allocate the list to recv into.
01071    int* recvBuf = new int[totalLength];
01072 
01073    //now we're ready to do the gather.
01074    MPI_Allgatherv(remoteInds, len, MPI_INT, recvBuf, lengths, offsets,
01075                   MPI_INT, comm);
01076 
01077    //now we'll run through the list and put block sizes into a list of
01078    //the same length as the total recvBuf list.
01079    int* blkSizes = new int[totalLength];
01080    int index;
01081 
01082    for(int j=0; j<totalLength; j++) {
01083       if (inUpdate(recvBuf[j], index)) {
01084          blkSizes[j] = Amat_->rpntr[index+1]-Amat_->rpntr[index];
01085       }
01086       else blkSizes[j] = 0;
01087    }
01088 
01089    //now we'll reduce this info back onto all processors. We'll use MPI_SUM.
01090    //Since the sizes we did NOT supply hold a 0, and each spot in the list
01091    //should only have a nonzero size from 1 processor, the result will be
01092    //that each spot in the result list has the correct value.
01093    int* recvSizes = new int[totalLength];
01094 
01095    MPI_Allreduce(blkSizes, recvSizes, totalLength, MPI_INT, MPI_SUM, comm);
01096 
01097    //and finally, we just need to run our section of the list of recv'd sizes,
01098    //and transfer them into the remoteBlkSizes list.
01099    int offset = offsets[thisProc];
01100    for(int k=0; k<len; k++) {
01101       remoteBlkSizes[k] = recvSizes[offset + k];
01102       if (recvSizes[offset+k] <= 0)
01103          messageAbort("getRemoteBlkSizes: recvd a size <= 0.");
01104    }
01105 
01106    delete [] lengths;
01107    delete [] offsets;
01108    delete [] recvBuf;
01109    delete [] blkSizes;
01110    delete [] recvSizes;
01111 #endif
01112 }
01113 
01114 //==============================================================================
01115 void AztecDVBR_Matrix::insertList(int item, int*& list, int& len) {
01116 //
01117 //insert 'item' in 'list', if it's not already in there,
01118 //and update the list's length, 'len'.
01119 //
01120 //We want to keep the list ordered, so we'll insert item in
01121 //the list after the biggest existing entry that's smaller, and
01122 //before the smallest existing entry that's bigger.
01123 //
01124 
01125    if (len <= 0) {
01126       list = new int[1];
01127       list[0] = item;
01128       len = 1;
01129       return;
01130    }
01131 
01132    int index = AZ_find_index(item, list, len);
01133 
01134    if (index >= 0) return;
01135 
01136    int* newList = new int[len+1];
01137 
01138    //bring over the contents of the old list, putting in the new
01139    //one at the appropriate point.
01140    int inserted = 0;
01141    for(int i=0; i<len; i++) {
01142       if (!inserted) {
01143          if (list[i] < item) newList[i] = list[i];
01144          else {
01145             newList[i] = item;
01146             inserted = 1;
01147          }
01148       }
01149       else newList[i] = list[i-1];
01150    }
01151 
01152    //now put in the last list entry
01153    if (inserted) newList[len] = list[len-1];
01154    else newList[len] = item;
01155 
01156    //delete the old memory and reset the pointer.
01157    if (len > 0) delete [] list;
01158    list = newList;
01159 
01160    //update the length.
01161    len++;
01162 }
01163 
01164 }//namespace fei_trilinos
01165 
01166 #endif
01167 //HAVE_FEI_AZTECOO
01168 
 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