GCC Code Coverage Report

Directory: ./
File: src/eps/tutorials/ex59.c
Date: 2026-07-29 03:58:07
Exec Total Coverage
Lines: 173 182 95.1%
Functions: 8 8 100.0%
Branches: 466 802 58.1%
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 self-consistent field (SCF) iteration.
30 */
31 typedef struct {
32 // DKS context
33 PetscInt n; // Spatial mesh size
34 PetscInt k; // Number of eigenvectors to compute
35 PetscReal alpha; // Parameter controlling the nonlinearity
36 Mat L; // 1D discrete Laplacian matrix
37 KSP ksp; // Linear solver (to apply L^-1)
38
39 Vec rho; // Vector to store the electronic density
40 Vec z; // Intermediate vector to store the result z = L^-1 * rho
41
42 // SCF context
43 EPS eps; // Solver for each H
44 Mat H; // Hamiltonian matrix
45 BV X; // Block of eigenvectors
46 } DKSContext;
47
48 /*
49 Initialize the context for the 1D discrete Kohn-Sham problem.
50
51 This function allocates the necessary memory, builds the 1D Laplacian operator,
52 and configures the linear solver (KSP) to apply L^-1.
53
54 Arguments:
55 comm - MPI communicator
56 n - Spatial mesh size
57 k - Number of eigenvectors to compute
58 alpha - Parameter controlling the nonlinearity
59 ctx_out - Output pointer where the created context will be stored
60 */
61 40 PetscErrorCode DKSCreate(MPI_Comm comm,PetscInt n,PetscInt k,PetscReal alpha,DKSContext **ctx_out)
62 {
63 40 DKSContext *ctx;
64 40 PetscInt Istart,Iend,i;
65 40 PC pc;
66
67
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 40 PetscFunctionBeginUser;
68
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));
69 40 ctx->n=n;
70 40 ctx->k=k;
71 40 ctx->alpha=alpha;
72
73 // 1. Create and fill the Laplacian L
74
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));
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(MatSetSizes(ctx->L,PETSC_DECIDE,PETSC_DECIDE,n,n));
76
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));
77
78
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));
79
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 240 for (i=Istart;i<Iend;i++) {
80
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));
81
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));
82
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));
83 }
84
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));
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(MatAssemblyEnd(ctx->L,MAT_FINAL_ASSEMBLY));
86
87 // 2. Configure the KSP
88
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));
89
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));
90
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));
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(KSPGetPC(ctx->ksp,&pc));
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(PCSetType(pc,PCCHOLESKY));
93
94
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 40 if (PetscDefined(HAVE_MUMPS)) PetscCall(PCFactorSetMatSolverType(pc,MATSOLVERMUMPS));
95
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(KSPSetFromOptions(ctx->ksp));
97
98 // 3. Create internal work vectors
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(MatCreateVecs(ctx->L,&ctx->rho,&ctx->z));
100
101 40 *ctx_out=ctx;
102
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);
103 }
104
105 /*
106 Configure the SLEPc objects necessary for the SCF iteration.
107
108 This function prepares the Hamiltonian matrix by cloning the structure of the Laplacian,
109 initializes the basis vectors (BV) block for the eigenvectors, and configures the eigenvalue
110 solver (EPS) by defining the problem type and adjusting its tolerance.
111
112 Arguments:
113 ctx - Pointer to the previously initialized DKS context
114 tol - Desired tolerance for the internal eigenvalue solver
115 */
116 40 PetscErrorCode DKSSetupSCF(DKSContext *ctx,PetscReal tol)
117 {
118 40 MPI_Comm comm;
119
120
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 40 PetscFunctionBeginUser;
121 40 comm = PetscObjectComm((PetscObject)ctx->L);
122 /* Prepare the H matrix */
123
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));
124
125 /* Configure the eigenvector block X */
126
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));
127
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));
128
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));
129
130 /* Configure SLEPc EPS */
131
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));
132
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));
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(EPSSetWhichEigenpairs(ctx->eps,EPS_SMALLEST_REAL));
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(EPSSetDimensions(ctx->eps,ctx->k,PETSC_DECIDE,PETSC_DECIDE));
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(EPSSetTolerances(ctx->eps,tol,PETSC_DECIDE));
136
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));
137
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);
138 }
139
140 /*
141 Generate the initial guess X0 using the exact eigenvectors of the 1D Laplacian.
142
143 This function fills the vector block X with the initial guess, which greatly
144 improves the convergence of the SCF loop.
145
146 Formula used: vv = [1:n]'/(n+1)*pi; X0 = sin(vv * [1:k])*sqrt(2/(n+1));
147
148 Arguments:
149 ctx - Pointer to the previously initialized DKS context
150 */
151 40 PetscErrorCode DKSGenerateInitialGuess(DKSContext *ctx)
152 {
153 40 PetscInt Istart,Iend,i,j;
154 40 PetscReal h_val,norm_factor,val;
155 40 Vec col;
156 40 PetscScalar *x_local;
157
158
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 40 PetscFunctionBeginUser;
159 // Mathematical constants of the formula
160 // X0 = sin( [1:n]' * [1:k] * pi/(n+1) ) * sqrt(2/(n+1));
161 40 h_val=PETSC_PI/(ctx->n+1.0);
162 40 norm_factor=PetscSqrtReal(2.0/(ctx->n+1.0));
163
164 // Iterate over each column (eigenvector)
165
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 160 for (j=0;j<ctx->k;j++) {
166
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));
167
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));
168
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));
169
170
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 720 for (i=Istart;i<Iend;i++) {
171 // i and j start at 0 in C, so we add 1 for the mathematical formula
172 600 val=PetscSinReal((i+1.0)*(j+1.0)*h_val)*norm_factor;
173 600 x_local[i-Istart]=val;
174 }
175
176
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));
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(BVRestoreColumn(ctx->X,j,&col));
178 }
179
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);
180 }
181
182 /*
183 Calculate the electronic density from the eigenvectors.
184
185 This function computes the density rho(X) as the sum of the squares of the
186 components of each eigenvector. The formula used is: rho(X) = diag(X * X').
187
188 Arguments:
189 ctx - Pointer to the initialized DKS context
190 rho_out - Pre-created vector where the computed density will be stored
191 */
192 530 PetscErrorCode DKSCalculateDensity(DKSContext *ctx,Vec rho_out)
193 {
194 530 PetscInt i,j,n_loc;
195 530 Vec col;
196 530 const PetscScalar *x_local;
197 530 PetscScalar *rho_local;
198
199
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 530 PetscFunctionBeginUser;
200
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.
 530 PetscCall(VecSet(rho_out,0.0));
201
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.
 530 PetscCall(VecGetLocalSize(rho_out,&n_loc));
202
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.
 530 PetscCall(VecGetArray(rho_out,&rho_local));
203
204
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 2120 for (j=0;j<ctx->k;j++) {
205
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.
 1590 PetscCall(BVGetColumn(ctx->X,j,&col));
206
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.
 1590 PetscCall(VecGetArrayRead(col,&x_local));
207
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 9540 for (i=0;i<n_loc;i++) {
208 7950 rho_local[i]+=x_local[i]*PetscConj(x_local[i]);
209 }
210
211
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.
 1590 PetscCall(VecRestoreArrayRead(col,&x_local));
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.
 1590 PetscCall(BVRestoreColumn(ctx->X,j,&col));
213 }
214
215
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.
 530 PetscCall(VecRestoreArray(rho_out,&rho_local));
216
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);
217 }
218
219 /*
220 Build the Hamiltonian from an input density.
221
222 This function solves the linear system L * z = rho_in using the configured
223 KSP solver. Then, it assembles the updated Hamiltonian matrix stored in the context
224 using the formula: H = L + alpha * Diag(z), where z = L^-1 * rho_in.
225
226 Arguments:
227 ctx - Pointer to the initialized DKS context (ctx->H and ctx->z are updated)
228 rho_in - Vector with the proposed input electronic density
229 */
230 490 PetscErrorCode DKSBuildHamiltonian(DKSContext *ctx,Vec rho_in)
231 {
232
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 490 PetscFunctionBeginUser;
233 // 1. Solve L * z = rho_in
234
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.
 490 PetscCall(KSPSolve(ctx->ksp,rho_in,ctx->z));
235
236 // 2. Build H = L + alpha * Diag(z)
237
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.
 490 PetscCall(MatCopy(ctx->L,ctx->H,SAME_NONZERO_PATTERN));
238
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.
 490 PetscCall(VecScale(ctx->z,ctx->alpha));
239
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.
 490 PetscCall(MatDiagonalSet(ctx->H,ctx->z,ADD_VALUES));
240
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);
241 }
242
243 /*
244 Evaluation function (Callback) for the SNES nonlinear solver.
245
246 In each SNES iteration, this function receives a proposed density (rho_in),
247 builds the Hamiltonian, solves the eigenvalue equation, computes the
248 resulting density, and returns the residual F = rho_out - rho_in.
249
250 Arguments:
251 snes - The nonlinear solver context
252 rho_in - Input electronic density proposed by SNES
253 F - Vector where the computed residual will be stored
254 ctx_void - Pointer to the user's DKS context (DKSContext)
255 */
256 490 PetscErrorCode DKSIterationSCF(SNES snes,Vec rho_in,Vec F,void *ctx_void)
257 {
258 490 DKSContext *ctx=(DKSContext*)ctx_void;
259 490 PetscInt j,nconv;
260 490 Vec xr,col;
261 490 MPI_Comm comm;
262
263
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 490 PetscFunctionBeginUser;
264 490 comm = PetscObjectComm((PetscObject)ctx->L);
265 /* 1. Build the physics (Hamiltonian) with the guess proposed by SNES */
266
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.
 490 PetscCall(DKSBuildHamiltonian(ctx,rho_in));
267
268 /* 2. Solve with the current H */
269 // Pass our newly built H matrix to SLEPc
270
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.
 490 PetscCall(EPSSetOperators(ctx->eps,ctx->H,NULL));
271
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.
 490 PetscCall(EPSSolve(ctx->eps));
272
273 // (Safety check: verify that SLEPc has not failed internally)
274
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.
 490 PetscCall(EPSGetConverged(ctx->eps,&nconv));
275
1/8
✗ Branch 0 not taken.
✓ Branch 1 taken 10 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.
 490 if (nconv<ctx->k) PetscCall(PetscPrintf(comm,"Warning: SLEPc only converged %" PetscInt_FMT " out of %" PetscInt_FMT " eigenvalues.\n",nconv,ctx->k));
276
277 /* 3. Extract the new eigenvectors and store them in ctx->X */
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.
 490 PetscCall(MatCreateVecs(ctx->H,&xr,NULL));
279
280
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 1960 for (j=0;j<ctx->k;j++) {
281
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.
 1470 PetscCall(EPSGetEigenvector(ctx->eps,j,xr,NULL)); // Extracts the j-th vector
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.
 1470 PetscCall(BVGetColumn(ctx->X,j,&col)); // Gets column j of X
283
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.
 1470 PetscCall(VecCopy(xr,col)); // Copies the data
284
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.
 1470 PetscCall(BVRestoreColumn(ctx->X,j,&col)); // Restores the column
285 }
286
287 /* 4. Calculate the new density generated by these electrons */
288 // We use ctx->rho as a temporary work vector to store rho_out
289
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.
 490 PetscCall(DKSCalculateDensity(ctx,ctx->rho));
290
291 /* 5. Calculate the nonlinear residual: F = rho_out - rho_in */
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.
 490 PetscCall(VecWAXPY(F,-1.0,ctx->rho,rho_in));
293
294 /* Cleanup of the temporary memory of this iteration */
295
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.
 490 PetscCall(VecDestroy(&xr));
296
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);
297 }
298
299 /*
300 Free the memory associated with the DKS context.
301
302 This function destroys all internal PETSc objects created during
303 initialization and frees the memory of the main structure.
304
305 Arguments:
306 ctx - Pointer to the DKS context (set to NULL upon completion)
307 */
308 40 PetscErrorCode DKSDestroy(DKSContext **ctx)
309 {
310
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 40 PetscFunctionBeginUser;
311
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);
312
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.
 40 PetscCall(MatDestroy(&(*ctx)->L));
314
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));
315
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));
316
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));
317
318
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));
319
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));
320
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));
321
322
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));
323
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);
324 }
325
326 40 int main(int argc,char **argv)
327 {
328 40 DKSContext *ctx;
329 40 Vec H_diag,rho_guess;
330 40 PetscScalar diagonal_sum,kr;
331 40 PetscInt n=5,k=3,maxit=100,its;
332 40 PetscInt nconv;
333 40 PetscReal alpha=0.5,rtol=SLEPC_DEFAULT_TOL,stol=1e-12;
334 40 SNES snes;
335 40 Vec F; // Vector to store the residual (F = rho_out - rho_in)
336 40 PetscBool verbose=PETSC_FALSE;
337 40 SNESConvergedReason reason;
338 40 SNESLineSearch linesearch;
339
340
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 40 PetscFunctionBeginUser;
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(SlepcInitialize(&argc,&argv,NULL,help));
342
343
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 SCF ---\n"));
344
345 /* 1. Read options from the command line (if provided by the user) */
346
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));
347
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));
348
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));
349
350 /* 2. Initialization of the DKS context */
351
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(DKSCreate(PETSC_COMM_WORLD,n,k,alpha,&ctx));
352
353 /* 3. Configure general objects for the SCF (Matrices, EPS, eigenvectors) */
354
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(DKSSetupSCF(ctx,rtol*0.1));
355
356 /* 4. Generate the initial guess for the eigenvectors and their density */
357
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));
358
359 /* 5. Calculate the initial density (rho_guess) from X0 */
360 // Since SNES works with densities, the density is extracted from our X0
361
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));
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(DKSCalculateDensity(ctx,rho_guess));
363
364 /* 6. Prepare the residual vector F by cloning the structure of rho */
365
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));
366
367 /* 7. Create and set up the nonlinear solver engine (SNES) */
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(SNESCreate(PETSC_COMM_WORLD,&snes));
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(SNESSetFunction(snes,F,DKSIterationSCF,ctx));
370
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));
371
372 // Default -> NRICHARDSON with step lambda = 1.0 to simulate basic SCF
373
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));
374
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));
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(SNESLineSearchSetType(linesearch,SNESLINESEARCHNONE));
376
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));
377
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(SNESSetFromOptions(snes));
379
380
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));
381
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));
382
383 /* 8. SCF Loop */
384
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 SCF iterations...\n"));
385
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));
386
387
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));
388
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));
389
390 /* 9. Result analysis */
391
1/2
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
 40 if (reason>0) {
392
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]));
393
394 /* Show the final result */
395
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(EPSGetConverged(ctx->eps,&nconv));
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(PetscPrintf(PETSC_COMM_WORLD,"\n--- Eigenvalues (%" PetscInt_FMT " found) ---\n",nconv));
398
399
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 240 for (PetscInt i=0;i<nconv;i++) {
400
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));
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.
 200 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Eigenvalue[%" PetscInt_FMT "] = %10.6f\n",i,(double)PetscRealPart(kr)));
402 }
403
404
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));
405
406
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
 40 if (verbose) {
407 ✗ PetscCall(MatCreateVecs(ctx->H,NULL,&H_diag));
408 ✗ PetscCall(MatGetDiagonal(ctx->H,H_diag));
409
410 ✗ PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Diagonal of the converged Hamiltonian:\n"));
411 ✗ PetscCall(VecView(H_diag,PETSC_VIEWER_STDOUT_WORLD));
412
413 ✗ PetscCall(VecSum(H_diag,&diagonal_sum));
414 ✗ PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Sum of the diagonal: %g\n",(double)PetscRealPart(diagonal_sum)));
415 ✗ PetscCall(VecDestroy(&H_diag));
416 }
417
418 } else {
419 ✗ PetscCall(SNESGetConvergedReason(snes,&reason));
420 ✗ PetscCall(PetscPrintf(PETSC_COMM_WORLD,"--> ERROR: SCF did not converge (Reason: %s)\n",SNESConvergedReasons[reason]));
421 }
422
423 /* 10. Clean up memory */
424
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));
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(VecDestroy(&rho_guess));
426
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));
427
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(DKSDestroy(&ctx));
428
429
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());
430 return 0;
431 }
432
433 /*TEST
434
435 testset:
436 filter: sed -e "s/1.364212/1.364211/" -e "s/4.864809/4.864808/" -e "s/1.921852/1.921853/" -e "s/1.921854/1.921853/" -e "s/2.931104/2.931103/" -e "s/3.957516/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/"
437 output_file: output/ex59_1.out
438 test:
439 suffix: 1
440 test:
441 suffix: 2
442 args: -snes_type anderson -snes_anderson_m 7
443 test:
444 suffix: 3
445 args: -snes_type ngmres -npc_snes_type nrichardson -snes_npc_side right
446 test:
447 suffix: 4
448 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
449
450 TEST*/
451