GCC Code Coverage Report

Directory: ./
File: src/sys/classes/ds/impls/ghiep/hz.c
Date: 2025-10-04 04:19:13
Exec Total Coverage
Lines: 204 217 94.0%
Functions: 4 4 100.0%
Branches: 287 468 61.3%
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 HZ iteration for generalized symmetric-indefinite eigenproblem.
12 Based on Matlab code from David Watkins.
13
14 References:
15
16 [1] D.S. Watkins, The Matrix Eigenvalue Problem: GR and Krylov Subspace
17 Methods, SIAM, 2007.
18
19 [2] M.A. Brebner, J. Grad, "Eigenvalues of Ax = lambda Bx for real
20 symmetric matrices A and B computed by reduction to pseudosymmetric
21 form and the HR process", Linear Alg. Appl. 43:99-118, 1982.
22 */
23
24 #include <slepc/private/dsimpl.h>
25 #include <slepcblaslapack.h>
26
27 /*
28 Sets up a 2-by-2 matrix to eliminate y in the vector [x y]'.
29 Transformation is rotator if sygn = 1 and hyperbolic if sygn = -1.
30 */
31 3550 static PetscErrorCode UnifiedRotation(PetscReal x,PetscReal y,PetscReal sygn,PetscReal *rot,PetscReal *rcond,PetscBool *swap)
32 {
33 3550 PetscReal nrm,c,s;
34
35
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 3550 PetscFunctionBegin;
36 3550 *swap = PETSC_FALSE;
37
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
 3550 if (y == 0) {
38 ✗ rot[0] = 1.0; rot[1] = 0.0; rot[2] = 0.0; rot[3] = 1.0;
39 ✗ *rcond = 1.0;
40 } else {
41
6/6
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
✓ Branch 2 taken 10 times.
✓ Branch 3 taken 10 times.
✓ Branch 4 taken 10 times.
✓ Branch 5 taken 10 times.
 3550 nrm = PetscMax(PetscAbs(x),PetscAbs(y));
42 3550 c = x/nrm; s = y/nrm;
43
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 3550 PetscCheck(sygn==1.0 || sygn==-1,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Value of sygn sent to transetup must be 1 or -1");
44
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 3550 if (sygn == 1.0) { /* set up a rotator */
45 2390 nrm = SlepcAbs(c,s);
46 2390 c = c/nrm; s = s/nrm;
47 /* rot = [c s; -s c]; */
48 2390 rot[0] = c; rot[1] = -s; rot[2] = s; rot[3] = c;
49 2390 *rcond = 1.0;
50 } else { /* sygn == -1, set up a hyperbolic transformation */
51 1160 nrm = c*c-s*s;
52
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 1160 if (nrm > 0) nrm = PetscSqrtReal(nrm);
53 else {
54
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 590 PetscCheck(nrm<0,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Breakdown in construction of hyperbolic transformation");
55 590 nrm = PetscSqrtReal(-nrm);
56 590 *swap = PETSC_TRUE;
57 }
58 1160 c = c/nrm; s = s/nrm;
59 /* rot = [c -s; -s c]; */
60 1160 rot[0] = c; rot[1] = -s; rot[2] = -s; rot[3] = c;
61
6/6
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
✓ Branch 2 taken 10 times.
✓ Branch 3 taken 10 times.
✓ Branch 4 taken 10 times.
✓ Branch 5 taken 10 times.
 1210 *rcond = PetscAbs(PetscAbs(s)-PetscAbs(c))/(PetscAbs(s)+PetscAbs(c));
62 }
63 }
64
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.
 710 PetscFunctionReturn(PETSC_SUCCESS);
65 }
66
67 370 static PetscErrorCode HZStep(PetscBLASInt ntop,PetscBLASInt nn,PetscReal tr,PetscReal dt,PetscReal *aa,PetscReal *bb,PetscReal *dd,PetscScalar *uu,PetscInt n,PetscInt ld,PetscBool *flag)
68 {
69 370 PetscBLASInt one=1;
70 370 PetscInt k,jj,ii;
71 370 PetscBLASInt n_;
72 370 PetscReal bulge10,bulge20,bulge30,bulge31,bulge41,bulge42;
73 370 PetscReal sygn,rcond=1.0,worstcond,rot[4],buf[2],t;
74 370 PetscScalar rtmp;
75 370 PetscBool swap;
76
77
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 370 PetscFunctionBegin;
78 370 *flag = PETSC_FALSE;
79 370 worstcond = 1.0;
80
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.
 370 PetscCall(PetscBLASIntCast(n,&n_));
81
82 /* Build initial bulge that sets step in motion */
83 370 bulge10 = dd[ntop+1]*(aa[ntop]*(aa[ntop] - dd[ntop]*tr) + dt*dd[ntop]*dd[ntop]) + dd[ntop]*bb[ntop]*bb[ntop];
84 370 bulge20 = bb[ntop]*(dd[ntop+1]*aa[ntop] + dd[ntop]*aa[ntop+1] - dd[ntop]*dd[ntop+1]*tr);
85 370 bulge30 = bb[ntop]*bb[ntop+1]*dd[ntop];
86 370 bulge31 = 0.0;
87 370 bulge41 = 0.0;
88 370 bulge42 = 0.0;
89
90 /* Chase the bulge */
91
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 2370 for (jj=ntop;jj<nn-1;jj++) {
92
93 /* Check for trivial bulge */
94
18/18
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
✓ Branch 2 taken 10 times.
✓ Branch 3 taken 10 times.
✓ Branch 4 taken 10 times.
✓ Branch 5 taken 10 times.
✓ Branch 6 taken 10 times.
✓ Branch 7 taken 10 times.
✓ Branch 8 taken 10 times.
✓ Branch 9 taken 10 times.
✓ Branch 10 taken 10 times.
✓ Branch 11 taken 10 times.
✓ Branch 12 taken 10 times.
✓ Branch 13 taken 10 times.
✓ Branch 14 taken 10 times.
✓ Branch 15 taken 10 times.
✓ Branch 16 taken 10 times.
✓ Branch 17 taken 10 times.
 4870 if (jj>ntop && PetscMax(PetscMax(PetscAbs(bulge10),PetscAbs(bulge20)),PetscAbs(bulge30))<PETSC_MACHINE_EPSILON*(PetscAbs(aa[jj]) + PetscAbs(aa[jj+1]))) {
95 80 bb[jj-1] = 0.0; /* deflate and move on */
96
97 } else { /* carry out the step */
98
99 /* Annihilate tip entry bulge30 */
100
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 1920 if (bulge30 != 0.0) {
101
102 /* Make an interchange if necessary to ensure that the
103 first transformation is othogonal, not hyperbolic. */
104
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 1630 if (dd[jj+1] != dd[jj+2]) { /* make an interchange */
105
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 450 if (dd[jj] != dd[jj+1]) { /* interchange 1st and 2nd */
106 200 buf[0] = bulge20; bulge20 = bulge10; bulge10 = buf[0];
107 200 buf[0] = aa[jj]; aa[jj] = aa[jj+1]; aa[jj+1] = buf[0];
108 200 buf[0] = bb[jj+1]; bb[jj+1] = bulge31; bulge31 = buf[0];
109 200 buf[0] = dd[jj]; dd[jj] = dd[jj+1]; dd[jj+1] = buf[0];
110
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 2200 for (k=0;k<n;k++) {
111 2000 rtmp = uu[k+jj*ld]; uu[k+jj*ld] = uu[k+(jj+1)*ld]; uu[k+(jj+1)*ld] = rtmp;
112 }
113 } else { /* interchange 1st and 3rd */
114 250 buf[0] = bulge30; bulge30 = bulge10; bulge10 = buf[0];
115 250 buf[0] = aa[jj]; aa[jj] = aa[jj+2]; aa[jj+2] = buf[0];
116 250 buf[0] = bb[jj]; bb[jj] = bb[jj+1]; bb[jj+1] = buf[0];
117 250 buf[0] = dd[jj]; dd[jj] = dd[jj+2]; dd[jj+2] = buf[0];
118
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 250 if (jj + 2 < nn-1) {
119 210 bulge41 = bb[jj+2];
120 210 bb[jj+2] = 0;
121 }
122
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 2750 for (k=0;k<n;k++) {
123 2500 rtmp = uu[k+jj*ld]; uu[k+jj*ld] = uu[k+(jj+2)*ld]; uu[k+(jj+2)*ld] = rtmp;
124 }
125 }
126 }
127
128 /* Set up transforming matrix rot. */
129
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 1630 PetscCall(UnifiedRotation(bulge20,bulge30,1,rot,&rcond,&swap));
130
131 /* Apply transforming matrix rot to T. */
132 1630 bulge20 = rot[0]*bulge20 + rot[2]*bulge30;
133 1630 buf[0] = rot[0]*bb[jj] + rot[2]*bulge31;
134 1630 buf[1] = rot[1]*bb[jj] + rot[3]*bulge31;
135 1630 bb[jj] = buf[0];
136 1630 bulge31 = buf[1];
137 1630 buf[0] = rot[0]*rot[0]*aa[jj+1] + 2.0*rot[0]*rot[2]*bb[jj+1] + rot[2]*rot[2]*aa[jj+2];
138 1630 buf[1] = rot[1]*rot[1]*aa[jj+1] + 2.0*rot[3]*rot[1]*bb[jj+1] + rot[3]*rot[3]*aa[jj+2];
139 1630 bb[jj+1] = rot[1]*rot[0]*aa[jj+1] + rot[3]*rot[2]*aa[jj+2] + (rot[3]*rot[0] + rot[1]*rot[2])*bb[jj+1];
140 1630 aa[jj+1] = buf[0];
141 1630 aa[jj+2] = buf[1];
142
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 1630 if (jj + 2 < nn-1) {
143 1260 bulge42 = bb[jj+2]*rot[2];
144 1260 bb[jj+2] = bb[jj+2]*rot[3];
145 }
146
147 /* Accumulate transforming matrix */
148
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.
 1630 PetscCallBLAS("BLASrot",BLASMIXEDrot_(&n_,uu+(jj+1)*ld,&one,uu+(jj+2)*ld,&one,&rot[0],&rot[2]));
149 }
150
151 /* Annihilate inner entry bulge20 */
152
1/2
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
 1920 if (bulge20 != 0.0) {
153
154 /* Begin by setting up transforming matrix rot */
155 1920 sygn = dd[jj]*dd[jj+1];
156
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.
 1920 PetscCall(UnifiedRotation(bulge10,bulge20,sygn,rot,&rcond,&swap));
157
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
 1920 if (rcond<PETSC_MACHINE_EPSILON) {
158 ✗ *flag = PETSC_TRUE;
159 ✗ PetscFunctionReturn(PETSC_SUCCESS);
160 }
161
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 1920 if (rcond < worstcond) worstcond = rcond;
162
163 /* Apply transforming matrix rot to T */
164
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 1920 if (jj > ntop) bb[jj-1] = rot[0]*bulge10 + rot[2]*bulge20;
165 1920 buf[0] = rot[0]*rot[0]*aa[jj] + 2*rot[0]*rot[2]*bb[jj] + rot[2]*rot[2]*aa[jj+1];
166 1920 buf[1] = rot[1]*rot[1]*aa[jj] + 2*rot[3]*rot[1]*bb[jj] + rot[3]*rot[3]*aa[jj+1];
167 1920 bb[jj] = rot[1]*rot[0]*aa[jj] + rot[3]*rot[2]*aa[jj+1] + (rot[3]*rot[0] + rot[1]*rot[2])*bb[jj];
168 1920 aa[jj] = buf[0];
169 1920 aa[jj+1] = buf[1];
170
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 1920 if (jj + 1 < nn-1) {
171 /* buf = [ bulge31 bb(jj+1) ] * rot' */
172 1630 buf[0] = rot[0]*bulge31 + rot[2]*bb[jj+1];
173 1630 buf[1] = rot[1]*bulge31 + rot[3]*bb[jj+1];
174 1630 bulge31 = buf[0];
175 1630 bb[jj+1] = buf[1];
176 }
177
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 1920 if (jj + 2 < nn-1) {
178 /* buf = [bulge41 bulge42] * rot' */
179 1260 buf[0] = rot[0]*bulge41 + rot[2]*bulge42;
180 1260 buf[1] = rot[1]*bulge41 + rot[3]*bulge42;
181 1260 bulge41 = buf[0];
182 1260 bulge42 = buf[1];
183 }
184
185 /* Apply transforming matrix rot to D */
186
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 1920 if (swap == 1) {
187 590 buf[0] = dd[jj]; dd[jj] = dd[jj+1]; dd[jj+1] = buf[0];
188 }
189
190 /* Accumulate transforming matrix, uu(jj:jj+1,:) = rot*uu(jj:jj+1,:) */
191
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 1920 if (sygn==1) {
192
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("BLASrot",BLASMIXEDrot_(&n_,uu+jj*ld,&one,uu+(jj+1)*ld,&one,&rot[0],&rot[2]));
193 } else {
194
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 1160 if (PetscAbsReal(rot[0])>PetscAbsReal(rot[1])) { /* Type I */
195 570 t = rot[1]/rot[0];
196
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 6270 for (ii=0;ii<n;ii++) {
197 5700 uu[jj*ld+ii] = rot[0]*uu[jj*ld+ii] + rot[1]*uu[(jj+1)*ld+ii];
198 5700 uu[(jj+1)*ld+ii] = t*uu[jj*ld+ii] + uu[(jj+1)*ld+ii]/rot[0];
199 }
200 } else { /* Type II */
201 590 t = rot[0]/rot[1];
202
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 6490 for (ii=0;ii<n;ii++) {
203 5900 rtmp = uu[jj*ld+ii];
204 5900 uu[jj*ld+ii] = rot[0]*uu[jj*ld+ii] + rot[1]*uu[(jj+1)*ld+ii];
205 5900 uu[(jj+1)*ld+ii] = t*uu[jj*ld+ii] + rtmp/rot[1];
206 }
207 }
208 }
209 }
210 }
211
212 /* Adjust bulge for next step */
213 2000 bulge10 = bb[jj];
214 2000 bulge20 = bulge31;
215 2000 bulge30 = bulge41;
216 2000 bulge31 = bulge42;
217 2000 bulge41 = 0.0;
218 2000 bulge42 = 0.0;
219 }
220
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.
 74 PetscFunctionReturn(PETSC_SUCCESS);
221 }
222
223 30 static PetscErrorCode HZIteration(PetscBLASInt nn,PetscBLASInt cgd,PetscReal *aa,PetscReal *bb,PetscReal *dd,PetscScalar *uu,PetscBLASInt ld)
224 {
225 30 PetscBLASInt j2,one=1,its,nits,nstop,jj,ntop,nbot,ntry;
226 30 PetscReal htr,det,dis,dif,tn,kt,c,s,tr,dt;
227 30 PetscBool flag=PETSC_FALSE;
228
229
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 30 PetscFunctionBegin;
230 30 its = 0;
231 30 nbot = nn-1;
232 30 nits = 0;
233 30 nstop = 40*(nn - cgd);
234
235
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 570 while (nbot >= cgd && nits < nstop) {
236
237 /* Check for zeros on the subdiagonal */
238 540 jj = nbot - 1;
239
10/10
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
✓ Branch 2 taken 10 times.
✓ Branch 3 taken 10 times.
✓ Branch 4 taken 10 times.
✓ Branch 5 taken 10 times.
✓ Branch 6 taken 10 times.
✓ Branch 7 taken 10 times.
✓ Branch 8 taken 10 times.
✓ Branch 9 taken 10 times.
 2650 while (jj>=cgd && PetscAbs(bb[jj])>PETSC_MACHINE_EPSILON*(PetscAbs(aa[jj])+PetscAbs(aa[jj+1]))) jj = jj-1;
240
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 540 if (jj>=cgd) bb[jj]=0;
241 540 ntop = jj + 1; /* starting point for step */
242
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 540 if (ntop == nbot) { /* isolate single eigenvalue */
243 nbot = ntop - 1;
244 its = 0;
245
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 480 } else if (ntop+1 == nbot) { /* isolate pair of eigenvalues */
246 110 htr = 0.5*(aa[ntop]*dd[ntop] + aa[nbot]*dd[nbot]);
247 110 det = dd[ntop]*dd[nbot]*(aa[ntop]*aa[nbot]-bb[ntop]*bb[ntop]);
248 110 dis = htr*htr - det;
249
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 110 if (dis > 0) { /* distinct real eigenvalues */
250
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 80 if (dd[ntop] == dd[nbot]) { /* separate the eigenvalues by a Jacobi rotator */
251 60 dif = aa[ntop]-aa[nbot];
252
4/4
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
✓ Branch 2 taken 10 times.
✓ Branch 3 taken 10 times.
 60 if (2.0*PetscAbs(bb[ntop])<=dif) {
253 30 tn = 2*bb[ntop]/dif;
254 30 tn = tn/(1.0 + PetscSqrtReal(1.0+tn*tn));
255 } else {
256 30 kt = dif/(2.0*bb[ntop]);
257
3/4
✓ Branch 0 taken 6 times.
✓ Branch 1 taken 10 times.
✓ Branch 2 taken 6 times.
✗ Branch 3 not taken.
 42 tn = PetscSign(kt)/(PetscAbsReal(kt)+PetscSqrtReal(1.0+kt*kt));
258 }
259 60 c = 1.0/PetscSqrtReal(1.0 + tn*tn);
260 60 s = c*tn;
261 60 aa[ntop] = aa[ntop] + tn*bb[ntop];
262 60 aa[nbot] = aa[nbot] - tn*bb[ntop];
263 60 bb[ntop] = 0;
264 60 j2 = nn-cgd;
265
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.
 60 PetscCallBLAS("BLASrot",BLASMIXEDrot_(&j2,uu+ntop*ld+cgd,&one,uu+nbot*ld+cgd,&one,&c,&s));
266 }
267 }
268 nbot = ntop - 1;
269 } else { /* Do an HZ iteration */
270 370 its = its + 1;
271 370 nits = nits + 1;
272 370 tr = aa[nbot-1]*dd[nbot-1] + aa[nbot]*dd[nbot];
273 370 dt = dd[nbot-1]*dd[nbot]*(aa[nbot-1]*aa[nbot]-bb[nbot-1]*bb[nbot-1]);
274
1/2
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
 370 for (ntry=1;ntry<=6;ntry++) {
275
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.
 370 PetscCall(HZStep(ntop,nbot+1,tr,dt,aa,bb,dd,uu,nn,ld,&flag));
276
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
 370 if (!flag) break;
277 ✗ PetscCheck(ntry<6,PETSC_COMM_SELF,PETSC_ERR_CONV_FAILED,"Unable to complete hz step after six tries");
278 ✗ tr = 0.9*tr; dt = 0.81*dt;
279 }
280 }
281 }
282
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.
 6 PetscFunctionReturn(PETSC_SUCCESS);
283 }
284
285 40 PetscErrorCode DSSolve_GHIEP_HZ(DS ds,PetscScalar *wr,PetscScalar *wi)
286 {
287 40 PetscInt i,off;
288 40 PetscBLASInt n,l,n1,ld = 0;
289 40 const PetscScalar *A,*B;
290 40 PetscScalar *Q;
291 40 PetscReal *d,*e,*s;
292
293
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 40 PetscFunctionBegin;
294 #if !defined(PETSC_USE_COMPLEX)
295
2/8
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 20 PetscAssertPointer(wi,3);
296 #endif
297
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(PetscBLASIntCast(ds->ld,&ld));
298
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(PetscBLASIntCast(ds->n-ds->l,&n1));
299 40 off = ds->l + ds->l*ld;
300
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(DSGetArrayReal(ds,DS_MAT_T,&d));
301
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(DSGetArrayReal(ds,DS_MAT_D,&s));
302 40 e = d + ld;
303 #if defined(PETSC_USE_DEBUG)
304 /* Check signature */
305
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.
 8 PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_B],&B));
306
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
 78 for (i=0;i<ds->n;i++) {
307
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
 70 PetscReal de = (ds->compact)?s[i]:PetscRealPart(B[i*ld+i]);
308
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 70 PetscCheck(de==1.0 || de==-1.0,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Diagonal elements of the signature matrix must be 1 or -1");
309 }
310
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.
 8 PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_B],&B));
311 #endif
312 /* Quick return */
313
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 40 if (n1 == 1) {
314
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.
 10 PetscCall(MatDenseGetArray(ds->omat[DS_MAT_Q],&Q));
315
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 60 for (i=0;i<=ds->l;i++) Q[i+i*ld] = 1.0;
316
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.
 10 PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_Q],&Q));
317
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.
 10 PetscCall(DSGHIEPComplexEigs(ds,0,ds->l,wr,wi));
318
1/2
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
 10 if (ds->compact) {
319 10 wr[ds->l] = d[ds->l]/s[ds->l]; wi[ds->l] = 0.0;
320 } else {
321 ✗ PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_A],&A));
322 ✗ PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_B],&B));
323 ✗ d[ds->l] = PetscRealPart(A[off]);
324 ✗ s[ds->l] = PetscRealPart(B[off]);
325 ✗ PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_A],&A));
326 ✗ PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_B],&B));
327 ✗ wr[ds->l] = d[ds->l]/s[ds->l]; wi[ds->l] = 0.0;
328 }
329
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.
 10 PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&d));
330
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.
 10 PetscCall(DSRestoreArrayReal(ds,DS_MAT_D,&s));
331
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.
 2 PetscFunctionReturn(PETSC_SUCCESS);
332 }
333 /* Reduce to pseudotriadiagonal form */
334
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 30 PetscCall(DSIntermediate_GHIEP(ds));
335
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.
 30 PetscCall(MatDenseGetArray(ds->omat[DS_MAT_Q],&Q));
336
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.
 30 PetscCall(PetscBLASIntCast(ds->n,&n));
337
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.
 30 PetscCall(PetscBLASIntCast(ds->l,&l));
338
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.
 30 PetscCall(HZIteration(n,l,d,e,s,Q,ld));
339
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.
 30 PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_Q],&Q));
340
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.
 30 PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&d));
341
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.
 30 PetscCall(DSRestoreArrayReal(ds,DS_MAT_D,&s));
342
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.
 30 if (!ds->compact) PetscCall(DSSwitchFormat_GHIEP(ds,PETSC_FALSE));
343 /* Undo from diagonal the blocks with real eigenvalues*/
344
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 30 PetscCall(DSGHIEPRealBlocks(ds));
345
346 /* Recover eigenvalues from diagonal */
347
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.
 30 PetscCall(DSGHIEPComplexEigs(ds,0,ds->n,wr,wi));
348 #if defined(PETSC_USE_COMPLEX)
349
1/2
✓ Branch 0 taken 5 times.
✗ Branch 1 not taken.
 15 if (wi) {
350
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
 155 for (i=ds->l;i<ds->n;i++) wi[i] = 0.0;
351 }
352 #endif
353 30 ds->t = ds->n;
354
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.
 30 PetscFunctionReturn(PETSC_SUCCESS);
355 }
356