Epetra Package Browser (Single Doxygen Collection) Development
Epetra_Import_Util.cpp
Go to the documentation of this file.
00001 /*
00002 //@HEADER
00003 // ************************************************************************
00004 // 
00005 //               Epetra: Linear Algebra Services Package 
00006 //                 Copyright 2011 Sandia Corporation
00007 // 
00008 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00009 // the U.S. Government retains certain rights in this software.
00010 //
00011 // Redistribution and use in source and binary forms, with or without
00012 // modification, are permitted provided that the following conditions are
00013 // met:
00014 //
00015 // 1. Redistributions of source code must retain the above copyright
00016 // notice, this list of conditions and the following disclaimer.
00017 //
00018 // 2. Redistributions in binary form must reproduce the above copyright
00019 // notice, this list of conditions and the following disclaimer in the
00020 // documentation and/or other materials provided with the distribution.
00021 //
00022 // 3. Neither the name of the Corporation nor the names of the
00023 // contributors may be used to endorse or promote products derived from
00024 // this software without specific prior written permission.
00025 //
00026 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00027 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00028 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00029 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00030 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00031 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00032 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00033 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00034 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00035 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00036 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00037 //
00038 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00039 // 
00040 // ************************************************************************
00041 //@HEADER
00042 */
00043 
00044 #include "Epetra_ConfigDefs.h"
00045 #include "Epetra_Import_Util.h"
00046 #include "Epetra_Util.h"
00047 #include "Epetra_Comm.h"
00048 #include "Epetra_BlockMap.h"
00049 #include "Epetra_Map.h"
00050 #include "Epetra_Import.h"
00051 #include "Epetra_CrsMatrix.h"
00052 #include "Epetra_HashTable.h"
00053 
00054 #include <stdexcept>
00055 
00056 
00057 namespace Epetra_Import_Util { 
00058 
00059 // =========================================================================
00060 // =========================================================================
00061 template<typename int_type>
00062 int TPackAndPrepareWithOwningPIDs(const Epetra_CrsMatrix & A, 
00063          int NumExportIDs,
00064          int * ExportLIDs,
00065          int & LenExports,
00066          char *& Exports,
00067          int & SizeOfPacket,
00068          int * Sizes,
00069          bool & VarSizes,
00070          std::vector<int>& pids)
00071 {
00072   int i, j;
00073 
00074   VarSizes = true; //enable variable block size data comm
00075 
00076   int TotalSendLength = 0;
00077   int * IntSizes = 0; 
00078   if( NumExportIDs>0 ) IntSizes = new int[NumExportIDs];
00079 
00080   int SizeofIntType = sizeof(int_type);
00081 
00082   for(i = 0; i < NumExportIDs; ++i) {
00083     int NumEntries;
00084     A.NumMyRowEntries( ExportLIDs[i], NumEntries );
00085     // Will have NumEntries doubles, 2*NumEntries +2 ints pack them interleaved     Sizes[i] = NumEntries;
00086     // NTS: We add the owning PID as the SECOND int of the pair for each entry
00087     Sizes[i] = NumEntries;
00088     // NOTE: Mixing and matching Int Types would be more efficient, BUT what about variable alignment?
00089     IntSizes[i] = 1 + (((2*NumEntries+2)*SizeofIntType)/(int)sizeof(double));
00090     TotalSendLength += (Sizes[i]+IntSizes[i]);
00091   }    
00092   
00093   double * DoubleExports = 0; 
00094   SizeOfPacket = (int)sizeof(double);
00095        
00096   //setup buffer locally for memory management by this object
00097   if( TotalSendLength*SizeOfPacket > LenExports ) {
00098     if( LenExports > 0 ) delete [] Exports;
00099     LenExports = TotalSendLength*SizeOfPacket;
00100     DoubleExports = new double[TotalSendLength];
00101     for( int i = 0; i < TotalSendLength; ++i ) DoubleExports[i] = 0.0;
00102     Exports = (char *) DoubleExports;
00103   } 
00104   
00105   int NumEntries;
00106   double * values;
00107   double * valptr, * dintptr; 
00108 
00109   // Each segment of Exports will be filled by a packed row of information for each row as follows:
00110   // 1st int: GRID of row where GRID is the global row ID for the source matrix
00111   // next int:  NumEntries, Number of indices in row
00112   // next 2*NumEntries: The actual indices and owning [1] PID each for the row in (GID,PID) pairs with the GID first.
00113 
00114   // [1] Owning is defined in the sense of "Who owns the GID in the DomainMap," aka, who sends the GID in the importer
00115 
00116   const Epetra_Map & rowMap = A.RowMap();
00117   const Epetra_Map & colMap = A.ColMap();
00118   
00119   if( NumExportIDs > 0 ) {
00120     int_type * Indices;
00121     int_type FromRow; 
00122     int_type * intptr;                         
00123     
00124     int maxNumEntries = A.MaxNumEntries();
00125     std::vector<int> MyIndices(maxNumEntries);
00126     dintptr = (double *) Exports;
00127     valptr = dintptr + IntSizes[0];
00128     intptr = (int_type *) dintptr;
00129     for (i=0; i<NumExportIDs; i++) {
00130       FromRow   = (int_type) rowMap.GID64(ExportLIDs[i]);
00131       intptr[0] = FromRow;
00132       values    = valptr;
00133       Indices   = intptr + 2;
00134       EPETRA_CHK_ERR(A.ExtractMyRowCopy(ExportLIDs[i], maxNumEntries, NumEntries, values, &MyIndices[0]));
00135       for (j=0; j<NumEntries; j++) {
00136   Indices[2*j]   = (int_type) colMap.GID64(MyIndices[j]);   // convert to GIDs
00137   Indices[2*j+1] = pids[MyIndices[j]];                      // PID owning the entry.
00138       }
00139       intptr[1] = NumEntries; // Load second slot of segment
00140       if( i < (NumExportIDs-1) ) {
00141   dintptr += (IntSizes[i]+Sizes[i]);
00142   valptr = dintptr + IntSizes[i+1];
00143   intptr = (int_type *) dintptr;
00144       } 
00145     }    
00146 
00147     for(i = 0; i < NumExportIDs; ++i )
00148       Sizes[i] += IntSizes[i];
00149   }
00150   
00151   if( IntSizes ) delete [] IntSizes;
00152   
00153   return(0);
00154 }
00155 
00156 
00157 // =========================================================================
00158 // =========================================================================
00159 int PackAndPrepareWithOwningPIDs(const Epetra_CrsMatrix & SourceMatrix, 
00160          int NumExportIDs,
00161          int * ExportLIDs,
00162          int & LenExports,
00163          char *& Exports,
00164          int & SizeOfPacket,
00165          int * Sizes,
00166          bool & VarSizes,
00167          std::vector<int>& SourcePids) 
00168 {
00169   if(SourceMatrix.RowMap().GlobalIndicesInt()) {
00170     return TPackAndPrepareWithOwningPIDs<int>(SourceMatrix,NumExportIDs,ExportLIDs,LenExports,Exports,SizeOfPacket,Sizes,VarSizes,SourcePids);
00171   }
00172   else if(SourceMatrix.RowMap().GlobalIndicesLongLong()) {
00173     return TPackAndPrepareWithOwningPIDs<long long>(SourceMatrix,NumExportIDs,ExportLIDs,LenExports,Exports,SizeOfPacket,Sizes,VarSizes,SourcePids);
00174   }
00175   else
00176     throw std::runtime_error("UnpackWithOwningPIDsCount: Unable to determine source global index type");
00177 }
00178 
00179 
00180 // =========================================================================
00181 // =========================================================================
00182 template<typename int_type>
00183 int TUnpackWithOwningPIDsCount(const Epetra_CrsMatrix& SourceMatrix, 
00184              int NumSameIDs,
00185              int NumRemoteIDs,
00186              const int * RemoteLIDs,
00187              int NumPermuteIDs,
00188              const int *PermuteToLIDs,
00189              const int *PermuteFromLIDs,
00190              int LenImports,
00191              char* Imports)
00192 {  
00193   int i,nnz=0;
00194   int SizeofIntType = (int)sizeof(int_type);
00195 
00196   // SameIDs: Always first, always in the same place
00197   for(i=0; i<NumSameIDs; i++)
00198     nnz+=SourceMatrix.NumMyEntries(i);
00199   
00200   // PermuteIDs: Still local, but reordered
00201   for(i=0; i<NumPermuteIDs; i++)
00202     nnz += SourceMatrix.NumMyEntries(PermuteFromLIDs[i]);
00203   
00204   // RemoteIDs:  RemoteLIDs tells us the ID, we need to look up the length the hard way.  See UnpackAndCombine for where this code came from
00205   if(NumRemoteIDs > 0) {
00206     double * dintptr = (double *) Imports;
00207     // General version
00208     int_type *  intptr = (int_type *) dintptr;
00209     int   NumEntries = (int) intptr[1];
00210     int      IntSize = 1 + (((2*NumEntries+2)*SizeofIntType)/(int)sizeof(double));
00211     for(i=0; i<NumRemoteIDs; i++) {
00212       nnz += NumEntries;
00213       
00214       if( i < (NumRemoteIDs-1) ) {
00215   dintptr   += IntSize + NumEntries;
00216   intptr     = (int_type *) dintptr;
00217   NumEntries = (int) intptr[1];
00218   IntSize    = 1 + (((2*NumEntries+2)*SizeofIntType)/(int)sizeof(double));
00219       }
00220     }
00221   }
00222 
00223   return nnz;
00224 }
00225 
00226   
00227 // =========================================================================
00228 // =========================================================================
00229 int UnpackWithOwningPIDsCount(const Epetra_CrsMatrix& SourceMatrix, 
00230             int NumSameIDs,
00231             int NumRemoteIDs,
00232             const int * RemoteLIDs,
00233             int NumPermuteIDs,
00234             const int *PermuteToLIDs,
00235             const int *PermuteFromLIDs,
00236             int LenImports,
00237             char* Imports)  
00238 {
00239   if(SourceMatrix.RowMap().GlobalIndicesInt()) {
00240     return TUnpackWithOwningPIDsCount<int>(SourceMatrix,NumSameIDs,NumRemoteIDs,RemoteLIDs,NumPermuteIDs,
00241              PermuteToLIDs,PermuteFromLIDs,LenImports,Imports);
00242   }
00243   else if(SourceMatrix.RowMap().GlobalIndicesLongLong()) {
00244    return TUnpackWithOwningPIDsCount<long long>(SourceMatrix,NumSameIDs,NumRemoteIDs,RemoteLIDs,NumPermuteIDs,
00245             PermuteToLIDs,PermuteFromLIDs,LenImports,Imports);
00246   }
00247   else
00248     throw std::runtime_error("UnpackWithOwningPIDsCount: Unable to determine source global index type");
00249 
00250 }
00251 
00252         
00253 // =========================================================================
00254 // =========================================================================
00255 template<typename int_type>
00256 int TUnpackAndCombineIntoCrsArrays(const Epetra_CrsMatrix& SourceMatrix, 
00257           int NumSameIDs,
00258           int NumRemoteIDs,
00259           const int * RemoteLIDs,
00260           int NumPermuteIDs,
00261           const int *PermuteToLIDs,
00262           const int *PermuteFromLIDs,
00263           int LenImports,
00264           char* Imports,
00265           int TargetNumRows,
00266           int TargetNumNonzeros,
00267           int MyTargetPID,
00268           int * CSR_rowptr,
00269           int_type * CSR_colind,
00270           double * CSR_vals,
00271           const std::vector<int> &SourcePids,
00272           std::vector<int> &TargetPids)
00273 {
00274   // What we really need to know is where in the CSR arrays each row should start (aka the rowptr).
00275   // We do that by (a) having each row record it's size in the rowptr (b) doing a cumulative sum to get the rowptr values correct and 
00276   // (c) Having each row copied into the right colind / values locations.
00277 
00278   // From Epetra_CrsMatrix UnpackAndCombine():
00279   // Each segment of Exports will be filled by a packed row of information for each row as follows:
00280   // 1st int: GRID of row where GRID is the global row ID for the source matrix
00281   // next int:  NumEntries, Number of indices in row.
00282   // next NumEntries: The actual indices for the row.
00283 
00284   int i,j,rv;
00285   int N = TargetNumRows;
00286   int mynnz = TargetNumNonzeros;
00287   // In the case of reduced communicators, the SourceMatrix won't have the right "MyPID", so thus we have to supply it.
00288   int MyPID = MyTargetPID;
00289 
00290   int SizeofIntType = sizeof(int_type);
00291 
00292   // Zero the rowptr
00293   for(i=0; i<N+1; i++) CSR_rowptr[i]=0;
00294 
00295   // SameIDs: Always first, always in the same place
00296   for(i=0; i<NumSameIDs; i++)
00297     CSR_rowptr[i]=SourceMatrix.NumMyEntries(i);
00298   
00299   // PermuteIDs: Still local, but reordered
00300   for(i=0; i<NumPermuteIDs; i++)
00301     CSR_rowptr[PermuteToLIDs[i]] = SourceMatrix.NumMyEntries(PermuteFromLIDs[i]);
00302   
00303   // RemoteIDs:  RemoteLIDs tells us the ID, we need to look up the length the hard way.  See UnpackAndCombine for where this code came from
00304   if(NumRemoteIDs > 0) {
00305     double * dintptr = (double *) Imports;    
00306     int_type    *  intptr = (int_type *) dintptr;
00307     int   NumEntries = (int) intptr[1];
00308     int      IntSize = 1 + (((2*NumEntries+2)*SizeofIntType)/(int)sizeof(double));
00309     for(i=0; i<NumRemoteIDs; i++) {
00310       CSR_rowptr[RemoteLIDs[i]] += NumEntries;
00311       
00312       if( i < (NumRemoteIDs-1) ) {
00313   dintptr   += IntSize + NumEntries;
00314   intptr     = (int_type *) dintptr;
00315   NumEntries = (int) intptr[1];
00316   IntSize    = 1 + (((2*NumEntries+2)*SizeofIntType)/(int)sizeof(double));
00317       }
00318     }
00319   }
00320 
00321   // If multiple procs contribute to a row;
00322   std::vector<int> NewStartRow(N+1);
00323  
00324   // Turn row length into a real CSR_rowptr
00325   int last_len = CSR_rowptr[0];
00326   CSR_rowptr[0] = 0;
00327   for(i=1; i<N+1; i++){
00328     int new_len    = CSR_rowptr[i];
00329     CSR_rowptr[i]  = last_len + CSR_rowptr[i-1];
00330     NewStartRow[i] = CSR_rowptr[i];
00331     last_len       = new_len;
00332   }
00333 
00334   // Preseed TargetPids with -1 for local
00335   if(TargetPids.size()!=(size_t)mynnz) TargetPids.resize(mynnz);
00336   TargetPids.assign(mynnz,-1);
00337  
00338   // Grab pointers for SourceMatrix
00339   int    * Source_rowptr, * Source_colind;
00340   double * Source_vals;
00341   rv=SourceMatrix.ExtractCrsDataPointers(Source_rowptr,Source_colind,Source_vals);
00342   if(rv) throw std::runtime_error("UnpackAndCombineIntoCrsArrays: failed in ExtractCrsDataPointers");
00343 
00344   // SameIDs: Copy the data over
00345   for(i=0; i<NumSameIDs; i++) {
00346     int FromRow = Source_rowptr[i];
00347     int ToRow   = CSR_rowptr[i];
00348     NewStartRow[i] += Source_rowptr[i+1]-Source_rowptr[i];
00349 
00350     for(j=Source_rowptr[i]; j<Source_rowptr[i+1]; j++) {
00351       CSR_vals[ToRow + j - FromRow]   = Source_vals[j];      
00352       CSR_colind[ToRow + j - FromRow] = (int_type) SourceMatrix.GCID64(Source_colind[j]);
00353       TargetPids[ToRow + j - FromRow] = (SourcePids[Source_colind[j]] != MyPID) ? SourcePids[Source_colind[j]] : -1;
00354     }
00355   }
00356 
00357   // PermuteIDs: Copy the data over
00358   for(i=0; i<NumPermuteIDs; i++) {
00359     int FromLID = PermuteFromLIDs[i];
00360     int FromRow = Source_rowptr[FromLID];
00361     int ToRow   = CSR_rowptr[PermuteToLIDs[i]];
00362 
00363     NewStartRow[PermuteToLIDs[i]] += Source_rowptr[FromLID+1]-Source_rowptr[FromLID];
00364 
00365     for(j=Source_rowptr[FromLID]; j<Source_rowptr[FromLID+1]; j++) {
00366       CSR_vals[ToRow + j - FromRow]   = Source_vals[j];      
00367       CSR_colind[ToRow + j - FromRow] = (int_type) SourceMatrix.GCID64(Source_colind[j]);
00368       TargetPids[ToRow + j - FromRow] = (SourcePids[Source_colind[j]] != MyPID) ? SourcePids[Source_colind[j]] : -1;
00369     }
00370   }
00371 
00372   // RemoteIDs: Loop structure following UnpackAndCombine  
00373   if(NumRemoteIDs > 0) {
00374     double * dintptr   = (double *) Imports;
00375     int_type * intptr  = (int_type *) dintptr;
00376     int NumEntries     = (int) intptr[1];
00377     int IntSize        = 1 + (((2*NumEntries+2)*SizeofIntType)/(int)sizeof(double));
00378     double* valptr     = dintptr + IntSize;
00379       
00380     for (i=0; i<NumRemoteIDs; i++) {
00381       int ToLID    = RemoteLIDs[i];
00382       int StartRow = NewStartRow[ToLID];
00383       NewStartRow[ToLID]+=NumEntries; 
00384       
00385       double * values    = valptr;
00386       int_type * Indices = intptr + 2;
00387       for(j=0; j<NumEntries; j++){
00388   CSR_vals[StartRow + j]   = values[j];   
00389   CSR_colind[StartRow + j] = Indices[2*j];
00390   TargetPids[StartRow + j] = (Indices[2*j+1] != MyPID) ? Indices[2*j+1] : -1;
00391       }
00392       if( i < (NumRemoteIDs-1) ) {
00393   dintptr   += IntSize + NumEntries;
00394   intptr     = (int_type *) dintptr;
00395   NumEntries = (int) intptr[1];
00396   IntSize    = 1 + (((2*NumEntries+2)*SizeofIntType)/(int)sizeof(double));
00397   valptr     = dintptr + IntSize;
00398       }
00399     }
00400   }
00401 
00402   return 0;
00403 }
00404 
00405 
00406 
00407 // =========================================================================
00408 // =========================================================================
00409 int UnpackAndCombineIntoCrsArrays(const Epetra_CrsMatrix& SourceMatrix, 
00410           int NumSameIDs,
00411           int NumRemoteIDs,
00412           const int * RemoteLIDs,
00413           int NumPermuteIDs,
00414           const int *PermuteToLIDs,
00415           const int *PermuteFromLIDs,
00416           int LenImports,
00417           char* Imports,
00418           int TargetNumRows,
00419           int TargetNumNonzeros,
00420           int MyTargetPID,
00421           int * CSR_rowptr,
00422           int * CSR_colind,
00423           double * CSR_values,
00424           const std::vector<int> &SourcePids,
00425           std::vector<int> &TargetPids)
00426 {
00427   if(SourceMatrix.RowMap().GlobalIndicesInt()) {
00428     return TUnpackAndCombineIntoCrsArrays<int>(SourceMatrix,NumSameIDs,NumRemoteIDs,RemoteLIDs,NumPermuteIDs,
00429                  PermuteToLIDs,PermuteFromLIDs,LenImports,Imports,TargetNumRows,TargetNumNonzeros,MyTargetPID,
00430                  CSR_rowptr,CSR_colind,CSR_values,SourcePids,TargetPids);
00431   }
00432   else 
00433     throw std::runtime_error("UnpackAndCombineIntoCrsArrays: int type not matched.");
00434 }
00435 
00436 // =========================================================================
00437 // =========================================================================
00438 int UnpackAndCombineIntoCrsArrays(const Epetra_CrsMatrix& SourceMatrix, 
00439           int NumSameIDs,
00440           int NumRemoteIDs,
00441           const int * RemoteLIDs,
00442           int NumPermuteIDs,
00443           const int *PermuteToLIDs,
00444           const int *PermuteFromLIDs,
00445           int LenImports,
00446           char* Imports,
00447           int TargetNumRows,
00448           int TargetNumNonzeros,
00449           int MyTargetPID,
00450           int * CSR_rowptr,
00451           long long * CSR_colind,
00452           double * CSR_values,
00453           const std::vector<int> &SourcePids,
00454           std::vector<int> &TargetPids)
00455 {
00456   if(SourceMatrix.RowMap().GlobalIndicesLongLong()) {
00457     return TUnpackAndCombineIntoCrsArrays<long long>(SourceMatrix,NumSameIDs,NumRemoteIDs,RemoteLIDs,NumPermuteIDs,
00458                  PermuteToLIDs,PermuteFromLIDs,LenImports,Imports,TargetNumRows,TargetNumNonzeros,MyTargetPID,
00459                  CSR_rowptr,CSR_colind,CSR_values,SourcePids,TargetPids);
00460   }
00461   else
00462     throw std::runtime_error("UnpackAndCombineIntoCrsArrays: int type not matched.");
00463 }
00464 
00465 
00466 // =========================================================================
00467 // =========================================================================
00468   template<typename int_type, class MapType1, class MapType2>
00469  int TLowCommunicationMakeColMapAndReindex(int N, const int * rowptr, int * colind_LID, const int_type *colind_GID, const Epetra_Map& domainMap, const int * owningPIDs, bool SortGhostsAssociatedWithEachProcessor, std::vector<int>& RemotePIDs, MapType1 & NewColMap)
00470    {
00471   int i,j;
00472 
00473 
00474   // Sanity checks
00475   bool UseLL;
00476   if(domainMap.GlobalIndicesLongLong()) UseLL=true;
00477   else if(domainMap.GlobalIndicesInt()) UseLL=false;
00478   else throw std::runtime_error("LowCommunicationMakeColMapAndReindex: cannot detect int type.");
00479 
00480   // Scan all column indices and sort into two groups: 
00481   // Local:  those whose GID matches a GID of the domain map on this processor and
00482   // Remote: All others.
00483   int numDomainElements = domainMap.NumMyElements();
00484   bool * LocalGIDs  = 0;
00485   if (numDomainElements>0) LocalGIDs  = new bool[numDomainElements];
00486   for (i=0; i<numDomainElements; i++) LocalGIDs[i] = false; // Assume domain GIDs are not local
00487 
00488   bool DoSizes = !domainMap.ConstantElementSize(); // If not constant element size, then error
00489   if(DoSizes) throw std::runtime_error("LowCommunicationMakeColMapAndReindex: cannot handle non-constant sized domainMap.");
00490 
00491 
00492   // In principle it is good to have RemoteGIDs and RemotGIDList be as long as the number of remote GIDs
00493   // on this processor, but this would require two passes through the column IDs, so we make it the max of 100
00494   // and the number of block rows.
00495   const int numMyBlockRows = N;
00496   int  hashsize = numMyBlockRows; if (hashsize < 100) hashsize = 100;
00497   Epetra_HashTable<int_type> RemoteGIDs(hashsize); 
00498   std::vector<int_type> RemoteGIDList; RemoteGIDList.reserve(hashsize);
00499   std::vector<int> PIDList;            PIDList.reserve(hashsize);
00500 
00501   // Here we start using the *int* colind array.  If int_type==int this clobbers the GIDs, if
00502   // int_type==long long, then this is the first use of the colind array.
00503   // For *local* GID's set colind with with their LID in the domainMap.  For *remote* GIDs, 
00504   // we set colind with (numDomainElements+NumRemoteColGIDs) before the increment of
00505   // the remote count.  These numberings will be separate because no local LID is greater 
00506   // than numDomainElements. 
00507 
00508   int NumLocalColGIDs = 0;
00509   int NumRemoteColGIDs = 0;
00510   for(i = 0; i < numMyBlockRows; i++) {
00511     for(j = rowptr[i]; j < rowptr[i+1]; j++) {
00512       int_type GID = colind_GID[j];
00513       // Check if GID matches a row GID
00514       int LID = domainMap.LID(GID);
00515       if(LID != -1) {
00516   bool alreadyFound = LocalGIDs[LID];
00517   if (!alreadyFound) {
00518           LocalGIDs[LID] = true; // There is a column in the graph associated with this domain map GID
00519           NumLocalColGIDs++;
00520   }
00521   colind_LID[j] = LID; 
00522       }
00523       else {
00524   int_type hash_value=RemoteGIDs.Get(GID);
00525   if(hash_value  == -1) { // This means its a new remote GID
00526     int PID = owningPIDs[j];
00527     if(PID==-1) throw std::runtime_error("LowCommunicationMakeColMapAndReindex: Cannot figure out if PID is owned.");
00528     colind_LID[j] = numDomainElements + NumRemoteColGIDs;
00529     RemoteGIDs.Add(GID, NumRemoteColGIDs);
00530     RemoteGIDList.push_back(GID);
00531     PIDList.push_back(PID);
00532     NumRemoteColGIDs++;
00533   }
00534   else
00535     colind_LID[j] = numDomainElements + hash_value;   
00536       }
00537     }
00538   }
00539 
00540   // Possible short-circuit:  If all domain map GIDs are present as column indices, then set ColMap=domainMap and quit
00541   if (domainMap.Comm().NumProc()==1) { 
00542     
00543     if (NumRemoteColGIDs!=0) {
00544       throw std::runtime_error("Some column IDs are not in domainMap.  If matrix is rectangular, you must pass in a domainMap"); 
00545       // Sanity test: When one processor,there can be no remoteGIDs
00546     }
00547     if (NumLocalColGIDs==numDomainElements) {
00548       if (LocalGIDs!=0) delete [] LocalGIDs; 
00549       // In this case, we just use the domainMap's indices, which is, not coincidently, what we clobbered colind with up above anyway. 
00550       // No further reindexing is needed.
00551       NewColMap = domainMap;
00552       return 0;
00553     }
00554   }
00555       
00556   // Now build integer array containing column GIDs
00557   // Build back end, containing remote GIDs, first
00558   int numMyBlockCols = NumLocalColGIDs + NumRemoteColGIDs;
00559   std::vector<int_type> ColIndices;
00560   int_type * RemoteColIndices=0;
00561   if(numMyBlockCols > 0) {
00562     ColIndices.resize(numMyBlockCols);
00563     if(NumLocalColGIDs!=numMyBlockCols) RemoteColIndices = &ColIndices[NumLocalColGIDs]; // Points to back end of ColIndices
00564     else RemoteColIndices=0;
00565   }
00566 
00567   for(i = 0; i < NumRemoteColGIDs; i++) 
00568     RemoteColIndices[i] = RemoteGIDList[i]; 
00569 
00570   // Build permute array for *remote* reindexing.
00571   std::vector<int> RemotePermuteIDs(NumRemoteColGIDs);
00572   for(i=0; i<NumRemoteColGIDs; i++) RemotePermuteIDs[i]=i;
00573 
00574   // Sort External column indices so that all columns coming from a given remote processor are contiguous
00575   int NumListsInt=0;
00576   int NumListsLL =0;
00577   int * IntSortLists[2];
00578   long long * LLSortLists[2];
00579   int * RemotePermuteIDs_ptr = RemotePermuteIDs.size() ? &RemotePermuteIDs[0] : 0;
00580   if(!UseLL) {
00581     // int version
00582     IntSortLists[0] = (int*) RemoteColIndices;
00583     IntSortLists[1] = RemotePermuteIDs_ptr;
00584     NumListsInt=2;
00585   }
00586   else {
00587     //LL version
00588     LLSortLists[0]  = (long long*) RemoteColIndices;
00589     IntSortLists[0] = RemotePermuteIDs_ptr;
00590     NumListsInt     = NumListsLL = 1;
00591   }
00592 
00593   int * PIDList_ptr = PIDList.size() ? &PIDList[0] : 0;
00594   Epetra_Util::Sort(true, NumRemoteColGIDs, PIDList_ptr, 0, 0, NumListsInt, IntSortLists,NumListsLL,LLSortLists);
00595 
00596   // Stash the RemotePIDs  
00597   PIDList.resize(NumRemoteColGIDs);
00598   RemotePIDs = PIDList;
00599 
00600   if (SortGhostsAssociatedWithEachProcessor) {
00601     // Sort external column indices so that columns from a given remote processor are not only contiguous
00602     // but also in ascending order. NOTE: I don't know if the number of externals associated
00603     // with a given remote processor is known at this point ... so I count them here.
00604 
00605     // NTS: Only sort the RemoteColIndices this time...
00606     int StartCurrent, StartNext;
00607     StartCurrent = 0; StartNext = 1;
00608     while ( StartNext < NumRemoteColGIDs ) {
00609       if (PIDList[StartNext]==PIDList[StartNext-1]) StartNext++;
00610       else {
00611   IntSortLists[0] =  &RemotePermuteIDs[StartCurrent];
00612   Epetra_Util::Sort(true,StartNext-StartCurrent, &(RemoteColIndices[StartCurrent]),0,0,1,IntSortLists,0,0);
00613         StartCurrent = StartNext; StartNext++;
00614       }
00615     }
00616     IntSortLists[0] =  &RemotePermuteIDs[StartCurrent];
00617     Epetra_Util::Sort(true, StartNext-StartCurrent, &(RemoteColIndices[StartCurrent]), 0, 0, 1,IntSortLists,0,0);
00618   }
00619 
00620   // Reverse the permutation to get the information we actually care about
00621   std::vector<int> ReverseRemotePermuteIDs(NumRemoteColGIDs);
00622   for(i=0; i<NumRemoteColGIDs; i++) ReverseRemotePermuteIDs[RemotePermuteIDs[i]]=i;
00623 
00624   // Build permute array for *local* reindexing.
00625   bool use_local_permute=false;
00626   std::vector<int> LocalPermuteIDs(numDomainElements);
00627 
00628   // Now fill front end. Two cases:
00629   // (1) If the number of Local column GIDs is the same as the number of Local domain GIDs, we
00630   //     can simply read the domain GIDs into the front part of ColIndices, otherwise 
00631   // (2) We step through the GIDs of the domainMap, checking to see if each domain GID is a column GID.
00632   //     we want to do this to maintain a consistent ordering of GIDs between the columns and the domain.
00633 
00634   if(NumLocalColGIDs == domainMap.NumMyElements()) {
00635     if(NumLocalColGIDs > 0) {
00636       domainMap.MyGlobalElements(&ColIndices[0]); // Load Global Indices into first numMyBlockCols elements column GID list
00637     }
00638   }
00639   else {
00640     int_type* MyGlobalElements = 0;
00641     domainMap.MyGlobalElementsPtr(MyGlobalElements);
00642 
00643     int NumLocalAgain = 0;
00644     use_local_permute = true;    
00645     for(i = 0; i < numDomainElements; i++) {
00646       if(LocalGIDs[i]) {
00647   LocalPermuteIDs[i] = NumLocalAgain;
00648   ColIndices[NumLocalAgain++] = MyGlobalElements[i];
00649       }
00650     }
00651     assert(NumLocalAgain==NumLocalColGIDs); // Sanity test
00652   }
00653 
00654   // Done with this array
00655   if (LocalGIDs!=0) delete [] LocalGIDs; 
00656 
00657   // Make Column map with same element sizes as Domain map 
00658   int_type * ColIndices_ptr  = ColIndices.size() ? &ColIndices[0] : 0;
00659   MapType2 temp((int_type)(-1), numMyBlockCols, ColIndices_ptr, (int_type)domainMap.IndexBase64(), domainMap.Comm());
00660   NewColMap = temp;
00661 
00662   // Low-cost reindex of the matrix
00663   for(i=0; i<numMyBlockRows; i++){
00664     for(j=rowptr[i]; j<rowptr[i+1]; j++){
00665       int ID=colind_LID[j];
00666       if(ID < numDomainElements){
00667   if(use_local_permute) colind_LID[j] = LocalPermuteIDs[colind_LID[j]];
00668   // In the case where use_local_permute==false, we just copy the DomainMap's ordering, which it so happens
00669   // is what we put in colind to begin with.
00670       }
00671       else
00672   colind_LID[j] =  NumLocalColGIDs + ReverseRemotePermuteIDs[colind_LID[j]-numDomainElements];
00673     }
00674   }
00675   
00676   return 0;
00677 }
00678 
00679 
00680 // =========================================================================
00681 // =========================================================================
00682 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00683 int LowCommunicationMakeColMapAndReindex(int N, const int * rowptr, int * colind, const Epetra_Map& domainMap, const int * owningPIDs, bool SortGhostsAssociatedWithEachProcessor, std::vector<int>& RemotePIDs, Epetra_BlockMap & NewColMap) {
00684 
00685   
00686   if(domainMap.GlobalIndicesInt()) 
00687     return TLowCommunicationMakeColMapAndReindex<int,Epetra_BlockMap,Epetra_Map>(N,rowptr,colind,colind,domainMap,owningPIDs,SortGhostsAssociatedWithEachProcessor,RemotePIDs,NewColMap); 
00688   else
00689     throw std::runtime_error("LowCommunicationMakeColMapAndReindex: GID type mismatch.");
00690   return -1;
00691 }
00692 #endif
00693 
00694 // =========================================================================
00695 // =========================================================================
00696 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00697 int LowCommunicationMakeColMapAndReindex(int N, const int * rowptr, int * colind_LID, long long * colind_GID, const Epetra_Map& domainMap, const int * owningPIDs, bool SortGhostsAssociatedWithEachProcessor, std::vector<int>& RemotePIDs, Epetra_BlockMap & NewColMap) {
00698 
00699   
00700   if(domainMap.GlobalIndicesLongLong()) 
00701     return TLowCommunicationMakeColMapAndReindex<long long,Epetra_BlockMap,Epetra_Map>(N,rowptr,colind_LID,colind_GID,domainMap,owningPIDs,SortGhostsAssociatedWithEachProcessor,RemotePIDs,NewColMap); 
00702   else
00703     throw std::runtime_error("LowCommunicationMakeColMapAndReindex: GID type mismatch.");
00704   return -1;
00705 }
00706 #endif
00707 
00708 
00709 }// end namespace Epetra_Import_Util
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines