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 |