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 terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00009 // license for use of this work by or on behalf of the U.S. Government.
00010 // 
00011 // This library is free software; you can redistribute it and/or modify
00012 // it under the terms of the GNU Lesser General Public License as
00013 // published by the Free Software Foundation; either version 2.1 of the
00014 // License, or (at your option) any later version.
00015 //  
00016 // This library is distributed in the hope that it will be useful, but
00017 // WITHOUT ANY WARRANTY; without even the implied warranty of
00018 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00019 // Lesser General Public License for more details.
00020 //  
00021 // You should have received a copy of the GNU Lesser General Public
00022 // License along with this library; if not, write to the Free Software
00023 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00024 // USA
00025 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00026 // 
00027 // ************************************************************************
00028 //@HEADER
00029 
00030 #include "Epetra_Util.h"
00031 #include "Epetra_Object.h"
00032 #include "Epetra_Comm.h"
00033 #include "Epetra_Directory.h"
00034 #include "Epetra_Map.h"
00035 #include "Epetra_LocalMap.h"
00036 #include "Epetra_CrsGraph.h"
00037 #include "Epetra_CrsMatrix.h"
00038 #include "Epetra_MultiVector.h"
00039 #include "Epetra_IntVector.h"
00040 #include "Epetra_Import.h"
00041 
00042 const double Epetra_Util::chopVal_ = 1.0e-15;
00043 
00044 //=========================================================================
00045 unsigned int Epetra_Util::RandomInt() {
00046 
00047   const int a = 16807;
00048   const int m = 2147483647;
00049   const int q = 127773;
00050   const int r = 2836;
00051 
00052   int hi = Seed_ / q;
00053   int lo = Seed_ % q;
00054   int test = a * lo - r * hi;
00055   if (test > 0)
00056     Seed_ = test;
00057   else
00058     Seed_ = test + m;
00059   
00060   return(Seed_);
00061 }
00062 
00063 //=========================================================================
00064 double Epetra_Util::RandomDouble() {
00065   const double Modulus = 2147483647.0;
00066   const double DbleOne = 1.0;
00067   const double DbleTwo = 2.0;
00068 
00069   double randdouble = RandomInt(); // implicit conversion from int to double
00070   randdouble = DbleTwo * (randdouble / Modulus) - DbleOne; // scale to (-1.0, 1.0)
00071 
00072   return(randdouble);
00073 }
00074 
00075 //=========================================================================
00076 unsigned int Epetra_Util::Seed() const {
00077   return(Seed_);
00078 }
00079 
00080 //=========================================================================
00081 int Epetra_Util::SetSeed(unsigned int Seed) {
00082   Seed_ = Seed;
00083   return(0);
00084 }
00085 
00086 //=============================================================================
00087 void Epetra_Util::Sort(bool SortAscending, int NumKeys, int * Keys, 
00088            int NumDoubleCompanions,double ** DoubleCompanions, 
00089            int NumIntCompanions, int ** IntCompanions)
00090 {
00091   int i;
00092 
00093   int n = NumKeys;
00094   int * const list = Keys;
00095   int m = n/2;
00096   
00097   while (m > 0) {
00098     int max = n - m;
00099     for (int j=0; j<max; j++)
00100       {
00101   for (int k=j; k>=0; k-=m)
00102     {
00103       if ((SortAscending && list[k+m] >= list[k]) || 
00104     ( !SortAscending && list[k+m] <= list[k]))
00105         break;
00106       int temp = list[k+m];
00107       list[k+m] = list[k];
00108       list[k] = temp;
00109       for (i=0; i<NumDoubleCompanions; i++) {
00110         double dtemp = DoubleCompanions[i][k+m];
00111       DoubleCompanions[i][k+m] = DoubleCompanions[i][k];
00112       DoubleCompanions[i][k] = dtemp;
00113       }
00114       for (i=0; i<NumIntCompanions; i++) {
00115         int itemp = IntCompanions[i][k+m];
00116       IntCompanions[i][k+m] = IntCompanions[i][k];
00117       IntCompanions[i][k] = itemp;
00118       }
00119     }
00120       }
00121     m = m/2;
00122   }
00123 }
00124 
00125 //----------------------------------------------------------------------------
00126 Epetra_Map
00127 Epetra_Util::Create_OneToOne_Map(const Epetra_Map& usermap,
00128          bool high_rank_proc_owns_shared)
00129 {
00130   //if usermap is already 1-to-1 then we'll just return a copy of it.
00131   if (usermap.IsOneToOne()) {
00132     Epetra_Map newmap(usermap);
00133     return(newmap);
00134   }
00135 
00136   int myPID = usermap.Comm().MyPID();
00137   Epetra_Directory* directory = usermap.Comm().CreateDirectory(usermap);
00138 
00139   int numMyElems = usermap.NumMyElements();
00140   const int* myElems = usermap.MyGlobalElements();
00141 
00142   int* owner_procs = new int[numMyElems];
00143 
00144   directory->GetDirectoryEntries(usermap, numMyElems, myElems, owner_procs,
00145          0, 0, high_rank_proc_owns_shared);
00146 
00147   //we'll fill a list of map-elements which belong on this processor
00148 
00149   int* myOwnedElems = new int[numMyElems];
00150   int numMyOwnedElems = 0;
00151 
00152   for(int i=0; i<numMyElems; ++i) {
00153     int GID = myElems[i];
00154     int owner = owner_procs[i];
00155 
00156     if (myPID == owner) {
00157       myOwnedElems[numMyOwnedElems++] = GID;
00158     }
00159   }
00160 
00161   Epetra_Map one_to_one_map(-1, numMyOwnedElems, myOwnedElems,
00162        usermap.IndexBase(), usermap.Comm());
00163 
00164   delete [] myOwnedElems;
00165   delete [] owner_procs;
00166   delete directory;
00167 
00168   return(one_to_one_map);
00169 }
00170 
00171 //----------------------------------------------------------------------------
00172 Epetra_Map
00173 Epetra_Util::Create_Root_Map(const Epetra_Map& usermap,
00174          int root)
00175 {
00176 
00177   int numProc = usermap.Comm().NumProc();
00178   if (numProc==1) {
00179     Epetra_Map newmap(usermap);
00180     return(newmap);
00181   }
00182 
00183   const Epetra_Comm & comm = usermap.Comm();
00184   bool isRoot = usermap.Comm().MyPID()==root;
00185 
00186   //if usermap is already completely owned by root then we'll just return a copy of it.
00187   int quickreturn = 0;
00188   int globalquickreturn = 0;
00189 
00190   if (isRoot) {
00191     if (usermap.NumMyElements()==usermap.NumGlobalElements()) quickreturn = 1;
00192   }
00193   else {
00194     if (usermap.NumMyElements()==0) quickreturn = 1;
00195   }
00196   usermap.Comm().MinAll(&quickreturn, &globalquickreturn, 1);
00197   
00198   if (globalquickreturn==1) {
00199     Epetra_Map newmap(usermap);
00200     return(newmap);
00201   }
00202   
00203   // Linear map: Simple case, just put all GIDs linearly on root processor
00204   if (usermap.LinearMap() && root!=-1) {
00205     int numMyElements = 0;
00206     if (isRoot) numMyElements = usermap.MaxAllGID()+1;
00207     Epetra_Map newmap(-1, numMyElements, usermap.IndexBase(), comm);
00208     return(newmap);
00209   }
00210 
00211   if (!usermap.UniqueGIDs()) 
00212     throw usermap.ReportError("usermap must have unique GIDs",-1);
00213 
00214   // General map
00215 
00216   // Build IntVector of the GIDs, then ship them to root processor
00217   int numMyElements = usermap.NumMyElements();
00218   Epetra_Map allGidsMap(-1, numMyElements, 0, comm);
00219   Epetra_IntVector allGids(allGidsMap);
00220   for (int i=0; i<numMyElements; i++) allGids[i] = usermap.GID(i);
00221   
00222   int numGlobalElements = usermap.NumGlobalElements();
00223   if (root!=-1) {
00224     int n1 = 0; if (isRoot) n1 = numGlobalElements;
00225     Epetra_Map allGidsOnRootMap(-1, n1, 0, comm);
00226     Epetra_Import importer(allGidsOnRootMap, allGidsMap);
00227     Epetra_IntVector allGidsOnRoot(allGidsOnRootMap);
00228     allGidsOnRoot.Import(allGids, importer, Insert);
00229     
00230     Epetra_Map rootMap(-1, allGidsOnRoot.MyLength(), allGidsOnRoot.Values(), usermap.IndexBase(), comm);
00231     return(rootMap);
00232   }
00233   else {
00234     int n1 = numGlobalElements;
00235     Epetra_LocalMap allGidsOnRootMap(n1, 0, comm);
00236     Epetra_Import importer(allGidsOnRootMap, allGidsMap);
00237     Epetra_IntVector allGidsOnRoot(allGidsOnRootMap);
00238     allGidsOnRoot.Import(allGids, importer, Insert);
00239     
00240     Epetra_Map rootMap(-1, allGidsOnRoot.MyLength(), allGidsOnRoot.Values(), usermap.IndexBase(), comm);
00241 
00242     return(rootMap);
00243   }
00244 }
00245 
00246 //----------------------------------------------------------------------------
00247 Epetra_BlockMap
00248 Epetra_Util::Create_OneToOne_BlockMap(const Epetra_BlockMap& usermap,
00249               bool high_rank_proc_owns_shared)
00250 {
00251   //if usermap is already 1-to-1 then we'll just return a copy of it.
00252   if (usermap.IsOneToOne()) {
00253     Epetra_BlockMap newmap(usermap);
00254     return(newmap);
00255   }
00256 
00257   int myPID = usermap.Comm().MyPID();
00258   Epetra_Directory* directory = usermap.Comm().CreateDirectory(usermap);
00259 
00260   int numMyElems = usermap.NumMyElements();
00261   const int* myElems = usermap.MyGlobalElements();
00262 
00263   int* owner_procs = new int[numMyElems*2];
00264   int* sizes = owner_procs+numMyElems;
00265 
00266   directory->GetDirectoryEntries(usermap, numMyElems, myElems, owner_procs,
00267          0, sizes, high_rank_proc_owns_shared);
00268 
00269   //we'll fill a list of map-elements which belong on this processor
00270 
00271   int* myOwnedElems = new int[numMyElems*2];
00272   int* ownedSizes = myOwnedElems+numMyElems;
00273   int numMyOwnedElems = 0;
00274 
00275   for(int i=0; i<numMyElems; ++i) {
00276     int GID = myElems[i];
00277     int owner = owner_procs[i];
00278 
00279     if (myPID == owner) {
00280       ownedSizes[numMyOwnedElems] = sizes[i];
00281       myOwnedElems[numMyOwnedElems++] = GID;
00282     }
00283   }
00284 
00285   Epetra_BlockMap one_to_one_map(-1, numMyOwnedElems, myOwnedElems,
00286          sizes, usermap.IndexBase(), usermap.Comm());
00287 
00288   delete [] myOwnedElems;
00289   delete [] owner_procs;
00290   delete directory;
00291 
00292   return(one_to_one_map);
00293 }
00294 
00295 //----------------------------------------------------------------------------
00296 int Epetra_Util_binary_search(int item,
00297                               const int* list,
00298                               int len,
00299                               int& insertPoint)
00300 {
00301   if (len < 1) {
00302     insertPoint = 0;
00303     return(-1);
00304   }
00305 
00306   unsigned start = 0, end = len - 1;
00307 
00308   while(end - start > 1) {
00309     unsigned mid = (start + end) >> 1;
00310     if (list[mid] < item) start = mid;
00311     else end = mid;
00312   }
00313 
00314   if (list[start] == item) return(start);
00315   if (list[end] == item) return(end);
00316 
00317   if (list[end] < item) {
00318     insertPoint = end+1;
00319     return(-1);
00320   }
00321 
00322   if (list[start] < item) insertPoint = end;
00323   else insertPoint = start;
00324 
00325   return(-1);
00326 }
00327 
00328 //=========================================================================
00329 int Epetra_Util_ExtractHbData(Epetra_CrsMatrix * A, Epetra_MultiVector * LHS,
00330             Epetra_MultiVector * RHS,
00331             int & M, int & N, int & nz, int * & ptr,
00332             int * & ind, double * & val, int & Nrhs,
00333             double * & rhs, int & ldrhs,
00334             double * & lhs, int & ldlhs) {
00335 
00336   int ierr = 0;
00337   if (A==0) EPETRA_CHK_ERR(-1); // This matrix is defined
00338   if (!A->IndicesAreContiguous()) { // Data must be contiguous for this to work
00339     EPETRA_CHK_ERR(A->MakeDataContiguous()); // Call MakeDataContiguous() method on the matrix
00340     ierr = 1; // Warn User that we changed the matrix
00341   }
00342   
00343   M = A->NumMyRows();
00344   N = A->NumMyCols();
00345   nz = A->NumMyNonzeros();
00346   val = (*A)[0];        // Dangerous, but cheap and effective way to access first element in 
00347   
00348   const Epetra_CrsGraph & Graph = A->Graph();
00349   ind = Graph[0];  // list of values and indices
00350   
00351   Nrhs = 0; // Assume no rhs, lhs
00352 
00353   if (RHS!=0) {
00354     Nrhs = RHS->NumVectors();
00355     if (Nrhs>1)
00356     if (!RHS->ConstantStride()) {EPETRA_CHK_ERR(-2)}; // Must have strided vectors
00357     ldrhs = RHS->Stride();
00358     rhs = (*RHS)[0]; // Dangerous but effective (again)
00359   }
00360   if (LHS!=0) {
00361     int Nlhs = LHS->NumVectors();
00362     if (Nlhs!=Nrhs) {EPETRA_CHK_ERR(-3)}; // Must have same number of rhs and lhs
00363     if (Nlhs>1)
00364     if (!LHS->ConstantStride()) {EPETRA_CHK_ERR(-4)}; // Must have strided vectors
00365   ldlhs = LHS->Stride();
00366   lhs = (*LHS)[0];
00367   }
00368   
00369   // Finally build ptr vector
00370   
00371   if (ptr==0) {
00372     ptr = new int[M+1];
00373     ptr[0] = 0;
00374     for (int i=0; i<M; i++) ptr[i+1] = ptr[i] + Graph.NumMyIndices(i);
00375   }
00376   EPETRA_CHK_ERR(ierr);
00377   return(0);
00378 }

Generated on Thu Sep 18 12:37:58 2008 for Epetra Package Browser (Single Doxygen Collection) by doxygen 1.3.9.1