examples/ShellMatrix_petsc_Lindbladian
examples/ShellMatrix_petsc_Lindbladian/src/main.cpp
1#include <iostream>
2#include <mpi.h>
3#include <complex>
4
5/* slepceps includes eigen value routines */
6#include <slepc.h>
7
8/* Including Eigen */
9#include <Eigen/Dense>
10#include <Eigen/Eigenvalues>
11
12/* Time measurement */
13#include <chrono>
14
15/* Maximal system size */
16#define MaxSites 64
17
18/* Hilbert space dimension */
19#define Q 2
20
21/* ScalarType */
22#define ScalarType std::complex<double>
23
24/* Monitoring function */
25PetscErrorCode monitor(EPS eps,PetscInt its,PetscInt nconv,PetscScalar *eigr,PetscScalar *eigi,PetscReal *errest,PetscInt nest,void *mctx){
26
27 PetscMPIInt myrank;
28 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &myrank));
29
30 if (myrank == 0){
31 std::cout << "[myrank: " << myrank << "] - [monitor] - iteration = " << its << " - nconv = " << nconv << " - errors: ";
32 for(PetscInt c_for = 0; c_for < nconv+1; c_for++){
33 std::cout << errest[c_for] << " ";
34 }
35 std::cout << std::endl;
36 }
37 PetscErrorCode ierr = PETSC_SUCCESS;
38 return ierr;
39}
40
41#include "DanceQ.h"
42using namespace danceq;
43
44int main(int argc, char *argv[]) {
45
46
47
48 PetscInitialize(&argc, &argv, NULL, NULL);
49
50 /* System size */
51 PetscInt L = 6;
52
53 /* Anisotropy */
54 PetscReal Jz = 1.;
55
56 /* Number of eigenvalues */
57 PetscInt nev = 2;
58
59 /* Maximal dimension of subspace */
60 PetscInt ncv = 100;
61
62 /* Coupling strength */
63 PetscReal gamma = 1;
64
65 /* Spectral shift improves the convergence significantly */
66 PetscReal shift = 100.;
67
68 /* Tolerance for the solver */
69 PetscReal tol = 1e-8;
70
71 PetscOptionsGetInt(NULL, NULL, "-L", &L, NULL);
72 PetscOptionsGetReal(NULL, NULL, "-Jz", &Jz, NULL);
73 PetscOptionsGetInt(NULL, NULL, "-nev", &nev, NULL);
74 PetscOptionsGetInt(NULL, NULL, "-ncv", &ncv, NULL);
75 PetscOptionsGetReal(NULL, NULL, "-shift", &shift, NULL);
76 PetscOptionsGetReal(NULL, NULL, "-tol", &tol, NULL);
77 PetscOptionsGetReal(NULL, NULL, "-gamma", &gamma, NULL);
78
79 /* Number of MPI ranks */
80 PetscMPIInt world_size;
81
82 /* Rank number */
83 PetscMPIInt myrank;
84
85 PetscCall(SlepcInitialize(&argc, &argv, (char *)0, (char *)0));
86 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &world_size));
87 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &myrank));
88
89 if(myrank == 0){
90 std::cout << "[myrank: " << myrank << "] - [main] - MPI processes: world_size = " << world_size << std::endl;
91 std::cout << "[myrank: " << myrank << "] - [main] - system size: L = " << L << std::endl;
92 std::cout << "[myrank: " << myrank << "] - [main] - interaction strength: Jz = " << Jz << std::endl;
93 std::cout << "[myrank: " << myrank << "] - [main] - coupling strength: gamma = " << gamma << std::endl;
94 std::cout << "[myrank: " << myrank << "] - [main] - number of requested eigen values: nev = " << nev << std::endl;
95 std::cout << "[myrank: " << myrank << "] - [main] - maximal subspace dimension: ncv = " << ncv << std::endl;
96 std::cout << "[myrank: " << myrank << "] - [main] - spectral shift to improve convergence: shift = " << shift << std::endl;
97 std::cout << "[myrank: " << myrank << "] - [main] - tolerance: tol = " << tol << std::endl;
98 }
99
100
101 /* Operator class */
102 Lindbladian_NoU1 Lind(L);
103
104 /* Simply XXZ-Heisenberg model with periodic boundary conditions */
105 for(uint64_t i = 0; i < L; i++){
106 Lind.add_operator(.5, {i,(i+1)%L}, {"S+","S-"});
107 Lind.add_operator(.5, {i,(i+1)%L}, {"S-","S+"});
108 Lind.add_operator(Jz, {i,(i+1)%L}, {"Sz","Sz"});
109 }
110
111 /* Non-hermitian jump operators */
112 Lind.add_jump_operator(gamma, {0}, {"S+"});
113 Lind.add_jump_operator(gamma, {0}, {"S-"});
114
115
116 if(myrank == 0){
117 Lind.info();
118 }
119
120 /* A simple spectral shift improves the convergence significantly */
121 Lind.add_diagonal(std::complex<double>(shift,0.));
122
123 /* Creating the matrix shell for matrix-free for multiplication */
124 auto Lind_matrix = Lind.create_PetscShellMatrix();
125
126 /* Setting up the eigen problem solver */
127 EPS eps;
128 PetscCall(EPSCreate(PETSC_COMM_WORLD, &eps));
129 PetscCall(EPSSetOperators(eps, Lind_matrix, NULL));
130 PetscCall(EPSSetProblemType(eps, EPS_NHEP));
131 PetscCall(EPSSetWhichEigenpairs(eps, EPS_LARGEST_REAL));
132 PetscCall(EPSSetType(eps,EPSKRYLOVSCHUR));
133 PetscCall(EPSSetDimensions(eps, nev, ncv, PETSC_DEFAULT));
134
135 /* Define tolerance, last argument refers to the maximal number of iterations */
136 PetscCall(EPSSetTolerances(eps,tol,10000));
137
138 /* Setting the monitoring function */
139 PetscCall(EPSMonitorSet(eps,&monitor,NULL,NULL));
140
141 /* Solve the system */
142 auto start = std::chrono::high_resolution_clock::now();
143 PetscCall(EPSSolve(eps));
144 auto end = std::chrono::high_resolution_clock::now();
145 double duration = std::chrono::duration_cast<std::chrono::seconds>(end-start).count();
146
147
148 /* Get information from the solver */
149 PetscInt nconv;
150 PetscCall(EPSGetConverged(eps, &nconv));
151 PetscInt iter_number;
152 EPSGetIterationNumber(eps,&iter_number);
153
154
155 if(myrank == 0){
156 for(uint64_t c_for = 0UL; c_for < static_cast<uint64_t>(nconv); c_for++){
157 PetscScalar E, E_imag;
158 PetscReal err;
159 EPSGetErrorEstimate(eps,static_cast<PetscInt>(c_for),&err);
160 PetscCall(EPSGetEigenvalue(eps,static_cast<PetscInt>(c_for), &E, &E_imag));
161
162 std::cout << "[myrank: " << myrank << "] - [main] - " << c_for << "th eigenvalue = " << E - shift << " with error = " << err << std::endl;
163 }
164
165 }
166
167 if(myrank == 0){
168 std::cout << "[myrank: " << myrank << "] - [main] - MPI processes: world_size = " << world_size << std::endl;
169 std::cout << "[myrank: " << myrank << "] - [main] - system size: L = " << L << std::endl;
170 std::cout << "[myrank: " << myrank << "] - [main] - interaction strength: Jz = " << Jz << std::endl;
171 std::cout << "[myrank: " << myrank << "] - [main] - coupling strength: gamma = " << gamma << std::endl;
172 std::cout << "[myrank: " << myrank << "] - [main] - number of requested eigen values: nev = " << nev << std::endl;
173 std::cout << "[myrank: " << myrank << "] - [main] - maximal subspace dimension: ncv = " << ncv << std::endl;
174 std::cout << "[myrank: " << myrank << "] - [main] - spectral shift to improve convergence: shift = " << shift << std::endl;
175 std::cout << "[myrank: " << myrank << "] - [main] - tolerance: tol = " << tol << std::endl;
176 std::cout << "[myrank: " << myrank << "] - [main] - number of iteration = " << iter_number << std::endl;
177 std::cout << "[myrank: " << myrank << "] - [main] - time = " << duration << "s" << std::endl;
178 std::cout << "[myrank: " << myrank << "] - [main] - time per iteration = " << duration/iter_number << "s" << std::endl;
179 }
180
181
182 PetscCall(MatDestroy(&Lind_matrix));
183 PetscCall(PetscFinalize());
184
185 return 0;
186}