GCC Code Coverage Report

Directory: ./
File: src/eps/tutorials/ex59.c
Date: 2026-09-29 04:16:25
Exec Total Coverage
Lines: 180 188 95.7%
Functions: 8 8 100.0%
Branches: 492 830 59.3%
Line Branch Exec Source
1 /*
2 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3 SLEPc - Scalable Library for Eigenvalue Problem Computations
4 Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
5
6 This file is part of SLEPc.
7 SLEPc is distributed under a 2-clause BSD license (see LICENSE).
8 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9 */
10
11 static char help[] = "1-D discrete Kohn-Sham model solved as a Nonlinear Eigenvalue Problem with eigenvector dependence (NEPv).\n\n"
12 "Implemented by embedding an EPS inside a nonlinear SNES loop.\n\n"
13 "The Kohn-Sham (KS) equations simplify a complex, interacting many-electron system by mapping it\n"
14 "onto an equivalent system of noninteracting particles. The goal is to find the single-particle\n"
15 "wavefunctions (eigenvectors) and orbital energies (eigenvalues) that self-consistently describe the ground state.\n\n"
16 "The nonlinear Hamiltonian matrix is defined as H(X) = L + alpha * Diag(L^-1 * rho(X)), where:\n"
17 " L = 1-D discrete Laplacian operator matrix.\n"
18 " rho = electronic density vector computed from the current block of eigenvectors X as rho = diag(X * X').\n\n"
19 "The command line options are:\n"
20 " -n <n>, where <n> = number of grid points.\n"
21 " -k <k>, where <k> = number of eigenvectors to compute.\n"
22 " -alpha <alpha>, where <alpha> = real scaling parameter controlling the nonlinearity.\n"
23 " -verbose, to print the converged Hamiltonian diagonal and its sum.\n\n";
24
25 #include <slepceps.h>
26
27 /*
28 Context structure to store the necessary data regarding the discrete Kohn-Sham (DKS) problem
29 as well as the objects needed for the nonlinear iteration.
30 */
31 typedef struct {
32 // DKS context
33
34 // Inputs
35 PetscInt n; // Spatial mesh size
36 PetscInt k; // Number of eigenvectors to compute
37 PetscReal alpha; // Parameter controlling the nonlinearity
38
39 // Constants
40 Mat L; // 1D discrete Laplacian matrix
41
42 // Working vectors and base solvers
43 KSP ksp; // Linear solver (to apply L^-1)
44 Vec rho; // Vector to store the electronic density
45 Vec z; // Intermediate vector to store the result z = L^-1 * rho
46 Vec xr; // Temporary vector to extract eigenvectors during each iteration
47
48 // Nonlinear iteration context
49 EPS eps; // Solver for each H
50 Mat H; // Hamiltonian matrix
51 BV X; // Block of eigenvectors
52 Vec *initial_space; // Array of vectors to store the EPS initial space
53 } DKSContext;
54
55 /*
56 Initialize the context for the 1D discrete Kohn-Sham problem.
57
58 This function allocates the necessary memory, builds the 1D Laplacian operator,
59 and configures the linear solver (KSP) to apply L^-1.
60
61 Arguments:
62 comm - MPI communicator
63 n - Spatial mesh size
64 k - Number of eigenvectors to compute
65 alpha - Parameter controlling the nonlinearity
66 ctx_out - Output pointer where the created context will be stored
67 */
68 40 PetscErrorCode DKSCreateContext(MPI_Comm comm,PetscInt n,PetscInt k,PetscReal alpha,DKSContext **ctx_out)
69 {
70 40 DKSContext *ctx;
71 40 PetscInt Istart,Iend,i;
72 40 PC pc;
73
74
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 40 PetscFunctionBeginUser;
75
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 40 PetscCall(PetscNew(&ctx));
76 40 ctx->n=n;
77 40 ctx->k=k;
78 40 ctx->alpha=alpha;
79
80 // 1. Create and fill the Laplacian L
81
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 40 PetscCall(MatCreate(comm,&ctx->L));
82
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 40 PetscCall(MatSetSizes(ctx->L,PETSC_DECIDE,PETSC_DECIDE,n,n));
83
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 40 PetscCall(MatSetFromOptions(ctx->L));
84
85
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 40 PetscCall(MatGetOwnershipRange(ctx->L,&Istart,&Iend));
86
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 240 for (i=Istart;i<Iend;i++) {
87
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 200 PetscCall(MatSetValue(ctx->L,i,i,2.0,INSERT_VALUES));
88
6/8
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 8 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
 200 if (i>0) PetscCall(MatSetValue(ctx->L,i,i-1,-1.0,INSERT_VALUES));
89
6/8
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 8 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
 200 if (i<n-1) PetscCall(MatSetValue(ctx->L,i,i+1,-1.0,INSERT_VALUES));
90 }
91
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 40 PetscCall(MatAssemblyBegin(ctx->L,MAT_FINAL_ASSEMBLY));
92
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 40 PetscCall(MatAssemblyEnd(ctx->L,MAT_FINAL_ASSEMBLY));
93
94 // 2. Configure the KSP
95
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 40 PetscCall(KSPCreate(comm,&ctx->ksp));
96
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 40 PetscCall(KSPSetOperators(ctx->ksp,ctx->L,ctx->L));
97
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 40 PetscCall(KSPSetType(ctx->ksp,KSPPREONLY));
98
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 40 PetscCall(KSPGetPC(ctx->ksp,&pc));
99
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 40 PetscCall(PCSetType(pc,PCCHOLESKY));
100
101
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 40 if (PetscDefined(HAVE_MUMPS)) PetscCall(PCFactorSetMatSolverType(pc,MATSOLVERMUMPS));
102
103
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 40 PetscCall(KSPSetFromOptions(ctx->ksp));
104
105 // 3. Create internal work vectors
106
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 40 PetscCall(MatCreateVecs(ctx->L,&ctx->rho,&ctx->z));
107
108 40 *ctx_out=ctx;
109
5/12
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 2 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
 40 PetscFunctionReturn(PETSC_SUCCESS);
110 }
111
112 /*
113 Configure the SLEPc objects necessary for the nonlinear iteration.
114
115 This function prepares the Hamiltonian matrix by cloning the structure of the Laplacian,
116 initializes the basis vectors (BV) block for the eigenvectors, and configures the eigenvalue
117 solver (EPS) by defining the problem type and adjusting its tolerance.
118
119 Arguments:
120 ctx - Pointer to the previously initialized DKS context
121 tol - Desired tolerance for the internal eigenvalue solver
122 */
123 40 PetscErrorCode DKSSetUpIterationWorkspace(DKSContext *ctx,PetscReal tol)
124 {
125 40 MPI_Comm comm;
126
127
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 40 PetscFunctionBeginUser;
128 40 comm = PetscObjectComm((PetscObject)ctx->L);
129 /* Prepare the H matrix */
130
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 40 PetscCall(MatDuplicate(ctx->L,MAT_DO_NOT_COPY_VALUES,&ctx->H));
131
132 /* Configure the eigenvector block X */
133
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 40 PetscCall(BVCreate(comm,&ctx->X));
134
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 40 PetscCall(BVSetSizesFromVec(ctx->X,ctx->rho,ctx->k));
135
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 40 PetscCall(BVSetFromOptions(ctx->X));
136
137 /* Configure SLEPc EPS */
138
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 40 PetscCall(EPSCreate(comm,&ctx->eps));
139
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 40 PetscCall(EPSSetProblemType(ctx->eps,EPS_HEP));
140
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 40 PetscCall(EPSSetWhichEigenpairs(ctx->eps,EPS_SMALLEST_REAL));
141
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 40 PetscCall(EPSSetDimensions(ctx->eps,ctx->k,PETSC_DECIDE,PETSC_DECIDE));
142
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 40 PetscCall(EPSSetTolerances(ctx->eps,tol,PETSC_DECIDE));
143
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 40 PetscCall(EPSSetFromOptions(ctx->eps));
144
145 /* Pre-allocate temporary vectors for the nonlinear iteration loop */
146
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 40 PetscCall(VecDuplicateVecs(ctx->rho,ctx->k,&ctx->initial_space));
147
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 40 PetscCall(MatCreateVecs(ctx->H,&ctx->xr,NULL));
148
5/12
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 2 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
 8 PetscFunctionReturn(PETSC_SUCCESS);
149 }
150
151 /*
152 Generate the initial guess X0 using the exact eigenvectors of the 1D Laplacian.
153
154 This function fills the vector block X with the initial guess, which greatly
155 improves the convergence of the nonlinear iteration loop.
156
157 Formula used: vv = [1:n]'/(n+1)*pi; X0 = sin(vv * [1:k])*sqrt(2/(n+1));
158
159 Arguments:
160 ctx - Pointer to the previously initialized DKS context
161 */
162 40 PetscErrorCode DKSGenerateInitialGuess(DKSContext *ctx)
163 {
164 40 PetscInt Istart,Iend,i,j;
165 40 PetscReal h_val,norm_factor,val;
166 40 Vec col;
167 40 PetscScalar *x_local;
168
169
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 40 PetscFunctionBeginUser;
170 // Mathematical constants of the formula
171 // X0 = sin( [1:n]' * [1:k] * pi/(n+1) ) * sqrt(2/(n+1));
172 40 h_val=PETSC_PI/(ctx->n+1.0);
173 40 norm_factor=PetscSqrtReal(2.0/(ctx->n+1.0));
174
175 // Iterate over each column (eigenvector)
176
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 160 for (j=0;j<ctx->k;j++) {
177
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 120 PetscCall(BVGetColumn(ctx->X,j,&col));
178
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 120 PetscCall(VecGetOwnershipRange(col,&Istart,&Iend));
179
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 120 PetscCall(VecGetArray(col,&x_local));
180
181
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 720 for (i=Istart;i<Iend;i++) {
182 // i and j start at 0 in C, so we add 1 for the mathematical formula
183 600 val=PetscSinReal((i+1.0)*(j+1.0)*h_val)*norm_factor;
184 600 x_local[i-Istart]=val;
185 }
186
187
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 120 PetscCall(VecRestoreArray(col,&x_local));
188
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 120 PetscCall(BVRestoreColumn(ctx->X,j,&col));
189 }
190
5/12
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 2 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
 8 PetscFunctionReturn(PETSC_SUCCESS);
191 }
192
193 /*
194 Compute the eigenvector-dependent term from the eigenvectors.
195
196 In the context of the Kohn-Sham model, this function computes the electronic
197 density rho(X) as the sum of the squares of the components of each eigenvector.
198 The formula used is: rho(X) = diag(X * X').
199
200 Arguments:
201 ctx - Pointer to the initialized DKS context
202 rho_out - Pre-created vector where the computed term will be stored
203 */
204 542 PetscErrorCode DKSComputeEigvecDependentTerm(DKSContext *ctx,Vec rho_out)
205 {
206 542 PetscInt i,j,n_loc;
207 542 Vec col;
208 542 const PetscScalar *x_local;
209 542 PetscScalar *rho_local;
210
211
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 542 PetscFunctionBeginUser;
212
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 542 PetscCall(VecSet(rho_out,0.0));
213
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 542 PetscCall(VecGetLocalSize(rho_out,&n_loc));
214
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 542 PetscCall(VecGetArray(rho_out,&rho_local));
215
216
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 2168 for (j=0;j<ctx->k;j++) {
217
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 1626 PetscCall(BVGetColumn(ctx->X,j,&col));
218
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 1626 PetscCall(VecGetArrayRead(col,&x_local));
219
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 9756 for (i=0;i<n_loc;i++) {
220 8130 rho_local[i]+=x_local[i]*PetscConj(x_local[i]);
221 }
222
223
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 1626 PetscCall(VecRestoreArrayRead(col,&x_local));
224
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 1626 PetscCall(BVRestoreColumn(ctx->X,j,&col));
225 }
226
227
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 542 PetscCall(VecRestoreArray(rho_out,&rho_local));
228
5/12
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 2 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
 106 PetscFunctionReturn(PETSC_SUCCESS);
229 }
230
231 /*
232 Build the Hamiltonian from an input eigenvector-dependent term.
233
234 This function solves the linear system L * z = rho_in using the configured
235 KSP solver. Then, it assembles the updated Hamiltonian matrix stored in the context
236 using the formula: H = L + alpha * Diag(z), where z = L^-1 * rho_in.
237
238 Arguments:
239 ctx - Pointer to the initialized DKS context (ctx->H and ctx->z are updated)
240 rho_in - Vector with the proposed input eigenvector-dependent term
241 */
242 502 PetscErrorCode DKSBuildHamiltonian(DKSContext *ctx,Vec rho_in)
243 {
244
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 502 PetscFunctionBeginUser;
245 // 1. Solve L * z = rho_in
246
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 502 PetscCall(KSPSolve(ctx->ksp,rho_in,ctx->z));
247
248 // 2. Build H = L + alpha * Diag(z)
249
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 502 PetscCall(MatCopy(ctx->L,ctx->H,SAME_NONZERO_PATTERN));
250
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 502 PetscCall(VecScale(ctx->z,ctx->alpha));
251
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 502 PetscCall(MatDiagonalSet(ctx->H,ctx->z,ADD_VALUES));
252
5/12
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 2 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
 98 PetscFunctionReturn(PETSC_SUCCESS);
253 }
254
255 /*
256 Evaluation function (Callback) for the SNES nonlinear solver.
257
258 In each SNES iteration, this function receives a proposed eigenvector-dependent
259 term (rho_in), builds the Hamiltonian, solves the eigenvalue equation, computes the
260 resulting term, and returns the residual F = rho_out - rho_in.
261
262 Arguments:
263 snes - The nonlinear solver context
264 rho_in - Input state (eigenvector-dependent term) proposed by SNES
265 F - Vector where the computed residual will be stored
266 ctx_void - Pointer to the user's DKS context (DKSContext)
267 */
268 502 PetscErrorCode DKSNonlinearIteration(SNES snes,Vec rho_in,Vec F,void *ctx_void)
269 {
270 502 DKSContext *ctx=(DKSContext*)ctx_void;
271 502 PetscInt j,nconv;
272 502 Vec col;
273 502 MPI_Comm comm;
274
275
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 502 PetscFunctionBeginUser;
276 502 comm = PetscObjectComm((PetscObject)ctx->L);
277 /* 1. Build the physics (Hamiltonian) with the guess proposed by SNES */
278
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 502 PetscCall(DKSBuildHamiltonian(ctx,rho_in));
279
280 /* 2. Solve with the current H */
281 // Pass our newly built H matrix to SLEPc
282
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 502 PetscCall(EPSSetOperators(ctx->eps,ctx->H,NULL));
283
284 // Extract the eigenvectors from the previous iteration to use as the initial space
285
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 2008 for (j=0; j<ctx->k; j++) {
286
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 1506 PetscCall(BVGetColumn(ctx->X,j,&col));
287
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 1506 PetscCall(VecCopy(col,ctx->initial_space[j]));
288
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 1506 PetscCall(BVRestoreColumn(ctx->X,j,&col));
289 }
290
291 // Inject the initial space into the EPS
292
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 502 PetscCall(EPSSetInitialSpace(ctx->eps,ctx->k,ctx->initial_space));
293
294
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 502 PetscCall(EPSSolve(ctx->eps));
295
296 // (Safety check: verify that SLEPc has not failed internally)
297
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 502 PetscCall(EPSGetConverged(ctx->eps,&nconv));
298
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 502 PetscCheck(nconv>=ctx->k,comm,PETSC_ERR_NOT_CONVERGED,"SLEPc only converged %" PetscInt_FMT " out of %" PetscInt_FMT " eigenvalues in this SNES iteration",nconv,ctx->k);
299
300 /* 3. Extract the new eigenvectors and store them in ctx->X */
301
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 2008 for (j=0;j<ctx->k;j++) {
302
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 1506 PetscCall(EPSGetEigenvector(ctx->eps,j,ctx->xr,NULL)); // Extracts the j-th vector
303
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 1506 PetscCall(BVGetColumn(ctx->X,j,&col)); // Gets column j of X
304
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 1506 PetscCall(VecCopy(ctx->xr,col)); // Copies the data
305
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 1506 PetscCall(BVRestoreColumn(ctx->X,j,&col)); // Restores the column
306 }
307
308 /* 4. Compute the new eigenvector-dependent term */
309 // We use ctx->rho as a temporary work vector to store rho_out
310
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 502 PetscCall(DKSComputeEigvecDependentTerm(ctx,ctx->rho));
311
312 /* 5. Calculate the nonlinear residual: F = rho_in - rho_out */
313
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 502 PetscCall(VecWAXPY(F,-1.0,ctx->rho,rho_in));
314
5/12
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 2 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
 98 PetscFunctionReturn(PETSC_SUCCESS);
315 }
316
317 /*
318 Free the memory associated with the DKS context.
319
320 This function destroys all internal PETSc objects created during
321 initialization and frees the memory of the main structure.
322
323 Arguments:
324 ctx - Pointer to the DKS context (set to NULL upon completion)
325 */
326 40 PetscErrorCode DKSDestroyContext(DKSContext **ctx)
327 {
328
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 40 PetscFunctionBeginUser;
329
2/14
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
 40 if (!*ctx) PetscFunctionReturn(PETSC_SUCCESS);
330
331
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 40 PetscCall(MatDestroy(&(*ctx)->L));
332
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 40 PetscCall(KSPDestroy(&(*ctx)->ksp));
333
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 40 PetscCall(VecDestroy(&(*ctx)->rho));
334
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 40 PetscCall(VecDestroy(&(*ctx)->z));
335
336
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 40 PetscCall(MatDestroy(&(*ctx)->H));
337
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 40 PetscCall(BVDestroy(&(*ctx)->X));
338
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 40 PetscCall(EPSDestroy(&(*ctx)->eps));
339
340
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 40 PetscCall(VecDestroy(&(*ctx)->xr));
341
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 40 PetscCall(VecDestroyVecs((*ctx)->k, &(*ctx)->initial_space));
342
343
5/8
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 8 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
 40 PetscCall(PetscFree(*ctx));
344
5/12
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 2 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
 8 PetscFunctionReturn(PETSC_SUCCESS);
345 }
346
347 40 int main(int argc,char **argv)
348 {
349 40 DKSContext *ctx;
350 40 Vec H_diag,rho_guess;
351 40 PetscScalar diagonal_sum,kr;
352 40 PetscInt n=5,k=3,maxit=100,its,i;
353 40 PetscInt nconv;
354 40 PetscReal alpha=0.5,rtol=SLEPC_DEFAULT_TOL,stol=1e-12;
355 40 SNES snes;
356 40 Vec F; // Vector to store the residual (F = rho_out - rho_in)
357 40 PetscBool verbose=PETSC_FALSE;
358 40 SNESConvergedReason reason;
359 40 SNESLineSearch linesearch;
360
361
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 40 PetscFunctionBeginUser;
362
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 40 PetscCall(SlepcInitialize(&argc,&argv,NULL,help));
363
364
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 40 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"--- DKS SNES NEPv ---\n"));
365
366 /* 1. Read options from the command line (if provided by the user) */
367
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 40 PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
368
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 40 PetscCall(PetscOptionsGetInt(NULL,NULL,"-k",&k,NULL));
369
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 40 PetscCall(PetscOptionsGetReal(NULL,NULL,"-alpha",&alpha,NULL));
370
371 /* 2. Initialization of the DKS context */
372
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 40 PetscCall(DKSCreateContext(PETSC_COMM_WORLD,n,k,alpha,&ctx));
373
374 /* 3. Configure general objects for the nonlinear iteration (Matrices, EPS, eigenvectors) */
375
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 40 PetscCall(DKSSetUpIterationWorkspace(ctx,rtol*0.1));
376
377 /* 4. Generate the initial guess for the eigenvectors and their density */
378
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 40 PetscCall(DKSGenerateInitialGuess(ctx));
379
380 /* 5. Calculate the initial density (rho_guess) from X0 */
381 // Since SNES works with densities, the density is extracted from our X0
382
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 40 PetscCall(VecDuplicate(ctx->rho,&rho_guess));
383
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 40 PetscCall(DKSComputeEigvecDependentTerm(ctx,rho_guess));
384
385 /* 6. Prepare the residual vector F by cloning the structure of rho */
386
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 40 PetscCall(VecDuplicate(ctx->rho,&F));
387
388 /* 7. Create and set up the nonlinear solver engine (SNES) */
389
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 40 PetscCall(SNESCreate(PETSC_COMM_WORLD,&snes));
390
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 40 PetscCall(SNESSetFunction(snes,F,DKSNonlinearIteration,ctx));
391
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 40 PetscCall(SNESSetTolerances(snes,PETSC_DETERMINE,rtol,stol,maxit,PETSC_DETERMINE));
392
393 // Default -> NRICHARDSON with step lambda = 1.0 to simulate basic SCF
394
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 40 PetscCall(SNESSetType(snes,SNESNRICHARDSON));
395
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 40 PetscCall(SNESGetLineSearch(snes,&linesearch));
396
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 40 PetscCall(SNESLineSearchSetType(linesearch,SNESLINESEARCHNONE));
397
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 40 PetscCall(SNESLineSearchSetDamping(linesearch,1.0));
398
399
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 40 PetscCall(SNESSetFromOptions(snes));
400
401
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 40 PetscCall(SNESGetTolerances(snes,NULL,&rtol,NULL,&maxit,NULL));
402
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 40 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Current parameters: n=%" PetscInt_FMT ", k=%" PetscInt_FMT ", alpha=%g, rtol=%g, maxit=%" PetscInt_FMT "\n",n,k,(double)alpha,(double)rtol,maxit));
403
404 /* 8. Nonlinear iteration loop */
405
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 40 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nStarting nonlinear iterations...\n"));
406
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 40 PetscCall(SNESSolve(snes,NULL,rho_guess));
407
408
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 40 PetscCall(SNESGetConvergedReason(snes,&reason));
409
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 40 PetscCall(SNESGetIterationNumber(snes,&its));
410
411 /* 9. Result analysis */
412
1/2
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
 40 if (reason>0) {
413
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 40 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"--> CONVERGENCE REACHED in %" PetscInt_FMT " iterations (Reason: %s).\n",its,SNESConvergedReasons[reason]));
414
415 /* Show the final result */
416
417
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 40 PetscCall(EPSGetConverged(ctx->eps,&nconv));
418
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 40 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n--- Eigenvalues (%" PetscInt_FMT " found) ---\n",nconv));
419
420
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 240 for (i=0;i<nconv;i++) {
421
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 200 PetscCall(EPSGetEigenvalue(ctx->eps,i,&kr,NULL));
422
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 200 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Eigenvalue[%" PetscInt_FMT "] = %10.6f\n",i,(double)PetscRealPart(kr)));
423 }
424
425
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 40 PetscCall(PetscOptionsHasName(NULL,NULL,"-verbose",&verbose));
426
427
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
 40 if (verbose) {
428 ✗ PetscCall(MatCreateVecs(ctx->H,NULL,&H_diag));
429 ✗ PetscCall(MatGetDiagonal(ctx->H,H_diag));
430
431 ✗ PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Diagonal of the converged Hamiltonian:\n"));
432 ✗ PetscCall(VecView(H_diag,PETSC_VIEWER_STDOUT_WORLD));
433
434 ✗ PetscCall(VecSum(H_diag,&diagonal_sum));
435 ✗ PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Sum of the diagonal: %g\n",(double)PetscRealPart(diagonal_sum)));
436 ✗ PetscCall(VecDestroy(&H_diag));
437 }
438
439 } else {
440 ✗ PetscCall(PetscPrintf(PETSC_COMM_WORLD,"--> ERROR: SNES did not converge (Reason: %s)\n",SNESConvergedReasons[reason]));
441 }
442
443 /* 10. Clean up memory */
444
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 40 PetscCall(VecDestroy(&F));
445
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 40 PetscCall(VecDestroy(&rho_guess));
446
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 40 PetscCall(SNESDestroy(&snes));
447
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 40 PetscCall(DKSDestroyContext(&ctx));
448
449
2/6
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 40 PetscCall(SlepcFinalize());
450 return 0;
451 }
452
453 /*TEST
454
455 testset:
456 filter: sed -e "s/1.364212/1.364211/" -e "s/4.864809/4.864808/" -e "s/1.92185[24]/1.921853/" -e "s/2.931104/2.931103/" -e "s/3.95751[46]/3.957515/" -e "s/rtol=1e-05/rtol=1e-08/" -e "s/rtol=1e-16/rtol=1e-08/" -e "s/[0-9]\{1,\} iterations/8 iterations/" -e "s/CONVERGED_SNORM_RELATIVE/CONVERGED_FNORM_RELATIVE/"
457 output_file: output/ex59_1.out
458 test:
459 suffix: 1
460 test:
461 suffix: 2
462 args: -snes_type anderson -snes_anderson_m 7
463 test:
464 suffix: 3
465 args: -snes_type ngmres -npc_snes_type nrichardson -snes_npc_side right
466 test:
467 suffix: 4
468 args: -snes_type composite -snes_composite_type additiveoptimal -snes_composite_sneses anderson,nrichardson -sub_0_snes_anderson_m 7 -sub_0_snes_anderson_beta 0.4
469
470 TEST*/
471