1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
|
#include <iostream>
#include <cmath>
#include <cassert>
#include <fstream>
#include <slepceps.h>
#include <petscmat.h>
int main(int argc, char ** argv)
{
SlepcInitialize(&argc, &argv, NULL, NULL);
std::string infile(argv[1]);
std::string outfile(argv[2]);
PetscMPIInt mpi_rank;
PetscMPIInt mpi_size;
MPI_Comm_rank(PETSC_COMM_WORLD, &mpi_rank);
MPI_Comm_size(PETSC_COMM_WORLD, &mpi_size);
PetscErrorCode err;
Mat L;
err = MatCreate(PETSC_COMM_WORLD, &L); CHKERRQ(err);
err = MatSetFromOptions(L); CHKERRQ(err);
PetscViewer viewer;
err = PetscViewerBinaryOpen(PETSC_COMM_WORLD, infile.c_str(), FILE_MODE_READ, &viewer); CHKERRQ(err);
err = MatLoad(L, viewer); CHKERRQ(err);
err = MatAssemblyBegin(L, MAT_FINAL_ASSEMBLY); CHKERRQ(err);
err = MatAssemblyEnd(L, MAT_FINAL_ASSEMBLY); CHKERRQ(err);
PetscInt m;
PetscInt n;
err = MatGetSize(L, &m, &n); CHKERRQ(err);
assert(m == n);
PetscPrintf(PETSC_COMM_WORLD, "Laplacian is size %d.\n", m);
EPS eps;
err = EPSCreate(PETSC_COMM_WORLD, &eps); CHKERRQ(err);
err = EPSSetOperators(eps, L, NULL); CHKERRQ(err);
err = EPSSetProblemType(eps, EPS_HEP); CHKERRQ(err);
err = EPSSetWhichEigenpairs(eps, EPS_SMALLEST_MAGNITUDE); CHKERRQ(err);
err = EPSSetFromOptions(eps); CHKERRQ(err);
PetscPrintf(PETSC_COMM_WORLD, "Solving...\n");
err = EPSSolve(eps); CHKERRQ(err);
PetscInt num_its;
err = EPSGetIterationNumber(eps, &num_its); CHKERRQ(err);
PetscPrintf(PETSC_COMM_WORLD, "Took %d iterations.\n", num_its);
PetscInt num_ev;
err = EPSGetConverged(eps, &num_ev); CHKERRQ(err);
PetscPrintf(PETSC_COMM_WORLD, "Have %d eigenpairs.\n", num_ev);
Vec vec;
err = MatCreateVecs(L, &vec, NULL); CHKERRQ(err);
PetscReal val;
PetscReal cutoff = 1e-6;
PetscInt estimated_betti = 0;
PetscInt vec_size;
err = VecGetSize(vec, &vec_size); CHKERRQ(err);
Vec vec_r0;
if (mpi_rank == 0)
{
err = VecCreateSeq(PETSC_COMM_SELF, vec_size, &vec_r0); CHKERRQ(err);
}
VecScatter scatter;
err = VecScatterCreateToZero(vec, &scatter, &vec_r0); CHKERRQ(err);
std::ofstream ofs;
if (mpi_rank == 0)
{
ofs.open(outfile, std::ios::out);
}
for (PetscInt i = 0; i < num_ev; i++)
{
err = EPSGetEigenpair(eps, i, &val, NULL, vec, NULL); CHKERRQ(err);
PetscPrintf(PETSC_COMM_WORLD, "%d: %.20f\n", i, val);
if (std::abs(val) < cutoff)
{
estimated_betti++;
err = VecScatterBegin(scatter, vec, vec_r0, INSERT_VALUES, SCATTER_FORWARD); CHKERRQ(err);
err = VecScatterEnd(scatter, vec, vec_r0, INSERT_VALUES, SCATTER_FORWARD); CHKERRQ(err);
if (mpi_rank == 0)
{
PetscReal * tmp;
err = VecGetArray(vec_r0, &tmp); CHKERRQ(err);
for (PetscInt j = 0; j < vec_size; j++)
{
ofs << tmp[j] << " ";
}
ofs << std::endl;
err = VecRestoreArray(vec_r0, &tmp); CHKERRQ(err);
}
}
}
PetscPrintf(PETSC_COMM_WORLD, "Estimated Betti number: %d.\n", estimated_betti);
if (mpi_rank == 0)
{
ofs.close();
}
err = VecScatterDestroy(&scatter); CHKERRQ(err);
err = VecDestroy(&vec); CHKERRQ(err);
err = VecDestroy(&vec_r0); CHKERRQ(err);
err = EPSDestroy(&eps); CHKERRQ(err);
err = MatDestroy(&L); CHKERRQ(err);
SlepcFinalize();
return 0;
}
|
