GCC Code Coverage Report

Directory: ./
File: src/sys/classes/ds/impls/hep/bdc/dsbtdc.c
Date: 2026-07-29 03:58:07
Exec Total Coverage
Lines: 132 155 85.2%
Functions: 1 1 100.0%
Branches: 158 366 43.2%
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 BDC - Block-divide and conquer (see description in README file)
12 */
13
14 #include <slepc/private/dsimpl.h>
15 #include <slepcblaslapack.h>
16
17 10 PetscErrorCode BDC_dsbtdc_(const char *jobz,const char *jobacc,PetscBLASInt n,
18 PetscBLASInt nblks,PetscBLASInt *ksizes,PetscReal *d,PetscBLASInt l1d,
19 PetscBLASInt l2d,PetscReal *e,PetscBLASInt l1e,PetscBLASInt l2e,PetscReal tol,
20 PetscReal tau1,PetscReal tau2,PetscReal *ev,PetscReal *z,PetscBLASInt ldz,
21 PetscReal *work,PetscBLASInt lwork,PetscBLASInt *iwork,PetscBLASInt liwork,
22 PetscReal *mingap,PetscBLASInt *mingapi,PetscBLASInt *oinfo,
23 PetscBLASInt jobz_len,PetscBLASInt jobacc_len)
24 {
25 /* -- Routine written in LAPACK Version 3.0 style -- */
26 /* *************************************************** */
27 /* Written by */
28 /* Michael Moldaschl and Wilfried Gansterer */
29 /* University of Vienna */
30 /* last modification: March 28, 2014 */
31
32 /* Small adaptations of original code written by */
33 /* Wilfried Gansterer and Bob Ward, */
34 /* Department of Computer Science, University of Tennessee */
35 /* see https://doi.org/10.1137/S1064827501399432 */
36 /* *************************************************** */
37
38 /* Purpose */
39 /* ======= */
40
41 /* DSBTDC computes approximations to all eigenvalues and eigenvectors */
42 /* of a symmetric block tridiagonal matrix using the divide and */
43 /* conquer method with lower rank approximations to the subdiagonal blocks. */
44
45 /* This code makes very mild assumptions about floating point */
46 /* arithmetic. It will work on machines with a guard digit in */
47 /* add/subtract, or on those binary machines without guard digits */
48 /* which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2. */
49 /* It could conceivably fail on hexadecimal or decimal machines */
50 /* without guard digits, but we know of none. See DLAED3M for details. */
51
52 /* Arguments */
53 /* ========= */
54
55 /* JOBZ (input) CHARACTER*1 */
56 /* = 'N': Compute eigenvalues only (not implemented); */
57 /* = 'D': Compute eigenvalues and eigenvectors. Eigenvectors */
58 /* are accumulated in the divide-and-conquer process. */
59
60 /* JOBACC (input) CHARACTER*1 */
61 /* = 'A' ("automatic"): The accuracy parameters TAU1 and TAU2 */
62 /* are determined automatically from the */
63 /* parameter TOL according to the analytical */
64 /* bounds. In that case the input values of */
65 /* TAU1 and TAU2 are irrelevant (ignored). */
66 /* = 'M' ("manual"): The input values of the accuracy parameters */
67 /* TAU1 and TAU2 are used. In that case the input */
68 /* value of the parameter TOL is irrelevant */
69 /* (ignored). */
70
71 /* N (input) INTEGER */
72 /* The dimension of the symmetric block tridiagonal matrix. */
73 /* N >= 1. */
74
75 /* NBLKS (input) INTEGER, 1 <= NBLKS <= N */
76 /* The number of diagonal blocks in the matrix. */
77
78 /* KSIZES (input) INTEGER array, dimension (NBLKS) */
79 /* The dimensions of the square diagonal blocks from top left */
80 /* to bottom right. KSIZES(I) >= 1 for all I, and the sum of */
81 /* KSIZES(I) for I = 1 to NBLKS has to be equal to N. */
82
83 /* D (input) DOUBLE PRECISION array, dimension (L1D,L2D,NBLKS) */
84 /* The lower triangular elements of the symmetric diagonal */
85 /* blocks of the block tridiagonal matrix. The elements of the top */
86 /* left diagonal block, which is of dimension KSIZES(1), have to */
87 /* be placed in D(*,*,1); the elements of the next diagonal */
88 /* block, which is of dimension KSIZES(2), have to be placed in */
89 /* D(*,*,2); etc. */
90
91 /* L1D (input) INTEGER */
92 /* The leading dimension of the array D. L1D >= max(3,KMAX), */
93 /* where KMAX is the dimension of the largest diagonal block, */
94 /* i.e., KMAX = max_I (KSIZES(I)). */
95
96 /* L2D (input) INTEGER */
97 /* The second dimension of the array D. L2D >= max(3,KMAX), */
98 /* where KMAX is as stated in L1D above. */
99
100 /* E (input) DOUBLE PRECISION array, dimension (L1E,L2E,NBLKS-1) */
101 /* The elements of the subdiagonal blocks of the */
102 /* block tridiagonal matrix. The elements of the top left */
103 /* subdiagonal block, which is KSIZES(2) x KSIZES(1), have to be */
104 /* placed in E(*,*,1); the elements of the next subdiagonal block, */
105 /* which is KSIZES(3) x KSIZES(2), have to be placed in E(*,*,2); etc. */
106 /* During runtime, the original contents of E(*,*,K) is */
107 /* overwritten by the singular vectors and singular values of */
108 /* the lower rank representation. */
109
110 /* L1E (input) INTEGER */
111 /* The leading dimension of the array E. L1E >= max(3,2*KMAX+1), */
112 /* where KMAX is as stated in L1D above. The size of L1E enables */
113 /* the storage of ALL singular vectors and singular values for */
114 /* the corresponding off-diagonal block in E(*,*,K) and therefore */
115 /* there are no restrictions on the rank of the approximation */
116 /* (only the "natural" restriction */
117 /* RANK(K) .LE. MIN(KSIZES(K),KSIZES(K+1))). */
118
119 /* L2E (input) INTEGER */
120 /* The second dimension of the array E. L2E >= max(3,2*KMAX+1), */
121 /* where KMAX is as stated in L1D above. The size of L2E enables */
122 /* the storage of ALL singular vectors and singular values for */
123 /* the corresponding off-diagonal block in E(*,*,K) and therefore */
124 /* there are no restrictions on the rank of the approximation */
125 /* (only the "natural" restriction */
126 /* RANK(K) .LE. MIN(KSIZES(K),KSIZES(K+1))). */
127
128 /* TOL (input) DOUBLE PRECISION, TOL.LE.TOLMAX */
129 /* User specified tolerance for the residuals of the computed */
130 /* eigenpairs. If (JOBACC.EQ.'A') then it is used to determine */
131 /* TAU1 and TAU2; ignored otherwise. */
132 /* If (TOL.LT.40*EPS .AND. JOBACC.EQ.'A') then TAU1 is set to machine */
133 /* epsilon and TAU2 is set to the standard deflation tolerance from */
134 /* LAPACK. */
135
136 /* TAU1 (input) DOUBLE PRECISION, TAU1.LE.TOLMAX/2 */
137 /* User specified tolerance for determining the rank of the */
138 /* lower rank approximations to the off-diagonal blocks. */
139 /* The rank for each off-diagonal block is determined such that */
140 /* the resulting absolute eigenvalue error is less than or equal */
141 /* to TAU1. */
142 /* If (JOBACC.EQ.'A') then TAU1 is determined automatically from */
143 /* TOL and the input value is ignored. */
144 /* If (JOBACC.EQ.'M' .AND. TAU1.LT.20*EPS) then TAU1 is set to */
145 /* machine epsilon. */
146
147 /* TAU2 (input) DOUBLE PRECISION, TAU2.LE.TOLMAX/2 */
148 /* User specified deflation tolerance for the routine DIBTDC. */
149 /* If (1.0D-1.GT.TAU2.GT.20*EPS) then TAU2 is used as */
150 /* the deflation tolerance in DSRTDF (EPS is the machine epsilon). */
151 /* If (TAU2.LE.20*EPS) then the standard deflation tolerance from */
152 /* LAPACK is used as the deflation tolerance in DSRTDF. */
153 /* If (JOBACC.EQ.'A') then TAU2 is determined automatically from */
154 /* TOL and the input value is ignored. */
155 /* If (JOBACC.EQ.'M' .AND. TAU2.LT.20*EPS) then TAU2 is set to */
156 /* the standard deflation tolerance from LAPACK. */
157
158 /* EV (output) DOUBLE PRECISION array, dimension (N) */
159 /* If INFO = 0, then EV contains the computed eigenvalues of the */
160 /* symmetric block tridiagonal matrix in ascending order. */
161
162 /* Z (output) DOUBLE PRECISION array, dimension (LDZ,N) */
163 /* If (JOBZ.EQ.'D' .AND. INFO = 0) */
164 /* then Z contains the orthonormal eigenvectors of the symmetric */
165 /* block tridiagonal matrix computed by the routine DIBTDC */
166 /* (accumulated in the divide-and-conquer process). */
167 /* If (-199 < INFO < -99) then Z contains the orthonormal */
168 /* eigenvectors of the symmetric block tridiagonal matrix, */
169 /* computed without divide-and-conquer (quick returns). */
170 /* Otherwise not referenced. */
171
172 /* LDZ (input) INTEGER */
173 /* The leading dimension of the array Z. LDZ >= max(1,N). */
174
175 /* WORK (workspace/output) DOUBLE PRECISION array, dimension (LWORK) */
176
177 /* LWORK (input) INTEGER */
178 /* The dimension of the array WORK. */
179 /* If NBLKS.EQ.1, then LWORK has to be at least 2N^2+6N+1 */
180 /* (for the call of DSYEVD). */
181 /* If NBLKS.GE.2 and (JOBZ.EQ.'D') then the absolute minimum */
182 /* required for DIBTDC is (N**2 + 3*N). This will not always */
183 /* suffice, though, the routine will return a corresponding */
184 /* error code and report how much work space was missing (see */
185 /* INFO). */
186 /* In order to guarantee correct results in all cases where */
187 /* NBLKS.GE.2, LWORK must be at least (2*N**2 + 3*N). */
188
189 /* IWORK (workspace/output) INTEGER array, dimension (LIWORK) */
190
191 /* LIWORK (input) INTEGER */
192 /* The dimension of the array IWORK. */
193 /* LIWORK must be at least (5*N + 5*NBLKS - 1) (for DIBTDC) */
194 /* Note that this should also suffice for the call of DSYEVD on a */
195 /* diagonal block which requires (5*KMAX + 3). */
196
197 /* MINGAP (output) DOUBLE PRECISION */
198 /* The minimum "gap" between the approximate eigenvalues */
199 /* computed, i.e., MIN( ABS(EV(I+1)-EV(I)) for I=1,2,..., N-1 */
200 /* IF (MINGAP.LE.TOL/10) THEN a warning flag is returned in INFO, */
201 /* because the computed eigenvectors may be unreliable individually */
202 /* (only the subspaces spanned are approximated reliably). */
203
204 /* MINGAPI (output) INTEGER */
205 /* Index I where the minimum gap in the spectrum occurred. */
206
207 /* INFO (output) INTEGER */
208 /* = 0: successful exit, no special cases occurred. */
209 /* < -200: not enough workspace. Space for ABS(INFO + 200) */
210 /* numbers is required in addition to the workspace provided, */
211 /* otherwise some of the computed eigenvectors will be incorrect. */
212 /* < -99, > -199: successful exit, but quick returns. */
213 /* if INFO = -100, successful exit, but the input matrix */
214 /* was the zero matrix and no */
215 /* divide-and-conquer was performed */
216 /* if INFO = -101, successful exit, but N was 1 and no */
217 /* divide-and-conquer was performed */
218 /* if INFO = -102, successful exit, but only a single */
219 /* dense block. Standard dense solver */
220 /* was called, no divide-and-conquer was */
221 /* performed */
222 /* if INFO = -103, successful exit, but warning that */
223 /* MINGAP.LE.TOL/10 and therefore the */
224 /* eigenvectors corresponding to close */
225 /* approximate eigenvalues may individually */
226 /* be unreliable (although taken together they */
227 /* do approximate the corresponding subspace to */
228 /* the desired accuracy) */
229 /* = -99: error in the preprocessing in DIBTDC (when determining */
230 /* the merging order). */
231 /* < 0, > -99: illegal arguments. */
232 /* if INFO = -i, the i-th argument had an illegal value. */
233 /* > 0: The algorithm failed to compute an eigenvalue while */
234 /* working on the submatrix lying in rows and columns */
235 /* INFO/(N+1) through mod(INFO,N+1). */
236
237 /* Further Details */
238 /* =============== */
239
240 /* Small modifications of code written by */
241 /* Wilfried Gansterer and Bob Ward, */
242 /* Department of Computer Science, University of Tennessee */
243 /* see https://doi.org/10.1137/S1064827501399432 */
244
245 /* Based on the design of the LAPACK code sstedc.f written by Jeff */
246 /* Rutter, Computer Science Division, University of California at */
247 /* Berkeley, and modified by Francoise Tisseur, University of Tennessee. */
248
249 /* ===================================================================== */
250
251 /* .. Parameters .. */
252
253 #define TOLMAX 0.1
254
255 /* TOLMAX .... upper bound for tolerances TOL, TAU1, TAU2 */
256 /* NOTE: in the routine DIBTDC, the value */
257 /* 1.D-1 is hardcoded for TOLMAX ! */
258
259 10 PetscBLASInt i, j, k, i1, iwspc, lwmin, start;
260 10 PetscBLASInt ii, ip, nk, rk, np, iu, rp1, ldu;
261 10 PetscBLASInt ksk, ivt, iend, kchk=0, kmax=0, one=1, zero=0;
262 10 PetscBLASInt ldvt, ksum=0, kskp1, spneed, nrblks, liwmin, isvals;
263 10 PetscReal p, d2, eps, dmax, emax, done = 1.0;
264 10 PetscReal dnrm, tiny, anorm, exdnrm=0, dropsv, absdiff;
265
266
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 10 PetscFunctionBegin;
267 /* Determine machine epsilon. */
268 10 eps = LAPACKlamch_("Epsilon");
269
270 10 *oinfo = 0;
271
272
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
 10 if (*(unsigned char *)jobz != 'N' && *(unsigned char *)jobz != 'D') *oinfo = -1;
273
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
 10 else if (*(unsigned char *)jobacc != 'A' && *(unsigned char *)jobacc != 'M') *oinfo = -2;
274
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
 10 else if (n < 1) *oinfo = -3;
275
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
 10 else if (nblks < 1 || nblks > n) *oinfo = -4;
276
1/2
✓ Branch 0 taken 5 times.
✗ Branch 1 not taken.
 10 if (*oinfo == 0) {
277
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
 55 for (k = 0; k < nblks; ++k) {
278 45 ksk = ksizes[k];
279 45 ksum += ksk;
280 45 if (ksk > kmax) kmax = ksk;
281
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
 45 if (ksk < 1) kchk = 1;
282 }
283
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
 10 if (nblks == 1) lwmin = 2*n*n + n*6 + 1;
284 5 else lwmin = n*n + n*3;
285 10 liwmin = n * 5 + nblks * 5 - 4;
286
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
 10 if (ksum != n || kchk == 1) *oinfo = -5;
287
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
 10 else if (l1d < PetscMax(3,kmax)) *oinfo = -7;
288
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
 10 else if (l2d < PetscMax(3,kmax)) *oinfo = -8;
289
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
 10 else if (l1e < PetscMax(3,2*kmax+1)) *oinfo = -10;
290
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
 10 else if (l2e < PetscMax(3,2*kmax+1)) *oinfo = -11;
291
2/4
✓ Branch 0 taken 5 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 5 times.
 10 else if (*(unsigned char *)jobacc == 'A' && tol > TOLMAX) *oinfo = -12;
292
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 10 else if (*(unsigned char *)jobacc == 'M' && tau1 > TOLMAX/2) *oinfo = -13;
293
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 10 else if (*(unsigned char *)jobacc == 'M' && tau2 > TOLMAX/2) *oinfo = -14;
294
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
 10 else if (ldz < PetscMax(1,n)) *oinfo = -17;
295
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
 10 else if (lwork < lwmin) *oinfo = -19;
296
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
 10 else if (liwork < liwmin) *oinfo = -21;
297 }
298
299
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 10 PetscCheck(!*oinfo,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Wrong argument %" PetscBLASInt_FMT " in DSBTDC",-(*oinfo));
300
301 /* Quick return if possible */
302
303
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
 10 if (n == 1) {
304 ✗ ev[0] = d[0]; z[0] = 1.;
305 ✗ *oinfo = -101;
306 ✗ PetscFunctionReturn(PETSC_SUCCESS);
307 }
308
309 /* If NBLKS is equal to 1, then solve the problem with standard */
310 /* dense solver (in this case KSIZES(1) = N). */
311
312
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
 10 if (nblks == 1) {
313
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
 20 for (i = 0; i < n; ++i) {
314
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
 45 for (j = 0; j <= i; ++j) {
315 30 z[i + j*ldz] = d[i + j*l1d];
316 }
317 }
318
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.
 5 PetscCallLAPACKInfo("LAPACKsyevd",LAPACKsyevd_("V", "L", &n, z, &ldz, ev, work, &lwork, iwork, &liwork, &info));
319 5 *oinfo = -102;
320
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.
 5 PetscFunctionReturn(PETSC_SUCCESS);
321 }
322
323 /* determine the accuracy parameters (if requested) */
324
325
1/2
✓ Branch 0 taken 5 times.
✗ Branch 1 not taken.
 5 if (*(unsigned char *)jobacc == 'A') {
326 5 tau1 = tol / 2;
327
1/2
✓ Branch 0 taken 5 times.
✗ Branch 1 not taken.
 5 if (tau1 < eps * 20) tau1 = eps;
328 tau2 = tol / 2;
329 }
330
331 /* Initialize Z as the identity matrix */
332
333
1/2
✓ Branch 0 taken 5 times.
✗ Branch 1 not taken.
 5 if (*(unsigned char *)jobz == 'D') {
334
4/4
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
✓ Branch 2 taken 5 times.
✓ Branch 3 taken 5 times.
 3005 for (j=0;j<n;j++) for (i=0;i<n;i++) z[i+j*ldz] = 0.0;
335
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
 125 for (i=0;i<n;i++) z[i+i*ldz] = 1.0;
336 }
337
338 /* Determine the off-diagonal ranks, form and store the lower rank */
339 /* approximations based on the tolerance parameters, the */
340 /* RANK(K) largest singular values and the associated singular */
341 /* vectors of each subdiagonal block. Also find the maximum norm of */
342 /* the subdiagonal blocks (in EMAX). */
343
344 /* Compute SVDs of the subdiagonal blocks.... */
345
346 /* EMAX .... maximum norm of the off-diagonal blocks */
347
348 emax = 0.;
349
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
 40 for (k = 0; k < nblks-1; ++k) {
350 35 ksk = ksizes[k];
351 35 kskp1 = ksizes[k+1];
352 35 isvals = 0;
353
354 /* Note that min(KSKP1,KSK).LE.N/2 (equality possible for */
355 /* NBLKS=2), and therefore storing the singular values requires */
356 /* at most N/2 entries of the * array WORK. */
357
358 35 iu = isvals + n / 2;
359 35 ivt = isvals + n / 2;
360
361 /* Call of DGESVD: The space for U is not referenced, since */
362 /* JOBU='O' and therefore this portion of the array WORK */
363 /* is not referenced for U. */
364
365 35 ldu = kskp1;
366 35 ldvt = PetscMin(kskp1,ksk);
367 35 iwspc = ivt + n * n / 2;
368
369 /* Note that the minimum workspace required for this call */
370 /* of DGESVD is: N/2 for storing the singular values + N**2/2 for */
371 /* storing V^T + 5*N/2 workspace = N**2/2 + 3*N. */
372
373 35 i1 = lwork - iwspc;
374
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.
 35 PetscCallLAPACKInfo("LAPACKgesvd",LAPACKgesvd_("O", "S", &kskp1, &ksk,
375 &e[k*l1e*l2e], &l1e, &work[isvals],
376 &work[iu], &ldu, &work[ivt], &ldvt, &work[iwspc], &i1, &info));
377
378 /* Note that after the return from DGESVD U is stored in */
379 /* E(*,*,K), and V^{\top} is stored in WORK(IVT, IVT+1, ....) */
380
381 /* determine the ranks RANK() for the approximations */
382
383 35 rk = PetscMin(ksk,kskp1);
384 35 L8:
385 35 dropsv = work[isvals - 1 + rk];
386
387
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
 35 if (dropsv * 2. <= tau1) {
388
389 /* the error caused by dropping singular value RK is */
390 /* small enough, try to reduce the rank by one more */
391
392 ✗ if (--rk > 0) goto L8;
393 ✗ else iwork[k] = 0;
394 } else {
395
396 /* the error caused by dropping singular value RK is */
397 /* too large already, RK is the rank required to achieve the */
398 /* desired accuracy */
399
400 35 iwork[k] = rk;
401 }
402
403 /* ************************************************************************** */
404
405 /* Store the first RANK(K) terms of the SVD of the current */
406 /* off-diagonal block. */
407 /* NOTE that here it is required that L1E, L2E >= 2*KMAX+1 in order */
408 /* to have enough space for storing singular vectors and values up */
409 /* to the full SVD of an off-diagonal block !!!! */
410
411 /* u1-u_RANK(K) is already contained in E(:,1:RANK(K),K) (as a */
412 /* result of the call of DGESVD !), the sigma1-sigmaK are to be */
413 /* stored in E(1:RANK(K),RANK(K)+1,K), and v1-v_RANK(K) are to be */
414 /* stored in E(:,RANK(K)+2:2*RANK(K)+1,K) */
415
416 35 rp1 = iwork[k];
417
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
 140 for (j = 0; j < iwork[k]; ++j) {
418
419 /* store sigma_J in E(J,RANK(K)+1,K) */
420
421 105 e[j + (rp1 + k*l2e)* l1e] = work[isvals + j];
422
423 /* update maximum norm of subdiagonal blocks */
424
425
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
 105 if (e[j + (rp1 + k*l2e)*l1e] > emax) {
426 5 emax = e[j + (rp1 + k*l2e)*l1e];
427 }
428
429 /* store v_J in E(:,RANK(K)+1+J,K) */
430 /* (note that WORK contains V^{\top} and therefore */
431 /* we need to read rowwise !) */
432
433
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
 420 for (i = 1; i <= ksk; ++i) {
434 315 e[i-1 + (rp1+j+1 + k*l2e)*l1e] = work[ivt+j + (i-1)*ldvt];
435 }
436 }
437
438 }
439
440 /* Compute the maximum norm of diagonal blocks and store the norm */
441 /* of each diagonal block in E(RP1,RP1,K) (after the singular values); */
442 /* store the norm of the last diagonal block in EXDNRM. */
443
444 /* DMAX .... maximum one-norm of the diagonal blocks */
445
446 dmax = 0.;
447
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
 45 for (k = 0; k < nblks; ++k) {
448 40 rp1 = iwork[k];
449
450 /* compute the one-norm of diagonal block K */
451
452 40 dnrm = LAPACKlansy_("1", "L", &ksizes[k], &d[k*l1d*l2d], &l1d, work);
453
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
 40 if (k+1 == nblks) exdnrm = dnrm;
454 35 else e[rp1 + (rp1 + k*l2e)*l1e] = dnrm;
455
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
 40 if (dnrm > dmax) dmax = dnrm;
456 }
457
458 /* Check for zero matrix. */
459
460
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
 5 if (emax == 0. && dmax == 0.) {
461 ✗ for (i = 0; i < n; ++i) ev[i] = 0.;
462 ✗ *oinfo = -100;
463 ✗ PetscFunctionReturn(PETSC_SUCCESS);
464 }
465
466 /* **************************************************************** */
467
468 /* ....Identify irreducible parts of the block tridiagonal matrix */
469 /* [while (START <= NBLKS)].... */
470
471 start = 0;
472 np = 0;
473 5 L10:
474
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
 10 if (start < nblks) {
475
476 /* Let IEND be the number of the next subdiagonal block such that */
477 /* its RANK is 0 or IEND = NBLKS if no such subdiagonal exists. */
478 /* The matrix identified by the elements between the diagonal block START */
479 /* and the diagonal block IEND constitutes an independent (irreducible) */
480 /* sub-problem. */
481
482 iend = start;
483
484 40 L20:
485
1/2
✓ Branch 0 taken 5 times.
✗ Branch 1 not taken.
 40 if (iend < nblks) {
486 40 rk = iwork[iend];
487
488 /* NOTE: if RANK(IEND).EQ.0 then decoupling happens due to */
489 /* reduced accuracy requirements ! (because in this case */
490 /* we would not merge the corresponding two diagonal blocks) */
491
492 /* NOTE: seems like any combination may potentially happen: */
493 /* (i) RANK = 0 but no decoupling due to small norm of */
494 /* off-diagonal block (corresponding diagonal blocks */
495 /* also have small norm) as well as */
496 /* (ii) RANK > 0 but decoupling due to small norm of */
497 /* off-diagonal block (corresponding diagonal blocks */
498 /* have very large norm) */
499 /* case (i) is ruled out by checking for RANK = 0 above */
500 /* (we decide to decouple all the time when the rank */
501 /* of an off-diagonal block is zero, independently of */
502 /* the norms of the corresponding diagonal blocks. */
503
504
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
 40 if (rk > 0) {
505
506 /* check for decoupling due to small norm of off-diagonal block */
507 /* (relative to the norms of the corresponding diagonal blocks) */
508
509
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
 35 if (iend == nblks-2) {
510 5 d2 = PetscSqrtReal(exdnrm);
511 } else {
512 30 d2 = PetscSqrtReal(e[iwork[iend+1] + (iwork[iend+1] + (iend+1)*l2e)*l1e]);
513 }
514
515 /* this definition of TINY is analogous to the definition */
516 /* in the tridiagonal divide&conquer (dstedc) */
517
518 35 tiny = eps * PetscSqrtReal(e[iwork[iend] + (iwork[iend] + iend*l2e)*l1e])*d2;
519
1/2
✓ Branch 0 taken 5 times.
✗ Branch 1 not taken.
 35 if (e[(iwork[iend] + iend*l2e)*l1e] > tiny) {
520
521 /* no decoupling due to small norm of off-diagonal block */
522
523 35 ++iend;
524 35 goto L20;
525 }
526 }
527 }
528
529 /* ....(Sub) Problem determined: between diagonal blocks */
530 /* START and IEND. Compute its size and solve it.... */
531
532 5 nrblks = iend - start + 1;
533
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
 5 if (nrblks == 1) {
534
535 /* Isolated problem is a single diagonal block */
536
537 ✗ nk = ksizes[start];
538
539 /* copy this isolated block into Z */
540
541 ✗ for (i = 0; i < nk; ++i) {
542 ✗ ip = np + i + 1;
543 ✗ for (j = 0; j <= i; ++j) z[ip + (np+j+1)*ldz] = d[i + (j + start*l2d)*l1d];
544 }
545
546 /* check whether there is enough workspace */
547
548 ✗ spneed = 2*nk*nk + nk * 6 + 1;
549 ✗ PetscCheck(spneed<=lwork,PETSC_COMM_SELF,PETSC_ERR_MEM,"dsbtdc: not enough workspace for DSYEVD, info = %" PetscBLASInt_FMT,lwork - 200 - spneed);
550
551 ✗ PetscCallLAPACKInfo("LAPACKsyevd",LAPACKsyevd_("V", "L", &nk,
552 &z[np + np*ldz], &ldz, &ev[np],
553 work, &lwork, &iwork[nblks-1], &liwork, &info));
554 ✗ start = iend + 1;
555 ✗ np += nk;
556
557 /* go to the next irreducible subproblem */
558
559 ✗ goto L10;
560 }
561
562 /* ....Isolated problem consists of more than one diagonal block. */
563 /* Start the divide and conquer algorithm.... */
564
565 /* Scale: Divide by the maximum of all norms of diagonal blocks */
566 /* and singular values of the subdiagonal blocks */
567
568 /* ....determine maximum of the norms of all diagonal and subdiagonal */
569 /* blocks.... */
570
571
1/2
✓ Branch 0 taken 5 times.
✗ Branch 1 not taken.
 5 if (iend == nblks-1) anorm = exdnrm;
572 ✗ else anorm = e[iwork[iend] + (iwork[iend] + iend*l2e)*l1e];
573
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
 40 for (k = start; k < iend; ++k) {
574 35 rp1 = iwork[k];
575
576 /* norm of diagonal block */
577
1/2
✓ Branch 0 taken 5 times.
✗ Branch 1 not taken.
 35 anorm = PetscMax(anorm,e[rp1 + (rp1 + k*l2e)*l1e]);
578
579 /* singular value of subdiagonal block */
580
1/2
✓ Branch 0 taken 5 times.
✗ Branch 1 not taken.
 70 anorm = PetscMax(anorm,e[(rp1 + k*l2e)*l1e]);
581 }
582
583 5 nk = 0;
584
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
 45 for (k = start; k < iend+1; ++k) {
585 40 ksk = ksizes[k];
586 40 nk += ksk;
587
588 /* scale the diagonal block */
589
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.
 40 PetscCallLAPACKInfo("LAPACKlascl",LAPACKlascl_("L", &zero, &zero,
590 &anorm, &done, &ksk, &ksk, &d[k*l2d*l1d], &l1d, &info));
591
592 /* scale the (approximated) off-diagonal block by dividing its */
593 /* singular values */
594
595
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
 40 if (k != iend) {
596
597 /* the last subdiagonal block has index IEND-1 !!!! */
598
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
 140 for (i = 0; i < iwork[k]; ++i) {
599 105 e[i + (iwork[k] + k*l2e)*l1e] /= anorm;
600 }
601 }
602 }
603
604 /* call the block-tridiagonal divide-and-conquer on the */
605 /* irreducible subproblem which has been identified */
606
607
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.
 5 PetscCall(BDC_dibtdc_(jobz, nk, nrblks, &ksizes[start], &d[start*l1d*l2d], l1d, l2d,
608 &e[start*l2e*l1e], &iwork[start], l1e, l2e, tau2, &ev[np],
609 &z[np + np*ldz], ldz, work, lwork, &iwork[nblks-1], liwork, oinfo, 1));
610
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 5 PetscCheck(!*oinfo,PETSC_COMM_SELF,PETSC_ERR_LIB,"dsbtdc: Error in DIBTDC, oinfo = %" PetscBLASInt_FMT,*oinfo);
611
612 /* ************************************************************************** */
613
614 /* Scale back the computed eigenvalues. */
615
616
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.
 5 PetscCallLAPACKInfo("LAPACKlascl",LAPACKlascl_("G", &zero, &zero, &done,
617 &anorm, &nk, &one, &ev[np], &nk, &info));
618
619 5 start = iend + 1;
620 5 np += nk;
621
622 /* Go to the next irreducible subproblem. */
623
624 5 goto L10;
625 }
626
627 /* ....If the problem split any number of times, then the eigenvalues */
628 /* will not be properly ordered. Here we permute the eigenvalues */
629 /* (and the associated eigenvectors) across the irreducible parts */
630 /* into ascending order.... */
631
632 /* IF(NRBLKS.LT.NBLKS)THEN */
633
634 /* Use Selection Sort to minimize swaps of eigenvectors */
635
636
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
 120 for (ii = 1; ii < n; ++ii) {
637 115 i = ii;
638 115 k = i;
639 115 p = ev[i];
640
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
 1495 for (j = ii; j < n; ++j) {
641
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
 1380 if (ev[j] < p) {
642 ✗ k = j;
643 ✗ p = ev[j];
644 }
645 }
646
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
 115 if (k != i) {
647 ✗ ev[k] = ev[i];
648 ✗ ev[i] = p;
649
0/20
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ 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.
 115 PetscCallBLAS("BLASswap",BLASswap_(&n, &z[i*ldz], &one, &z[k*ldz], &one));
650 }
651 }
652
653 /* ...Compute MINGAP (minimum difference between neighboring eigenvalue */
654 /* approximations).............................................. */
655
656 5 *mingap = ev[1] - ev[0];
657
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 5 PetscCheck(*mingap>=0.,PETSC_COMM_SELF,PETSC_ERR_LIB,"dsbtdc: Eigenvalue approximations are not ordered properly. Approximation 1 is larger than approximation 2.");
658 5 *mingapi = 1;
659
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
 115 for (i = 2; i < n; ++i) {
660 110 absdiff = ev[i] - ev[i-1];
661
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 110 PetscCheck(absdiff>=0.,PETSC_COMM_SELF,PETSC_ERR_LIB,"dsbtdc: Eigenvalue approximations are not ordered properly. Approximation %" PetscBLASInt_FMT " is larger than approximation %" PetscBLASInt_FMT ".",i,i+1);
662
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
 110 if (absdiff < *mingap) {
663 5 *mingap = absdiff;
664 5 *mingapi = i;
665 }
666 }
667
668 /* check whether the minimum gap between eigenvalue approximations */
669 /* may indicate severe inaccuracies in the eigenvector approximations */
670
671
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
 5 if (*mingap <= tol / 10) *oinfo = -103;
672
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.
 1 PetscFunctionReturn(PETSC_SUCCESS);
673 }
674