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