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_Util.h"
00044 #include "Epetra_Object.h"
00045 #include "Epetra_Comm.h"
00046 #include "Epetra_Directory.h"
00047 #include "Epetra_Map.h"
00048 #include "Epetra_LocalMap.h"
00049 #include "Epetra_CrsGraph.h"
00050 #include "Epetra_CrsMatrix.h"
00051 #include "Epetra_MultiVector.h"
00052 #include "Epetra_IntVector.h"
00053 #include "Epetra_Import.h"
00054 
00055 const double Epetra_Util::chopVal_ = 1.0e-15;
00056 
00057 //=========================================================================
00058 double Epetra_Util::Chop(const double & Value){
00059   if (std::abs(Value) < chopVal_) return 0;
00060   return Value;
00061 }
00062 
00063 //=========================================================================
00064 unsigned int Epetra_Util::RandomInt() {
00065 
00066   const int a = 16807;
00067   const int m = 2147483647;
00068   const int q = 127773;
00069   const int r = 2836;
00070 
00071   int hi = Seed_ / q;
00072   int lo = Seed_ % q;
00073   int test = a * lo - r * hi;
00074   if (test > 0)
00075     Seed_ = test;
00076   else
00077     Seed_ = test + m;
00078   
00079   return(Seed_);
00080 }
00081 
00082 //=========================================================================
00083 double Epetra_Util::RandomDouble() {
00084   const double Modulus = 2147483647.0;
00085   const double DbleOne = 1.0;
00086   const double DbleTwo = 2.0;
00087 
00088   double randdouble = RandomInt(); // implicit conversion from int to double
00089   randdouble = DbleTwo * (randdouble / Modulus) - DbleOne; // scale to (-1.0, 1.0)
00090 
00091   return(randdouble);
00092 }
00093 
00094 //=========================================================================
00095 unsigned int Epetra_Util::Seed() const {
00096   return(Seed_);
00097 }
00098 
00099 //=========================================================================
00100 int Epetra_Util::SetSeed(unsigned int Seed_in) {
00101   Seed_ = Seed_in;
00102   return(0);
00103 }
00104 
00105 //=============================================================================
00106 void Epetra_Util::Sort(bool SortAscending, int NumKeys, int * Keys, 
00107            int NumDoubleCompanions,double ** DoubleCompanions, 
00108            int NumIntCompanions, int ** IntCompanions)
00109 {
00110   int i;
00111 
00112   int n = NumKeys;
00113   int * const list = Keys;
00114   int m = n/2;
00115   
00116   while (m > 0) {
00117     int max = n - m;
00118     for (int j=0; j<max; j++)
00119       {
00120   for (int k=j; k>=0; k-=m)
00121     {
00122       if ((SortAscending && list[k+m] >= list[k]) || 
00123     ( !SortAscending && list[k+m] <= list[k]))
00124         break;
00125       int temp = list[k+m];
00126       list[k+m] = list[k];
00127       list[k] = temp;
00128       for (i=0; i<NumDoubleCompanions; i++) {
00129         double dtemp = DoubleCompanions[i][k+m];
00130       DoubleCompanions[i][k+m] = DoubleCompanions[i][k];
00131       DoubleCompanions[i][k] = dtemp;
00132       }
00133       for (i=0; i<NumIntCompanions; i++) {
00134         int itemp = IntCompanions[i][k+m];
00135       IntCompanions[i][k+m] = IntCompanions[i][k];
00136       IntCompanions[i][k] = itemp;
00137       }
00138     }
00139       }
00140     m = m/2;
00141   }
00142 }
00143 
00144 //----------------------------------------------------------------------------
00145 Epetra_Map
00146 Epetra_Util::Create_OneToOne_Map(const Epetra_Map& usermap,
00147          bool high_rank_proc_owns_shared)
00148 {
00149   //if usermap is already 1-to-1 then we'll just return a copy of it.
00150   if (usermap.IsOneToOne()) {
00151     Epetra_Map newmap(usermap);
00152     return(newmap);
00153   }
00154 
00155   int myPID = usermap.Comm().MyPID();
00156   Epetra_Directory* directory = usermap.Comm().CreateDirectory(usermap);
00157 
00158   int numMyElems = usermap.NumMyElements();
00159   const int* myElems = usermap.MyGlobalElements();
00160 
00161   int* owner_procs = new int[numMyElems];
00162 
00163   directory->GetDirectoryEntries(usermap, numMyElems, myElems, owner_procs,
00164          0, 0, high_rank_proc_owns_shared);
00165 
00166   //we'll fill a list of map-elements which belong on this processor
00167 
00168   int* myOwnedElems = new int[numMyElems];
00169   int numMyOwnedElems = 0;
00170 
00171   for(int i=0; i<numMyElems; ++i) {
00172     int GID = myElems[i];
00173     int owner = owner_procs[i];
00174 
00175     if (myPID == owner) {
00176       myOwnedElems[numMyOwnedElems++] = GID;
00177     }
00178   }
00179 
00180   Epetra_Map one_to_one_map(-1, numMyOwnedElems, myOwnedElems,
00181        usermap.IndexBase(), usermap.Comm());
00182 
00183   delete [] myOwnedElems;
00184   delete [] owner_procs;
00185   delete directory;
00186 
00187   return(one_to_one_map);
00188 }
00189 
00190 //----------------------------------------------------------------------------
00191 Epetra_Map
00192 Epetra_Util::Create_Root_Map(const Epetra_Map& usermap,
00193          int root)
00194 {
00195 
00196   int numProc = usermap.Comm().NumProc();
00197   if (numProc==1) {
00198     Epetra_Map newmap(usermap);
00199     return(newmap);
00200   }
00201 
00202   const Epetra_Comm & comm = usermap.Comm();
00203   bool isRoot = usermap.Comm().MyPID()==root;
00204 
00205   //if usermap is already completely owned by root then we'll just return a copy of it.
00206   int quickreturn = 0;
00207   int globalquickreturn = 0;
00208 
00209   if (isRoot) {
00210     if (usermap.NumMyElements()==usermap.NumGlobalElements()) quickreturn = 1;
00211   }
00212   else {
00213     if (usermap.NumMyElements()==0) quickreturn = 1;
00214   }
00215   usermap.Comm().MinAll(&quickreturn, &globalquickreturn, 1);
00216   
00217   if (globalquickreturn==1) {
00218     Epetra_Map newmap(usermap);
00219     return(newmap);
00220   }
00221   
00222   // Linear map: Simple case, just put all GIDs linearly on root processor
00223   if (usermap.LinearMap() && root!=-1) {
00224     int numMyElements = 0;
00225     if (isRoot) numMyElements = usermap.MaxAllGID()+1;
00226     Epetra_Map newmap(-1, numMyElements, usermap.IndexBase(), comm);
00227     return(newmap);
00228   }
00229 
00230   if (!usermap.UniqueGIDs()) 
00231     throw usermap.ReportError("usermap must have unique GIDs",-1);
00232 
00233   // General map
00234 
00235   // Build IntVector of the GIDs, then ship them to root processor
00236   int numMyElements = usermap.NumMyElements();
00237   Epetra_Map allGidsMap(-1, numMyElements, 0, comm);
00238   Epetra_IntVector allGids(allGidsMap);
00239   for (int i=0; i<numMyElements; i++) allGids[i] = usermap.GID(i);
00240   
00241   int numGlobalElements = usermap.NumGlobalElements();
00242   if (root!=-1) {
00243     int n1 = 0; if (isRoot) n1 = numGlobalElements;
00244     Epetra_Map allGidsOnRootMap(-1, n1, 0, comm);
00245     Epetra_Import importer(allGidsOnRootMap, allGidsMap);
00246     Epetra_IntVector allGidsOnRoot(allGidsOnRootMap);
00247     allGidsOnRoot.Import(allGids, importer, Insert);
00248     
00249     Epetra_Map rootMap(-1, allGidsOnRoot.MyLength(), allGidsOnRoot.Values(), usermap.IndexBase(), comm);
00250     return(rootMap);
00251   }
00252   else {
00253     int n1 = numGlobalElements;
00254     Epetra_LocalMap allGidsOnRootMap(n1, 0, comm);
00255     Epetra_Import importer(allGidsOnRootMap, allGidsMap);
00256     Epetra_IntVector allGidsOnRoot(allGidsOnRootMap);
00257     allGidsOnRoot.Import(allGids, importer, Insert);
00258     
00259     Epetra_Map rootMap(-1, allGidsOnRoot.MyLength(), allGidsOnRoot.Values(), usermap.IndexBase(), comm);
00260 
00261     return(rootMap);
00262   }
00263 }
00264 
00265 //----------------------------------------------------------------------------
00266 Epetra_BlockMap
00267 Epetra_Util::Create_OneToOne_BlockMap(const Epetra_BlockMap& usermap,
00268               bool high_rank_proc_owns_shared)
00269 {
00270   //if usermap is already 1-to-1 then we'll just return a copy of it.
00271   if (usermap.IsOneToOne()) {
00272     Epetra_BlockMap newmap(usermap);
00273     return(newmap);
00274   }
00275 
00276   int myPID = usermap.Comm().MyPID();
00277   Epetra_Directory* directory = usermap.Comm().CreateDirectory(usermap);
00278 
00279   int numMyElems = usermap.NumMyElements();
00280   const int* myElems = usermap.MyGlobalElements();
00281 
00282   int* owner_procs = new int[numMyElems*2];
00283   int* sizes = owner_procs+numMyElems;
00284 
00285   directory->GetDirectoryEntries(usermap, numMyElems, myElems, owner_procs,
00286          0, sizes, high_rank_proc_owns_shared);
00287 
00288   //we'll fill a list of map-elements which belong on this processor
00289 
00290   int* myOwnedElems = new int[numMyElems*2];
00291   int* ownedSizes = myOwnedElems+numMyElems;
00292   int numMyOwnedElems = 0;
00293 
00294   for(int i=0; i<numMyElems; ++i) {
00295     int GID = myElems[i];
00296     int owner = owner_procs[i];
00297 
00298     if (myPID == owner) {
00299       ownedSizes[numMyOwnedElems] = sizes[i];
00300       myOwnedElems[numMyOwnedElems++] = GID;
00301     }
00302   }
00303 
00304   Epetra_BlockMap one_to_one_map(-1, numMyOwnedElems, myOwnedElems,
00305          sizes, usermap.IndexBase(), usermap.Comm());
00306 
00307   delete [] myOwnedElems;
00308   delete [] owner_procs;
00309   delete directory;
00310 
00311   return(one_to_one_map);
00312 }
00313 
00314 //----------------------------------------------------------------------------
00315 int Epetra_Util_binary_search(int item,
00316                               const int* list,
00317                               int len,
00318                               int& insertPoint)
00319 {
00320   if (len < 1) {
00321     insertPoint = 0;
00322     return(-1);
00323   }
00324 
00325   unsigned start = 0, end = len - 1;
00326 
00327   while(end - start > 1) {
00328     unsigned mid = (start + end) >> 1;
00329     if (list[mid] < item) start = mid;
00330     else end = mid;
00331   }
00332 
00333   if (list[start] == item) return(start);
00334   if (list[end] == item) return(end);
00335 
00336   if (list[end] < item) {
00337     insertPoint = end+1;
00338     return(-1);
00339   }
00340 
00341   if (list[start] < item) insertPoint = end;
00342   else insertPoint = start;
00343 
00344   return(-1);
00345 }
00346 
00347 //=========================================================================
00348 int Epetra_Util_ExtractHbData(Epetra_CrsMatrix * A, Epetra_MultiVector * LHS,
00349             Epetra_MultiVector * RHS,
00350             int & M, int & N, int & nz, int * & ptr,
00351             int * & ind, double * & val, int & Nrhs,
00352             double * & rhs, int & ldrhs,
00353             double * & lhs, int & ldlhs) {
00354 
00355   int ierr = 0;
00356   if (A==0) EPETRA_CHK_ERR(-1); // This matrix is defined
00357   if (!A->IndicesAreContiguous()) { // Data must be contiguous for this to work
00358     EPETRA_CHK_ERR(A->MakeDataContiguous()); // Call MakeDataContiguous() method on the matrix
00359     ierr = 1; // Warn User that we changed the matrix
00360   }
00361   
00362   M = A->NumMyRows();
00363   N = A->NumMyCols();
00364   nz = A->NumMyNonzeros();
00365   val = (*A)[0];        // Dangerous, but cheap and effective way to access first element in 
00366   
00367   const Epetra_CrsGraph & Graph = A->Graph();
00368   ind = Graph[0];  // list of values and indices
00369   
00370   Nrhs = 0; // Assume no rhs, lhs
00371 
00372   if (RHS!=0) {
00373     Nrhs = RHS->NumVectors();
00374     if (Nrhs>1)
00375     if (!RHS->ConstantStride()) {EPETRA_CHK_ERR(-2)}; // Must have strided vectors
00376     ldrhs = RHS->Stride();
00377     rhs = (*RHS)[0]; // Dangerous but effective (again)
00378   }
00379   if (LHS!=0) {
00380     int Nlhs = LHS->NumVectors();
00381     if (Nlhs!=Nrhs) {EPETRA_CHK_ERR(-3)}; // Must have same number of rhs and lhs
00382     if (Nlhs>1)
00383     if (!LHS->ConstantStride()) {EPETRA_CHK_ERR(-4)}; // Must have strided vectors
00384   ldlhs = LHS->Stride();
00385   lhs = (*LHS)[0];
00386   }
00387   
00388   // Finally build ptr vector
00389   
00390   if (ptr==0) {
00391     ptr = new int[M+1];
00392     ptr[0] = 0;
00393     for (int i=0; i<M; i++) ptr[i+1] = ptr[i] + Graph.NumMyIndices(i);
00394   }
00395   EPETRA_CHK_ERR(ierr);
00396   return(0);
00397 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines