Chaste  Release::3.4
PetscTools.cpp
1 /*
2 
3 Copyright (c) 2005-2016, University of Oxford.
4 All rights reserved.
5 
6 University of Oxford means the Chancellor, Masters and Scholars of the
7 University of Oxford, having an administrative office at Wellington
8 Square, Oxford OX1 2JD, UK.
9 
10 This file is part of Chaste.
11 
12 Redistribution and use in source and binary forms, with or without
13 modification, are permitted provided that the following conditions are met:
14  * Redistributions of source code must retain the above copyright notice,
15  this list of conditions and the following disclaimer.
16  * Redistributions in binary form must reproduce the above copyright notice,
17  this list of conditions and the following disclaimer in the documentation
18  and/or other materials provided with the distribution.
19  * Neither the name of the University of Oxford nor the names of its
20  contributors may be used to endorse or promote products derived from this
21  software without specific prior written permission.
22 
23 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
24 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
25 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
26 ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
27 LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
28 CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
29 GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
30 HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
31 LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
32 OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
33 
34 */
35 
36 #include "PetscTools.hpp"
37 #include "Exception.hpp"
38 #include "Warnings.hpp"
39 #include <iostream>
40 #include <sstream>
41 #include <cassert>
42 #include <cstring> // For strcmp etc. Needed in gcc-4.3
43 
45 unsigned PetscTools::mNumProcessors = 0;
46 unsigned PetscTools::mRank = 0;
48 
49 #ifndef NDEBUG
50 // Uncomment this to trace calls to PetscTools::Barrier
51 //#define DEBUG_BARRIERS
52 #endif
53 
54 #ifdef DEBUG_BARRIERS
55 static unsigned mNumBarriers = 0u;
56 #endif
57 
59 {
60  PetscBool is_there;
61  PetscInitialized(&is_there);
62  if (is_there)
63  {
64  mPetscIsInitialised = true;
65 
66  PetscInt num_procs;
67  MPI_Comm_size(PETSC_COMM_WORLD, &num_procs);
68  mNumProcessors = (unsigned) num_procs;
69 
70  PetscInt my_rank;
71  MPI_Comm_rank(PETSC_COMM_WORLD, &my_rank);
72  mRank = (unsigned) my_rank;
73  }
74  else
75  {
76  // No PETSc
77  mPetscIsInitialised = false;
78  mNumProcessors = 1;
79  mRank = 0;
80  }
81 }
82 
83 // Information methods
84 
86 {
87  CheckCache();
88  return mPetscIsInitialised;
89 }
90 
92 {
93  CheckCache();
94  return (mIsolateProcesses || mNumProcessors == 1);
95 }
96 
98 {
99  CheckCache();
100  return (mNumProcessors > 1);
101 }
102 
104 {
105  return mIsolateProcesses;
106 }
107 
109 {
110  CheckCache();
111  return mNumProcessors;
112 }
113 
115 {
116  CheckCache();
117  return mRank;
118 }
119 
121 {
122  CheckCache();
123  return (mRank == MASTER_RANK || mIsolateProcesses);
124 }
125 
127 {
128  CheckCache();
129  return (mRank == mNumProcessors - 1 || mIsolateProcesses);
130 }
131 
132 // Little utility methods
133 
134 void PetscTools::Barrier(const std::string callerId)
135 {
136  CheckCache();
138  {
139 #ifdef DEBUG_BARRIERS
140  // "Before" is alphabetically before "Post" so that one can sort the output on process/event/barrier
141  std::cout << "DEBUG: proc " << PetscTools::GetMyRank() << ": Before " << "Barrier " << mNumBarriers << " \""<< callerId << "\"." << std::endl << std::flush;
142 #endif
143  PetscBarrier(PETSC_NULL);
144 #ifdef DEBUG_BARRIERS
145  std::cout << "DEBUG: proc " << PetscTools::GetMyRank() << ": Post " << "Barrier " << mNumBarriers++ << " \""<< callerId << "\"." << std::endl << std::flush;
146 #endif
147  }
148 }
149 
151 {
152  Barrier("PetscTools::RoundRobin"); // We want barriers both before all and after all, just in case
153  const unsigned my_rank = GetMyRank();
154  for (unsigned turn=0; turn<my_rank; turn++)
155  {
156  Barrier("PetscTools::RoundRobin");
157  }
158 }
159 
161 {
162  const unsigned num_procs = GetNumProcs();
163  for (unsigned turn=GetMyRank(); turn<num_procs; turn++)
164  {
165  Barrier("PetscTools::RoundRobin");
166  }
167 }
168 
170 {
171  mIsolateProcesses = isolate;
172 }
173 
175 {
176  if (mIsolateProcesses)
177  {
178  return PETSC_COMM_SELF;
179  }
180  else
181  {
182  return PETSC_COMM_WORLD;
183  }
184 }
185 
187 {
188  CheckCache();
189  unsigned my_flag = (unsigned) flag;
190  unsigned anyones_flag_is_true = my_flag;
192  {
193  MPI_Allreduce(&my_flag, &anyones_flag_is_true, 1, MPI_UNSIGNED, MPI_MAX, PETSC_COMM_WORLD);
194  }
195  return (anyones_flag_is_true == 1);
196 }
197 
199 {
200  bool anyones_error=ReplicateBool(flag);
201  if (flag)
202  {
203  // Return control to exception thrower
204  return;
205  }
206  if (anyones_error)
207  {
208  EXCEPTION("Another process threw an exception; bailing out.");
209  }
210 }
211 
212 // Vector & matrix creation routines
213 
214 Vec PetscTools::CreateVec(int size, int localSize, bool ignoreOffProcEntries)
215 {
216  assert(size >= 0); // There is one test where we create a zero-sized vector
217  Vec ret;
218  VecCreate(PETSC_COMM_WORLD, &ret);
219  VecSetSizes(ret, localSize, size); // localSize usually defaults to PETSC_DECIDE
220  VecSetFromOptions(ret);
221 
222  if (ignoreOffProcEntries)
223  {
224 #if (PETSC_VERSION_MAJOR == 3) //PETSc 3.x.x
225  VecSetOption(ret, VEC_IGNORE_OFF_PROC_ENTRIES, PETSC_TRUE);
226 #else
227  VecSetOption(ret, VEC_IGNORE_OFF_PROC_ENTRIES);
228 #endif
229  }
230 
231  return ret;
232 }
233 
234 Vec PetscTools::CreateVec(std::vector<double> data)
235 {
236  assert(data.size() > 0);
237  Vec ret = CreateVec(data.size());
238 
239  double* p_ret;
240  VecGetArray(ret, &p_ret);
241  int lo, hi;
242  VecGetOwnershipRange(ret, &lo, &hi);
243 
244  for (int global_index=lo; global_index<hi; global_index++)
245  {
246  int local_index = global_index - lo;
247  p_ret[local_index] = data[global_index];
248  }
249  VecRestoreArray(ret, &p_ret);
250 
251  return ret;
252 }
253 
254 Vec PetscTools::CreateAndSetVec(int size, double value)
255 {
256  assert(size > 0);
257  Vec ret = CreateVec(size);
258 
259 #if (PETSC_VERSION_MAJOR == 2 && PETSC_VERSION_MINOR == 2) //PETSc 2.2
260  VecSet(&value, ret);
261 #else
262  VecSet(ret, value);
263 #endif
264 
265  return ret;
266 }
267 
268 void PetscTools::SetupMat(Mat& rMat, int numRows, int numColumns,
269  unsigned rowPreallocation,
270  int numLocalRows,
271  int numLocalColumns,
272  bool ignoreOffProcEntries,
273  bool newAllocationError)
274 {
275  assert(numRows > 0);
276  assert(numColumns > 0);
277  if ((int) rowPreallocation>numColumns)
278  {
279  WARNING("Preallocation failure: requested number of nonzeros per row greater than number of columns");//+rowPreallocation+">"+numColumns);
280  rowPreallocation = numColumns;
281  }
282 
283 #if (PETSC_VERSION_MAJOR == 2 && PETSC_VERSION_MINOR == 2) //PETSc 2.2
284  MatCreate(PETSC_COMM_WORLD,numLocalRows,numLocalColumns,numRows,numColumns,&rMat);
285 #else //New API
286  MatCreate(PETSC_COMM_WORLD,&rMat);
287  MatSetSizes(rMat,numLocalRows,numLocalColumns,numRows,numColumns);
288 #endif
289 
291  {
292  MatSetType(rMat, MATSEQAIJ);
293  if (rowPreallocation > 0)
294  {
295  MatSeqAIJSetPreallocation(rMat, rowPreallocation, PETSC_NULL);
296  }
297  }
298  else
299  {
300  MatSetType(rMat, MATMPIAIJ);
301  if (rowPreallocation > 0)
302  {
303  MatMPIAIJSetPreallocation(rMat, rowPreallocation, PETSC_NULL, rowPreallocation, PETSC_NULL);
304  }
305  }
306 
307  MatSetFromOptions(rMat);
308 
309  if (ignoreOffProcEntries)//&& IsParallel())
310  {
311  if (rowPreallocation == 0)
312  {
313  // We aren't allowed to do non-zero allocation after setting MAT_IGNORE_OFF_PROC_ENTRIES
314  WARNING("Ignoring MAT_IGNORE_OFF_PROC_ENTRIES flag because we might set non-zeroes later");
315  }
316  else
317  {
318 #if (PETSC_VERSION_MAJOR == 3) //PETSc 3.x.x
319  MatSetOption(rMat, MAT_IGNORE_OFF_PROC_ENTRIES, PETSC_TRUE);
320 #else
321  MatSetOption(rMat, MAT_IGNORE_OFF_PROC_ENTRIES);
322 #endif
323  }
324  }
325 #if (PETSC_VERSION_MAJOR == 3 && PETSC_VERSION_MINOR >= 3) //PETSc 3.3 or later
326  if (newAllocationError == false)
327  {
328  MatSetOption(rMat, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
329  }
330 #endif
331 }
332 
333 void PetscTools::DumpPetscObject(const Mat& rMat, const std::string& rOutputFileFullPath)
334 {
335  PetscViewer view;
336 #if (PETSC_VERSION_MAJOR == 2 && PETSC_VERSION_MINOR == 2) //PETSc 2.2
337  PetscViewerFileType type = PETSC_FILE_CREATE;
338 #else
339  PetscFileMode type = FILE_MODE_WRITE;
340 #endif
341 
342  PetscViewerBinaryOpen(PETSC_COMM_WORLD, rOutputFileFullPath.c_str(), type, &view);
343  MatView(rMat, view);
344  PetscViewerDestroy(PETSC_DESTROY_PARAM(view));
345 }
346 
347 void PetscTools::DumpPetscObject(const Vec& rVec, const std::string& rOutputFileFullPath)
348 {
349  PetscViewer view;
350 #if (PETSC_VERSION_MAJOR == 2 && PETSC_VERSION_MINOR == 2) //PETSc 2.2
351  PetscViewerFileType type = PETSC_FILE_CREATE;
352 #else
353  PetscFileMode type = FILE_MODE_WRITE;
354 #endif
355 
356  PetscViewerBinaryOpen(PETSC_COMM_WORLD, rOutputFileFullPath.c_str(), type, &view);
357  VecView(rVec, view);
358  PetscViewerDestroy(PETSC_DESTROY_PARAM(view));
359 }
360 
361 void PetscTools::ReadPetscObject(Mat& rMat, const std::string& rOutputFileFullPath, Vec rParallelLayout)
362 {
363  /*
364  * PETSc (as of 3.1) doesn't provide any method for loading a Mat object with a user-defined parallel
365  * layout, i.e. there's no equivalent to VecLoadIntoVector for Mat's.
366  *
367  * It seems to be in their future plans though: http://lists.mcs.anl.gov/pipermail/petsc-users/2010-March/006062.html
368  */
369 
370  PetscViewer view;
371 #if (PETSC_VERSION_MAJOR == 2 && PETSC_VERSION_MINOR == 2) //PETSc 2.2
372  PetscViewerFileType type = PETSC_FILE_RDONLY;
373 #else
374  PetscFileMode type = FILE_MODE_READ;
375 #endif
376 
377  PetscViewerBinaryOpen(PETSC_COMM_WORLD, rOutputFileFullPath.c_str(),
378  type, &view);
379 
380 #if (PETSC_VERSION_MAJOR == 3 && PETSC_VERSION_MINOR >= 2) //PETSc 3.2 or later
381  MatCreate(PETSC_COMM_WORLD,&rMat);
382  MatSetType(rMat,MATMPIAIJ);
383  MatLoad(rMat,view);
384 #else
385  MatLoad(view, MATMPIAIJ, &rMat);
386 #endif
387 
388  PetscViewerDestroy(PETSC_DESTROY_PARAM(view));
389 
390  if (rParallelLayout != NULL)
391  {
392  /*
393  * The idea is to copy rMat into a matrix that has the appropriate
394  * parallel layout. Inefficient...
395  */
396  PetscInt num_rows, num_local_rows;
397 
398  VecGetSize(rParallelLayout, &num_rows);
399  VecGetLocalSize(rParallelLayout, &num_local_rows);
400 
401  Mat temp_mat;
403  PetscTools::SetupMat(temp_mat, num_rows, num_rows, 100, num_local_rows, num_local_rows, false);
404 
405  MatCopy(rMat, temp_mat, DIFFERENT_NONZERO_PATTERN);
406 
407  PetscTools::Destroy(rMat);
408  rMat = temp_mat;
409  }
410 }
411 
412 void PetscTools::ReadPetscObject(Vec& rVec, const std::string& rOutputFileFullPath, Vec rParallelLayout)
413 {
414  PetscViewer view;
415 #if (PETSC_VERSION_MAJOR == 2 && PETSC_VERSION_MINOR == 2) //PETSc 2.2
416  PetscViewerFileType type = PETSC_FILE_RDONLY;
417 #else
418  PetscFileMode type = FILE_MODE_READ;
419 #endif
420 
421  PetscViewerBinaryOpen(PETSC_COMM_WORLD, rOutputFileFullPath.c_str(),
422  type, &view);
423  if (rParallelLayout == NULL)
424  {
425 #if (PETSC_VERSION_MAJOR == 3 && PETSC_VERSION_MINOR >= 2) //PETSc 3.2 or later
426  VecCreate(PETSC_COMM_WORLD,&rVec);
427  VecSetType(rVec,VECMPI);
428  VecLoad(rVec,view);
429 #else
430  VecLoad(view, VECMPI, &rVec);
431 #endif
432  }
433  else
434  {
435  VecDuplicate(rParallelLayout, &rVec);
436 #if (PETSC_VERSION_MAJOR == 3 && PETSC_VERSION_MINOR >= 2) //PETSc 3.2 or later
437  VecLoad(rVec,view);
438 #else
439  VecLoadIntoVector(view, rVec);
440 #endif
441  }
442  PetscViewerDestroy(PETSC_DESTROY_PARAM(view));
443 }
444 
446 {
447 #ifdef __INTEL_COMPILER
448  //Old versions of the intel compiler can result in a PETSC ERROR and the program aborting if parmetis is checked for before
449  //some other calls to Petsc are made. This nasty hack ensures that the HasParMetis method always works on Intel
450  Mat mat;
451  PetscTools::SetupMat(mat, 1, 1, 1);
452  PetscTools::Destroy(mat);
453 #endif
454 
455  MatPartitioning part;
456  MatPartitioningCreate(PETSC_COMM_WORLD, &part);
457 
458 
459  // We are expecting an error from PETSC on systems that don't have the interface, so suppress it
460  // in case it aborts
461  PetscPushErrorHandler(PetscIgnoreErrorHandler, NULL);
462 
463 #if (PETSC_VERSION_MAJOR == 2 || (PETSC_VERSION_MAJOR == 3 && PETSC_VERSION_MINOR < 2))
464  PetscErrorCode parmetis_installed_error = MatPartitioningSetType(part,MAT_PARTITIONING_PARMETIS);
465 #else
466  PetscErrorCode parmetis_installed_error = MatPartitioningSetType(part,MATPARTITIONINGPARMETIS);
467 #endif
468 
469  // Stop supressing error
470  PetscPopErrorHandler();
471 
472  // Note that this method probably leaks memory inside PETSc because if MatPartitioningCreate fails
473  // then there isn't a proper handle to destroy.
474  MatPartitioningDestroy(PETSC_DESTROY_PARAM(part));
475 
476  // Get out of jail free card for Windows where the latest configuration of the integration machine shows that our implementation doesn't work as expected.
477 #ifdef _MSC_VER
478 //\todo #2016 (or similar ticket). The method NodePartitioner::PetscMatrixPartitioning is not working in parallel
479  if (parmetis_installed_error == 0)
480  {
481  WARN_ONCE_ONLY("The PETSc/parMETIS interface is correctly installed but does not yet work in Windows so matrix-based partitioning will be turned off.");
482  }
483  return false;
484 #endif
485 
486  return (parmetis_installed_error == 0);
487 }
static bool ReplicateBool(bool flag)
Definition: PetscTools.cpp:186
static void Barrier(const std::string callerId="")
Definition: PetscTools.cpp:134
static void ReadPetscObject(Mat &rMat, const std::string &rOutputFileFullPath, Vec rParallelLayout=NULL)
Definition: PetscTools.cpp:361
static void DumpPetscObject(const Mat &rMat, const std::string &rOutputFileFullPath)
Definition: PetscTools.cpp:333
#define PETSC_DESTROY_PARAM(x)
Definition: PetscTools.hpp:69
#define EXCEPTION(message)
Definition: Exception.hpp:143
static bool AmMaster()
Definition: PetscTools.cpp:120
static Vec CreateVec(int size, int localSize=PETSC_DECIDE, bool ignoreOffProcEntries=true)
Definition: PetscTools.cpp:214
static MPI_Comm GetWorld()
Definition: PetscTools.cpp:174
static unsigned mNumProcessors
Definition: PetscTools.hpp:118
static bool mPetscIsInitialised
Definition: PetscTools.hpp:115
static unsigned mRank
Definition: PetscTools.hpp:121
static bool IsSequential()
Definition: PetscTools.cpp:91
static void EndRoundRobin()
Definition: PetscTools.cpp:160
static bool HasParMetis()
Definition: PetscTools.cpp:445
static Vec CreateAndSetVec(int size, double value)
Definition: PetscTools.cpp:254
static void ReplicateException(bool flag)
Definition: PetscTools.cpp:198
static void BeginRoundRobin()
Definition: PetscTools.cpp:150
static void SetupMat(Mat &rMat, int numRows, int numColumns, unsigned rowPreallocation, int numLocalRows=PETSC_DECIDE, int numLocalColumns=PETSC_DECIDE, bool ignoreOffProcEntries=true, bool newAllocationError=true)
Definition: PetscTools.cpp:268
static void Destroy(Vec &rVec)
Definition: PetscTools.hpp:351
static bool IsParallel()
Definition: PetscTools.cpp:97
static bool IsIsolated()
Definition: PetscTools.cpp:103
static void CheckCache()
Definition: PetscTools.hpp:127
static bool IsInitialised()
Definition: PetscTools.cpp:85
static bool AmTopMost()
Definition: PetscTools.cpp:126
static void IsolateProcesses(bool isolate=true)
Definition: PetscTools.cpp:169
static unsigned GetNumProcs()
Definition: PetscTools.cpp:108
static bool mIsolateProcesses
Definition: PetscTools.hpp:124
static unsigned GetMyRank()
Definition: PetscTools.cpp:114
static void ResetCache()
Definition: PetscTools.cpp:58
static const unsigned MASTER_RANK
Definition: PetscTools.hpp:140