GCC Code Coverage Report


Directory: ./
File: src/sys/classes/ds/impls/hep/dshep.c
Date: 2026-07-06 04:10:30
Exec Total Coverage
Lines: 517 545 94.9%
Functions: 20 20 100.0%
Branches: 1373 2461 55.8%

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 #include <slepc/private/dsimpl.h>
12 #include <slepcblaslapack.h>
13
14 3809 static PetscErrorCode DSAllocate_HEP(DS ds,PetscInt ld)
15 {
16
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
3809 PetscFunctionBegin;
17
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.
3809 if (!ds->compact) PetscCall(DSAllocateMat_Private(ds,DS_MAT_A));
18
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.
3809 PetscCall(DSAllocateMat_Private(ds,DS_MAT_Q));
19
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.
3809 PetscCall(DSAllocateMat_Private(ds,DS_MAT_T));
20
5/8
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 8 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
3809 PetscCall(PetscFree(ds->perm));
21
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.
3809 PetscCall(PetscMalloc1(ld,&ds->perm));
22
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.
717 PetscFunctionReturn(PETSC_SUCCESS);
23 }
24
25 /* 0 l k n-1
26 -----------------------------------------
27 |* . . |
28 | * . . |
29 | * . . |
30 | * . . |
31 |. . . . o o |
32 | o o |
33 | o o |
34 | o o |
35 | o o |
36 | o o |
37 |. . . . o o o o o o o x |
38 | x x x |
39 | x x x |
40 | x x x |
41 | x x x |
42 | x x x |
43 | x x x |
44 | x x x |
45 | x x x|
46 | x x|
47 -----------------------------------------
48 */
49
50 1431 static PetscErrorCode DSSwitchFormat_HEP(DS ds)
51 {
52 1431 PetscReal *T;
53 1431 PetscScalar *A;
54 1431 PetscInt i,n=ds->n,k=ds->k,ld=ds->ld;
55
56
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
1431 PetscFunctionBegin;
57 /* switch from compact (arrow) to dense storage */
58
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.
1431 PetscCall(MatDenseGetArrayWrite(ds->omat[DS_MAT_A],&A));
59
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.
1431 PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
60
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.
1431 PetscCall(PetscArrayzero(A,ld*ld));
61
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
1431 for (i=0;i<k;i++) {
62 A[i+i*ld] = T[i];
63 A[k+i*ld] = T[i+ld];
64 A[i+k*ld] = T[i+ld];
65 }
66 1431 A[k+k*ld] = T[k];
67
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
25223 for (i=k+1;i<n;i++) {
68 23792 A[i+i*ld] = T[i];
69 23792 A[i-1+i*ld] = T[i-1+ld];
70 23792 A[i+(i-1)*ld] = T[i-1+ld];
71 }
72
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
1431 if (ds->extrarow) A[n+(n-1)*ld] = T[n-1+ld];
73
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.
1431 PetscCall(MatDenseRestoreArrayWrite(ds->omat[DS_MAT_A],&A));
74
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.
1431 PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
75
6/12
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 2 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 2 times.
257 PetscFunctionReturn(PETSC_SUCCESS);
76 }
77
78 208 static PetscErrorCode DSView_HEP(DS ds,PetscViewer viewer)
79 {
80 208 PetscViewerFormat format;
81 208 PetscInt i,j,r,c,rows;
82 208 PetscReal *T,value;
83 208 const char *methodname[] = {
84 "Implicit QR method (_steqr)",
85 "Relatively Robust Representations (_stevr)",
86 "Divide and Conquer method (_stedc)",
87 "Block Divide and Conquer method (dsbtdc)"
88 };
89 208 const int nmeth=PETSC_STATIC_ARRAY_LENGTH(methodname);
90
91
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
208 PetscFunctionBegin;
92
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.
208 PetscCall(PetscViewerGetFormat(viewer,&format));
93
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
208 if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
94
6/8
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 10 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 1 times.
148 if (ds->bs>1) PetscCall(PetscViewerASCIIPrintf(viewer,"block size: %" PetscInt_FMT "\n",ds->bs));
95
5/8
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 8 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
148 if (ds->method<nmeth) PetscCall(PetscViewerASCIIPrintf(viewer,"solving the problem with: %s\n",methodname[ds->method]));
96
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.
148 PetscFunctionReturn(PETSC_SUCCESS);
97 }
98
1/2
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
60 if (ds->compact) {
99
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
60 PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
100
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
60 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE));
101
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
60 rows = ds->extrarow? ds->n+1: ds->n;
102
1/2
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
60 if (format == PETSC_VIEWER_ASCII_MATLAB) {
103
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
60 PetscCall(PetscViewerASCIIPrintf(viewer,"%% Size = %" PetscInt_FMT " %" PetscInt_FMT "\n",rows,ds->n));
104
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
60 PetscCall(PetscViewerASCIIPrintf(viewer,"zzz = zeros(%" PetscInt_FMT ",3);\n",3*ds->n));
105
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
60 PetscCall(PetscViewerASCIIPrintf(viewer,"zzz = [\n"));
106
7/8
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 10 times.
✓ Branch 3 taken 8 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
✓ Branch 6 taken 2 times.
✓ Branch 7 taken 2 times.
600 for (i=0;i<ds->n;i++) PetscCall(PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %18.16e\n",i+1,i+1,(double)T[i]));
107
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
570 for (i=0;i<rows-1;i++) {
108
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
510 r = PetscMax(i+2,ds->k+1);
109 510 c = i+1;
110
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.
510 PetscCall(PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %18.16e\n",r,c,(double)T[i+ds->ld]));
111
3/4
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
✓ Branch 2 taken 10 times.
✗ Branch 3 not taken.
510 if (i<ds->n-1 && ds->k<ds->n) { /* do not print vertical arrow when k=n */
112
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.
480 PetscCall(PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %18.16e\n",c,r,(double)T[i+ds->ld]));
113 }
114 }
115
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
60 PetscCall(PetscViewerASCIIPrintf(viewer,"];\n%s = spconvert(zzz);\n",DSMatName[DS_MAT_T]));
116 } else {
117 for (i=0;i<rows;i++) {
118 for (j=0;j<ds->n;j++) {
119 if (i==j) value = T[i];
120 else if ((i<ds->k && j==ds->k) || (i==ds->k && j<ds->k)) value = T[PetscMin(i,j)+ds->ld];
121 else if (i==j+1 && i>ds->k) value = T[i-1+ds->ld];
122 else if (i+1==j && j>ds->k) value = T[j-1+ds->ld];
123 else value = 0.0;
124 PetscCall(PetscViewerASCIIPrintf(viewer," %18.16e ",(double)value));
125 }
126 PetscCall(PetscViewerASCIIPrintf(viewer,"\n"));
127 }
128 }
129
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
60 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE));
130
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
60 PetscCall(PetscViewerFlush(viewer));
131
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
60 PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
132 } else PetscCall(DSViewMat(ds,viewer,DS_MAT_A));
133
1/8
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
60 if (ds->state>DS_STATE_INTERMEDIATE) PetscCall(DSViewMat(ds,viewer,DS_MAT_Q));
134
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.
12 PetscFunctionReturn(PETSC_SUCCESS);
135 }
136
137 245378 static PetscErrorCode DSVectors_HEP(DS ds,DSMatType mat,PetscInt *j,PetscReal *rnorm)
138 {
139 245378 PetscScalar *Z;
140 245378 const PetscScalar *Q;
141 245378 PetscInt ld = ds->ld;
142
143
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
245378 PetscFunctionBegin;
144
1/3
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
245378 switch (mat) {
145 245378 case DS_MAT_X:
146 case DS_MAT_Y:
147
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
245378 if (j) {
148
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.
244766 PetscCall(MatDenseGetArray(ds->omat[mat],&Z));
149
1/2
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
244766 if (ds->state>=DS_STATE_CONDENSED) {
150
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.
244766 PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_Q],&Q));
151
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.
244766 PetscCall(PetscArraycpy(Z+(*j)*ld,Q+(*j)*ld,ld));
152
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
244766 if (rnorm) *rnorm = PetscAbsScalar(Q[ds->n-1+(*j)*ld]);
153
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.
244766 PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_Q],&Q));
154 } else {
155 PetscCall(PetscArrayzero(Z+(*j)*ld,ld));
156 Z[(*j)+(*j)*ld] = 1.0;
157 if (rnorm) *rnorm = 0.0;
158 }
159
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
244766 PetscCall(MatDenseRestoreArray(ds->omat[mat],&Z));
160 } else {
161
5/8
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 8 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
612 if (ds->state>=DS_STATE_CONDENSED) PetscCall(MatCopy(ds->omat[DS_MAT_Q],ds->omat[mat],SAME_NONZERO_PATTERN));
162 else PetscCall(DSSetIdentity(ds,mat));
163 }
164 47169 break;
165 case DS_MAT_U:
166 case DS_MAT_V:
167 SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"Not implemented yet");
168 default:
169 SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Invalid mat parameter");
170 }
171
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.
47169 PetscFunctionReturn(PETSC_SUCCESS);
172 }
173
174 /*
175 ARROWTRIDIAG reduces a symmetric arrowhead matrix of the form
176
177 [ d 0 0 0 e ]
178 [ 0 d 0 0 e ]
179 A = [ 0 0 d 0 e ]
180 [ 0 0 0 d e ]
181 [ e e e e d ]
182
183 to tridiagonal form
184
185 [ d e 0 0 0 ]
186 [ e d e 0 0 ]
187 T = Q'*A*Q = [ 0 e d e 0 ],
188 [ 0 0 e d e ]
189 [ 0 0 0 e d ]
190
191 where Q is an orthogonal matrix. Rutishauser's algorithm is used to
192 perform the reduction, which requires O(n**2) flops. The accumulation
193 of the orthogonal factor Q, however, requires O(n**3) flops.
194
195 Arguments
196 =========
197
198 N (input) INTEGER
199 The order of the matrix A. N >= 0.
200
201 D (input/output) DOUBLE PRECISION array, dimension (N)
202 On entry, the diagonal entries of the matrix A to be
203 reduced.
204 On exit, the diagonal entries of the reduced matrix T.
205
206 E (input/output) DOUBLE PRECISION array, dimension (N-1)
207 On entry, the off-diagonal entries of the matrix A to be
208 reduced.
209 On exit, the subdiagonal entries of the reduced matrix T.
210
211 Q (input/output) DOUBLE PRECISION array, dimension (LDQ, N)
212 On exit, the orthogonal matrix Q.
213
214 LDQ (input) INTEGER
215 The leading dimension of the array Q.
216
217 Note
218 ====
219 Based on Fortran code contributed by Daniel Kressner
220 */
221 26256 PetscErrorCode DSArrowTridiag(PetscBLASInt n,PetscReal *d,PetscReal *e,PetscScalar *Q,PetscBLASInt ld)
222 {
223 26256 PetscBLASInt i,j,j2,one=1;
224 26256 PetscReal c,s,p,off,temp;
225
226
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
26256 PetscFunctionBegin;
227
8/14
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 2 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 2 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✓ Branch 13 taken 2 times.
26256 if (n<=2) PetscFunctionReturn(PETSC_SUCCESS);
228
229
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
206865 for (j=0;j<n-2;j++) {
230
231 /* Eliminate entry e(j) by a rotation in the planes (j,j+1) */
232 180732 temp = e[j+1];
233
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.
180732 PetscCallBLAS("LAPACKlartg",LAPACKREALlartg_(&temp,&e[j],&c,&s,&e[j+1]));
234 180732 s = -s;
235
236 /* Apply rotation to diagonal elements */
237 180732 temp = d[j+1];
238 180732 e[j] = c*s*(temp-d[j]);
239 180732 d[j+1] = s*s*d[j] + c*c*temp;
240 180732 d[j] = c*c*d[j] + s*s*temp;
241
242 /* Apply rotation to Q */
243 180732 j2 = j+2;
244
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.
180732 PetscCallBLAS("BLASrot",BLASMIXEDrot_(&j2,Q+j*ld,&one,Q+(j+1)*ld,&one,&c,&s));
245
246 /* Chase newly introduced off-diagonal entry to the top left corner */
247
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
1111586 for (i=j-1;i>=0;i--) {
248 930854 off = -s*e[i];
249 930854 e[i] = c*e[i];
250 930854 temp = e[i+1];
251
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.
930854 PetscCallBLAS("LAPACKlartg",LAPACKREALlartg_(&temp,&off,&c,&s,&e[i+1]));
252 930854 s = -s;
253 930854 temp = (d[i]-d[i+1])*s - 2.0*c*e[i];
254 930854 p = s*temp;
255 930854 d[i+1] += p;
256 930854 d[i] -= p;
257 930854 e[i] = -e[i] - c*temp;
258
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.
930854 PetscCallBLAS("BLASrot",BLASMIXEDrot_(&j2,Q+i*ld,&one,Q+(i+1)*ld,&one,&c,&s));
259 }
260 }
261
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.
5158 PetscFunctionReturn(PETSC_SUCCESS);
262 }
263
264 /*
265 Reduce to tridiagonal form by means of DSArrowTridiag.
266 */
267 107201 static PetscErrorCode DSIntermediate_HEP(DS ds)
268 {
269 107201 PetscInt i;
270 107201 PetscBLASInt n1 = 0,n2,lwork,l = 0,n = 0,ld,off;
271 107201 PetscScalar *Q,*work,*tau;
272 107201 const PetscScalar *A;
273 107201 PetscReal *d,*e;
274 107201 Mat At,Qt; /* trailing submatrices */
275
276
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
107201 PetscFunctionBegin;
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.
107201 PetscCall(PetscBLASIntCast(ds->n,&n));
278
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.
107201 PetscCall(PetscBLASIntCast(ds->l,&l));
279
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.
107201 PetscCall(PetscBLASIntCast(ds->ld,&ld));
280
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.
107201 PetscCall(PetscBLASIntCast(PetscMax(0,ds->k-l+1),&n1)); /* size of leading block, excl. locked */
281 107201 n2 = n-l; /* n2 = n1 + size of trailing block */
282 107201 off = l+l*ld;
283
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.
107201 PetscCall(DSGetArrayReal(ds,DS_MAT_T,&d));
284 107201 e = d+ld;
285
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.
107201 PetscCall(DSSetIdentity(ds,DS_MAT_Q));
286
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.
107201 PetscCall(MatDenseGetArray(ds->omat[DS_MAT_Q],&Q));
287
288
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
107201 if (ds->compact) {
289
290
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.
36204 if (ds->state<DS_STATE_INTERMEDIATE) PetscCall(DSArrowTridiag(n1,d+l,e+l,Q+off,ld));
291
292 } else {
293
294
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.
70997 PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_A],&A));
295
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
73459 for (i=0;i<l;i++) { d[i] = PetscRealPart(A[i+i*ld]); e[i] = 0.0; }
296
297
1/2
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
70997 if (ds->state<DS_STATE_INTERMEDIATE) {
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.
70997 PetscCall(MatDenseGetSubMatrix(ds->omat[DS_MAT_A],ds->l,ds->n,ds->l,ds->n,&At));
299
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.
70997 PetscCall(MatDenseGetSubMatrix(ds->omat[DS_MAT_Q],ds->l,ds->n,ds->l,ds->n,&Qt));
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.
70997 PetscCall(MatCopy(At,Qt,SAME_NONZERO_PATTERN));
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.
70997 PetscCall(MatDenseRestoreSubMatrix(ds->omat[DS_MAT_A],&At));
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.
70997 PetscCall(MatDenseRestoreSubMatrix(ds->omat[DS_MAT_Q],&Qt));
303
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.
70997 PetscCall(DSAllocateWork_Private(ds,ld+ld*ld,0,0));
304 70997 tau = ds->work;
305 70997 work = ds->work+ld;
306 70997 lwork = ld*ld;
307
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.
70997 PetscCallLAPACKInfo("LAPACKsytrd",LAPACKsytrd_("L",&n2,Q+off,&ld,d+l,e+l,tau,work,&lwork,&info));
308
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.
70997 PetscCallLAPACKInfo("LAPACKorgtr",LAPACKorgtr_("L",&n2,Q+off,&ld,tau,work,&lwork,&info));
309 } else {
310 /* copy tridiagonal to d,e */
311 for (i=l;i<n;i++) d[i] = PetscRealPart(A[i+i*ld]);
312 for (i=l;i<n-1;i++) e[i] = PetscRealPart(A[(i+1)+i*ld]);
313 }
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.
70997 PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_A],&A));
315 }
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.
107201 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.
107201 PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&d));
318
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.
20703 PetscFunctionReturn(PETSC_SUCCESS);
319 }
320
321 174037 static PetscErrorCode DSSort_HEP(DS ds,PetscScalar *wr,PetscScalar *wi,PetscScalar *rr,PetscScalar *ri,PetscInt *k)
322 {
323 174037 PetscInt n,l,i,*perm,ld=ds->ld;
324 174037 PetscScalar *A;
325 174037 PetscReal *d;
326
327
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
174037 PetscFunctionBegin;
328
2/14
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
174037 if (!ds->sc) PetscFunctionReturn(PETSC_SUCCESS);
329 174037 n = ds->n;
330 174037 l = ds->l;
331
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.
174037 PetscCall(DSGetArrayReal(ds,DS_MAT_T,&d));
332 174037 perm = ds->perm;
333
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.
174037 if (!rr) PetscCall(DSSortEigenvaluesReal_Private(ds,d,perm));
334
4/6
✓ Branch 0 taken 9 times.
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
135340 else PetscCall(DSSortEigenvalues_Private(ds,rr,ri,perm,PETSC_FALSE));
335
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
2306716 for (i=l;i<n;i++) wr[i] = d[perm[i]];
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.
174037 PetscCall(DSPermuteColumns_Private(ds,l,n,n,DS_MAT_Q,perm));
337
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
2306716 for (i=l;i<n;i++) d[i] = PetscRealPart(wr[i]);
338
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
174037 if (!ds->compact) {
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.
137833 PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
340
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
1712379 for (i=l;i<n;i++) A[i+i*ld] = wr[i];
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.
137833 PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
342 }
343
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.
174037 PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&d));
344
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.
33599 PetscFunctionReturn(PETSC_SUCCESS);
345 }
346
347 26079 static PetscErrorCode DSUpdateExtraRow_HEP(DS ds)
348 {
349 26079 PetscInt i;
350 26079 PetscBLASInt n,ld,incx=1;
351 26079 PetscScalar *A,*x,*y,one=1.0,zero=0.0;
352 26079 PetscReal *T,*e,beta;
353 26079 const PetscScalar *Q;
354
355
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
26079 PetscFunctionBegin;
356
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.
26079 PetscCall(PetscBLASIntCast(ds->n,&n));
357
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.
26079 PetscCall(PetscBLASIntCast(ds->ld,&ld));
358
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.
26079 PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_Q],&Q));
359
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
26079 if (ds->compact) {
360
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.
26049 PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
361 26049 e = T+ld;
362 26049 beta = e[n-1]; /* in compact, we assume all entries are zero except the last one */
363
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
462769 for (i=0;i<n;i++) e[i] = PetscRealPart(beta*Q[n-1+i*ld]);
364
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.
26049 PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
365 26049 ds->k = n;
366 } else {
367
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_A],&A));
368
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(DSAllocateWork_Private(ds,2*ld,0,0));
369 30 x = ds->work;
370 30 y = ds->work+ld;
371
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
390 for (i=0;i<n;i++) x[i] = PetscConj(A[n+i*ld]);
372
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.
30 PetscCallBLAS("BLASgemv",BLASgemv_("C",&n,&n,&one,Q,&ld,x,&incx,&zero,y,&incx));
373
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
390 for (i=0;i<n;i++) A[n+i*ld] = PetscConj(y[i]);
374 30 ds->k = n;
375
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_A],&A));
376 }
377
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.
26079 PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_Q],&Q));
378
6/12
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 2 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 2 times.
5029 PetscFunctionReturn(PETSC_SUCCESS);
379 }
380
381 107121 static PetscErrorCode DSSolve_HEP_QR(DS ds,PetscScalar *wr,PetscScalar *wi)
382 {
383 107121 PetscInt i;
384 107121 PetscBLASInt n1,l = 0,n = 0,ld,off;
385 107121 PetscScalar *Q,*A;
386 107121 PetscReal *d,*e;
387
388
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
107121 PetscFunctionBegin;
389
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
107121 PetscCheck(ds->bs==1,PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"This method is not prepared for bs>1");
390
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.
107121 PetscCall(PetscBLASIntCast(ds->n,&n));
391
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.
107121 PetscCall(PetscBLASIntCast(ds->l,&l));
392
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.
107121 PetscCall(PetscBLASIntCast(ds->ld,&ld));
393 107121 n1 = n-l; /* n1 = size of leading block, excl. locked + size of trailing block */
394 107121 off = l+l*ld;
395
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.
107121 PetscCall(DSGetArrayReal(ds,DS_MAT_T,&d));
396 107121 e = d+ld;
397
398 /* Reduce to tridiagonal form */
399
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.
107121 PetscCall(DSIntermediate_HEP(ds));
400
401 /* Solve the tridiagonal eigenproblem */
402
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
185862 for (i=0;i<l;i++) wr[i] = d[i];
403
404
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.
107121 PetscCall(DSAllocateWork_Private(ds,0,2*ld,0));
405
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.
107121 PetscCall(MatDenseGetArray(ds->omat[DS_MAT_Q],&Q));
406
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.
107121 PetscCallLAPACKInfo("LAPACKsteqr",LAPACKsteqr_("V",&n1,d+l,e+l,Q+off,&ld,ds->rwork,&info));
407
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.
107121 PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_Q],&Q));
408
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
1491714 for (i=l;i<n;i++) wr[i] = d[i];
409
410 /* Create diagonal matrix as a result */
411
6/8
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
✓ Branch 2 taken 5 times.
✓ Branch 3 taken 5 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
107121 if (ds->compact) PetscCall(PetscArrayzero(e,n-1));
412 else {
413
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.
70957 PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
414
7/8
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
✓ Branch 2 taken 10 times.
✓ Branch 3 taken 8 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
✓ Branch 6 taken 2 times.
✓ Branch 7 taken 2 times.
897697 for (i=l;i<n;i++) PetscCall(PetscArrayzero(A+l+i*ld,n-l));
415
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
897697 for (i=l;i<n;i++) A[i+i*ld] = d[i];
416
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.
70957 PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
417 }
418
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.
107121 PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&d));
419
420 /* Set zero wi */
421
4/4
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
✓ Branch 2 taken 10 times.
✓ Branch 3 taken 10 times.
1318459 if (wi) for (i=l;i<n;i++) wi[i] = 0.0;
422
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.
20687 PetscFunctionReturn(PETSC_SUCCESS);
423 }
424
425 40 static PetscErrorCode DSSolve_HEP_MRRR(DS ds,PetscScalar *wr,PetscScalar *wi)
426 {
427 40 Mat At,Qt; /* trailing submatrices */
428 40 PetscInt i;
429 40 PetscBLASInt n1 = 0,n2 = 0,n3,lrwork,liwork,l = 0,n = 0,m = 0,ld,off,il,iu,*isuppz;
430 40 PetscScalar *A,*Q,*W=NULL,one=1.0,zero=0.0;
431 40 PetscReal *d,*e,abstol=0.0,vl,vu;
432 #if defined(PETSC_USE_COMPLEX)
433 20 PetscInt j;
434 20 PetscReal *Qr,*ritz;
435 #endif
436
437
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
40 PetscFunctionBegin;
438
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
40 PetscCheck(ds->bs==1,PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"This method is not prepared for bs>1");
439
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
40 PetscCall(PetscBLASIntCast(ds->n,&n));
440
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
40 PetscCall(PetscBLASIntCast(ds->l,&l));
441
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));
442
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->k-l+1,&n1)); /* size of leading block, excl. locked */
443
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(n-ds->k-1,&n2)); /* size of trailing block */
444 40 n3 = n1+n2;
445 40 off = l+l*ld;
446
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));
447 40 e = d+ld;
448
449 /* Reduce to tridiagonal form */
450
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(DSIntermediate_HEP(ds));
451
452 /* Solve the tridiagonal eigenproblem */
453
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
80 for (i=0;i<l;i++) wr[i] = d[i];
454
455
1/2
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
40 if (ds->state<DS_STATE_INTERMEDIATE) { /* Q contains useful info */
456
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(DSAllocateMat_Private(ds,DS_MAT_W));
457
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(MatCopy(ds->omat[DS_MAT_Q],ds->omat[DS_MAT_W],SAME_NONZERO_PATTERN));
458 }
459
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(MatDenseGetArray(ds->omat[DS_MAT_Q],&Q));
460 40 lrwork = 20*ld;
461 40 liwork = 10*ld;
462 #if defined(PETSC_USE_COMPLEX)
463
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
20 PetscCall(DSAllocateWork_Private(ds,0,lrwork+ld+ld*ld,liwork+2*ld));
464 #else
465
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
20 PetscCall(DSAllocateWork_Private(ds,0,lrwork+ld,liwork+2*ld));
466 #endif
467 40 isuppz = ds->iwork+liwork;
468 #if defined(PETSC_USE_COMPLEX)
469 20 ritz = ds->rwork+lrwork;
470 20 Qr = ds->rwork+lrwork+ld;
471
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.
20 PetscCallLAPACKInfo("LAPACKstevr",LAPACKstevr_("V","A",&n3,d+l,e+l,&vl,&vu,&il,&iu,&abstol,&m,ritz+l,Qr+off,&ld,isuppz,ds->rwork,&lrwork,ds->iwork,&liwork,&info));
472
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
210 for (i=l;i<n;i++) wr[i] = ritz[i];
473 #else
474
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.
20 PetscCallLAPACKInfo("LAPACKstevr",LAPACKstevr_("V","A",&n3,d+l,e+l,&vl,&vu,&il,&iu,&abstol,&m,wr+l,Q+off,&ld,isuppz,ds->rwork,&lrwork,ds->iwork,&liwork,&info));
475 #endif
476 #if defined(PETSC_USE_COMPLEX)
477
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
210 for (i=l;i<n;i++)
478
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
2120 for (j=l;j<n;j++)
479 1930 Q[i+j*ld] = Qr[i+j*ld];
480 #endif
481
1/2
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
40 if (ds->state<DS_STATE_INTERMEDIATE) { /* accumulate previous Q */
482
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.
40 if (ds->compact) PetscCall(DSAllocateMat_Private(ds,DS_MAT_A));
483
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(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
484
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(MatDenseGetArray(ds->omat[DS_MAT_W],&W));
485
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.
40 PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n3,&n3,&n3,&one,W+off,&ld,Q+off,&ld,&zero,A+off,&ld));
486
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(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
487
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(MatDenseRestoreArray(ds->omat[DS_MAT_W],&W));
488
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(MatDenseGetSubMatrix(ds->omat[DS_MAT_A],ds->l,ds->n,ds->l,ds->n,&At));
489
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
40 PetscCall(MatDenseGetSubMatrix(ds->omat[DS_MAT_Q],ds->l,ds->n,ds->l,ds->n,&Qt));
490
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
40 PetscCall(MatCopy(At,Qt,SAME_NONZERO_PATTERN));
491
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
40 PetscCall(MatDenseRestoreSubMatrix(ds->omat[DS_MAT_A],&At));
492
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
40 PetscCall(MatDenseRestoreSubMatrix(ds->omat[DS_MAT_Q],&Qt));
493 }
494
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(MatDenseRestoreArray(ds->omat[DS_MAT_Q],&Q));
495
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
420 for (i=l;i<n;i++) d[i] = PetscRealPart(wr[i]);
496
497 /* Create diagonal matrix as a result */
498
6/8
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
✓ Branch 2 taken 5 times.
✓ Branch 3 taken 5 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
40 if (ds->compact) PetscCall(PetscArrayzero(e,n-1));
499 else {
500
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
20 PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
501
7/8
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
✓ Branch 2 taken 10 times.
✓ Branch 3 taken 8 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
✓ Branch 6 taken 2 times.
✓ Branch 7 taken 2 times.
260 for (i=l;i<n;i++) PetscCall(PetscArrayzero(A+l+i*ld,n-l));
502
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
260 for (i=l;i<n;i++) A[i+i*ld] = d[i];
503
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
20 PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
504 }
505
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(DSRestoreArrayReal(ds,DS_MAT_T,&d));
506
507 /* Set zero wi */
508
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
40 if (wi) for (i=l;i<n;i++) wi[i] = 0.0;
509
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.
8 PetscFunctionReturn(PETSC_SUCCESS);
510 }
511
512 40 static PetscErrorCode DSSolve_HEP_DC(DS ds,PetscScalar *wr,PetscScalar *wi)
513 {
514 40 PetscInt i;
515 40 PetscBLASInt n1,l = 0,ld,off,lrwork,liwork;
516 40 PetscScalar *Q,*A;
517 40 PetscReal *d,*e;
518 #if defined(PETSC_USE_COMPLEX)
519 20 PetscBLASInt lwork;
520 20 PetscInt j;
521 #endif
522
523
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
40 PetscFunctionBegin;
524
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
40 PetscCheck(ds->bs==1,PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"This method is not prepared for bs>1");
525
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->l,&l));
526
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));
527
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));
528 40 off = l+l*ld;
529
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));
530 40 e = d+ld;
531
532 /* Reduce to tridiagonal form */
533
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(DSIntermediate_HEP(ds));
534
535 /* Solve the tridiagonal eigenproblem */
536
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
80 for (i=0;i<l;i++) wr[i] = d[i];
537
538 40 lrwork = 5*n1*n1+3*n1+1;
539 40 liwork = 5*n1*n1+6*n1+6;
540
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(MatDenseGetArray(ds->omat[DS_MAT_Q],&Q));
541 #if !defined(PETSC_USE_COMPLEX)
542
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
20 PetscCall(DSAllocateWork_Private(ds,0,lrwork,liwork));
543
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.
20 PetscCallLAPACKInfo("LAPACKstedc",LAPACKstedc_("V",&n1,d+l,e+l,Q+off,&ld,ds->rwork,&lrwork,ds->iwork,&liwork,&info));
544 #else
545 20 lwork = ld*ld;
546
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
20 PetscCall(DSAllocateWork_Private(ds,lwork,lrwork,liwork));
547
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.
20 PetscCallLAPACKInfo("LAPACKstedc",LAPACKstedc_("V",&n1,d+l,e+l,Q+off,&ld,ds->work,&lwork,ds->rwork,&lrwork,ds->iwork,&liwork,&info));
548 /* Fixing Lapack bug*/
549
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
210 for (j=ds->l;j<ds->n;j++)
550
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
330 for (i=0;i<ds->l;i++) Q[i+j*ld] = 0.0;
551 #endif
552
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(MatDenseRestoreArray(ds->omat[DS_MAT_Q],&Q));
553
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
420 for (i=l;i<ds->n;i++) wr[i] = d[i];
554
555 /* Create diagonal matrix as a result */
556
6/8
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
✓ Branch 2 taken 5 times.
✓ Branch 3 taken 5 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
40 if (ds->compact) PetscCall(PetscArrayzero(e,ds->n-1));
557 else {
558
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
20 PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
559
7/8
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
✓ Branch 2 taken 10 times.
✓ Branch 3 taken 8 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
✓ Branch 6 taken 2 times.
✓ Branch 7 taken 2 times.
260 for (i=l;i<ds->n;i++) PetscCall(PetscArrayzero(A+l+i*ld,ds->n-l));
560
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
260 for (i=l;i<ds->n;i++) A[i+i*ld] = d[i];
561
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
20 PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
562 }
563
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(DSRestoreArrayReal(ds,DS_MAT_T,&d));
564
565 /* Set zero wi */
566
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
40 if (wi) for (i=l;i<ds->n;i++) wi[i] = 0.0;
567
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.
8 PetscFunctionReturn(PETSC_SUCCESS);
568 }
569
570 #if !defined(PETSC_USE_COMPLEX)
571 10 static PetscErrorCode DSSolve_HEP_BDC(DS ds,PetscScalar *wr,PetscScalar *wi)
572 {
573 10 PetscBLASInt i,j,k,m,n = 0,info,nblks,bs = 0,ld = 0,lde,lrwork,liwork,*ksizes,*iwork,mingapi;
574 10 PetscScalar *Q,*A;
575 10 PetscReal *D,*E,*d,*e,tol=PETSC_MACHINE_EPSILON/2,tau1=1e-16,tau2=1e-18,*rwork,mingap;
576
577
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
10 PetscFunctionBegin;
578
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
10 PetscCheck(ds->l==0,PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"This method is not prepared for l>1");
579
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
10 PetscCheck(!ds->compact,PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"Not implemented for compact storage");
580
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.
10 PetscCall(PetscBLASIntCast(ds->ld,&ld));
581
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.
10 PetscCall(PetscBLASIntCast(ds->bs,&bs));
582
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.
10 PetscCall(PetscBLASIntCast(ds->n,&n));
583 10 nblks = n/bs;
584
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.
10 PetscCall(DSGetArrayReal(ds,DS_MAT_T,&d));
585 10 e = d+ld;
586 10 lrwork = 4*n*n+60*n+1;
587 10 liwork = 5*n+5*nblks-1;
588 10 lde = 2*bs+1;
589
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.
10 PetscCall(DSAllocateWork_Private(ds,bs*n+lde*lde*(nblks-1),lrwork,nblks+liwork));
590 10 D = ds->work;
591 10 E = ds->work+bs*n;
592 10 rwork = ds->rwork;
593 10 ksizes = ds->iwork;
594 10 iwork = ds->iwork+nblks;
595
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 3 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
10 PetscCall(PetscArrayzero(iwork,liwork));
596
597 /* Copy matrix to block tridiagonal format */
598
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.
10 PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
599 j=0;
600
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
55 for (i=0;i<nblks;i++) {
601 45 ksizes[i]=bs;
602
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
180 for (k=0;k<bs;k++)
603
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
540 for (m=0;m<bs;m++)
604 405 D[k+m*bs+i*bs*bs] = PetscRealPart(A[j+k+(j+m)*n]);
605 45 j = j + bs;
606 }
607 j=0;
608
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
45 for (i=0;i<nblks-1;i++) {
609
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
140 for (k=0;k<bs;k++)
610
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
420 for (m=0;m<bs;m++)
611 315 E[k+m*lde+i*lde*lde] = PetscRealPart(A[j+bs+k+(j+m)*n]);
612 35 j = j + bs;
613 }
614
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.
10 PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
615
616 /* Solve the block tridiagonal eigenproblem */
617
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.
10 PetscCall(MatDenseGetArray(ds->omat[DS_MAT_Q],&Q));
618
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.
10 PetscCall(BDC_dsbtdc_("D","A",n,nblks,ksizes,D,bs,bs,E,lde,lde,tol,tau1,tau2,d,Q,n,rwork,lrwork,iwork,liwork,&mingap,&mingapi,&info,1,1));
619
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.
10 PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_Q],&Q));
620
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
145 for (i=0;i<ds->n;i++) wr[i] = d[i];
621
622 /* Create diagonal matrix as a result */
623
1/8
✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
10 if (ds->compact) PetscCall(PetscArrayzero(e,ds->n-1));
624 else {
625
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.
10 PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
626
7/8
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 3 times.
✓ Branch 2 taken 5 times.
✓ Branch 3 taken 4 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
✓ Branch 6 taken 1 times.
✓ Branch 7 taken 1 times.
145 for (i=0;i<ds->n;i++) PetscCall(PetscArrayzero(A+i*ld,ds->n));
627
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
145 for (i=0;i<ds->n;i++) A[i+i*ld] = wr[i];
628
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.
10 PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
629 }
630
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.
10 PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&d));
631
632 /* Set zero wi */
633
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
10 if (wi) for (i=0;i<ds->n;i++) wi[i] = 0.0;
634
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.
2 PetscFunctionReturn(PETSC_SUCCESS);
635 }
636 #endif
637
638 26180 static PetscErrorCode DSTruncate_HEP(DS ds,PetscInt n,PetscBool trim)
639 {
640 26180 PetscInt i,ld=ds->ld,l=ds->l;
641 26180 PetscScalar *A;
642
643
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
26180 PetscFunctionBegin;
644
3/10
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 10 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
26180 if (!ds->compact && ds->extrarow) PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
645
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
26180 if (trim) {
646
3/4
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 10 times.
2926 if (!ds->compact && ds->extrarow) { /* clean extra row */
647 for (i=l;i<ds->n;i++) A[ds->n+i*ld] = 0.0;
648 }
649 2926 ds->l = 0;
650 2926 ds->k = 0;
651 2926 ds->n = n;
652 2926 ds->t = ds->n; /* truncated length equal to the new dimension */
653 } else {
654
1/6
✗ 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.
23254 if (!ds->compact && ds->extrarow && ds->k==ds->n) {
655 /* copy entries of extra row to the new position, then clean last row */
656 for (i=l;i<n;i++) A[n+i*ld] = A[ds->n+i*ld];
657 for (i=l;i<ds->n;i++) A[ds->n+i*ld] = 0.0;
658 }
659
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
23254 ds->k = (ds->extrarow)? n: 0;
660 23254 ds->t = ds->n; /* truncated length equal to previous dimension */
661 23254 ds->n = n;
662 }
663
3/10
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 10 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
26180 if (!ds->compact && ds->extrarow) PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
664
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.
5046 PetscFunctionReturn(PETSC_SUCCESS);
665 }
666
667 #if !defined(PETSC_HAVE_MPIUNI)
668 20 static PetscErrorCode DSSynchronize_HEP(DS ds,PetscScalar eigr[],PetscScalar eigi[])
669 {
670 20 PetscInt ld=ds->ld,l=ds->l,k=0,kr=0;
671 20 PetscMPIInt n,rank,off=0,size,ldn,ld3;
672 20 PetscScalar *A,*Q;
673 20 PetscReal *T;
674
675
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
20 PetscFunctionBegin;
676
1/2
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
20 if (ds->compact) kr = 3*ld;
677 else k = (ds->n-l)*ld;
678
1/2
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
20 if (ds->state>DS_STATE_RAW) k += (ds->n-l)*ld;
679
1/2
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
20 if (eigr) k += (ds->n-l);
680
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
20 PetscCall(DSAllocateWork_Private(ds,k+kr,0,0));
681
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
20 PetscCall(PetscMPIIntCast(k*sizeof(PetscScalar)+kr*sizeof(PetscReal),&size));
682
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
20 PetscCall(PetscMPIIntCast(ds->n-l,&n));
683
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
20 PetscCall(PetscMPIIntCast(ld*(ds->n-l),&ldn));
684
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
20 PetscCall(PetscMPIIntCast(ld*3,&ld3));
685
5/8
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 8 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
20 if (ds->compact) PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
686 else PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
687
5/8
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 8 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
20 if (ds->state>DS_STATE_RAW) PetscCall(MatDenseGetArray(ds->omat[DS_MAT_Q],&Q));
688
14/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 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.
✓ Branch 12 taken 2 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✓ Branch 15 taken 2 times.
✓ Branch 16 taken 2 times.
✗ Branch 17 not taken.
✓ Branch 18 taken 2 times.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✓ Branch 21 taken 2 times.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✓ Branch 25 taken 2 times.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
20 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)ds),&rank));
689
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
20 if (!rank) {
690
15/30
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 8 times.
✓ Branch 4 taken 2 times.
✓ Branch 5 taken 8 times.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 2 times.
✗ Branch 9 not taken.
✓ 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 not taken.
✓ Branch 17 taken 2 times.
✓ Branch 18 taken 2 times.
✗ Branch 19 not taken.
✓ Branch 20 taken 2 times.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✓ Branch 23 taken 2 times.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✓ Branch 27 taken 2 times.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
10 if (ds->compact) PetscCallMPI(MPI_Pack(T,ld3,MPIU_REAL,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
691 else PetscCallMPI(MPI_Pack(A+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
692
15/30
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 8 times.
✓ Branch 4 taken 2 times.
✓ Branch 5 taken 8 times.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 2 times.
✗ Branch 9 not taken.
✓ 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 not taken.
✓ Branch 17 taken 2 times.
✓ Branch 18 taken 2 times.
✗ Branch 19 not taken.
✓ Branch 20 taken 2 times.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✓ Branch 23 taken 2 times.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✓ Branch 27 taken 2 times.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
10 if (ds->state>DS_STATE_RAW) PetscCallMPI(MPI_Pack(Q+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
693
15/30
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 8 times.
✓ Branch 4 taken 2 times.
✓ Branch 5 taken 8 times.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 2 times.
✗ Branch 9 not taken.
✓ 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 not taken.
✓ Branch 17 taken 2 times.
✓ Branch 18 taken 2 times.
✗ Branch 19 not taken.
✓ Branch 20 taken 2 times.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✓ Branch 23 taken 2 times.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✓ Branch 27 taken 2 times.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
10 if (eigr) PetscCallMPI(MPI_Pack(eigr+l,n,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
694 }
695
15/30
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 8 times.
✓ Branch 4 taken 2 times.
✓ Branch 5 taken 8 times.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 2 times.
✗ Branch 9 not taken.
✓ 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 not taken.
✓ Branch 17 taken 2 times.
✓ Branch 18 taken 2 times.
✗ Branch 19 not taken.
✓ Branch 20 taken 2 times.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✓ Branch 23 taken 2 times.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✓ Branch 27 taken 2 times.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
40 PetscCallMPI(MPI_Bcast(ds->work,size,MPI_BYTE,0,PetscObjectComm((PetscObject)ds)));
696
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
20 if (rank) {
697
15/30
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 8 times.
✓ Branch 4 taken 2 times.
✓ Branch 5 taken 8 times.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 2 times.
✗ Branch 9 not taken.
✓ 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 not taken.
✓ Branch 17 taken 2 times.
✓ Branch 18 taken 2 times.
✗ Branch 19 not taken.
✓ Branch 20 taken 2 times.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✓ Branch 23 taken 2 times.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✓ Branch 27 taken 2 times.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
10 if (ds->compact) PetscCallMPI(MPI_Unpack(ds->work,size,&off,T,ld3,MPIU_REAL,PetscObjectComm((PetscObject)ds)));
698 else PetscCallMPI(MPI_Unpack(ds->work,size,&off,A+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
699
15/30
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 8 times.
✓ Branch 4 taken 2 times.
✓ Branch 5 taken 8 times.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 2 times.
✗ Branch 9 not taken.
✓ 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 not taken.
✓ Branch 17 taken 2 times.
✓ Branch 18 taken 2 times.
✗ Branch 19 not taken.
✓ Branch 20 taken 2 times.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✓ Branch 23 taken 2 times.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✓ Branch 27 taken 2 times.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
10 if (ds->state>DS_STATE_RAW) PetscCallMPI(MPI_Unpack(ds->work,size,&off,Q+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
700
15/30
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 8 times.
✓ Branch 4 taken 2 times.
✓ Branch 5 taken 8 times.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 2 times.
✗ Branch 9 not taken.
✓ 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 not taken.
✓ Branch 17 taken 2 times.
✓ Branch 18 taken 2 times.
✗ Branch 19 not taken.
✓ Branch 20 taken 2 times.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✓ Branch 23 taken 2 times.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✓ Branch 27 taken 2 times.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
10 if (eigr) PetscCallMPI(MPI_Unpack(ds->work,size,&off,eigr+l,n,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
701 }
702
5/8
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 8 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
20 if (ds->compact) PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
703 else PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
704
5/8
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 8 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
20 if (ds->state>DS_STATE_RAW) PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_Q],&Q));
705
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.
4 PetscFunctionReturn(PETSC_SUCCESS);
706 }
707 #endif
708
709 1431 static PetscErrorCode DSCond_HEP(DS ds,PetscReal *cond)
710 {
711 1431 PetscScalar *work;
712 1431 PetscReal *rwork;
713 1431 PetscBLASInt *ipiv;
714 1431 PetscBLASInt lwork,n,ld;
715 1431 PetscReal hn,hin;
716 1431 PetscScalar *A;
717
718
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
1431 PetscFunctionBegin;
719
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.
1431 PetscCall(PetscBLASIntCast(ds->n,&n));
720
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.
1431 PetscCall(PetscBLASIntCast(ds->ld,&ld));
721 1431 lwork = 8*ld;
722
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
1431 PetscCall(DSAllocateWork_Private(ds,lwork,ld,ld));
723 1431 work = ds->work;
724 1431 rwork = ds->rwork;
725 1431 ipiv = ds->iwork;
726
1/8
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
1431 if (ds->compact) PetscCall(DSAllocateMat_Private(ds,DS_MAT_A));
727
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.
1431 PetscCall(DSSwitchFormat_HEP(ds));
728
729 /* use workspace matrix W to avoid overwriting A */
730
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.
1431 PetscCall(DSAllocateMat_Private(ds,DS_MAT_W));
731
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.
1431 PetscCall(MatCopy(ds->omat[DS_MAT_A],ds->omat[DS_MAT_W],SAME_NONZERO_PATTERN));
732
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.
1431 PetscCall(MatDenseGetArray(ds->omat[DS_MAT_W],&A));
733
734 /* norm of A */
735 1431 hn = LAPACKlange_("I",&n,&n,A,&ld,rwork);
736
737 /* norm of inv(A) */
738
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.
1431 PetscCallLAPACKInfo("LAPACKgetrf",LAPACKgetrf_(&n,&n,A,&ld,ipiv,&info));
739
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.
1431 PetscCallLAPACKInfo("LAPACKgetri",LAPACKgetri_(&n,A,&ld,ipiv,work,&lwork,&info));
740 1431 hin = LAPACKlange_("I",&n,&n,A,&ld,rwork);
741
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.
1431 PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_W],&A));
742
743 1431 *cond = hn*hin;
744
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.
1431 PetscFunctionReturn(PETSC_SUCCESS);
745 }
746
747 126 static PetscErrorCode DSTranslateRKS_HEP(DS ds,PetscScalar alpha)
748 {
749 126 PetscInt i,j,k=ds->k;
750 126 PetscScalar *Q,*A,*R,*tau,*work;
751 126 PetscBLASInt ld,n1,n0,lwork;
752
753
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
126 PetscFunctionBegin;
754
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.
126 PetscCall(PetscBLASIntCast(ds->ld,&ld));
755
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.
126 PetscCall(DSAllocateWork_Private(ds,ld*ld,0,0));
756 126 tau = ds->work;
757 126 work = ds->work+ld;
758
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.
126 PetscCall(PetscBLASIntCast(ld*(ld-1),&lwork));
759
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.
126 PetscCall(DSAllocateMat_Private(ds,DS_MAT_W));
760
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.
126 PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
761
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
126 PetscCall(MatDenseGetArrayWrite(ds->omat[DS_MAT_Q],&Q));
762
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.
126 PetscCall(MatDenseGetArrayWrite(ds->omat[DS_MAT_W],&R));
763
764 /* copy I+alpha*A */
765
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.
126 PetscCall(PetscArrayzero(Q,ld*ld));
766
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.
126 PetscCall(PetscArrayzero(R,ld*ld));
767
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
1863 for (i=0;i<k;i++) {
768 1737 Q[i+i*ld] = 1.0 + alpha*A[i+i*ld];
769 1737 Q[k+i*ld] = alpha*A[k+i*ld];
770 }
771
772 /* compute qr */
773
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.
126 PetscCall(PetscBLASIntCast(k+1,&n1));
774
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.
126 PetscCall(PetscBLASIntCast(k,&n0));
775
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.
126 PetscCallLAPACKInfo("LAPACKgeqrf",LAPACKgeqrf_(&n1,&n0,Q,&ld,tau,work,&lwork,&info));
776
777 /* copy R from Q */
778
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
1863 for (j=0;j<k;j++)
779
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
14872 for (i=0;i<=j;i++)
780 13135 R[i+j*ld] = Q[i+j*ld];
781
782 /* compute orthogonal matrix in Q */
783
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.
126 PetscCallLAPACKInfo("LAPACKorgqr",LAPACKorgqr_(&n1,&n1,&n0,Q,&ld,tau,work,&lwork,&info));
784
785 /* compute the updated matrix of projected problem */
786
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
1863 for (j=0;j<k;j++)
787
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
28007 for (i=0;i<k+1;i++)
788 26270 A[j*ld+i] = Q[i*ld+j];
789 126 alpha = -1.0/alpha;
790
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.
126 PetscCallBLAS("BLAStrsm",BLAStrsm_("R","U","N","N",&n1,&n0,&alpha,R,&ld,A,&ld));
791
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
1863 for (i=0;i<k;i++)
792 1737 A[ld*i+i] -= alpha;
793
794
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.
126 PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
795
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.
126 PetscCall(MatDenseRestoreArrayWrite(ds->omat[DS_MAT_Q],&Q));
796
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.
126 PetscCall(MatDenseRestoreArrayWrite(ds->omat[DS_MAT_W],&R));
797
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.
23 PetscFunctionReturn(PETSC_SUCCESS);
798 }
799
800 237667 static PetscErrorCode DSHermitian_HEP(DS ds,DSMatType m,PetscBool *flg)
801 {
802
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
237667 PetscFunctionBegin;
803
4/4
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
✓ Branch 2 taken 10 times.
✓ Branch 3 taken 10 times.
237667 if (m==DS_MAT_A && !ds->extrarow) *flg = PETSC_TRUE;
804 165832 else *flg = PETSC_FALSE;
805
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.
237667 PetscFunctionReturn(PETSC_SUCCESS);
806 }
807
808 688 static PetscErrorCode DSSetCompact_HEP(DS ds,PetscBool comp)
809 {
810
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
688 PetscFunctionBegin;
811
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.
688 if (!comp) PetscCall(DSAllocateMat_Private(ds,DS_MAT_A));
812
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.
130 PetscFunctionReturn(PETSC_SUCCESS);
813 }
814
815 63 static PetscErrorCode DSReallocate_HEP(DS ds,PetscInt ld)
816 {
817 63 PetscInt i,*perm=ds->perm;
818
819
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
63 PetscFunctionBegin;
820
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
1449 for (i=0;i<DS_NUM_MAT;i++) {
821
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
1386 if (!ds->compact && i==DS_MAT_A) continue;
822
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.
1386 if (i!=DS_MAT_Q && i!=DS_MAT_T) PetscCall(MatDestroy(&ds->omat[i]));
823 }
824
825
1/8
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
63 if (!ds->compact) PetscCall(DSReallocateMat_Private(ds,DS_MAT_A,ld));
826
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
63 PetscCall(DSReallocateMat_Private(ds,DS_MAT_Q,ld));
827
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.
63 PetscCall(DSReallocateMat_Private(ds,DS_MAT_T,ld));
828
829
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.
63 PetscCall(PetscMalloc1(ld,&ds->perm));
830
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.
63 PetscCall(PetscArraycpy(ds->perm,perm,ds->ld));
831
6/8
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 8 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
63 PetscCall(PetscFree(perm));
832
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.
13 PetscFunctionReturn(PETSC_SUCCESS);
833 }
834
835 /*MC
836 DSHEP - Dense Hermitian Eigenvalue Problem.
837
838 Notes:
839 The problem is expressed as $AX = X\Lambda$, where $A$ is real symmetric
840 (or complex Hermitian). $\Lambda$ is a diagonal matrix whose diagonal
841 elements are the arguments of `DSSolve()`. After solve, $A$ is overwritten
842 with $\Lambda$.
843
844 In the intermediate state $A$ is reduced to tridiagonal form. In compact
845 storage format, the symmetric tridiagonal matrix is stored in $T$.
846
847 Used DS matrices:
848 + `DS_MAT_A` - problem matrix (used only if `compact=PETSC_FALSE`)
849 . `DS_MAT_T` - symmetric tridiagonal matrix
850 - `DS_MAT_Q` - orthogonal/unitary transformation that reduces to tridiagonal form
851 (intermediate step) or matrix of orthogonal eigenvectors, which is equal to $X$
852
853 Implemented methods:
854 + 0 - Implicit QR (`_steqr`)
855 . 1 - Multiple Relatively Robust Representations (`_stevr`)
856 . 2 - Divide and Conquer (`_stedc`)
857 - 3 - Block Divide and Conquer (real scalars only)
858
859 Level: beginner
860
861 .seealso: [](sec:ds), `DSCreate()`, `DSSetType()`, `DSType`, `DSSetCompact()`
862 M*/
863 3711 SLEPC_EXTERN PetscErrorCode DSCreate_HEP(DS ds)
864 {
865
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
3711 PetscFunctionBegin;
866 3711 ds->ops->allocate = DSAllocate_HEP;
867 3711 ds->ops->view = DSView_HEP;
868 3711 ds->ops->vectors = DSVectors_HEP;
869 3711 ds->ops->solve[0] = DSSolve_HEP_QR;
870 3711 ds->ops->solve[1] = DSSolve_HEP_MRRR;
871 3711 ds->ops->solve[2] = DSSolve_HEP_DC;
872 #if !defined(PETSC_USE_COMPLEX)
873 1947 ds->ops->solve[3] = DSSolve_HEP_BDC;
874 #endif
875 3711 ds->ops->sort = DSSort_HEP;
876 3711 ds->ops->truncate = DSTruncate_HEP;
877 3711 ds->ops->update = DSUpdateExtraRow_HEP;
878 3711 ds->ops->cond = DSCond_HEP;
879 3711 ds->ops->transrks = DSTranslateRKS_HEP;
880 3711 ds->ops->hermitian = DSHermitian_HEP;
881 #if !defined(PETSC_HAVE_MPIUNI)
882 3711 ds->ops->synchronize = DSSynchronize_HEP;
883 #endif
884 3711 ds->ops->setcompact = DSSetCompact_HEP;
885 3711 ds->ops->reallocate = DSReallocate_HEP;
886
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.
3711 PetscFunctionReturn(PETSC_SUCCESS);
887 }
888