MiscTools.cpp
1 #include <petscsys.h>
2 #include <slepceps.h>
3 
4 #include <vector>
5 #include <string>
6 #include <iostream>
7 #include <sstream>
8 #include <fstream>
9 #include <map>
10 #include <iomanip>
11 
12 #include <sys/stat.h>
13 #include <sys/types.h>
14 #include <dirent.h>
15 #include "MiscTools.hpp"
16 
17 
18 /* Obtained from: https://gist.github.com/orlp/3551590 */
19 PETSC_EXTERN int64_t ipow(int64_t base, uint8_t exp) {
20  static const uint8_t highest_bit_set[] = {
21  0, 1, 2, 2, 3, 3, 3, 3,
22  4, 4, 4, 4, 4, 4, 4, 4,
23  5, 5, 5, 5, 5, 5, 5, 5,
24  5, 5, 5, 5, 5, 5, 5, 5,
25  6, 6, 6, 6, 6, 6, 6, 6,
26  6, 6, 6, 6, 6, 6, 6, 6,
27  6, 6, 6, 6, 6, 6, 6, 6,
28  6, 6, 6, 6, 6, 6, 6, 255, // anything past 63 is a guaranteed overflow with base > 1
29  255, 255, 255, 255, 255, 255, 255, 255,
30  255, 255, 255, 255, 255, 255, 255, 255,
31  255, 255, 255, 255, 255, 255, 255, 255,
32  255, 255, 255, 255, 255, 255, 255, 255,
33  255, 255, 255, 255, 255, 255, 255, 255,
34  255, 255, 255, 255, 255, 255, 255, 255,
35  255, 255, 255, 255, 255, 255, 255, 255,
36  255, 255, 255, 255, 255, 255, 255, 255,
37  255, 255, 255, 255, 255, 255, 255, 255,
38  255, 255, 255, 255, 255, 255, 255, 255,
39  255, 255, 255, 255, 255, 255, 255, 255,
40  255, 255, 255, 255, 255, 255, 255, 255,
41  255, 255, 255, 255, 255, 255, 255, 255,
42  255, 255, 255, 255, 255, 255, 255, 255,
43  255, 255, 255, 255, 255, 255, 255, 255,
44  255, 255, 255, 255, 255, 255, 255, 255,
45  255, 255, 255, 255, 255, 255, 255, 255,
46  255, 255, 255, 255, 255, 255, 255, 255,
47  255, 255, 255, 255, 255, 255, 255, 255,
48  255, 255, 255, 255, 255, 255, 255, 255,
49  255, 255, 255, 255, 255, 255, 255, 255,
50  255, 255, 255, 255, 255, 255, 255, 255,
51  255, 255, 255, 255, 255, 255, 255, 255,
52  255, 255, 255, 255, 255, 255, 255, 255,
53  };
54 
55  uint64_t result = 1;
56 
57  switch (highest_bit_set[exp]) {
58  case 255: // we use 255 as an overflow marker and return 0 on overflow/underflow
59  if (base == 1) {
60  return 1;
61  }
62 
63  if (base == -1) {
64  return 1 - 2 * (exp & 1);
65  }
66 
67  return 0;
68  case 6:
69  if (exp & 1) result *= base;
70  exp >>= 1;
71  base *= base;
72  case 5:
73  if (exp & 1) result *= base;
74  exp >>= 1;
75  base *= base;
76  case 4:
77  if (exp & 1) result *= base;
78  exp >>= 1;
79  base *= base;
80  case 3:
81  if (exp & 1) result *= base;
82  exp >>= 1;
83  base *= base;
84  case 2:
85  if (exp & 1) result *= base;
86  exp >>= 1;
87  base *= base;
88  case 1:
89  if (exp & 1) result *= base;
90  default:
91  return result;
92  }
93 }
94 
95 
96 PetscErrorCode PreSplitOwnership(const MPI_Comm comm, const PetscInt N, PetscInt& locrows, PetscInt& Istart)
97 {
98  PetscErrorCode ierr = 0;
99 
100 #if 0
101  /* The petsc way */
102  PetscInt Nsize = N;
103  PetscInt Lrows = PETSC_DECIDE;
104  ierr = PetscSplitOwnership(comm, &Lrows, &Nsize); CHKERRQ(ierr);
105  Istart = 0;
106  ierr = MPI_Exscan(&Lrows, &Istart, 1, MPIU_INT, MPI_SUM, comm); CHKERRQ(ierr);
107  locrows = Lrows;
108 #else
109  /* Calculate predictively */
110  PetscMPIInt nprocs,rank;
111  ierr = MPI_Comm_size(comm, &nprocs); CHKERRQ(ierr);
112  ierr = MPI_Comm_rank(comm, &rank); CHKERRQ(ierr);
113  const PetscInt remrows = N % nprocs;
114  locrows = N / nprocs + PetscInt(rank < remrows);
115  Istart = N / nprocs * rank + (rank < remrows ? rank : remrows);
116 #endif
117 
118  return ierr;
119 }
120 
121 
122 PetscErrorCode SplitOwnership(
123  const PetscMPIInt& rank,
124  const PetscMPIInt& nprocs ,
125  const PetscInt N,
126  PetscInt& locrows,
127  PetscInt& Istart)
128 {
129  const PetscInt remrows = N % nprocs;
130  locrows = N / nprocs + PetscInt(rank < remrows);
131  Istart = N / nprocs * rank + (rank < remrows ? rank : remrows);
132  return 0;
133 }
134 
135 
136 PETSC_EXTERN PetscErrorCode InitSingleSiteOperator(const MPI_Comm& comm, const PetscInt dim, Mat* mat)
137 {
138  PetscErrorCode ierr = 0;
139 
140  if(*mat) SETERRQ(comm, 1, "Matrix was previously initialized. First, destroy the matrix and set to NULL.");
141 
142  ierr = MatCreate(comm, mat); CHKERRQ(ierr);
143  ierr = MatSetSizes(*mat, PETSC_DECIDE, PETSC_DECIDE, dim, dim); CHKERRQ(ierr);
144  ierr = MatSetType(*mat, MATMPIAIJ); CHKERRQ(ierr);
145  ierr = MatSetFromOptions(*mat); CHKERRQ(ierr);
146  ierr = MatSetUp(*mat); CHKERRQ(ierr);
147 
148  ierr = MatSetOption(*mat, MAT_NO_OFF_PROC_ENTRIES , PETSC_TRUE);
149  ierr = MatSetOption(*mat, MAT_NO_OFF_PROC_ZERO_ROWS , PETSC_TRUE);
150  ierr = MatSetOption(*mat, MAT_IGNORE_OFF_PROC_ENTRIES , PETSC_TRUE);
151  ierr = MatSetOption(*mat, MAT_IGNORE_ZERO_ENTRIES , PETSC_TRUE);
152 
153  return ierr;
154 }
155 
156 
157 PETSC_EXTERN PetscErrorCode MatSetOption_MultipleMats(
158  const std::vector<Mat>& matrices,
159  const std::vector<MatOption>& options,
160  const std::vector<PetscBool>& flgs)
161 {
162  PetscErrorCode ierr = 0;
163 
164  if(flgs.size() == 1)
165  {
166  for(const Mat& mat: matrices){
167  for(const MatOption& op: options){
168  ierr = MatSetOption(mat, op, flgs[0]); CHKERRQ(ierr);
169  }
170  }
171  }
172  else if(flgs.size() == options.size())
173  {
174  for(const Mat& mat: matrices){
175  size_t i = 0;
176  for(const MatOption& op: options){
177  ierr = MatSetOption(mat, op, flgs[i++]); CHKERRQ(ierr);
178  }
179  }
180  }
181  else
182  {
183  SETERRQ(PETSC_COMM_WORLD, 1, "Incorrect input. Either flgs.size() == 1 or flgs.size() == options.size()");
184  }
185 
186  return ierr;
187 }
188 
189 
190 PETSC_EXTERN PetscErrorCode MatSetOption_MultipleMatGroups(
191  const std::vector<std::vector<Mat>>& matgroups,
192  const std::vector<MatOption>& options,
193  const std::vector<PetscBool>& flgs)
194 {
195  PetscErrorCode ierr = 0;
196 
197  for(const std::vector<Mat>& matrices: matgroups){
198  ierr = MatSetOption_MultipleMats(matrices, options, flgs); CHKERRQ(ierr);
199  }
200 
201  return ierr;
202 }
203 
204 
205 /*----- Spin-1/2 functions -----*/
206 
207 PETSC_EXTERN PetscErrorCode MatSpinOneHalfSzCreate(const MPI_Comm& comm, Mat& Sz)
208 {
209  PetscErrorCode ierr = 0;
210 
211  PetscInt loc_dim = 2;
212  ierr = InitSingleSiteOperator(comm, loc_dim, &Sz); CHKERRQ(ierr);
213 
214  PetscInt locrows = 0, Istart = 0;
215  ierr = PreSplitOwnership(comm, loc_dim, locrows, Istart); CHKERRQ(ierr);
216  PetscInt Iend = Istart + locrows;
217 
218  if (Istart <= 0 && 0 < Iend){
219  ierr = MatSetValue(Sz, 0, 0, +0.5, INSERT_VALUES); CHKERRQ(ierr);
220  }
221  if (Istart <= 1 && 1 < Iend){
222  ierr = MatSetValue(Sz, 1, 1, -0.5, INSERT_VALUES); CHKERRQ(ierr);
223  }
224 
225  ierr = MatAssemblyBegin(Sz, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
226  ierr = MatAssemblyEnd(Sz, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
227  return ierr;
228 }
229 
230 
231 PETSC_EXTERN PetscErrorCode MatSpinOneHalfSpCreate(const MPI_Comm& comm, Mat& Sp)
232 {
233  PetscErrorCode ierr = 0;
234 
235  PetscInt loc_dim = 2;
236  ierr = InitSingleSiteOperator(comm, loc_dim, &Sp); CHKERRQ(ierr);
237 
238  PetscInt locrows = 0, Istart = 0;
239  ierr = PreSplitOwnership(comm, loc_dim, locrows, Istart); CHKERRQ(ierr);
240  PetscInt Iend = Istart + locrows;
241 
242  if (Istart <= 0 && 0 < Iend){
243  ierr = MatSetValue(Sp, 0, 1, +1.0, INSERT_VALUES); CHKERRQ(ierr);
244  }
245 
246  ierr = MatAssemblyBegin(Sp, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
247  ierr = MatAssemblyEnd(Sp, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
248  return ierr;
249 }
250 
251 
252 PETSC_EXTERN PetscErrorCode MatEnsureAssembled(const Mat& matin)
253 {
254  PetscErrorCode ierr = 0;
255 
256  PetscBool assembled;
257  ierr = MatAssembled(matin, &assembled); CHKERRQ(ierr);
258  if(!assembled)
259  {
260  ierr = MatAssemblyBegin(matin, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
261  ierr = MatAssemblyEnd(matin, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
262  }
263 
264  return ierr;
265 }
266 
267 PETSC_EXTERN PetscErrorCode MatEnsureAssembled_MultipleMats(const std::vector<Mat>& matrices)
268 {
269  PetscErrorCode ierr = 0;
270 
271  for(const Mat& mat: matrices){
272  ierr = MatEnsureAssembled(mat); CHKERRQ(ierr);
273  }
274 
275  return ierr;
276 }
277 
278 
279 PETSC_EXTERN PetscErrorCode MatEnsureAssembled_MultipleMatGroups(const std::vector<std::vector<Mat>>& matgroups)
280 {
281  PetscErrorCode ierr = 0;
282 
283  for(const std::vector<Mat>& matrices: matgroups){
284  ierr = MatEnsureAssembled_MultipleMats(matrices); CHKERRQ(ierr);
285  }
286 
287  return ierr;
288 }
289 
290 
291 PETSC_EXTERN PetscErrorCode MatEyeCreate(const MPI_Comm& comm, const PetscInt& dim, Mat& eye)
292 {
293  PetscErrorCode ierr = 0;
294 
295  ierr = InitSingleSiteOperator(comm, dim, &eye); CHKERRQ(ierr);
296  ierr = MatEnsureAssembled(eye); CHKERRQ(ierr);
297  ierr = MatShift(eye, 1.00); CHKERRQ(ierr);
298  ierr = MatEnsureAssembled(eye); CHKERRQ(ierr);
299 
300  return ierr;
301 }
302 
303 
304 /*----- Utlities -----*/
305 
306 PETSC_EXTERN PetscErrorCode Makedir(const std::string& dir_name)
307 {
308  PetscErrorCode ierr;
309  DIR *dir = opendir(dir_name.c_str());
310  if(!dir){
311  /* Info on mode_t: https://jameshfisher.com/2017/02/24/what-is-mode_t.html */
312  ierr = mkdir(dir_name.c_str(), S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
313  if(ierr){
314  DIR *dir = opendir(dir_name.c_str());
315  if(!dir){
316  if(ierr) PetscPrintf(PETSC_COMM_SELF,"mkdir error code: %d for dir: %s\n",
317  ierr, dir_name.c_str());
318  CHKERRQ(ierr);
319  }
320  closedir(dir);
321  }
322  } else {
323  closedir(dir);
324  }
325  return(0);
326 }
327 
328 
329 PetscErrorCode SetOptionsFromFile(
330  MPI_Comm& mpi_comm,
331  const std::string& filename
332  )
333 {
334  PetscErrorCode ierr;
335  PetscMPIInt mpi_rank;
336  ierr = MPI_Comm_rank(mpi_comm, &mpi_rank); CHKERRQ(ierr);
337  std::map< std::string, std::string > infomap;
338  ierr = RetrieveInfoFile<std::string>(mpi_comm,filename,infomap); CHKERRQ(ierr);
339  for(auto it: infomap) {
340  ierr = PetscOptionsSetValue(NULL,(it.first).c_str(),(it.second).c_str()); CHKERRQ(ierr);
341  }
342  if(!mpi_rank) {
343  std::cout << "-----------------------------------------" << std::endl;
344  std::cout << "WARNING:\n"
345  "The following directives from " << filename << "\n" <<
346  "will override command-line arguments:" << std::endl;
347  for(auto it: infomap)
348  std::cout << " " << std::setw(20) << it.first << " " << it.second << std::endl;
349  }
350  return(0);
351 }
Makedir
PETSC_EXTERN PetscErrorCode Makedir(const std::string &dir_name)
Creates a directory named dir_name.
Definition: MiscTools.cpp:306
SetOptionsFromFile
PetscErrorCode SetOptionsFromFile(MPI_Comm &mpi_comm, const std::string &filename)
Reads a file containing keys and values and sets them as command line arguments.
Definition: MiscTools.cpp:329
This site was generated by Sphinx using Doxygen with a customized theme from doxygen-bootstrapped.