GCC Code Coverage Report


Directory: ./
File: src/sys/classes/st/impls/filter/chebyshev.c
Date: 2026-02-22 03:58:10
Exec Total Coverage
Lines: 169 173 97.7%
Functions: 9 9 100.0%
Branches: 292 471 62.0%

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 A polynomial filter based on a truncated Chebyshev series.
12 */
13
14 #include <slepc/private/stimpl.h>
15 #include "filter.h"
16
17 #define DAMPING_LANCZOS_POW 2
18
19 /*
20 Gateway to Chebyshev for evaluating y=p(A)*x
21 (Clenshaw with a vector)
22 */
23 5235 static PetscErrorCode MatMult_Chebyshev(Mat A,Vec x,Vec y)
24 {
25 5235 ST st;
26 5235 ST_FILTER *ctx;
27 5235 CHEBYSHEV_CTX cheby;
28 5235 PetscInt degree,i,w0=0,w1=1;
29 5235 PetscReal s1,s2;
30
31
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
5235 PetscFunctionBegin;
32
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
5235 PetscCall(MatShellGetContext(A,&st));
33 5235 ctx = (ST_FILTER*)st->data;
34 5235 cheby = (CHEBYSHEV_CTX)ctx->data;
35 5235 degree = ctx->polyDegree;
36 5235 s1 = 4/(ctx->right - ctx->left);
37 5235 s2 = -2*(ctx->right+ctx->left)/(ctx->right-ctx->left);
38
39
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
5235 PetscCall(VecSet(st->work[w0],0.0));
40
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
5235 PetscCall(VecCopy(x,y));
41
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
5235 PetscCall(VecScale(y,cheby->damping_coeffs[degree]*cheby->coeffs[degree]));
42
43
2/2
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 8 times.
644235 for (i=0;i<degree;i++) {
44 /* swap work vectors, save a copy of current y */
45 639000 w0 = 1-w0;
46 639000 w1 = 1-w1;
47
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
639000 PetscCall(VecCopy(y,st->work[w0]));
48 /* compute next y */
49
2/2
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 8 times.
639000 if (i == degree-1) {
50 5235 s1 = s1/2;
51 5235 s2 = s2/2;
52 }
53
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
639000 PetscCall(VecAXPBYPCZ(y,cheby->damping_coeffs[degree-i-1]*cheby->coeffs[degree-i-1],-1,s2,x,st->work[w1]));
54
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
639000 PetscCall(MatMult(ctx->T,st->work[w0],st->work[w1]));
55
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
639000 PetscCall(VecAXPY(y,s1,st->work[w1]));
56 }
57
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.
1263 PetscFunctionReturn(PETSC_SUCCESS);
58 }
59
60 /*
61 Gateway to Chebyshev for evaluating C=p(A)*B
62 (Clenshaw with a matrix)
63 */
64 366 static PetscErrorCode MatMatMult_Chebyshev(Mat A,Mat B,Mat C,void *pctx)
65 {
66 366 ST st;
67 366 ST_FILTER *ctx;
68 366 CHEBYSHEV_CTX cheby;
69 366 PetscInt degree,i,m1,m2,w0=0,w1=1;
70 366 PetscReal s1,s2;
71
72
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
366 PetscFunctionBegin;
73
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
366 PetscCall(MatShellGetContext(A,&st));
74 366 ctx = (ST_FILTER*)st->data;
75 366 cheby = (CHEBYSHEV_CTX)ctx->data;
76 366 degree = ctx->polyDegree;
77 366 s1 = 4/(ctx->right - ctx->left);
78 366 s2 = -2*(ctx->right+ctx->left)/(ctx->right-ctx->left);
79
80
2/2
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 8 times.
366 if (ctx->nW) { /* check if work matrices must be resized */
81
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
356 PetscCall(MatGetSize(B,NULL,&m1));
82
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
356 PetscCall(MatGetSize(ctx->W[0],NULL,&m2));
83
2/2
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 8 times.
356 if (m1!=m2) {
84
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
41 PetscCall(MatDestroyMatrices(ctx->nW,&ctx->W));
85 41 ctx->nW = 0;
86 }
87 }
88
2/2
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 8 times.
366 if (!ctx->nW) { /* allocate work matrices */
89 51 ctx->nW = 2;
90
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
51 PetscCall(PetscMalloc1(ctx->nW,&ctx->W));
91
7/8
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 8 times.
✓ Branch 3 taken 6 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
✓ Branch 6 taken 2 times.
✓ Branch 7 taken 2 times.
153 for (i=0;i<ctx->nW;i++) PetscCall(MatDuplicate(B,MAT_DO_NOT_COPY_VALUES,&ctx->W[i]));
92 }
93
94
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
366 PetscCall(MatScale(ctx->W[w0],0));
95
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
366 PetscCall(MatCopy(B,C,SAME_NONZERO_PATTERN));
96
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
366 PetscCall(MatScale(C,cheby->damping_coeffs[degree]*cheby->coeffs[degree]));
97
98
2/2
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 8 times.
73566 for (i=0;i<degree;i++) {
99 /* swap work matrices, save a copy of current C */
100 73200 w0 = 1-w0;
101 73200 w1 = 1-w1;
102
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
73200 PetscCall(MatCopy(C,ctx->W[w0],SAME_NONZERO_PATTERN));
103
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 8 times.
73200 if (i == degree-1) {
104 366 s1 = s1/2;
105 366 s2 = s2/2;
106 }
107
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
73200 PetscCall(MatScale(C,s2));
108
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
73200 PetscCall(MatAXPY(C,-1,ctx->W[w1],SAME_NONZERO_PATTERN));
109
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
73200 PetscCall(MatAXPY(C,cheby->damping_coeffs[degree-i-1]*cheby->coeffs[degree-i-1],B,SAME_NONZERO_PATTERN));
110
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
73200 PetscCall(MatMatMult(ctx->T,ctx->W[w0],MAT_REUSE_MATRIX,PETSC_CURRENT,&ctx->W[w1]));
111
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
73200 PetscCall(MatAXPY(C,s1,ctx->W[w1],SAME_NONZERO_PATTERN));
112 }
113
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.
73 PetscFunctionReturn(PETSC_SUCCESS);
114 }
115
116 76 static PetscErrorCode Chebyshev_Compute_Coefficients(ST st)
117 {
118 76 ST_FILTER *ctx = (ST_FILTER*)st->data;
119 76 CHEBYSHEV_CTX cheby = (CHEBYSHEV_CTX)ctx->data;
120 76 PetscInt i,degree;
121 76 PetscReal s_inta,s_intb,acosa,acosb,t1,t2,s;
122
123
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
76 PetscFunctionBegin;
124 76 degree = ctx->polyDegree;
125
5/8
✓ Branch 0 taken 8 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 6 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
76 PetscCall(PetscFree(cheby->coeffs));
126
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
76 PetscCall(PetscMalloc1(degree+1,&cheby->coeffs));
127
1/2
✓ Branch 0 taken 8 times.
✗ Branch 1 not taken.
76 s_inta = PetscMax((2.0*ctx->inta - (ctx->right + ctx->left))/(ctx->right - ctx->left),-1.0);
128
1/2
✓ Branch 0 taken 8 times.
✗ Branch 1 not taken.
76 s_intb = PetscMin((2.0*ctx->intb - (ctx->right + ctx->left))/(ctx->right - ctx->left),1.0);
129 76 acosa = PetscAcosReal(s_inta);
130 76 acosb = PetscAcosReal(s_intb);
131 76 t1 = (acosa - acosb)/PETSC_PI;
132 76 t2 = t1;
133 76 cheby->coeffs[0] = t1;
134
2/2
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 8 times.
10076 for (i=0;i<degree;i++) {
135 10000 cheby->coeffs[i+1] = 2.0*(PetscSinReal((i+1)*acosa) - PetscSinReal((i+1)*acosb)) / ((i+1)*PETSC_PI);
136 10000 t1 = t1 + cheby->coeffs[i+1]*PetscCosReal((i+1)*acosa);
137 10000 t2 = t2 + cheby->coeffs[i+1]*PetscCosReal((i+1)*acosb);
138 }
139
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
76 if (s_inta == -1.0) s = 0.5/PetscAbs(t1);
140
1/8
✗ Branch 0 not taken.
✓ Branch 1 taken 8 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.
76 else if (s_intb == 1.0) s = 0.5/PetscMin(PetscAbs(t1),PetscAbs(t2));
141
4/6
✗ Branch 0 not taken.
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 8 times.
✓ Branch 4 taken 8 times.
✓ Branch 5 taken 8 times.
132 else s = 0.5/PetscMin(PetscAbs(t1),PetscAbs(t2));
142
2/2
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 8 times.
10152 for (i=0;i<degree+1;i++) cheby->coeffs[i] = s*cheby->coeffs[i];
143
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.
18 PetscFunctionReturn(PETSC_SUCCESS);
144 }
145
146 76 static PetscErrorCode Chebyshev_Compute_Damping_Coefficients(ST st)
147 {
148 76 ST_FILTER *ctx = (ST_FILTER*)st->data;
149 76 CHEBYSHEV_CTX cheby = (CHEBYSHEV_CTX)ctx->data;
150 76 PetscInt i,degree;
151 76 PetscReal invdeg;
152
153
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
76 PetscFunctionBegin;
154 76 degree = ctx->polyDegree;
155
5/8
✓ Branch 0 taken 8 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 6 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
76 PetscCall(PetscFree(cheby->damping_coeffs));
156
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
76 PetscCall(PetscMalloc1(degree+1,&cheby->damping_coeffs));
157 76 cheby->damping_coeffs[0] = 1;
158
4/5
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 8 times.
✓ Branch 3 taken 8 times.
✗ Branch 4 not taken.
76 switch (ctx->damping) {
159 case ST_FILTER_DAMPING_NONE:
160
2/2
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 8 times.
6036 for (i=1;i<=degree;i++) cheby->damping_coeffs[i] = 1;
161 break;
162 case ST_FILTER_DAMPING_LANCZOS:
163
2/2
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 8 times.
1616 for (i=1;i<=degree;i++) cheby->damping_coeffs[i] = PetscPowRealInt(PetscSinReal(i*PETSC_PI/(degree+1))/(i*PETSC_PI/(degree+1)),DAMPING_LANCZOS_POW);
164 break;
165 8 case ST_FILTER_DAMPING_JACKSON:
166 8 invdeg = 1/((PetscReal)degree+2);
167
2/2
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 8 times.
808 for (i=1;i<=degree;i++) cheby->damping_coeffs[i] = invdeg*((degree+2-i)*PetscCosReal(PETSC_PI*i*invdeg)+PetscSinReal(PETSC_PI*i*invdeg)/PetscTanReal(PETSC_PI*invdeg));
168 break;
169 case ST_FILTER_DAMPING_FEJER:
170
2/2
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 8 times.
1616 for (i=1;i<=degree;i++) cheby->damping_coeffs[i] = 1-i/(PetscReal)(degree+1);
171 break;
172 default:
173 SETERRQ(PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_WRONG,"Invalid Chebyshev filter damping type");
174 }
175
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.
18 PetscFunctionReturn(PETSC_SUCCESS);
176 }
177
178 /*
179 Computes the coefficients and damping coefficients for the Chebyshev filter
180
181 Creates the shifted (and scaled) matrix and the filter P(z).
182 G is a shell matrix whose MatMult() and MatMatMult() applies the filter.
183 */
184 76 static PetscErrorCode STComputeOperator_Filter_Chebyshev(ST st,Mat *G)
185 {
186 76 ST_FILTER *ctx = (ST_FILTER*)st->data;
187 76 PetscInt n,m,N,M;
188
189
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
76 PetscFunctionBegin;
190
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
76 PetscCall(STSetWorkVecs(st,2));
191
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
76 PetscCall(Chebyshev_Compute_Coefficients(st));
192
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
76 PetscCall(Chebyshev_Compute_Damping_Coefficients(st));
193 /* create shell matrix*/
194
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
76 PetscCall(MatDestroy(&ctx->T));
195
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
76 PetscCall(MatDuplicate(st->A[0],MAT_COPY_VALUES,&ctx->T));
196
2/2
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 8 times.
76 if (!*G) {
197
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
60 PetscCall(MatGetSize(ctx->T,&N,&M));
198
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
60 PetscCall(MatGetLocalSize(ctx->T,&n,&m));
199
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
60 PetscCall(MatCreateShell(PetscObjectComm((PetscObject)st),n,m,N,M,st,G));
200
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
60 PetscCall(MatShellSetOperation(*G,MATOP_MULT,(PetscErrorCodeFn*)MatMult_Chebyshev));
201
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
60 PetscCall(MatShellSetMatProductOperation(*G,MATPRODUCT_AB,NULL,MatMatMult_Chebyshev,NULL,MATDENSE,MATDENSE));
202
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
60 PetscCall(MatShellSetMatProductOperation(*G,MATPRODUCT_AB,NULL,MatMatMult_Chebyshev,NULL,MATDENSECUDA,MATDENSECUDA));
203
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
60 PetscCall(MatShellSetMatProductOperation(*G,MATPRODUCT_AB,NULL,MatMatMult_Chebyshev,NULL,MATDENSEHIP,MATDENSEHIP));
204 }
205
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.
18 PetscFunctionReturn(PETSC_SUCCESS);
206 }
207
208 246 static PetscErrorCode Chebyshev_Clenshaw_Damping_Real(ST st,PetscReal x,PetscReal *px)
209 {
210 246 ST_FILTER *ctx = (ST_FILTER*)st->data;
211 246 CHEBYSHEV_CTX cheby = (CHEBYSHEV_CTX)ctx->data;
212 246 PetscReal t1=0,t2=0;
213 246 PetscInt i,degree;
214
215
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
246 PetscFunctionBegin;
216 246 degree = ctx->polyDegree;
217
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
246 PetscCheck(cheby->coeffs,PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_NULL,"Chebyshev coefficients are not computed");
218
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
246 PetscCheck(cheby->damping_coeffs,PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_NULL,"Chebyshev damping coefficients are not computed");
219 246 x = 2*x;
220
2/2
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 8 times.
22046 for (i=degree;i>1;i=i-2) {
221 21800 t2 = cheby->damping_coeffs[i]*cheby->coeffs[i] + x*t1 - t2;
222 21800 t1 = cheby->damping_coeffs[i-1]*cheby->coeffs[i-1] + x*t2 - t1;
223 }
224
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 8 times.
246 if (degree % 2) {
225 t2 = cheby->damping_coeffs[1]*cheby->coeffs[1] + x*t1 - t2;
226 t1 = t2;
227 }
228 246 *px = cheby->damping_coeffs[0]*cheby->coeffs[0] + .5*x*t1 - t2;
229
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.
246 PetscFunctionReturn(PETSC_SUCCESS);
230 }
231
232 /*
233 Computes the threshold for the Chebyshev filter
234
235 The threshold is P(scale(ctx->inta)), where ctx->inta is the lower bound of the interval.
236 It should be the same value if the upper bound was used instead.
237 */
238 246 static PetscErrorCode STFilterGetThreshold_Filter_Chebyshev(ST st,PetscReal *gamma)
239 {
240 246 ST_FILTER *ctx = (ST_FILTER*)st->data;
241 246 PetscReal c,e;
242
243
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
246 PetscFunctionBegin;
244 /* scale ctx->inta */
245 246 c = (ctx->right + ctx->left) / 2.0;
246 246 e = (ctx->right - ctx->left) / 2.0;
247
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
246 PetscCall(Chebyshev_Clenshaw_Damping_Real(st,(ctx->inta-c)/e,gamma));
248
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.
52 PetscFunctionReturn(PETSC_SUCCESS);
249 }
250
251 60 static PetscErrorCode STDestroy_Filter_Chebyshev(ST st)
252 {
253 60 ST_FILTER *ctx = (ST_FILTER*)st->data;
254 60 CHEBYSHEV_CTX cheby = (CHEBYSHEV_CTX)ctx->data;
255
256
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
60 PetscFunctionBegin;
257
5/8
✓ Branch 0 taken 8 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 6 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
60 PetscCall(PetscFree(cheby->coeffs));
258
5/8
✓ Branch 0 taken 8 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 6 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
60 PetscCall(PetscFree(cheby->damping_coeffs));
259
6/8
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 6 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
60 PetscCall(PetscFree(cheby));
260
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.
14 PetscFunctionReturn(PETSC_SUCCESS);
261 }
262
263 60 PetscErrorCode STCreate_Filter_Chebyshev(ST st)
264 {
265 60 ST_FILTER *ctx = (ST_FILTER*)st->data;
266 60 CHEBYSHEV_CTX cheby;
267
268
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
60 PetscFunctionBegin;
269
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
60 PetscCall(PetscNew(&cheby));
270 60 ctx->data = (void*)cheby;
271
272 60 ctx->computeoperator = STComputeOperator_Filter_Chebyshev;
273 60 ctx->getthreshold = STFilterGetThreshold_Filter_Chebyshev;
274 60 ctx->destroy = STDestroy_Filter_Chebyshev;
275
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.
60 PetscFunctionReturn(PETSC_SUCCESS);
276 }
277