GCC Code Coverage Report


Directory: ./
File: src/sys/classes/fn/impls/exp/fnexp.c
Date: 2026-07-06 04:10:30
Exec Total Coverage
Lines: 976 1196 81.6%
Functions: 16 17 94.1%
Branches: 1821 4490 40.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 Exponential function exp(x)
12 */
13
14 #include <slepc/private/fnimpl.h> /*I "slepcfn.h" I*/
15 #include <slepcblaslapack.h>
16
17 28278 static PetscErrorCode FNEvaluateFunction_Exp(FN fn,PetscScalar x,PetscScalar *y)
18 {
19
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
28278 PetscFunctionBegin;
20 28278 *y = PetscExpScalar(x);
21
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.
28278 PetscFunctionReturn(PETSC_SUCCESS);
22 }
23
24 8907 static PetscErrorCode FNEvaluateDerivative_Exp(FN fn,PetscScalar x,PetscScalar *y)
25 {
26
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
8907 PetscFunctionBegin;
27 8907 *y = PetscExpScalar(x);
28
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.
8907 PetscFunctionReturn(PETSC_SUCCESS);
29 }
30
31 #define MAX_PADE 6
32
33 160 static PetscErrorCode FNEvaluateFunctionMat_Exp_Pade(FN fn,Mat A,Mat B)
34 {
35 160 PetscBLASInt n=0,ld,ld2,*ipiv,inc=1;
36 160 PetscInt m,j,k,sexp;
37 160 PetscBool odd;
38 160 const PetscInt p=MAX_PADE;
39 160 PetscReal c[MAX_PADE+1],s,*rwork;
40 160 PetscScalar scale,mone=-1.0,one=1.0,two=2.0,zero=0.0;
41 160 PetscScalar *Ba,*As,*A2,*Q,*P,*W,*aux;
42 160 const PetscScalar *Aa;
43
44
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
160 PetscFunctionBegin;
45
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(MatDenseGetArrayRead(A,&Aa));
46
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(MatDenseGetArray(B,&Ba));
47
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(A,&m,NULL));
48
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(PetscBLASIntCast(m,&n));
49 160 ld = n;
50 160 ld2 = ld*ld;
51 160 P = Ba;
52
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(PetscMalloc6(m*m,&Q,m*m,&W,m*m,&As,m*m,&A2,ld,&rwork,ld,&ipiv));
53
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(PetscArraycpy(As,Aa,ld2));
54
55 /* Pade' coefficients */
56 160 c[0] = 1.0;
57
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
1120 for (k=1;k<=p;k++) c[k] = c[k-1]*(p+1-k)/(k*(2*p+1-k));
58
59 /* Scaling */
60 160 s = LAPACKlange_("I",&n,&n,As,&ld,rwork);
61
3/6
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
160 PetscCall(PetscLogFlops(1.0*n*n));
62
1/2
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
160 if (s>0.5) {
63
1/2
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
160 sexp = PetscMax(0,(int)(PetscLogReal(s)/PetscLogReal(2.0))+2);
64 160 scale = PetscPowRealInt(2.0,-sexp);
65
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.
160 PetscCallBLAS("BLASscal",BLASscal_(&ld2,&scale,As,&inc));
66
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.
160 PetscCall(PetscLogFlops(1.0*n*n));
67 } else sexp = 0;
68
69 /* Horner evaluation */
70
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.
160 PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n,&n,&n,&one,As,&ld,As,&ld,&zero,A2,&ld));
71
3/6
✓ Branch 0 taken 5 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
160 PetscCall(PetscLogFlops(2.0*n*n*n));
72
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.
160 PetscCall(PetscArrayzero(Q,ld2));
73
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.
160 PetscCall(PetscArrayzero(P,ld2));
74
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
6160 for (j=0;j<n;j++) {
75 6000 Q[j+j*ld] = c[p];
76 6000 P[j+j*ld] = c[p-1];
77 }
78
79 odd = PETSC_TRUE;
80
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
960 for (k=p-1;k>0;k--) {
81
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
800 if (odd) {
82
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.
480 PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n,&n,&n,&one,Q,&ld,A2,&ld,&zero,W,&ld));
83 480 SlepcSwap(Q,W,aux);
84
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
18480 for (j=0;j<n;j++) Q[j+j*ld] += c[k-1];
85 odd = PETSC_FALSE;
86 } else {
87
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("BLASgemm",BLASgemm_("N","N",&n,&n,&n,&one,P,&ld,A2,&ld,&zero,W,&ld));
88 320 SlepcSwap(P,W,aux);
89
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
12320 for (j=0;j<n;j++) P[j+j*ld] += c[k-1];
90 odd = PETSC_TRUE;
91 }
92
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.
800 PetscCall(PetscLogFlops(2.0*n*n*n));
93 }
94 /*if (odd) {
95 PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n,&n,&n,&one,Q,&ld,As,&ld,&zero,W,&ld));
96 SlepcSwap(Q,W,aux);
97 PetscCallBLAS("BLASaxpy",BLASaxpy_(&ld2,&mone,P,&inc,Q,&inc));
98 PetscCallLAPACKInfo("LAPACKgesv",LAPACKgesv_(&n,&n,Q,&ld,ipiv,P,&ld,&info));
99 PetscCallBLAS("BLASscal",BLASscal_(&ld2,&two,P,&inc));
100 for (j=0;j<n;j++) P[j+j*ld] += 1.0;
101 PetscCallBLAS("BLASscal",BLASscal_(&ld2,&mone,P,&inc));
102 } else {*/
103
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.
160 PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n,&n,&n,&one,P,&ld,As,&ld,&zero,W,&ld));
104 160 SlepcSwap(P,W,aux);
105
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.
160 PetscCallBLAS("BLASaxpy",BLASaxpy_(&ld2,&mone,P,&inc,Q,&inc));
106
13/28
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 8 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 10 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.
✗ Branch 20 not taken.
✓ Branch 21 taken 2 times.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✓ Branch 25 taken 2 times.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
160 PetscCallLAPACKInfo("LAPACKgesv",LAPACKgesv_(&n,&n,Q,&ld,ipiv,P,&ld,&info));
107
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.
160 PetscCallBLAS("BLASscal",BLASscal_(&ld2,&two,P,&inc));
108
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
6160 for (j=0;j<n;j++) P[j+j*ld] += 1.0;
109 /*}*/
110
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.
160 PetscCall(PetscLogFlops(2.0*n*n*n+2.0*n*n*n/3.0+4.0*n*n));
111
112
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
680 for (k=1;k<=sexp;k++) {
113
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.
520 PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n,&n,&n,&one,P,&ld,P,&ld,&zero,W,&ld));
114
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.
520 PetscCall(PetscArraycpy(P,W,ld2));
115 }
116
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.
160 if (P!=Ba) PetscCall(PetscArraycpy(Ba,P,ld2));
117
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.
160 PetscCall(PetscLogFlops(2.0*n*n*n*sexp));
118
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.
160 PetscCall(PetscFree6(Q,W,As,A2,rwork,ipiv));
120
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(MatDenseRestoreArrayRead(A,&Aa));
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.
160 PetscCall(MatDenseRestoreArray(B,&Ba));
122
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.
32 PetscFunctionReturn(PETSC_SUCCESS);
123 }
124
125 #if defined(PETSC_HAVE_COMPLEX)
126 /*
127 * Set scaling factor (s) and Pade degree (k,m)
128 */
129 552 static PetscErrorCode sexpm_params(PetscReal nrm,PetscInt *s,PetscInt *k,PetscInt *m)
130 {
131
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
552 PetscFunctionBegin;
132
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
552 if (nrm>1) {
133
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
352 if (nrm<200) {*s = 4; *k = 5; *m = *k-1;}
134
1/2
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
40 else if (nrm<1e4) {*s = 4; *k = 4; *m = *k+1;}
135 else if (nrm<1e6) {*s = 4; *k = 3; *m = *k+1;}
136 else if (nrm<1e9) {*s = 3; *k = 3; *m = *k+1;}
137 else if (nrm<1e11) {*s = 2; *k = 3; *m = *k+1;}
138 else if (nrm<1e12) {*s = 2; *k = 2; *m = *k+1;}
139 else if (nrm<1e14) {*s = 2; *k = 1; *m = *k+1;}
140 else {*s = 1; *k = 1; *m = *k+1;}
141 } else { /* nrm<1 */
142
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
200 if (nrm>0.5) {*s = 4; *k = 4; *m = *k-1;}
143
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
160 else if (nrm>0.3) {*s = 3; *k = 4; *m = *k-1;}
144
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
120 else if (nrm>0.15) {*s = 2; *k = 4; *m = *k-1;}
145
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
120 else if (nrm>0.07) {*s = 1; *k = 4; *m = *k-1;}
146
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
120 else if (nrm>0.01) {*s = 0; *k = 4; *m = *k-1;}
147
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
120 else if (nrm>3e-4) {*s = 0; *k = 3; *m = *k-1;}
148
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
120 else if (nrm>1e-5) {*s = 0; *k = 3; *m = 0;}
149
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
120 else if (nrm>1e-8) {*s = 0; *k = 2; *m = 0;}
150 120 else {*s = 0; *k = 1; *m = 0;}
151 }
152
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.
552 PetscFunctionReturn(PETSC_SUCCESS);
153 }
154
155 /*
156 * Partial fraction form coefficients.
157 * If query, the function returns the size necessary to store the coefficients.
158 */
159 432 static PetscErrorCode getcoeffs(PetscInt k,PetscInt m,PetscComplex *r,PetscComplex *q,PetscComplex *remain,PetscBool query)
160 {
161 432 PetscInt i;
162 432 const PetscComplex /* m == k+1 */
163 432 p1r4[5] = {-1.582680186458572e+01 - 2.412564578224361e+01*PETSC_i,
164 432 -1.582680186458572e+01 + 2.412564578224361e+01*PETSC_i,
165 432 1.499984465975511e+02 + 6.804227952202417e+01*PETSC_i,
166 432 1.499984465975511e+02 - 6.804227952202417e+01*PETSC_i,
167 -2.733432894659307e+02 },
168 432 p1q4[5] = { 3.655694325463550e+00 + 6.543736899360086e+00*PETSC_i,
169 432 3.655694325463550e+00 - 6.543736899360086e+00*PETSC_i,
170 432 5.700953298671832e+00 + 3.210265600308496e+00*PETSC_i,
171 432 5.700953298671832e+00 - 3.210265600308496e+00*PETSC_i,
172 6.286704751729261e+00 },
173 432 p1r3[4] = {-1.130153999597152e+01 + 1.247167585025031e+01*PETSC_i,
174 432 -1.130153999597152e+01 - 1.247167585025031e+01*PETSC_i,
175 432 1.330153999597152e+01 - 6.007173273704750e+01*PETSC_i,
176 432 1.330153999597152e+01 + 6.007173273704750e+01*PETSC_i},
177 432 p1q3[4] = { 3.212806896871536e+00 + 4.773087433276636e+00*PETSC_i,
178 432 3.212806896871536e+00 - 4.773087433276636e+00*PETSC_i,
179 432 4.787193103128464e+00 + 1.567476416895212e+00*PETSC_i,
180 432 4.787193103128464e+00 - 1.567476416895212e+00*PETSC_i},
181 432 p1r2[3] = { 7.648749087422928e+00 + 4.171640244747463e+00*PETSC_i,
182 432 7.648749087422928e+00 - 4.171640244747463e+00*PETSC_i,
183 -1.829749817484586e+01 },
184 432 p1q2[3] = { 2.681082873627756e+00 + 3.050430199247411e+00*PETSC_i,
185 432 2.681082873627756e+00 - 3.050430199247411e+00*PETSC_i,
186 3.637834252744491e+00 },
187 432 p1r1[2] = { 1.000000000000000e+00 - 3.535533905932738e+00*PETSC_i,
188 432 1.000000000000000e+00 + 3.535533905932738e+00*PETSC_i},
189 432 p1q1[2] = { 2.000000000000000e+00 + 1.414213562373095e+00*PETSC_i,
190 432 2.000000000000000e+00 - 1.414213562373095e+00*PETSC_i};
191 432 const PetscComplex /* m == k-1 */
192 432 m1r5[4] = {-1.423367961376821e+02 - 1.385465094833037e+01*PETSC_i,
193 432 -1.423367961376821e+02 + 1.385465094833037e+01*PETSC_i,
194 432 2.647367961376822e+02 - 4.814394493714596e+02*PETSC_i,
195 432 2.647367961376822e+02 + 4.814394493714596e+02*PETSC_i},
196 432 m1q5[4] = { 5.203941240131764e+00 + 5.805856841805367e+00*PETSC_i,
197 432 5.203941240131764e+00 - 5.805856841805367e+00*PETSC_i,
198 432 6.796058759868242e+00 + 1.886649260140217e+00*PETSC_i,
199 432 6.796058759868242e+00 - 1.886649260140217e+00*PETSC_i},
200 432 m1r4[3] = { 2.484269593165883e+01 + 7.460342395992306e+01*PETSC_i,
201 432 2.484269593165883e+01 - 7.460342395992306e+01*PETSC_i,
202 -1.734353918633177e+02 },
203 432 m1q4[3] = { 4.675757014491557e+00 + 3.913489560603711e+00*PETSC_i,
204 432 4.675757014491557e+00 - 3.913489560603711e+00*PETSC_i,
205 5.648485971016893e+00 },
206 432 m1r3[2] = { 2.533333333333333e+01 - 2.733333333333333e+01*PETSC_i,
207 432 2.533333333333333e+01 + 2.733333333333333e+01*PETSC_i},
208 432 m1q3[2] = { 4.000000000000000e+00 + 2.000000000000000e+00*PETSC_i,
209 432 4.000000000000000e+00 - 2.000000000000000e+00*PETSC_i};
210 432 const PetscScalar /* m == k-1 */
211 432 m1remain5[2] = { 2.000000000000000e-01, 9.800000000000000e+00},
212 432 m1remain4[2] = {-2.500000000000000e-01, -7.750000000000000e+00},
213 432 m1remain3[2] = { 3.333333333333333e-01, 5.666666666666667e+00},
214 432 m1remain2[2] = {-0.5, -3.5},
215 432 remain3[4] = {1.0/6.0, 1.0/2.0, 1, 1},
216 432 remain2[3] = {1.0/2.0, 1, 1};
217
218
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
432 PetscFunctionBegin;
219
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
432 if (query) { /* query about buffer's size */
220
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
216 if (m==k+1) {
221 20 *remain = 0;
222 20 *r = *q = k+1;
223
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); /* quick return */
224 }
225
1/2
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
196 if (m==k-1) {
226 196 *remain = 2;
227
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
196 if (k==5) *r = *q = 4;
228
1/2
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
40 else if (k==4) *r = *q = 3;
229 else if (k==3) *r = *q = 2;
230 else if (k==2) *r = *q = 1;
231 }
232
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
196 if (m==0) {
233 *r = *q = 0;
234 *remain = k+1;
235 }
236 } else {
237
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
216 if (m==k+1) {
238
1/2
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
20 if (k==4) {
239
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
120 for (i=0;i<5;i++) { r[i] = p1r4[i]; q[i] = p1q4[i]; }
240 } else if (k==3) {
241 for (i=0;i<4;i++) { r[i] = p1r3[i]; q[i] = p1q3[i]; }
242 } else if (k==2) {
243 for (i=0;i<3;i++) { r[i] = p1r2[i]; q[i] = p1q2[i]; }
244 } else if (k==1) {
245 for (i=0;i<2;i++) { r[i] = p1r1[i]; q[i] = p1q1[i]; }
246 }
247
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); /* quick return */
248 }
249
1/2
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
196 if (m==k-1) {
250
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
196 if (k==5) {
251
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
780 for (i=0;i<4;i++) { r[i] = m1r5[i]; q[i] = m1q5[i]; }
252
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
468 for (i=0;i<2;i++) remain[i] = m1remain5[i];
253
1/2
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
40 } else if (k==4) {
254
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
160 for (i=0;i<3;i++) { r[i] = m1r4[i]; q[i] = m1q4[i]; }
255
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
120 for (i=0;i<2;i++) remain[i] = m1remain4[i];
256 } else if (k==3) {
257 for (i=0;i<2;i++) { r[i] = m1r3[i]; q[i] = m1q3[i]; remain[i] = m1remain3[i]; }
258 } else if (k==2) {
259 r[0] = -13.5; q[0] = 3;
260 for (i=0;i<2;i++) remain[i] = m1remain2[i];
261 }
262 }
263
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
196 if (m==0) {
264 r = q = NULL;
265 if (k==3) {
266 for (i=0;i<4;i++) remain[i] = remain3[i];
267 } else if (k==2) {
268 for (i=0;i<3;i++) remain[i] = remain2[i];
269 }
270 }
271 }
272
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.
72 PetscFunctionReturn(PETSC_SUCCESS);
273 }
274
275 /*
276 * Product form coefficients.
277 * If query, the function returns the size necessary to store the coefficients.
278 */
279 432 static PetscErrorCode getcoeffsproduct(PetscInt k,PetscInt m,PetscComplex *p,PetscComplex *q,PetscComplex *mult,PetscBool query)
280 {
281 432 PetscInt i;
282 432 const PetscComplex /* m == k+1 */
283 432 p1p4[4] = {-5.203941240131764e+00 + 5.805856841805367e+00*PETSC_i,
284 432 -5.203941240131764e+00 - 5.805856841805367e+00*PETSC_i,
285 432 -6.796058759868242e+00 + 1.886649260140217e+00*PETSC_i,
286 432 -6.796058759868242e+00 - 1.886649260140217e+00*PETSC_i},
287 432 p1q4[5] = { 3.655694325463550e+00 + 6.543736899360086e+00*PETSC_i,
288 432 3.655694325463550e+00 - 6.543736899360086e+00*PETSC_i,
289 6.286704751729261e+00 ,
290 432 5.700953298671832e+00 + 3.210265600308496e+00*PETSC_i,
291 432 5.700953298671832e+00 - 3.210265600308496e+00*PETSC_i},
292 432 p1p3[3] = {-4.675757014491557e+00 + 3.913489560603711e+00*PETSC_i,
293 432 -4.675757014491557e+00 - 3.913489560603711e+00*PETSC_i,
294 -5.648485971016893e+00 },
295 432 p1q3[4] = { 3.212806896871536e+00 + 4.773087433276636e+00*PETSC_i,
296 432 3.212806896871536e+00 - 4.773087433276636e+00*PETSC_i,
297 432 4.787193103128464e+00 + 1.567476416895212e+00*PETSC_i,
298 432 4.787193103128464e+00 - 1.567476416895212e+00*PETSC_i},
299 432 p1p2[2] = {-4.00000000000000e+00 + 2.000000000000000e+00*PETSC_i,
300 432 -4.00000000000000e+00 - 2.000000000000000e+00*PETSC_i},
301 432 p1q2[3] = { 2.681082873627756e+00 + 3.050430199247411e+00*PETSC_i,
302 432 2.681082873627756e+00 - 3.050430199247411e+00*PETSC_i,
303 3.637834252744491e+00 },
304 432 p1q1[2] = { 2.000000000000000e+00 + 1.414213562373095e+00*PETSC_i,
305 432 2.000000000000000e+00 - 1.414213562373095e+00*PETSC_i};
306 432 const PetscComplex /* m == k-1 */
307 432 m1p5[5] = {-3.655694325463550e+00 + 6.543736899360086e+00*PETSC_i,
308 432 -3.655694325463550e+00 - 6.543736899360086e+00*PETSC_i,
309 -6.286704751729261e+00 ,
310 432 -5.700953298671832e+00 + 3.210265600308496e+00*PETSC_i,
311 432 -5.700953298671832e+00 - 3.210265600308496e+00*PETSC_i},
312 432 m1q5[4] = { 5.203941240131764e+00 + 5.805856841805367e+00*PETSC_i,
313 432 5.203941240131764e+00 - 5.805856841805367e+00*PETSC_i,
314 432 6.796058759868242e+00 + 1.886649260140217e+00*PETSC_i,
315 432 6.796058759868242e+00 - 1.886649260140217e+00*PETSC_i},
316 432 m1p4[4] = {-3.212806896871536e+00 + 4.773087433276636e+00*PETSC_i,
317 432 -3.212806896871536e+00 - 4.773087433276636e+00*PETSC_i,
318 432 -4.787193103128464e+00 + 1.567476416895212e+00*PETSC_i,
319 432 -4.787193103128464e+00 - 1.567476416895212e+00*PETSC_i},
320 432 m1q4[3] = { 4.675757014491557e+00 + 3.913489560603711e+00*PETSC_i,
321 432 4.675757014491557e+00 - 3.913489560603711e+00*PETSC_i,
322 5.648485971016893e+00 },
323 432 m1p3[3] = {-2.681082873627756e+00 + 3.050430199247411e+00*PETSC_i,
324 432 -2.681082873627756e+00 - 3.050430199247411e+00*PETSC_i,
325 -3.637834252744491e+00 },
326 432 m1q3[2] = { 4.000000000000000e+00 + 2.000000000000000e+00*PETSC_i,
327 432 4.000000000000000e+00 - 2.000000000000001e+00*PETSC_i},
328 432 m1p2[2] = {-2.000000000000000e+00 + 1.414213562373095e+00*PETSC_i,
329 432 -2.000000000000000e+00 - 1.414213562373095e+00*PETSC_i};
330
331
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
432 PetscFunctionBegin;
332
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
432 if (query) {
333
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
216 if (m == k+1) {
334 20 *mult = 1;
335 20 *p = k;
336 20 *q = k+1;
337
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);
338 }
339
1/2
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
196 if (m==k-1) {
340 196 *mult = 1;
341 196 *p = k;
342 196 *q = k-1;
343 }
344 } else {
345
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
216 if (m == k+1) {
346 20 *mult = PetscPowInt(-1,m);
347 20 *mult *= m;
348
1/2
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
20 if (k==4) {
349
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
100 for (i=0;i<4;i++) { p[i] = p1p4[i]; q[i] = p1q4[i]; }
350 20 q[4] = p1q4[4];
351 } else if (k==3) {
352 for (i=0;i<3;i++) { p[i] = p1p3[i]; q[i] = p1q3[i]; }
353 q[3] = p1q3[3];
354 } else if (k==2) {
355 for (i=0;i<2;i++) { p[i] = p1p2[i]; q[i] = p1q2[i]; }
356 q[2] = p1q2[2];
357 } else if (k==1) {
358 p[0] = -3;
359 for (i=0;i<2;i++) q[i] = p1q1[i];
360 }
361
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);
362 }
363
1/2
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
196 if (m==k-1) {
364 196 *mult = PetscPowInt(-1,m);
365 196 *mult /= k;
366
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
196 if (k==5) {
367
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
780 for (i=0;i<4;i++) { p[i] = m1p5[i]; q[i] = m1q5[i]; }
368 156 p[4] = m1p5[4];
369
1/2
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
40 } else if (k==4) {
370
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
160 for (i=0;i<3;i++) { p[i] = m1p4[i]; q[i] = m1q4[i]; }
371 40 p[3] = m1p4[3];
372 } else if (k==3) {
373 for (i=0;i<2;i++) { p[i] = m1p3[i]; q[i] = m1q3[i]; }
374 p[2] = m1p3[2];
375 } else if (k==2) {
376 for (i=0;i<2;i++) p[i] = m1p2[i];
377 q[0] = 3;
378 }
379 }
380 }
381
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.
72 PetscFunctionReturn(PETSC_SUCCESS);
382 }
383 #endif /* PETSC_HAVE_COMPLEX */
384
385 #if defined(PETSC_USE_COMPLEX)
386 100 static PetscErrorCode getisreal(PetscInt n,PetscComplex *a,PetscBool *result)
387 {
388 100 PetscInt i;
389
390
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
100 PetscFunctionBegin;
391 100 *result=PETSC_TRUE;
392
3/4
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
✓ Branch 2 taken 5 times.
✗ Branch 3 not taken.
296100 for (i=0;i<n&&*result;i++) {
393
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
296000 if (PetscImaginaryPartComplex(a[i])) *result=PETSC_FALSE;
394 }
395
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.
100 PetscFunctionReturn(PETSC_SUCCESS);
396 }
397 #endif
398
399 /*
400 * Matrix exponential implementation based on algorithm and matlab code by Stefan Guettel
401 * and Yuji Nakatsukasa
402 *
403 * Stefan Guettel and Yuji Nakatsukasa, "Scaled and Squared Subdiagonal Pade
404 * Approximation for the Matrix Exponential",
405 * SIAM J. Matrix Anal. Appl. 37(1):145-170, 2016.
406 * https://doi.org/10.1137/15M1027553
407 */
408 520 static PetscErrorCode FNEvaluateFunctionMat_Exp_GuettelNakatsukasa(FN fn,Mat A,Mat B)
409 {
410 #if !defined(PETSC_HAVE_COMPLEX)
411 PetscFunctionBegin;
412 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"This function requires C99 or C++ complex support");
413 #else
414 520 PetscInt i,j,n_,s,k,m,mod;
415 520 PetscBLASInt n=0,n2=0,irsize=0,rsizediv2,ipsize=0,iremainsize=0,*piv,minlen,lwork=0,one=1;
416 520 PetscReal nrm,shift=0.0;
417 #if defined(PETSC_USE_COMPLEX)
418 260 PetscReal *rwork=NULL;
419 #endif
420 520 PetscComplex *As,*RR,*RR2,*expmA,*expmA2,*Maux,*Maux2,rsize,*r,psize,*p,remainsize,*remainterm,*rootp,*rootq,mult=0.0,scale,cone=1.0,czero=0.0,*aux;
421 520 PetscScalar *Ba,*Ba2,*sMaux,*wr,*wi,expshift,sone=1.0,szero=0.0,*saux;
422 520 const PetscScalar *Aa;
423 520 PetscBool isreal,flg;
424 520 PetscBLASInt query=-1;
425 520 PetscScalar work1,*work;
426
427
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
520 PetscFunctionBegin;
428
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.
520 PetscCall(MatGetSize(A,&n_,NULL));
429
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.
520 PetscCall(PetscBLASIntCast(n_,&n));
430
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.
520 PetscCall(MatDenseGetArrayRead(A,&Aa));
431
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.
520 PetscCall(MatDenseGetArray(B,&Ba));
432 520 Ba2 = Ba;
433
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.
520 PetscCall(PetscBLASIntCast(n*n,&n2));
434
435
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.
520 PetscCall(PetscMalloc2(n2,&sMaux,n2,&Maux));
436 520 Maux2 = Maux;
437
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.
520 PetscCall(PetscOptionsGetReal(NULL,NULL,"-fn_expm_estimated_eig",&shift,&flg));
438
1/2
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
520 if (!flg) {
439
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.
520 PetscCall(PetscMalloc2(n,&wr,n,&wi));
440
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.
520 PetscCall(PetscArraycpy(sMaux,Aa,n2));
441 /* estimate rightmost eigenvalue and shift A with it */
442 #if !defined(PETSC_USE_COMPLEX)
443
13/28
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ 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 5 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.
✗ Branch 20 not taken.
✓ Branch 21 taken 1 times.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✓ Branch 25 taken 1 times.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
260 PetscCallLAPACKInfo("LAPACKgeev",LAPACKgeev_("N","N",&n,sMaux,&n,wr,wi,NULL,&n,NULL,&n,&work1,&query,&info));
444
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.
260 PetscCall(PetscBLASIntCast((PetscInt)work1,&lwork));
445
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.
260 PetscCall(PetscMalloc1(lwork,&work));
446
13/28
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ 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 5 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.
✗ Branch 20 not taken.
✓ Branch 21 taken 1 times.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✓ Branch 25 taken 1 times.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
260 PetscCallLAPACKInfo("LAPACKgeev",LAPACKgeev_("N","N",&n,sMaux,&n,wr,wi,NULL,&n,NULL,&n,work,&lwork,&info));
447
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.
260 PetscCall(PetscFree(work));
448 #else
449
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.
260 PetscCall(PetscArraycpy(Maux,Aa,n2));
450
13/28
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ 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 5 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.
✗ Branch 20 not taken.
✓ Branch 21 taken 1 times.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✓ Branch 25 taken 1 times.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
260 PetscCallLAPACKInfo("LAPACKgeev",LAPACKgeev_("N","N",&n,Maux,&n,wr,NULL,&n,NULL,&n,&work1,&query,rwork,&info));
451
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.
260 PetscCall(PetscBLASIntCast((PetscInt)PetscRealPart(work1),&lwork));
452
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.
260 PetscCall(PetscMalloc2(2*n,&rwork,lwork,&work));
453
13/28
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ 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 5 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.
✗ Branch 20 not taken.
✓ Branch 21 taken 1 times.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✓ Branch 25 taken 1 times.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
260 PetscCallLAPACKInfo("LAPACKgeev",LAPACKgeev_("N","N",&n,Maux,&n,wr,NULL,&n,NULL,&n,work,&lwork,rwork,&info));
454
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.
260 PetscCall(PetscFree2(rwork,work));
455 #endif
456
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.
520 PetscCall(PetscLogFlops(25.0*n*n*n+(n*n*n)/3.0+1.0*n*n*n));
457
458 520 shift = PetscRealPart(wr[0]);
459
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
14000 for (i=1;i<n;i++) {
460
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
13480 if (PetscRealPart(wr[i]) > shift) shift = PetscRealPart(wr[i]);
461 }
462
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.
520 PetscCall(PetscFree2(wr,wi));
463 }
464 /* shift so that largest real part is (about) 0 */
465
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.
520 PetscCall(PetscArraycpy(sMaux,Aa,n2));
466
1/2
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
520 if (shift) {
467
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
14520 for (i=0;i<n;i++) sMaux[i+i*n] -= shift;
468
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.
520 PetscCall(PetscLogFlops(1.0*n));
469 }
470 #if defined(PETSC_USE_COMPLEX)
471
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.
260 PetscCall(PetscArraycpy(Maux,Aa,n2));
472
1/2
✓ Branch 0 taken 5 times.
✗ Branch 1 not taken.
260 if (shift) {
473
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
7260 for (i=0;i<n;i++) Maux[i+i*n] -= shift;
474
3/6
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
260 PetscCall(PetscLogFlops(1.0*n));
475 }
476 #endif
477
478 /* estimate norm(A) and select the scaling factor */
479 520 nrm = LAPACKlange_("O",&n,&n,sMaux,&n,NULL);
480
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.
520 PetscCall(PetscLogFlops(1.0*n*n));
481
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.
520 PetscCall(sexpm_params(nrm,&s,&k,&m));
482
4/6
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
✓ Branch 2 taken 10 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 10 times.
✗ Branch 5 not taken.
520 if (s==0 && k==1 && m==0) { /* exp(A) = I+A to eps! */
483
1/2
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
120 if (shift) expshift = PetscExpReal(shift);
484
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
1320 for (i=0;i<n;i++) sMaux[i+i*n] += 1.0;
485
1/2
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
120 if (shift) {
486
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.
120 PetscCallBLAS("BLASscal",BLASscal_(&n2,&expshift,sMaux,&one));
487
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.
120 PetscCall(PetscLogFlops(1.0*(n+n2)));
488 } else PetscCall(PetscLogFlops(1.0*n));
489
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(PetscArraycpy(Ba,sMaux,n2));
490
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(PetscFree2(sMaux,Maux));
491
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(MatDenseRestoreArrayRead(A,&Aa));
492
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(MatDenseRestoreArray(B,&Ba));
493
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.
24 PetscFunctionReturn(PETSC_SUCCESS); /* quick return */
494 }
495
496
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.
400 PetscCall(PetscMalloc4(n2,&expmA,n2,&As,n2,&RR,n,&piv));
497 400 expmA2 = expmA; RR2 = RR;
498 /* scale matrix */
499 #if !defined(PETSC_USE_COMPLEX)
500
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
592200 for (i=0;i<n2;i++) {
501 592000 As[i] = sMaux[i];
502 }
503 #else
504
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.
200 PetscCall(PetscArraycpy(As,sMaux,n2));
505 #endif
506 400 scale = 1.0/PetscPowRealInt(2.0,s);
507
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.
400 PetscCallBLAS("BLASCOMPLEXscal",BLASCOMPLEXscal_(&n2,&scale,As,&one));
508
4/6
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
400 PetscCall(SlepcLogFlopsComplex(1.0*n2));
509
510 /* evaluate Pade approximant (partial fraction or product form) */
511
3/4
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 10 times.
400 if (fn->method==3 || !m) { /* partial fraction */
512
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
200 PetscCall(getcoeffs(k,m,&rsize,&psize,&remainsize,PETSC_TRUE));
513
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
200 PetscCall(PetscBLASIntCast((PetscInt)PetscRealPartComplex(rsize),&irsize));
514
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
200 PetscCall(PetscBLASIntCast((PetscInt)PetscRealPartComplex(psize),&ipsize));
515
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
200 PetscCall(PetscBLASIntCast((PetscInt)PetscRealPartComplex(remainsize),&iremainsize));
516
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
200 PetscCall(PetscMalloc3(irsize,&r,ipsize,&p,iremainsize,&remainterm));
517
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
200 PetscCall(getcoeffs(k,m,r,p,remainterm,PETSC_FALSE));
518
519
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.
200 PetscCall(PetscArrayzero(expmA,n2));
520 #if !defined(PETSC_USE_COMPLEX)
521 100 isreal = PETSC_TRUE;
522 #else
523
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.
100 PetscCall(getisreal(n2,Maux,&isreal));
524 #endif
525
1/2
✓ Branch 0 taken 5 times.
✗ Branch 1 not taken.
200 if (isreal) {
526 200 rsizediv2 = irsize/2;
527
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
560 for (i=0;i<rsizediv2;i++) { /* use partial fraction to get R(As) */
528
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.
360 PetscCall(PetscArraycpy(Maux,As,n2));
529
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.
360 PetscCall(PetscArrayzero(RR,n2));
530
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
12760 for (j=0;j<n;j++) {
531 12400 Maux[j+j*n] -= p[2*i];
532 12400 RR[j+j*n] = r[2*i];
533 }
534
13/28
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 8 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 10 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.
✗ Branch 20 not taken.
✓ Branch 21 taken 2 times.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✓ Branch 25 taken 2 times.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
360 PetscCallLAPACKInfo("LAPACKCOMPLEXgesv",LAPACKCOMPLEXgesv_(&n,&n,Maux,&n,piv,RR,&n,&info));
535
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
1180360 for (j=0;j<n2;j++) {
536 1180000 expmA[j] += RR[j] + PetscConj(RR[j]);
537 }
538 /* loop(n) + gesv + loop(n2) */
539
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.
360 PetscCall(SlepcLogFlopsComplex(1.0*n+(2.0*n*n*n/3.0+2.0*n*n*n)+2.0*n2));
540 }
541
542 200 mod = ipsize % 2;
543
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
200 if (mod) {
544
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.
60 PetscCall(PetscArraycpy(Maux,As,n2));
545
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.
60 PetscCall(PetscArrayzero(RR,n2));
546
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
660 for (j=0;j<n;j++) {
547 600 Maux[j+j*n] -= p[ipsize-1];
548 600 RR[j+j*n] = r[irsize-1];
549 }
550
13/28
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 8 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 10 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.
✗ Branch 20 not taken.
✓ Branch 21 taken 2 times.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✓ Branch 25 taken 2 times.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
60 PetscCallLAPACKInfo("LAPACKCOMPLEXgesv",LAPACKCOMPLEXgesv_(&n,&n,Maux,&n,piv,RR,&n,&info));
551
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
6060 for (j=0;j<n2;j++) {
552 6000 expmA[j] += RR[j];
553 }
554
4/6
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 times.
✓ Branch 5 taken 1 times.
60 PetscCall(SlepcLogFlopsComplex(1.0*n+(2.0*n*n*n/3.0+2.0*n*n*n)+1.0*n2));
555 }
556 } else { /* complex */
557 for (i=0;i<irsize;i++) { /* use partial fraction to get R(As) */
558 PetscCall(PetscArraycpy(Maux,As,n2));
559 PetscCall(PetscArrayzero(RR,n2));
560 for (j=0;j<n;j++) {
561 Maux[j+j*n] -= p[i];
562 RR[j+j*n] = r[i];
563 }
564 PetscCallLAPACKInfo("LAPACKCOMPLEXgesv",LAPACKCOMPLEXgesv_(&n,&n,Maux,&n,piv,RR,&n,&info));
565 for (j=0;j<n2;j++) {
566 expmA[j] += RR[j];
567 }
568 PetscCall(SlepcLogFlopsComplex(1.0*n+(2.0*n*n*n/3.0+2.0*n*n*n)+1.0*n2));
569 }
570 }
571
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
560 for (i=0;i<iremainsize;i++) {
572
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
360 if (!i) {
573
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.
180 PetscCall(PetscArrayzero(RR,n2));
574
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
6380 for (j=0;j<n;j++) {
575 6200 RR[j+j*n] = remainterm[iremainsize-1];
576 }
577 } else {
578
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.
180 PetscCall(PetscArraycpy(RR,As,n2));
579
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
180 for (j=1;j<i;j++) {
580 PetscCallBLAS("BLASCOMPLEXgemm",BLASCOMPLEXgemm_("N","N",&n,&n,&n,&cone,RR,&n,RR,&n,&czero,Maux,&n));
581 SlepcSwap(RR,Maux,aux);
582 PetscCall(SlepcLogFlopsComplex(2.0*n*n*n));
583 }
584
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.
180 PetscCallBLAS("BLASCOMPLEXscal",BLASCOMPLEXscal_(&n2,&remainterm[iremainsize-1-i],RR,&one));
585
3/6
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
180 PetscCall(SlepcLogFlopsComplex(1.0*n2));
586 }
587
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
1180360 for (j=0;j<n2;j++) {
588 1180000 expmA[j] += RR[j];
589 }
590
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.
360 PetscCall(SlepcLogFlopsComplex(1.0*n2));
591 }
592
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
200 PetscCall(PetscFree3(r,p,remainterm));
593 } else { /* product form, default */
594
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
200 PetscCall(getcoeffsproduct(k,m,&rsize,&psize,&mult,PETSC_TRUE));
595
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
200 PetscCall(PetscBLASIntCast((PetscInt)PetscRealPartComplex(rsize),&irsize));
596
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
200 PetscCall(PetscBLASIntCast((PetscInt)PetscRealPartComplex(psize),&ipsize));
597
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
200 PetscCall(PetscMalloc2(irsize,&rootp,ipsize,&rootq));
598
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
200 PetscCall(getcoeffsproduct(k,m,rootp,rootq,&mult,PETSC_FALSE));
599
600
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.
200 PetscCall(PetscArrayzero(expmA,n2));
601
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
6600 for (i=0;i<n;i++) { /* initialize */
602 6400 expmA[i+i*n] = 1.0;
603 }
604 200 minlen = PetscMin(irsize,ipsize);
605
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
960 for (i=0;i<minlen;i++) {
606
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.
760 PetscCall(PetscArraycpy(RR,As,n2));
607
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
25960 for (j=0;j<n;j++) {
608 25200 RR[j+j*n] -= rootp[i];
609 }
610
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.
760 PetscCallBLAS("BLASCOMPLEXgemm",BLASCOMPLEXgemm_("N","N",&n,&n,&n,&cone,RR,&n,expmA,&n,&czero,Maux,&n));
611 760 SlepcSwap(expmA,Maux,aux);
612
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.
760 PetscCall(PetscArraycpy(RR,As,n2));
613
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
25960 for (j=0;j<n;j++) {
614 25200 RR[j+j*n] -= rootq[i];
615 }
616
13/28
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 8 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 10 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.
✗ Branch 20 not taken.
✓ Branch 21 taken 2 times.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✓ Branch 25 taken 2 times.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
760 PetscCallLAPACKInfo("LAPACKCOMPLEXgesv",LAPACKCOMPLEXgesv_(&n,&n,RR,&n,piv,expmA,&n,&info));
617 /* loop(n) + gemm + loop(n) + gesv */
618
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.
760 PetscCall(SlepcLogFlopsComplex(1.0*n+(2.0*n*n*n)+1.0*n+(2.0*n*n*n/3.0+2.0*n*n*n)));
619 }
620 /* extra numerator */
621
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
380 for (i=minlen;i<irsize;i++) {
622
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.
180 PetscCall(PetscArraycpy(RR,As,n2));
623
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
6380 for (j=0;j<n;j++) {
624 6200 RR[j+j*n] -= rootp[i];
625 }
626
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.
180 PetscCallBLAS("BLASCOMPLEXgemm",BLASCOMPLEXgemm_("N","N",&n,&n,&n,&cone,RR,&n,expmA,&n,&czero,Maux,&n));
627 180 SlepcSwap(expmA,Maux,aux);
628
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.
180 PetscCall(SlepcLogFlopsComplex(1.0*n+2.0*n*n*n));
629 }
630 /* extra denominator */
631
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
220 for (i=minlen;i<ipsize;i++) {
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.
20 PetscCall(PetscArraycpy(RR,As,n2));
633
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
220 for (j=0;j<n;j++) RR[j+j*n] -= rootq[i];
634
13/28
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 8 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 10 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.
✗ Branch 20 not taken.
✓ Branch 21 taken 2 times.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✓ Branch 25 taken 2 times.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
20 PetscCallLAPACKInfo("LAPACKCOMPLEXgesv",LAPACKCOMPLEXgesv_(&n,&n,RR,&n,piv,expmA,&n,&info));
635
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.
20 PetscCall(SlepcLogFlopsComplex(1.0*n+(2.0*n*n*n/3.0+2.0*n*n*n)));
636 }
637
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.
200 PetscCallBLAS("BLASCOMPLEXscal",BLASCOMPLEXscal_(&n2,&mult,expmA,&one));
638
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.
200 PetscCall(SlepcLogFlopsComplex(1.0*n2));
639
5/6
✓ Branch 0 taken 6 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 times.
✓ Branch 5 taken 1 times.
200 PetscCall(PetscFree2(rootp,rootq));
640 }
641
642 #if !defined(PETSC_USE_COMPLEX)
643
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
592200 for (i=0;i<n2;i++) {
644 592000 Ba2[i] = PetscRealPartComplex(expmA[i]);
645 }
646 #else
647
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.
200 PetscCall(PetscArraycpy(Ba2,expmA,n2));
648 #endif
649
650 /* perform repeated squaring */
651
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
1960 for (i=0;i<s;i++) { /* final squaring */
652
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.
1560 PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n,&n,&n,&sone,Ba2,&n,Ba2,&n,&szero,sMaux,&n));
653 1560 SlepcSwap(Ba2,sMaux,saux);
654
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.
1560 PetscCall(PetscLogFlops(2.0*n*n*n));
655 }
656
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
400 if (Ba2!=Ba) {
657
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
40 PetscCall(PetscArraycpy(Ba,Ba2,n2));
658 40 sMaux = Ba2;
659 }
660
1/2
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
400 if (shift) {
661 400 expshift = PetscExpReal(shift);
662
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.
400 PetscCallBLAS("BLASscal",BLASscal_(&n2,&expshift,Ba,&one));
663
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.
400 PetscCall(PetscLogFlops(1.0*n2));
664 }
665
666 /* restore pointers */
667 400 Maux = Maux2; expmA = expmA2; RR = RR2;
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.
400 PetscCall(PetscFree2(sMaux,Maux));
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.
400 PetscCall(PetscFree4(expmA,As,RR,piv));
670
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.
400 PetscCall(MatDenseRestoreArrayRead(A,&Aa));
671
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.
400 PetscCall(MatDenseRestoreArray(B,&Ba));
672
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.
80 PetscFunctionReturn(PETSC_SUCCESS);
673 #endif
674 }
675
676 #define SMALLN 100
677
678 /*
679 * Function needed to compute optimal parameters (required workspace is 3*n*n)
680 */
681 4600 static PetscInt ell(PetscBLASInt n,PetscScalar *A,PetscReal coeff,PetscInt m,PetscScalar *work,PetscRandom rand)
682 {
683 4600 PetscScalar *Ascaled=work;
684 4600 PetscReal nrm,alpha,beta,rwork[1];
685 4600 PetscInt t;
686 4600 PetscBLASInt i,j;
687
688
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
4600 PetscFunctionBegin;
689 4600 beta = PetscPowReal(coeff,1.0/(2*m+1));
690
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
116392 for (i=0;i<n;i++)
691
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
6671592 for (j=0;j<n;j++)
692 6559800 Ascaled[i+j*n] = beta*PetscAbsScalar(A[i+j*n]);
693 4600 nrm = LAPACKlange_("O",&n,&n,A,&n,rwork);
694
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.
4600 PetscCall(PetscLogFlops(2.0*n*n));
695
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.
4600 PetscCall(SlepcNormAm(n,Ascaled,2*m+1,work+n*n,rand,&alpha));
696 4600 alpha /= nrm;
697
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
4600 t = PetscMax((PetscInt)PetscCeilReal(PetscLogReal(2.0*alpha/PETSC_MACHINE_EPSILON)/PetscLogReal(2.0)/(2*m)),0);
698
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.
904 PetscFunctionReturn(t);
699 }
700
701 /*
702 * Compute scaling parameter (s) and order of Pade approximant (m) (required workspace is 4*n*n)
703 */
704 4598 static PetscErrorCode expm_params(PetscInt n,PetscScalar **Apowers,PetscInt *s,PetscInt *m,PetscScalar *work)
705 {
706 4598 PetscScalar sfactor,sone=1.0,szero=0.0,*A=Apowers[0],*Ascaled;
707 4598 PetscReal d4,d6,d8,d10,eta1,eta3,eta4,eta5,rwork[1];
708 4598 PetscBLASInt n_=0,n2,one=1;
709 4598 PetscRandom rand;
710 4598 const PetscReal coeff[5] = { 9.92063492063492e-06, 9.94131285136576e-11, /* backward error function */
711 2.22819456055356e-16, 1.69079293431187e-22, 8.82996160201868e-36 };
712 4598 const PetscReal theta[5] = { 1.495585217958292e-002, /* m = 3 */
713 2.539398330063230e-001, /* m = 5 */
714 9.504178996162932e-001, /* m = 7 */
715 2.097847961257068e+000, /* m = 9 */
716 5.371920351148152e+000 }; /* m = 13 */
717
718
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
4598 PetscFunctionBegin;
719 4598 *s = 0;
720 4598 *m = 13;
721
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.
4598 PetscCall(PetscBLASIntCast(n,&n_));
722
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.
4598 PetscCall(PetscRandomCreate(PETSC_COMM_SELF,&rand));
723 4598 d4 = PetscPowReal(LAPACKlange_("O",&n_,&n_,Apowers[2],&n_,rwork),1.0/4.0);
724
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
4598 if (d4==0.0) { /* safeguard for the case A = 0 */
725 *m = 3;
726 goto done;
727 }
728 4598 d6 = PetscPowReal(LAPACKlange_("O",&n_,&n_,Apowers[3],&n_,rwork),1.0/6.0);
729
4/6
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
4598 PetscCall(PetscLogFlops(2.0*n*n));
730
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
4598 eta1 = PetscMax(d4,d6);
731
3/4
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
✓ Branch 2 taken 10 times.
✗ Branch 3 not taken.
4598 if (eta1<=theta[0] && !ell(n_,A,coeff[0],3,work,rand)) {
732 1115 *m = 3;
733 1115 goto done;
734 }
735
3/4
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
✓ Branch 2 taken 10 times.
✗ Branch 3 not taken.
3483 if (eta1<=theta[1] && !ell(n_,A,coeff[1],5,work,rand)) {
736 1000 *m = 5;
737 1000 goto done;
738 }
739
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
2483 if (n<SMALLN) {
740
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.
2343 PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n_,&n_,&n_,&sone,Apowers[2],&n_,Apowers[2],&n_,&szero,work,&n_));
741 2343 d8 = PetscPowReal(LAPACKlange_("O",&n_,&n_,work,&n_,rwork),1.0/8.0);
742
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.
2343 PetscCall(PetscLogFlops(2.0*n*n*n+1.0*n*n));
743 } else {
744
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.
140 PetscCall(SlepcNormAm(n_,Apowers[2],2,work,rand,&d8));
745 140 d8 = PetscPowReal(d8,1.0/8.0);
746 }
747
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
2483 eta3 = PetscMax(d6,d8);
748
3/4
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
✓ Branch 2 taken 10 times.
✗ Branch 3 not taken.
2483 if (eta3<=theta[2] && !ell(n_,A,coeff[2],7,work,rand)) {
749 114 *m = 7;
750 114 goto done;
751 }
752
4/4
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
✓ Branch 2 taken 10 times.
✓ Branch 3 taken 2 times.
2369 if (eta3<=theta[3] && !ell(n_,A,coeff[3],9,work,rand)) {
753 177 *m = 9;
754 177 goto done;
755 }
756
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
2192 if (n<SMALLN) {
757
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.
2052 PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n_,&n_,&n_,&sone,Apowers[2],&n_,Apowers[3],&n_,&szero,work,&n_));
758 2052 d10 = PetscPowReal(LAPACKlange_("O",&n_,&n_,work,&n_,rwork),1.0/10.0);
759
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.
2052 PetscCall(PetscLogFlops(2.0*n*n*n+1.0*n*n));
760 } else {
761
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.
140 PetscCall(SlepcNormAm(n_,Apowers[1],5,work,rand,&d10));
762 140 d10 = PetscPowReal(d10,1.0/10.0);
763 }
764
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
2192 eta4 = PetscMax(d8,d10);
765
1/2
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
2192 eta5 = PetscMin(eta3,eta4);
766
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
2192 *s = PetscMax((PetscInt)PetscCeilReal(PetscLogReal(eta5/theta[4])/PetscLogReal(2.0)),0);
767
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
2192 if (*s) {
768 753 Ascaled = work+3*n*n;
769 753 n2 = n_*n_;
770
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.
753 PetscCallBLAS("BLAScopy",BLAScopy_(&n2,A,&one,Ascaled,&one));
771 753 sfactor = PetscPowRealInt(2.0,-(*s));
772
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.
753 PetscCallBLAS("BLASscal",BLASscal_(&n2,&sfactor,Ascaled,&one));
773
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.
753 PetscCall(PetscLogFlops(1.0*n*n));
774 } else Ascaled = A;
775 2192 *s += ell(n_,Ascaled,coeff[4],13,work,rand);
776 4598 done:
777
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.
4598 PetscCall(PetscRandomDestroy(&rand));
778
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.
904 PetscFunctionReturn(PETSC_SUCCESS);
779 }
780
781 /*
782 * Matrix exponential implementation based on algorithm and matlab code by N. Higham and co-authors
783 *
784 * N. J. Higham, "The scaling and squaring method for the matrix exponential
785 * revisited", SIAM J. Matrix Anal. Appl. 26(4):1179-1193, 2005.
786 */
787 4550 PetscErrorCode FNEvaluateFunctionMat_Exp_Higham(FN fn,Mat A,Mat B)
788 {
789 4550 PetscBLASInt n_=0,n2,*ipiv,one=1;
790 4550 PetscInt n,m,j,s;
791 4550 PetscScalar scale,smone=-1.0,sone=1.0,stwo=2.0,szero=0.0;
792 4550 PetscScalar *Ba,*Apowers[5],*Q,*P,*W,*work,*aux;
793 4550 const PetscScalar *Aa,*c;
794 4550 const PetscScalar c3[4] = { 120, 60, 12, 1 };
795 4550 const PetscScalar c5[6] = { 30240, 15120, 3360, 420, 30, 1 };
796 4550 const PetscScalar c7[8] = { 17297280, 8648640, 1995840, 277200, 25200, 1512, 56, 1 };
797 4550 const PetscScalar c9[10] = { 17643225600.0, 8821612800.0, 2075673600, 302702400, 30270240,
798 2162160, 110880, 3960, 90, 1 };
799 4550 const PetscScalar c13[14] = { 64764752532480000.0, 32382376266240000.0, 7771770303897600.0,
800 1187353796428800.0, 129060195264000.0, 10559470521600.0,
801 670442572800.0, 33522128640.0, 1323241920.0,
802 40840800, 960960, 16380, 182, 1 };
803
804
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
4550 PetscFunctionBegin;
805
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.
4550 PetscCall(MatDenseGetArrayRead(A,&Aa));
806
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.
4550 PetscCall(MatDenseGetArray(B,&Ba));
807
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.
4550 PetscCall(MatGetSize(A,&n,NULL));
808
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.
4550 PetscCall(PetscBLASIntCast(n,&n_));
809 4550 n2 = n_*n_;
810
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.
4550 PetscCall(PetscMalloc2(8*n*n,&work,n,&ipiv));
811
812 /* Matrix powers */
813 4550 Apowers[0] = work; /* Apowers[0] = A */
814 4550 Apowers[1] = Apowers[0] + n*n; /* Apowers[1] = A^2 */
815 4550 Apowers[2] = Apowers[1] + n*n; /* Apowers[2] = A^4 */
816 4550 Apowers[3] = Apowers[2] + n*n; /* Apowers[3] = A^6 */
817 4550 Apowers[4] = Apowers[3] + n*n; /* Apowers[4] = A^8 */
818
819
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.
4550 PetscCall(PetscArraycpy(Apowers[0],Aa,n2));
820
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.
4550 PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n_,&n_,&n_,&sone,Apowers[0],&n_,Apowers[0],&n_,&szero,Apowers[1],&n_));
821
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.
4550 PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n_,&n_,&n_,&sone,Apowers[1],&n_,Apowers[1],&n_,&szero,Apowers[2],&n_));
822
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.
4550 PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n_,&n_,&n_,&sone,Apowers[1],&n_,Apowers[2],&n_,&szero,Apowers[3],&n_));
823
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.
4550 PetscCall(PetscLogFlops(6.0*n*n*n));
824
825 /* Compute scaling parameter and order of Pade approximant */
826
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.
4550 PetscCall(expm_params(n,Apowers,&s,&m,Apowers[4]));
827
828
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
4550 if (s) { /* rescale */
829
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
3925 for (j=0;j<4;j++) {
830
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
3140 scale = PetscPowRealInt(2.0,-PetscMax(2*j,1)*s);
831
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.
3140 PetscCallBLAS("BLASscal",BLASscal_(&n2,&scale,Apowers[j],&one));
832 }
833
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.
785 PetscCall(PetscLogFlops(4.0*n*n));
834 }
835
836 /* Evaluate the Pade approximant */
837
5/6
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
✓ Branch 2 taken 10 times.
✓ Branch 3 taken 10 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 10 times.
4550 switch (m) {
838 case 3: c = c3; break;
839 1000 case 5: c = c5; break;
840 98 case 7: c = c7; break;
841 177 case 9: c = c9; break;
842 2160 case 13: c = c13; break;
843 default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Wrong value of m %" PetscInt_FMT,m);
844 }
845 4550 P = Ba;
846 4550 Q = Apowers[4] + n*n;
847 4550 W = Q + n*n;
848
2/3
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
4550 switch (m) {
849 2390 case 3:
850 case 5:
851 case 7:
852 case 9:
853
12/22
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 8 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 2 times.
✓ Branch 10 taken 2 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✓ Branch 13 taken 2 times.
✓ Branch 14 taken 2 times.
✗ Branch 15 not taken.
✓ Branch 16 taken 2 times.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✓ Branch 19 taken 2 times.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
2390 if (m==9) PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n_,&n_,&n_,&sone,Apowers[1],&n_,Apowers[3],&n_,&szero,Apowers[4],&n_));
854
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.
2390 PetscCall(PetscArrayzero(P,n2));
855
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.
2390 PetscCall(PetscArrayzero(Q,n2));
856
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
25329 for (j=0;j<n;j++) {
857 22939 P[j+j*n] = c[1];
858 22939 Q[j+j*n] = c[0];
859 }
860
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
6507 for (j=m;j>=3;j-=2) {
861
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.
4117 PetscCallBLAS("BLASaxpy",BLASaxpy_(&n2,&c[j],Apowers[(j+1)/2-1],&one,P,&one));
862
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.
4117 PetscCallBLAS("BLASaxpy",BLASaxpy_(&n2,&c[j-1],Apowers[(j+1)/2-1],&one,Q,&one));
863
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.
4117 PetscCall(PetscLogFlops(4.0*n*n));
864 }
865
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.
2390 PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n_,&n_,&n_,&sone,Apowers[0],&n_,P,&n_,&szero,W,&n_));
866
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.
2390 PetscCall(PetscLogFlops(2.0*n*n*n));
867 2822 SlepcSwap(P,W,aux);
868 1918 break;
869 2160 case 13:
870 /* P = A*(Apowers[3]*(c[13]*Apowers[3] + c[11]*Apowers[2] + c[9]*Apowers[1])
871 + c[7]*Apowers[3] + c[5]*Apowers[2] + c[3]*Apowers[1] + c[1]*I) */
872
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.
2160 PetscCallBLAS("BLAScopy",BLAScopy_(&n2,Apowers[3],&one,P,&one));
873
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.
2160 PetscCallBLAS("BLASscal",BLASscal_(&n2,&c[13],P,&one));
874
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.
2160 PetscCallBLAS("BLASaxpy",BLASaxpy_(&n2,&c[11],Apowers[2],&one,P,&one));
875
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.
2160 PetscCallBLAS("BLASaxpy",BLASaxpy_(&n2,&c[9],Apowers[1],&one,P,&one));
876
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.
2160 PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n_,&n_,&n_,&sone,Apowers[3],&n_,P,&n_,&szero,W,&n_));
877
3/6
✓ Branch 0 taken 5 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
2160 PetscCall(PetscLogFlops(5.0*n*n+2.0*n*n*n));
878
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.
2160 PetscCall(PetscArrayzero(P,n2));
879
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
90473 for (j=0;j<n;j++) P[j+j*n] = c[1];
880
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.
2160 PetscCallBLAS("BLASaxpy",BLASaxpy_(&n2,&c[7],Apowers[3],&one,P,&one));
881
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.
2160 PetscCallBLAS("BLASaxpy",BLASaxpy_(&n2,&c[5],Apowers[2],&one,P,&one));
882
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.
2160 PetscCallBLAS("BLASaxpy",BLASaxpy_(&n2,&c[3],Apowers[1],&one,P,&one));
883
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.
2160 PetscCallBLAS("BLASaxpy",BLASaxpy_(&n2,&sone,P,&one,W,&one));
884
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.
2160 PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n_,&n_,&n_,&sone,Apowers[0],&n_,W,&n_,&szero,P,&n_));
885
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.
2160 PetscCall(PetscLogFlops(7.0*n*n+2.0*n*n*n));
886 /* Q = Apowers[3]*(c[12]*Apowers[3] + c[10]*Apowers[2] + c[8]*Apowers[1])
887 + c[6]*Apowers[3] + c[4]*Apowers[2] + c[2]*Apowers[1] + c[0]*I */
888
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.
2160 PetscCallBLAS("BLAScopy",BLAScopy_(&n2,Apowers[3],&one,Q,&one));
889
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.
2160 PetscCallBLAS("BLASscal",BLASscal_(&n2,&c[12],Q,&one));
890
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.
2160 PetscCallBLAS("BLASaxpy",BLASaxpy_(&n2,&c[10],Apowers[2],&one,Q,&one));
891
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.
2160 PetscCallBLAS("BLASaxpy",BLASaxpy_(&n2,&c[8],Apowers[1],&one,Q,&one));
892
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.
2160 PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n_,&n_,&n_,&sone,Apowers[3],&n_,Q,&n_,&szero,W,&n_));
893
3/6
✓ Branch 0 taken 5 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
2160 PetscCall(PetscLogFlops(5.0*n*n+2.0*n*n*n));
894
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.
2160 PetscCall(PetscArrayzero(Q,n2));
895
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
90473 for (j=0;j<n;j++) Q[j+j*n] = c[0];
896
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.
2160 PetscCallBLAS("BLASaxpy",BLASaxpy_(&n2,&c[6],Apowers[3],&one,Q,&one));
897
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.
2160 PetscCallBLAS("BLASaxpy",BLASaxpy_(&n2,&c[4],Apowers[2],&one,Q,&one));
898
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.
2160 PetscCallBLAS("BLASaxpy",BLASaxpy_(&n2,&c[2],Apowers[1],&one,Q,&one));
899
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.
2160 PetscCallBLAS("BLASaxpy",BLASaxpy_(&n2,&sone,W,&one,Q,&one));
900
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.
2160 PetscCall(PetscLogFlops(7.0*n*n));
901 break;
902 default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Wrong value of m %" PetscInt_FMT,m);
903 }
904
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.
4550 PetscCallBLAS("BLASaxpy",BLASaxpy_(&n2,&smone,P,&one,Q,&one));
905
13/28
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 8 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 10 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.
✗ Branch 20 not taken.
✓ Branch 21 taken 2 times.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✓ Branch 25 taken 2 times.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
4550 PetscCallLAPACKInfo("LAPACKgesv",LAPACKgesv_(&n_,&n_,Q,&n_,ipiv,P,&n_,&info));
906
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.
4550 PetscCallBLAS("BLASscal",BLASscal_(&n2,&stwo,P,&one));
907
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
115802 for (j=0;j<n;j++) P[j+j*n] += 1.0;
908
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.
4550 PetscCall(PetscLogFlops(2.0*n*n*n/3.0+4.0*n*n));
909
910 /* Squaring */
911
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
11170 for (j=1;j<=s;j++) {
912
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.
6620 PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n_,&n_,&n_,&sone,P,&n_,P,&n_,&szero,W,&n_));
913 6620 SlepcSwap(P,W,aux);
914 }
915
6/8
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 8 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
4550 if (P!=Ba) PetscCall(PetscArraycpy(Ba,P,n2));
916
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.
4550 PetscCall(PetscLogFlops(2.0*n*n*n*s));
917
918
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.
4550 PetscCall(PetscFree2(work,ipiv));
919
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.
4550 PetscCall(MatDenseRestoreArrayRead(A,&Aa));
920
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.
4550 PetscCall(MatDenseRestoreArray(B,&Ba));
921
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.
904 PetscFunctionReturn(PETSC_SUCCESS);
922 }
923
924 #if defined(PETSC_HAVE_CUDA)
925 #include "../src/sys/classes/fn/impls/cuda/fnutilcuda.h"
926 #include <slepccupmblas.h>
927
928 PetscErrorCode FNEvaluateFunctionMat_Exp_Pade_CUDA(FN fn,Mat A,Mat B)
929 {
930 PetscBLASInt n=0,ld,ld2,*d_ipiv,*d_info,info,one=1;
931 PetscInt m,k,sexp;
932 PetscBool odd;
933 const PetscInt p=MAX_PADE;
934 PetscReal c[MAX_PADE+1],s;
935 PetscScalar scale,smone=-1.0,sone=1.0,stwo=2.0,szero=0.0;
936 const PetscScalar *Aa;
937 PetscScalar *d_Ba,*d_As,*d_A2,*d_Q,*d_P,*d_W,*aux,**ppP,**d_ppP,**ppQ,**d_ppQ;
938 cublasHandle_t cublasv2handle;
939
940 PetscFunctionBegin;
941 PetscCall(PetscDeviceInitialize(PETSC_DEVICE_CUDA)); /* For CUDA event timers */
942 PetscCall(PetscCUBLASGetHandle(&cublasv2handle));
943 PetscCall(MatGetSize(A,&m,NULL));
944 PetscCall(PetscBLASIntCast(m,&n));
945 ld = n;
946 ld2 = ld*ld;
947 if (A==B) {
948 PetscCallCUDA(cudaMalloc((void **)&d_As,sizeof(PetscScalar)*m*m));
949 PetscCall(MatDenseCUDAGetArrayRead(A,&Aa));
950 PetscCallCUDA(cudaMemcpy(d_As,Aa,sizeof(PetscScalar)*ld2,cudaMemcpyDeviceToDevice));
951 PetscCall(MatDenseCUDARestoreArrayRead(A,&Aa));
952 } else PetscCall(MatDenseCUDAGetArrayRead(A,(const PetscScalar**)&d_As));
953 PetscCall(MatDenseCUDAGetArrayWrite(B,&d_Ba));
954
955 PetscCallCUDA(cudaMalloc((void **)&d_Q,sizeof(PetscScalar)*m*m));
956 PetscCallCUDA(cudaMalloc((void **)&d_W,sizeof(PetscScalar)*m*m));
957 PetscCallCUDA(cudaMalloc((void **)&d_A2,sizeof(PetscScalar)*m*m));
958 PetscCallCUDA(cudaMalloc((void **)&d_ipiv,sizeof(PetscBLASInt)*ld));
959 PetscCallCUDA(cudaMalloc((void **)&d_info,sizeof(PetscBLASInt)));
960 PetscCallCUDA(cudaMalloc((void **)&d_ppP,sizeof(PetscScalar*)));
961 PetscCallCUDA(cudaMalloc((void **)&d_ppQ,sizeof(PetscScalar*)));
962
963 PetscCall(PetscMalloc1(1,&ppP));
964 PetscCall(PetscMalloc1(1,&ppQ));
965
966 d_P = d_Ba;
967 PetscCall(PetscLogGpuTimeBegin());
968
969 /* Pade' coefficients */
970 c[0] = 1.0;
971 for (k=1;k<=p;k++) c[k] = c[k-1]*(p+1-k)/(k*(2*p+1-k));
972
973 /* Scaling */
974 PetscCallCUBLAS(cublasXnrm2(cublasv2handle,ld2,d_As,one,&s));
975 if (s>0.5) {
976 sexp = PetscMax(0,(int)(PetscLogReal(s)/PetscLogReal(2.0))+2);
977 scale = PetscPowRealInt(2.0,-sexp);
978 PetscCallCUBLAS(cublasXscal(cublasv2handle,ld2,&scale,d_As,one));
979 PetscCall(PetscLogGpuFlops(1.0*n*n));
980 } else sexp = 0;
981
982 /* Horner evaluation */
983 PetscCallCUBLAS(cublasXgemm(cublasv2handle,CUBLAS_OP_N,CUBLAS_OP_N,n,n,n,&sone,d_As,ld,d_As,ld,&szero,d_A2,ld));
984 PetscCall(PetscLogGpuFlops(2.0*n*n*n));
985 PetscCallCUDA(cudaMemset(d_Q,0,sizeof(PetscScalar)*ld2));
986 PetscCallCUDA(cudaMemset(d_P,0,sizeof(PetscScalar)*ld2));
987 PetscCall(set_diagonal(n,d_Q,ld,c[p]));
988 PetscCall(set_diagonal(n,d_P,ld,c[p-1]));
989
990 odd = PETSC_TRUE;
991 for (k=p-1;k>0;k--) {
992 if (odd) {
993 PetscCallCUBLAS(cublasXgemm(cublasv2handle,CUBLAS_OP_N,CUBLAS_OP_N,n,n,n,&sone,d_Q,ld,d_A2,ld,&szero,d_W,ld));
994 SlepcSwap(d_Q,d_W,aux);
995 PetscCall(shift_diagonal(n,d_Q,ld,c[k-1]));
996 odd = PETSC_FALSE;
997 } else {
998 PetscCallCUBLAS(cublasXgemm(cublasv2handle,CUBLAS_OP_N,CUBLAS_OP_N,n,n,n,&sone,d_P,ld,d_A2,ld,&szero,d_W,ld));
999 SlepcSwap(d_P,d_W,aux);
1000 PetscCall(shift_diagonal(n,d_P,ld,c[k-1]));
1001 odd = PETSC_TRUE;
1002 }
1003 PetscCall(PetscLogGpuFlops(2.0*n*n*n));
1004 }
1005 if (odd) {
1006 PetscCallCUBLAS(cublasXgemm(cublasv2handle,CUBLAS_OP_N,CUBLAS_OP_N,n,n,n,&sone,d_Q,ld,d_As,ld,&szero,d_W,ld));
1007 SlepcSwap(d_Q,d_W,aux);
1008 PetscCallCUBLAS(cublasXaxpy(cublasv2handle,ld2,&smone,d_P,one,d_Q,one));
1009
1010 ppQ[0] = d_Q;
1011 ppP[0] = d_P;
1012 PetscCallCUDA(cudaMemcpy(d_ppQ,ppQ,sizeof(PetscScalar*),cudaMemcpyHostToDevice));
1013 PetscCallCUDA(cudaMemcpy(d_ppP,ppP,sizeof(PetscScalar*),cudaMemcpyHostToDevice));
1014
1015 PetscCallCUBLAS(cublasXgetrfBatched(cublasv2handle,n,d_ppQ,ld,d_ipiv,d_info,one));
1016 PetscCallCUDA(cudaMemcpy(&info,d_info,sizeof(PetscBLASInt),cudaMemcpyDeviceToHost));
1017 PetscCheck(info>=0,PETSC_COMM_SELF,PETSC_ERR_LIB,"LAPACKgetrf: Illegal value on argument %" PetscBLASInt_FMT,PetscAbsInt(info));
1018 PetscCheck(info<=0,PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"LAPACKgetrf: Matrix is singular. U(%" PetscBLASInt_FMT ",%" PetscBLASInt_FMT ") is zero",info,info);
1019 PetscCallCUBLAS(cublasXgetrsBatched(cublasv2handle,CUBLAS_OP_N,n,n,(const PetscScalar **)d_ppQ,ld,d_ipiv,d_ppP,ld,&info,one));
1020 PetscCheck(info>=0,PETSC_COMM_SELF,PETSC_ERR_LIB,"LAPACKgetri: Illegal value on argument %" PetscBLASInt_FMT,PetscAbsInt(info));
1021 PetscCheck(info<=0,PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"LAPACKgetri: Matrix is singular. U(%" PetscBLASInt_FMT ",%" PetscBLASInt_FMT ") is zero",info,info);
1022 PetscCallCUBLAS(cublasXscal(cublasv2handle,ld2,&stwo,d_P,one));
1023 PetscCall(shift_diagonal(n,d_P,ld,sone));
1024 PetscCallCUBLAS(cublasXscal(cublasv2handle,ld2,&smone,d_P,one));
1025 } else {
1026 PetscCallCUBLAS(cublasXgemm(cublasv2handle,CUBLAS_OP_N,CUBLAS_OP_N,n,n,n,&sone,d_P,ld,d_As,ld,&szero,d_W,ld));
1027 SlepcSwap(d_P,d_W,aux);
1028 PetscCallCUBLAS(cublasXaxpy(cublasv2handle,ld2,&smone,d_P,one,d_Q,one));
1029
1030 ppQ[0] = d_Q;
1031 ppP[0] = d_P;
1032 PetscCallCUDA(cudaMemcpy(d_ppQ,ppQ,sizeof(PetscScalar*),cudaMemcpyHostToDevice));
1033 PetscCallCUDA(cudaMemcpy(d_ppP,ppP,sizeof(PetscScalar*),cudaMemcpyHostToDevice));
1034
1035 PetscCallCUBLAS(cublasXgetrfBatched(cublasv2handle,n,d_ppQ,ld,d_ipiv,d_info,one));
1036 PetscCallCUDA(cudaMemcpy(&info,d_info,sizeof(PetscBLASInt),cudaMemcpyDeviceToHost));
1037 PetscCheck(info>=0,PETSC_COMM_SELF,PETSC_ERR_LIB,"LAPACKgetrf: Illegal value on argument %" PetscBLASInt_FMT,PetscAbsInt(info));
1038 PetscCheck(info<=0,PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"LAPACKgetrf: Matrix is singular. U(%" PetscBLASInt_FMT ",%" PetscBLASInt_FMT ") is zero",info,info);
1039 PetscCallCUBLAS(cublasXgetrsBatched(cublasv2handle,CUBLAS_OP_N,n,n,(const PetscScalar **)d_ppQ,ld,d_ipiv,d_ppP,ld,&info,one));
1040 PetscCheck(info>=0,PETSC_COMM_SELF,PETSC_ERR_LIB,"LAPACKgetri: Illegal value on argument %" PetscBLASInt_FMT,PetscAbsInt(info));
1041 PetscCheck(info<=0,PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"LAPACKgetri: Matrix is singular. U(%" PetscBLASInt_FMT ",%" PetscBLASInt_FMT ") is zero",info,info);
1042 PetscCallCUBLAS(cublasXscal(cublasv2handle,ld2,&stwo,d_P,one));
1043 PetscCall(shift_diagonal(n,d_P,ld,sone));
1044 }
1045 PetscCall(PetscLogGpuFlops(2.0*n*n*n+2.0*n*n*n/3.0+4.0*n*n));
1046
1047 for (k=1;k<=sexp;k++) {
1048 PetscCallCUBLAS(cublasXgemm(cublasv2handle,CUBLAS_OP_N,CUBLAS_OP_N,n,n,n,&sone,d_P,ld,d_P,ld,&szero,d_W,ld));
1049 PetscCallCUDA(cudaMemcpy(d_P,d_W,sizeof(PetscScalar)*ld2,cudaMemcpyDeviceToDevice));
1050 }
1051 PetscCall(PetscLogGpuFlops(2.0*n*n*n*sexp));
1052
1053 PetscCall(PetscLogGpuTimeEnd());
1054 PetscCallCUDA(cudaFree(d_Q));
1055 PetscCallCUDA(cudaFree(d_W));
1056 PetscCallCUDA(cudaFree(d_A2));
1057 PetscCallCUDA(cudaFree(d_ipiv));
1058 PetscCallCUDA(cudaFree(d_info));
1059 PetscCallCUDA(cudaFree(d_ppP));
1060 PetscCallCUDA(cudaFree(d_ppQ));
1061
1062 PetscCall(PetscFree(ppP));
1063 PetscCall(PetscFree(ppQ));
1064
1065 if (d_P!=d_Ba) PetscCallCUDA(cudaMemcpy(d_Ba,d_P,sizeof(PetscScalar)*ld2,cudaMemcpyDeviceToDevice));
1066 if (A!=B) {
1067 if (s>0.5) { /* undo scaling */
1068 scale = 1.0/scale;
1069 PetscCallCUBLAS(cublasXscal(cublasv2handle,ld2,&scale,d_As,one));
1070 }
1071 PetscCall(MatDenseCUDARestoreArrayRead(A,(const PetscScalar**)&d_As));
1072 } else PetscCallCUDA(cudaFree(d_As));
1073 PetscCall(MatDenseCUDARestoreArrayWrite(B,&d_Ba));
1074 PetscFunctionReturn(PETSC_SUCCESS);
1075 }
1076
1077 #if defined(PETSC_HAVE_MAGMA)
1078 #include <slepcmagma.h>
1079
1080 16 PetscErrorCode FNEvaluateFunctionMat_Exp_Pade_CUDAm(FN fn,Mat A,Mat B)
1081 {
1082 16 PetscBLASInt n=0,ld,ld2,*piv,one=1;
1083 16 PetscInt m,k,sexp;
1084 16 PetscBool odd;
1085 16 const PetscInt p=MAX_PADE;
1086 16 PetscReal c[MAX_PADE+1],s;
1087 16 PetscScalar scale,smone=-1.0,sone=1.0,stwo=2.0,szero=0.0;
1088 16 const PetscScalar *Aa;
1089 16 PetscScalar *d_Ba,*d_As,*d_A2,*d_Q,*d_P,*d_W,*aux;
1090 16 cublasHandle_t cublasv2handle;
1091
1092 16 PetscFunctionBegin;
1093
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
16 PetscCall(PetscDeviceInitialize(PETSC_DEVICE_CUDA)); /* For CUDA event timers */
1094
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
16 PetscCall(PetscCUBLASGetHandle(&cublasv2handle));
1095
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
16 PetscCall(SlepcMagmaInit());
1096
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
16 PetscCall(MatGetSize(A,&m,NULL));
1097
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
16 PetscCall(PetscBLASIntCast(m,&n));
1098 16 ld = n;
1099 16 ld2 = ld*ld;
1100
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
16 if (A==B) {
1101
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
4 PetscCallCUDA(cudaMalloc((void **)&d_As,sizeof(PetscScalar)*m*m));
1102
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
4 PetscCall(MatDenseCUDAGetArrayRead(A,&Aa));
1103
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
4 PetscCallCUDA(cudaMemcpy(d_As,Aa,sizeof(PetscScalar)*ld2,cudaMemcpyDeviceToDevice));
1104
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
4 PetscCall(MatDenseCUDARestoreArrayRead(A,&Aa));
1105
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
12 } else PetscCall(MatDenseCUDAGetArrayRead(A,(const PetscScalar**)&d_As));
1106
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
16 PetscCall(MatDenseCUDAGetArrayWrite(B,&d_Ba));
1107
1108
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
16 PetscCallCUDA(cudaMalloc((void **)&d_Q,sizeof(PetscScalar)*m*m));
1109
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
16 PetscCallCUDA(cudaMalloc((void **)&d_W,sizeof(PetscScalar)*m*m));
1110
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
16 PetscCallCUDA(cudaMalloc((void **)&d_A2,sizeof(PetscScalar)*m*m));
1111
1112
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
16 PetscCall(PetscMalloc1(n,&piv));
1113
1114 16 d_P = d_Ba;
1115
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
16 PetscCall(PetscLogGpuTimeBegin());
1116
1117 /* Pade' coefficients */
1118 16 c[0] = 1.0;
1119
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
112 for (k=1;k<=p;k++) c[k] = c[k-1]*(p+1-k)/(k*(2*p+1-k));
1120
1121 /* Scaling */
1122
1/10
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
16 PetscCallCUBLAS(cublasXnrm2(cublasv2handle,ld2,d_As,one,&s));
1123
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
16 PetscCall(PetscLogGpuFlops(1.0*n*n));
1124
1125
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
16 if (s>0.5) {
1126
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
16 sexp = PetscMax(0,(int)(PetscLogReal(s)/PetscLogReal(2.0))+2);
1127 16 scale = PetscPowRealInt(2.0,-sexp);
1128
1/10
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
16 PetscCallCUBLAS(cublasXscal(cublasv2handle,ld2,&scale,d_As,one));
1129
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
16 PetscCall(PetscLogGpuFlops(1.0*n*n));
1130 } else sexp = 0;
1131
1132 /* Horner evaluation */
1133
1/10
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
16 PetscCallCUBLAS(cublasXgemm(cublasv2handle,CUBLAS_OP_N,CUBLAS_OP_N,n,n,n,&sone,d_As,ld,d_As,ld,&szero,d_A2,ld));
1134
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
16 PetscCall(PetscLogGpuFlops(2.0*n*n*n));
1135
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
16 PetscCallCUDA(cudaMemset(d_Q,0,sizeof(PetscScalar)*ld2));
1136
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
16 PetscCallCUDA(cudaMemset(d_P,0,sizeof(PetscScalar)*ld2));
1137
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
16 PetscCall(set_diagonal(n,d_Q,ld,c[p]));
1138
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
16 PetscCall(set_diagonal(n,d_P,ld,c[p-1]));
1139
1140 odd = PETSC_TRUE;
1141
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
96 for (k=p-1;k>0;k--) {
1142
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
80 if (odd) {
1143
1/10
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
48 PetscCallCUBLAS(cublasXgemm(cublasv2handle,CUBLAS_OP_N,CUBLAS_OP_N,n,n,n,&sone,d_Q,ld,d_A2,ld,&szero,d_W,ld));
1144 48 SlepcSwap(d_Q,d_W,aux);
1145
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
48 PetscCall(shift_diagonal(n,d_Q,ld,c[k-1]));
1146 odd = PETSC_FALSE;
1147 } else {
1148
1/10
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
32 PetscCallCUBLAS(cublasXgemm(cublasv2handle,CUBLAS_OP_N,CUBLAS_OP_N,n,n,n,&sone,d_P,ld,d_A2,ld,&szero,d_W,ld));
1149 32 SlepcSwap(d_P,d_W,aux);
1150
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
32 PetscCall(shift_diagonal(n,d_P,ld,c[k-1]));
1151 odd = PETSC_TRUE;
1152 }
1153
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
80 PetscCall(PetscLogGpuFlops(2.0*n*n*n));
1154 }
1155
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
16 if (odd) {
1156 PetscCallCUBLAS(cublasXgemm(cublasv2handle,CUBLAS_OP_N,CUBLAS_OP_N,n,n,n,&sone,d_Q,ld,d_As,ld,&szero,d_W,ld));
1157 SlepcSwap(d_Q,d_W,aux);
1158 PetscCallCUBLAS(cublasXaxpy(cublasv2handle,ld2,&smone,d_P,one,d_Q,one));
1159 PetscCallMAGMA(magma_xgesv_gpu,n,n,d_Q,ld,piv,d_P,ld);
1160 PetscCallCUBLAS(cublasXscal(cublasv2handle,ld2,&stwo,d_P,one));
1161 PetscCall(shift_diagonal(n,d_P,ld,sone));
1162 PetscCallCUBLAS(cublasXscal(cublasv2handle,ld2,&smone,d_P,one));
1163 } else {
1164
1/10
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
16 PetscCallCUBLAS(cublasXgemm(cublasv2handle,CUBLAS_OP_N,CUBLAS_OP_N,n,n,n,&sone,d_P,ld,d_As,ld,&szero,d_W,ld));
1165 16 SlepcSwap(d_P,d_W,aux);
1166
1/10
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
16 PetscCallCUBLAS(cublasXaxpy(cublasv2handle,ld2,&smone,d_P,one,d_Q,one));
1167
2/6
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
16 PetscCallMAGMA(magma_xgesv_gpu,n,n,d_Q,ld,piv,d_P,ld);
1168
1/10
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
16 PetscCallCUBLAS(cublasXscal(cublasv2handle,ld2,&stwo,d_P,one));
1169
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
16 PetscCall(shift_diagonal(n,d_P,ld,sone));
1170 }
1171
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
16 PetscCall(PetscLogGpuFlops(2.0*n*n*n+2.0*n*n*n/3.0+4.0*n*n));
1172
1173
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
96 for (k=1;k<=sexp;k++) {
1174
1/10
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
80 PetscCallCUBLAS(cublasXgemm(cublasv2handle,CUBLAS_OP_N,CUBLAS_OP_N,n,n,n,&sone,d_P,ld,d_P,ld,&szero,d_W,ld));
1175
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
80 PetscCallCUDA(cudaMemcpy(d_P,d_W,sizeof(PetscScalar)*ld2,cudaMemcpyDeviceToDevice));
1176 }
1177
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
16 PetscCall(PetscLogGpuFlops(2.0*n*n*n*sexp));
1178
1179
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
16 PetscCall(PetscLogGpuTimeEnd());
1180
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
16 PetscCallCUDA(cudaFree(d_Q));
1181
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
16 PetscCallCUDA(cudaFree(d_W));
1182
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
16 PetscCallCUDA(cudaFree(d_A2));
1183
2/4
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
16 PetscCall(PetscFree(piv));
1184
1185
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
16 if (d_P!=d_Ba) PetscCallCUDA(cudaMemcpy(d_Ba,d_P,sizeof(PetscScalar)*ld2,cudaMemcpyDeviceToDevice));
1186
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
16 if (A!=B) {
1187
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
12 if (s>0.5) { /* undo scaling */
1188 12 scale = 1.0/scale;
1189
1/10
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
12 PetscCallCUBLAS(cublasXscal(cublasv2handle,ld2,&scale,d_As,one));
1190 }
1191
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
12 PetscCall(MatDenseCUDARestoreArrayRead(A,(const PetscScalar**)&d_As));
1192
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
4 } else PetscCallCUDA(cudaFree(d_As));
1193
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
16 PetscCall(MatDenseCUDARestoreArrayWrite(B,&d_Ba));
1194 PetscFunctionReturn(PETSC_SUCCESS);
1195 }
1196
1197 /*
1198 * Matrix exponential implementation based on algorithm and matlab code by N. Higham and co-authors
1199 *
1200 * N. J. Higham, "The scaling and squaring method for the matrix exponential
1201 * revisited", SIAM J. Matrix Anal. Appl. 26(4):1179-1193, 2005.
1202 */
1203 48 PetscErrorCode FNEvaluateFunctionMat_Exp_Higham_CUDAm(FN fn,Mat A,Mat B)
1204 {
1205 48 PetscBLASInt n_=0,n2,*ipiv,one=1;
1206 48 PetscInt n,m,j,s;
1207 48 PetscScalar scale,smone=-1.0,sone=1.0,stwo=2.0,szero=0.0;
1208 48 PetscScalar *d_Ba,*Apowers[5],*d_Apowers[5],*d_Q,*d_P,*d_W,*work,*d_work,*aux;
1209 48 const PetscScalar *Aa,*c;
1210 48 const PetscScalar c3[4] = { 120, 60, 12, 1 };
1211 48 const PetscScalar c5[6] = { 30240, 15120, 3360, 420, 30, 1 };
1212 48 const PetscScalar c7[8] = { 17297280, 8648640, 1995840, 277200, 25200, 1512, 56, 1 };
1213 48 const PetscScalar c9[10] = { 17643225600, 8821612800, 2075673600, 302702400, 30270240,
1214 2162160, 110880, 3960, 90, 1 };
1215 48 const PetscScalar c13[14] = { 64764752532480000, 32382376266240000, 7771770303897600,
1216 1187353796428800, 129060195264000, 10559470521600,
1217 670442572800, 33522128640, 1323241920,
1218 40840800, 960960, 16380, 182, 1 };
1219 48 cublasHandle_t cublasv2handle;
1220
1221 48 PetscFunctionBegin;
1222
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
48 PetscCall(PetscDeviceInitialize(PETSC_DEVICE_CUDA)); /* For CUDA event timers */
1223
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
48 PetscCall(PetscCUBLASGetHandle(&cublasv2handle));
1224
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
48 PetscCall(SlepcMagmaInit());
1225
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
48 PetscCall(MatGetSize(A,&n,NULL));
1226
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
48 PetscCall(PetscBLASIntCast(n,&n_));
1227 48 n2 = n_*n_;
1228
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
48 PetscCall(PetscMalloc2(8*n*n,&work,n,&ipiv));
1229 /* Matrix powers */
1230 48 Apowers[0] = work; /* Apowers[0] = A */
1231 48 Apowers[1] = Apowers[0] + n*n; /* Apowers[1] = A^2 */
1232 48 Apowers[2] = Apowers[1] + n*n; /* Apowers[2] = A^4 */
1233 48 Apowers[3] = Apowers[2] + n*n; /* Apowers[3] = A^6 */
1234 48 Apowers[4] = Apowers[3] + n*n; /* Apowers[4] = A^8 */
1235
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
48 if (A==B) {
1236
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
4 PetscCallCUDA(cudaMalloc((void**)&d_work,7*n*n*sizeof(PetscScalar)));
1237 4 d_Apowers[0] = d_work; /* d_Apowers[0] = A */
1238 4 d_Apowers[1] = d_Apowers[0] + n*n; /* d_Apowers[1] = A^2 */
1239
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
4 PetscCall(MatDenseCUDAGetArrayRead(A,&Aa));
1240
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
4 PetscCallCUDA(cudaMemcpy(d_Apowers[0],Aa,n2*sizeof(PetscScalar),cudaMemcpyDeviceToDevice));
1241
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
4 PetscCall(MatDenseCUDARestoreArrayRead(A,&Aa));
1242 } else {
1243
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
44 PetscCallCUDA(cudaMalloc((void**)&d_work,6*n*n*sizeof(PetscScalar)));
1244
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
44 PetscCall(MatDenseCUDAGetArrayRead(A,(const PetscScalar**)&d_Apowers[0]));
1245 44 d_Apowers[1] = d_work; /* d_Apowers[1] = A^2 */
1246 }
1247
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
48 PetscCall(MatDenseCUDAGetArrayWrite(B,&d_Ba));
1248 48 d_Apowers[2] = d_Apowers[1] + n*n; /* d_Apowers[2] = A^4 */
1249 48 d_Apowers[3] = d_Apowers[2] + n*n; /* d_Apowers[3] = A^6 */
1250 48 d_Apowers[4] = d_Apowers[3] + n*n; /* d_Apowers[4] = A^8 */
1251 48 d_Q = d_Apowers[4] + n*n;
1252 48 d_W = d_Q + n*n;
1253
1254
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
48 PetscCall(PetscLogGpuTimeBegin());
1255
1256
1/10
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
48 PetscCallCUBLAS(cublasXgemm(cublasv2handle,CUBLAS_OP_N,CUBLAS_OP_N,n_,n_,n_,&sone,d_Apowers[0],n_,d_Apowers[0],n_,&szero,d_Apowers[1],n_));
1257
1/10
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
48 PetscCallCUBLAS(cublasXgemm(cublasv2handle,CUBLAS_OP_N,CUBLAS_OP_N,n_,n_,n_,&sone,d_Apowers[1],n_,d_Apowers[1],n_,&szero,d_Apowers[2],n_));
1258
1/10
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
48 PetscCallCUBLAS(cublasXgemm(cublasv2handle,CUBLAS_OP_N,CUBLAS_OP_N,n_,n_,n_,&sone,d_Apowers[1],n_,d_Apowers[2],n_,&szero,d_Apowers[3],n_));
1259
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
48 PetscCall(PetscLogGpuFlops(6.0*n*n*n));
1260
1261
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
48 PetscCallCUDA(cudaMemcpy(Apowers[0],d_Apowers[0],n2*sizeof(PetscScalar),cudaMemcpyDeviceToHost));
1262
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
48 PetscCallCUDA(cudaMemcpy(Apowers[1],d_Apowers[1],3*n2*sizeof(PetscScalar),cudaMemcpyDeviceToHost));
1263
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
48 PetscCall(PetscLogGpuToCpu(4*n2*sizeof(PetscScalar)));
1264 /* Compute scaling parameter and order of Pade approximant */
1265
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
48 PetscCall(expm_params(n,Apowers,&s,&m,Apowers[4]));
1266
1267
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
48 if (s) { /* rescale */
1268
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
80 for (j=0;j<4;j++) {
1269
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
64 scale = PetscPowRealInt(2.0,-PetscMax(2*j,1)*s);
1270
1/10
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
64 PetscCallCUBLAS(cublasXscal(cublasv2handle,n2,&scale,d_Apowers[j],one));
1271 }
1272
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
16 PetscCall(PetscLogGpuFlops(4.0*n*n));
1273 }
1274
1275 /* Evaluate the Pade approximant */
1276
2/6
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
48 switch (m) {
1277 case 3: c = c3; break;
1278 case 5: c = c5; break;
1279 16 case 7: c = c7; break;
1280 case 9: c = c9; break;
1281 32 case 13: c = c13; break;
1282 default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Wrong value of m %" PetscInt_FMT,m);
1283 }
1284 48 d_P = d_Ba;
1285
2/3
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
48 switch (m) {
1286 16 case 3:
1287 case 5:
1288 case 7:
1289 case 9:
1290
1/12
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
16 if (m==9) PetscCallCUBLAS(cublasXgemm(cublasv2handle,CUBLAS_OP_N,CUBLAS_OP_N,n_,n_,n_,&sone,d_Apowers[1],n_,d_Apowers[3],n_,&szero,d_Apowers[4],n_));
1291
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
16 PetscCallCUDA(cudaMemset(d_P,0,sizeof(PetscScalar)*n2));
1292
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
16 PetscCallCUDA(cudaMemset(d_Q,0,sizeof(PetscScalar)*n2));
1293
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
16 PetscCall(set_diagonal(n,d_P,n,c[1]));
1294
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
16 PetscCall(set_diagonal(n,d_Q,n,c[0]));
1295
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
64 for (j=m;j>=3;j-=2) {
1296
1/10
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
48 PetscCallCUBLAS(cublasXaxpy(cublasv2handle,n2,&c[j],d_Apowers[(j+1)/2-1],one,d_P,one));
1297
1/10
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
48 PetscCallCUBLAS(cublasXaxpy(cublasv2handle,n2,&c[j-1],d_Apowers[(j+1)/2-1],one,d_Q,one));
1298
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
48 PetscCall(PetscLogGpuFlops(4.0*n*n));
1299 }
1300
1/10
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
16 PetscCallCUBLAS(cublasXgemm(cublasv2handle,CUBLAS_OP_N,CUBLAS_OP_N,n_,n_,n_,&sone,d_Apowers[0],n_,d_P,n_,&szero,d_W,n_));
1301
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
16 PetscCall(PetscLogGpuFlops(2.0*n*n*n));
1302 48 SlepcSwap(d_P,d_W,aux);
1303 break;
1304 32 case 13:
1305 /* P = A*(Apowers[3]*(c[13]*Apowers[3] + c[11]*Apowers[2] + c[9]*Apowers[1])
1306 + c[7]*Apowers[3] + c[5]*Apowers[2] + c[3]*Apowers[1] + c[1]*I) */
1307
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
32 PetscCallCUDA(cudaMemcpy(d_P,d_Apowers[3],n2*sizeof(PetscScalar),cudaMemcpyDeviceToDevice));
1308
1/10
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
32 PetscCallCUBLAS(cublasXscal(cublasv2handle,n2,&c[13],d_P,one));
1309
1/10
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
32 PetscCallCUBLAS(cublasXaxpy(cublasv2handle,n2,&c[11],d_Apowers[2],one,d_P,one));
1310
1/10
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
32 PetscCallCUBLAS(cublasXaxpy(cublasv2handle,n2,&c[9],d_Apowers[1],one,d_P,one));
1311
1/10
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
32 PetscCallCUBLAS(cublasXgemm(cublasv2handle,CUBLAS_OP_N,CUBLAS_OP_N,n_,n_,n_,&sone,d_Apowers[3],n_,d_P,n_,&szero,d_W,n_));
1312
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
32 PetscCall(PetscLogGpuFlops(5.0*n*n+2.0*n*n*n));
1313
1314
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
32 PetscCallCUDA(cudaMemset(d_P,0,sizeof(PetscScalar)*n2));
1315
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
32 PetscCall(set_diagonal(n,d_P,n,c[1]));
1316
1/10
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
32 PetscCallCUBLAS(cublasXaxpy(cublasv2handle,n2,&c[7],d_Apowers[3],one,d_P,one));
1317
1/10
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
32 PetscCallCUBLAS(cublasXaxpy(cublasv2handle,n2,&c[5],d_Apowers[2],one,d_P,one));
1318
1/10
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
32 PetscCallCUBLAS(cublasXaxpy(cublasv2handle,n2,&c[3],d_Apowers[1],one,d_P,one));
1319
1/10
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
32 PetscCallCUBLAS(cublasXaxpy(cublasv2handle,n2,&sone,d_P,one,d_W,one));
1320
1/10
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
32 PetscCallCUBLAS(cublasXgemm(cublasv2handle,CUBLAS_OP_N,CUBLAS_OP_N,n_,n_,n_,&sone,d_Apowers[0],n_,d_W,n_,&szero,d_P,n_));
1321
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
32 PetscCall(PetscLogGpuFlops(7.0*n*n+2.0*n*n*n));
1322 /* Q = Apowers[3]*(c[12]*Apowers[3] + c[10]*Apowers[2] + c[8]*Apowers[1])
1323 + c[6]*Apowers[3] + c[4]*Apowers[2] + c[2]*Apowers[1] + c[0]*I */
1324
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
32 PetscCallCUDA(cudaMemcpy(d_Q,d_Apowers[3],n2*sizeof(PetscScalar),cudaMemcpyDeviceToDevice));
1325
1/10
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
32 PetscCallCUBLAS(cublasXscal(cublasv2handle,n2,&c[12],d_Q,one));
1326
1/10
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
32 PetscCallCUBLAS(cublasXaxpy(cublasv2handle,n2,&c[10],d_Apowers[2],one,d_Q,one));
1327
1/10
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
32 PetscCallCUBLAS(cublasXaxpy(cublasv2handle,n2,&c[8],d_Apowers[1],one,d_Q,one));
1328
1/10
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
32 PetscCallCUBLAS(cublasXgemm(cublasv2handle,CUBLAS_OP_N,CUBLAS_OP_N,n_,n_,n_,&sone,d_Apowers[3],n_,d_Q,n_,&szero,d_W,n_));
1329
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
32 PetscCall(PetscLogGpuFlops(5.0*n*n+2.0*n*n*n));
1330
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
32 PetscCallCUDA(cudaMemset(d_Q,0,sizeof(PetscScalar)*n2));
1331
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
32 PetscCall(set_diagonal(n,d_Q,n,c[0]));
1332
1/10
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
32 PetscCallCUBLAS(cublasXaxpy(cublasv2handle,n2,&c[6],d_Apowers[3],one,d_Q,one));
1333
1/10
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
32 PetscCallCUBLAS(cublasXaxpy(cublasv2handle,n2,&c[4],d_Apowers[2],one,d_Q,one));
1334
1/10
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
32 PetscCallCUBLAS(cublasXaxpy(cublasv2handle,n2,&c[2],d_Apowers[1],one,d_Q,one));
1335
1/10
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
32 PetscCallCUBLAS(cublasXaxpy(cublasv2handle,n2,&sone,d_W,one,d_Q,one));
1336
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
32 PetscCall(PetscLogGpuFlops(7.0*n*n));
1337 break;
1338 default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Wrong value of m %" PetscInt_FMT,m);
1339 }
1340
1/10
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
48 PetscCallCUBLAS(cublasXaxpy(cublasv2handle,n2,&smone,d_P,one,d_Q,one));
1341
1342
2/6
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
48 PetscCallMAGMA(magma_xgesv_gpu,n_,n_,d_Q,n_,ipiv,d_P,n_);
1343
1344
1/10
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
48 PetscCallCUBLAS(cublasXscal(cublasv2handle,n2,&stwo,d_P,one));
1345
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
48 PetscCall(shift_diagonal(n,d_P,n,sone));
1346
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
48 PetscCall(PetscLogGpuFlops(2.0*n*n*n/3.0+4.0*n*n));
1347
1348 /* Squaring */
1349
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
64 for (j=1;j<=s;j++) {
1350
1/10
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
16 PetscCallCUBLAS(cublasXgemm(cublasv2handle,CUBLAS_OP_N,CUBLAS_OP_N,n_,n_,n_,&sone,d_P,n_,d_P,n_,&szero,d_W,n_));
1351 16 SlepcSwap(d_P,d_W,aux);
1352 }
1353
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
48 PetscCall(PetscLogGpuFlops(2.0*n*n*n*s));
1354
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
48 PetscCall(PetscLogGpuTimeEnd());
1355
1356
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
48 PetscCall(PetscFree2(work,ipiv));
1357
3/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
48 if (d_P!=d_Ba) PetscCallCUDA(cudaMemcpy(d_Ba,d_P,n2*sizeof(PetscScalar),cudaMemcpyDeviceToDevice));
1358
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
48 if (A!=B) {
1359
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
44 if (s>0.5) { /* undo scaling */
1360 12 scale = 1.0/PetscPowRealInt(2.0,-s);
1361
1/10
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
12 PetscCallCUBLAS(cublasXscal(cublasv2handle,n2,&scale,d_Apowers[0],one));
1362 }
1363
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
44 PetscCall(MatDenseCUDARestoreArrayRead(A,(const PetscScalar**)&d_Apowers[0]));
1364 }
1365
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
48 PetscCall(MatDenseCUDARestoreArrayWrite(B,&d_Ba));
1366
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
48 PetscCallCUDA(cudaFree(d_work));
1367 PetscFunctionReturn(PETSC_SUCCESS);
1368 }
1369
1370 /*
1371 * Matrix exponential implementation based on algorithm and matlab code by Stefan Guettel
1372 * and Yuji Nakatsukasa
1373 *
1374 * Stefan Guettel and Yuji Nakatsukasa, "Scaled and Squared Subdiagonal Pade'
1375 * Approximation for the Matrix Exponential",
1376 * SIAM J. Matrix Anal. Appl. 37(1):145-170, 2016.
1377 * https://doi.org/10.1137/15M1027553
1378 */
1379 32 PetscErrorCode FNEvaluateFunctionMat_Exp_GuettelNakatsukasa_CUDAm(FN fn,Mat A,Mat B)
1380 {
1381 32 PetscInt i,j,n_,s,k,m,mod;
1382 32 PetscBLASInt n=0,n2=0,irsize=0,rsizediv2,ipsize=0,iremainsize=0,query=-1,*piv,minlen,lwork=0,one=1;
1383 32 PetscReal nrm,shift=0.0,rone=1.0,rzero=0.0;
1384 #if defined(PETSC_USE_COMPLEX)
1385 16 PetscReal *rwork=NULL;
1386 #endif
1387 32 PetscComplex *d_As,*d_RR,*d_RR2,*d_expmA,*d_expmA2,*d_Maux,*d_Maux2,rsize,*r,psize,*p,remainsize,*remainterm,*rootp,*rootq,mult=0.0,scale,cone=1.0,czero=0.0,*aux;
1388 32 PetscScalar *d_Aa,*d_Ba,*d_Ba2,*Maux,*d_sMaux,*wr,*wi,expshift,sone=1.0,szero=0.0,*work,work1,*saux;
1389 32 const PetscScalar *Aa;
1390 32 PetscBool isreal,*d_isreal,flg;
1391 32 cublasHandle_t cublasv2handle;
1392
1393 32 PetscFunctionBegin;
1394
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
32 PetscCall(PetscDeviceInitialize(PETSC_DEVICE_CUDA)); /* For CUDA event timers */
1395
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
32 PetscCall(PetscCUBLASGetHandle(&cublasv2handle));
1396
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
32 PetscCall(SlepcMagmaInit());
1397
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
32 PetscCall(MatGetSize(A,&n_,NULL));
1398
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
32 PetscCall(PetscBLASIntCast(n_,&n));
1399
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
32 PetscCall(PetscBLASIntCast(n*n,&n2));
1400
1401
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
32 if (A==B) {
1402
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
8 PetscCallCUDA(cudaMalloc((void **)&d_Aa,sizeof(PetscScalar)*n2));
1403
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
8 PetscCall(MatDenseCUDAGetArrayRead(A,&Aa));
1404
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
8 PetscCallCUDA(cudaMemcpy(d_Aa,Aa,sizeof(PetscScalar)*n2,cudaMemcpyDeviceToDevice));
1405
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
8 PetscCall(MatDenseCUDARestoreArrayRead(A,&Aa));
1406
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
24 } else PetscCall(MatDenseCUDAGetArrayRead(A,(const PetscScalar**)&d_Aa));
1407
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
32 PetscCall(MatDenseCUDAGetArrayWrite(B,&d_Ba));
1408 32 d_Ba2 = d_Ba;
1409
1410
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
32 PetscCallCUDA(cudaMalloc((void **)&d_isreal,sizeof(PetscBool)));
1411
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
32 PetscCallCUDA(cudaMalloc((void **)&d_sMaux,sizeof(PetscScalar)*n2));
1412
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
32 PetscCallCUDA(cudaMalloc((void **)&d_Maux,sizeof(PetscComplex)*n2));
1413
1414
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
32 PetscCall(PetscLogGpuTimeBegin());
1415 32 d_Maux2 = d_Maux;
1416
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
32 PetscCall(PetscOptionsGetReal(NULL,NULL,"-fn_expm_estimated_eig",&shift,&flg));
1417
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
32 if (!flg) {
1418
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
32 PetscCall(PetscMalloc2(n,&wr,n,&wi));
1419 /* estimate rightmost eigenvalue and shift A with it */
1420
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
32 PetscCall(PetscMalloc1(n2,&Maux));
1421
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
32 PetscCall(MatDenseGetArrayRead(A,&Aa));
1422
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
32 PetscCall(PetscArraycpy(Maux,Aa,n2));
1423
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
32 PetscCall(MatDenseRestoreArrayRead(A,&Aa));
1424 #if !defined(PETSC_USE_COMPLEX)
1425
2/6
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
16 PetscCallMAGMA(magma_xgeev,MagmaNoVec,MagmaNoVec,n,Maux,n,wr,wi,NULL,n,NULL,n,&work1,query);
1426
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
16 PetscCall(PetscBLASIntCast((PetscInt)PetscRealPart(work1),&lwork));
1427
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
16 PetscCall(PetscMalloc1(lwork,&work));
1428
2/6
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
16 PetscCallMAGMA(magma_xgeev,MagmaNoVec,MagmaNoVec,n,Maux,n,wr,wi,NULL,n,NULL,n,work,lwork);
1429
2/4
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
16 PetscCall(PetscFree(work));
1430 #else
1431
2/6
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
16 PetscCallMAGMA(magma_xgeev,MagmaNoVec,MagmaNoVec,n,Maux,n,wr,NULL,n,NULL,n,&work1,query,rwork);
1432
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
16 PetscCall(PetscBLASIntCast((PetscInt)PetscRealPart(work1),&lwork));
1433
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
16 PetscCall(PetscMalloc2(2*n,&rwork,lwork,&work));
1434
2/6
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
16 PetscCallMAGMA(magma_xgeev,MagmaNoVec,MagmaNoVec,n,Maux,n,wr,NULL,n,NULL,n,work,lwork,rwork);
1435
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
16 PetscCall(PetscFree2(rwork,work));
1436 #endif
1437
2/4
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
32 PetscCall(PetscFree(Maux));
1438
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
32 PetscCall(PetscLogGpuFlops(25.0*n*n*n+(n*n*n)/3.0+1.0*n*n*n));
1439
1440 32 shift = PetscRealPart(wr[0]);
1441
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
320 for (i=1;i<n;i++) {
1442
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 2 times.
288 if (PetscRealPart(wr[i]) > shift) shift = PetscRealPart(wr[i]);
1443 }
1444
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
32 PetscCall(PetscFree2(wr,wi));
1445 }
1446 /* shift so that largest real part is (about) 0 */
1447
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
32 PetscCallCUDA(cudaMemcpy(d_sMaux,d_Aa,sizeof(PetscScalar)*n2,cudaMemcpyDeviceToDevice));
1448
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
32 if (shift) {
1449
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
32 PetscCall(shift_diagonal(n,d_sMaux,n,-shift));
1450
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
32 PetscCall(PetscLogGpuFlops(1.0*n));
1451 }
1452 #if defined(PETSC_USE_COMPLEX)
1453
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
16 PetscCallCUDA(cudaMemcpy(d_Maux,d_Aa,sizeof(PetscScalar)*n2,cudaMemcpyDeviceToDevice));
1454
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
16 if (shift) {
1455
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
16 PetscCall(shift_diagonal(n,d_Maux,n,-shift));
1456
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
16 PetscCall(PetscLogGpuFlops(1.0*n));
1457 }
1458 #endif
1459
3/4
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
32 if (A!=B) PetscCall(MatDenseCUDARestoreArrayRead(A,(const PetscScalar**)&d_Aa));
1460
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
8 else PetscCallCUDA(cudaFree(d_Aa));
1461
1462 /* estimate norm(A) and select the scaling factor */
1463
1/10
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
32 PetscCallCUBLAS(cublasXnrm2(cublasv2handle,n2,d_sMaux,one,&nrm));
1464
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
32 PetscCall(PetscLogGpuFlops(2.0*n*n));
1465
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
32 PetscCall(sexpm_params(nrm,&s,&k,&m));
1466
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
32 if (s==0 && k==1 && m==0) { /* exp(A) = I+A to eps! */
1467 if (shift) expshift = PetscExpReal(shift);
1468 PetscCall(shift_Cdiagonal(n,d_Maux,n,rone,rzero));
1469 if (shift) {
1470 PetscCallCUBLAS(cublasXscal(cublasv2handle,n2,&expshift,d_sMaux,one));
1471 PetscCall(PetscLogGpuFlops(1.0*(n+n2)));
1472 } else PetscCall(PetscLogGpuFlops(1.0*n));
1473 PetscCallCUDA(cudaMemcpy(d_Ba,d_sMaux,sizeof(PetscScalar)*n2,cudaMemcpyDeviceToDevice));
1474 PetscCallCUDA(cudaFree(d_isreal));
1475 PetscCallCUDA(cudaFree(d_sMaux));
1476 PetscCallCUDA(cudaFree(d_Maux));
1477 PetscCall(MatDenseCUDARestoreArrayWrite(B,&d_Ba));
1478 PetscFunctionReturn(PETSC_SUCCESS); /* quick return */
1479 }
1480
1481
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
32 PetscCallCUDA(cudaMalloc((void **)&d_expmA,sizeof(PetscComplex)*n2));
1482
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
32 PetscCallCUDA(cudaMalloc((void **)&d_As,sizeof(PetscComplex)*n2));
1483
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
32 PetscCallCUDA(cudaMalloc((void **)&d_RR,sizeof(PetscComplex)*n2));
1484 32 d_expmA2 = d_expmA; d_RR2 = d_RR;
1485
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
32 PetscCall(PetscMalloc1(n,&piv));
1486 /* scale matrix */
1487 #if !defined(PETSC_USE_COMPLEX)
1488
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
16 PetscCall(copy_array2D_S2C(n,n,d_As,n,d_sMaux,n));
1489 #else
1490
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
16 PetscCallCUDA(cudaMemcpy(d_As,d_sMaux,sizeof(PetscScalar)*n2,cudaMemcpyDeviceToDevice));
1491 #endif
1492 32 scale = 1.0/PetscPowRealInt(2.0,s);
1493
1/10
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
32 PetscCallCUBLAS(cublasXCscal(cublasv2handle,n2,(const cuComplex *)&scale,(cuComplex *)d_As,one));
1494
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
32 PetscCall(SlepcLogGpuFlopsComplex(1.0*n2));
1495
1496 /* evaluate Pade approximant (partial fraction or product form) */
1497
3/4
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
32 if (fn->method==3 || !m) { /* partial fraction */
1498
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
16 PetscCall(getcoeffs(k,m,&rsize,&psize,&remainsize,PETSC_TRUE));
1499
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
16 PetscCall(PetscBLASIntCast((PetscInt)PetscRealPartComplex(rsize),&irsize));
1500
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
16 PetscCall(PetscBLASIntCast((PetscInt)PetscRealPartComplex(psize),&ipsize));
1501
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
16 PetscCall(PetscBLASIntCast((PetscInt)PetscRealPartComplex(remainsize),&iremainsize));
1502
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
16 PetscCall(PetscMalloc3(irsize,&r,ipsize,&p,iremainsize,&remainterm));
1503
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
16 PetscCall(getcoeffs(k,m,r,p,remainterm,PETSC_FALSE));
1504
1505
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
16 PetscCallCUDA(cudaMemset(d_expmA,0,sizeof(PetscComplex)*n2));
1506 #if !defined(PETSC_USE_COMPLEX)
1507 8 isreal = PETSC_TRUE;
1508 #else
1509
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
8 PetscCall(getisreal_array2D(n,n,d_Maux,n,d_isreal));
1510
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
8 PetscCallCUDA(cudaMemcpy(&isreal,d_isreal,sizeof(PetscBool),cudaMemcpyDeviceToHost));
1511 #endif
1512
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
16 if (isreal) {
1513 16 rsizediv2 = irsize/2;
1514
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
48 for (i=0;i<rsizediv2;i++) { /* use partial fraction to get R(As) */
1515
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
32 PetscCallCUDA(cudaMemcpy(d_Maux,d_As,sizeof(PetscComplex)*n2,cudaMemcpyDeviceToDevice));
1516
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
32 PetscCallCUDA(cudaMemset(d_RR,0,sizeof(PetscComplex)*n2));
1517
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
32 PetscCall(shift_Cdiagonal(n,d_Maux,n,-PetscRealPartComplex(p[2*i]),-PetscImaginaryPartComplex(p[2*i])));
1518
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
32 PetscCall(set_Cdiagonal(n,d_RR,n,PetscRealPartComplex(r[2*i]),PetscImaginaryPartComplex(r[2*i])));
1519
2/6
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
32 PetscCallMAGMA(magma_Cgesv_gpu,n,n,d_Maux,n,piv,d_RR,n);
1520
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
32 PetscCall(add_array2D_Conj(n,n,d_RR,n));
1521
1/10
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
32 PetscCallCUBLAS(cublasXCaxpy(cublasv2handle,n2,&cone,d_RR,one,d_expmA,one));
1522 /* shift(n) + gesv + axpy(n2) */
1523
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
32 PetscCall(SlepcLogGpuFlopsComplex(1.0*n+(2.0*n*n*n/3.0+2.0*n*n*n)+2.0*n2));
1524 }
1525
1526 16 mod = ipsize % 2;
1527
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 1 times.
16 if (mod) {
1528 PetscCallCUDA(cudaMemcpy(d_Maux,d_As,sizeof(PetscComplex)*n2,cudaMemcpyDeviceToDevice));
1529 PetscCallCUDA(cudaMemset(d_RR,0,sizeof(PetscComplex)*n2));
1530 PetscCall(shift_Cdiagonal(n,d_Maux,n,-PetscRealPartComplex(p[ipsize-1]),-PetscImaginaryPartComplex(p[ipsize-1])));
1531 PetscCall(set_Cdiagonal(n,d_RR,n,PetscRealPartComplex(r[irsize-1]),PetscImaginaryPartComplex(r[irsize-1])));
1532 PetscCallMAGMA(magma_Cgesv_gpu,n,n,d_Maux,n,piv,d_RR,n);
1533 PetscCallCUBLAS(cublasXCaxpy(cublasv2handle,n2,&cone,d_RR,one,d_expmA,one));
1534 PetscCall(SlepcLogGpuFlopsComplex(1.0*n+(2.0*n*n*n/3.0+2.0*n*n*n)+1.0*n2));
1535 }
1536 } else { /* complex */
1537 for (i=0;i<irsize;i++) { /* use partial fraction to get R(As) */
1538 PetscCallCUDA(cudaMemcpy(d_Maux,d_As,sizeof(PetscComplex)*n2,cudaMemcpyDeviceToDevice));
1539 PetscCallCUDA(cudaMemset(d_RR,0,sizeof(PetscComplex)*n2));
1540 PetscCall(shift_Cdiagonal(n,d_Maux,n,-PetscRealPartComplex(p[i]),-PetscImaginaryPartComplex(p[i])));
1541 PetscCall(set_Cdiagonal(n,d_RR,n,PetscRealPartComplex(r[i]),PetscImaginaryPartComplex(r[i])));
1542 PetscCallMAGMA(magma_Cgesv_gpu,n,n,d_Maux,n,piv,d_RR,n);
1543 PetscCallCUBLAS(cublasXCaxpy(cublasv2handle,n2,&cone,d_RR,one,d_expmA,one));
1544 PetscCall(SlepcLogGpuFlopsComplex(1.0*n+(2.0*n*n*n/3.0+2.0*n*n*n)+1.0*n2));
1545 }
1546 }
1547
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
48 for (i=0;i<iremainsize;i++) {
1548
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
32 if (!i) {
1549
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
16 PetscCallCUDA(cudaMemset(d_RR,0,sizeof(PetscComplex)*n2));
1550
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
16 PetscCall(set_Cdiagonal(n,d_RR,n,PetscRealPartComplex(remainterm[iremainsize-1]),PetscImaginaryPartComplex(remainterm[iremainsize-1])));
1551 } else {
1552
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
16 PetscCallCUDA(cudaMemcpy(d_RR,d_As,sizeof(PetscComplex)*n2,cudaMemcpyDeviceToDevice));
1553
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
16 for (j=1;j<i;j++) {
1554 PetscCallCUBLAS(cublasXCgemm(cublasv2handle,CUBLAS_OP_N,CUBLAS_OP_N,n,n,n,&cone,d_RR,n,d_RR,n,&czero,d_Maux,n));
1555 SlepcSwap(d_RR,d_Maux,aux);
1556 PetscCall(SlepcLogGpuFlopsComplex(2.0*n*n*n));
1557 }
1558
1/10
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
16 PetscCallCUBLAS(cublasXCscal(cublasv2handle,n2,&remainterm[iremainsize-1-i],d_RR,one));
1559
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
16 PetscCall(SlepcLogGpuFlopsComplex(1.0*n2));
1560 }
1561
1/10
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
32 PetscCallCUBLAS(cublasXCaxpy(cublasv2handle,n2,&cone,d_RR,one,d_expmA,one));
1562
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
32 PetscCall(SlepcLogGpuFlopsComplex(1.0*n2));
1563 }
1564
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
16 PetscCall(PetscFree3(r,p,remainterm));
1565 } else { /* product form, default */
1566
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
16 PetscCall(getcoeffsproduct(k,m,&rsize,&psize,&mult,PETSC_TRUE));
1567
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
16 PetscCall(PetscBLASIntCast((PetscInt)PetscRealPartComplex(rsize),&irsize));
1568
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
16 PetscCall(PetscBLASIntCast((PetscInt)PetscRealPartComplex(psize),&ipsize));
1569
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
16 PetscCall(PetscMalloc2(irsize,&rootp,ipsize,&rootq));
1570
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
16 PetscCall(getcoeffsproduct(k,m,rootp,rootq,&mult,PETSC_FALSE));
1571
1572
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
16 PetscCallCUDA(cudaMemset(d_expmA,0,sizeof(PetscComplex)*n2));
1573
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
16 PetscCall(set_Cdiagonal(n,d_expmA,n,rone,rzero)); /* initialize */
1574 16 minlen = PetscMin(irsize,ipsize);
1575
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
80 for (i=0;i<minlen;i++) {
1576
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
64 PetscCallCUDA(cudaMemcpy(d_RR,d_As,sizeof(PetscComplex)*n2,cudaMemcpyDeviceToDevice));
1577
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
64 PetscCall(shift_Cdiagonal(n,d_RR,n,-PetscRealPartComplex(rootp[i]),-PetscImaginaryPartComplex(rootp[i])));
1578
1/10
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
64 PetscCallCUBLAS(cublasXCgemm(cublasv2handle,CUBLAS_OP_N,CUBLAS_OP_N,n,n,n,&cone,d_RR,n,d_expmA,n,&czero,d_Maux,n));
1579 64 SlepcSwap(d_expmA,d_Maux,aux);
1580
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
64 PetscCallCUDA(cudaMemcpy(d_RR,d_As,sizeof(PetscComplex)*n2,cudaMemcpyDeviceToDevice));
1581
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
64 PetscCall(shift_Cdiagonal(n,d_RR,n,-PetscRealPartComplex(rootq[i]),-PetscImaginaryPartComplex(rootq[i])));
1582
2/6
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
64 PetscCallMAGMA(magma_Cgesv_gpu,n,n,d_RR,n,piv,d_expmA,n);
1583 /* shift(n) + gemm + shift(n) + gesv */
1584
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
64 PetscCall(SlepcLogGpuFlopsComplex(1.0*n+(2.0*n*n*n)+1.0*n+(2.0*n*n*n/3.0+2.0*n*n*n)));
1585 }
1586 /* extra enumerator */
1587
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
32 for (i=minlen;i<irsize;i++) {
1588
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
16 PetscCallCUDA(cudaMemcpy(d_RR,d_As,sizeof(PetscComplex)*n2,cudaMemcpyDeviceToDevice));
1589
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
16 PetscCall(shift_Cdiagonal(n,d_RR,n,-PetscRealPartComplex(rootp[i]),-PetscImaginaryPartComplex(rootp[i])));
1590
1/10
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
16 PetscCallCUBLAS(cublasXCgemm(cublasv2handle,CUBLAS_OP_N,CUBLAS_OP_N,n,n,n,&cone,d_RR,n,d_expmA,n,&czero,d_Maux,n));
1591 16 SlepcSwap(d_expmA,d_Maux,aux);
1592
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
16 PetscCall(SlepcLogGpuFlopsComplex(1.0*n+2.0*n*n*n));
1593 }
1594 /* extra denominator */
1595
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
16 for (i=minlen;i<ipsize;i++) {
1596 PetscCallCUDA(cudaMemcpy(d_RR,d_As,sizeof(PetscComplex)*n2,cudaMemcpyDeviceToDevice));
1597 PetscCall(shift_Cdiagonal(n,d_RR,n,-PetscRealPartComplex(rootq[i]),-PetscImaginaryPartComplex(rootq[i])));
1598 PetscCallMAGMA(magma_Cgesv_gpu,n,n,d_RR,n,piv,d_expmA,n);
1599 PetscCall(SlepcLogGpuFlopsComplex(1.0*n+(2.0*n*n*n/3.0+2.0*n*n*n)));
1600 }
1601
1/10
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
16 PetscCallCUBLAS(cublasXCscal(cublasv2handle,n2,&mult,d_expmA,one));
1602
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
16 PetscCall(SlepcLogGpuFlopsComplex(1.0*n2));
1603
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
16 PetscCall(PetscFree2(rootp,rootq));
1604 }
1605
1606 #if !defined(PETSC_USE_COMPLEX)
1607
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
16 PetscCall(copy_array2D_C2S(n,n,d_Ba2,n,d_expmA,n));
1608 #else
1609
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
16 PetscCallCUDA(cudaMemcpy(d_Ba2,d_expmA,sizeof(PetscScalar)*n2,cudaMemcpyDeviceToDevice));
1610 #endif
1611
1612 /* perform repeated squaring */
1613
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
160 for (i=0;i<s;i++) { /* final squaring */
1614
1/10
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
128 PetscCallCUBLAS(cublasXgemm(cublasv2handle,CUBLAS_OP_N,CUBLAS_OP_N,n,n,n,&sone,d_Ba2,n,d_Ba2,n,&szero,d_sMaux,n));
1615 128 SlepcSwap(d_Ba2,d_sMaux,saux);
1616
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
128 PetscCall(PetscLogGpuFlops(2.0*n*n*n));
1617 }
1618
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
32 if (d_Ba2!=d_Ba) {
1619 PetscCallCUDA(cudaMemcpy(d_Ba,d_Ba2,sizeof(PetscScalar)*n2,cudaMemcpyDeviceToDevice));
1620 d_sMaux = d_Ba2;
1621 }
1622
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
32 if (shift) {
1623 32 expshift = PetscExpReal(shift);
1624
1/10
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
32 PetscCallCUBLAS(cublasXscal(cublasv2handle,n2,&expshift,d_Ba,one));
1625
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
32 PetscCall(PetscLogGpuFlops(1.0*n2));
1626 }
1627
1628
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
32 PetscCall(PetscLogGpuTimeEnd());
1629
1630 /* restore pointers */
1631 32 d_Maux = d_Maux2; d_expmA = d_expmA2; d_RR = d_RR2;
1632
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
32 PetscCall(MatDenseCUDARestoreArrayWrite(B,&d_Ba));
1633
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
32 PetscCallCUDA(cudaFree(d_isreal));
1634
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
32 PetscCallCUDA(cudaFree(d_sMaux));
1635
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
32 PetscCallCUDA(cudaFree(d_Maux));
1636
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
32 PetscCallCUDA(cudaFree(d_expmA));
1637
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
32 PetscCallCUDA(cudaFree(d_As));
1638
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
32 PetscCallCUDA(cudaFree(d_RR));
1639
2/4
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
32 PetscCall(PetscFree(piv));
1640 PetscFunctionReturn(PETSC_SUCCESS);
1641 }
1642 #endif /* PETSC_HAVE_MAGMA */
1643 #endif /* PETSC_HAVE_CUDA */
1644
1645 319 static PetscErrorCode FNView_Exp(FN fn,PetscViewer viewer)
1646 {
1647 319 PetscBool isascii;
1648 319 char str[50];
1649 319 const char *methodname[] = {
1650 "scaling & squaring, [m/m] Pade approximant (Higham)",
1651 "scaling & squaring, [6/6] Pade approximant",
1652 "scaling & squaring, subdiagonal Pade approximant (product form)",
1653 "scaling & squaring, subdiagonal Pade approximant (partial fraction)"
1654 };
1655 319 const int nmeth=PETSC_STATIC_ARRAY_LENGTH(methodname);
1656
1657
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
319 PetscFunctionBegin;
1658
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.
319 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
1659
1/2
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
319 if (isascii) {
1660
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
319 if (fn->beta==(PetscScalar)1.0) {
1661
6/8
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 8 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
269 if (fn->alpha==(PetscScalar)1.0) PetscCall(PetscViewerASCIIPrintf(viewer," exponential: exp(x)\n"));
1662 else {
1663
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.
105 PetscCall(SlepcSNPrintfScalar(str,sizeof(str),fn->alpha,PETSC_TRUE));
1664
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.
105 PetscCall(PetscViewerASCIIPrintf(viewer," exponential: exp(%s*x)\n",str));
1665 }
1666 } else {
1667
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(SlepcSNPrintfScalar(str,sizeof(str),fn->beta,PETSC_TRUE));
1668
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 (fn->alpha==(PetscScalar)1.0) PetscCall(PetscViewerASCIIPrintf(viewer," exponential: %s*exp(x)\n",str));
1669 else {
1670
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(PetscViewerASCIIPrintf(viewer," exponential: %s",str));
1671
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(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE));
1672
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(SlepcSNPrintfScalar(str,sizeof(str),fn->alpha,PETSC_TRUE));
1673
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(PetscViewerASCIIPrintf(viewer,"*exp(%s*x)\n",str));
1674
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(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE));
1675 }
1676 }
1677
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.
319 if (fn->method<nmeth) PetscCall(PetscViewerASCIIPrintf(viewer," computing matrix functions with: %s\n",methodname[fn->method]));
1678 }
1679
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.
59 PetscFunctionReturn(PETSC_SUCCESS);
1680 }
1681
1682 /*MC
1683 FNEXP - FNEXP = "exp" - The exponential function $f(x)=e^x$.
1684
1685 Level: beginner
1686
1687 .seealso: [](sec:fn), `FN`, `FNType`, `FNSetType()`
1688 M*/
1689
1690 1160 SLEPC_EXTERN PetscErrorCode FNCreate_Exp(FN fn)
1691 {
1692
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
1160 PetscFunctionBegin;
1693 1160 fn->ops->evaluatefunction = FNEvaluateFunction_Exp;
1694 1160 fn->ops->evaluatederivative = FNEvaluateDerivative_Exp;
1695 1160 fn->ops->evaluatefunctionmat[0] = FNEvaluateFunctionMat_Exp_Higham;
1696 1160 fn->ops->evaluatefunctionmat[1] = FNEvaluateFunctionMat_Exp_Pade;
1697 1160 fn->ops->evaluatefunctionmat[2] = FNEvaluateFunctionMat_Exp_GuettelNakatsukasa; /* product form */
1698 1160 fn->ops->evaluatefunctionmat[3] = FNEvaluateFunctionMat_Exp_GuettelNakatsukasa; /* partial fraction */
1699 #if defined(PETSC_HAVE_CUDA)
1700 255 fn->ops->evaluatefunctionmatcuda[1] = FNEvaluateFunctionMat_Exp_Pade_CUDA;
1701 #if defined(PETSC_HAVE_MAGMA)
1702 255 fn->ops->evaluatefunctionmatcuda[0] = FNEvaluateFunctionMat_Exp_Higham_CUDAm;
1703 255 fn->ops->evaluatefunctionmatcuda[1] = FNEvaluateFunctionMat_Exp_Pade_CUDAm;
1704 255 fn->ops->evaluatefunctionmatcuda[2] = FNEvaluateFunctionMat_Exp_GuettelNakatsukasa_CUDAm; /* product form */
1705 255 fn->ops->evaluatefunctionmatcuda[3] = FNEvaluateFunctionMat_Exp_GuettelNakatsukasa_CUDAm; /* partial fraction */
1706 #endif
1707 #endif
1708 1160 fn->ops->view = FNView_Exp;
1709
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.
1160 PetscFunctionReturn(PETSC_SUCCESS);
1710 }
1711