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