GCC Code Coverage Report

Directory: ./
File: src/sys/classes/fn/impls/fnutil.c
Date: 2026-07-06 04:10:30
Exec Total Coverage
Lines: 388 391 99.2%
Functions: 8 8 100.0%
Branches: 762 1602 47.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 Utility subroutines common to several impls
12 */
13
14 #include <slepc/private/fnimpl.h> /*I "slepcfn.h" I*/
15 #include <slepcblaslapack.h>
16
17 /*
18 Compute the square root of an upper quasi-triangular matrix T,
19 using Higham's algorithm (LAA 88, 1987). T is overwritten with sqrtm(T).
20 */
21 1360 static PetscErrorCode SlepcMatDenseSqrt(PetscBLASInt n,PetscScalar *T,PetscBLASInt ld)
22 {
23 1360 PetscScalar one=1.0,mone=-1.0;
24 1360 PetscReal scal;
25 1360 PetscBLASInt i,j,si,sj,r,ione=1;
26 #if !defined(PETSC_USE_COMPLEX)
27 609 PetscReal alpha,theta,mu,mu2;
28 #endif
29
30
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 1360 PetscFunctionBegin;
31
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 47296 for (j=0;j<n;j++) {
32 #if defined(PETSC_USE_COMPLEX)
33 27644 sj = 1;
34 27644 T[j+j*ld] = PetscSqrtScalar(T[j+j*ld]);
35 #else
36
4/4
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
✓ Branch 2 taken 5 times.
✓ Branch 3 taken 5 times.
 18292 sj = (j==n-1 || T[j+1+j*ld] == 0.0)? 1: 2;
37
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
 18292 if (sj==1) {
38
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 14292 PetscCheck(T[j+j*ld]>=0.0,PETSC_COMM_SELF,PETSC_ERR_USER_INPUT,"Matrix has a real negative eigenvalue, no real primary square root exists");
39 14292 T[j+j*ld] = PetscSqrtReal(T[j+j*ld]);
40 } else {
41 /* square root of 2x2 block */
42 4000 theta = (T[j+j*ld]+T[j+1+(j+1)*ld])/2.0;
43 4000 mu = (T[j+j*ld]-T[j+1+(j+1)*ld])/2.0;
44 4000 mu2 = -mu*mu-T[j+1+j*ld]*T[j+(j+1)*ld];
45 4000 mu = PetscSqrtReal(mu2);
46
1/2
✓ Branch 0 taken 5 times.
✗ Branch 1 not taken.
 4000 if (theta>0.0) alpha = PetscSqrtReal((theta+PetscSqrtReal(theta*theta+mu2))/2.0);
47 ✗ else alpha = mu/PetscSqrtReal(2.0*(-theta+PetscSqrtReal(theta*theta+mu2)));
48 4000 T[j+j*ld] /= 2.0*alpha;
49 4000 T[j+1+(j+1)*ld] /= 2.0*alpha;
50 4000 T[j+(j+1)*ld] /= 2.0*alpha;
51 4000 T[j+1+j*ld] /= 2.0*alpha;
52 4000 T[j+j*ld] += alpha-theta/(2.0*alpha);
53 4000 T[j+1+(j+1)*ld] += alpha-theta/(2.0*alpha);
54 }
55 #endif
56
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 1178366 for (i=j-1;i>=0;i--) {
57 #if defined(PETSC_USE_COMPLEX)
58 721688 si = 1;
59 #else
60
4/4
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
✓ Branch 2 taken 5 times.
✓ Branch 3 taken 5 times.
 410742 si = (i==0 || T[i+(i-1)*ld] == 0.0)? 1: 2;
61
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
 410742 if (si==2) i--;
62 #endif
63 /* solve Sylvester equation of order si x sj */
64 1132430 r = j-i-si;
65
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.
 1132430 if (r) PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&si,&sj,&r,&mone,T+i+(i+si)*ld,&ld,T+i+si+j*ld,&ld,&one,T+i+j*ld,&ld));
66
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.
 1132430 PetscCallLAPACKInfo("LAPACKtrsyl",LAPACKtrsyl_("N","N",&ione,&si,&sj,T+i+i*ld,&ld,T+j+j*ld,&ld,T+i+j*ld,&ld,&scal,&info));
67
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 1132430 PetscCheck(scal==1.0,PETSC_COMM_SELF,PETSC_ERR_SUP,"Current implementation cannot handle scale factor %g",(double)scal);
68 }
69
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 10 times.
 45936 if (sj==2) j++;
70 }
71
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.
 266 PetscFunctionReturn(PETSC_SUCCESS);
72 }
73
74 #define BLOCKSIZE 64
75
76 /*
77 Schur method for the square root of an upper quasi-triangular matrix T.
78 T is overwritten with sqrtm(T).
79 If firstonly then only the first column of T will contain relevant values.
80 */
81 814 PetscErrorCode FNSqrtmSchur(FN fn,PetscBLASInt n,PetscScalar *T,PetscBLASInt ld,PetscBool firstonly)
82 {
83 814 PetscBLASInt i,j,k,r,ione=1,sdim,lwork,*s,*p,bs=BLOCKSIZE;
84 814 PetscScalar *wr,*W,*Q,*work,one=1.0,zero=0.0,mone=-1.0;
85 814 PetscInt m,nblk;
86 814 PetscReal scal;
87 #if defined(PETSC_USE_COMPLEX)
88 444 PetscReal *rwork;
89 #else
90 370 PetscReal *wi;
91 #endif
92
93
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 814 PetscFunctionBegin;
94 814 m = n;
95 814 nblk = (m+bs-1)/bs;
96 814 lwork = 5*n;
97
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 814 k = firstonly? 1: n;
98
99 /* compute Schur decomposition A*Q = Q*T */
100 #if !defined(PETSC_USE_COMPLEX)
101
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.
 370 PetscCall(PetscMalloc7(m,&wr,m,&wi,m*k,&W,m*m,&Q,lwork,&work,nblk,&s,nblk,&p));
102
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.
 370 PetscCallLAPACKInfo("LAPACKgees",LAPACKgees_("V","N",NULL,&n,T,&ld,&sdim,wr,wi,Q,&ld,work,&lwork,NULL,&info));
103 #else
104
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.
 444 PetscCall(PetscMalloc7(m,&wr,m,&rwork,m*k,&W,m*m,&Q,lwork,&work,nblk,&s,nblk,&p));
105
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.
 444 PetscCallLAPACKInfo("LAPACKgees",LAPACKgees_("V","N",NULL,&n,T,&ld,&sdim,wr,Q,&ld,work,&lwork,rwork,NULL,&info));
106 #endif
107
108 /* determine block sizes and positions, to avoid cutting 2x2 blocks */
109 814 j = 0;
110 814 p[j] = 0;
111 1844 do {
112 1360 s[j] = PetscMin(bs,n-p[j]);
113 #if !defined(PETSC_USE_COMPLEX)
114
4/4
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
✓ Branch 2 taken 5 times.
✓ Branch 3 taken 5 times.
 609 if (p[j]+s[j]!=n && T[p[j]+s[j]+(p[j]+s[j]-1)*ld]!=0.0) s[j]++;
115 #endif
116
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 1360 if (p[j]+s[j]==n) break;
117 546 j++;
118 546 p[j] = p[j-1]+s[j-1];
119 62 } while (1);
120 2174 nblk = j+1;
121
122
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 2174 for (j=0;j<nblk;j++) {
123 /* evaluate f(T_jj) */
124
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 1360 PetscCall(SlepcMatDenseSqrt(s[j],T+p[j]+p[j]*ld,ld));
125
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 1906 for (i=j-1;i>=0;i--) {
126 /* solve Sylvester equation for block (i,j) */
127 546 r = p[j]-p[i]-s[i];
128
1/22
✗ 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.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
 546 if (r) PetscCallBLAS("BLASgemm",BLASgemm_("N","N",s+i,s+j,&r,&mone,T+p[i]+(p[i]+s[i])*ld,&ld,T+p[i]+s[i]+p[j]*ld,&ld,&one,T+p[i]+p[j]*ld,&ld));
129
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.
 546 PetscCallLAPACKInfo("LAPACKtrsyl",LAPACKtrsyl_("N","N",&ione,s+i,s+j,T+p[i]+p[i]*ld,&ld,T+p[j]+p[j]*ld,&ld,T+p[i]+p[j]*ld,&ld,&scal,&info));
130
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 546 PetscCheck(scal==1.0,PETSC_COMM_SELF,PETSC_ERR_SUP,"Current implementation cannot handle scale factor %g",(double)scal);
131 }
132 }
133
134 /* backtransform B = Q*T*Q' */
135
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.
 814 PetscCallBLAS("BLASgemm",BLASgemm_("N","C",&n,&k,&n,&one,T,&ld,Q,&ld,&zero,W,&ld));
136
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.
 814 PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n,&k,&n,&one,Q,&ld,W,&ld,&zero,T,&ld));
137
138 /* flop count: Schur decomposition, triangular square root, and backtransform */
139
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.
 814 PetscCall(PetscLogFlops(25.0*n*n*n+n*n*n/3.0+4.0*n*n*k));
140
141 #if !defined(PETSC_USE_COMPLEX)
142
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.
 370 PetscCall(PetscFree7(wr,wi,W,Q,work,s,p));
143 #else
144
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.
 444 PetscCall(PetscFree7(wr,rwork,W,Q,work,s,p));
145 #endif
146
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.
 158 PetscFunctionReturn(PETSC_SUCCESS);
147 }
148
149 #define DBMAXIT 25
150
151 /*
152 Computes the principal square root of the matrix T using the product form
153 of the Denman-Beavers iteration.
154 T is overwritten with sqrtm(T) or inv(sqrtm(T)) depending on flag inv.
155 */
156 252 PetscErrorCode FNSqrtmDenmanBeavers(FN fn,PetscBLASInt n,PetscScalar *T,PetscBLASInt ld,PetscBool inv)
157 {
158 252 PetscScalar *Told,*M=NULL,*invM,*work,work1,prod,alpha;
159 252 PetscScalar szero=0.0,sone=1.0,smone=-1.0,spfive=0.5,sp25=0.25;
160 252 PetscReal tol,Mres=0.0,detM,g,reldiff,fnormdiff,fnormT,rwork[1];
161 252 PetscBLASInt N,i,it,*piv=NULL,query=-1,lwork;
162 252 const PetscBLASInt one=1;
163 252 PetscBool converged=PETSC_FALSE,scale;
164 252 unsigned int ftz;
165
166
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 252 PetscFunctionBegin;
167 252 N = n*n;
168 252 tol = PetscSqrtReal((PetscReal)n)*PETSC_MACHINE_EPSILON/2;
169 252 scale = PetscDefined(USE_REAL_SINGLE)? PETSC_FALSE: PETSC_TRUE;
170
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.
 252 PetscCall(SlepcSetFlushToZero(&ftz));
171
172 /* query work size */
173
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.
 252 PetscCallLAPACKInfo("LAPACKgetri",LAPACKgetri_(&n,M,&ld,piv,&work1,&query,&info));
174
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 252 PetscCall(PetscBLASIntCast((PetscInt)PetscRealPart(work1),&lwork));
175
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 252 PetscCall(PetscMalloc5(lwork,&work,n,&piv,n*n,&Told,n*n,&M,n*n,&invM));
176
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 252 PetscCall(PetscArraycpy(M,T,n*n));
177
178
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 252 if (inv) { /* start recurrence with I instead of A */
179
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.
 120 PetscCall(PetscArrayzero(T,n*n));
180
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 1320 for (i=0;i<n;i++) T[i+i*ld] += 1.0;
181 }
182
183
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 1512 for (it=0;it<DBMAXIT && !converged;it++) {
184
185
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 1260 if (scale) { /* g = (abs(det(M)))^(-1/(2*n)) */
186
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.
 756 PetscCall(PetscArraycpy(invM,M,n*n));
187
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.
 756 PetscCallLAPACKInfo("LAPACKgetrf",LAPACKgetrf_(&n,&n,invM,&ld,piv,&info));
188 756 prod = invM[0];
189
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 40284 for (i=1;i<n;i++) prod *= invM[i+i*ld];
190 756 detM = PetscAbsScalar(prod);
191
1/2
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
 756 g = (detM>PETSC_MAX_REAL)? 0.5: PetscPowReal(detM,-1.0/(2.0*n));
192 756 alpha = g;
193
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.
 756 PetscCallBLAS("BLASscal",BLASscal_(&N,&alpha,T,&one));
194 756 alpha = g*g;
195
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.
 756 PetscCallBLAS("BLASscal",BLASscal_(&N,&alpha,M,&one));
196
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.
 756 PetscCall(PetscLogFlops(2.0*n*n*n/3.0+2.0*n*n));
197 }
198
199
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.
 1260 PetscCall(PetscArraycpy(Told,T,n*n));
200
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.
 1260 PetscCall(PetscArraycpy(invM,M,n*n));
201
202
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.
 1260 PetscCallLAPACKInfo("LAPACKgetrf",LAPACKgetrf_(&n,&n,invM,&ld,piv,&info));
203
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.
 1260 PetscCallLAPACKInfo("LAPACKgetri",LAPACKgetri_(&n,invM,&ld,piv,work,&lwork,&info));
204
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.
 1260 PetscCall(PetscLogFlops(2.0*n*n*n/3.0+4.0*n*n*n/3.0));
205
206
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 72000 for (i=0;i<n;i++) invM[i+i*ld] += 1.0;
207
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.
 1260 PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n,&n,&n,&spfive,Told,&ld,invM,&ld,&szero,T,&ld));
208
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 72000 for (i=0;i<n;i++) invM[i+i*ld] -= 1.0;
209
210
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.
 1260 PetscCallBLAS("BLASaxpy",BLASaxpy_(&N,&sone,invM,&one,M,&one));
211
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.
 1260 PetscCallBLAS("BLASscal",BLASscal_(&N,&sp25,M,&one));
212
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 72000 for (i=0;i<n;i++) M[i+i*ld] -= 0.5;
213
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.
 1260 PetscCall(PetscLogFlops(2.0*n*n*n+2.0*n*n));
214
215 1260 Mres = LAPACKlange_("F",&n,&n,M,&n,rwork);
216
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 73260 for (i=0;i<n;i++) M[i+i*ld] += 1.0;
217
218
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 1260 if (scale) {
219 /* reldiff = norm(T - Told,'fro')/norm(T,'fro') */
220
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.
 756 PetscCallBLAS("BLASaxpy",BLASaxpy_(&N,&smone,T,&one,Told,&one));
221 756 fnormdiff = LAPACKlange_("F",&n,&n,Told,&n,rwork);
222 756 fnormT = LAPACKlange_("F",&n,&n,T,&n,rwork);
223
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.
 756 PetscCall(PetscLogFlops(7.0*n*n));
224 756 reldiff = fnormdiff/fnormT;
225
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.
 756 PetscCall(PetscInfo(fn,"it: %" PetscBLASInt_FMT " reldiff: %g scale: %g tol*scale: %g\n",it,(double)reldiff,(double)g,(double)(tol*g)));
226
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 756 if (reldiff<1e-2) scale = PETSC_FALSE; /* Switch off scaling */
227 }
228
229
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 1260 if (Mres<=tol) converged = PETSC_TRUE;
230 }
231
232
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 252 PetscCheck(Mres<=tol,PETSC_COMM_SELF,PETSC_ERR_LIB,"SQRTM not converged after %d iterations",DBMAXIT);
233
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.
 252 PetscCall(PetscFree5(work,piv,Told,M,invM));
234
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.
 74 PetscCall(SlepcResetFlushToZero(&ftz));
235
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.
 48 PetscFunctionReturn(PETSC_SUCCESS);
236 }
237
238 #define NSMAXIT 50
239
240 /*
241 Computes the principal square root of the matrix A using the Newton-Schulz iteration.
242 T is overwritten with sqrtm(T) or inv(sqrtm(T)) depending on flag inv.
243 */
244 252 PetscErrorCode FNSqrtmNewtonSchulz(FN fn,PetscBLASInt n,PetscScalar *A,PetscBLASInt ld,PetscBool inv)
245 {
246 252 PetscScalar *Y=A,*Yold,*Z,*Zold,*M;
247 252 PetscScalar szero=0.0,sone=1.0,smone=-1.0,spfive=0.5,sthree=3.0;
248 252 PetscReal sqrtnrm,tol,Yres=0.0,nrm,rwork[1],done=1.0;
249 252 PetscBLASInt i,it,N,one=1,zero=0;
250 252 PetscBool converged=PETSC_FALSE;
251 252 unsigned int ftz;
252
253
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 252 PetscFunctionBegin;
254 252 N = n*n;
255 252 tol = PetscSqrtReal((PetscReal)n)*PETSC_MACHINE_EPSILON/2;
256
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.
 252 PetscCall(SlepcSetFlushToZero(&ftz));
257
258
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.
 252 PetscCall(PetscMalloc4(N,&Yold,N,&Z,N,&Zold,N,&M));
259
260 /* scale */
261
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.
 252 PetscCall(PetscArraycpy(Z,A,N));
262
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 13680 for (i=0;i<n;i++) Z[i+i*ld] -= 1.0;
263 252 nrm = LAPACKlange_("fro",&n,&n,Z,&n,rwork);
264 252 sqrtnrm = PetscSqrtReal(nrm);
265
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.
 252 PetscCallLAPACKInfo("LAPACKlascl",LAPACKlascl_("G",&zero,&zero,&nrm,&done,&N,&one,A,&N,&info));
266 252 tol *= nrm;
267
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.
 252 PetscCall(PetscInfo(fn,"||I-A||_F = %g, new tol: %g\n",(double)nrm,(double)tol));
268
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.
 252 PetscCall(PetscLogFlops(2.0*n*n));
269
270 /* Z = I */
271
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.
 252 PetscCall(PetscArrayzero(Z,N));
272
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 13680 for (i=0;i<n;i++) Z[i+i*ld] = 1.0;
273
274
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 2956 for (it=0;it<NSMAXIT && !converged;it++) {
275 /* Yold = Y, Zold = Z */
276
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.
 2704 PetscCall(PetscArraycpy(Yold,Y,N));
277
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.
 2704 PetscCall(PetscArraycpy(Zold,Z,N));
278
279 /* M = (3*I-Zold*Yold) */
280
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.
 2704 PetscCall(PetscArrayzero(M,N));
281
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 171440 for (i=0;i<n;i++) M[i+i*ld] = sthree;
282
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.
 2704 PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n,&n,&n,&smone,Zold,&ld,Yold,&ld,&sone,M,&ld));
283
284 /* Y = (1/2)*Yold*M, Z = (1/2)*M*Zold */
285
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.
 2704 PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n,&n,&n,&spfive,Yold,&ld,M,&ld,&szero,Y,&ld));
286
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.
 2704 PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n,&n,&n,&spfive,M,&ld,Zold,&ld,&szero,Z,&ld));
287
288 /* reldiff = norm(Y-Yold,'fro')/norm(Y,'fro') */
289
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.
 2704 PetscCallBLAS("BLASaxpy",BLASaxpy_(&N,&smone,Y,&one,Yold,&one));
290 2704 Yres = LAPACKlange_("fro",&n,&n,Yold,&n,rwork);
291
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 2704 PetscCheck(!PetscIsNanReal(Yres),PETSC_COMM_SELF,PETSC_ERR_FP,"The computed norm is not-a-number");
292
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 2704 if (Yres<=tol) converged = PETSC_TRUE;
293
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.
 2704 PetscCall(PetscInfo(fn,"it: %" PetscBLASInt_FMT " res: %g\n",it,(double)Yres));
294
295
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.
 2704 PetscCall(PetscLogFlops(6.0*n*n*n+2.0*n*n));
296 }
297
298
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 252 PetscCheck(Yres<=tol,PETSC_COMM_SELF,PETSC_ERR_LIB,"SQRTM not converged after %d iterations",NSMAXIT);
299
300 /* undo scaling */
301
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 252 if (inv) {
302
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(A,Z,N));
303
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.
 120 PetscCallLAPACKInfo("LAPACKlascl",LAPACKlascl_("G",&zero,&zero,&sqrtnrm,&done,&N,&one,A,&N,&info));
304
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.
 132 } else PetscCallLAPACKInfo("LAPACKlascl",LAPACKlascl_("G",&zero,&zero,&done,&sqrtnrm,&N,&one,A,&N,&info));
305
306
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.
 252 PetscCall(PetscFree4(Yold,Z,Zold,M));
307
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.
 74 PetscCall(SlepcResetFlushToZero(&ftz));
308
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.
 48 PetscFunctionReturn(PETSC_SUCCESS);
309 }
310
311 #if defined(PETSC_HAVE_CUDA)
312 #include "../src/sys/classes/fn/impls/cuda/fnutilcuda.h"
313 #include <slepccupmblas.h>
314
315 /*
316 * Matrix square root by Newton-Schulz iteration. CUDA version.
317 * Computes the principal square root of the matrix A using the
318 * Newton-Schulz iteration. A is overwritten with sqrtm(A).
319 */
320 48 PetscErrorCode FNSqrtmNewtonSchulz_CUDA(FN fn,PetscBLASInt n,PetscScalar *d_A,PetscBLASInt ld,PetscBool inv)
321 {
322 48 PetscScalar *d_Yold,*d_Z,*d_Zold,*d_M,alpha;
323 48 PetscReal nrm,sqrtnrm,tol,Yres=0.0;
324 48 const PetscScalar szero=0.0,sone=1.0,smone=-1.0,spfive=0.5,sthree=3.0;
325 48 PetscInt it;
326 48 PetscBLASInt N;
327 48 const PetscBLASInt one=1;
328 48 PetscBool converged=PETSC_FALSE;
329 48 cublasHandle_t cublasv2handle;
330
331 48 PetscFunctionBegin;
332
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 48 PetscCall(PetscDeviceInitialize(PETSC_DEVICE_CUDA)); /* For CUDA event timers */
333
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 48 PetscCall(PetscCUBLASGetHandle(&cublasv2handle));
334 48 N = n*n;
335 48 tol = PetscSqrtReal((PetscReal)n)*PETSC_MACHINE_EPSILON/2;
336
337
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 48 PetscCallCUDA(cudaMalloc((void **)&d_Yold,sizeof(PetscScalar)*N));
338
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 48 PetscCallCUDA(cudaMalloc((void **)&d_Z,sizeof(PetscScalar)*N));
339
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 48 PetscCallCUDA(cudaMalloc((void **)&d_Zold,sizeof(PetscScalar)*N));
340
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 48 PetscCallCUDA(cudaMalloc((void **)&d_M,sizeof(PetscScalar)*N));
341
342
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 48 PetscCall(PetscLogGpuTimeBegin());
343
344 /* Z = I; */
345
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 48 PetscCallCUDA(cudaMemset(d_Z,0,sizeof(PetscScalar)*N));
346
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 48 PetscCall(set_diagonal(n,d_Z,ld,sone));
347
348 /* scale */
349
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,N,&smone,d_A,one,d_Z,one));
350
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(cublasXnrm2(cublasv2handle,N,d_Z,one,&nrm));
351 48 sqrtnrm = PetscSqrtReal(nrm);
352 48 alpha = 1.0/nrm;
353
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,N,&alpha,d_A,one));
354 48 tol *= nrm;
355
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 48 PetscCall(PetscInfo(fn,"||I-A||_F = %g, new tol: %g\n",(double)nrm,(double)tol));
356
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 48 PetscCall(PetscLogGpuFlops(2.0*n*n));
357
358 /* Z = I; */
359
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 48 PetscCallCUDA(cudaMemset(d_Z,0,sizeof(PetscScalar)*N));
360
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 48 PetscCall(set_diagonal(n,d_Z,ld,sone));
361
362
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
 560 for (it=0;it<NSMAXIT && !converged;it++) {
363 /* Yold = Y, Zold = Z */
364
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 512 PetscCallCUDA(cudaMemcpy(d_Yold,d_A,sizeof(PetscScalar)*N,cudaMemcpyDeviceToDevice));
365
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 512 PetscCallCUDA(cudaMemcpy(d_Zold,d_Z,sizeof(PetscScalar)*N,cudaMemcpyDeviceToDevice));
366
367 /* M = (3*I - Zold*Yold) */
368
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 512 PetscCallCUDA(cudaMemset(d_M,0,sizeof(PetscScalar)*N));
369
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 512 PetscCall(set_diagonal(n,d_M,ld,sthree));
370
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.
 512 PetscCallCUBLAS(cublasXgemm(cublasv2handle,CUBLAS_OP_N,CUBLAS_OP_N,n,n,n,&smone,d_Zold,ld,d_Yold,ld,&sone,d_M,ld));
371
372 /* Y = (1/2) * Yold * M, Z = (1/2) * M * Zold */
373
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.
 512 PetscCallCUBLAS(cublasXgemm(cublasv2handle,CUBLAS_OP_N,CUBLAS_OP_N,n,n,n,&spfive,d_Yold,ld,d_M,ld,&szero,d_A,ld));
374
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.
 512 PetscCallCUBLAS(cublasXgemm(cublasv2handle,CUBLAS_OP_N,CUBLAS_OP_N,n,n,n,&spfive,d_M,ld,d_Zold,ld,&szero,d_Z,ld));
375
376 /* reldiff = norm(Y-Yold,'fro')/norm(Y,'fro') */
377
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.
 512 PetscCallCUBLAS(cublasXaxpy(cublasv2handle,N,&smone,d_A,one,d_Yold,one));
378
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.
 512 PetscCallCUBLAS(cublasXnrm2(cublasv2handle,N,d_Yold,one,&Yres));
379
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 512 PetscCheck(!PetscIsNanReal(Yres),PETSC_COMM_SELF,PETSC_ERR_FP,"The computed norm is not-a-number");
380
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
 512 if (Yres<=tol) converged = PETSC_TRUE;
381
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 512 PetscCall(PetscInfo(fn,"it: %" PetscInt_FMT " res: %g\n",it,(double)Yres));
382
383
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 512 PetscCall(PetscLogGpuFlops(6.0*n*n*n+2.0*n*n));
384 }
385
386
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 48 PetscCheck(Yres<=tol,PETSC_COMM_SELF,PETSC_ERR_LIB,"SQRTM not converged after %d iterations", NSMAXIT);
387
388 /* undo scaling */
389
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
 48 if (inv) {
390 24 alpha = 1.0/sqrtnrm;
391
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.
 24 PetscCallCUBLAS(cublasXscal(cublasv2handle,N,&alpha,d_Z,one));
392
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 24 PetscCallCUDA(cudaMemcpy(d_A,d_Z,sizeof(PetscScalar)*N,cudaMemcpyDeviceToDevice));
393 } else {
394 24 alpha = sqrtnrm;
395
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.
 24 PetscCallCUBLAS(cublasXscal(cublasv2handle,N,&alpha,d_A,one));
396 }
397
398
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 48 PetscCall(PetscLogGpuTimeEnd());
399
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 48 PetscCallCUDA(cudaFree(d_Yold));
400
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 48 PetscCallCUDA(cudaFree(d_Z));
401
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 48 PetscCallCUDA(cudaFree(d_Zold));
402
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 48 PetscCallCUDA(cudaFree(d_M));
403 PetscFunctionReturn(PETSC_SUCCESS);
404 }
405
406 #if defined(PETSC_HAVE_MAGMA)
407 #include <slepcmagma.h>
408
409 /*
410 * Matrix square root by product form of Denman-Beavers iteration. CUDA version.
411 * Computes the principal square root of the matrix T using the product form
412 * of the Denman-Beavers iteration. T is overwritten with sqrtm(T).
413 */
414 48 PetscErrorCode FNSqrtmDenmanBeavers_CUDAm(FN fn,PetscBLASInt n,PetscScalar *d_T,PetscBLASInt ld,PetscBool inv)
415 {
416 48 PetscScalar *d_Told,*d_M,*d_invM,*d_work,prod,szero=0.0,sone=1.0,smone=-1.0,spfive=0.5,sneg_pfive=-0.5,sp25=0.25,alpha;
417 48 PetscReal tol,Mres=0.0,detM,g,reldiff,fnormdiff,fnormT;
418 48 PetscInt it,lwork,nb;
419 48 PetscBLASInt N,one=1,*piv=NULL;
420 48 PetscBool converged=PETSC_FALSE,scale;
421 48 cublasHandle_t cublasv2handle;
422
423 48 PetscFunctionBegin;
424
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 48 PetscCall(PetscDeviceInitialize(PETSC_DEVICE_CUDA)); /* For CUDA event timers */
425
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 48 PetscCall(PetscCUBLASGetHandle(&cublasv2handle));
426
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 48 PetscCall(SlepcMagmaInit());
427 48 N = n*n;
428 48 scale = PetscDefined(USE_REAL_SINGLE)? PETSC_FALSE: PETSC_TRUE;
429 48 tol = PetscSqrtReal((PetscReal)n)*PETSC_MACHINE_EPSILON/2;
430
431 /* query work size */
432 48 nb = magma_get_xgetri_nb(n);
433 48 lwork = nb*n;
434
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 48 PetscCall(PetscMalloc1(n,&piv));
435
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 48 PetscCallCUDA(cudaMalloc((void **)&d_work,sizeof(PetscScalar)*lwork));
436
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 48 PetscCallCUDA(cudaMalloc((void **)&d_Told,sizeof(PetscScalar)*N));
437
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 48 PetscCallCUDA(cudaMalloc((void **)&d_M,sizeof(PetscScalar)*N));
438
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 48 PetscCallCUDA(cudaMalloc((void **)&d_invM,sizeof(PetscScalar)*N));
439
440
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 48 PetscCall(PetscLogGpuTimeBegin());
441
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 48 PetscCallCUDA(cudaMemcpy(d_M,d_T,sizeof(PetscScalar)*N,cudaMemcpyDeviceToDevice));
442
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
 48 if (inv) { /* start recurrence with I instead of A */
443
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 24 PetscCallCUDA(cudaMemset(d_T,0,sizeof(PetscScalar)*N));
444
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 24 PetscCall(set_diagonal(n,d_T,ld,1.0));
445 }
446
447
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
 304 for (it=0;it<DBMAXIT && !converged;it++) {
448
449
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
 256 if (scale) { /* g = (abs(det(M)))^(-1/(2*n)); */
450
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 160 PetscCallCUDA(cudaMemcpy(d_invM,d_M,sizeof(PetscScalar)*N,cudaMemcpyDeviceToDevice));
451
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.
 160 PetscCallMAGMA(magma_xgetrf_gpu,n,n,d_invM,ld,piv);
452
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 160 PetscCall(mult_diagonal(n,d_invM,ld,&prod));
453 160 detM = PetscAbsScalar(prod);
454
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 160 g = (detM>PETSC_MAX_REAL)? 0.5: PetscPowReal(detM,-1.0/(2.0*n));
455 160 alpha = g;
456
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.
 160 PetscCallCUBLAS(cublasXscal(cublasv2handle,N,&alpha,d_T,one));
457 160 alpha = g*g;
458
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.
 160 PetscCallCUBLAS(cublasXscal(cublasv2handle,N,&alpha,d_M,one));
459
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 160 PetscCall(PetscLogGpuFlops(2.0*n*n*n/3.0+2.0*n*n));
460 }
461
462
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 256 PetscCallCUDA(cudaMemcpy(d_Told,d_T,sizeof(PetscScalar)*N,cudaMemcpyDeviceToDevice));
463
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 256 PetscCallCUDA(cudaMemcpy(d_invM,d_M,sizeof(PetscScalar)*N,cudaMemcpyDeviceToDevice));
464
465
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.
 256 PetscCallMAGMA(magma_xgetrf_gpu,n,n,d_invM,ld,piv);
466
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.
 256 PetscCallMAGMA(magma_xgetri_gpu,n,d_invM,ld,piv,d_work,lwork);
467
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 256 PetscCall(PetscLogGpuFlops(2.0*n*n*n/3.0+4.0*n*n*n/3.0));
468
469
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 256 PetscCall(shift_diagonal(n,d_invM,ld,sone));
470
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.
 256 PetscCallCUBLAS(cublasXgemm(cublasv2handle,CUBLAS_OP_N,CUBLAS_OP_N,n,n,n,&spfive,d_Told,ld,d_invM,ld,&szero,d_T,ld));
471
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 256 PetscCall(shift_diagonal(n,d_invM,ld,smone));
472
473
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.
 256 PetscCallCUBLAS(cublasXaxpy(cublasv2handle,N,&sone,d_invM,one,d_M,one));
474
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.
 256 PetscCallCUBLAS(cublasXscal(cublasv2handle,N,&sp25,d_M,one));
475
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 256 PetscCall(shift_diagonal(n,d_M,ld,sneg_pfive));
476
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 256 PetscCall(PetscLogGpuFlops(2.0*n*n*n+2.0*n*n));
477
478
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.
 256 PetscCallCUBLAS(cublasXnrm2(cublasv2handle,N,d_M,one,&Mres));
479
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 256 PetscCall(shift_diagonal(n,d_M,ld,sone));
480
481
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
 256 if (scale) {
482 /* reldiff = norm(T - Told,'fro')/norm(T,'fro'); */
483
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.
 160 PetscCallCUBLAS(cublasXaxpy(cublasv2handle,N,&smone,d_T,one,d_Told,one));
484
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.
 160 PetscCallCUBLAS(cublasXnrm2(cublasv2handle,N,d_Told,one,&fnormdiff));
485
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.
 160 PetscCallCUBLAS(cublasXnrm2(cublasv2handle,N,d_T,one,&fnormT));
486
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 160 PetscCall(PetscLogGpuFlops(7.0*n*n));
487 160 reldiff = fnormdiff/fnormT;
488
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 160 PetscCall(PetscInfo(fn,"it: %" PetscInt_FMT " reldiff: %g scale: %g tol*scale: %g\n",it,(double)reldiff,(double)g,(double)tol*g));
489
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
 160 if (reldiff<1e-2) scale = PETSC_FALSE; /* Switch to no scaling. */
490 }
491
492
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 256 PetscCall(PetscInfo(fn,"it: %" PetscInt_FMT " Mres: %g\n",it,(double)Mres));
493
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
 256 if (Mres<=tol) converged = PETSC_TRUE;
494 }
495
496
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 48 PetscCheck(Mres<=tol,PETSC_COMM_SELF,PETSC_ERR_LIB,"SQRTM not converged after %d iterations", DBMAXIT);
497
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 48 PetscCall(PetscLogGpuTimeEnd());
498
2/4
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
 48 PetscCall(PetscFree(piv));
499
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 48 PetscCallCUDA(cudaFree(d_work));
500
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 48 PetscCallCUDA(cudaFree(d_Told));
501
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 48 PetscCallCUDA(cudaFree(d_M));
502
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 48 PetscCallCUDA(cudaFree(d_invM));
503 PetscFunctionReturn(PETSC_SUCCESS);
504 }
505 #endif /* PETSC_HAVE_MAGMA */
506
507 #endif /* PETSC_HAVE_CUDA */
508
509 #define ITMAX 5
510
511 /*
512 Estimate norm(A^m,1) by block 1-norm power method (required workspace is 11*n)
513 */
514 280 static PetscErrorCode SlepcNormEst1(PetscBLASInt n,PetscScalar *A,PetscInt m,PetscScalar *work,PetscRandom rand,PetscReal *nrm)
515 {
516 280 PetscScalar *X,*Y,*Z,*S,*S_old,*aux,val,sone=1.0,szero=0.0;
517 280 PetscReal est=0.0,est_old,vals[2]={0.0,0.0},*zvals,maxzval[2],raux;
518 280 PetscBLASInt i,j,t=2,it=0,ind[2],est_j=0,m1;
519
520
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 280 PetscFunctionBegin;
521 280 X = work;
522 280 Y = work + 2*n;
523 280 Z = work + 4*n;
524 280 S = work + 6*n;
525 280 S_old = work + 8*n;
526 280 zvals = (PetscReal*)(work + 10*n);
527
528
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 33980 for (i=0;i<n;i++) { /* X has columns of unit 1-norm */
529 33700 X[i] = 1.0/n;
530
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.
 33700 PetscCall(PetscRandomGetValue(rand,&val));
531
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 33700 if (PetscRealPart(val) < 0.5) X[i+n] = -1.0/n;
532 16595 else X[i+n] = 1.0/n;
533 }
534
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 67680 for (i=0;i<t*n;i++) S[i] = 0.0;
535 280 ind[0] = 0; ind[1] = 0;
536 280 est_old = 0;
537 612 while (1) {
538 612 it++;
539
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 2676 for (j=0;j<m;j++) { /* Y = A^m*X */
540
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.
 2064 PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n,&t,&n,&sone,A,&n,X,&n,&szero,Y,&n));
541
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 2064 if (j<m-1) SlepcSwap(X,Y,aux);
542 }
543
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 1836 for (j=0;j<t;j++) { /* vals[j] = norm(Y(:,j),1) */
544 1224 vals[j] = 0.0;
545
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 148536 for (i=0;i<n;i++) vals[j] += PetscAbsScalar(Y[i+j*n]);
546 }
547
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 612 if (vals[0]<vals[1]) {
548 340 SlepcSwap(vals[0],vals[1],raux);
549 340 m1 = 1;
550 } else m1 = 0;
551 612 est = vals[0];
552
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 612 if (est>est_old || it==2) est_j = ind[m1];
553
3/4
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
✓ Branch 2 taken 10 times.
✗ Branch 3 not taken.
 612 if (it>=2 && est<=est_old) {
554 est = est_old;
555 break;
556 }
557 134748 est_old = est;
558
1/2
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
 360 if (it>ITMAX) break;
559 147924 SlepcSwap(S,S_old,aux);
560
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 147924 for (i=0;i<t*n;i++) { /* S = sign(Y) */
561
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 256856 S[i] = (PetscRealPart(Y[i]) < 0.0)? -1.0: 1.0;
562 }
563
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 2676 for (j=0;j<m;j++) { /* Z = (A^T)^m*S */
564
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.
 2064 PetscCallBLAS("BLASgemm",BLASgemm_("C","N",&n,&t,&n,&sone,A,&n,S,&n,&szero,Z,&n));
565
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 2064 if (j<m-1) SlepcSwap(S,Z,aux);
566 }
567 612 maxzval[0] = -1; maxzval[1] = -1;
568 612 ind[0] = 0; ind[1] = 0;
569
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 74268 for (i=0;i<n;i++) { /* zvals[i] = norm(Z(i,:),inf) */
570
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 73656 zvals[i] = PetscMax(PetscAbsScalar(Z[i+0*n]),PetscAbsScalar(Z[i+1*n]));
571
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 73656 if (zvals[i]>maxzval[0]) {
572 7251 maxzval[0] = zvals[i];
573 7251 ind[0] = i;
574
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 66405 } else if (zvals[i]>maxzval[1]) {
575 3362 maxzval[1] = zvals[i];
576 3362 ind[1] = i;
577 }
578 }
579
4/4
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
✓ Branch 2 taken 5 times.
✓ Branch 3 taken 10 times.
 612 if (it>=2 && maxzval[0]==zvals[est_j]) break;
580
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 80244 for (i=0;i<t*n;i++) X[i] = 0.0;
581
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 996 for (j=0;j<t;j++) X[ind[j]+j*n] = 1.0;
582 }
583 280 *nrm = est;
584 /* Flop count is roughly (it * 2*m * t*gemv) = 4*its*m*t*n*n */
585
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.
 280 PetscCall(PetscLogFlops(4.0*it*m*t*n*n));
586
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.
 280 PetscFunctionReturn(PETSC_SUCCESS);
587 }
588
589 #define SMALLN 100
590
591 /*
592 Estimate norm(A^m,1) (required workspace is 2*n*n)
593 */
594 5390 PetscErrorCode SlepcNormAm(PetscBLASInt n,PetscScalar *A,PetscInt m,PetscScalar *work,PetscRandom rand,PetscReal *nrm)
595 {
596 5390 PetscScalar *v=work,*w=work+n*n,*aux,sone=1.0,szero=0.0;
597 5390 PetscReal rwork[1],tmp;
598 5390 PetscBLASInt i,j,one=1;
599 5390 PetscBool isrealpos=PETSC_TRUE;
600
601
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 5390 PetscFunctionBegin;
602
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 5390 if (n<SMALLN) { /* compute matrix power explicitly */
603
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
 4930 if (m==1) {
604 ✗ *nrm = LAPACKlange_("O",&n,&n,A,&n,rwork);
605 ✗ PetscCall(PetscLogFlops(1.0*n*n));
606 } else { /* m>=2 */
607
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.
 4930 PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n,&n,&n,&sone,A,&n,A,&n,&szero,v,&n));
608
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 75830 for (j=0;j<m-2;j++) {
609
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.
 70900 PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n,&n,&n,&sone,A,&n,v,&n,&szero,w,&n));
610 70900 SlepcSwap(v,w,aux);
611 }
612 4930 *nrm = LAPACKlange_("O",&n,&n,v,&n,rwork);
613
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.
 4930 PetscCall(PetscLogFlops(2.0*n*n*n*(m-1)+1.0*n*n));
614 }
615 } else {
616
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 55010 for (i=0;i<n;i++)
617
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 4421260 for (j=0;j<n;j++)
618 #if defined(PETSC_USE_COMPLEX)
619
3/4
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
✓ Branch 2 taken 5 times.
✗ Branch 3 not taken.
 2200145 if (PetscRealPart(A[i+j*n])<0.0 || PetscImaginaryPart(A[i+j*n])!=0.0) { isrealpos = PETSC_FALSE; break; }
620 #else
621
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
 2200145 if (A[i+j*n]<0.0) { isrealpos = PETSC_FALSE; break; }
622 #endif
623
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 460 if (isrealpos) { /* for positive matrices only */
624
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 21030 for (i=0;i<n;i++) v[i] = 1.0;
625
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 4380 for (j=0;j<m;j++) { /* w = A'*v */
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.
 4200 PetscCallBLAS("BLASgemv",BLASgemv_("C",&n,&n,&sone,A,&n,v,&one,&szero,w,&one));
627 4200 SlepcSwap(v,w,aux);
628 }
629
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(PetscLogFlops(2.0*n*n*m));
630 180 *nrm = 0.0;
631
4/4
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
✓ Branch 2 taken 10 times.
✓ Branch 3 taken 10 times.
 21030 for (i=0;i<n;i++) if ((tmp = PetscAbsScalar(v[i])) > *nrm) *nrm = tmp; /* norm(v,inf) */
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.
 280 } else PetscCall(SlepcNormEst1(n,A,m,work,rand,nrm));
633 }
634
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.
 1062 PetscFunctionReturn(PETSC_SUCCESS);
635 }
636