Epetra Package Browser (Single Doxygen Collection) Development
Epetra_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 #include "Epetra_ConfigDefs.h"
00044 #include "Epetra_Util.h"
00045 #include "Epetra_Object.h"
00046 #include "Epetra_Comm.h"
00047 #include "Epetra_Directory.h"
00048 #include "Epetra_Map.h"
00049 #include "Epetra_LocalMap.h"
00050 #include "Epetra_CrsGraph.h"
00051 #include "Epetra_CrsMatrix.h"
00052 #include "Epetra_MultiVector.h"
00053 #include "Epetra_IntVector.h"
00054 #include "Epetra_GIDTypeVector.h"
00055 #include "Epetra_Import.h"
00056 #include <limits>
00057 
00058 #ifdef HAVE_MPI
00059 #include "Epetra_MpiDistributor.h"
00060 #endif
00061 
00062 const double Epetra_Util::chopVal_ = 1.0e-15;
00063 
00064 //=========================================================================
00065 double Epetra_Util::Chop(const double & Value){
00066   if (std::abs(Value) < chopVal_) return 0;
00067   return Value;
00068 }
00069 
00070 //=========================================================================
00071 unsigned int Epetra_Util::RandomInt() {
00072 
00073   const int a = 16807;
00074   const int m = 2147483647;
00075   const int q = 127773;
00076   const int r = 2836;
00077 
00078   int hi = Seed_ / q;
00079   int lo = Seed_ % q;
00080   int test = a * lo - r * hi;
00081   if (test > 0)
00082     Seed_ = test;
00083   else
00084     Seed_ = test + m;
00085 
00086   return(Seed_);
00087 }
00088 
00089 //=========================================================================
00090 double Epetra_Util::RandomDouble() {
00091   const double Modulus = 2147483647.0;
00092   const double DbleOne = 1.0;
00093   const double DbleTwo = 2.0;
00094 
00095   double randdouble = RandomInt(); // implicit conversion from int to double
00096   randdouble = DbleTwo * (randdouble / Modulus) - DbleOne; // scale to (-1.0, 1.0)
00097 
00098   return(randdouble);
00099 }
00100 
00101 //=========================================================================
00102 unsigned int Epetra_Util::Seed() const {
00103   return(Seed_);
00104 }
00105 
00106 //=========================================================================
00107 int Epetra_Util::SetSeed(unsigned int Seed_in) {
00108   Seed_ = Seed_in;
00109   return(0);
00110 }
00111 
00112 //=============================================================================
00113 template<typename T>
00114 void Epetra_Util::Sort(bool SortAscending, int NumKeys, T * Keys,
00115            int NumDoubleCompanions,double ** DoubleCompanions,
00116            int NumIntCompanions, int ** IntCompanions,
00117            int NumLongLongCompanions, long long ** LongLongCompanions)
00118 {
00119   int i;
00120 
00121   int n = NumKeys;
00122   T * const list = Keys;
00123   int m = n/2;
00124 
00125   while (m > 0) {
00126     int max = n - m;
00127     for (int j=0; j<max; j++)
00128       {
00129   for (int k=j; k>=0; k-=m)
00130     {
00131       if ((SortAscending && list[k+m] >= list[k]) ||
00132     ( !SortAscending && list[k+m] <= list[k]))
00133         break;
00134       T temp = list[k+m];
00135       list[k+m] = list[k];
00136       list[k] = temp;
00137       for (i=0; i<NumDoubleCompanions; i++) {
00138         double dtemp = DoubleCompanions[i][k+m];
00139       DoubleCompanions[i][k+m] = DoubleCompanions[i][k];
00140       DoubleCompanions[i][k] = dtemp;
00141       }
00142       for (i=0; i<NumIntCompanions; i++) {
00143         int itemp = IntCompanions[i][k+m];
00144       IntCompanions[i][k+m] = IntCompanions[i][k];
00145       IntCompanions[i][k] = itemp;
00146       }
00147       for (i=0; i<NumLongLongCompanions; i++) {
00148         long long LLtemp = LongLongCompanions[i][k+m];
00149       LongLongCompanions[i][k+m] = LongLongCompanions[i][k];
00150       LongLongCompanions[i][k] = LLtemp;
00151       }
00152     }
00153       }
00154     m = m/2;
00155   }
00156 }
00157 
00158 void Epetra_Util::Sort(bool SortAscending, int NumKeys, int * Keys,
00159            int NumDoubleCompanions,double ** DoubleCompanions,
00160            int NumIntCompanions, int ** IntCompanions,
00161            int NumLongLongCompanions, long long ** LongLongCompanions)
00162 {
00163   Sort<int>(SortAscending, NumKeys, Keys,
00164            NumDoubleCompanions, DoubleCompanions,
00165            NumIntCompanions, IntCompanions,
00166            NumLongLongCompanions, LongLongCompanions);
00167 }
00168 
00169 void Epetra_Util::Sort(bool SortAscending, int NumKeys, long long * Keys,
00170            int NumDoubleCompanions,double ** DoubleCompanions,
00171            int NumIntCompanions, int ** IntCompanions,
00172            int NumLongLongCompanions, long long ** LongLongCompanions)
00173 {
00174   Sort<long long>(SortAscending, NumKeys, Keys,
00175            NumDoubleCompanions, DoubleCompanions,
00176            NumIntCompanions, IntCompanions,
00177            NumLongLongCompanions, LongLongCompanions);
00178 }
00179 
00180 void Epetra_Util::Sort(bool SortAscending, int NumKeys, double * Keys,
00181            int NumDoubleCompanions,double ** DoubleCompanions,
00182            int NumIntCompanions, int ** IntCompanions,
00183            int NumLongLongCompanions, long long ** LongLongCompanions)
00184 {
00185   Sort<double>(SortAscending, NumKeys, Keys,
00186          NumDoubleCompanions, DoubleCompanions,
00187          NumIntCompanions, IntCompanions,
00188          NumLongLongCompanions, LongLongCompanions);
00189 }
00190 
00191 
00192 void Epetra_Util::Sort(bool SortAscending, int NumKeys, int * Keys,
00193            int NumDoubleCompanions,double ** DoubleCompanions,
00194            int NumIntCompanions, int ** IntCompanions)
00195 {
00196   int i;
00197 
00198   int n = NumKeys;
00199   int * const list = Keys;
00200   int m = n/2;
00201 
00202   while (m > 0) {
00203     int max = n - m;
00204     for (int j=0; j<max; j++)
00205       {
00206   for (int k=j; k>=0; k-=m)
00207     {
00208       if ((SortAscending && list[k+m] >= list[k]) ||
00209     ( !SortAscending && list[k+m] <= list[k]))
00210         break;
00211       int temp = list[k+m];
00212       list[k+m] = list[k];
00213       list[k] = temp;
00214       for (i=0; i<NumDoubleCompanions; i++) {
00215         double dtemp = DoubleCompanions[i][k+m];
00216       DoubleCompanions[i][k+m] = DoubleCompanions[i][k];
00217       DoubleCompanions[i][k] = dtemp;
00218       }
00219       for (i=0; i<NumIntCompanions; i++) {
00220         int itemp = IntCompanions[i][k+m];
00221       IntCompanions[i][k+m] = IntCompanions[i][k];
00222       IntCompanions[i][k] = itemp;
00223       }
00224     }
00225       }
00226     m = m/2;
00227   }
00228 }
00229 
00230 //----------------------------------------------------------------------------
00231 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES // FIXME
00232 // FIXME long long
00233 Epetra_Map
00234 Epetra_Util::Create_OneToOne_Map(const Epetra_Map& usermap,
00235          bool high_rank_proc_owns_shared)
00236 {
00237   //if usermap is already 1-to-1 then we'll just return a copy of it.
00238   if (usermap.IsOneToOne()) {
00239     Epetra_Map newmap(usermap);
00240     return(newmap);
00241   }
00242 
00243   int myPID = usermap.Comm().MyPID();
00244   Epetra_Directory* directory = usermap.Comm().CreateDirectory(usermap);
00245 
00246   int numMyElems = usermap.NumMyElements();
00247   const int* myElems = usermap.MyGlobalElements();
00248 
00249   int* owner_procs = new int[numMyElems];
00250 
00251   directory->GetDirectoryEntries(usermap, numMyElems, myElems, owner_procs,
00252          0, 0, high_rank_proc_owns_shared);
00253 
00254   //we'll fill a list of map-elements which belong on this processor
00255 
00256   int* myOwnedElems = new int[numMyElems];
00257   int numMyOwnedElems = 0;
00258 
00259   for(int i=0; i<numMyElems; ++i) {
00260     int GID = myElems[i];
00261     int owner = owner_procs[i];
00262 
00263     if (myPID == owner) {
00264       myOwnedElems[numMyOwnedElems++] = GID;
00265     }
00266   }
00267 
00268   Epetra_Map one_to_one_map(-1, numMyOwnedElems, myOwnedElems,
00269        usermap.IndexBase(), usermap.Comm()); // CJ TODO FIXME long long
00270 
00271   delete [] myOwnedElems;
00272   delete [] owner_procs;
00273   delete directory;
00274 
00275   return(one_to_one_map);
00276 }
00277 #endif // EPETRA_NO_32BIT_GLOBAL_INDICES
00278 //----------------------------------------------------------------------------
00279 
00280 template<typename int_type>
00281 static Epetra_Map TCreate_Root_Map(const Epetra_Map& usermap,
00282          int root)
00283 {
00284   int numProc = usermap.Comm().NumProc();
00285   if (numProc==1) {
00286     Epetra_Map newmap(usermap);
00287     return(newmap);
00288   }
00289 
00290   const Epetra_Comm & comm = usermap.Comm();
00291   bool isRoot = usermap.Comm().MyPID()==root;
00292 
00293   //if usermap is already completely owned by root then we'll just return a copy of it.
00294   int quickreturn = 0;
00295   int globalquickreturn = 0;
00296 
00297   if (isRoot) {
00298     if (usermap.NumMyElements()==usermap.NumGlobalElements64()) quickreturn = 1;
00299   }
00300   else {
00301     if (usermap.NumMyElements()==0) quickreturn = 1;
00302   }
00303   usermap.Comm().MinAll(&quickreturn, &globalquickreturn, 1);
00304 
00305   if (globalquickreturn==1) {
00306     Epetra_Map newmap(usermap);
00307     return(newmap);
00308   }
00309 
00310   // Linear map: Simple case, just put all GIDs linearly on root processor
00311   if (usermap.LinearMap() && root!=-1) {
00312     int numMyElements = 0;
00313     if(usermap.MaxAllGID64()+1 > std::numeric_limits<int>::max())
00314       throw "Epetra_Util::Create_Root_Map: cannot fit all gids in int";
00315     if (isRoot) numMyElements = (int)(usermap.MaxAllGID64()+1);
00316     Epetra_Map newmap((int_type) -1, numMyElements, (int_type)usermap.IndexBase64(), comm);
00317     return(newmap);
00318   }
00319 
00320   if (!usermap.UniqueGIDs())
00321     throw usermap.ReportError("usermap must have unique GIDs",-1);
00322 
00323   // General map
00324 
00325   // Build IntVector of the GIDs, then ship them to root processor
00326   int numMyElements = usermap.NumMyElements();
00327   Epetra_Map allGidsMap((int_type) -1, numMyElements, (int_type) 0, comm);
00328   typename Epetra_GIDTypeVector<int_type>::impl allGids(allGidsMap);
00329   for (int i=0; i<numMyElements; i++) allGids[i] = (int_type) usermap.GID64(i);
00330 
00331   if(usermap.MaxAllGID64() > std::numeric_limits<int>::max())
00332     throw "Epetra_Util::Create_Root_Map: cannot fit all gids in int";
00333   int numGlobalElements = (int) usermap.NumGlobalElements64();
00334   if (root!=-1) {
00335     int n1 = 0; if (isRoot) n1 = numGlobalElements;
00336     Epetra_Map allGidsOnRootMap((int_type) -1, n1, (int_type) 0, comm);
00337     Epetra_Import importer(allGidsOnRootMap, allGidsMap);
00338     typename Epetra_GIDTypeVector<int_type>::impl allGidsOnRoot(allGidsOnRootMap);
00339     allGidsOnRoot.Import(allGids, importer, Insert);
00340 
00341     Epetra_Map rootMap((int_type)-1, allGidsOnRoot.MyLength(), allGidsOnRoot.Values(), (int_type)usermap.IndexBase64(), comm);
00342     return(rootMap);
00343   }
00344   else {
00345     int n1 = numGlobalElements;
00346     Epetra_LocalMap allGidsOnRootMap((int_type) n1, (int_type) 0, comm);
00347     Epetra_Import importer(allGidsOnRootMap, allGidsMap);
00348     typename Epetra_GIDTypeVector<int_type>::impl allGidsOnRoot(allGidsOnRootMap);
00349     allGidsOnRoot.Import(allGids, importer, Insert);
00350 
00351     Epetra_Map rootMap((int_type) -1, allGidsOnRoot.MyLength(), allGidsOnRoot.Values(), (int_type)usermap.IndexBase64(), comm);
00352 
00353     return(rootMap);
00354   }
00355 }
00356 
00357 Epetra_Map Epetra_Util::Create_Root_Map(const Epetra_Map& usermap,
00358          int root)
00359 {
00360 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00361   if(usermap.GlobalIndicesInt()) {
00362     return TCreate_Root_Map<int>(usermap, root);
00363   }
00364   else
00365 #endif
00366 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00367   if(usermap.GlobalIndicesLongLong()) {
00368     return TCreate_Root_Map<long long>(usermap, root);
00369   }
00370   else
00371 #endif
00372     throw "Epetra_Util::Create_Root_Map: GlobalIndices type unknown";
00373 }
00374 
00375 //----------------------------------------------------------------------------
00376 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES // FIXME
00377 // FIXME long long
00378 Epetra_BlockMap
00379 Epetra_Util::Create_OneToOne_BlockMap(const Epetra_BlockMap& usermap,
00380               bool high_rank_proc_owns_shared)
00381 {
00382 // FIXME long long
00383 
00384   //if usermap is already 1-to-1 then we'll just return a copy of it.
00385   if (usermap.IsOneToOne()) {
00386     Epetra_BlockMap newmap(usermap);
00387     return(newmap);
00388   }
00389 
00390   int myPID = usermap.Comm().MyPID();
00391   Epetra_Directory* directory = usermap.Comm().CreateDirectory(usermap);
00392 
00393   int numMyElems = usermap.NumMyElements();
00394   const int* myElems = usermap.MyGlobalElements();
00395 
00396   int* owner_procs = new int[numMyElems*2];
00397   int* sizes = owner_procs+numMyElems;
00398 
00399   directory->GetDirectoryEntries(usermap, numMyElems, myElems, owner_procs,
00400          0, sizes, high_rank_proc_owns_shared);
00401 
00402   //we'll fill a list of map-elements which belong on this processor
00403 
00404   int* myOwnedElems = new int[numMyElems*2];
00405   int* ownedSizes = myOwnedElems+numMyElems;
00406   int numMyOwnedElems = 0;
00407 
00408   for(int i=0; i<numMyElems; ++i) {
00409     int GID = myElems[i];
00410     int owner = owner_procs[i];
00411 
00412     if (myPID == owner) {
00413       ownedSizes[numMyOwnedElems] = sizes[i];
00414       myOwnedElems[numMyOwnedElems++] = GID;
00415     }
00416   }
00417 
00418   Epetra_BlockMap one_to_one_map(-1, numMyOwnedElems, myOwnedElems,
00419          sizes, usermap.IndexBase(), usermap.Comm()); // CJ TODO FIXME long long
00420 
00421   delete [] myOwnedElems;
00422   delete [] owner_procs;
00423   delete directory;
00424 
00425   return(one_to_one_map);
00426 }
00427 #endif // EPETRA_NO_32BIT_GLOBAL_INDICES
00428 
00429 
00430 //----------------------------------------------------------------------------
00431 int Epetra_Util::SortCrsEntries(int NumRows, const int *CRS_rowptr, int *CRS_colind, double *CRS_vals){
00432   // For each row, sort column entries from smallest to largest.
00433   // Use shell sort. Stable sort so it is fast if indices are already sorted.
00434   // Code copied from  Epetra_CrsMatrix::SortEntries()
00435   int nnz = CRS_rowptr[NumRows];
00436 
00437   for(int i = 0; i < NumRows; i++){
00438     int start=CRS_rowptr[i];
00439     if(start >= nnz) continue;
00440 
00441     double* locValues = &CRS_vals[start];
00442     int NumEntries    = CRS_rowptr[i+1] - start;
00443     int* locIndices   = &CRS_colind[start];
00444     
00445     int n = NumEntries;
00446     int m = n/2;
00447 
00448     while(m > 0) {
00449       int max = n - m;
00450       for(int j = 0; j < max; j++) {
00451   for(int k = j; k >= 0; k-=m) {
00452     if(locIndices[k+m] >= locIndices[k])
00453       break;
00454     double dtemp = locValues[k+m];
00455     locValues[k+m] = locValues[k];
00456     locValues[k] = dtemp;
00457     int itemp = locIndices[k+m];
00458     locIndices[k+m] = locIndices[k];
00459     locIndices[k] = itemp;
00460   }
00461       }
00462       m = m/2;
00463     }
00464   }
00465   return(0);
00466 }
00467 
00468 //----------------------------------------------------------------------------
00469 int Epetra_Util::SortAndMergeCrsEntries(int NumRows, int *CRS_rowptr, int *CRS_colind, double *CRS_vals){
00470   // For each row, sort column entries from smallest to largest, merging column ids that are identical by adding values.
00471   // Use shell sort. Stable sort so it is fast if indices are already sorted.
00472   // Code copied from  Epetra_CrsMatrix::SortEntries()
00473 
00474   int nnz = CRS_rowptr[NumRows];
00475   int new_curr=CRS_rowptr[0], old_curr=CRS_rowptr[0];
00476 
00477   for(int i = 0; i < NumRows; i++){
00478     int start=CRS_rowptr[i];
00479     if(start >= nnz) continue;
00480 
00481     double* locValues = &CRS_vals[start];
00482     int NumEntries    = CRS_rowptr[i+1] - start;
00483     int* locIndices   = &CRS_colind[start];
00484 
00485     // Sort phase
00486     int n = NumEntries;
00487     int m = n/2;
00488 
00489     while(m > 0) {
00490       int max = n - m;
00491       for(int j = 0; j < max; j++) {
00492         for(int k = j; k >= 0; k-=m) {
00493           if(locIndices[k+m] >= locIndices[k])
00494             break;
00495           double dtemp = locValues[k+m];
00496           locValues[k+m] = locValues[k];
00497           locValues[k] = dtemp;
00498           int itemp = locIndices[k+m];
00499           locIndices[k+m] = locIndices[k];
00500           locIndices[k] = itemp;
00501         }
00502       }
00503       m = m/2;
00504     }
00505 
00506     // Merge & shrink
00507     for(int j=CRS_rowptr[i]; j < CRS_rowptr[i+1]; j++) {
00508       if(j > CRS_rowptr[i] && CRS_colind[j]==CRS_colind[new_curr-1]) {
00509   CRS_vals[new_curr-1] += CRS_vals[j];
00510       }
00511       else if(new_curr==j) {
00512   new_curr++;
00513       }
00514       else {
00515   CRS_colind[new_curr] = CRS_colind[j];
00516   CRS_vals[new_curr]   = CRS_vals[j];
00517   new_curr++;
00518       }
00519     }
00520 
00521     CRS_rowptr[i] = old_curr;
00522     old_curr=new_curr;
00523   }
00524 
00525   CRS_rowptr[NumRows] = new_curr;
00526   return (0);
00527 }
00528 
00529 
00530 //----------------------------------------------------------------------------
00531 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00532 int Epetra_Util::GetPidGidPairs(const Epetra_Import & Importer,std::vector< std::pair<int,int> > & gpids, bool use_minus_one_for_local){
00533   // Put the (PID,GID) pair in member of Importer.TargetMap() in gpids.  If use_minus_one_for_local==true, put in -1 instead of MyPID.
00534   // This only works if we have an MpiDistributor in our Importer.  Otherwise return an error.
00535 #ifdef HAVE_MPI
00536   Epetra_MpiDistributor *D=dynamic_cast<Epetra_MpiDistributor*>(&Importer.Distributor());
00537   if(!D) EPETRA_CHK_ERR(-2);
00538 
00539   int i,j,k;
00540   int mypid=Importer.TargetMap().Comm().MyPID();
00541   int N=Importer.TargetMap().NumMyElements();
00542 
00543   // Get the importer's data
00544   const int *RemoteLIDs  = Importer.RemoteLIDs();
00545 
00546   // Get the distributor's data
00547   int NumReceives        = D->NumReceives();
00548   const int *ProcsFrom   = D->ProcsFrom();
00549   const int *LengthsFrom = D->LengthsFrom();
00550 
00551   // Resize the outgoing data structure
00552   gpids.resize(N);
00553 
00554   // Start by claiming that I own all the data
00555   if(use_minus_one_for_local)
00556     for(i=0;i <N; i++) gpids[i]=std::make_pair(-1,Importer.TargetMap().GID(i));
00557   else
00558     for(i=0;i <N; i++) gpids[i]=std::make_pair(mypid,Importer.TargetMap().GID(i));
00559 
00560   // Now, for each remote ID, record who actually owns it.  This loop follows the operation order in the
00561   // MpiDistributor so it ought to duplicate that effect.
00562   for(i=0,j=0;i<NumReceives;i++){
00563     int pid=ProcsFrom[i];
00564     for(k=0;k<LengthsFrom[i];k++){
00565       if(pid!=mypid) gpids[RemoteLIDs[j]].first=pid;
00566       j++;
00567     }
00568   }
00569   return 0;
00570 #else
00571   EPETRA_CHK_ERR(-10);
00572 #endif
00573 }
00574 #endif
00575 
00576 //----------------------------------------------------------------------------
00577 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00578 int Epetra_Util::GetPidGidPairs(const Epetra_Import & Importer,std::vector< std::pair<int,long long> > & gpids, bool use_minus_one_for_local){
00579   // Put the (PID,GID) pair in member of Importer.TargetMap() in gpids.  If use_minus_one_for_local==true, put in -1 instead of MyPID.
00580   // This only works if we have an MpiDistributor in our Importer.  Otheriwise return an error.
00581 #ifdef HAVE_MPI
00582   Epetra_MpiDistributor *D=dynamic_cast<Epetra_MpiDistributor*>(&Importer.Distributor());
00583   if(!D) EPETRA_CHK_ERR(-2);
00584 
00585   int i,j,k;
00586   int mypid=Importer.TargetMap().Comm().MyPID();
00587   int N=Importer.TargetMap().NumMyElements();
00588 
00589   // Get the importer's data
00590   const int *RemoteLIDs  = Importer.RemoteLIDs();
00591 
00592   // Get the distributor's data
00593   int NumReceives        = D->NumReceives();
00594   const int *ProcsFrom   = D->ProcsFrom();
00595   const int *LengthsFrom = D->LengthsFrom();
00596 
00597   // Resize the outgoing data structure
00598   gpids.resize(N);
00599 
00600   // Start by claiming that I own all the data
00601   if(use_minus_one_for_local)
00602     for(i=0;i <N; i++) gpids[i]=std::make_pair(-1,Importer.TargetMap().GID64(i));
00603   else
00604     for(i=0;i <N; i++) gpids[i]=std::make_pair(mypid,Importer.TargetMap().GID64(i));
00605 
00606   // Now, for each remote ID, record who actually owns it.  This loop follows the operation order in the
00607   // MpiDistributor so it ought to duplicate that effect.
00608   for(i=0,j=0;i<NumReceives;i++){
00609     int pid=ProcsFrom[i];
00610     for(k=0;k<LengthsFrom[i];k++){
00611       if(pid!=mypid) gpids[RemoteLIDs[j]].first=pid;
00612       j++;
00613     }
00614   }
00615   return 0;
00616 #else
00617   EPETRA_CHK_ERR(-10);
00618 #endif
00619 }
00620 #endif
00621 
00622 
00623 //----------------------------------------------------------------------------
00624 int Epetra_Util::GetPids(const Epetra_Import & Importer, std::vector<int> &pids, bool use_minus_one_for_local){
00625 #ifdef HAVE_MPI
00626   Epetra_MpiDistributor *D=dynamic_cast<Epetra_MpiDistributor*>(&Importer.Distributor());
00627   if(!D) EPETRA_CHK_ERR(-2);
00628 
00629   int i,j,k;
00630   int mypid=Importer.TargetMap().Comm().MyPID();
00631   int N=Importer.TargetMap().NumMyElements();
00632 
00633   // Get the importer's data
00634   const int *RemoteLIDs  = Importer.RemoteLIDs();
00635 
00636   // Get the distributor's data
00637   int NumReceives        = D->NumReceives();
00638   const int *ProcsFrom   = D->ProcsFrom();
00639   const int *LengthsFrom = D->LengthsFrom();
00640 
00641   // Resize the outgoing data structure
00642   pids.resize(N);
00643 
00644   // Start by claiming that I own all the data
00645   if(use_minus_one_for_local)
00646     for(i=0; i<N; i++) pids[i]=-1;
00647   else
00648     for(i=0; i<N; i++) pids[i]=mypid;
00649 
00650   // Now, for each remote ID, record who actually owns it.  This loop follows the operation order in the
00651   // MpiDistributor so it ought to duplicate that effect.
00652   for(i=0,j=0;i<NumReceives;i++){
00653     int pid=ProcsFrom[i];
00654     for(k=0;k<LengthsFrom[i];k++){
00655       if(pid!=mypid) pids[RemoteLIDs[j]]=pid;
00656       j++;
00657     }
00658   }
00659   return 0;
00660 #else
00661   EPETRA_CHK_ERR(-10);
00662 #endif
00663 }
00664 
00665 //----------------------------------------------------------------------------
00666 int Epetra_Util::GetRemotePIDs(const Epetra_Import & Importer, std::vector<int> &RemotePIDs){
00667 #ifdef HAVE_MPI
00668   Epetra_MpiDistributor *D=dynamic_cast<Epetra_MpiDistributor*>(&Importer.Distributor());
00669   if(!D) {
00670     RemotePIDs.resize(0);
00671     return 0;
00672   }
00673 
00674   int i,j,k;
00675 
00676   // Get the distributor's data
00677   int NumReceives        = D->NumReceives();
00678   const int *ProcsFrom   = D->ProcsFrom();
00679   const int *LengthsFrom = D->LengthsFrom();
00680 
00681   // Resize the outgoing data structure
00682   RemotePIDs.resize(Importer.NumRemoteIDs());
00683 
00684   // Now, for each remote ID, record who actually owns it.  This loop follows the operation order in the
00685   // MpiDistributor so it ought to duplicate that effect.
00686   for(i=0,j=0;i<NumReceives;i++){
00687     int pid=ProcsFrom[i];
00688     for(k=0;k<LengthsFrom[i];k++){
00689       RemotePIDs[j]=pid;
00690       j++;
00691     }
00692   }
00693   return 0;
00694 #else
00695   RemotePIDs.resize(0);
00696   return 0;
00697 #endif
00698 }
00699 
00700 //----------------------------------------------------------------------------
00701 template<typename T>
00702 int Epetra_Util_binary_search(T item,
00703                               const T* list,
00704                               int len,
00705                               int& insertPoint)
00706 {
00707   if (len < 1) {
00708     insertPoint = 0;
00709     return(-1);
00710   }
00711 
00712   unsigned start = 0, end = len - 1;
00713 
00714   while(end - start > 1) {
00715     unsigned mid = (start + end) >> 1;
00716     if (list[mid] < item) start = mid;
00717     else end = mid;
00718   }
00719 
00720   if (list[start] == item) return(start);
00721   if (list[end] == item) return(end);
00722 
00723   if (list[end] < item) {
00724     insertPoint = end+1;
00725     return(-1);
00726   }
00727 
00728   if (list[start] < item) insertPoint = end;
00729   else insertPoint = start;
00730 
00731   return(-1);
00732 }
00733 
00734 int Epetra_Util_binary_search(int item,
00735                               const int* list,
00736                               int len,
00737                               int& insertPoint)
00738 {
00739   return Epetra_Util_binary_search<int>(item, list, len, insertPoint);
00740 }
00741 
00742 //----------------------------------------------------------------------------
00743 int Epetra_Util_binary_search(long long item,
00744                               const long long* list,
00745                               int len,
00746                               int& insertPoint)
00747 {
00748   return Epetra_Util_binary_search<long long>(item, list, len, insertPoint);
00749 }
00750 
00751 //----------------------------------------------------------------------------
00752 template<typename T>
00753 int Epetra_Util_binary_search_aux(T item,
00754                               const int* list,
00755                               const T* aux_list,
00756                               int len,
00757                               int& insertPoint)
00758 {
00759   if (len < 1) {
00760     insertPoint = 0;
00761     return(-1);
00762   }
00763 
00764   unsigned start = 0, end = len - 1;
00765 
00766   while(end - start > 1) {
00767     unsigned mid = (start + end) >> 1;
00768     if (aux_list[list[mid]] < item) start = mid;
00769     else end = mid;
00770   }
00771 
00772   if (aux_list[list[start]] == item) return(start);
00773   if (aux_list[list[end]] == item) return(end);
00774 
00775   if (aux_list[list[end]] < item) {
00776     insertPoint = end+1;
00777     return(-1);
00778   }
00779 
00780   if (aux_list[list[start]] < item) insertPoint = end;
00781   else insertPoint = start;
00782 
00783   return(-1);
00784 }
00785 
00786 //----------------------------------------------------------------------------
00787 int Epetra_Util_binary_search_aux(int item,
00788                               const int* list,
00789                               const int* aux_list,
00790                               int len,
00791                               int& insertPoint)
00792 {
00793   return Epetra_Util_binary_search_aux<int>(item, list, aux_list, len, insertPoint);
00794 }
00795 
00796 //----------------------------------------------------------------------------
00797 int Epetra_Util_binary_search_aux(long long item,
00798                               const int* list,
00799                               const long long* aux_list,
00800                               int len,
00801                               int& insertPoint)
00802 {
00803   return Epetra_Util_binary_search_aux<long long>(item, list, aux_list, len, insertPoint);
00804 }
00805 
00806 
00807 
00808 //=========================================================================
00809 int Epetra_Util_ExtractHbData(Epetra_CrsMatrix * A, Epetra_MultiVector * LHS,
00810             Epetra_MultiVector * RHS,
00811             int & M, int & N, int & nz, int * & ptr,
00812             int * & ind, double * & val, int & Nrhs,
00813             double * & rhs, int & ldrhs,
00814             double * & lhs, int & ldlhs) {
00815 
00816   int ierr = 0;
00817   if (A==0) EPETRA_CHK_ERR(-1); // This matrix is defined
00818   if (!A->IndicesAreContiguous()) { // Data must be contiguous for this to work
00819     EPETRA_CHK_ERR(A->MakeDataContiguous()); // Call MakeDataContiguous() method on the matrix
00820     ierr = 1; // Warn User that we changed the matrix
00821   }
00822 
00823   M = A->NumMyRows();
00824   N = A->NumMyCols();
00825   nz = A->NumMyNonzeros();
00826   val = (*A)[0];        // Dangerous, but cheap and effective way to access first element in
00827 
00828   const Epetra_CrsGraph & Graph = A->Graph();
00829   ind = Graph[0];  // list of values and indices
00830 
00831   Nrhs = 0; // Assume no rhs, lhs
00832 
00833   if (RHS!=0) {
00834     Nrhs = RHS->NumVectors();
00835     if (Nrhs>1)
00836     if (!RHS->ConstantStride()) {EPETRA_CHK_ERR(-2)}; // Must have strided vectors
00837     ldrhs = RHS->Stride();
00838     rhs = (*RHS)[0]; // Dangerous but effective (again)
00839   }
00840   if (LHS!=0) {
00841     int Nlhs = LHS->NumVectors();
00842     if (Nlhs!=Nrhs) {EPETRA_CHK_ERR(-3)}; // Must have same number of rhs and lhs
00843     if (Nlhs>1)
00844     if (!LHS->ConstantStride()) {EPETRA_CHK_ERR(-4)}; // Must have strided vectors
00845   ldlhs = LHS->Stride();
00846   lhs = (*LHS)[0];
00847   }
00848 
00849   // Finally build ptr vector
00850 
00851   if (ptr==0) {
00852     ptr = new int[M+1];
00853     ptr[0] = 0;
00854     for (int i=0; i<M; i++) ptr[i+1] = ptr[i] + Graph.NumMyIndices(i);
00855   }
00856   EPETRA_CHK_ERR(ierr);
00857   return(0);
00858 }
00859 
00860 // Explicitly instantiate these two even though a call to them is made.
00861 // Possible fix for a bug reported by Kenneth Belcourt.
00862 template void Epetra_Util::Sort<int>      (bool, int, int *,       int, double **, int, int **, int, long long **);
00863 template void Epetra_Util::Sort<long long>(bool, int, long long *, int, double **, int, int **, int, long long **);
00864 
00865 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines