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