GCC Code Coverage Report


Directory: ./
File: src/mfn/impls/krylov/mfnkrylov.c
Date: 2025-10-03 04:28:47
Exec Total Coverage
Lines: 70 70 100.0%
Functions: 3 3 100.0%
Branches: 178 262 67.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: "krylov"
12
13 Method: Arnoldi with Eiermann-Ernst restart
14
15 Algorithm:
16
17 Build Arnoldi approximations using f(H) for the Hessenberg matrix H,
18 restart by discarding the Krylov basis but keeping H.
19
20 References:
21
22 [1] M. Eiermann and O. Ernst, "A restarted Krylov subspace method
23 for the evaluation of matrix functions", SIAM J. Numer. Anal.
24 44(6):2481-2504, 2006.
25 */
26
27 #include <slepc/private/mfnimpl.h>
28 #include <slepcblaslapack.h>
29
30 160 static PetscErrorCode MFNSetUp_Krylov(MFN mfn)
31 {
32 160 PetscInt N;
33
34
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
160 PetscFunctionBegin;
35
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.
160 PetscCall(MatGetSize(mfn->A,&N,NULL));
36
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
160 if (mfn->ncv==PETSC_DETERMINE) mfn->ncv = PetscMin(30,N);
37
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
160 if (mfn->max_it==PETSC_DETERMINE) mfn->max_it = 100;
38
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.
160 PetscCall(MFNAllocateSolution(mfn,1));
39
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.
29 PetscFunctionReturn(PETSC_SUCCESS);
40 }
41
42 964 static PetscErrorCode MFNSolve_Krylov(MFN mfn,Vec b,Vec x)
43 {
44 964 PetscInt n=0,m,ld,ldh,j;
45 964 PetscBLASInt m_,inc=1;
46 964 Mat M,G=NULL,H=NULL;
47 964 Vec F=NULL;
48 964 PetscScalar *marray,*farray,*harray;
49 964 const PetscScalar *garray;
50 964 PetscReal beta,betaold=0.0,nrm=1.0;
51 964 PetscBool breakdown;
52
53
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
964 PetscFunctionBegin;
54 964 m = mfn->ncv;
55 964 ld = m+1;
56
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.
964 PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,ld,m,NULL,&M));
57
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.
964 PetscCall(MatDenseGetArray(M,&marray));
58
59 /* set initial vector to b/||b|| */
60
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.
964 PetscCall(BVInsertVec(mfn->V,0,b));
61
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.
964 PetscCall(BVScaleColumn(mfn->V,0,1.0/mfn->bnorm));
62
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.
964 PetscCall(VecSet(x,0.0));
63
64 /* Restart loop */
65
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
3262 while (mfn->reason == MFN_CONVERGED_ITERATING) {
66 2298 mfn->its++;
67
68 /* compute Arnoldi factorization */
69
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.
2298 PetscCall(BVMatArnoldi(mfn->V,mfn->transpose_solve?mfn->AT:mfn->A,M,0,&m,&beta,&breakdown));
70
71 /* save previous Hessenberg matrix in G; allocate new storage for H and f(H) */
72
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
2298 if (mfn->its>1) { G = H; H = NULL; }
73 2298 ldh = n+m;
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.
2298 PetscCall(MFN_CreateVec(ldh,&F));
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.
2298 PetscCall(MFN_CreateDenseMat(ldh,&H));
76
77 /* glue together the previous H and the new H obtained with Arnoldi */
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.
2298 PetscCall(MatDenseGetArray(H,&harray));
79
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.
51270 for (j=0;j<m;j++) PetscCall(PetscArraycpy(harray+n+(j+n)*ldh,marray+j*ld,m));
80
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
2298 if (mfn->its>1) {
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.
1334 PetscCall(MatDenseGetArrayRead(G,&garray));
82
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.
44270 for (j=0;j<n;j++) PetscCall(PetscArraycpy(harray+j*ldh,garray+j*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.
1334 PetscCall(MatDenseRestoreArrayRead(G,&garray));
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.
1334 PetscCall(MatDestroy(&G));
85 1334 harray[n+(n-1)*ldh] = betaold;
86 }
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.
2298 PetscCall(MatDenseRestoreArray(H,&harray));
88
89
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
2298 if (mfn->its==1) {
90 /* set symmetry flag of H from A */
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.
964 PetscCall(MatPropagateSymmetryOptions(mfn->A,H));
92 }
93
94 /* evaluate f(H) */
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.
2298 PetscCall(FNEvaluateFunctionMatVec(mfn->fn,H,F));
96
97 /* x += ||b||*V*f(H)*e_1 */
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.
2298 PetscCall(VecGetArray(F,&farray));
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.
2298 PetscCall(PetscBLASIntCast(m,&m_));
100 2298 nrm = BLASnrm2_(&m_,farray+n,&inc); /* relative norm of the update ||u||/||b|| */
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.
2298 PetscCall(MFNMonitor(mfn,mfn->its,nrm));
102
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
51270 for (j=0;j<m;j++) farray[j+n] *= mfn->bnorm;
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.
2298 PetscCall(BVSetActiveColumns(mfn->V,0,m));
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.
2298 PetscCall(BVMultVec(mfn->V,1.0,1.0,x,farray+n));
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.
2298 PetscCall(VecRestoreArray(F,&farray));
106
107 /* check convergence */
108
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
2298 if (mfn->its >= mfn->max_it) mfn->reason = MFN_DIVERGED_ITS;
109
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
2298 if (mfn->its>1) {
110
5/8
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 10 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 10 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 10 times.
✓ Branch 7 taken 10 times.
1334 if (m<mfn->ncv || breakdown || beta==0.0 || nrm<mfn->tol) mfn->reason = MFN_CONVERGED_TOL;
111 }
112
113 /* restart with vector v_{m+1} */
114
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
2298 if (mfn->reason == MFN_CONVERGED_ITERATING) {
115
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.
1334 PetscCall(BVCopyColumn(mfn->V,m,0));
116 1334 n += m;
117 1334 betaold = beta;
118 }
119 }
120
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.
964 PetscCall(MatDestroy(&H));
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.
964 PetscCall(MatDestroy(&G));
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.
964 PetscCall(VecDestroy(&F));
124
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.
964 PetscCall(MatDenseRestoreArray(M,&marray));
125
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.
964 PetscCall(MatDestroy(&M));
126
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.
189 PetscFunctionReturn(PETSC_SUCCESS);
127 }
128
129 174 SLEPC_EXTERN PetscErrorCode MFNCreate_Krylov(MFN mfn)
130 {
131
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
174 PetscFunctionBegin;
132 174 mfn->ops->solve = MFNSolve_Krylov;
133 174 mfn->ops->setup = MFNSetUp_Krylov;
134
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.
174 PetscFunctionReturn(PETSC_SUCCESS);
135 }
136