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 | Utility subroutines common to several impls | ||
12 | */ | ||
13 | |||
14 | #include <slepc/private/dsimpl.h> /*I "slepcds.h" I*/ | ||
15 | #include <slepcblaslapack.h> | ||
16 | |||
17 | /* | ||
18 | Compute the (real) Schur form of A. At the end, A is (quasi-)triangular and Q | ||
19 | contains the unitary matrix of Schur vectors. Eigenvalues are returned in wr,wi | ||
20 | */ | ||
21 | 39823 | PetscErrorCode DSSolve_NHEP_Private(DS ds,DSMatType mA,DSMatType mQ,PetscScalar *wr,PetscScalar *wi) | |
22 | { | ||
23 | 39823 | PetscScalar *work,*tau,*A,*Q; | |
24 | 39823 | PetscInt i,j; | |
25 | 39823 | PetscBLASInt ilo,lwork,info,n,k,ld; | |
26 | |||
27 |
1/2✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
|
39823 | PetscFunctionBegin; |
28 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
39823 | PetscCall(MatDenseGetArray(ds->omat[mA],&A)); |
29 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
39823 | PetscCall(MatDenseGetArray(ds->omat[mQ],&Q)); |
30 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
39823 | PetscCall(PetscBLASIntCast(ds->n,&n)); |
31 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
39823 | PetscCall(PetscBLASIntCast(ds->ld,&ld)); |
32 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
39823 | PetscCall(PetscBLASIntCast(ds->l+1,&ilo)); |
33 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
39823 | PetscCall(PetscBLASIntCast(ds->k,&k)); |
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.
|
39823 | PetscCall(DSAllocateWork_Private(ds,ld+6*ld,0,0)); |
35 | 39823 | tau = ds->work; | |
36 | 39823 | work = ds->work+ld; | |
37 | 39823 | lwork = 6*ld; | |
38 | |||
39 | /* initialize orthogonal matrix */ | ||
40 |
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.
|
39823 | PetscCall(PetscArrayzero(Q,ld*ld)); |
41 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
764051 | for (i=0;i<n;i++) Q[i+i*ld] = 1.0; |
42 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
39823 | if (n==1) { /* quick return */ |
43 | 460 | wr[0] = A[0]; | |
44 |
1/2✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
|
460 | if (wi) wi[0] = 0.0; |
45 |
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.
|
460 | PetscFunctionReturn(PETSC_SUCCESS); |
46 | } | ||
47 | |||
48 | /* reduce to upper Hessenberg form */ | ||
49 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
39363 | if (ds->state<DS_STATE_INTERMEDIATE) { |
50 |
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.
|
26777 | PetscCallBLAS("LAPACKgehrd",LAPACKgehrd_(&n,&ilo,&n,A,&ld,tau,work,&lwork,&info)); |
51 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
26777 | SlepcCheckLapackInfo("gehrd",info); |
52 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
497946 | for (j=0;j<n-1;j++) { |
53 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
5950127 | for (i=j+2;i<n;i++) { |
54 | 5478958 | Q[i+j*ld] = A[i+j*ld]; | |
55 | 5478958 | A[i+j*ld] = 0.0; | |
56 | } | ||
57 | } | ||
58 |
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.
|
26777 | PetscCallBLAS("LAPACKorghr",LAPACKorghr_(&n,&ilo,&n,Q,&ld,tau,work,&lwork,&info)); |
59 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
26777 | SlepcCheckLapackInfo("orghr",info); |
60 | } | ||
61 | |||
62 | /* compute the (real) Schur form */ | ||
63 | #if !defined(PETSC_USE_COMPLEX) | ||
64 |
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.
|
19729 | PetscCallBLAS("LAPACKhseqr",LAPACKhseqr_("S","V",&n,&ilo,&n,A,&ld,wr,wi,Q,&ld,work,&lwork,&info)); |
65 |
2/2✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
|
38633 | for (j=0;j<ds->l;j++) { |
66 |
3/4✓ Branch 0 taken 5 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 5 times.
✓ Branch 3 taken 5 times.
|
18904 | if (j==n-1 || A[j+1+j*ld] == 0.0) { |
67 | /* real eigenvalue */ | ||
68 | 17159 | wr[j] = A[j+j*ld]; | |
69 | 17159 | wi[j] = 0.0; | |
70 | } else { | ||
71 | /* complex eigenvalue */ | ||
72 | 1745 | wr[j] = A[j+j*ld]; | |
73 | 1745 | wr[j+1] = A[j+j*ld]; | |
74 | 1745 | wi[j] = PetscSqrtReal(PetscAbsReal(A[j+1+j*ld]))*PetscSqrtReal(PetscAbsReal(A[j+(j+1)*ld])); | |
75 | 1745 | wi[j+1] = -wi[j]; | |
76 | 1745 | j++; | |
77 | } | ||
78 | } | ||
79 | #else | ||
80 |
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.
|
19634 | PetscCallBLAS("LAPACKhseqr",LAPACKhseqr_("S","V",&n,&ilo,&n,A,&ld,wr,Q,&ld,work,&lwork,&info)); |
81 |
4/4✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
✓ Branch 2 taken 5 times.
✓ Branch 3 taken 5 times.
|
347403 | if (wi) for (i=ds->l;i<n;i++) wi[i] = 0.0; |
82 | #endif | ||
83 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
39363 | SlepcCheckLapackInfo("hseqr",info); |
84 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
39363 | PetscCall(MatDenseRestoreArray(ds->omat[mA],&A)); |
85 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
39363 | PetscCall(MatDenseRestoreArray(ds->omat[mQ],&Q)); |
86 |
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.
|
7578 | PetscFunctionReturn(PETSC_SUCCESS); |
87 | } | ||
88 | |||
89 | /* | ||
90 | Sort a Schur form represented by the (quasi-)triangular matrix T and | ||
91 | the unitary matrix Q, and return the sorted eigenvalues in wr,wi | ||
92 | */ | ||
93 | 43678 | PetscErrorCode DSSort_NHEP_Total(DS ds,DSMatType mT,DSMatType mQ,PetscScalar *wr,PetscScalar *wi) | |
94 | { | ||
95 | 43678 | PetscScalar re,*T,*Q; | |
96 | 43678 | PetscInt i,j,pos,result; | |
97 | 43678 | PetscBLASInt ifst,ilst,info,n,ld; | |
98 | #if !defined(PETSC_USE_COMPLEX) | ||
99 | 21911 | PetscScalar *work,im; | |
100 | #endif | ||
101 | |||
102 |
1/2✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
|
43678 | PetscFunctionBegin; |
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.
|
43678 | PetscCall(MatDenseGetArray(ds->omat[mT],&T)); |
104 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
43678 | PetscCall(MatDenseGetArray(ds->omat[mQ],&Q)); |
105 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
43678 | PetscCall(PetscBLASIntCast(ds->n,&n)); |
106 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
43678 | PetscCall(PetscBLASIntCast(ds->ld,&ld)); |
107 | #if !defined(PETSC_USE_COMPLEX) | ||
108 |
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.
|
21911 | PetscCall(DSAllocateWork_Private(ds,ld,0,0)); |
109 | 21911 | work = ds->work; | |
110 | #endif | ||
111 | /* selection sort */ | ||
112 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
686116 | for (i=ds->l;i<n-1;i++) { |
113 | 642438 | re = wr[i]; | |
114 | #if !defined(PETSC_USE_COMPLEX) | ||
115 | 300077 | im = wi[i]; | |
116 | #endif | ||
117 | 642438 | pos = 0; | |
118 | 642438 | j=i+1; /* j points to the next eigenvalue */ | |
119 | #if !defined(PETSC_USE_COMPLEX) | ||
120 |
2/2✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
|
300077 | if (im != 0) j=i+2; |
121 | #endif | ||
122 | /* find minimum eigenvalue */ | ||
123 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
7557684 | for (;j<n;j++) { |
124 | #if !defined(PETSC_USE_COMPLEX) | ||
125 |
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.
|
3131940 | PetscCall(SlepcSCCompare(ds->sc,re,im,wr[j],wi[j],&result)); |
126 | #else | ||
127 |
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.
|
3783306 | PetscCall(SlepcSCCompare(ds->sc,re,0.0,wr[j],0.0,&result)); |
128 | #endif | ||
129 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
6915246 | if (result > 0) { |
130 | 2610035 | re = wr[j]; | |
131 | #if !defined(PETSC_USE_COMPLEX) | ||
132 | 1222242 | im = wi[j]; | |
133 | #endif | ||
134 | 2610035 | pos = j; | |
135 | } | ||
136 | #if !defined(PETSC_USE_COMPLEX) | ||
137 |
2/2✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
|
3131940 | if (wi[j] != 0) j++; |
138 | #endif | ||
139 | } | ||
140 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
642438 | if (pos) { |
141 | /* interchange blocks */ | ||
142 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
436931 | PetscCall(PetscBLASIntCast(pos+1,&ifst)); |
143 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
436931 | PetscCall(PetscBLASIntCast(i+1,&ilst)); |
144 | #if !defined(PETSC_USE_COMPLEX) | ||
145 |
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.
|
201607 | PetscCallBLAS("LAPACKtrexc",LAPACKtrexc_("V",&n,T,&ld,Q,&ld,&ifst,&ilst,work,&info)); |
146 | #else | ||
147 |
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.
|
235324 | PetscCallBLAS("LAPACKtrexc",LAPACKtrexc_("V",&n,T,&ld,Q,&ld,&ifst,&ilst,&info)); |
148 | #endif | ||
149 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
436931 | SlepcCheckLapackInfo("trexc",info); |
150 | /* recover original eigenvalues from T matrix */ | ||
151 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
5798775 | for (j=i;j<n;j++) { |
152 | 5361844 | wr[j] = T[j+j*ld]; | |
153 | #if !defined(PETSC_USE_COMPLEX) | ||
154 |
4/4✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
✓ Branch 2 taken 5 times.
✓ Branch 3 taken 5 times.
|
2397764 | if (j<n-1 && T[j+1+j*ld] != 0.0) { |
155 | /* complex conjugate eigenvalue */ | ||
156 | 262242 | wi[j] = PetscSqrtReal(PetscAbsReal(T[j+1+j*ld]))*PetscSqrtReal(PetscAbsReal(T[j+(j+1)*ld])); | |
157 | 262242 | wr[j+1] = wr[j]; | |
158 | 262242 | wi[j+1] = -wi[j]; | |
159 | 262242 | j++; | |
160 | 2135522 | } else wi[j] = 0.0; | |
161 | #endif | ||
162 | } | ||
163 | } | ||
164 | #if !defined(PETSC_USE_COMPLEX) | ||
165 |
2/2✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
|
300077 | if (wi[i] != 0) i++; |
166 | #endif | ||
167 | } | ||
168 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
43678 | PetscCall(MatDenseRestoreArray(ds->omat[mT],&T)); |
169 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
43678 | PetscCall(MatDenseRestoreArray(ds->omat[mQ],&Q)); |
170 |
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.
|
8432 | PetscFunctionReturn(PETSC_SUCCESS); |
171 | } | ||
172 | |||
173 | /* | ||
174 | Reorder a Schur form represented by T,Q according to a permutation perm, | ||
175 | and return the sorted eigenvalues in wr,wi | ||
176 | */ | ||
177 | 187 | PetscErrorCode DSSortWithPermutation_NHEP_Private(DS ds,PetscInt *perm,DSMatType mT,DSMatType mQ,PetscScalar *wr,PetscScalar *wi) | |
178 | { | ||
179 | 187 | PetscInt i,j,pos,inc=1; | |
180 | 187 | PetscBLASInt ifst,ilst,info,n,ld; | |
181 | 187 | PetscScalar *T,*Q; | |
182 | #if !defined(PETSC_USE_COMPLEX) | ||
183 | 75 | PetscScalar *work; | |
184 | #endif | ||
185 | |||
186 |
1/2✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
|
187 | PetscFunctionBegin; |
187 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
187 | PetscCall(MatDenseGetArray(ds->omat[mT],&T)); |
188 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
187 | PetscCall(MatDenseGetArray(ds->omat[mQ],&Q)); |
189 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
187 | PetscCall(PetscBLASIntCast(ds->n,&n)); |
190 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
187 | PetscCall(PetscBLASIntCast(ds->ld,&ld)); |
191 | #if !defined(PETSC_USE_COMPLEX) | ||
192 |
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.
|
75 | PetscCall(DSAllocateWork_Private(ds,ld,0,0)); |
193 | 75 | work = ds->work; | |
194 | #endif | ||
195 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
4223 | for (i=ds->l;i<n-1;i++) { |
196 | 4036 | pos = perm[i]; | |
197 | #if !defined(PETSC_USE_COMPLEX) | ||
198 |
3/4✓ Branch 0 taken 5 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 5 times.
✓ Branch 3 taken 5 times.
|
2100 | inc = (pos<n-1 && T[pos+1+pos*ld] != 0.0)? 2: 1; |
199 | #endif | ||
200 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
4036 | if (pos!=i) { |
201 | #if !defined(PETSC_USE_COMPLEX) | ||
202 |
5/12✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 5 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 5 times.
✓ Branch 7 taken 5 times.
✓ Branch 8 taken 5 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
|
445 | PetscCheck((T[pos+(pos-1)*ld]==0.0 || perm[i+1]==pos-1) && (pos==n-1 || (T[pos+1+pos*ld]==0.0 || perm[i+1]==pos+1)),PETSC_COMM_SELF,PETSC_ERR_FP,"Invalid permutation due to a 2x2 block at position %" PetscInt_FMT,pos); |
203 | #endif | ||
204 | /* interchange blocks */ | ||
205 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
1352 | PetscCall(PetscBLASIntCast(pos+1,&ifst)); |
206 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
1352 | PetscCall(PetscBLASIntCast(i+1,&ilst)); |
207 | #if !defined(PETSC_USE_COMPLEX) | ||
208 |
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.
|
445 | PetscCallBLAS("LAPACKtrexc",LAPACKtrexc_("V",&n,T,&ld,Q,&ld,&ifst,&ilst,work,&info)); |
209 | #else | ||
210 |
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.
|
907 | PetscCallBLAS("LAPACKtrexc",LAPACKtrexc_("V",&n,T,&ld,Q,&ld,&ifst,&ilst,&info)); |
211 | #endif | ||
212 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
1352 | SlepcCheckLapackInfo("trexc",info); |
213 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
13584 | for (j=i+1;j<n;j++) { |
214 |
4/4✓ Branch 0 taken 10 times.
✓ Branch 1 taken 5 times.
✓ Branch 2 taken 10 times.
✓ Branch 3 taken 10 times.
|
12232 | if (perm[j]>=i && perm[j]<pos) perm[j]+=inc; |
215 | } | ||
216 | 1352 | perm[i] = i; | |
217 |
2/2✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
|
1352 | if (inc==2) perm[i+1] = i+1; |
218 | } | ||
219 |
2/2✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
|
3646 | if (inc==2) i++; |
220 | } | ||
221 | /* recover original eigenvalues from T matrix */ | ||
222 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
4365 | for (j=ds->l;j<n;j++) { |
223 | 4178 | wr[j] = T[j+j*ld]; | |
224 | #if !defined(PETSC_USE_COMPLEX) | ||
225 |
4/4✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
✓ Branch 2 taken 5 times.
✓ Branch 3 taken 5 times.
|
2130 | if (j<n-1 && T[j+1+j*ld] != 0.0) { |
226 | /* complex conjugate eigenvalue */ | ||
227 | 140 | wi[j] = PetscSqrtReal(PetscAbsReal(T[j+1+j*ld]))*PetscSqrtReal(PetscAbsReal(T[j+(j+1)*ld])); | |
228 | 140 | wr[j+1] = wr[j]; | |
229 | 140 | wi[j+1] = -wi[j]; | |
230 | 140 | j++; | |
231 | 1990 | } else wi[j] = 0.0; | |
232 | #endif | ||
233 | } | ||
234 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
187 | PetscCall(MatDenseRestoreArray(ds->omat[mT],&T)); |
235 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
187 | PetscCall(MatDenseRestoreArray(ds->omat[mQ],&Q)); |
236 |
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.
|
48 | PetscFunctionReturn(PETSC_SUCCESS); |
237 | } | ||
238 |