GCC Code Coverage Report


Directory: ./
File: src/sys/classes/ds/impls/ghiep/dsghiep.c
Date: 2025-10-03 04:28:47
Exec Total Coverage
Lines: 706 783 90.2%
Functions: 17 17 100.0%
Branches: 1520 2769 54.9%

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 579 static PetscErrorCode DSAllocate_GHIEP(DS ds,PetscInt ld)
15 {
16
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
579 PetscFunctionBegin;
17
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.
579 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.
579 PetscCall(DSAllocateMat_Private(ds,DS_MAT_B));
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.
579 PetscCall(DSAllocateMat_Private(ds,DS_MAT_Q));
20
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.
579 PetscCall(DSAllocateMat_Private(ds,DS_MAT_T));
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.
579 PetscCall(DSAllocateMat_Private(ds,DS_MAT_D));
22
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.
579 PetscCall(PetscFree(ds->perm));
23
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.
579 PetscCall(PetscMalloc1(ld,&ds->perm));
24
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.
121 PetscFunctionReturn(PETSC_SUCCESS);
25 }
26
27 10181 PetscErrorCode DSSwitchFormat_GHIEP(DS ds,PetscBool tocompact)
28 {
29 10181 PetscReal *T,*S;
30 10181 PetscScalar *A,*B;
31 10181 PetscInt i,n,ld;
32
33
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
10181 PetscFunctionBegin;
34
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.
10181 PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
35
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.
10181 PetscCall(MatDenseGetArray(ds->omat[DS_MAT_B],&B));
36
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.
10181 PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
37
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.
10181 PetscCall(DSGetArrayReal(ds,DS_MAT_D,&S));
38 10181 n = ds->n;
39 10181 ld = ds->ld;
40
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
10181 if (tocompact) { /* switch from dense (arrow) to compact storage */
41
4/6
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
60 PetscCall(PetscArrayzero(T,n));
42
4/6
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
60 PetscCall(PetscArrayzero(T+ld,n));
43
4/6
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
60 PetscCall(PetscArrayzero(T+2*ld,n));
44
4/6
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
60 PetscCall(PetscArrayzero(S,n));
45
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
600 for (i=0;i<n-1;i++) {
46 540 T[i] = PetscRealPart(A[i+i*ld]);
47 540 T[ld+i] = PetscRealPart(A[i+1+i*ld]);
48 540 S[i] = PetscRealPart(B[i+i*ld]);
49 }
50 60 T[n-1] = PetscRealPart(A[n-1+(n-1)*ld]);
51 60 S[n-1] = PetscRealPart(B[n-1+(n-1)*ld]);
52
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
60 for (i=ds->l;i<ds->k;i++) T[2*ld+i] = PetscRealPart(A[ds->k+i*ld]);
53 } else { /* switch from compact (arrow) to dense storage */
54
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
339958 for (i=0;i<n;i++) {
55
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.
329837 PetscCall(PetscArrayzero(A+i*ld,n));
56
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.
401921 PetscCall(PetscArrayzero(B+i*ld,n));
57 }
58
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
329837 for (i=0;i<n-1;i++) {
59 319716 A[i+i*ld] = T[i];
60 319716 A[i+1+i*ld] = T[ld+i];
61 319716 A[i+(i+1)*ld] = T[ld+i];
62 319716 B[i+i*ld] = S[i];
63 }
64 10121 A[n-1+(n-1)*ld] = T[n-1];
65 10121 B[n-1+(n-1)*ld] = S[n-1];
66
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
219727 for (i=ds->l;i<ds->k;i++) {
67 209606 A[ds->k+i*ld] = T[2*ld+i];
68 209606 A[i+ds->k*ld] = T[2*ld+i];
69 }
70 }
71
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.
10181 PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
72
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.
10181 PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_B],&B));
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.
10181 PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
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.
10181 PetscCall(DSRestoreArrayReal(ds,DS_MAT_D,&S));
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.
2019 PetscFunctionReturn(PETSC_SUCCESS);
76 }
77
78 180 static PetscErrorCode DSView_GHIEP(DS ds,PetscViewer viewer)
79 {
80 180 PetscViewerFormat format;
81 180 PetscInt i,j;
82 180 PetscReal *T,*S,value;
83 180 const char *methodname[] = {
84 "QR + Inverse Iteration",
85 "HZ method",
86 "QR"
87 };
88 180 const int nmeth=PETSC_STATIC_ARRAY_LENGTH(methodname);
89
90
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
180 PetscFunctionBegin;
91
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
180 PetscCall(PetscViewerGetFormat(viewer,&format));
92
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
180 if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
93
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.
120 if (ds->method<nmeth) PetscCall(PetscViewerASCIIPrintf(viewer,"solving the problem with: %s\n",methodname[ds->method]));
94
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.
120 PetscFunctionReturn(PETSC_SUCCESS);
95 }
96
1/2
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
60 if (ds->compact) {
97
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));
98
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_D,&S));
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(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE));
100
1/2
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
60 if (format == PETSC_VIEWER_ASCII_MATLAB) {
101
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",ds->n,ds->n));
102
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));
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,"zzz = [\n"));
104
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.
510 for (i=0;i<ds->n;i++) PetscCall(PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %18.16e\n",i+1,i+1,(double)T[i]));
105
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
450 for (i=0;i<ds->n-1;i++) {
106
3/4
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
✓ Branch 2 taken 10 times.
✗ Branch 3 not taken.
390 if (T[i+ds->ld] !=0 && i!=ds->k-1) {
107
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
120 PetscCall(PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %18.16e\n",i+2,i+1,(double)T[i+ds->ld]));
108
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.
390 PetscCall(PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %18.16e\n",i+1,i+2,(double)T[i+ds->ld]));
109 }
110 }
111
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
150 for (i = ds->l;i<ds->k;i++) {
112
1/2
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
90 if (T[i+2*ds->ld]) {
113
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.
90 PetscCall(PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %18.16e\n",ds->k+1,i+1,(double)T[i+2*ds->ld]));
114
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
90 PetscCall(PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %18.16e\n",i+1,ds->k+1,(double)T[i+2*ds->ld]));
115 }
116 }
117
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_A]));
118
119
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
60 PetscCall(PetscViewerASCIIPrintf(viewer,"%% Size = %" PetscInt_FMT " %" PetscInt_FMT "\n",ds->n,ds->n));
120
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
60 PetscCall(PetscViewerASCIIPrintf(viewer,"omega = zeros(%" PetscInt_FMT ",3);\n",3*ds->n));
121
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
60 PetscCall(PetscViewerASCIIPrintf(viewer,"omega = [\n"));
122
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.
510 for (i=0;i<ds->n;i++) PetscCall(PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %18.16e\n",i+1,i+1,(double)S[i]));
123
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(omega);\n",DSMatName[DS_MAT_B]));
124
125 } else {
126 PetscCall(PetscViewerASCIIPrintf(viewer,"T\n"));
127 for (i=0;i<ds->n;i++) {
128 for (j=0;j<ds->n;j++) {
129 if (i==j) value = T[i];
130 else if (i==j+1 || j==i+1) value = T[PetscMin(i,j)+ds->ld];
131 else if ((i<ds->k && j==ds->k) || (i==ds->k && j<ds->k)) value = T[PetscMin(i,j)+2*ds->ld];
132 else value = 0.0;
133 PetscCall(PetscViewerASCIIPrintf(viewer," %18.16e ",(double)value));
134 }
135 PetscCall(PetscViewerASCIIPrintf(viewer,"\n"));
136 }
137 PetscCall(PetscViewerASCIIPrintf(viewer,"omega\n"));
138 for (i=0;i<ds->n;i++) {
139 for (j=0;j<ds->n;j++) {
140 if (i==j) value = S[i];
141 else value = 0.0;
142 PetscCall(PetscViewerASCIIPrintf(viewer," %18.16e ",(double)value));
143 }
144 PetscCall(PetscViewerASCIIPrintf(viewer,"\n"));
145 }
146 }
147
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));
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.
60 PetscCall(PetscViewerFlush(viewer));
149
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));
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.
60 PetscCall(DSRestoreArrayReal(ds,DS_MAT_D,&S));
151 } else {
152 PetscCall(DSViewMat(ds,viewer,DS_MAT_A));
153 PetscCall(DSViewMat(ds,viewer,DS_MAT_B));
154 }
155
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));
156
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);
157 }
158
159 17642 static PetscErrorCode DSVectors_GHIEP_Eigen_Some(DS ds,PetscInt *idx,PetscReal *rnorm)
160 {
161 17642 PetscReal b[4],M[4],*T,*S,d1,d2,s1,s2,e,scal1,scal2,wr1,wr2,wi,ep,norm;
162 17642 PetscScalar *X,Y[4],alpha,szero=0.0;
163 17642 const PetscScalar *A,*B,*Q;
164 17642 PetscInt k;
165 17642 PetscBLASInt two=2,n_,ld,one=1;
166 #if !defined(PETSC_USE_COMPLEX)
167 12519 PetscBLASInt four=4;
168 #endif
169
170
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
17642 PetscFunctionBegin;
171
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.
17642 PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_A],&A));
172
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.
17642 PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_B],&B));
173
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.
17642 PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_Q],&Q));
174
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
17642 PetscCall(MatDenseGetArray(ds->omat[DS_MAT_X],&X));
175
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
17642 PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
176
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
17642 PetscCall(DSGetArrayReal(ds,DS_MAT_D,&S));
177 17642 k = *idx;
178
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.
17642 PetscCall(PetscBLASIntCast(ds->n,&n_));
179
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.
17642 PetscCall(PetscBLASIntCast(ds->ld,&ld));
180
2/4
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 10 times.
✗ Branch 3 not taken.
17642 if (k < ds->n-1) e = (ds->compact)?T[k+ld]:PetscRealPart(A[(k+1)+ld*k]);
181 else e = 0.0;
182
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
17642 if (e == 0.0) { /* Real */
183
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.
16840 if (ds->state>=DS_STATE_CONDENSED) PetscCall(PetscArraycpy(X+k*ld,Q+k*ld,ld));
184 else {
185 PetscCall(PetscArrayzero(X+k*ds->ld,ds->ld));
186 X[k+k*ds->ld] = 1.0;
187 }
188
1/2
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
16840 if (rnorm) *rnorm = PetscAbsScalar(X[ds->n-1+k*ld]);
189 } else { /* 2x2 block */
190
1/2
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
802 if (ds->compact) {
191 802 s1 = S[k];
192 802 d1 = T[k];
193 802 s2 = S[k+1];
194 802 d2 = T[k+1];
195 } else {
196 s1 = PetscRealPart(B[k*ld+k]);
197 d1 = PetscRealPart(A[k+k*ld]);
198 s2 = PetscRealPart(B[(k+1)*ld+k+1]);
199 d2 = PetscRealPart(A[k+1+(k+1)*ld]);
200 }
201 802 M[0] = d1; M[1] = e; M[2] = e; M[3]= d2;
202 802 b[0] = s1; b[1] = 0.0; b[2] = 0.0; b[3] = s2;
203 802 ep = LAPACKlamch_("S");
204 /* Compute eigenvalues of the block */
205
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.
802 PetscCallBLAS("LAPACKlag2",LAPACKlag2_(M,&two,b,&two,&ep,&scal1,&scal2,&wr1,&wr2,&wi));
206
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
802 PetscCheck(wi!=0.0,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Real block in DSVectors_GHIEP");
207 /* Complex eigenvalues */
208
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
802 PetscCheck(scal1>=ep,PETSC_COMM_SELF,PETSC_ERR_FP,"Nearly infinite eigenvalue");
209 802 wr1 /= scal1;
210 802 wi /= scal1;
211 #if !defined(PETSC_USE_COMPLEX)
212
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
582 if (SlepcAbs(s1*d1-wr1,wi)<SlepcAbs(s2*d2-wr1,wi)) {
213 84 Y[0] = wr1-s2*d2; Y[1] = s2*e; Y[2] = wi; Y[3] = 0.0;
214 } else {
215 498 Y[0] = s1*e; Y[1] = wr1-s1*d1; Y[2] = 0.0; Y[3] = wi;
216 }
217 582 norm = BLASnrm2_(&four,Y,&one);
218 582 norm = 1.0/norm;
219
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
582 if (ds->state >= DS_STATE_CONDENSED) {
220 455 alpha = norm;
221
10/20
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 1 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.
455 PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n_,&two,&two,&alpha,Q+k*ld,&ld,Y,&two,&szero,X+k*ld,&ld));
222
1/2
✓ Branch 0 taken 5 times.
✗ Branch 1 not taken.
455 if (rnorm) *rnorm = SlepcAbsEigenvalue(X[ds->n-1+k*ld],X[ds->n-1+(k+1)*ld]);
223 } else {
224
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.
127 PetscCall(PetscArrayzero(X+k*ld,2*ld));
225 127 X[k*ld+k] = Y[0]*norm;
226 127 X[k*ld+k+1] = Y[1]*norm;
227 127 X[(k+1)*ld+k] = Y[2]*norm;
228 127 X[(k+1)*ld+k+1] = Y[3]*norm;
229 }
230 #else
231
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 5 times.
220 if (SlepcAbs(s1*d1-wr1,wi)<SlepcAbs(s2*d2-wr1,wi)) {
232 22 Y[0] = PetscCMPLX(wr1-s2*d2,wi);
233 22 Y[1] = s2*e;
234 } else {
235 198 Y[0] = s1*e;
236 198 Y[1] = PetscCMPLX(wr1-s1*d1,wi);
237 }
238 220 norm = BLASnrm2_(&two,Y,&one);
239 220 norm = 1.0/norm;
240
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
220 if (ds->state >= DS_STATE_CONDENSED) {
241 153 alpha = norm;
242
10/20
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 1 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.
153 PetscCallBLAS("BLASgemv",BLASgemv_("N",&n_,&two,&alpha,Q+k*ld,&ld,Y,&one,&szero,X+k*ld,&one));
243
1/2
✓ Branch 0 taken 5 times.
✗ Branch 1 not taken.
153 if (rnorm) *rnorm = PetscAbsScalar(X[ds->n-1+k*ld]);
244 } else {
245
4/6
✓ Branch 0 taken 3 times.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
67 PetscCall(PetscArrayzero(X+k*ld,2*ld));
246 67 X[k*ld+k] = Y[0]*norm;
247 67 X[k*ld+k+1] = Y[1]*norm;
248 }
249 220 X[(k+1)*ld+k] = PetscConj(X[k*ld+k]);
250 220 X[(k+1)*ld+k+1] = PetscConj(X[k*ld+k+1]);
251 #endif
252 802 (*idx)++;
253 }
254
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.
17642 PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_A],&A));
255
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.
17642 PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_B],&B));
256
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.
17642 PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_Q],&Q));
257
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.
17642 PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_X],&X));
258
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
17642 PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
259
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.
17642 PetscCall(DSRestoreArrayReal(ds,DS_MAT_D,&S));
260
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.
3883 PetscFunctionReturn(PETSC_SUCCESS);
261 }
262
263 17734 static PetscErrorCode DSVectors_GHIEP(DS ds,DSMatType mat,PetscInt *k,PetscReal *rnorm)
264 {
265 17734 PetscScalar *Z;
266 17734 const PetscScalar *A,*Q;
267 17734 PetscInt i;
268 17734 PetscReal e,*T,*S;
269
270
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
17734 PetscFunctionBegin;
271
1/3
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
17734 switch (mat) {
272 17734 case DS_MAT_X:
273 case DS_MAT_Y:
274
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.
17734 if (k) PetscCall(DSVectors_GHIEP_Eigen_Some(ds,k,rnorm));
275 else {
276
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
286 PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_A],&A));
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.
286 PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_Q],&Q));
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.
286 PetscCall(MatDenseGetArray(ds->omat[mat],&Z));
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.
286 PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
280
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
4096 for (i=0; i<ds->n; i++) {
281
1/2
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
3810 e = (ds->compact)?T[i+ds->ld]:PetscRealPart(A[(i+1)+ds->ld*i]);
282
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
3810 if (e == 0.0) { /* real */
283
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.
3616 if (ds->state >= DS_STATE_CONDENSED) PetscCall(PetscArraycpy(Z+i*ds->ld,Q+i*ds->ld,ds->ld));
284 else {
285
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.
3616 PetscCall(PetscArrayzero(Z+i*ds->ld,ds->ld));
286 3616 Z[i+i*ds->ld] = 1.0;
287 }
288 } else {
289
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.
194 PetscCall(DSVectors_GHIEP_Eigen_Some(ds,&i,rnorm));
290
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.
194 PetscCall(DSGetArrayReal(ds,DS_MAT_D,&S));
291 194 S[i]=S[i-1]=1.0;
292
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.
3810 PetscCall(DSRestoreArrayReal(ds,DS_MAT_D,&S));
293 }
294 }
295
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.
286 PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_A],&A));
296
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.
286 PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_Q],&Q));
297
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
286 PetscCall(MatDenseRestoreArray(ds->omat[mat],&Z));
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.
286 PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
299 }
300 3903 break;
301 case DS_MAT_U:
302 case DS_MAT_V:
303 SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"Not implemented yet");
304 default:
305 SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Invalid mat parameter");
306 }
307
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.
3903 PetscFunctionReturn(PETSC_SUCCESS);
308 }
309
310 /*
311 Extract the eigenvalues contained in the block-diagonal of the indefinite problem.
312 Only the index range n0..n1 is processed.
313 */
314 10560 PetscErrorCode DSGHIEPComplexEigs(DS ds,PetscInt n0,PetscInt n1,PetscScalar *wr,PetscScalar *wi)
315 {
316 10560 PetscInt k,ld;
317 10560 PetscBLASInt two=2;
318 10560 const PetscScalar *A,*B;
319 10560 PetscReal *D,*T,b[4],M[4],d1,d2,s1,s2,e,scal1,scal2,ep,wr1,wr2,wi1;
320
321
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
10560 PetscFunctionBegin;
322 10560 ld = ds->ld;
323
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.
10560 PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_A],&A));
324
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.
10560 PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_B],&B));
325
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.
10560 PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
326
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.
10560 PetscCall(DSGetArrayReal(ds,DS_MAT_D,&D));
327
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
22581 for (k=n0;k<n1;k++) {
328
4/4
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
✓ Branch 2 taken 10 times.
✓ Branch 3 taken 10 times.
12021 if (k < n1-1) e = (ds->compact)?T[ld+k]:PetscRealPart(A[(k+1)+ld*k]);
329 else e = 0.0;
330
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
11048 if (e==0.0) { /* real eigenvalue */
331
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
11966 wr[k] = (ds->compact)?T[k]/D[k]:A[k+k*ld]/B[k+k*ld];
332 #if !defined(PETSC_USE_COMPLEX)
333 6563 wi[k] = 0.0 ;
334 #endif
335 } else { /* diagonal block */
336
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
55 if (ds->compact) {
337 35 s1 = D[k];
338 35 d1 = T[k];
339 35 s2 = D[k+1];
340 35 d2 = T[k+1];
341 } else {
342 20 s1 = PetscRealPart(B[k*ld+k]);
343 20 d1 = PetscRealPart(A[k+k*ld]);
344 20 s2 = PetscRealPart(B[(k+1)*ld+k+1]);
345 20 d2 = PetscRealPart(A[k+1+(k+1)*ld]);
346 }
347 55 M[0] = d1; M[1] = e; M[2] = e; M[3]= d2;
348 55 b[0] = s1; b[1] = 0.0; b[2] = 0.0; b[3] = s2;
349 55 ep = LAPACKlamch_("S");
350 /* Compute eigenvalues of the block */
351
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.
55 PetscCallBLAS("LAPACKlag2",LAPACKlag2_(M,&two,b,&two,&ep,&scal1,&scal2,&wr1,&wr2,&wi1));
352
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
55 PetscCheck(scal1>=ep,PETSC_COMM_SELF,PETSC_ERR_FP,"Nearly infinite eigenvalue");
353
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
55 if (wi1==0.0) { /* Real eigenvalues */
354 PetscCheck(scal2>=ep,PETSC_COMM_SELF,PETSC_ERR_FP,"Nearly infinite eigenvalue");
355 wr[k] = wr1/scal1; wr[k+1] = wr2/scal2;
356 #if !defined(PETSC_USE_COMPLEX)
357 wi[k] = wi[k+1] = 0.0;
358 #endif
359 } else { /* Complex eigenvalues */
360 #if !defined(PETSC_USE_COMPLEX)
361 30 wr[k] = wr1/scal1;
362 30 wr[k+1] = wr[k];
363 30 wi[k] = wi1/scal1;
364 30 wi[k+1] = -wi[k];
365 #else
366 25 wr[k] = PetscCMPLX(wr1,wi1)/scal1;
367 25 wr[k+1] = PetscConj(wr[k]);
368 #endif
369 }
370 55 k++;
371 }
372 }
373 #if defined(PETSC_USE_COMPLEX)
374
1/2
✓ Branch 0 taken 5 times.
✗ Branch 1 not taken.
1816 if (wi) {
375
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
7269 for (k=n0;k<n1;k++) wi[k] = 0.0;
376 }
377 #endif
378
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.
10560 PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_A],&A));
379
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.
10560 PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_B],&B));
380
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.
10560 PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
381
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.
10560 PetscCall(DSRestoreArrayReal(ds,DS_MAT_D,&D));
382
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.
2100 PetscFunctionReturn(PETSC_SUCCESS);
383 }
384
385 10560 static PetscErrorCode DSSort_GHIEP(DS ds,PetscScalar *wr,PetscScalar *wi,PetscScalar *rr,PetscScalar *ri,PetscInt *k)
386 {
387 10560 PetscInt n,i,*perm;
388 10560 PetscReal *d,*e,*s;
389
390
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
10560 PetscFunctionBegin;
391 #if !defined(PETSC_USE_COMPLEX)
392
2/8
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
8744 PetscAssertPointer(wi,3);
393 #endif
394 10560 n = ds->n;
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.
10560 PetscCall(DSGetArrayReal(ds,DS_MAT_T,&d));
396 10560 e = d + ds->ld;
397
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.
10560 PetscCall(DSGetArrayReal(ds,DS_MAT_D,&s));
398
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.
10560 PetscCall(DSAllocateWork_Private(ds,ds->ld,ds->ld,0));
399 10560 perm = ds->perm;
400
1/2
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
10560 if (!rr) {
401 10560 rr = wr;
402 10560 ri = wi;
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.
10560 PetscCall(DSSortEigenvalues_Private(ds,rr,ri,perm,PETSC_TRUE));
405
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.
10560 if (!ds->compact) PetscCall(DSSwitchFormat_GHIEP(ds,PETSC_TRUE));
406
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.
10560 PetscCall(PetscArraycpy(ds->work,wr,n));
407
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
346556 for (i=ds->l;i<n;i++) wr[i] = *(ds->work+perm[i]);
408 #if !defined(PETSC_USE_COMPLEX)
409
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.
8744 PetscCall(PetscArraycpy(ds->work,wi,n));
410
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
296405 for (i=ds->l;i<n;i++) wi[i] = *(ds->work+perm[i]);
411 #endif
412
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.
10560 PetscCall(PetscArraycpy(ds->rwork,s,n));
413
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
346556 for (i=ds->l;i<n;i++) s[i] = *(ds->rwork+perm[i]);
414
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.
10560 PetscCall(PetscArraycpy(ds->rwork,d,n));
415
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
346556 for (i=ds->l;i<n;i++) d[i] = *(ds->rwork+perm[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.
10560 PetscCall(PetscArraycpy(ds->rwork,e,n-1));
417
4/6
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
10560 PetscCall(PetscArrayzero(e+ds->l,n-1-ds->l));
418
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
335996 for (i=ds->l;i<n-1;i++) {
419
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
325436 if (perm[i]<n-1) e[i] = *(ds->rwork+perm[i]);
420 }
421
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.
10560 if (!ds->compact) PetscCall(DSSwitchFormat_GHIEP(ds,PETSC_FALSE));
422
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.
10560 PetscCall(DSPermuteColumns_Private(ds,ds->l,n,n,DS_MAT_Q,perm));
423
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.
10560 PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&d));
424
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.
10560 PetscCall(DSRestoreArrayReal(ds,DS_MAT_D,&s));
425
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.
2100 PetscFunctionReturn(PETSC_SUCCESS);
426 }
427
428 10470 static PetscErrorCode DSUpdateExtraRow_GHIEP(DS ds)
429 {
430 10470 PetscInt i;
431 10470 PetscBLASInt n,ld,incx=1;
432 10470 PetscScalar *A,*x,*y,one=1.0,zero=0.0;
433 10470 const PetscScalar *Q;
434 10470 PetscReal *T,*b,*r,beta;
435
436
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
10470 PetscFunctionBegin;
437
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
10470 PetscCall(PetscBLASIntCast(ds->n,&n));
438
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.
10470 PetscCall(PetscBLASIntCast(ds->ld,&ld));
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.
10470 PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_Q],&Q));
440
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
10470 if (ds->compact) {
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.
10440 PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
442 10440 b = T+ld;
443 10440 r = T+2*ld;
444 10440 beta = b[n-1]; /* in compact, we assume all entries are zero except the last one */
445
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
357182 for (i=0;i<n;i++) r[i] = PetscRealPart(beta*Q[n-1+i*ld]);
446 10440 ds->k = n;
447
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.
10440 PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
448 } else {
449
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));
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.
30 PetscCall(DSAllocateWork_Private(ds,2*ld,0,0));
451 30 x = ds->work;
452 30 y = ds->work+ld;
453
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
330 for (i=0;i<n;i++) x[i] = PetscConj(A[n+i*ld]);
454
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));
455
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
330 for (i=0;i<n;i++) A[n+i*ld] = PetscConj(y[i]);
456 30 ds->k = n;
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.
30 PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
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.
10470 PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_Q],&Q));
460
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.
2082 PetscFunctionReturn(PETSC_SUCCESS);
461 }
462
463 /*
464 Get eigenvectors with inverse iteration.
465 The system matrix is in Hessenberg form.
466 */
467 10470 PetscErrorCode DSGHIEPInverseIteration(DS ds,PetscScalar *wr,PetscScalar *wi)
468 {
469 10470 PetscInt i,off;
470 10470 PetscBLASInt *select,*infoC,ld,n1,mout,info;
471 10470 const PetscScalar *A,*B;
472 10470 PetscScalar *H,*X;
473 10470 PetscReal *s,*d,*e;
474 #if defined(PETSC_USE_COMPLEX)
475 1771 PetscInt j;
476 #endif
477
478
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
10470 PetscFunctionBegin;
479
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.
10470 PetscCall(PetscBLASIntCast(ds->ld,&ld));
480
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.
10470 PetscCall(PetscBLASIntCast(ds->n-ds->l,&n1));
481
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
10470 PetscCall(DSAllocateWork_Private(ds,ld*ld+2*ld,ld,2*ld));
482
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.
10470 PetscCall(DSAllocateMat_Private(ds,DS_MAT_W));
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.
10470 PetscCall(MatDenseGetArrayRead(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.
10470 PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_B],&B));
485
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.
10470 PetscCall(MatDenseGetArrayWrite(ds->omat[DS_MAT_W],&H));
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.
10470 PetscCall(DSGetArrayReal(ds,DS_MAT_T,&d));
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.
10470 PetscCall(DSGetArrayReal(ds,DS_MAT_D,&s));
488 10470 e = d + ld;
489 10470 select = ds->iwork;
490 10470 infoC = ds->iwork + ld;
491 10470 off = ds->l+ds->l*ld;
492
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
10470 if (ds->compact) {
493 10450 H[off] = d[ds->l]*s[ds->l];
494 10450 H[off+ld] = e[ds->l]*s[ds->l];
495
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
324756 for (i=ds->l+1;i<ds->n-1;i++) {
496 314306 H[i+(i-1)*ld] = e[i-1]*s[i];
497 314306 H[i+i*ld] = d[i]*s[i];
498 314306 H[i+(i+1)*ld] = e[i]*s[i];
499 }
500 10450 H[ds->n-1+(ds->n-2)*ld] = e[ds->n-2]*s[ds->n-1];
501 10450 H[ds->n-1+(ds->n-1)*ld] = d[ds->n-1]*s[ds->n-1];
502 } else {
503 20 s[ds->l] = PetscRealPart(B[off]);
504 20 H[off] = A[off]*s[ds->l];
505 20 H[off+ld] = A[off+ld]*s[ds->l];
506
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
180 for (i=ds->l+1;i<ds->n-1;i++) {
507 160 s[i] = PetscRealPart(B[i+i*ld]);
508 160 H[i+(i-1)*ld] = A[i+(i-1)*ld]*s[i];
509 160 H[i+i*ld] = A[i+i*ld]*s[i];
510 160 H[i+(i+1)*ld] = A[i+(i+1)*ld]*s[i];
511 }
512 20 s[ds->n-1] = PetscRealPart(B[ds->n-1+(ds->n-1)*ld]);
513 20 H[ds->n-1+(ds->n-2)*ld] = A[ds->n-1+(ds->n-2)*ld]*s[ds->n-1];
514 20 H[ds->n-1+(ds->n-1)*ld] = A[ds->n-1+(ds->n-1)*ld]*s[ds->n-1];
515 }
516
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
10470 PetscCall(DSAllocateMat_Private(ds,DS_MAT_X));
517
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
10470 PetscCall(MatDenseGetArray(ds->omat[DS_MAT_X],&X));
518
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
345876 for (i=0;i<n1;i++) select[i] = 1;
519 #if !defined(PETSC_USE_COMPLEX)
520
10/20
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 1 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.
8699 PetscCallBLAS("LAPACKhsein",LAPACKhsein_("R","N","N",select,&n1,H+off,&ld,wr+ds->l,wi+ds->l,NULL,&ld,X+off,&ld,&n1,&mout,ds->work,NULL,infoC,&info));
521 #else
522
10/20
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 1 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.
1771 PetscCallBLAS("LAPACKhsein",LAPACKhsein_("R","N","N",select,&n1,H+off,&ld,wr+ds->l,NULL,&ld,X+off,&ld,&n1,&mout,ds->work,ds->rwork,NULL,infoC,&info));
523
524 /* Separate real and imaginary part of complex eigenvectors */
525
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
48652 for (j=ds->l;j<ds->n;j++) {
526
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
46881 if (PetscAbsReal(PetscImaginaryPart(wr[j])) > PetscAbsScalar(wr[j])*PETSC_SQRT_MACHINE_EPSILON) {
527
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
65916 for (i=ds->l;i<ds->n;i++) {
528 64757 X[i+(j+1)*ds->ld] = PetscImaginaryPart(X[i+j*ds->ld]);
529 64757 X[i+j*ds->ld] = PetscRealPart(X[i+j*ds->ld]);
530 }
531 1159 j++;
532 }
533 }
534 #endif
535
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
10470 SlepcCheckLapackInfo("hsein",info);
536
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.
10470 PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_A],&A));
537
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.
10470 PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_B],&B));
538
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.
10470 PetscCall(MatDenseRestoreArrayWrite(ds->omat[DS_MAT_W],&H));
539
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.
10470 PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_X],&X));
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.
10470 PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&d));
541
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.
10470 PetscCall(DSRestoreArrayReal(ds,DS_MAT_D,&s));
542
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.
10470 PetscCall(DSGHIEPOrthogEigenv(ds,DS_MAT_X,wr,wi,PETSC_TRUE));
543
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.
2082 PetscFunctionReturn(PETSC_SUCCESS);
544 }
545
546 /*
547 Undo 2x2 blocks that have real eigenvalues.
548 */
549 30 PetscErrorCode DSGHIEPRealBlocks(DS ds)
550 {
551 30 PetscInt i;
552 30 PetscReal e,d1,d2,s1,s2,ss1,ss2,t,dd,ss;
553 30 PetscReal maxy,ep,scal1,scal2,snorm;
554 30 PetscReal *T,*D,b[4],M[4],wr1,wr2,wi;
555 30 PetscScalar *A,*B,*Q,Y[4],sone=1.0,szero=0.0;
556 30 PetscBLASInt m,two=2,ld;
557 30 PetscBool isreal;
558
559
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
30 PetscFunctionBegin;
560
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
30 PetscCall(PetscBLASIntCast(ds->ld,&ld));
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.
30 PetscCall(PetscBLASIntCast(ds->n-ds->l,&m));
562
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));
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.
30 PetscCall(MatDenseGetArray(ds->omat[DS_MAT_B],&B));
564
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
30 PetscCall(MatDenseGetArray(ds->omat[DS_MAT_Q],&Q));
565
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(DSGetArrayReal(ds,DS_MAT_T,&T));
566
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(DSGetArrayReal(ds,DS_MAT_D,&D));
567
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*m,0,0));
568
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
230 for (i=ds->l;i<ds->n-1;i++) {
569
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
200 e = (ds->compact)?T[ld+i]:PetscRealPart(A[(i+1)+ld*i]);
570
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
200 if (e != 0.0) { /* 2x2 block */
571
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
50 if (ds->compact) {
572 10 s1 = D[i];
573 10 d1 = T[i];
574 10 s2 = D[i+1];
575 10 d2 = T[i+1];
576 } else {
577 40 s1 = PetscRealPart(B[i*ld+i]);
578 40 d1 = PetscRealPart(A[i*ld+i]);
579 40 s2 = PetscRealPart(B[(i+1)*ld+i+1]);
580 40 d2 = PetscRealPart(A[(i+1)*ld+i+1]);
581 }
582 50 isreal = PETSC_FALSE;
583
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
50 if (s1==s2) { /* apply a Jacobi rotation to compute the eigendecomposition */
584 dd = d1-d2;
585 if (2*PetscAbsReal(e) <= dd) {
586 t = 2*e/dd;
587 t = t/(1 + PetscSqrtReal(1+t*t));
588 } else {
589 t = dd/(2*e);
590 ss = (t>=0)?1.0:-1.0;
591 t = ss/(PetscAbsReal(t)+PetscSqrtReal(1+t*t));
592 }
593 Y[0] = 1/PetscSqrtReal(1 + t*t); Y[3] = Y[0]; /* c */
594 Y[1] = Y[0]*t; Y[2] = -Y[1]; /* s */
595 wr1 = d1+t*e; wr2 = d2-t*e;
596 ss1 = s1; ss2 = s2;
597 isreal = PETSC_TRUE;
598 } else {
599 50 ss1 = 1.0; ss2 = 1.0,
600 50 M[0] = d1; M[1] = e; M[2] = e; M[3]= d2;
601 50 b[0] = s1; b[1] = 0.0; b[2] = 0.0; b[3] = s2;
602 50 ep = LAPACKlamch_("S");
603
604 /* Compute eigenvalues of the block */
605
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.
50 PetscCallBLAS("LAPACKlag2",LAPACKlag2_(M,&two,b,&two,&ep,&scal1,&scal2,&wr1,&wr2,&wi));
606
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
50 if (wi==0.0) { /* Real eigenvalues */
607 20 isreal = PETSC_TRUE;
608
2/6
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 10 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
20 PetscCheck(scal1>=ep && scal2>=ep,PETSC_COMM_SELF,PETSC_ERR_FP,"Nearly infinite eigenvalue");
609 20 wr1 /= scal1;
610 20 wr2 /= scal2;
611
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
20 if (PetscAbsReal(s1*d1-wr1)<PetscAbsReal(s2*d2-wr1)) {
612 Y[0] = wr1-s2*d2;
613 Y[1] = s2*e;
614 } else {
615 20 Y[0] = s1*e;
616 20 Y[1] = wr1-s1*d1;
617 }
618 /* normalize with a signature*/
619
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
20 maxy = PetscMax(PetscAbsScalar(Y[0]),PetscAbsScalar(Y[1]));
620 20 scal1 = PetscRealPart(Y[0])/maxy;
621 20 scal2 = PetscRealPart(Y[1])/maxy;
622 20 snorm = scal1*scal1*s1 + scal2*scal2*s2;
623
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
20 if (snorm<0) { ss1 = -1.0; snorm = -snorm; }
624 20 snorm = maxy*PetscSqrtReal(snorm);
625 20 Y[0] = Y[0]/snorm;
626 20 Y[1] = Y[1]/snorm;
627
1/2
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
20 if (PetscAbsReal(s1*d1-wr2)<PetscAbsReal(s2*d2-wr2)) {
628 20 Y[2] = wr2-s2*d2;
629 20 Y[3] = s2*e;
630 } else {
631 Y[2] = s1*e;
632 Y[3] = wr2-s1*d1;
633 }
634
1/2
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
20 maxy = PetscMax(PetscAbsScalar(Y[2]),PetscAbsScalar(Y[3]));
635 20 scal1 = PetscRealPart(Y[2])/maxy;
636 20 scal2 = PetscRealPart(Y[3])/maxy;
637 20 snorm = scal1*scal1*s1 + scal2*scal2*s2;
638
1/2
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
20 if (snorm<0) { ss2 = -1.0; snorm = -snorm; }
639 20 snorm = maxy*PetscSqrtReal(snorm); Y[2] = Y[2]/snorm; Y[3] = Y[3]/snorm;
640 }
641 50 wr1 *= ss1; wr2 *= ss2;
642 }
643
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
50 if (isreal) {
644
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
20 if (ds->compact) {
645 D[i] = ss1;
646 T[i] = wr1;
647 D[i+1] = ss2;
648 T[i+1] = wr2;
649 T[ld+i] = 0.0;
650 } else {
651 20 B[i*ld+i] = ss1;
652 20 A[i*ld+i] = wr1;
653 20 B[(i+1)*ld+i+1] = ss2;
654 20 A[(i+1)*ld+i+1] = wr2;
655 20 A[(i+1)+ld*i] = 0.0;
656 20 A[i+ld*(i+1)] = 0.0;
657 }
658
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.
20 PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&m,&two,&two,&sone,Q+ds->l+i*ld,&ld,Y,&two,&szero,ds->work,&m));
659
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
20 PetscCall(PetscArraycpy(Q+ds->l+i*ld,ds->work,m));
660
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
20 PetscCall(PetscArraycpy(Q+ds->l+(i+1)*ld,ds->work+m,m));
661 }
662 50 i++;
663 }
664 }
665
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));
666
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_B],&B));
667
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
30 PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_Q],&Q));
668
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
30 PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
669
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
30 PetscCall(DSRestoreArrayReal(ds,DS_MAT_D,&D));
670
6/12
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 2 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 2 times.
6 PetscFunctionReturn(PETSC_SUCCESS);
671 }
672
673 10480 static PetscErrorCode DSSolve_GHIEP_QR_II(DS ds,PetscScalar *wr,PetscScalar *wi)
674 {
675 10480 PetscInt i,off;
676 10480 PetscBLASInt n1,ld,one=1,info,lwork;
677 10480 const PetscScalar *A,*B;
678 10480 PetscScalar *H,*Q;
679 10480 PetscReal *d,*e,*s;
680 #if defined(PETSC_USE_COMPLEX)
681 1776 PetscInt j;
682 #endif
683
684
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
10480 PetscFunctionBegin;
685 #if !defined(PETSC_USE_COMPLEX)
686
2/8
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
8704 PetscAssertPointer(wi,3);
687 #endif
688
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.
10480 PetscCall(PetscBLASIntCast(ds->n-ds->l,&n1));
689
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.
10480 PetscCall(PetscBLASIntCast(ds->ld,&ld));
690 10480 off = ds->l + ds->l*ld;
691
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.
10480 PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_A],&A));
692
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.
10480 PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_B],&B));
693
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.
10480 PetscCall(DSGetArrayReal(ds,DS_MAT_T,&d));
694
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.
10480 PetscCall(DSGetArrayReal(ds,DS_MAT_D,&s));
695 10480 e = d + ld;
696 #if defined(PETSC_USE_DEBUG)
697 /* Check signature */
698
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
72320 for (i=0;i<ds->n;i++) {
699
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
70236 PetscReal de = (ds->compact)?s[i]:PetscRealPart(B[i*ld+i]);
700
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
70236 PetscCheck(de==1.0 || de==-1.0,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Diagonal elements of the signature matrix must be 1 or -1");
701 }
702 #endif
703
704 /* Quick return if possible */
705
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
10480 if (n1 == 1) {
706
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
10 PetscCall(MatDenseGetArray(ds->omat[DS_MAT_Q],&Q));
707
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
60 for (i=0;i<=ds->l;i++) Q[i+i*ld] = 1.0;
708
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
10 PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_Q],&Q));
709
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
10 PetscCall(DSGHIEPComplexEigs(ds,0,ds->l,wr,wi));
710
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
10 if (!ds->compact) {
711 d[ds->l] = PetscRealPart(A[off]);
712 s[ds->l] = PetscRealPart(B[off]);
713 }
714
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
10 PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_A],&A));
715
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
10 PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_B],&B));
716 10 wr[ds->l] = d[ds->l]/s[ds->l];
717
1/2
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
10 if (wi) wi[ds->l] = 0.0;
718
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
10 PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&d));
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.
10 PetscCall(DSRestoreArrayReal(ds,DS_MAT_D,&s));
720
6/12
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 2 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 2 times.
2 PetscFunctionReturn(PETSC_SUCCESS);
721 }
722
723
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.
10470 PetscCall(DSAllocateWork_Private(ds,ld*ld,2*ld,ld*2));
724 10470 lwork = ld*ld;
725
726 /* Reduce to pseudotriadiagonal form */
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.
10470 PetscCall(DSIntermediate_GHIEP(ds));
728
729 /* Compute Eigenvalues (QR) */
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.
10470 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.
10470 PetscCall(MatDenseGetArrayWrite(ds->omat[DS_MAT_W],&H));
732
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
10470 if (ds->compact) {
733 10450 H[off] = d[ds->l]*s[ds->l];
734 10450 H[off+ld] = e[ds->l]*s[ds->l];
735
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
324756 for (i=ds->l+1;i<ds->n-1;i++) {
736 314306 H[i+(i-1)*ld] = e[i-1]*s[i];
737 314306 H[i+i*ld] = d[i]*s[i];
738 314306 H[i+(i+1)*ld] = e[i]*s[i];
739 }
740 10450 H[ds->n-1+(ds->n-2)*ld] = e[ds->n-2]*s[ds->n-1];
741 10450 H[ds->n-1+(ds->n-1)*ld] = d[ds->n-1]*s[ds->n-1];
742 } else {
743 20 s[ds->l] = PetscRealPart(B[off]);
744 20 H[off] = A[off]*s[ds->l];
745 20 H[off+ld] = A[off+ld]*s[ds->l];
746
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
180 for (i=ds->l+1;i<ds->n-1;i++) {
747 160 s[i] = PetscRealPart(B[i+i*ld]);
748 160 H[i+(i-1)*ld] = A[i+(i-1)*ld]*s[i];
749 160 H[i+i*ld] = A[i+i*ld]*s[i];
750 160 H[i+(i+1)*ld] = A[i+(i+1)*ld]*s[i];
751 }
752 20 s[ds->n-1] = PetscRealPart(B[ds->n-1+(ds->n-1)*ld]);
753 20 H[ds->n-1+(ds->n-2)*ld] = A[ds->n-1+(ds->n-2)*ld]*s[ds->n-1];
754 20 H[ds->n-1+(ds->n-1)*ld] = A[ds->n-1+(ds->n-1)*ld]*s[ds->n-1];
755 }
756
757 #if !defined(PETSC_USE_COMPLEX)
758
10/20
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 1 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.
8699 PetscCallBLAS("LAPACKhseqr",LAPACKhseqr_("E","N",&n1,&one,&n1,H+off,&ld,wr+ds->l,wi+ds->l,NULL,&ld,ds->work,&lwork,&info));
759 #else
760
10/20
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 1 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.
1771 PetscCallBLAS("LAPACKhseqr",LAPACKhseqr_("E","N",&n1,&one,&n1,H+off,&ld,wr+ds->l,NULL,&ld,ds->work,&lwork,&info));
761
4/4
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
✓ Branch 2 taken 5 times.
✓ Branch 3 taken 5 times.
49811 for (i=ds->l;i<ds->n;i++) if (PetscAbsReal(PetscImaginaryPart(wr[i]))<10*PETSC_MACHINE_EPSILON) wr[i] = PetscRealPart(wr[i]);
762 /* Sort to have consecutive conjugate pairs */
763
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
48601 for (i=ds->l;i<ds->n;i++) {
764 46830 j=i+1;
765
4/4
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
✓ Branch 2 taken 5 times.
✓ Branch 3 taken 5 times.
1127556 while (j<ds->n && (PetscAbsScalar(wr[i]-PetscConj(wr[j]))>PetscAbsScalar(wr[i])*PETSC_SQRT_MACHINE_EPSILON)) j++;
766
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
46830 if (j==ds->n) {
767
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
45620 PetscCheck(PetscAbsReal(PetscImaginaryPart(wr[i]))<PetscAbsScalar(wr[i])*PETSC_SQRT_MACHINE_EPSILON,PETSC_COMM_SELF,PETSC_ERR_LIB,"Found complex without conjugate pair");
768 45620 wr[i] = PetscRealPart(wr[i]);
769 } else { /* complex eigenvalue */
770 1210 wr[j] = wr[i+1];
771
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
1210 if (PetscImaginaryPart(wr[i])<0) wr[i] = PetscConj(wr[i]);
772 1210 wr[i+1] = PetscConj(wr[i]);
773 1210 i++;
774 }
775 }
776 #endif
777
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
10470 SlepcCheckLapackInfo("hseqr",info);
778
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.
10470 PetscCall(MatDenseRestoreArrayWrite(ds->omat[DS_MAT_W],&H));
779
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.
10470 PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_A],&A));
780
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.
10470 PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_B],&B));
781
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.
10470 PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&d));
782
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.
10470 PetscCall(DSRestoreArrayReal(ds,DS_MAT_D,&s));
783
784 /* Compute Eigenvectors with Inverse Iteration */
785
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.
10470 PetscCall(DSGHIEPInverseIteration(ds,wr,wi));
786
787 /* Recover eigenvalues from diagonal */
788
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.
10470 PetscCall(DSGHIEPComplexEigs(ds,0,ds->l,wr,wi));
789 #if defined(PETSC_USE_COMPLEX)
790
1/2
✓ Branch 0 taken 5 times.
✗ Branch 1 not taken.
1771 if (wi) {
791
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
49811 for (i=ds->l;i<ds->n;i++) wi[i] = 0.0;
792 }
793 #endif
794
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.
2082 PetscFunctionReturn(PETSC_SUCCESS);
795 }
796
797 40 static PetscErrorCode DSSolve_GHIEP_QR(DS ds,PetscScalar *wr,PetscScalar *wi)
798 {
799 40 PetscInt i,j,off,nwu=0,n,lw,lwr,nwru=0;
800 40 PetscBLASInt n_,ld,info,lwork,ilo,ihi;
801 40 const PetscScalar *A,*B;
802 40 PetscScalar *H,*Q,*X;
803 40 PetscReal *d,*s,*scale,nrm,*rcde,*rcdv;
804 #if defined(PETSC_USE_COMPLEX)
805 20 PetscInt k;
806 #endif
807
808
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
40 PetscFunctionBegin;
809 #if !defined(PETSC_USE_COMPLEX)
810
2/8
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
20 PetscAssertPointer(wi,3);
811 #endif
812 40 n = ds->n-ds->l;
813
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,&n_));
814
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));
815 40 off = ds->l + ds->l*ld;
816
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(MatDenseGetArrayRead(ds->omat[DS_MAT_A],&A));
817
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(MatDenseGetArrayRead(ds->omat[DS_MAT_B],&B));
818
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));
819
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
40 PetscCall(DSGetArrayReal(ds,DS_MAT_D,&s));
820 #if defined(PETSC_USE_DEBUG)
821 /* Check signature */
822
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
78 for (i=0;i<ds->n;i++) {
823
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
70 PetscReal de = (ds->compact)?s[i]:PetscRealPart(B[i*ld+i]);
824
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
70 PetscCheck(de==1.0 || de==-1.0,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Diagonal elements of the signature matrix must be 1 or -1");
825 }
826 #endif
827
828 /* Quick return if possible */
829
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
40 if (n_ == 1) {
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.
10 PetscCall(MatDenseGetArray(ds->omat[DS_MAT_Q],&Q));
831
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
60 for (i=0;i<=ds->l;i++) Q[i+i*ld] = 1.0;
832
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
10 PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_Q],&Q));
833
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
10 PetscCall(DSGHIEPComplexEigs(ds,0,ds->l,wr,wi));
834
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
10 if (!ds->compact) {
835 d[ds->l] = PetscRealPart(A[off]);
836 s[ds->l] = PetscRealPart(B[off]);
837 }
838
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
10 PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_A],&A));
839
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
10 PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_B],&B));
840 10 wr[ds->l] = d[ds->l]/s[ds->l];
841
1/2
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
10 if (wi) wi[ds->l] = 0.0;
842
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
10 PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&d));
843
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
10 PetscCall(DSRestoreArrayReal(ds,DS_MAT_D,&s));
844
6/12
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 2 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 2 times.
2 PetscFunctionReturn(PETSC_SUCCESS);
845 }
846
847 30 lw = 14*ld+ld*ld;
848 30 lwr = 7*ld;
849
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,lw,lwr,0));
850 30 scale = ds->rwork+nwru;
851 30 nwru += ld;
852 30 rcde = ds->rwork+nwru;
853 30 nwru += ld;
854 30 rcdv = ds->rwork+nwru;
855
856 /* Form pseudo-symmetric matrix */
857 30 H = ds->work+nwu;
858 30 nwu += n*n;
859
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.
30 PetscCall(PetscArrayzero(H,n*n));
860
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
30 if (ds->compact) {
861
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
80 for (i=0;i<n-1;i++) {
862 70 H[i+i*n] = s[ds->l+i]*d[ds->l+i];
863 70 H[i+1+i*n] = s[ds->l+i+1]*d[ld+ds->l+i];
864 70 H[i+(i+1)*n] = s[ds->l+i]*d[ld+ds->l+i];
865 }
866 10 H[n-1+(n-1)*n] = s[ds->l+n-1]*d[ds->l+n-1];
867
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
40 for (i=0;i<ds->k-ds->l;i++) {
868 30 H[ds->k-ds->l+i*n] = s[ds->k]*d[2*ld+ds->l+i];
869 30 H[i+(ds->k-ds->l)*n] = s[i+ds->l]*d[2*ld+ds->l+i];
870 }
871 } else {
872
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
220 for (j=0;j<n;j++) {
873
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
2200 for (i=0;i<n;i++) H[i+j*n] = B[off+i+i*ld]*A[off+i+j*ld];
874 }
875 }
876
877 /* Compute eigenpairs */
878
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
30 PetscCall(PetscBLASIntCast(lw-nwu,&lwork));
879
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(DSAllocateMat_Private(ds,DS_MAT_X));
880
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(MatDenseGetArrayWrite(ds->omat[DS_MAT_X],&X));
881 #if !defined(PETSC_USE_COMPLEX)
882
10/20
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 1 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.
15 PetscCallBLAS("LAPACKgeevx",LAPACKgeevx_("B","N","V","N",&n_,H,&n_,wr+ds->l,wi+ds->l,NULL,&ld,X+off,&ld,&ilo,&ihi,scale,&nrm,rcde,rcdv,ds->work+nwu,&lwork,NULL,&info));
883 #else
884
10/20
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 1 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.
15 PetscCallBLAS("LAPACKgeevx",LAPACKgeevx_("B","N","V","N",&n_,H,&n_,wr+ds->l,NULL,&ld,X+off,&ld,&ilo,&ihi,scale,&nrm,rcde,rcdv,ds->work+nwu,&lwork,ds->rwork+nwru,&info));
885
886 /* Sort to have consecutive conjugate pairs
887 Separate real and imaginary part of complex eigenvectors*/
888
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
140 for (i=ds->l;i<ds->n;i++) {
889 125 j=i+1;
890
4/4
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
✓ Branch 2 taken 5 times.
✓ Branch 3 taken 5 times.
480 while (j<ds->n && (PetscAbsScalar(wr[i]-PetscConj(wr[j]))>PetscAbsScalar(wr[i])*PETSC_SQRT_MACHINE_EPSILON)) j++;
891
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
125 if (j==ds->n) {
892
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
110 PetscCheck(PetscAbsReal(PetscImaginaryPart(wr[i]))<PetscAbsScalar(wr[i])*PETSC_SQRT_MACHINE_EPSILON,PETSC_COMM_SELF,PETSC_ERR_LIB,"Found complex without conjugate pair");
893 110 wr[i]=PetscRealPart(wr[i]); /* real eigenvalue */
894
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
1150 for (k=ds->l;k<ds->n;k++) {
895 1040 X[k+i*ds->ld] = PetscRealPart(X[k+i*ds->ld]);
896 }
897 } else { /* complex eigenvalue */
898
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
15 if (j!=i+1) {
899 wr[j] = wr[i+1];
900 PetscCall(PetscArraycpy(X+j*ds->ld,X+(i+1)*ds->ld,ds->ld));
901 }
902
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
15 if (PetscImaginaryPart(wr[i])<0) {
903 wr[i] = PetscConj(wr[i]);
904 for (k=ds->l;k<ds->n;k++) {
905 X[k+(i+1)*ds->ld] = -PetscImaginaryPart(X[k+i*ds->ld]);
906 X[k+i*ds->ld] = PetscRealPart(X[k+i*ds->ld]);
907 }
908 } else {
909
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
155 for (k=ds->l;k<ds->n;k++) {
910 140 X[k+(i+1)*ds->ld] = PetscImaginaryPart(X[k+i*ds->ld]);
911 140 X[k+i*ds->ld] = PetscRealPart(X[k+i*ds->ld]);
912 }
913 }
914 15 wr[i+1] = PetscConj(wr[i]);
915 15 i++;
916 }
917 }
918 #endif
919
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
30 SlepcCheckLapackInfo("geevx",info);
920
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
30 PetscCall(MatDenseRestoreArrayWrite(ds->omat[DS_MAT_X],&X));
921
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(MatDenseRestoreArrayRead(ds->omat[DS_MAT_A],&A));
922
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(MatDenseRestoreArrayRead(ds->omat[DS_MAT_B],&B));
923
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
30 PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&d));
924
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
30 PetscCall(DSRestoreArrayReal(ds,DS_MAT_D,&s));
925
926 /* Compute real s-orthonormal basis */
927
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(DSGHIEPOrthogEigenv(ds,DS_MAT_X,wr,wi,PETSC_FALSE));
928
929 /* Recover eigenvalues from diagonal */
930
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
30 PetscCall(DSGHIEPComplexEigs(ds,0,ds->l,wr,wi));
931 #if defined(PETSC_USE_COMPLEX)
932
1/2
✓ Branch 0 taken 5 times.
✗ Branch 1 not taken.
15 if (wi) {
933
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
155 for (i=ds->l;i<ds->n;i++) wi[i] = 0.0;
934 }
935 #endif
936
6/12
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 2 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 2 times.
6 PetscFunctionReturn(PETSC_SUCCESS);
937 }
938
939 9944 static PetscErrorCode DSGetTruncateSize_GHIEP(DS ds,PetscInt l,PetscInt n,PetscInt *k)
940 {
941 9944 PetscReal *T;
942
943
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
9944 PetscFunctionBegin;
944
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.
9944 PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
945
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
9944 if (T[l+(*k)-1+ds->ld] !=0.0) {
946
1/2
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
5965 if (l+(*k)<n-1) (*k)++;
947 else (*k)--;
948 }
949
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.
9944 PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
950
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.
1972 PetscFunctionReturn(PETSC_SUCCESS);
951 }
952
953 10443 static PetscErrorCode DSTruncate_GHIEP(DS ds,PetscInt n,PetscBool trim)
954 {
955 10443 PetscInt i,ld=ds->ld,l=ds->l;
956 10443 PetscScalar *A;
957 10443 PetscReal *T,*b,*r,*omega;
958
959
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
10443 PetscFunctionBegin;
960
1/2
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
10443 if (ds->compact) {
961
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.
10443 PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
962
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.
10443 PetscCall(DSGetArrayReal(ds,DS_MAT_D,&omega));
963 #if defined(PETSC_USE_DEBUG)
964 /* make sure diagonal 2x2 block is not broken */
965
6/10
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 2 times.
✓ Branch 4 taken 2 times.
✓ Branch 5 taken 2 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
2077 PetscCheck(ds->state<DS_STATE_CONDENSED || n==0 || n==ds->n || T[n-1+ld]==0.0,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"The given size would break a 2x2 block, call DSGetTruncateSize() first");
966 #endif
967 }
968
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
10443 if (trim) {
969
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
499 if (!ds->compact && ds->extrarow) { /* clean extra row */
970 PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
971 for (i=l;i<ds->n;i++) A[ds->n+i*ld] = 0.0;
972 PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
973 }
974 499 ds->l = 0;
975 499 ds->k = 0;
976 499 ds->n = n;
977 499 ds->t = ds->n; /* truncated length equal to the new dimension */
978 } else {
979
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.
9944 if (!ds->compact && ds->extrarow && ds->k==ds->n) {
980 /* copy entries of extra row to the new position, then clean last row */
981 PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
982 for (i=l;i<n;i++) A[n+i*ld] = A[ds->n+i*ld];
983 for (i=l;i<ds->n;i++) A[ds->n+i*ld] = 0.0;
984 PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
985 }
986
1/2
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
9944 if (ds->compact) {
987 9944 b = T+ld;
988 9944 r = T+2*ld;
989 9944 b[n-1] = r[n-1];
990 9944 b[n] = b[ds->n];
991 9944 omega[n] = omega[ds->n];
992 }
993
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
9944 ds->k = (ds->extrarow)? n: 0;
994 9944 ds->t = ds->n; /* truncated length equal to previous dimension */
995 9944 ds->n = n;
996 }
997
1/2
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
10443 if (ds->compact) {
998
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.
10443 PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
999
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.
10443 PetscCall(DSRestoreArrayReal(ds,DS_MAT_D,&omega));
1000 }
1001
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.
2077 PetscFunctionReturn(PETSC_SUCCESS);
1002 }
1003
1004 #if !defined(PETSC_HAVE_MPIUNI)
1005 20 static PetscErrorCode DSSynchronize_GHIEP(DS ds,PetscScalar eigr[],PetscScalar eigi[])
1006 {
1007 20 PetscScalar *A,*B,*Q;
1008 20 PetscReal *T,*D;
1009 20 PetscInt ld=ds->ld,l=ds->l,k=0,kr=0;
1010 20 PetscMPIInt n,rank,off=0,size,ldn,ld3,ld_;
1011
1012
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
20 PetscFunctionBegin;
1013
1/2
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
20 if (ds->compact) kr = 4*ld;
1014 else k = 2*(ds->n-l)*ld;
1015
1/2
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
20 if (ds->state>DS_STATE_RAW) k += (ds->n-l)*ld;
1016
1/2
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
20 if (eigr) k += (ds->n-l);
1017
1/2
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
20 if (eigi) k += (ds->n-l);
1018
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));
1019
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));
1020
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));
1021
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));
1022
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));
1023
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,&ld_));
1024
1/2
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
20 if (ds->compact) {
1025
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(DSGetArrayReal(ds,DS_MAT_T,&T));
1026
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(DSGetArrayReal(ds,DS_MAT_D,&D));
1027 } else {
1028 PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
1029 PetscCall(MatDenseGetArray(ds->omat[DS_MAT_B],&B));
1030 }
1031
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));
1032
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));
1033
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
20 if (!rank) {
1034
1/2
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
10 if (ds->compact) {
1035
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.
10 PetscCallMPI(MPI_Pack(T,ld3,MPIU_REAL,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
1036
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.
10 PetscCallMPI(MPI_Pack(D,ld_,MPIU_REAL,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
1037 } else {
1038 PetscCallMPI(MPI_Pack(A+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
1039 PetscCallMPI(MPI_Pack(B+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
1040 }
1041
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)));
1042
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)));
1043 #if !defined(PETSC_USE_COMPLEX)
1044
15/30
✓ Branch 0 taken 5 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1 times.
✓ Branch 3 taken 4 times.
✓ Branch 4 taken 1 times.
✓ Branch 5 taken 4 times.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 1 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✓ Branch 13 taken 1 times.
✓ Branch 14 taken 1 times.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✓ Branch 17 taken 1 times.
✓ Branch 18 taken 1 times.
✗ Branch 19 not taken.
✓ Branch 20 taken 1 times.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✓ Branch 23 taken 1 times.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✓ Branch 27 taken 1 times.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
5 if (eigi) PetscCallMPI(MPI_Pack(eigi+l,n,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
1045 #endif
1046 }
1047
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)));
1048
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
20 if (rank) {
1049
1/2
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
10 if (ds->compact) {
1050
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.
10 PetscCallMPI(MPI_Unpack(ds->work,size,&off,T,ld3,MPIU_REAL,PetscObjectComm((PetscObject)ds)));
1051
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.
10 PetscCallMPI(MPI_Unpack(ds->work,size,&off,D,ld_,MPIU_REAL,PetscObjectComm((PetscObject)ds)));
1052 } else {
1053 PetscCallMPI(MPI_Unpack(ds->work,size,&off,A+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
1054 PetscCallMPI(MPI_Unpack(ds->work,size,&off,B+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
1055 }
1056
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)));
1057
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)));
1058 #if !defined(PETSC_USE_COMPLEX)
1059
15/30
✓ Branch 0 taken 5 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1 times.
✓ Branch 3 taken 4 times.
✓ Branch 4 taken 1 times.
✓ Branch 5 taken 4 times.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 1 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✓ Branch 13 taken 1 times.
✓ Branch 14 taken 1 times.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✓ Branch 17 taken 1 times.
✓ Branch 18 taken 1 times.
✗ Branch 19 not taken.
✓ Branch 20 taken 1 times.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✓ Branch 23 taken 1 times.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✓ Branch 27 taken 1 times.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
5 if (eigi) PetscCallMPI(MPI_Unpack(ds->work,size,&off,eigi+l,n,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
1060 #endif
1061 }
1062
1/2
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
20 if (ds->compact) {
1063
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(DSRestoreArrayReal(ds,DS_MAT_T,&T));
1064
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(DSRestoreArrayReal(ds,DS_MAT_D,&D));
1065 } else {
1066 PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
1067 PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_B],&B));
1068 }
1069
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));
1070
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);
1071 }
1072 #endif
1073
1074 21789 static PetscErrorCode DSHermitian_GHIEP(DS ds,DSMatType m,PetscBool *flg)
1075 {
1076
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
21789 PetscFunctionBegin;
1077
2/6
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 10 times.
21789 if ((m==DS_MAT_A && !ds->extrarow) || m==DS_MAT_B) *flg = PETSC_TRUE;
1078 21789 else *flg = PETSC_FALSE;
1079
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.
21789 PetscFunctionReturn(PETSC_SUCCESS);
1080 }
1081
1082 /*MC
1083 DSGHIEP - Dense Generalized Hermitian Indefinite Eigenvalue Problem.
1084
1085 Level: beginner
1086
1087 Notes:
1088 The problem is expressed as A*X = B*X*Lambda, where both A and B are
1089 real symmetric (or complex Hermitian) and possibly indefinite. Lambda
1090 is a diagonal matrix whose diagonal elements are the arguments of DSSolve().
1091 After solve, A is overwritten with Lambda. Note that in the case of real
1092 scalars, A is overwritten with a real representation of Lambda, i.e.,
1093 complex conjugate eigenvalue pairs are stored as a 2x2 block in the
1094 quasi-diagonal matrix.
1095
1096 In the intermediate state A is reduced to tridiagonal form and B is
1097 transformed into a signature matrix. In compact storage format, these
1098 matrices are stored in T and D, respectively.
1099
1100 Used DS matrices:
1101 + DS_MAT_A - first problem matrix
1102 . DS_MAT_B - second problem matrix
1103 . DS_MAT_T - symmetric tridiagonal matrix of the reduced pencil
1104 . DS_MAT_D - diagonal matrix (signature) of the reduced pencil
1105 - DS_MAT_Q - pseudo-orthogonal transformation that reduces (A,B) to
1106 tridiagonal-diagonal form (intermediate step) or a real basis of eigenvectors
1107
1108 Implemented methods:
1109 + 0 - QR iteration plus inverse iteration for the eigenvectors
1110 . 1 - HZ iteration
1111 - 2 - QR iteration plus pseudo-orthogonalization for the eigenvectors
1112
1113 References:
1114 . 1. - C. Campos and J. E. Roman, "Restarted Q-Arnoldi-type methods exploiting
1115 symmetry in quadratic eigenvalue problems", BIT Numer. Math. 56(4):1213-1236, 2016.
1116
1117 .seealso: DSCreate(), DSSetType(), DSType
1118 M*/
1119 476 SLEPC_EXTERN PetscErrorCode DSCreate_GHIEP(DS ds)
1120 {
1121
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
476 PetscFunctionBegin;
1122 476 ds->ops->allocate = DSAllocate_GHIEP;
1123 476 ds->ops->view = DSView_GHIEP;
1124 476 ds->ops->vectors = DSVectors_GHIEP;
1125 476 ds->ops->solve[0] = DSSolve_GHIEP_QR_II;
1126 476 ds->ops->solve[1] = DSSolve_GHIEP_HZ;
1127 476 ds->ops->solve[2] = DSSolve_GHIEP_QR;
1128 476 ds->ops->sort = DSSort_GHIEP;
1129 #if !defined(PETSC_HAVE_MPIUNI)
1130 476 ds->ops->synchronize = DSSynchronize_GHIEP;
1131 #endif
1132 476 ds->ops->gettruncatesize = DSGetTruncateSize_GHIEP;
1133 476 ds->ops->truncate = DSTruncate_GHIEP;
1134 476 ds->ops->update = DSUpdateExtraRow_GHIEP;
1135 476 ds->ops->hermitian = DSHermitian_GHIEP;
1136
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.
476 PetscFunctionReturn(PETSC_SUCCESS);
1137 }
1138