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, int * Keys, 
00181            int NumDoubleCompanions,double ** DoubleCompanions, 
00182            int NumIntCompanions, int ** IntCompanions)
00183 {
00184   int i;
00185 
00186   int n = NumKeys;
00187   int * const list = Keys;
00188   int m = n/2;
00189   
00190   while (m > 0) {
00191     int max = n - m;
00192     for (int j=0; j<max; j++)
00193       {
00194   for (int k=j; k>=0; k-=m)
00195     {
00196       if ((SortAscending && list[k+m] >= list[k]) || 
00197     ( !SortAscending && list[k+m] <= list[k]))
00198         break;
00199       int temp = list[k+m];
00200       list[k+m] = list[k];
00201       list[k] = temp;
00202       for (i=0; i<NumDoubleCompanions; i++) {
00203         double dtemp = DoubleCompanions[i][k+m];
00204       DoubleCompanions[i][k+m] = DoubleCompanions[i][k];
00205       DoubleCompanions[i][k] = dtemp;
00206       }
00207       for (i=0; i<NumIntCompanions; i++) {
00208         int itemp = IntCompanions[i][k+m];
00209       IntCompanions[i][k+m] = IntCompanions[i][k];
00210       IntCompanions[i][k] = itemp;
00211       }
00212     }
00213       }
00214     m = m/2;
00215   }
00216 }
00217 
00218 //----------------------------------------------------------------------------
00219 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES // FIXME
00220 // FIXME long long
00221 Epetra_Map
00222 Epetra_Util::Create_OneToOne_Map(const Epetra_Map& usermap,
00223          bool high_rank_proc_owns_shared)
00224 {
00225   //if usermap is already 1-to-1 then we'll just return a copy of it.
00226   if (usermap.IsOneToOne()) {
00227     Epetra_Map newmap(usermap);
00228     return(newmap);
00229   }
00230 
00231   int myPID = usermap.Comm().MyPID();
00232   Epetra_Directory* directory = usermap.Comm().CreateDirectory(usermap);
00233 
00234   int numMyElems = usermap.NumMyElements();
00235   const int* myElems = usermap.MyGlobalElements();
00236 
00237   int* owner_procs = new int[numMyElems];
00238 
00239   directory->GetDirectoryEntries(usermap, numMyElems, myElems, owner_procs,
00240          0, 0, high_rank_proc_owns_shared);
00241 
00242   //we'll fill a list of map-elements which belong on this processor
00243 
00244   int* myOwnedElems = new int[numMyElems];
00245   int numMyOwnedElems = 0;
00246 
00247   for(int i=0; i<numMyElems; ++i) {
00248     int GID = myElems[i];
00249     int owner = owner_procs[i];
00250 
00251     if (myPID == owner) {
00252       myOwnedElems[numMyOwnedElems++] = GID;
00253     }
00254   }
00255 
00256   Epetra_Map one_to_one_map(-1, numMyOwnedElems, myOwnedElems,
00257        usermap.IndexBase(), usermap.Comm()); // CJ TODO FIXME long long
00258 
00259   delete [] myOwnedElems;
00260   delete [] owner_procs;
00261   delete directory;
00262 
00263   return(one_to_one_map);
00264 }
00265 #endif // EPETRA_NO_32BIT_GLOBAL_INDICES
00266 //----------------------------------------------------------------------------
00267 
00268 template<typename int_type>
00269 static Epetra_Map TCreate_Root_Map(const Epetra_Map& usermap,
00270          int root)
00271 {
00272   int numProc = usermap.Comm().NumProc();
00273   if (numProc==1) {
00274     Epetra_Map newmap(usermap);
00275     return(newmap);
00276   }
00277 
00278   const Epetra_Comm & comm = usermap.Comm();
00279   bool isRoot = usermap.Comm().MyPID()==root;
00280 
00281   //if usermap is already completely owned by root then we'll just return a copy of it.
00282   int quickreturn = 0;
00283   int globalquickreturn = 0;
00284 
00285   if (isRoot) {
00286     if (usermap.NumMyElements()==usermap.NumGlobalElements64()) quickreturn = 1;
00287   }
00288   else {
00289     if (usermap.NumMyElements()==0) quickreturn = 1;
00290   }
00291   usermap.Comm().MinAll(&quickreturn, &globalquickreturn, 1);
00292   
00293   if (globalquickreturn==1) {
00294     Epetra_Map newmap(usermap);
00295     return(newmap);
00296   }
00297   
00298   // Linear map: Simple case, just put all GIDs linearly on root processor
00299   if (usermap.LinearMap() && root!=-1) {
00300     int numMyElements = 0;
00301     if(usermap.MaxAllGID64()+1 > std::numeric_limits<int>::max())
00302       throw "Epetra_Util::Create_Root_Map: cannot fit all gids in int";
00303     if (isRoot) numMyElements = (int)(usermap.MaxAllGID64()+1);
00304     Epetra_Map newmap((int_type) -1, numMyElements, (int_type)usermap.IndexBase64(), comm);
00305     return(newmap);
00306   }
00307 
00308   if (!usermap.UniqueGIDs()) 
00309     throw usermap.ReportError("usermap must have unique GIDs",-1);
00310 
00311   // General map
00312 
00313   // Build IntVector of the GIDs, then ship them to root processor
00314   int numMyElements = usermap.NumMyElements();
00315   Epetra_Map allGidsMap((int_type) -1, numMyElements, (int_type) 0, comm);
00316   typename Epetra_GIDTypeVector<int_type>::impl allGids(allGidsMap);
00317   for (int i=0; i<numMyElements; i++) allGids[i] = (int_type) usermap.GID64(i);
00318   
00319   if(usermap.MaxAllGID64() > std::numeric_limits<int>::max())
00320     throw "Epetra_Util::Create_Root_Map: cannot fit all gids in int";
00321   int numGlobalElements = (int) usermap.NumGlobalElements64();
00322   if (root!=-1) {
00323     int n1 = 0; if (isRoot) n1 = numGlobalElements;
00324     Epetra_Map allGidsOnRootMap((int_type) -1, n1, (int_type) 0, comm);
00325     Epetra_Import importer(allGidsOnRootMap, allGidsMap);
00326     typename Epetra_GIDTypeVector<int_type>::impl allGidsOnRoot(allGidsOnRootMap);
00327     allGidsOnRoot.Import(allGids, importer, Insert);
00328     
00329     Epetra_Map rootMap((int_type)-1, allGidsOnRoot.MyLength(), allGidsOnRoot.Values(), (int_type)usermap.IndexBase64(), comm);
00330     return(rootMap);
00331   }
00332   else {
00333     int n1 = numGlobalElements;
00334     Epetra_LocalMap allGidsOnRootMap((int_type) n1, (int_type) 0, comm);
00335     Epetra_Import importer(allGidsOnRootMap, allGidsMap);
00336     typename Epetra_GIDTypeVector<int_type>::impl allGidsOnRoot(allGidsOnRootMap);
00337     allGidsOnRoot.Import(allGids, importer, Insert);
00338     
00339     Epetra_Map rootMap((int_type) -1, allGidsOnRoot.MyLength(), allGidsOnRoot.Values(), (int_type)usermap.IndexBase64(), comm);
00340 
00341     return(rootMap);
00342   }
00343 }
00344 
00345 Epetra_Map Epetra_Util::Create_Root_Map(const Epetra_Map& usermap,
00346          int root)
00347 {
00348 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00349   if(usermap.GlobalIndicesInt()) {
00350     return TCreate_Root_Map<int>(usermap, root);
00351   }
00352   else
00353 #endif
00354 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00355   if(usermap.GlobalIndicesLongLong()) {
00356     return TCreate_Root_Map<long long>(usermap, root);
00357   }
00358   else
00359 #endif
00360     throw "Epetra_Util::Create_Root_Map: GlobalIndices type unknown";
00361 }
00362 
00363 //----------------------------------------------------------------------------
00364 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES // FIXME
00365 // FIXME long long
00366 Epetra_BlockMap
00367 Epetra_Util::Create_OneToOne_BlockMap(const Epetra_BlockMap& usermap,
00368               bool high_rank_proc_owns_shared)
00369 {
00370 // FIXME long long
00371 
00372   //if usermap is already 1-to-1 then we'll just return a copy of it.
00373   if (usermap.IsOneToOne()) {
00374     Epetra_BlockMap newmap(usermap);
00375     return(newmap);
00376   }
00377 
00378   int myPID = usermap.Comm().MyPID();
00379   Epetra_Directory* directory = usermap.Comm().CreateDirectory(usermap);
00380 
00381   int numMyElems = usermap.NumMyElements();
00382   const int* myElems = usermap.MyGlobalElements();
00383 
00384   int* owner_procs = new int[numMyElems*2];
00385   int* sizes = owner_procs+numMyElems;
00386 
00387   directory->GetDirectoryEntries(usermap, numMyElems, myElems, owner_procs,
00388          0, sizes, high_rank_proc_owns_shared);
00389 
00390   //we'll fill a list of map-elements which belong on this processor
00391 
00392   int* myOwnedElems = new int[numMyElems*2];
00393   int* ownedSizes = myOwnedElems+numMyElems;
00394   int numMyOwnedElems = 0;
00395 
00396   for(int i=0; i<numMyElems; ++i) {
00397     int GID = myElems[i];
00398     int owner = owner_procs[i];
00399 
00400     if (myPID == owner) {
00401       ownedSizes[numMyOwnedElems] = sizes[i];
00402       myOwnedElems[numMyOwnedElems++] = GID;
00403     }
00404   }
00405 
00406   Epetra_BlockMap one_to_one_map(-1, numMyOwnedElems, myOwnedElems,
00407          sizes, usermap.IndexBase(), usermap.Comm()); // CJ TODO FIXME long long
00408 
00409   delete [] myOwnedElems;
00410   delete [] owner_procs;
00411   delete directory;
00412 
00413   return(one_to_one_map);
00414 }
00415 #endif // EPETRA_NO_32BIT_GLOBAL_INDICES
00416 
00417 
00418 //----------------------------------------------------------------------------
00419 int Epetra_Util::SortCrsEntries(int NumRows, const int *CRS_rowptr, int *CRS_colind, double *CRS_vals){
00420   // For each row, sort column entries from smallest to largest.
00421   // Use shell sort. Stable sort so it is fast if indices are already sorted.
00422   // Code copied from  Epetra_CrsMatrix::SortEntries() 
00423   int nnz = CRS_rowptr[NumRows];
00424 
00425   for(int i = 0; i < NumRows; i++){
00426     int start=CRS_rowptr[i];
00427     if(start >= nnz) continue;
00428 
00429     double* locValues = &CRS_vals[start];
00430     int NumEntries    = CRS_rowptr[i+1] - start;
00431     int* locIndices   = &CRS_colind[start];
00432     
00433     int n = NumEntries;
00434     int m = n/2;
00435     
00436     while(m > 0) {
00437       int max = n - m;
00438       for(int j = 0; j < max; j++) {
00439   for(int k = j; k >= 0; k-=m) {
00440     if(locIndices[k+m] >= locIndices[k])
00441       break;
00442     double dtemp = locValues[k+m];
00443     locValues[k+m] = locValues[k];
00444     locValues[k] = dtemp;
00445     int itemp = locIndices[k+m];
00446     locIndices[k+m] = locIndices[k];
00447     locIndices[k] = itemp;
00448   }
00449       }
00450       m = m/2;
00451     }
00452   }
00453   return(0);
00454 }
00455 
00456 
00457 //----------------------------------------------------------------------------
00458 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00459 int Epetra_Util::GetPidGidPairs(const Epetra_Import & Importer,std::vector< std::pair<int,int> > & gpids, bool use_minus_one_for_local){
00460   // 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.
00461   // This only works if we have an MpiDistributor in our Importer.  Otherwise return an error.
00462 #ifdef HAVE_MPI
00463   Epetra_MpiDistributor *D=dynamic_cast<Epetra_MpiDistributor*>(&Importer.Distributor());
00464   if(!D) EPETRA_CHK_ERR(-2);
00465 
00466   int i,j,k;
00467   int mypid=Importer.TargetMap().Comm().MyPID();
00468   int N=Importer.TargetMap().NumMyElements();
00469 
00470   // Get the importer's data
00471   const int *RemoteLIDs  = Importer.RemoteLIDs();
00472 
00473   // Get the distributor's data
00474   int NumReceives        = D->NumReceives();
00475   const int *ProcsFrom   = D->ProcsFrom();
00476   const int *LengthsFrom = D->LengthsFrom();
00477 
00478   // Resize the outgoing data structure
00479   gpids.resize(N);
00480 
00481   // Start by claiming that I own all the data
00482   if(use_minus_one_for_local)
00483     for(i=0;i <N; i++) gpids[i]=std::make_pair(-1,Importer.TargetMap().GID(i));
00484   else
00485     for(i=0;i <N; i++) gpids[i]=std::make_pair(mypid,Importer.TargetMap().GID(i));
00486 
00487   // Now, for each remote ID, record who actually owns it.  This loop follows the operation order in the
00488   // MpiDistributor so it ought to duplicate that effect.
00489   for(i=0,j=0;i<NumReceives;i++){
00490     int pid=ProcsFrom[i];
00491     for(k=0;k<LengthsFrom[i];k++){
00492       if(pid!=mypid) gpids[RemoteLIDs[j]].first=pid;
00493       j++;
00494     }    
00495   }
00496   return 0;
00497 #else
00498   EPETRA_CHK_ERR(-10);
00499 #endif
00500 }
00501 #endif
00502 
00503 //----------------------------------------------------------------------------
00504 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00505 int Epetra_Util::GetPidGidPairs(const Epetra_Import & Importer,std::vector< std::pair<int,long long> > & gpids, bool use_minus_one_for_local){
00506   // 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.
00507   // This only works if we have an MpiDistributor in our Importer.  Otheriwise return an error.
00508 #ifdef HAVE_MPI
00509   Epetra_MpiDistributor *D=dynamic_cast<Epetra_MpiDistributor*>(&Importer.Distributor());
00510   if(!D) EPETRA_CHK_ERR(-2);
00511 
00512   int i,j,k;
00513   int mypid=Importer.TargetMap().Comm().MyPID();
00514   int N=Importer.TargetMap().NumMyElements();
00515 
00516   // Get the importer's data
00517   const int *RemoteLIDs  = Importer.RemoteLIDs();
00518 
00519   // Get the distributor's data
00520   int NumReceives        = D->NumReceives();
00521   const int *ProcsFrom   = D->ProcsFrom();
00522   const int *LengthsFrom = D->LengthsFrom();
00523 
00524   // Resize the outgoing data structure
00525   gpids.resize(N);
00526 
00527   // Start by claiming that I own all the data
00528   if(use_minus_one_for_local)
00529     for(i=0;i <N; i++) gpids[i]=std::make_pair(-1,Importer.TargetMap().GID64(i));
00530   else
00531     for(i=0;i <N; i++) gpids[i]=std::make_pair(mypid,Importer.TargetMap().GID64(i));
00532 
00533   // Now, for each remote ID, record who actually owns it.  This loop follows the operation order in the
00534   // MpiDistributor so it ought to duplicate that effect.
00535   for(i=0,j=0;i<NumReceives;i++){
00536     int pid=ProcsFrom[i];
00537     for(k=0;k<LengthsFrom[i];k++){
00538       if(pid!=mypid) gpids[RemoteLIDs[j]].first=pid;
00539       j++;
00540     }    
00541   }
00542   return 0;
00543 #else
00544   EPETRA_CHK_ERR(-10);
00545 #endif
00546 }
00547 #endif
00548 
00549 
00550 //----------------------------------------------------------------------------
00551 int Epetra_Util::GetPids(const Epetra_Import & Importer, std::vector<int> &pids, bool use_minus_one_for_local){
00552 #ifdef HAVE_MPI
00553   Epetra_MpiDistributor *D=dynamic_cast<Epetra_MpiDistributor*>(&Importer.Distributor());
00554   if(!D) EPETRA_CHK_ERR(-2);
00555 
00556   int i,j,k;
00557   int mypid=Importer.TargetMap().Comm().MyPID();
00558   int N=Importer.TargetMap().NumMyElements();
00559 
00560   // Get the importer's data
00561   const int *RemoteLIDs  = Importer.RemoteLIDs();
00562 
00563   // Get the distributor's data
00564   int NumReceives        = D->NumReceives();
00565   const int *ProcsFrom   = D->ProcsFrom();
00566   const int *LengthsFrom = D->LengthsFrom();
00567   
00568   // Resize the outgoing data structure
00569   pids.resize(N);
00570 
00571   // Start by claiming that I own all the data
00572   if(use_minus_one_for_local)
00573     for(i=0; i<N; i++) pids[i]=-1;
00574   else
00575     for(i=0; i<N; i++) pids[i]=mypid;
00576 
00577   // Now, for each remote ID, record who actually owns it.  This loop follows the operation order in the
00578   // MpiDistributor so it ought to duplicate that effect.
00579   for(i=0,j=0;i<NumReceives;i++){
00580     int pid=ProcsFrom[i];
00581     for(k=0;k<LengthsFrom[i];k++){
00582       if(pid!=mypid) pids[RemoteLIDs[j]]=pid;
00583       j++;
00584     }    
00585   }
00586   return 0;
00587 #else
00588   EPETRA_CHK_ERR(-10);
00589 #endif
00590 }
00591 
00592 //----------------------------------------------------------------------------
00593 int Epetra_Util::GetRemotePIDs(const Epetra_Import & Importer, std::vector<int> &RemotePIDs){
00594 #ifdef HAVE_MPI
00595   Epetra_MpiDistributor *D=dynamic_cast<Epetra_MpiDistributor*>(&Importer.Distributor());
00596   if(!D) {
00597     RemotePIDs.resize(0); 
00598     return 0;
00599   }
00600 
00601   int i,j,k;
00602 
00603   // Get the distributor's data
00604   int NumReceives        = D->NumReceives();
00605   const int *ProcsFrom   = D->ProcsFrom();
00606   const int *LengthsFrom = D->LengthsFrom();
00607   
00608   // Resize the outgoing data structure
00609   RemotePIDs.resize(Importer.NumRemoteIDs());
00610 
00611   // Now, for each remote ID, record who actually owns it.  This loop follows the operation order in the
00612   // MpiDistributor so it ought to duplicate that effect.
00613   for(i=0,j=0;i<NumReceives;i++){
00614     int pid=ProcsFrom[i];    
00615     for(k=0;k<LengthsFrom[i];k++){
00616       RemotePIDs[j]=pid;
00617       j++;
00618     }    
00619   }
00620   return 0;
00621 #else
00622   RemotePIDs.resize(0);
00623   return 0;
00624 #endif
00625 }
00626 
00627 //----------------------------------------------------------------------------
00628 template<typename T>
00629 int Epetra_Util_binary_search(T item,
00630                               const T* list,
00631                               int len,
00632                               int& insertPoint)
00633 {
00634   if (len < 1) {
00635     insertPoint = 0;
00636     return(-1);
00637   }
00638 
00639   unsigned start = 0, end = len - 1;
00640 
00641   while(end - start > 1) {
00642     unsigned mid = (start + end) >> 1;
00643     if (list[mid] < item) start = mid;
00644     else end = mid;
00645   }
00646 
00647   if (list[start] == item) return(start);
00648   if (list[end] == item) return(end);
00649 
00650   if (list[end] < item) {
00651     insertPoint = end+1;
00652     return(-1);
00653   }
00654 
00655   if (list[start] < item) insertPoint = end;
00656   else insertPoint = start;
00657 
00658   return(-1);
00659 }
00660 
00661 int Epetra_Util_binary_search(int item,
00662                               const int* list,
00663                               int len,
00664                               int& insertPoint)
00665 {
00666   return Epetra_Util_binary_search<int>(item, list, len, insertPoint);
00667 }
00668 
00669 //----------------------------------------------------------------------------
00670 int Epetra_Util_binary_search(long long item,
00671                               const long long* list,
00672                               int len,
00673                               int& insertPoint)
00674 {
00675   return Epetra_Util_binary_search<long long>(item, list, len, insertPoint);
00676 }
00677 
00678 //----------------------------------------------------------------------------
00679 template<typename T>
00680 int Epetra_Util_binary_search_aux(T item,
00681                               const int* list,
00682                               const T* aux_list,
00683                               int len,
00684                               int& insertPoint)
00685 {
00686   if (len < 1) {
00687     insertPoint = 0;
00688     return(-1);
00689   }
00690 
00691   unsigned start = 0, end = len - 1;
00692 
00693   while(end - start > 1) {
00694     unsigned mid = (start + end) >> 1;
00695     if (aux_list[list[mid]] < item) start = mid;
00696     else end = mid;
00697   }
00698 
00699   if (aux_list[list[start]] == item) return(start);
00700   if (aux_list[list[end]] == item) return(end);
00701 
00702   if (aux_list[list[end]] < item) {
00703     insertPoint = end+1;
00704     return(-1);
00705   }
00706 
00707   if (aux_list[list[start]] < item) insertPoint = end;
00708   else insertPoint = start;
00709 
00710   return(-1);
00711 }
00712 
00713 //----------------------------------------------------------------------------
00714 int Epetra_Util_binary_search_aux(int item,
00715                               const int* list,
00716                               const int* aux_list,
00717                               int len,
00718                               int& insertPoint)
00719 {
00720   return Epetra_Util_binary_search_aux<int>(item, list, aux_list, len, insertPoint);
00721 }
00722 
00723 //----------------------------------------------------------------------------
00724 int Epetra_Util_binary_search_aux(long long item,
00725                               const int* list,
00726                               const long long* aux_list,
00727                               int len,
00728                               int& insertPoint)
00729 {
00730   return Epetra_Util_binary_search_aux<long long>(item, list, aux_list, len, insertPoint);
00731 }
00732 
00733 
00734 
00735 //=========================================================================
00736 int Epetra_Util_ExtractHbData(Epetra_CrsMatrix * A, Epetra_MultiVector * LHS,
00737             Epetra_MultiVector * RHS,
00738             int & M, int & N, int & nz, int * & ptr,
00739             int * & ind, double * & val, int & Nrhs,
00740             double * & rhs, int & ldrhs,
00741             double * & lhs, int & ldlhs) {
00742 
00743   int ierr = 0;
00744   if (A==0) EPETRA_CHK_ERR(-1); // This matrix is defined
00745   if (!A->IndicesAreContiguous()) { // Data must be contiguous for this to work
00746     EPETRA_CHK_ERR(A->MakeDataContiguous()); // Call MakeDataContiguous() method on the matrix
00747     ierr = 1; // Warn User that we changed the matrix
00748   }
00749   
00750   M = A->NumMyRows();
00751   N = A->NumMyCols();
00752   nz = A->NumMyNonzeros();
00753   val = (*A)[0];        // Dangerous, but cheap and effective way to access first element in 
00754   
00755   const Epetra_CrsGraph & Graph = A->Graph();
00756   ind = Graph[0];  // list of values and indices
00757   
00758   Nrhs = 0; // Assume no rhs, lhs
00759 
00760   if (RHS!=0) {
00761     Nrhs = RHS->NumVectors();
00762     if (Nrhs>1)
00763     if (!RHS->ConstantStride()) {EPETRA_CHK_ERR(-2)}; // Must have strided vectors
00764     ldrhs = RHS->Stride();
00765     rhs = (*RHS)[0]; // Dangerous but effective (again)
00766   }
00767   if (LHS!=0) {
00768     int Nlhs = LHS->NumVectors();
00769     if (Nlhs!=Nrhs) {EPETRA_CHK_ERR(-3)}; // Must have same number of rhs and lhs
00770     if (Nlhs>1)
00771     if (!LHS->ConstantStride()) {EPETRA_CHK_ERR(-4)}; // Must have strided vectors
00772   ldlhs = LHS->Stride();
00773   lhs = (*LHS)[0];
00774   }
00775   
00776   // Finally build ptr vector
00777   
00778   if (ptr==0) {
00779     ptr = new int[M+1];
00780     ptr[0] = 0;
00781     for (int i=0; i<M; i++) ptr[i+1] = ptr[i] + Graph.NumMyIndices(i);
00782   }
00783   EPETRA_CHK_ERR(ierr);
00784   return(0);
00785 }
00786 
00787 // Explicitly instantiate these two even though a call to them is made.
00788 // Possible fix for a bug reported by Kenneth Belcourt.
00789 template void Epetra_Util::Sort<int>      (bool, int, int *,       int, double **, int, int **, int, long long **);
00790 template void Epetra_Util::Sort<long long>(bool, int, long long *, int, double **, int, int **, int, long long **);
00791 
00792 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines