| 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 | struct HRtr | ||
| 15 | { | ||
| 16 | PetscScalar *data; | ||
| 17 | PetscInt m; | ||
| 18 | PetscInt idx[2]; | ||
| 19 | PetscInt n[2]; | ||
| 20 | PetscScalar tau[2]; | ||
| 21 | PetscReal alpha; | ||
| 22 | PetscReal cs; | ||
| 23 | PetscReal sn; | ||
| 24 | PetscInt type; | ||
| 25 | }; | ||
| 26 | |||
| 27 | /* | ||
| 28 | Generates a hyperbolic rotation | ||
| 29 | if x1*x1 - x2*x2 != 0 | ||
| 30 | r = sqrt(|x1*x1 - x2*x2|) | ||
| 31 | c = x1/r s = x2/r | ||
| 32 | |||
| 33 | | c -s||x1| |d*r| | ||
| 34 | |-s c||x2| = | 0 | | ||
| 35 | where d = 1 for type==1 and -1 for type==2 | ||
| 36 | Returns the condition number of the reduction | ||
| 37 | */ | ||
| 38 | 219908 | static PetscErrorCode HRGen(PetscReal x1,PetscReal x2,PetscInt *type,PetscReal *c,PetscReal *s,PetscReal *r,PetscReal *cond) | |
| 39 | { | ||
| 40 | 219908 | PetscReal t,n2,xa,xb; | |
| 41 | 219908 | PetscInt type_; | |
| 42 | |||
| 43 |
1/2✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
|
219908 | PetscFunctionBegin; |
| 44 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
|
219908 | if (x2==0.0) { |
| 45 | ✗ | *r = PetscAbsReal(x1); *c = (x1>=0.0)?1.0:-1.0; *s = 0.0; | |
| 46 | ✗ | if (type) *type = 1; | |
| 47 | ✗ | PetscFunctionReturn(PETSC_SUCCESS); | |
| 48 | } | ||
| 49 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
|
219908 | if (PetscAbsReal(x1) == PetscAbsReal(x2)) { |
| 50 | /* hyperbolic rotation doesn't exist */ | ||
| 51 | ✗ | *c = *s = *r = 0.0; | |
| 52 | ✗ | if (type) *type = 0; | |
| 53 | ✗ | *cond = PETSC_MAX_REAL; | |
| 54 | ✗ | PetscFunctionReturn(PETSC_SUCCESS); | |
| 55 | } | ||
| 56 | |||
| 57 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
219908 | if (PetscAbsReal(x1)>PetscAbsReal(x2)) { |
| 58 | xa = x1; xb = x2; type_ = 1; | ||
| 59 | } else { | ||
| 60 | 39317 | xa = x2; xb = x1; type_ = 2; | |
| 61 | } | ||
| 62 | 219908 | t = xb/xa; | |
| 63 | 219908 | n2 = PetscAbsReal(1 - t*t); | |
| 64 | 219908 | *r = PetscSqrtReal(n2)*PetscAbsReal(xa); | |
| 65 | 219908 | *c = x1/(*r); | |
| 66 | 219908 | *s = x2/(*r); | |
| 67 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
219908 | if (type_ == 2) *r *= -1; |
| 68 |
1/2✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
|
219908 | if (type) *type = type_; |
| 69 |
1/2✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
|
219908 | if (cond) *cond = (PetscAbsReal(*c) + PetscAbsReal(*s))/PetscAbsReal(PetscAbsReal(*c) - PetscAbsReal(*s)); |
| 70 |
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.
|
45433 | PetscFunctionReturn(PETSC_SUCCESS); |
| 71 | } | ||
| 72 | |||
| 73 | /* | ||
| 74 | |c s| | ||
| 75 | Applies an hyperbolic rotator |s c| | ||
| 76 | |c s| | ||
| 77 | [x1 x2]|s c| | ||
| 78 | */ | ||
| 79 | 532453 | static PetscErrorCode HRApply(PetscInt n,PetscScalar *x1,PetscInt inc1,PetscScalar *x2,PetscInt inc2,PetscReal c,PetscReal s) | |
| 80 | { | ||
| 81 | 532453 | PetscInt i; | |
| 82 | 532453 | PetscReal t; | |
| 83 | 532453 | PetscScalar tmp; | |
| 84 | |||
| 85 |
1/2✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
|
532453 | PetscFunctionBegin; |
| 86 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
532453 | if (PetscAbsReal(c)>PetscAbsReal(s)) { /* Type I */ |
| 87 | 455680 | t = s/c; | |
| 88 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
14047996 | for (i=0;i<n;i++) { |
| 89 | 13592316 | x1[i*inc1] = c*x1[i*inc1] + s*x2[i*inc2]; | |
| 90 | 13592316 | x2[i*inc2] = t*x1[i*inc1] + x2[i*inc2]/c; | |
| 91 | } | ||
| 92 | } else { /* Type II */ | ||
| 93 | 76773 | t = c/s; | |
| 94 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
1695381 | for (i=0;i<n;i++) { |
| 95 | 1618608 | tmp = x1[i*inc1]; | |
| 96 | 1618608 | x1[i*inc1] = c*x1[i*inc1] + s*x2[i*inc2]; | |
| 97 | 1618608 | x2[i*inc2] = t*x1[i*inc1] + tmp/s; | |
| 98 | } | ||
| 99 | } | ||
| 100 |
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.
|
532453 | PetscFunctionReturn(PETSC_SUCCESS); |
| 101 | } | ||
| 102 | |||
| 103 | /* | ||
| 104 | Reduction to tridiagonal-diagonal form (see F. Tisseur, SIMAX 26(1), 2004). | ||
| 105 | |||
| 106 | Input: | ||
| 107 | A symmetric (only lower triangular part is referred) | ||
| 108 | s vector +1 and -1 (signature matrix) | ||
| 109 | Output: | ||
| 110 | d,e | ||
| 111 | s | ||
| 112 | Q s-orthogonal matrix with Q^T*A*Q = T (symmetric tridiagonal matrix) | ||
| 113 | */ | ||
| 114 | 10001 | static PetscErrorCode TridiagDiag_HHR(PetscInt n,PetscScalar *A,PetscInt lda,PetscReal *s,PetscScalar* Q,PetscInt ldq,PetscBool flip,PetscReal *d,PetscReal *e,PetscInt *perm_,PetscScalar *work,PetscReal *rwork,PetscBLASInt *iwork) | |
| 115 | { | ||
| 116 | 10001 | PetscInt i,j,k,*ii,*jj,i0=0,ik=0,type; | |
| 117 | 10001 | PetscInt nwu=0; | |
| 118 | 10001 | PetscReal *ss,cond=1.0,cs,sn,r; | |
| 119 | 10001 | PetscScalar tau,t,*AA; | |
| 120 | 10001 | PetscBLASInt n0,n1,ni,inc=1,m,n_,lda_,ldq_,*perm,tmp; | |
| 121 | 10001 | PetscBool breakdown = PETSC_TRUE; | |
| 122 | |||
| 123 |
1/2✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
|
10001 | PetscFunctionBegin; |
| 124 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
|
10001 | if (n<3) { |
| 125 | ✗ | if (n==1) Q[0]=1; | |
| 126 | ✗ | if (n==2) { | |
| 127 | ✗ | Q[0] = Q[1+ldq] = 1; | |
| 128 | ✗ | Q[1] = Q[ldq] = 0; | |
| 129 | } | ||
| 130 | ✗ | PetscFunctionReturn(PETSC_SUCCESS); | |
| 131 | } | ||
| 132 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10001 | PetscCall(PetscBLASIntCast(lda,&lda_)); |
| 133 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10001 | PetscCall(PetscBLASIntCast(n,&n_)); |
| 134 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10001 | PetscCall(PetscBLASIntCast(ldq,&ldq_)); |
| 135 | 10001 | ss = rwork; | |
| 136 | 10001 | perm = iwork; | |
| 137 | 10001 | AA = work; | |
| 138 |
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.
|
229968 | for (i=0;i<n;i++) PetscCall(PetscArraycpy(AA+i*n,A+i*lda,n)); |
| 139 | 10001 | nwu += n*n; | |
| 140 | 10001 | k=0; | |
| 141 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
20002 | while (breakdown && k<n) { |
| 142 | 10001 | breakdown = PETSC_FALSE; | |
| 143 | /* Classify (and flip) A and s according to sign */ | ||
| 144 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
10001 | if (flip) { |
| 145 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
229528 | for (i=0;i<n;i++) { |
| 146 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
219567 | PetscCall(PetscBLASIntCast(n-1-perm_[i],&perm[i])); |
| 147 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
219567 | if (perm[i]==0) i0 = i; |
| 148 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
219567 | if (perm[i]==k) ik = i; |
| 149 | } | ||
| 150 | } else { | ||
| 151 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
440 | for (i=0;i<n;i++) { |
| 152 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
400 | PetscCall(PetscBLASIntCast(perm_[i],&perm[i])); |
| 153 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
400 | if (perm[i]==0) i0 = i; |
| 154 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
400 | if (perm[i]==k) ik = i; |
| 155 | } | ||
| 156 | } | ||
| 157 | 10001 | perm[ik] = 0; | |
| 158 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10001 | PetscCall(PetscBLASIntCast(k,&perm[i0])); |
| 159 | 10001 | i=1; | |
| 160 |
4/4✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
✓ Branch 2 taken 10 times.
✓ Branch 3 taken 10 times.
|
198755 | while (i<n-1 && s[perm[i-1]]==s[perm[0]]) { |
| 161 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
188754 | if (s[perm[i]]!=s[perm[0]]) { |
| 162 | 162033 | j=i+1; | |
| 163 |
4/4✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
✓ Branch 2 taken 10 times.
✓ Branch 3 taken 10 times.
|
177855 | while (j<n-1 && s[perm[j]]!=s[perm[0]])j++; |
| 164 | 162033 | tmp = perm[i]; perm[i] = perm[j]; perm[j] = tmp; | |
| 165 | } | ||
| 166 | 188754 | i++; | |
| 167 | } | ||
| 168 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
229968 | for (i=0;i<n;i++) { |
| 169 | 219967 | ss[i] = s[perm[i]]; | |
| 170 | } | ||
| 171 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
10001 | if (flip) { |
| 172 | ii = &j; | ||
| 173 | jj = &i; | ||
| 174 | } else { | ||
| 175 | 40 | ii = &i; | |
| 176 | 40 | jj = &j; | |
| 177 | } | ||
| 178 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
229968 | for (i=0;i<n;i++) |
| 179 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
5840550 | for (j=0;j<n;j++) |
| 180 | 5620583 | A[i+j*lda] = AA[perm[*ii]+perm[*jj]*n]; | |
| 181 | /* Initialize Q */ | ||
| 182 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
229968 | for (i=0;i<n;i++) { |
| 183 |
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.
|
219967 | PetscCall(PetscArrayzero(Q+i*ldq,n)); |
| 184 | 219967 | Q[perm[i]+i*ldq] = 1.0; | |
| 185 | } | ||
| 186 |
4/4✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
✓ Branch 2 taken 10 times.
✓ Branch 3 taken 10 times.
|
199038 | for (ni=1;ni<n && ss[ni]==ss[0]; ni++); |
| 187 | 10001 | n0 = ni-1; | |
| 188 | 10001 | n1 = n_-ni; | |
| 189 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
209966 | for (j=0;j<n-2;j++) { |
| 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.
|
199965 | PetscCall(PetscBLASIntCast(n-j-1,&m)); |
| 191 | /* Forming and applying reflectors */ | ||
| 192 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
199965 | if (n0 > 1) { |
| 193 |
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.
|
187336 | PetscCallBLAS("LAPACKlarfg",LAPACKlarfg_(&n0,A+ni-n0+j*lda,A+ni-n0+j*lda+1,&inc,&tau)); |
| 194 | /* Apply reflector */ | ||
| 195 |
1/2✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
|
187336 | if (PetscAbsScalar(tau) != 0.0) { |
| 196 | 187336 | t=*(A+ni-n0+j*lda); *(A+ni-n0+j*lda)=1.0; | |
| 197 |
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.
|
187336 | PetscCallBLAS("LAPACKlarf",LAPACKlarf_("R",&m,&n0,A+ni-n0+j*lda,&inc,&tau,A+j+1+(j+1)*lda,&lda_,work+nwu)); |
| 198 |
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.
|
187336 | PetscCallBLAS("LAPACKlarf",LAPACKlarf_("L",&n0,&m,A+ni-n0+j*lda,&inc,&tau,A+j+1+(j+1)*lda,&lda_,work+nwu)); |
| 199 | /* Update Q */ | ||
| 200 |
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.
|
187336 | PetscCallBLAS("LAPACKlarf",LAPACKlarf_("R",&n_,&n0,A+ni-n0+j*lda,&inc,&tau,Q+(j+1)*ldq,&ldq_,work+nwu)); |
| 201 | 187336 | *(A+ni-n0+j*lda) = t; | |
| 202 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
2490743 | for (i=1;i<n0;i++) { |
| 203 | 2303407 | *(A+ni-n0+j*lda+i) = 0.0; *(A+j+(ni-n0+i)*lda) = 0.0; | |
| 204 | } | ||
| 205 | 187336 | *(A+j+(ni-n0)*lda) = *(A+ni-n0+j*lda); | |
| 206 | } | ||
| 207 | } | ||
| 208 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
199965 | if (n1 > 1) { |
| 209 |
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.
|
14014 | PetscCallBLAS("LAPACKlarfg",LAPACKlarfg_(&n1,A+n-n1+j*lda,A+n-n1+j*lda+1,&inc,&tau)); |
| 210 | /* Apply reflector */ | ||
| 211 |
1/2✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
|
14014 | if (PetscAbsScalar(tau) != 0.0) { |
| 212 | 14014 | t=*(A+n-n1+j*lda); *(A+n-n1+j*lda)=1.0; | |
| 213 |
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.
|
14014 | PetscCallBLAS("LAPACKlarf",LAPACKlarf_("R",&m,&n1,A+n-n1+j*lda,&inc,&tau,A+j+1+(n-n1)*lda,&lda_,work+nwu)); |
| 214 |
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.
|
14014 | PetscCallBLAS("LAPACKlarf",LAPACKlarf_("L",&n1,&m,A+n-n1+j*lda,&inc,&tau,A+n-n1+(j+1)*lda,&lda_,work+nwu)); |
| 215 | /* Update Q */ | ||
| 216 |
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.
|
14014 | PetscCallBLAS("LAPACKlarf",LAPACKlarf_("R",&n_,&n1,A+n-n1+j*lda,&inc,&tau,Q+(n-n1)*ldq,&ldq_,work+nwu)); |
| 217 | 14014 | *(A+n-n1+j*lda) = t; | |
| 218 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
96072 | for (i=1;i<n1;i++) { |
| 219 | 82058 | *(A+n-n1+i+j*lda) = 0.0; *(A+j+(n-n1+i)*lda) = 0.0; | |
| 220 | } | ||
| 221 | 14014 | *(A+j+(n-n1)*lda) = *(A+n-n1+j*lda); | |
| 222 | } | ||
| 223 | } | ||
| 224 | /* Hyperbolic rotation */ | ||
| 225 |
4/4✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
✓ Branch 2 taken 10 times.
✓ Branch 3 taken 10 times.
|
199965 | if (n0 > 0 && n1 > 0) { |
| 226 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
104877 | PetscCall(HRGen(PetscRealPart(A[ni-n0+j*lda]),PetscRealPart(A[n-n1+j*lda]),&type,&cs,&sn,&r,&cond)); |
| 227 | /* Check condition number */ | ||
| 228 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
|
104877 | if (cond > 1.0/(10*PETSC_SQRT_MACHINE_EPSILON)) { |
| 229 | ✗ | breakdown = PETSC_TRUE; | |
| 230 | ✗ | k++; | |
| 231 | ✗ | PetscCheck(k<n && !flip,PETSC_COMM_SELF,PETSC_ERR_SUP,"Breakdown in construction of hyperbolic transformation"); | |
| 232 | break; | ||
| 233 | } | ||
| 234 | 104877 | A[ni-n0+j*lda] = r; A[n-n1+j*lda] = 0.0; | |
| 235 | 104877 | A[j+(ni-n0)*lda] = r; A[j+(n-n1)*lda] = 0.0; | |
| 236 | /* Apply to A */ | ||
| 237 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
104877 | PetscCall(HRApply(m,A+j+1+(ni-n0)*lda,1,A+j+1+(n-n1)*lda,1,cs,-sn)); |
| 238 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
104877 | PetscCall(HRApply(m,A+ni-n0+(j+1)*lda,lda,A+n-n1+(j+1)*lda,lda,cs,-sn)); |
| 239 | |||
| 240 | /* Update Q */ | ||
| 241 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
104877 | PetscCall(HRApply(n,Q+(ni-n0)*ldq,1,Q+(n-n1)*ldq,1,cs,-sn)); |
| 242 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
104877 | if (type==2) { |
| 243 | 8472 | ss[ni-n0] = -ss[ni-n0]; ss[n-n1] = -ss[n-n1]; | |
| 244 | 8472 | n0++;ni++;n1--; | |
| 245 | } | ||
| 246 | } | ||
| 247 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
199965 | if (n0>0) n0--; |
| 248 | 11164 | else n1--; | |
| 249 | } | ||
| 250 | } | ||
| 251 | |||
| 252 | /* flip matrices */ | ||
| 253 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
10001 | if (flip) { |
| 254 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
219567 | for (i=0;i<n-1;i++) { |
| 255 | 209606 | d[i] = PetscRealPart(A[n-i-1+(n-i-1)*lda]); | |
| 256 | 209606 | e[i] = PetscRealPart(A[n-i-1+(n-i-2)*lda]); | |
| 257 | 209606 | s[i] = ss[n-i-1]; | |
| 258 | } | ||
| 259 | 9961 | s[n-1] = ss[0]; | |
| 260 | 9961 | d[n-1] = PetscRealPart(A[0]); | |
| 261 |
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.
|
229528 | for (i=0;i<n;i++) PetscCall(PetscArraycpy(work+i*n,Q+i*ldq,n)); |
| 262 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
229528 | for (i=0;i<n;i++) |
| 263 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
5836150 | for (j=0;j<n;j++) |
| 264 | 5616583 | Q[i+j*ldq] = work[i+(n-j-1)*n]; | |
| 265 | } else { | ||
| 266 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
400 | for (i=0;i<n-1;i++) { |
| 267 | 360 | d[i] = PetscRealPart(A[i+i*lda]); | |
| 268 | 360 | e[i] = PetscRealPart(A[i+1+i*lda]); | |
| 269 | 360 | s[i] = ss[i]; | |
| 270 | } | ||
| 271 | 40 | s[n-1] = ss[n-1]; | |
| 272 | 40 | d[n-1] = PetscRealPart(A[n-1 + (n-1)*lda]); | |
| 273 | } | ||
| 274 |
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.
|
1983 | PetscFunctionReturn(PETSC_SUCCESS); |
| 275 | } | ||
| 276 | |||
| 277 | 335846 | static PetscErrorCode MadeHRtr(PetscInt sz,PetscInt n,PetscInt idx0,PetscInt n0,PetscInt idx1,PetscInt n1,struct HRtr *tr1,struct HRtr *tr2,PetscReal *ncond,PetscScalar *work) | |
| 278 | { | ||
| 279 | 335846 | PetscScalar *x,*y; | |
| 280 | 335846 | PetscReal ncond2=1.0; | |
| 281 | 335846 | PetscBLASInt n0_,n1_,inc=1; | |
| 282 | |||
| 283 |
1/2✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
|
335846 | PetscFunctionBegin; |
| 284 | /* Hyperbolic transformation to make zeros in x */ | ||
| 285 | 335846 | x = tr1->data; | |
| 286 | 335846 | tr1->n[0] = n0; | |
| 287 | 335846 | tr1->n[1] = n1; | |
| 288 | 335846 | tr1->idx[0] = idx0; | |
| 289 | 335846 | tr1->idx[1] = idx1; | |
| 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.
|
335846 | PetscCall(PetscBLASIntCast(tr1->n[0],&n0_)); |
| 291 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
335846 | PetscCall(PetscBLASIntCast(tr1->n[1],&n1_)); |
| 292 |
12/22✓ 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 taken 2 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 2 times.
✓ 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 taken 2 times.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✓ Branch 19 taken 2 times.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
|
335846 | if (tr1->n[0] > 1) PetscCallBLAS("LAPACKlarfg",LAPACKlarfg_(&n0_,x+tr1->idx[0],x+tr1->idx[0]+1,&inc,tr1->tau)); |
| 293 |
12/22✓ 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 taken 2 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 2 times.
✓ 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 taken 2 times.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✓ Branch 19 taken 2 times.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
|
335846 | if (tr1->n[1]> 1) PetscCallBLAS("LAPACKlarfg",LAPACKlarfg_(&n1_,x+tr1->idx[1],x+tr1->idx[1]+1,&inc,tr1->tau+1)); |
| 294 |
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.
|
335846 | if (tr1->idx[0]<tr1->idx[1]) PetscCall(HRGen(PetscRealPart(x[tr1->idx[0]]),PetscRealPart(x[tr1->idx[1]]),&tr1->type,&tr1->cs,&tr1->sn,&tr1->alpha,ncond)); |
| 295 | else { | ||
| 296 | 235037 | tr1->alpha = PetscRealPart(x[tr1->idx[0]]); | |
| 297 | 235037 | *ncond = 1.0; | |
| 298 | } | ||
| 299 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
335846 | if (sz==2) { |
| 300 | 21826 | y = tr2->data; | |
| 301 | /* Apply first transformation to second column */ | ||
| 302 |
3/4✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
✓ Branch 2 taken 10 times.
✗ Branch 3 not taken.
|
21826 | if (tr1->n[0] > 1 && PetscAbsScalar(tr1->tau[0])!=0.0) { |
| 303 | 20728 | x[tr1->idx[0]] = 1.0; | |
| 304 |
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.
|
20728 | PetscCallBLAS("LAPACKlarf",LAPACKlarf_("L",&n0_,&inc,x+tr1->idx[0],&inc,tr1->tau,y+tr1->idx[0],&n0_,work)); |
| 305 | } | ||
| 306 |
3/4✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
✓ Branch 2 taken 10 times.
✗ Branch 3 not taken.
|
21826 | if (tr1->n[1] > 1 && PetscAbsScalar(tr1->tau[1])!=0.0) { |
| 307 | 7716 | x[tr1->idx[1]] = 1.0; | |
| 308 |
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.
|
7716 | PetscCallBLAS("LAPACKlarf",LAPACKlarf_("L",&n1_,&inc,x+tr1->idx[1],&inc,tr1->tau+1,y+tr1->idx[1],&n1_,work)); |
| 309 | } | ||
| 310 |
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.
|
21826 | if (tr1->idx[0]<tr1->idx[1]) PetscCall(HRApply(1,y+tr1->idx[0],1,y+tr1->idx[1],1,tr1->cs,-tr1->sn)); |
| 311 | 21826 | tr2->n[0] = tr1->n[0]; | |
| 312 | 21826 | tr2->n[1] = tr1->n[1]; | |
| 313 | 21826 | tr2->idx[0] = tr1->idx[0]; | |
| 314 | 21826 | tr2->idx[1] = tr1->idx[1]; | |
| 315 |
3/4✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 10 times.
✓ Branch 3 taken 10 times.
|
21826 | if (tr1->idx[0]<tr1->idx[1] && tr1->type==2) { |
| 316 | 10913 | tr2->idx[1]++; tr2->n[1]--; tr2->n[0]++; | |
| 317 | } | ||
| 318 |
1/2✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
|
21826 | if (tr2->n[0]>0) { |
| 319 | 21826 | tr2->n[0]--; tr2->idx[0]++; | |
| 320 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
21826 | if (tr2->n[1]==0) tr2->idx[1] = tr2->idx[0]; |
| 321 | } else { | ||
| 322 | ✗ | tr2->n[1]--; tr2->idx[1]++; tr2->idx[0] = tr2->idx[1]; | |
| 323 | } | ||
| 324 | /* Hyperbolic transformation to make zeros in y */ | ||
| 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.
|
21826 | PetscCall(PetscBLASIntCast(tr2->n[0],&n0_)); |
| 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.
|
21826 | PetscCall(PetscBLASIntCast(tr2->n[1],&n1_)); |
| 327 |
12/22✓ 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 taken 2 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 2 times.
✓ 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 taken 2 times.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✓ Branch 19 taken 2 times.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
|
21826 | if (tr2->n[0] > 1) PetscCallBLAS("LAPACKlarfg",LAPACKlarfg_(&n0_,y+tr2->idx[0],y+tr2->idx[0]+1,&inc,tr2->tau)); |
| 328 |
12/22✓ 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 taken 2 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 2 times.
✓ 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 taken 2 times.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✓ Branch 19 taken 2 times.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
|
21826 | if (tr2->n[1]> 1) PetscCallBLAS("LAPACKlarfg",LAPACKlarfg_(&n1_,y+tr2->idx[1],y+tr2->idx[1]+1,&inc,tr2->tau+1)); |
| 329 |
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.
|
21826 | if (tr2->idx[0]<tr2->idx[1]) PetscCall(HRGen(PetscRealPart(y[tr2->idx[0]]),PetscRealPart(y[tr2->idx[1]]),&tr2->type,&tr2->cs,&tr2->sn,&tr2->alpha,&ncond2)); |
| 330 | else { | ||
| 331 | 7604 | tr2->alpha = PetscRealPart(y[tr2->idx[0]]); | |
| 332 | 7604 | ncond2 = 1.0; | |
| 333 | } | ||
| 334 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
21826 | if (ncond2>*ncond) *ncond = ncond2; |
| 335 | } | ||
| 336 |
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.
|
67713 | PetscFunctionReturn(PETSC_SUCCESS); |
| 337 | } | ||
| 338 | |||
| 339 | /* | ||
| 340 | Auxiliary function to try perform one iteration of hr routine, | ||
| 341 | checking condition number. If it is < tolD, apply the | ||
| 342 | transformation to H and R, if not, ok=false and it do nothing | ||
| 343 | tolE, tolerance to exchange complex pairs to improve conditioning | ||
| 344 | */ | ||
| 345 | 324933 | static PetscErrorCode TryHRIt(PetscInt n,PetscInt j,PetscInt sz,PetscScalar *H,PetscInt ldh,PetscScalar *R,PetscInt ldr,PetscReal *s,PetscBool *exg,PetscBool *ok,PetscInt *n0,PetscInt *n1,PetscInt *idx0,PetscInt *idx1,PetscReal *cond,PetscScalar *work) | |
| 346 | { | ||
| 347 | 324933 | struct HRtr *tr1,*tr2,tr1_t,tr2_t,tr1_te,tr2_te; | |
| 348 | 324933 | PetscScalar *x,*y; | |
| 349 | 324933 | PetscReal ncond=0,ncond_e; | |
| 350 | 324933 | PetscInt nwu=0,i,d=1; | |
| 351 | 324933 | PetscBLASInt n0_,n1_,inc=1,mh,mr,n_,ldr_,ldh_; | |
| 352 | 324933 | PetscReal tolD = 1e+5; | |
| 353 | |||
| 354 |
1/2✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
|
324933 | PetscFunctionBegin; |
| 355 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
|
324933 | if (cond) *cond = 1.0; |
| 356 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
324933 | PetscCall(PetscBLASIntCast(n,&n_)); |
| 357 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
324933 | PetscCall(PetscBLASIntCast(ldr,&ldr_)); |
| 358 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
324933 | PetscCall(PetscBLASIntCast(ldh,&ldh_)); |
| 359 | 324933 | x = work+nwu; | |
| 360 | 324933 | nwu += n; | |
| 361 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
324933 | PetscCall(PetscArraycpy(x,R+j*ldr,n)); |
| 362 | 324933 | *exg = PETSC_FALSE; | |
| 363 | 324933 | *ok = PETSC_TRUE; | |
| 364 | 324933 | tr1_t.data = x; | |
| 365 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
324933 | if (sz==1) { |
| 366 | /* Hyperbolic transformation to make zeros in x */ | ||
| 367 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
314020 | PetscCall(MadeHRtr(sz,n,*idx0,*n0,*idx1,*n1,&tr1_t,NULL,&ncond,work+nwu)); |
| 368 | /* Check condition number to single column*/ | ||
| 369 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
|
314020 | if (ncond>tolD) *ok = PETSC_FALSE; |
| 370 | tr1 = &tr1_t; | ||
| 371 | tr2 = &tr2_t; | ||
| 372 | } else { | ||
| 373 | 10913 | y = work+nwu; | |
| 374 | 10913 | nwu += n; | |
| 375 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10913 | PetscCall(PetscArraycpy(y,R+(j+1)*ldr,n)); |
| 376 | 10913 | tr2_t.data = y; | |
| 377 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10913 | PetscCall(MadeHRtr(sz,n,*idx0,*n0,*idx1,*n1,&tr1_t,&tr2_t,&ncond,work+nwu)); |
| 378 | /* Computing hyperbolic transformations also for exchanged vectors */ | ||
| 379 | 10913 | tr1_te.data = work+nwu; | |
| 380 | 10913 | nwu += n; | |
| 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.
|
10913 | PetscCall(PetscArraycpy(tr1_te.data,R+(j+1)*ldr,n)); |
| 382 | 10913 | tr2_te.data = work+nwu; | |
| 383 | 10913 | nwu += n; | |
| 384 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10913 | PetscCall(PetscArraycpy(tr2_te.data,R+j*ldr,n)); |
| 385 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10913 | PetscCall(MadeHRtr(sz,n,*idx0,*n0,*idx1,*n1,&tr1_te,&tr2_te,&ncond_e,work+nwu)); |
| 386 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
10913 | if (ncond > d*ncond_e) { |
| 387 | 6204 | *exg = PETSC_TRUE; | |
| 388 | 6204 | tr1 = &tr1_te; | |
| 389 | 6204 | tr2 = &tr2_te; | |
| 390 | 6204 | ncond = ncond_e; | |
| 391 | } else { | ||
| 392 | tr1 = &tr1_t; | ||
| 393 | tr2 = &tr2_t; | ||
| 394 | } | ||
| 395 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
|
10913 | if (ncond>tolD) *ok = PETSC_FALSE; |
| 396 | } | ||
| 397 |
1/2✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
|
324933 | if (*ok) { |
| 398 | /* Everything is OK, apply transformations to R and H */ | ||
| 399 | /* First column */ | ||
| 400 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
324933 | if (cond && *cond<ncond) *cond = ncond; |
| 401 | 324933 | x = tr1->data; | |
| 402 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
324933 | PetscCall(PetscBLASIntCast(tr1->n[0],&n0_)); |
| 403 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
324933 | PetscCall(PetscBLASIntCast(tr1->n[1],&n1_)); |
| 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.
|
324933 | PetscCall(PetscBLASIntCast(n-j-sz,&mr)); |
| 405 |
3/4✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
✓ Branch 2 taken 10 times.
✗ Branch 3 not taken.
|
324933 | if (tr1->n[0] > 1 && PetscAbsScalar(tr1->tau[0])!=0.0) { |
| 406 | 302394 | x[tr1->idx[0]] = 1.0; | |
| 407 |
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.
|
302394 | PetscCallBLAS("LAPACKlarf",LAPACKlarf_("L",&n0_,&mr,x+tr1->idx[0],&inc,tr1->tau,R+(j+sz)*ldr+tr1->idx[0],&ldr_,work+nwu)); |
| 408 |
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.
|
302394 | PetscCallBLAS("LAPACKlarf",LAPACKlarf_("R",&n_,&n0_,x+tr1->idx[0],&inc,tr1->tau,H+(tr1->idx[0])*ldh,&ldh_,work+nwu)); |
| 409 | } | ||
| 410 |
3/4✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
✓ Branch 2 taken 10 times.
✗ Branch 3 not taken.
|
324933 | if (tr1->n[1] > 1 && PetscAbsScalar(tr1->tau[1])!=0.0) { |
| 411 | 50965 | x[tr1->idx[1]] = 1.0; | |
| 412 |
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.
|
50965 | PetscCallBLAS("LAPACKlarf",LAPACKlarf_("L",&n1_,&mr,x+tr1->idx[1],&inc,tr1->tau+1,R+(j+sz)*ldr+tr1->idx[1],&ldr_,work+nwu)); |
| 413 |
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.
|
50965 | PetscCallBLAS("LAPACKlarf",LAPACKlarf_("R",&n_,&n1_,x+tr1->idx[1],&inc,tr1->tau+1,H+(tr1->idx[1])*ldh,&ldh_,work+nwu)); |
| 414 | } | ||
| 415 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
324933 | if (tr1->idx[0]<tr1->idx[1]) { |
| 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.
|
89896 | PetscCall(HRApply(mr,R+(j+sz)*ldr+tr1->idx[0],ldr,R+(j+sz)*ldr+tr1->idx[1],ldr,tr1->cs,-tr1->sn)); |
| 417 |
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.
|
89896 | if (tr1->type==1) PetscCall(HRApply(n,H+(tr1->idx[0])*ldh,1,H+(tr1->idx[1])*ldh,1,tr1->cs,tr1->sn)); |
| 418 | else { | ||
| 419 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
13879 | PetscCall(HRApply(n,H+(tr1->idx[0])*ldh,1,H+(tr1->idx[1])*ldh,1,-tr1->cs,-tr1->sn)); |
| 420 | 13879 | s[tr1->idx[0]] = -s[tr1->idx[0]]; | |
| 421 | 13879 | s[tr1->idx[1]] = -s[tr1->idx[1]]; | |
| 422 | } | ||
| 423 | } | ||
| 424 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
6703128 | for (i=0;i<tr1->idx[0];i++) *(R+j*ldr+i) = x[i]; |
| 425 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
6567374 | for (i=tr1->idx[0]+1;i<n;i++) *(R+j*ldr+i) = 0.0; |
| 426 | 324933 | *(R+j*ldr+tr1->idx[0]) = tr1->alpha; | |
| 427 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
324933 | if (sz==2) { |
| 428 | 10913 | y = tr2->data; | |
| 429 | /* Second column */ | ||
| 430 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10913 | PetscCall(PetscBLASIntCast(tr2->n[0],&n0_)); |
| 431 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10913 | PetscCall(PetscBLASIntCast(tr2->n[1],&n1_)); |
| 432 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10913 | PetscCall(PetscBLASIntCast(n-j-sz,&mr)); |
| 433 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10913 | PetscCall(PetscBLASIntCast(n-tr2->idx[0],&mh)); |
| 434 |
3/4✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
✓ Branch 2 taken 10 times.
✗ Branch 3 not taken.
|
10913 | if (tr2->n[0] > 1 && PetscAbsScalar(tr2->tau[0])!=0.0) { |
| 435 | 10101 | y[tr2->idx[0]] = 1.0; | |
| 436 |
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.
|
10101 | PetscCallBLAS("LAPACKlarf",LAPACKlarf_("L",&n0_,&mr,y+tr2->idx[0],&inc,tr2->tau,R+(j+2)*ldr+tr2->idx[0],&ldr_,work+nwu)); |
| 437 |
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.
|
10101 | PetscCallBLAS("LAPACKlarf",LAPACKlarf_("R",&n_,&n0_,y+tr2->idx[0],&inc,tr2->tau,H+(tr2->idx[0])*ldh,&ldh_,work+nwu)); |
| 438 | } | ||
| 439 |
3/4✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
✓ Branch 2 taken 10 times.
✗ Branch 3 not taken.
|
10913 | if (tr2->n[1] > 1 && PetscAbsScalar(tr2->tau[1])!=0.0) { |
| 440 | 3508 | y[tr2->idx[1]] = 1.0; | |
| 441 |
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.
|
3508 | PetscCallBLAS("LAPACKlarf",LAPACKlarf_("L",&n1_,&mr,y+tr2->idx[1],&inc,tr2->tau+1,R+(j+2)*ldr+tr2->idx[1],&ldr_,work+nwu)); |
| 442 |
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.
|
3508 | PetscCallBLAS("LAPACKlarf",LAPACKlarf_("R",&n_,&n1_,y+tr2->idx[1],&inc,tr2->tau+1,H+(tr2->idx[1])*ldh,&ldh_,work+nwu)); |
| 443 | } | ||
| 444 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
10913 | if (tr2->idx[0]<tr2->idx[1]) { |
| 445 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
8102 | PetscCall(HRApply(mr,R+(j+2)*ldr+tr2->idx[0],ldr,R+(j+2)*ldr+tr2->idx[1],ldr,tr2->cs,-tr2->sn)); |
| 446 |
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.
|
8102 | if (tr2->type==1) PetscCall(HRApply(n,H+(tr2->idx[0])*ldh,1,H+(tr2->idx[1])*ldh,1,tr2->cs,tr2->sn)); |
| 447 | else { | ||
| 448 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
6343 | PetscCall(HRApply(n,H+(tr2->idx[0])*ldh,1,H+(tr2->idx[1])*ldh,1,-tr2->cs,-tr2->sn)); |
| 449 | 6343 | s[tr2->idx[0]] = -s[tr2->idx[0]]; | |
| 450 | 6343 | s[tr2->idx[1]] = -s[tr2->idx[1]]; | |
| 451 | } | ||
| 452 | } | ||
| 453 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
122246 | for (i=0;i<tr2->idx[0]-1;i++) *(R+(j+1)*ldr+i) = y[i]; |
| 454 | 10913 | *(R+(j+1)*ldr+tr2->idx[0]-1) = y[tr2->idx[0]-1]; | |
| 455 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
268913 | for (i=tr2->idx[0]+1;i<n;i++) *(R+(j+1)*ldr+i) = 0.0; |
| 456 | 10913 | *(R+(j+1)*ldr+tr2->idx[0]) = tr2->alpha; | |
| 457 | 10913 | *n0 = tr2->n[0]; | |
| 458 | 10913 | *n1 = tr2->n[1]; | |
| 459 | 10913 | *idx0 = tr2->idx[0]; | |
| 460 | 10913 | *idx1 = tr2->idx[1]; | |
| 461 |
4/4✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
✓ Branch 2 taken 10 times.
✓ Branch 3 taken 10 times.
|
10913 | if (tr2->idx[0]<tr2->idx[1] && tr2->type==2) { |
| 462 | 6343 | (*idx1)++; (*n1)--; (*n0)++; | |
| 463 | } | ||
| 464 | } else { | ||
| 465 | 314020 | *n0 = tr1->n[0]; | |
| 466 | 314020 | *n1 = tr1->n[1]; | |
| 467 | 314020 | *idx0 = tr1->idx[0]; | |
| 468 | 314020 | *idx1 = tr1->idx[1]; | |
| 469 |
4/4✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
✓ Branch 2 taken 10 times.
✓ Branch 3 taken 10 times.
|
314020 | if (tr1->idx[0]<tr1->idx[1] && tr1->type==2) { |
| 470 | 9568 | (*idx1)++; (*n1)--; (*n0)++; | |
| 471 | } | ||
| 472 | } | ||
| 473 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
324933 | if (*n0>0) { |
| 474 | 313560 | (*n0)--; (*idx0)++; | |
| 475 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
313560 | if (*n1==0) *idx1 = *idx0; |
| 476 | } else { | ||
| 477 | 11373 | (*n1)--; (*idx1)++; *idx0 = *idx1; | |
| 478 | } | ||
| 479 | } | ||
| 480 |
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.
|
65465 | PetscFunctionReturn(PETSC_SUCCESS); |
| 481 | } | ||
| 482 | |||
| 483 | /* | ||
| 484 | compute V = HR whit H s-orthogonal and R upper triangular | ||
| 485 | (work: work space of minimum size 6*nv) | ||
| 486 | */ | ||
| 487 | 10520 | PetscErrorCode DSPseudoOrthog_HR(PetscInt *nv,PetscScalar *V,PetscInt ldv,PetscReal *s,PetscScalar *R,PetscInt ldr,PetscBLASInt *perm,PetscBLASInt *cmplxEig,PetscBool *breakdown,PetscScalar *work) | |
| 488 | { | ||
| 489 | 10520 | PetscInt i,j,n,n0,n1,np,idx0,idx1,sz=1,k=0,nwu=0; | |
| 490 | 10520 | PetscBLASInt t1,t2; | |
| 491 | 10520 | PetscScalar *col1,*col2; | |
| 492 | 10520 | PetscBool exg=PETSC_FALSE,ok=PETSC_FALSE; | |
| 493 | |||
| 494 |
1/2✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
|
10520 | PetscFunctionBegin; |
| 495 | 10520 | n = *nv; | |
| 496 | 10520 | col1 = work+nwu; | |
| 497 | 10520 | nwu += n; | |
| 498 | 10520 | col2 = work+nwu; | |
| 499 | 10520 | nwu += n; | |
| 500 | /* Sort R and s according to sing(s) */ | ||
| 501 | 10520 | np = 0; | |
| 502 |
4/4✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
✓ Branch 2 taken 10 times.
✓ Branch 3 taken 10 times.
|
346366 | for (i=0;i<n;i++) if (s[i]>0) np++; |
| 503 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
10520 | if (s[0]>0) n1 = np; |
| 504 | 994 | else n1 = n-np; | |
| 505 | 10520 | n0 = 0; | |
| 506 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
346366 | for (i=0;i<n;i++) { |
| 507 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
335846 | if (s[i]==s[0]) { |
| 508 | 304251 | s[n0] = s[0]; | |
| 509 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
304251 | PetscCall(PetscBLASIntCast(i,&perm[n0++])); |
| 510 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
335846 | } else PetscCall(PetscBLASIntCast(i,&perm[n1++])); |
| 511 | } | ||
| 512 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
42115 | for (i=n0;i<n;i++) s[i] = -s[0]; |
| 513 | 10520 | n1 -= n0; | |
| 514 | 10520 | idx0 = 0; | |
| 515 | 10520 | idx1 = n0; | |
| 516 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
10520 | if (idx1==n) idx1=idx0; |
| 517 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
346366 | for (i=0;i<n;i++) { |
| 518 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
13672574 | for (j=0;j<n;j++) R[j*ldr+i] = V[j*ldv+perm[i]]; |
| 519 | } | ||
| 520 | /* Initialize H */ | ||
| 521 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
346366 | for (i=0;i<n;i++) { |
| 522 |
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.
|
335846 | PetscCall(PetscArrayzero(V+i*ldv,n)); |
| 523 | 335846 | V[perm[i]+i*ldv] = 1.0; | |
| 524 | } | ||
| 525 |
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.
|
346366 | for (i=0;i<n;i++) PetscCall(PetscBLASIntCast(i,&perm[i])); |
| 526 | j = 0; | ||
| 527 | 335453 | while (j<n-k) { | |
| 528 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
324933 | if (cmplxEig[j]==0) sz=1; |
| 529 | 10913 | else sz=2; | |
| 530 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
324933 | PetscCall(TryHRIt(n,j,sz,V,ldv,R,ldr,s,&exg,&ok,&n0,&n1,&idx0,&idx1,NULL,work+nwu)); |
| 531 |
1/2✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
|
324933 | if (ok) { |
| 532 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
324933 | if (exg) cmplxEig[j] = -cmplxEig[j]; |
| 533 | 324933 | j = j+sz; | |
| 534 | } else { /* to be discarded */ | ||
| 535 | ✗ | k = k+1; | |
| 536 | ✗ | if (cmplxEig[j]==0) { | |
| 537 | ✗ | if (j<n) { | |
| 538 | ✗ | t1 = perm[j]; | |
| 539 | ✗ | for (i=j;i<n-1;i++) perm[i] = perm[i+1]; | |
| 540 | ✗ | perm[n-1] = t1; | |
| 541 | ✗ | t1 = cmplxEig[j]; | |
| 542 | ✗ | for (i=j;i<n-1;i++) cmplxEig[i] = cmplxEig[i+1]; | |
| 543 | ✗ | cmplxEig[n-1] = t1; | |
| 544 | ✗ | PetscCall(PetscArraycpy(col1,R+j*ldr,n)); | |
| 545 | ✗ | for (i=j;i<n-1;i++) PetscCall(PetscArraycpy(R+i*ldr,R+(i+1)*ldr,n)); | |
| 546 | ✗ | PetscCall(PetscArraycpy(R+(n-1)*ldr,col1,n)); | |
| 547 | } | ||
| 548 | } else { | ||
| 549 | ✗ | k = k+1; | |
| 550 | ✗ | if (j<n-1) { | |
| 551 | ✗ | t1 = perm[j]; t2 = perm[j+1]; | |
| 552 | ✗ | for (i=j;i<n-2;i++) perm[i] = perm[i+2]; | |
| 553 | ✗ | perm[n-2] = t1; perm[n-1] = t2; | |
| 554 | ✗ | t1 = cmplxEig[j]; t2 = cmplxEig[j+1]; | |
| 555 | ✗ | for (i=j;i<n-2;i++) cmplxEig[i] = cmplxEig[i+2]; | |
| 556 | ✗ | cmplxEig[n-2] = t1; cmplxEig[n-1] = t2; | |
| 557 | ✗ | PetscCall(PetscArraycpy(col1,R+j*ldr,n)); | |
| 558 | ✗ | PetscCall(PetscArraycpy(col2,R+(j+1)*ldr,n)); | |
| 559 | ✗ | for (i=j;i<n-2;i++) PetscCall(PetscArraycpy(R+i*ldr,R+(i+2)*ldr,n)); | |
| 560 | ✗ | PetscCall(PetscArraycpy(R+(n-2)*ldr,col1,n)); | |
| 561 |
4/8✗ Branch 0 not taken.
✗ Branch 1 not taken.
✓ Branch 2 taken 8 times.
✓ Branch 3 taken 8 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✓ Branch 6 taken 2 times.
✓ Branch 7 taken 2 times.
|
335453 | PetscCall(PetscArraycpy(R+(n-1)*ldr,col2,n)); |
| 562 | } | ||
| 563 | } | ||
| 564 | } | ||
| 565 | } | ||
| 566 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
|
10520 | if (k!=0) { |
| 567 | ✗ | if (breakdown) *breakdown = PETSC_TRUE; | |
| 568 | ✗ | *nv = n-k; | |
| 569 | } | ||
| 570 |
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.
|
2092 | PetscFunctionReturn(PETSC_SUCCESS); |
| 571 | } | ||
| 572 | |||
| 573 | 10500 | PetscErrorCode DSGHIEPOrthogEigenv(DS ds,DSMatType mat,PetscScalar *wr,PetscScalar *wi,PetscBool accum) | |
| 574 | { | ||
| 575 | 10500 | PetscInt lws,nwus=0,nwui=0,lwi,off,n,nv,ld,i,ldr,l; | |
| 576 | 10500 | const PetscScalar *B,*W; | |
| 577 | 10500 | PetscScalar *Q,*X,*R,*ts,szero=0.0,sone=1.0; | |
| 578 | 10500 | PetscReal *s,vi,vr,tr,*d,*e; | |
| 579 | 10500 | PetscBLASInt ld_,n_,nv_,*perm,*cmplxEig; | |
| 580 | |||
| 581 |
1/2✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
|
10500 | PetscFunctionBegin; |
| 582 | 10500 | l = ds->l; | |
| 583 | 10500 | n = ds->n-l; | |
| 584 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10500 | PetscCall(PetscBLASIntCast(n,&n_)); |
| 585 | 10500 | ld = ds->ld; | |
| 586 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10500 | PetscCall(PetscBLASIntCast(ld,&ld_)); |
| 587 | 10500 | off = l*ld+l; | |
| 588 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10500 | PetscCall(DSGetArrayReal(ds,DS_MAT_D,&s)); |
| 589 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
10500 | if (!ds->compact) { |
| 590 |
4/6✓ Branch 0 taken 2 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)); |
| 591 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
440 | for (i=l;i<ds->n;i++) s[i] = PetscRealPart(B[i*ld+i]); |
| 592 |
4/6✓ Branch 0 taken 2 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(MatDenseRestoreArrayRead(ds->omat[DS_MAT_B],&B)); |
| 593 | } | ||
| 594 | 10500 | lws = n*n+7*n; | |
| 595 | 10500 | lwi = 2*n; | |
| 596 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10500 | PetscCall(DSAllocateWork_Private(ds,lws,0,lwi)); |
| 597 | 10500 | R = ds->work+nwus; | |
| 598 | 10500 | nwus += n*n; | |
| 599 | 10500 | ldr = n; | |
| 600 | 10500 | perm = ds->iwork + nwui; | |
| 601 | 10500 | nwui += n; | |
| 602 | 10500 | cmplxEig = ds->iwork+nwui; | |
| 603 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10500 | PetscCall(MatDenseGetArray(ds->omat[mat],&X)); |
| 604 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
335273 | for (i=0;i<n;i++) { |
| 605 | #if defined(PETSC_USE_COMPLEX) | ||
| 606 | 47006 | vi = PetscImaginaryPart(wr[l+i]); | |
| 607 | #else | ||
| 608 |
1/2✓ Branch 0 taken 5 times.
✗ Branch 1 not taken.
|
277767 | vi = wi?PetscRealPart(wi[l+i]):0.0; |
| 609 | #endif | ||
| 610 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
324773 | if (vi!=0) { |
| 611 | 10913 | cmplxEig[i] = 1; | |
| 612 | 10913 | cmplxEig[i+1] = 2; | |
| 613 | 10913 | i++; | |
| 614 | 313860 | } else cmplxEig[i] = 0; | |
| 615 | } | ||
| 616 | 10500 | nv = n; | |
| 617 | |||
| 618 | /* Perform HR decomposition */ | ||
| 619 | /* Hyperbolic rotators */ | ||
| 620 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10500 | PetscCall(DSPseudoOrthog_HR(&nv,X+off,ld,s+l,R,ldr,perm,cmplxEig,NULL,ds->work+nwus)); |
| 621 | /* Sort wr,wi perm */ | ||
| 622 | 10500 | ts = ds->work+nwus; | |
| 623 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10500 | PetscCall(PetscArraycpy(ts,wr+l,n)); |
| 624 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
346186 | for (i=0;i<n;i++) wr[i+l] = ts[perm[i]]; |
| 625 | #if !defined(PETSC_USE_COMPLEX) | ||
| 626 |
1/2✓ Branch 0 taken 5 times.
✗ Branch 1 not taken.
|
8714 | if (wi) { |
| 627 |
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.
|
8714 | PetscCall(PetscArraycpy(ts,wi+l,n)); |
| 628 |
2/2✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
|
296220 | for (i=0;i<n;i++) wi[i+l] = ts[perm[i]]; |
| 629 | } | ||
| 630 | #endif | ||
| 631 | /* Projected Matrix */ | ||
| 632 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10500 | PetscCall(DSGetArrayReal(ds,DS_MAT_T,&d)); |
| 633 |
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.
|
10500 | PetscCall(PetscArrayzero(d+2*ld,ld)); |
| 634 | 10500 | e = d+ld; | |
| 635 | 10500 | d[l+nv-1] = PetscRealPart(wr[l+nv-1]*s[l+nv-1]); | |
| 636 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
325110 | for (i=0;i<nv-1;i++) { |
| 637 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
314610 | if (cmplxEig[i]==0) { /* Real */ |
| 638 | 303697 | d[l+i] = PetscRealPart(wr[l+i]*s[l+i]); | |
| 639 | 303697 | e[l+i] = 0.0; | |
| 640 | } else { | ||
| 641 | 10913 | vr = PetscRealPart(wr[l+i]); | |
| 642 | #if defined(PETSC_USE_COMPLEX) | ||
| 643 | 1174 | vi = PetscImaginaryPart(wr[l+i]); | |
| 644 | #else | ||
| 645 |
1/2✓ Branch 0 taken 5 times.
✗ Branch 1 not taken.
|
9739 | vi = wi?PetscRealPart(wi[l+i]):0.0; |
| 646 | #endif | ||
| 647 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
10913 | if (cmplxEig[i]==-1) vi = -vi; |
| 648 | 10913 | tr = PetscRealPart((R[i+(i+1)*ldr]/R[i+i*ldr]))*vi; | |
| 649 | 10913 | d[l+i] = (vr-tr)*s[l+i]; | |
| 650 | 10913 | d[l+i+1] = (vr+tr)*s[l+i+1]; | |
| 651 | 10913 | e[l+i] = PetscRealPart(s[l+i]*(R[(i+1)+(i+1)*ldr]/R[i+i*ldr])*vi); | |
| 652 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
10913 | if (i<nv-2) e[l+i+1] = 0.0; |
| 653 | i++; | ||
| 654 | } | ||
| 655 | } | ||
| 656 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10500 | PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&d)); |
| 657 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10500 | PetscCall(DSRestoreArrayReal(ds,DS_MAT_D,&s)); |
| 658 | /* accumulate previous Q */ | ||
| 659 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
10500 | if (accum) { |
| 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.
|
10470 | PetscCall(MatDenseGetArray(ds->omat[DS_MAT_Q],&Q)); |
| 661 |
4/6✓ Branch 0 taken 2 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(nv,&nv_)); |
| 662 |
4/6✓ Branch 0 taken 2 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)); |
| 663 |
4/6✓ Branch 0 taken 2 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(MatCopy(ds->omat[DS_MAT_Q],ds->omat[DS_MAT_W],SAME_NONZERO_PATTERN)); |
| 664 |
4/6✓ Branch 0 taken 2 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_W],&W)); |
| 665 |
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.
|
10470 | PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n_,&nv_,&n_,&sone,W+off,&ld_,X+off,&ld_,&szero,Q+off,&ld_)); |
| 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.
|
10470 | PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_W],&W)); |
| 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.
|
10470 | PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_Q],&Q)); |
| 668 | } | ||
| 669 | 10500 | ds->t = nv+l; | |
| 670 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10500 | PetscCall(MatDenseRestoreArray(ds->omat[mat],&X)); |
| 671 |
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.
|
10500 | if (!ds->compact) PetscCall(DSSwitchFormat_GHIEP(ds,PETSC_FALSE)); |
| 672 |
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.
|
2088 | PetscFunctionReturn(PETSC_SUCCESS); |
| 673 | } | ||
| 674 | |||
| 675 | /* | ||
| 676 | Reduce to tridiagonal-diagonal pair by means of TridiagDiag_HHR. | ||
| 677 | */ | ||
| 678 | 10500 | PetscErrorCode DSIntermediate_GHIEP(DS ds) | |
| 679 | { | ||
| 680 | 10500 | PetscInt i,ld,off; | |
| 681 | 10500 | PetscInt nwall,nwallr,nwalli; | |
| 682 | 10500 | PetscScalar *A,*B,*Q; | |
| 683 | 10500 | PetscReal *d,*e,*s; | |
| 684 | |||
| 685 |
1/2✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
|
10500 | PetscFunctionBegin; |
| 686 | 10500 | ld = ds->ld; | |
| 687 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10500 | PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A)); |
| 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.
|
10500 | PetscCall(MatDenseGetArray(ds->omat[DS_MAT_B],&B)); |
| 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.
|
10500 | PetscCall(MatDenseGetArray(ds->omat[DS_MAT_Q],&Q)); |
| 690 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10500 | PetscCall(DSGetArrayReal(ds,DS_MAT_T,&d)); |
| 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.
|
10500 | PetscCall(DSGetArrayReal(ds,DS_MAT_D,&s)); |
| 692 | 10500 | e = d+ld; | |
| 693 | 10500 | off = ds->l+ds->l*ld; | |
| 694 |
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.
|
10500 | PetscCall(PetscArrayzero(Q,ld*ld)); |
| 695 | 10500 | nwall = ld*ld+ld; | |
| 696 | 10500 | nwallr = ld; | |
| 697 | 10500 | nwalli = ld; | |
| 698 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10500 | PetscCall(DSAllocateWork_Private(ds,nwall,nwallr,nwalli)); |
| 699 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
357842 | for (i=0;i<ds->n;i++) Q[i+i*ld]=1.0; |
| 700 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
346186 | for (i=0;i<ds->n-ds->l;i++) *(ds->perm+i)=i; |
| 701 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
10500 | if (ds->compact) { |
| 702 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
10460 | if (ds->state < DS_STATE_INTERMEDIATE) { |
| 703 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
9961 | PetscCall(DSSwitchFormat_GHIEP(ds,PETSC_FALSE)); |
| 704 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
9961 | PetscCall(TridiagDiag_HHR(ds->k-ds->l+1,A+off,ld,s+ds->l,Q+off,ld,PETSC_TRUE,d+ds->l,e+ds->l,ds->perm,ds->work,ds->rwork,ds->iwork)); |
| 705 | 9961 | ds->k = ds->l; | |
| 706 |
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.
|
9961 | PetscCall(PetscArrayzero(d+2*ld+ds->l,ds->n-ds->l)); |
| 707 | } | ||
| 708 | } else { | ||
| 709 |
1/2✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
|
40 | if (ds->state < DS_STATE_INTERMEDIATE) { |
| 710 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
440 | for (i=0;i<ds->n;i++) s[i] = PetscRealPart(B[i+i*ld]); |
| 711 |
4/6✓ Branch 0 taken 2 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(TridiagDiag_HHR(ds->n-ds->l,A+off,ld,s+ds->l,Q+off,ld,PETSC_FALSE,d+ds->l,e+ds->l,ds->perm,ds->work,ds->rwork,ds->iwork)); |
| 712 |
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.
|
40 | PetscCall(PetscArrayzero(d+2*ld,ds->n)); |
| 713 | 40 | ds->k = ds->l; | |
| 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.
|
40 | PetscCall(DSSwitchFormat_GHIEP(ds,PETSC_FALSE)); |
| 715 | ✗ | } else PetscCall(DSSwitchFormat_GHIEP(ds,PETSC_TRUE)); | |
| 716 | } | ||
| 717 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10500 | PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A)); |
| 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.
|
10500 | PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_B],&B)); |
| 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.
|
10500 | PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_Q],&Q)); |
| 720 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10500 | PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&d)); |
| 721 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10500 | PetscCall(DSRestoreArrayReal(ds,DS_MAT_D,&s)); |
| 722 |
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.
|
2088 | PetscFunctionReturn(PETSC_SUCCESS); |
| 723 | } | ||
| 724 |