GCC Code Coverage Report

Directory: ./
File: src/mfn/impls/expokit/mfnexpokit.c
Date: 2025-12-10 04:20:18
Exec Total Coverage
Lines: 112 121 92.6%
Functions: 3 3 100.0%
Branches: 224 340 65.9%
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 SLEPc matrix function solver: "expokit"
12
13 Method: Arnoldi method tailored for the matrix exponential
14
15 Algorithm:
16
17 Uses Arnoldi relations to compute exp(t_step*A)*v_last for
18 several time steps.
19
20 References:
21
22 [1] R. Sidje, "Expokit: a software package for computing matrix
23 exponentials", ACM Trans. Math. Softw. 24(1):130-156, 1998.
24 */
25
26 #include <slepc/private/mfnimpl.h>
27
28 30 static PetscErrorCode MFNSetUp_Expokit(MFN mfn)
29 {
30 30 PetscInt N;
31 30 PetscBool isexp;
32
33
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 30 PetscFunctionBegin;
34
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.
 30 PetscCall(MatGetSize(mfn->A,&N,NULL));
35
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 30 if (mfn->ncv==PETSC_DETERMINE) mfn->ncv = PetscMin(30,N);
36
1/2
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
 30 if (mfn->max_it==PETSC_DETERMINE) mfn->max_it = 100;
37
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.
 30 PetscCall(MFNAllocateSolution(mfn,2));
38
39
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.
 30 PetscCall(PetscObjectTypeCompare((PetscObject)mfn->fn,FNEXP,&isexp));
40
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 30 PetscCheck(isexp,PETSC_COMM_SELF,PETSC_ERR_SUP,"This solver only supports the exponential function");
41
6/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 taken 2 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 2 times.
 6 PetscFunctionReturn(PETSC_SUCCESS);
42 }
43
44 84 static PetscErrorCode MFNSolve_Expokit(MFN mfn,Vec b,Vec x)
45 {
46 84 PetscInt mxstep,mxrej,m,mb,ld,i,j,ireject,mx,k1;
47 84 Vec v,r;
48 84 Mat H,M=NULL,K=NULL;
49 84 FN fn;
50 84 PetscScalar *Harray,*B,*F,*betaF,t,sgn,sfactor;
51 84 const PetscScalar *pK;
52 84 PetscReal anorm,avnorm,tol,err_loc,rndoff,t_out,t_new,t_now,t_step;
53 84 PetscReal xm,fact,s,p1,p2,beta,beta2,gamma,delta;
54 84 PetscBool breakdown;
55
56
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 84 PetscFunctionBegin;
57 84 m = mfn->ncv;
58 84 tol = mfn->tol;
59 84 mxstep = mfn->max_it;
60 84 mxrej = 10;
61 84 gamma = 0.9;
62 84 delta = 1.2;
63 84 mb = m;
64
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.
 84 PetscCall(FNGetScale(mfn->fn,&t,&sfactor));
65
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.
 84 PetscCall(FNDuplicate(mfn->fn,PetscObjectComm((PetscObject)mfn->fn),&fn));
66
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.
 84 PetscCall(FNSetScale(fn,1.0,1.0));
67 84 t_out = PetscAbsScalar(t);
68 84 t_now = 0.0;
69
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.
 84 PetscCall(MatNorm(mfn->A,NORM_INFINITY,&anorm));
70 84 rndoff = anorm*PETSC_MACHINE_EPSILON;
71
72 84 k1 = 2;
73 84 xm = 1.0/(PetscReal)m;
74 84 beta = mfn->bnorm;
75 84 fact = PetscPowRealInt((m+1)/2.72,m+1)*PetscSqrtReal(2.0*PETSC_PI*(m+1));
76 84 t_new = (1.0/anorm)*PetscPowReal((fact*tol)/(4.0*beta*anorm),xm);
77 84 s = PetscPowReal(10.0,PetscFloorReal(PetscLog10Real(t_new))-1);
78 84 t_new = PetscCeilReal(t_new/s)*s;
79 84 sgn = t/PetscAbsScalar(t);
80
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.
 84 PetscCall(VecCopy(b,x));
82 84 ld = m+2;
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.
 84 PetscCall(PetscCalloc2(m+1,&betaF,ld*ld,&B));
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.
 84 PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,ld,ld,NULL,&H));
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.
 84 PetscCall(MatDenseGetArray(H,&Harray));
86
87 178 while (mfn->reason == MFN_CONVERGED_ITERATING) {
88 94 mfn->its++;
89
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
 94 if (PetscIsInfOrNanReal(t_new)) t_new = PETSC_MAX_REAL;
90
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 94 t_step = PetscMin(t_out-t_now,t_new);
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.
 94 PetscCall(BVInsertVec(mfn->V,0,x));
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.
 94 PetscCall(BVScaleColumn(mfn->V,0,1.0/beta));
93
7/8
✓ Branch 0 taken 6 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 8 times.
✓ Branch 4 taken 2 times.
✓ Branch 5 taken 2 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
 94 PetscCall(BVMatArnoldi(mfn->V,mfn->transpose_solve?mfn->AT:mfn->A,H,0,&mb,&beta2,&breakdown));
94
1/2
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
 94 if (breakdown) {
95 k1 = 0;
96 t_step = t_out-t_now;
97 }
98
1/2
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
 94 if (k1!=0) {
99 94 Harray[m+1+ld*m] = 1.0;
100
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.
 94 PetscCall(BVGetColumn(mfn->V,m,&v));
101
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.
 94 PetscCall(BVGetColumn(mfn->V,m+1,&r));
102
7/8
✓ Branch 0 taken 6 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 8 times.
✓ Branch 4 taken 2 times.
✓ Branch 5 taken 2 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
 94 PetscCall(MatMult(mfn->transpose_solve?mfn->AT:mfn->A,v,r));
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.
 94 PetscCall(BVRestoreColumn(mfn->V,m,&v));
104
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.
 94 PetscCall(BVRestoreColumn(mfn->V,m+1,&r));
105
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.
 94 PetscCall(BVNormColumn(mfn->V,m+1,NORM_2,&avnorm));
106 }
107
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.
 94 PetscCall(PetscArraycpy(B,Harray,ld*ld));
108
109 ireject = 0;
110
1/2
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
 94 while (ireject <= mxrej) {
111 94 mx = mb + k1;
112
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 2862 for (i=0;i<mx;i++) {
113
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 85104 for (j=0;j<mx;j++) {
114 82336 Harray[i+j*ld] = sgn*B[i+j*ld]*t_step;
115 }
116 }
117
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.
 94 PetscCall(MFN_CreateDenseMat(mx,&M));
118
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.
 94 PetscCall(MFN_CreateDenseMat(mx,&K));
119
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.
 94 PetscCall(MatDenseGetArray(M,&F));
120
7/8
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 10 times.
✓ Branch 3 taken 8 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
✓ Branch 6 taken 2 times.
✓ Branch 7 taken 2 times.
 2862 for (j=0;j<mx;j++) PetscCall(PetscArraycpy(F+j*mx,Harray+j*ld,mx));
121
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.
 94 PetscCall(MatDenseRestoreArray(M,&F));
122
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.
 94 PetscCall(FNEvaluateFunctionMat(fn,M,K));
123
124
1/2
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
 94 if (k1==0) {
125 err_loc = tol;
126 break;
127 } else {
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.
 94 PetscCall(MatDenseGetArrayRead(K,&pK));
129 94 p1 = PetscAbsScalar(beta*pK[m]);
130 94 p2 = PetscAbsScalar(beta*pK[m+1]*avnorm);
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.
 94 PetscCall(MatDenseRestoreArrayRead(K,&pK));
132
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
 94 if (p1 > 10*p2) {
133 err_loc = p2;
134 xm = 1.0/(PetscReal)m;
135 ✗ } else if (p1 > p2) {
136 ✗ err_loc = (p1*p2)/(p1-p2);
137 ✗ xm = 1.0/(PetscReal)m;
138 } else {
139 ✗ err_loc = p1;
140 ✗ xm = 1.0/(PetscReal)(m-1);
141 }
142 }
143
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
 94 if (err_loc <= delta*t_step*tol) break;
144 else {
145 ✗ t_step = gamma*t_step*PetscPowReal(t_step*tol/err_loc,xm);
146 ✗ s = PetscPowReal(10.0,PetscFloorReal(PetscLog10Real(t_step))-1);
147 ✗ t_step = PetscCeilReal(t_step/s)*s;
148 ✗ ireject = ireject+1;
149 }
150 }
151
152 94 mx = mb + PetscMax(0,k1-1);
153
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.
 94 PetscCall(MatDenseGetArrayRead(K,&pK));
154
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 2768 for (j=0;j<mx;j++) betaF[j] = beta*pK[j];
155
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.
 94 PetscCall(MatDenseRestoreArrayRead(K,&pK));
156
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.
 94 PetscCall(BVSetActiveColumns(mfn->V,0,mx));
157
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.
 94 PetscCall(BVMultVec(mfn->V,1.0,0.0,x,betaF));
158
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.
 94 PetscCall(VecNorm(x,NORM_2,&beta));
159
160 94 t_now = t_now+t_step;
161
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 94 if (t_now>=t_out) mfn->reason = MFN_CONVERGED_TOL;
162 else {
163 10 t_new = gamma*t_step*PetscPowReal((t_step*tol)/err_loc,xm);
164 10 s = PetscPowReal(10.0,PetscFloorReal(PetscLog10Real(t_new))-1);
165 10 t_new = PetscCeilReal(t_new/s)*s;
166 }
167
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
 94 err_loc = PetscMax(err_loc,rndoff);
168
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
 94 if (mfn->its==mxstep) mfn->reason = MFN_DIVERGED_ITS;
169
7/8
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 10 times.
✓ Branch 3 taken 8 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
✓ Branch 6 taken 2 times.
✓ Branch 7 taken 2 times.
 178 PetscCall(MFNMonitor(mfn,mfn->its,err_loc));
170 }
171
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.
 84 PetscCall(VecScale(x,sfactor));
172
173
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.
 84 PetscCall(MatDestroy(&M));
174
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.
 84 PetscCall(MatDestroy(&K));
175
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.
 84 PetscCall(FNDestroy(&fn));
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.
 84 PetscCall(MatDenseRestoreArray(H,&Harray));
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.
 84 PetscCall(MatDestroy(&H));
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.
 84 PetscCall(PetscFree2(betaF,B));
179
6/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 taken 2 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 2 times.
 16 PetscFunctionReturn(PETSC_SUCCESS);
180 }
181
182 /*MC
183 MFNEXPOKIT - MFNEXPOKIT = "expokit" - An Arnoldi method tailored for the
184 matrix exponential, as implemented in the Expokit package.
185
186 Notes:
187 This solver implements the method of Expokit natively in SLEPc, rather than
188 interfacing to the original Expokit software {cite:p}`Sid98`.
189
190 The method uses Arnoldi relations to compute $\exp(\alpha A)v_\mathrm{last}$
191 for several time steps, where $\alpha$ is the time step size.
192
193 Level: beginner
194
195 .seealso: [](ch:mfn), `MFN`, `MFNType`, `MFNSetType()`
196 M*/
197 40 SLEPC_EXTERN PetscErrorCode MFNCreate_Expokit(MFN mfn)
198 {
199
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 40 PetscFunctionBegin;
200 40 mfn->ops->solve = MFNSolve_Expokit;
201 40 mfn->ops->setup = MFNSetUp_Expokit;
202
6/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 taken 2 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 2 times.
 40 PetscFunctionReturn(PETSC_SUCCESS);
203 }
204