GCC Code Coverage Report

Directory: ./
File: src/sys/classes/fn/impls/log/fnlog.c
Date: 2025-10-04 04:19:13
Exec Total Coverage
Lines: 368 394 93.4%
Functions: 17 17 100.0%
Branches: 603 1046 57.6%
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 Logarithm function log(x)
12 */
13
14 #include <slepc/private/fnimpl.h> /*I "slepcfn.h" I*/
15 #include <slepcblaslapack.h>
16
17 20 static PetscErrorCode FNEvaluateFunction_Log(FN fn,PetscScalar x,PetscScalar *y)
18 {
19
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 20 PetscFunctionBegin;
20 #if !defined(PETSC_USE_COMPLEX)
21
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 10 PetscCheck(x>=0.0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Function not defined in the requested value");
22 #endif
23 20 *y = PetscLogScalar(x);
24
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.
 20 PetscFunctionReturn(PETSC_SUCCESS);
25 }
26
27 20 static PetscErrorCode FNEvaluateDerivative_Log(FN fn,PetscScalar x,PetscScalar *y)
28 {
29
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 20 PetscFunctionBegin;
30
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 20 PetscCheck(x!=0.0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Derivative not defined in the requested value");
31 #if !defined(PETSC_USE_COMPLEX)
32
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 10 PetscCheck(x>0.0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Derivative not defined in the requested value");
33 #endif
34 20 *y = 1.0/x;
35
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.
 20 PetscFunctionReturn(PETSC_SUCCESS);
36 }
37
38 /*
39 Block structure of a quasitriangular matrix T. Returns a list of n-1 numbers, where
40 structure(j) encodes the block type of the j:j+1,j:j+1 diagonal block as one of:
41 0 - not the start of a block
42 1 - start of a 2x2 triangular block
43 2 - start of a 2x2 quasi-triangular (full) block
44 */
45 50 static PetscErrorCode qtri_struct(PetscBLASInt n,PetscScalar *T,PetscBLASInt ld,PetscInt *structure)
46 {
47 50 PetscBLASInt j;
48
49
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 50 PetscFunctionBegin;
50 #if defined(PETSC_USE_COMPLEX)
51
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
 2250 for (j=0;j<n-1;j++) structure[j] = 1;
52 #else
53
2/14
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 1 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.
 20 if (n==1) PetscFunctionReturn(PETSC_SUCCESS);
54
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
 20 else if (n==2) {
55 ✗ structure[0] = (T[1]==0.0)? 1: 2;
56 ✗ PetscFunctionReturn(PETSC_SUCCESS);
57 }
58 j = 0;
59
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
 1120 while (j<n-2) {
60
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
 1100 if (T[j+1+j*ld] != 0.0) { /* Start of a 2x2 full block */
61 360 structure[j++] = 2;
62 360 structure[j++] = 0;
63
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
 740 } else if (T[j+1+j*ld] == 0.0 && T[j+2+(j+1)*ld] == 0.0) { /* Start of a 2x2 triangular block */
64 730 structure[j++] = 1;
65 } else { /* Next block must start a 2x2 full block */
66 10 structure[j++] = 0;
67 }
68 }
69
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
 20 if (T[n-1+(n-2)*ld] != 0.0) { /* 2x2 full block at the end */
70 10 structure[n-2] = 2;
71
1/2
✓ Branch 0 taken 5 times.
✗ Branch 1 not taken.
 10 } else if (structure[n-3] == 0 || structure[n-3] == 1) {
72 10 structure[n-2] = 1;
73 }
74 #endif
75
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.
 34 PetscFunctionReturn(PETSC_SUCCESS);
76 }
77
78 /*
79 Compute scaling parameter (s) and order of Pade approximant (m).
80 wr,wi is overwritten. Required workspace is 3*n*n.
81 On output, Troot contains the sth square root of T.
82 */
83 50 static PetscErrorCode FNlogm_params(FN fn,PetscBLASInt n,PetscScalar *T,PetscBLASInt ld,PetscScalar *wr,PetscScalar *wi,PetscInt maxroots,PetscInt *s,PetscInt *m,PetscScalar *Troot,PetscScalar *work)
84 {
85 50 PetscInt i,j,k,p,s0;
86 50 PetscReal inrm,eta,a2,a3,a4,d2,d3,d4,d5;
87 50 PetscScalar *TrootmI=work+2*n*n;
88 50 PetscBool foundm=PETSC_FALSE,more;
89 50 PetscRandom rand;
90 50 const PetscReal xvals[] = { 1.586970738772063e-005, 2.313807884242979e-003, 1.938179313533253e-002,
91 6.209171588994762e-002, 1.276404810806775e-001, 2.060962623452836e-001, 2.879093714241194e-001 };
92 50 const PetscInt mmax=PETSC_STATIC_ARRAY_LENGTH(xvals);
93
94
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 50 PetscFunctionBegin;
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.
 50 PetscCall(PetscRandomCreate(PETSC_COMM_SELF,&rand));
96 /* get initial s0 so that T^(1/2^s0) < xvals(mmax) */
97 50 *s = 0;
98 220 do {
99 220 inrm = SlepcAbsEigenvalue(wr[0]-PetscRealConstant(1.0),wi[0]);
100
4/4
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
✓ Branch 2 taken 10 times.
✓ Branch 3 taken 10 times.
 16580 for (i=1;i<n;i++) inrm = PetscMax(inrm,SlepcAbsEigenvalue(wr[i]-PetscRealConstant(1.0),wi[i]));
101
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 220 if (inrm < xvals[mmax-1]) break;
102
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 12920 for (i=0;i<n;i++) {
103 #if defined(PETSC_USE_COMPLEX)
104 8250 wr[i] = PetscSqrtScalar(wr[i]);
105 #else
106 #if defined(PETSC_HAVE_COMPLEX)
107 4500 PetscComplex z = PetscSqrtComplex(PetscCMPLX(wr[i],wi[i]));
108 4500 wr[i] = PetscRealPartComplex(z);
109 4500 wi[i] = PetscImaginaryPartComplex(z);
110 #else
111 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"This operation requires a compiler with C99 or C++ complex support");
112 #endif
113 #endif
114 }
115 170 *s = *s + 1;
116
1/2
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
 170 } while (*s<maxroots);
117 50 s0 = *s;
118
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.
 50 if (*s == maxroots) PetscCall(PetscInfo(fn,"Too many matrix square roots\n"));
119
120 /* Troot = T */
121
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.
 3800 for (j=0;j<n;j++) PetscCall(PetscArraycpy(Troot+j*ld,T+j*ld,PetscMin(j+2,n)));
122
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.
 220 for (k=1;k<=PetscMin(*s,maxroots);k++) PetscCall(FNSqrtmSchur(fn,n,Troot,ld,PETSC_FALSE));
123 /* Compute value of s and m needed */
124 /* TrootmI = Troot - I */
125
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 3800 for (j=0;j<n;j++) {
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.
 3750 PetscCall(PetscArraycpy(TrootmI+j*ld,Troot+j*ld,PetscMin(j+2,n)));
127 3750 TrootmI[j+j*ld] -= 1.0;
128 }
129
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 50 PetscCall(SlepcNormAm(n,TrootmI,2,work,rand,&d2));
130 50 d2 = PetscPowReal(d2,1.0/2.0);
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.
 50 PetscCall(SlepcNormAm(n,TrootmI,3,work,rand,&d3));
132 50 d3 = PetscPowReal(d3,1.0/3.0);
133
1/2
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
 50 a2 = PetscMax(d2,d3);
134
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
 50 if (a2 <= xvals[1]) {
135 ✗ *m = (a2 <= xvals[0])? 1: 2;
136 ✗ foundm = PETSC_TRUE;
137 }
138 50 p = 0;
139
1/2
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
 120 while (!foundm) {
140 120 more = PETSC_FALSE; /* More norm checks needed */
141
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 120 if (*s > s0) {
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.
 70 PetscCall(SlepcNormAm(n,TrootmI,3,work,rand,&d3));
143 70 d3 = PetscPowReal(d3,1.0/3.0);
144 }
145
4/6
✓ Branch 0 taken 2 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(SlepcNormAm(n,TrootmI,4,work,rand,&d4));
146 120 d4 = PetscPowReal(d4,1.0/4.0);
147
1/2
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
 120 a3 = PetscMax(d3,d4);
148
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 120 if (a3 <= xvals[mmax-1]) {
149
3/4
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
✓ Branch 2 taken 10 times.
✗ Branch 3 not taken.
 330 for (j=2;j<mmax;j++) if (a3 <= xvals[j]) break;
150
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 70 if (j <= 5) {
151 20 *m = j+1;
152 20 break;
153
3/4
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
✓ Branch 2 taken 10 times.
✗ Branch 3 not taken.
 50 } else if (a3/2.0 <= xvals[4] && p < 2) {
154 20 more = PETSC_TRUE;
155 20 p = p + 1;
156 }
157 }
158
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 100 if (!more) {
159
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 80 PetscCall(SlepcNormAm(n,TrootmI,5,work,rand,&d5));
160 80 d5 = PetscPowReal(d5,1.0/5.0);
161
1/2
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
 80 a4 = PetscMax(d4,d5);
162
1/2
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
 80 eta = PetscMin(a3,a4);
163
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 80 if (eta <= xvals[mmax-1]) {
164
3/4
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
✓ Branch 2 taken 10 times.
✗ Branch 3 not taken.
 50 for (j=5;j<mmax;j++) if (eta <= xvals[j]) break;
165 30 *m = j + 1;
166 30 break;
167 }
168 }
169
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
 70 if (*s == maxroots) {
170 ✗ PetscCall(PetscInfo(fn,"Too many matrix square roots\n"));
171 ✗ *m = mmax; /* No good value found so take largest */
172 ✗ break;
173 }
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.
 70 PetscCall(FNSqrtmSchur(fn,n,Troot,ld,PETSC_FALSE));
175 /* TrootmI = Troot - I */
176
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 5320 for (j=0;j<n;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.
 5250 PetscCall(PetscArraycpy(TrootmI+j*ld,Troot+j*ld,PetscMin(j+2,n)));
178 5250 TrootmI[j+j*ld] -= 1.0;
179 }
180 70 *s = *s + 1;
181 }
182
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 50 PetscCall(PetscRandomDestroy(&rand));
183
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.
 10 PetscFunctionReturn(PETSC_SUCCESS);
184 }
185
186 #if !defined(PETSC_USE_COMPLEX)
187 /*
188 Computes a^(1/2^s) - 1 accurately, avoiding subtractive cancellation
189 */
190 10 static PetscScalar sqrt_obo(PetscScalar a,PetscInt s)
191 {
192 10 PetscScalar val,z0,r;
193 10 PetscReal angle;
194 10 PetscInt i,n0;
195
196
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 10 PetscFunctionBegin;
197
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
 10 if (s == 0) val = a-1.0;
198 else {
199 10 n0 = s;
200 10 angle = PetscAtan2Real(PetscImaginaryPart(a),PetscRealPart(a));
201
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
 10 if (angle >= PETSC_PI/2.0) {
202 ✗ a = PetscSqrtScalar(a);
203 ✗ n0 = s - 1;
204 }
205 10 z0 = a - 1.0;
206 10 a = PetscSqrtScalar(a);
207 10 r = 1.0 + a;
208
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
 40 for (i=0;i<n0-1;i++) {
209 30 a = PetscSqrtScalar(a);
210 30 r = r*(1.0+a);
211 }
212 10 val = z0/r;
213 }
214
6/12
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 1 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 1 times.
 10 PetscFunctionReturn(val);
215 }
216 #endif
217
218 /*
219 Square root of 2x2 matrix T from block diagonal of Schur form. Overwrites T.
220 */
221 16280 static PetscErrorCode sqrtm_tbt(PetscScalar *T)
222 {
223 16280 PetscScalar t11,t12,t21,t22,r11,r22;
224 #if !defined(PETSC_USE_COMPLEX)
225 5180 PetscScalar mu;
226 #endif
227
228
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 16280 PetscFunctionBegin;
229 16280 t11 = T[0]; t21 = T[1]; t12 = T[2]; t22 = T[3];
230
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 10 times.
 16280 if (t21 != 0.0) {
231 /* Compute square root of 2x2 quasitriangular block */
232 /* The algorithm assumes the special structure of real Schur form */
233 #if defined(PETSC_USE_COMPLEX)
234 ✗ SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Should not reach this line in complex scalars");
235 #else
236 1480 mu = PetscSqrtReal(-t21*t12);
237
1/2
✓ Branch 0 taken 5 times.
✗ Branch 1 not taken.
 1480 if (t11 > 0.0) r11 = PetscSqrtReal((t11+SlepcAbsEigenvalue(t11,mu))/2.0);
238 ✗ else r11 = mu / PetscSqrtReal(2.0*(-t11+SlepcAbsEigenvalue(t11,mu)));
239 1480 T[0] = r11;
240 1480 T[1] = t21/(2.0*r11);
241 1480 T[2] = t12/(2.0*r11);
242 1480 T[3] = r11;
243 #endif
244 } else {
245 /* Compute square root of 2x2 upper triangular block */
246 14800 r11 = PetscSqrtScalar(t11);
247 14800 r22 = PetscSqrtScalar(t22);
248 14800 T[0] = r11;
249 14800 T[2] = t12/(r11+r22);
250 14800 T[3] = r22;
251 }
252
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.
 16280 PetscFunctionReturn(PETSC_SUCCESS);
253 }
254
255 #if defined(PETSC_USE_COMPLEX)
256 /*
257 Unwinding number of z
258 */
259 2460 static inline PetscReal unwinding(PetscScalar z)
260 {
261 2460 PetscReal u;
262
263
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 2460 PetscFunctionBegin;
264 2460 u = PetscCeilReal((PetscImaginaryPart(z)-PETSC_PI)/(2.0*PETSC_PI));
265
6/12
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 1 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 1 times.
 2460 PetscFunctionReturn(u);
266 }
267 #endif
268
269 /*
270 Power of 2-by-2 upper triangular matrix A.
271 Returns the (1,2) element of the pth power of A, where p is an arbitrary real number
272 */
273 2460 static PetscScalar powerm2by2(PetscScalar A11,PetscScalar A22,PetscScalar A12,PetscReal p)
274 {
275 2460 PetscScalar a1,a2,a1p,a2p,loga1,loga2,w,dd,x12;
276
277
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 2460 PetscFunctionBegin;
278 2460 a1 = A11;
279 2460 a2 = A22;
280
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 2460 if (a1 == a2) {
281 1460 x12 = p*A12*PetscPowScalarReal(a1,p-1);
282
3/4
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 5 times.
✓ Branch 3 taken 10 times.
 1000 } else if (PetscAbsScalar(a1) < 0.5*PetscAbsScalar(a2) || PetscAbsScalar(a2) < 0.5*PetscAbsScalar(a1)) {
283 10 a1p = PetscPowScalarReal(a1,p);
284 10 a2p = PetscPowScalarReal(a2,p);
285 10 x12 = A12*(a2p-a1p)/(a2-a1);
286 } else { /* Close eigenvalues */
287 990 loga1 = PetscLogScalar(a1);
288 990 loga2 = PetscLogScalar(a2);
289 990 w = PetscAtanhScalar((a2-a1)/(a2+a1));
290 #if defined(PETSC_USE_COMPLEX)
291 980 w += PETSC_i*PETSC_PI*unwinding(loga2-loga1);
292 #endif
293 990 dd = 2.0*PetscExpScalar((loga1+loga2)*p/2.0)*PetscSinhScalar(p*w)/(a2-a1);
294 990 x12 = A12*dd;
295 }
296
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.
 2460 PetscFunctionReturn(x12);
297 }
298
299 /*
300 Recomputes diagonal blocks of T = X^(1/2^s) - 1 more accurately
301 */
302 50 static PetscErrorCode recompute_diag_blocks_sqrt(PetscBLASInt n,PetscScalar *Troot,PetscScalar *T,PetscBLASInt ld,PetscInt *blockStruct,PetscInt s)
303 {
304 50 PetscScalar A[4],P[4],M[4],Z0[4],det;
305 50 PetscInt i,j;
306 #if !defined(PETSC_USE_COMPLEX)
307 20 PetscInt last_block=0;
308 20 PetscScalar a;
309 #endif
310
311
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 50 PetscFunctionBegin;
312
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 3750 for (j=0;j<n-1;j++) {
313 #if !defined(PETSC_USE_COMPLEX)
314
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
 1480 switch (blockStruct[j]) {
315 370 case 0: /* Not start of a block */
316
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
 370 if (last_block != 0) {
317 last_block = 0;
318 } else { /* In a 1x1 block */
319 10 a = T[j+j*ld];
320 10 Troot[j+j*ld] = sqrt_obo(a,s);
321 }
322 break;
323 1110 default: /* In a 2x2 block */
324 1110 last_block = blockStruct[j];
325 #endif
326
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
 3330 if (s == 0) {
327 ✗ Troot[j+j*ld] = T[j+j*ld]-1.0;
328 ✗ Troot[j+1+j*ld] = T[j+1+j*ld];
329 ✗ Troot[j+(j+1)*ld] = T[j+(j+1)*ld];
330 ✗ Troot[j+1+(j+1)*ld] = T[j+1+(j+1)*ld]-1.0;
331 ✗ continue;
332 }
333 3330 A[0] = T[j+j*ld]; A[1] = T[j+1+j*ld]; A[2] = T[j+(j+1)*ld]; A[3] = T[j+1+(j+1)*ld];
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.
 3330 PetscCall(sqrtm_tbt(A));
335 /* Z0 = A - I */
336 3330 Z0[0] = A[0]-1.0; Z0[1] = A[1]; Z0[2] = A[2]; Z0[3] = A[3]-1.0;
337
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
 3330 if (s == 1) {
338 ✗ Troot[j+j*ld] = Z0[0];
339 ✗ Troot[j+1+j*ld] = Z0[1];
340 ✗ Troot[j+(j+1)*ld] = Z0[2];
341 ✗ Troot[j+1+(j+1)*ld] = Z0[3];
342 ✗ continue;
343 }
344
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 3330 PetscCall(sqrtm_tbt(A));
345 /* P = A + I */
346 3330 P[0] = A[0]+1.0; P[1] = A[1]; P[2] = A[2]; P[3] = A[3]+1.0;
347
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 12950 for (i=0;i<s-2;i++) {
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.
 9620 PetscCall(sqrtm_tbt(A));
349 /* P = P*(I + A) */
350 9620 M[0] = P[0]*(A[0]+1.0)+P[2]*A[1];
351 9620 M[1] = P[1]*(A[0]+1.0)+P[3]*A[1];
352 9620 M[2] = P[0]*A[2]+P[2]*(A[3]+1.0);
353 9620 M[3] = P[1]*A[2]+P[3]*(A[3]+1.0);
354 9620 P[0] = M[0]; P[1] = M[1]; P[2] = M[2]; P[3] = M[3];
355 }
356 /* Troot(j:j+1,j:j+1) = Z0 / P (via Cramer) */
357 3330 det = P[0]*P[3]-P[1]*P[2];
358 3330 Troot[j+j*ld] = (Z0[0]*P[3]-P[1]*Z0[2])/det;
359 3330 Troot[j+(j+1)*ld] = (P[0]*Z0[2]-Z0[0]*P[2])/det;
360 3330 Troot[j+1+j*ld] = (Z0[1]*P[3]-P[1]*Z0[3])/det;
361 3330 Troot[j+1+(j+1)*ld] = (P[0]*Z0[3]-Z0[1]*P[2])/det;
362 /* If block is upper triangular recompute the (1,2) element.
363 Skip when T(j,j) or T(j+1,j+1) < 0 since the implementation of atanh is nonstandard */
364
6/6
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 5 times.
✓ Branch 2 taken 10 times.
✓ Branch 3 taken 5 times.
✓ Branch 4 taken 10 times.
✓ Branch 5 taken 5 times.
 3330 if (T[j+1+j*ld]==0.0 && PetscRealPart(T[j+j*ld])>=0.0 && PetscRealPart(T[j+1+(j+1)*ld])>=0.0) {
365 2460 Troot[j+(j+1)*ld] = powerm2by2(T[j+j*ld],T[j+1+(j+1)*ld],T[j+(j+1)*ld],1.0/PetscPowInt(2,s));
366 }
367 #if !defined(PETSC_USE_COMPLEX)
368 }
369 #endif
370 }
371 #if !defined(PETSC_USE_COMPLEX)
372 /* If last diagonal entry is not in a block it will have been missed */
373
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
 20 if (blockStruct[n-2] == 0) {
374 ✗ a = T[n-1+(n-1)*ld];
375 ✗ Troot[n-1+(n-1)*ld] = sqrt_obo(a,s);
376 }
377 #endif
378
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.
 10 PetscFunctionReturn(PETSC_SUCCESS);
379 }
380
381 /*
382 Nodes x and weights w for n-point Gauss-Legendre quadrature (Q is n*n workspace)
383
384 G. H. Golub and J. H. Welsch, Calculation of Gauss quadrature
385 rules, Math. Comp., 23(106):221-230, 1969.
386 */
387 50 static PetscErrorCode gauss_legendre(PetscBLASInt n,PetscScalar *x,PetscScalar *w,PetscScalar *Q)
388 {
389 50 PetscScalar v,a,*work;
390 50 PetscReal *eig,dummy;
391 50 PetscBLASInt k,ld=n,lwork,info;
392 #if defined(PETSC_USE_COMPLEX)
393 30 PetscReal *rwork,rdummy;
394 #endif
395
396
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 50 PetscFunctionBegin;
397
4/6
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 50 PetscCall(PetscArrayzero(Q,n*n));
398
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 320 for (k=1;k<n;k++) {
399 270 v = k/PetscSqrtReal(4.0*k*k-1.0);
400 270 Q[k+(k-1)*n] = v;
401 270 Q[(k-1)+k*n] = v;
402 }
403
404 /* workspace query and memory allocation */
405 50 lwork = -1;
406 #if defined(PETSC_USE_COMPLEX)
407
10/20
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 1 times.
✓ Branch 8 taken 1 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 1 times.
✓ Branch 12 taken 1 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 1 times.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✓ Branch 17 taken 1 times.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
 30 PetscCallBLAS("LAPACKsyev",LAPACKsyev_("V","L",&n,Q,&ld,&dummy,&a,&lwork,&rdummy,&info));
408
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 30 PetscCall(PetscBLASIntCast((PetscInt)PetscRealPart(a),&lwork));
409
5/8
✓ Branch 0 taken 5 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1 times.
✓ Branch 3 taken 4 times.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 1 times.
 30 PetscCall(PetscMalloc3(n,&eig,lwork,&work,PetscMax(1,3*n-2),&rwork));
410 #else
411
10/20
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 1 times.
✓ Branch 8 taken 1 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 1 times.
✓ Branch 12 taken 1 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 1 times.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✓ Branch 17 taken 1 times.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
 20 PetscCallBLAS("LAPACKsyev",LAPACKsyev_("V","L",&n,Q,&ld,&dummy,&a,&lwork,&info));
412
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 20 PetscCall(PetscBLASIntCast((PetscInt)a,&lwork));
413
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 20 PetscCall(PetscMalloc2(n,&eig,lwork,&work));
414 #endif
415
416 /* compute eigendecomposition */
417 #if defined(PETSC_USE_COMPLEX)
418
10/20
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 1 times.
✓ Branch 8 taken 1 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 1 times.
✓ Branch 12 taken 1 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 1 times.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✓ Branch 17 taken 1 times.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
 30 PetscCallBLAS("LAPACKsyev",LAPACKsyev_("V","L",&n,Q,&ld,eig,work,&lwork,rwork,&info));
419 #else
420
10/20
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 1 times.
✓ Branch 8 taken 1 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 1 times.
✓ Branch 12 taken 1 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 1 times.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✓ Branch 17 taken 1 times.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
 20 PetscCallBLAS("LAPACKsyev",LAPACKsyev_("V","L",&n,Q,&ld,eig,work,&lwork,&info));
421 #endif
422
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 50 SlepcCheckLapackInfo("syev",info);
423
424
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 370 for (k=0;k<n;k++) {
425 320 x[k] = eig[k];
426 320 w[k] = 2.0*Q[k*n]*Q[k*n];
427 }
428 #if defined(PETSC_USE_COMPLEX)
429
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 30 PetscCall(PetscFree3(eig,work,rwork));
430 #else
431
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 20 PetscCall(PetscFree2(eig,work));
432 #endif
433
3/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 taken 2 times.
 50 PetscCall(PetscLogFlops(9.0*n*n*n+2.0*n*n*n));
434
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.
 50 PetscFunctionReturn(PETSC_SUCCESS);
435 }
436
437 /*
438 Pade approximation to log(1 + T) via partial fractions
439 */
440 50 static PetscErrorCode pade_approx(PetscBLASInt n,PetscScalar *T,PetscScalar *L,PetscBLASInt ld,PetscInt m,PetscScalar *work)
441 {
442 50 PetscScalar *K,*W,*nodes,*wts;
443 50 PetscBLASInt *ipiv,info,mm;
444 50 PetscInt i,j,k;
445
446
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 50 PetscFunctionBegin;
447 50 K = work;
448 50 W = work+n*n;
449 50 nodes = work+2*n*n;
450 50 wts = work+2*n*n+m;
451 50 ipiv = (PetscBLASInt*)(work+2*n*n+2*m);
452
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 50 PetscCall(PetscBLASIntCast(m,&mm));
453
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 50 PetscCall(gauss_legendre(mm,nodes,wts,L));
454 /* Convert from [-1,1] to [0,1] */
455
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 370 for (i=0;i<m;i++) {
456 320 nodes[i] = (nodes[i]+1.0)/2.0;
457 320 wts[i] = wts[i]/2.0;
458 }
459
4/6
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 50 PetscCall(PetscArrayzero(L,n*n));
460
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 370 for (k=0;k<m;k++) {
461
4/4
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
✓ Branch 2 taken 10 times.
✓ Branch 3 taken 10 times.
 1824320 for (i=0;i<n;i++) for (j=0;j<n;j++) K[i+j*ld] = nodes[k]*T[i+j*ld];
462
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 24320 for (i=0;i<n;i++) K[i+i*ld] += 1.0;
463
4/4
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
✓ Branch 2 taken 10 times.
✓ Branch 3 taken 10 times.
 1824320 for (i=0;i<n;i++) for (j=0;j<n;j++) W[i+j*ld] = T[i+j*ld];
464
10/20
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
✓ Branch 8 taken 2 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 2 times.
✓ Branch 12 taken 2 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 2 times.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✓ Branch 17 taken 2 times.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
 320 PetscCallBLAS("LAPACKgesv",LAPACKgesv_(&n,&n,K,&n,ipiv,W,&n,&info));
465
4/4
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
✓ Branch 2 taken 10 times.
✓ Branch 3 taken 10 times.
 1824320 for (i=0;i<n;i++) for (j=0;j<n;j++) L[i+j*ld] += wts[k]*W[i+j*ld];
466 }
467
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.
 10 PetscFunctionReturn(PETSC_SUCCESS);
468 }
469
470 /*
471 Recomputes diagonal blocks of T = X^(1/2^s) - 1 more accurately
472 */
473 50 static PetscErrorCode recompute_diag_blocks_log(PetscBLASInt n,PetscScalar *L,PetscScalar *T,PetscBLASInt ld,PetscInt *blockStruct)
474 {
475 50 PetscScalar a1,a2,a12,loga1,loga2,z,dd;
476 50 PetscInt j;
477 #if !defined(PETSC_USE_COMPLEX)
478 20 PetscInt last_block=0;
479 20 PetscScalar f,t;
480 #endif
481
482
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 50 PetscFunctionBegin;
483
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 3750 for (j=0;j<n-1;j++) {
484 #if !defined(PETSC_USE_COMPLEX)
485
3/4
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
✓ Branch 2 taken 5 times.
✗ Branch 3 not taken.
 1480 switch (blockStruct[j]) {
486 370 case 0: /* Not start of a block */
487
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
 370 if (last_block != 0) {
488 last_block = 0;
489 } else { /* In a 1x1 block */
490 10 L[j+j*ld] = PetscLogScalar(T[j+j*ld]);
491 }
492 break;
493 740 case 1: /* Start of upper-tri block */
494 740 last_block = 1;
495 #endif
496 2960 a1 = T[j+j*ld];
497 2960 a2 = T[j+1+(j+1)*ld];
498 2960 loga1 = PetscLogScalar(a1);
499 2960 loga2 = PetscLogScalar(a2);
500 2960 L[j+j*ld] = loga1;
501 2960 L[j+1+(j+1)*ld] = loga2;
502
5/6
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 5 times.
✓ Branch 2 taken 5 times.
✓ Branch 3 taken 10 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 5 times.
 2960 if ((PetscRealPart(a1)<0.0 && PetscImaginaryPart(a1)==0.0) || (PetscRealPart(a2)<0.0 && PetscImaginaryPart(a1)==0.0)) {
503 /* Problems with 2 x 2 formula for (1,2) block
504 since atanh is nonstandard, just redo diagonal part */
505 ✗ continue;
506 }
507
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 2960 if (a1 == a2) {
508 1460 a12 = T[j+(j+1)*ld]/a1;
509
3/4
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 5 times.
✓ Branch 3 taken 10 times.
 1500 } else if (PetscAbsScalar(a1)<0.5*PetscAbsScalar(a2) || PetscAbsScalar(a2)<0.5*PetscAbsScalar(a1)) {
510 10 a12 = T[j+(j+1)*ld]*(loga2-loga1)/(a2-a1);
511 } else { /* Close eigenvalues */
512 1490 z = (a2-a1)/(a2+a1);
513 1490 dd = 2.0*PetscAtanhScalar(z);
514 #if defined(PETSC_USE_COMPLEX)
515 1480 dd += 2.0*PETSC_i*PETSC_PI*unwinding(loga2-loga1);
516 #endif
517 1490 dd /= (a2-a1);
518 1490 a12 = T[j+(j+1)*ld]*dd;
519 }
520 2960 L[j+(j+1)*ld] = a12;
521 #if !defined(PETSC_USE_COMPLEX)
522 740 break;
523 370 case 2: /* Start of quasi-tri block */
524 370 last_block = 2;
525 370 f = 0.5*PetscLogScalar(T[j+j*ld]*T[j+j*ld]-T[j+(j+1)*ld]*T[j+1+j*ld]);
526 370 t = PetscAtan2Real(PetscSqrtScalar(-T[j+(j+1)*ld]*T[j+1+j*ld]),T[j+j*ld])/PetscSqrtScalar(-T[j+(j+1)*ld]*T[j+1+j*ld]);
527 370 L[j+j*ld] = f;
528 370 L[j+1+j*ld] = t*T[j+1+j*ld];
529 370 L[j+(j+1)*ld] = t*T[j+(j+1)*ld];
530 370 L[j+1+(j+1)*ld] = f;
531 }
532 #endif
533 }
534
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.
 50 PetscFunctionReturn(PETSC_SUCCESS);
535 }
536 /*
537 * Matrix logarithm implementation based on algorithm and matlab code by N. Higham and co-authors
538 *
539 * H. Al-Mohy and N. J. Higham, "Improved inverse scaling and squaring
540 * algorithms for the matrix logarithm", SIAM J. Sci. Comput. 34(4):C153-C169, 2012.
541 */
542 50 static PetscErrorCode FNLogmPade(FN fn,PetscBLASInt n,PetscScalar *T,PetscBLASInt ld,PetscBool firstonly)
543 {
544 50 PetscBLASInt k,sdim,lwork,info;
545 50 PetscScalar *wr,*wi=NULL,*W,*Q,*Troot,*L,*work,one=1.0,zero=0.0,alpha;
546 50 PetscInt i,j,s=0,m=0,*blockformat;
547 #if defined(PETSC_USE_COMPLEX)
548 30 PetscReal *rwork;
549 #endif
550
551
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 50 PetscFunctionBegin;
552 50 lwork = 3*n*n; /* gees needs only 5*n, but work is also passed to FNlogm_params */
553
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 50 k = firstonly? 1: n;
554
555 /* compute Schur decomposition A*Q = Q*T */
556
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 50 PetscCall(PetscCalloc7(n,&wr,n*k,&W,n*n,&Q,n*n,&Troot,n*n,&L,lwork,&work,n-1,&blockformat));
557 #if !defined(PETSC_USE_COMPLEX)
558
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 20 PetscCall(PetscMalloc1(n,&wi));
559
10/20
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 1 times.
✓ Branch 8 taken 1 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 1 times.
✓ Branch 12 taken 1 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 1 times.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✓ Branch 17 taken 1 times.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
 20 PetscCallBLAS("LAPACKgees",LAPACKgees_("V","N",NULL,&n,T,&ld,&sdim,wr,wi,Q,&ld,work,&lwork,NULL,&info));
560 #else
561
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 30 PetscCall(PetscMalloc1(n,&rwork));
562
10/20
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 1 times.
✓ Branch 8 taken 1 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 1 times.
✓ Branch 12 taken 1 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 1 times.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✓ Branch 17 taken 1 times.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
 30 PetscCallBLAS("LAPACKgees",LAPACKgees_("V","N",NULL,&n,T,&ld,&sdim,wr,Q,&ld,work,&lwork,rwork,NULL,&info));
563 #endif
564
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 50 SlepcCheckLapackInfo("gees",info);
565
566 #if !defined(PETSC_USE_COMPLEX)
567 /* check for negative real eigenvalues */
568
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
 1520 for (i=0;i<n;i++) {
569
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 1500 PetscCheck(wr[i]>=0.0 || wi[i]!=0.0,PETSC_COMM_SELF,PETSC_ERR_SUP,"Matrix has negative real eigenvalue; rerun with complex scalars");
570 }
571 #endif
572
573 /* get block structure of Schur factor */
574
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 50 PetscCall(qtri_struct(n,T,ld,blockformat));
575
576 /* get parameters */
577
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 50 PetscCall(FNlogm_params(fn,n,T,ld,wr,wi,100,&s,&m,Troot,work));
578
579 /* compute Troot - I = T(1/2^s) - I more accurately */
580
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 50 PetscCall(recompute_diag_blocks_sqrt(n,Troot,T,ld,blockformat,s));
581
582 /* compute Pade approximant */
583
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 50 PetscCall(pade_approx(n,Troot,L,ld,m,work));
584
585 /* scale back up, L = 2^s * L */
586 50 alpha = PetscPowInt(2,s);
587
4/4
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
✓ Branch 2 taken 10 times.
✓ Branch 3 taken 10 times.
 285050 for (i=0;i<n;i++) for (j=0;j<n;j++) L[i+j*ld] *= alpha;
588
589 /* recompute diagonal blocks */
590
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 50 PetscCall(recompute_diag_blocks_log(n,L,T,ld,blockformat));
591
592 /* backtransform B = Q*L*Q' */
593
10/20
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
✓ Branch 8 taken 2 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 2 times.
✓ Branch 12 taken 2 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 2 times.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✓ Branch 17 taken 2 times.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
 50 PetscCallBLAS("BLASgemm",BLASgemm_("N","C",&n,&k,&n,&one,L,&ld,Q,&ld,&zero,W,&ld));
594
10/20
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
✓ Branch 8 taken 2 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 2 times.
✓ Branch 12 taken 2 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 2 times.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✓ Branch 17 taken 2 times.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
 50 PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n,&k,&n,&one,Q,&ld,W,&ld,&zero,T,&ld));
595
596 /* flop count: Schur decomposition, and backtransform */
597
3/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 taken 2 times.
 50 PetscCall(PetscLogFlops(25.0*n*n*n+4.0*n*n*k));
598
599
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 50 PetscCall(PetscFree7(wr,W,Q,Troot,L,work,blockformat));
600 #if !defined(PETSC_USE_COMPLEX)
601
5/8
✓ Branch 0 taken 5 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1 times.
✓ Branch 3 taken 4 times.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 1 times.
 20 PetscCall(PetscFree(wi));
602 #else
603
5/8
✓ Branch 0 taken 5 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1 times.
✓ Branch 3 taken 4 times.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 1 times.
 30 PetscCall(PetscFree(rwork));
604 #endif
605
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.
 10 PetscFunctionReturn(PETSC_SUCCESS);
606 }
607
608 25 static PetscErrorCode FNEvaluateFunctionMat_Log_Higham(FN fn,Mat A,Mat B)
609 {
610 25 PetscBLASInt n = 0;
611 25 PetscScalar *T;
612 25 PetscInt m;
613
614
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 25 PetscFunctionBegin;
615
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.
 25 if (A!=B) PetscCall(MatCopy(A,B,SAME_NONZERO_PATTERN));
616
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 25 PetscCall(MatDenseGetArray(B,&T));
617
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 25 PetscCall(MatGetSize(A,&m,NULL));
618
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 25 PetscCall(PetscBLASIntCast(m,&n));
619
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 25 PetscCall(FNLogmPade(fn,n,T,n,PETSC_FALSE));
620
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 25 PetscCall(MatDenseRestoreArray(B,&T));
621
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.
 5 PetscFunctionReturn(PETSC_SUCCESS);
622 }
623
624 25 static PetscErrorCode FNEvaluateFunctionMatVec_Log_Higham(FN fn,Mat A,Vec v)
625 {
626 25 PetscBLASInt n = 0;
627 25 PetscScalar *T;
628 25 PetscInt m;
629 25 Mat B;
630
631
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 25 PetscFunctionBegin;
632
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 25 PetscCall(FN_AllocateWorkMat(fn,A,&B));
633
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 25 PetscCall(MatDenseGetArray(B,&T));
634
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 25 PetscCall(MatGetSize(A,&m,NULL));
635
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 25 PetscCall(PetscBLASIntCast(m,&n));
636
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 25 PetscCall(FNLogmPade(fn,n,T,n,PETSC_TRUE));
637
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 25 PetscCall(MatDenseRestoreArray(B,&T));
638
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 25 PetscCall(MatGetColumnVector(B,v,0));
639
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 25 PetscCall(FN_FreeWorkMat(fn,&B));
640
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.
 5 PetscFunctionReturn(PETSC_SUCCESS);
641 }
642
643 45 static PetscErrorCode FNView_Log(FN fn,PetscViewer viewer)
644 {
645 45 PetscBool isascii;
646 45 char str[50];
647 45 const char *methodname[] = {
648 "scaling & squaring, [m/m] Pade approximant (Higham)"
649 };
650 45 const int nmeth=PETSC_STATIC_ARRAY_LENGTH(methodname);
651
652
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 45 PetscFunctionBegin;
653
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 45 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
654
1/2
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
 45 if (isascii) {
655
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 45 if (fn->beta==(PetscScalar)1.0) {
656
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.
 10 if (fn->alpha==(PetscScalar)1.0) PetscCall(PetscViewerASCIIPrintf(viewer," logarithm: log(x)\n"));
657 else {
658 ✗ PetscCall(SlepcSNPrintfScalar(str,sizeof(str),fn->alpha,PETSC_TRUE));
659 ✗ PetscCall(PetscViewerASCIIPrintf(viewer," logarithm: log(%s*x)\n",str));
660 }
661 } else {
662
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 35 PetscCall(SlepcSNPrintfScalar(str,sizeof(str),fn->beta,PETSC_TRUE));
663
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.
 35 if (fn->alpha==(PetscScalar)1.0) PetscCall(PetscViewerASCIIPrintf(viewer," logarithm: %s*log(x)\n",str));
664 else {
665
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 35 PetscCall(PetscViewerASCIIPrintf(viewer," logarithm: %s",str));
666
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 35 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE));
667
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 35 PetscCall(SlepcSNPrintfScalar(str,sizeof(str),fn->alpha,PETSC_TRUE));
668
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 35 PetscCall(PetscViewerASCIIPrintf(viewer,"*log(%s*x)\n",str));
669
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 35 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE));
670 }
671 }
672
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.
 45 if (fn->method<nmeth) PetscCall(PetscViewerASCIIPrintf(viewer," computing matrix functions with: %s\n",methodname[fn->method]));
673 }
674
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.
 9 PetscFunctionReturn(PETSC_SUCCESS);
675 }
676
677 35 SLEPC_EXTERN PetscErrorCode FNCreate_Log(FN fn)
678 {
679
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 35 PetscFunctionBegin;
680 35 fn->ops->evaluatefunction = FNEvaluateFunction_Log;
681 35 fn->ops->evaluatederivative = FNEvaluateDerivative_Log;
682 35 fn->ops->evaluatefunctionmat[0] = FNEvaluateFunctionMat_Log_Higham;
683 35 fn->ops->evaluatefunctionmatvec[0] = FNEvaluateFunctionMatVec_Log_Higham;
684 35 fn->ops->view = FNView_Log;
685
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.
 35 PetscFunctionReturn(PETSC_SUCCESS);
686 }
687