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

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