| 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/slepccontour.h> | ||
| 12 | #include <slepcblaslapack.h> | ||
| 13 | |||
| 14 | /* | ||
| 15 | SlepcContourDataCreate - Create a contour data structure. | ||
| 16 | |||
| 17 | Input Parameters: | ||
| 18 | n - the number of integration points | ||
| 19 | npart - number of partitions for the subcommunicator | ||
| 20 | parent - parent object | ||
| 21 | */ | ||
| 22 | 447 | PetscErrorCode SlepcContourDataCreate(PetscInt n,PetscInt npart,PetscObject parent,SlepcContourData *contour) | |
| 23 | { | ||
| 24 |
1/2✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
|
447 | PetscFunctionBegin; |
| 25 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
447 | PetscCall(PetscNew(contour)); |
| 26 | 447 | (*contour)->parent = parent; | |
| 27 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
447 | PetscCall(PetscSubcommCreate(PetscObjectComm(parent),&(*contour)->subcomm)); |
| 28 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
447 | PetscCall(PetscSubcommSetNumber((*contour)->subcomm,npart)); |
| 29 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
447 | PetscCall(PetscSubcommSetType((*contour)->subcomm,PETSC_SUBCOMM_INTERLACED)); |
| 30 | 447 | (*contour)->npoints = n / npart; | |
| 31 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 8 times.
|
447 | if (n%npart > (*contour)->subcomm->color) (*contour)->npoints++; |
| 32 |
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.
|
116 | PetscFunctionReturn(PETSC_SUCCESS); |
| 33 | } | ||
| 34 | |||
| 35 | /* | ||
| 36 | SlepcContourDataReset - Resets the KSP objects in a contour data structure, | ||
| 37 | and destroys any objects whose size depends on the problem size. | ||
| 38 | */ | ||
| 39 | 337 | PetscErrorCode SlepcContourDataReset(SlepcContourData contour) | |
| 40 | { | ||
| 41 | 337 | PetscInt i; | |
| 42 | |||
| 43 |
1/2✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
|
337 | PetscFunctionBegin; |
| 44 |
1/2✓ Branch 0 taken 8 times.
✗ Branch 1 not taken.
|
337 | if (contour->ksp) { |
| 45 |
7/8✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 8 times.
✓ Branch 3 taken 6 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
✓ Branch 6 taken 2 times.
✓ Branch 7 taken 2 times.
|
7745 | for (i=0;i<contour->npoints;i++) PetscCall(KSPReset(contour->ksp[i])); |
| 46 | } | ||
| 47 |
2/2✓ Branch 0 taken 8 times.
✓ Branch 1 taken 8 times.
|
337 | if (contour->pA) { |
| 48 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
140 | PetscCall(MatDestroyMatrices(contour->nmat,&contour->pA)); |
| 49 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
140 | PetscCall(MatDestroyMatrices(contour->nmat,&contour->pP)); |
| 50 | 140 | contour->nmat = 0; | |
| 51 | } | ||
| 52 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
337 | PetscCall(VecScatterDestroy(&contour->scatterin)); |
| 53 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
337 | PetscCall(VecDestroy(&contour->xsub)); |
| 54 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
337 | PetscCall(VecDestroy(&contour->xdup)); |
| 55 |
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.
|
88 | PetscFunctionReturn(PETSC_SUCCESS); |
| 56 | } | ||
| 57 | |||
| 58 | /* | ||
| 59 | SlepcContourDataDestroy - Destroys the contour data structure. | ||
| 60 | */ | ||
| 61 | 657 | PetscErrorCode SlepcContourDataDestroy(SlepcContourData *contour) | |
| 62 | { | ||
| 63 | 657 | PetscInt i; | |
| 64 | |||
| 65 |
1/2✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
|
657 | PetscFunctionBegin; |
| 66 |
8/14✓ Branch 0 taken 8 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 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.
|
657 | if (!(*contour)) PetscFunctionReturn(PETSC_SUCCESS); |
| 67 |
1/2✓ Branch 0 taken 8 times.
✗ Branch 1 not taken.
|
447 | if ((*contour)->ksp) { |
| 68 |
7/8✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 8 times.
✓ Branch 3 taken 6 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
✓ Branch 6 taken 2 times.
✓ Branch 7 taken 2 times.
|
10799 | for (i=0;i<(*contour)->npoints;i++) PetscCall(KSPDestroy(&(*contour)->ksp[i])); |
| 69 |
5/8✓ Branch 0 taken 8 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 6 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
|
447 | PetscCall(PetscFree((*contour)->ksp)); |
| 70 | } | ||
| 71 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
447 | PetscCall(PetscSubcommDestroy(&(*contour)->subcomm)); |
| 72 |
5/8✓ Branch 0 taken 8 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 6 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
|
447 | PetscCall(PetscFree((*contour))); |
| 73 | 447 | *contour = NULL; | |
| 74 |
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.
|
447 | PetscFunctionReturn(PETSC_SUCCESS); |
| 75 | } | ||
| 76 | |||
| 77 | /* | ||
| 78 | SlepcContourRedundantMat - Creates redundant copies of the passed matrices in the subcomm. | ||
| 79 | |||
| 80 | Input Parameters: | ||
| 81 | nmat - the number of matrices | ||
| 82 | A - array of matrices | ||
| 83 | P - array of matrices (preconditioner) | ||
| 84 | */ | ||
| 85 | 431 | PetscErrorCode SlepcContourRedundantMat(SlepcContourData contour,PetscInt nmat,Mat *A,Mat *P) | |
| 86 | { | ||
| 87 | 431 | PetscInt i; | |
| 88 | 431 | MPI_Comm child; | |
| 89 | |||
| 90 |
1/2✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
|
431 | PetscFunctionBegin; |
| 91 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 8 times.
|
431 | if (contour->pA) { |
| 92 | ✗ | PetscCall(MatDestroyMatrices(contour->nmat,&contour->pA)); | |
| 93 | ✗ | PetscCall(MatDestroyMatrices(contour->nmat,&contour->pP)); | |
| 94 | ✗ | contour->nmat = 0; | |
| 95 | } | ||
| 96 |
3/4✓ Branch 0 taken 8 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 8 times.
✓ Branch 3 taken 8 times.
|
431 | if (contour->subcomm && contour->subcomm->n != 1) { |
| 97 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
140 | PetscCall(PetscSubcommGetChild(contour->subcomm,&child)); |
| 98 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
140 | PetscCall(PetscCalloc1(nmat,&contour->pA)); |
| 99 |
7/8✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 8 times.
✓ Branch 3 taken 6 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
✓ Branch 6 taken 2 times.
✓ Branch 7 taken 2 times.
|
440 | for (i=0;i<nmat;i++) PetscCall(MatCreateRedundantMatrix(A[i],contour->subcomm->n,child,MAT_INITIAL_MATRIX,&contour->pA[i])); |
| 100 |
2/2✓ Branch 0 taken 8 times.
✓ Branch 1 taken 8 times.
|
140 | if (P) { |
| 101 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
32 | PetscCall(PetscCalloc1(nmat,&contour->pP)); |
| 102 |
7/8✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 8 times.
✓ Branch 3 taken 6 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
✓ Branch 6 taken 2 times.
✓ Branch 7 taken 2 times.
|
112 | for (i=0;i<nmat;i++) PetscCall(MatCreateRedundantMatrix(P[i],contour->subcomm->n,child,MAT_INITIAL_MATRIX,&contour->pP[i])); |
| 103 | } | ||
| 104 | 140 | contour->nmat = nmat; | |
| 105 | } | ||
| 106 |
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.
|
112 | PetscFunctionReturn(PETSC_SUCCESS); |
| 107 | } | ||
| 108 | |||
| 109 | /* | ||
| 110 | SlepcContourScatterCreate - Creates a scatter context to communicate between a | ||
| 111 | regular vector and a vector xdup that can hold one duplicate per each subcommunicator | ||
| 112 | on the contiguous parent communicator. Also creates auxiliary vectors xdup and xsub | ||
| 113 | (the latter with the same layout as the redundant matrices in the subcommunicator). | ||
| 114 | |||
| 115 | Input Parameters: | ||
| 116 | v - the regular vector from which dimensions are taken | ||
| 117 | */ | ||
| 118 | 140 | PetscErrorCode SlepcContourScatterCreate(SlepcContourData contour,Vec v) | |
| 119 | { | ||
| 120 | 140 | IS is1,is2; | |
| 121 | 140 | PetscInt i,j,k,m,mstart,mend,mlocal; | |
| 122 | 140 | PetscInt *idx1,*idx2,mloc_sub; | |
| 123 | 140 | MPI_Comm contpar,parent; | |
| 124 | |||
| 125 |
1/2✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
|
140 | PetscFunctionBegin; |
| 126 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
140 | PetscCall(VecDestroy(&contour->xsub)); |
| 127 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
140 | PetscCall(MatCreateVecsEmpty(contour->pA[0],&contour->xsub,NULL)); |
| 128 | |||
| 129 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
140 | PetscCall(VecDestroy(&contour->xdup)); |
| 130 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
140 | PetscCall(MatGetLocalSize(contour->pA[0],&mloc_sub,NULL)); |
| 131 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
140 | PetscCall(PetscSubcommGetContiguousParent(contour->subcomm,&contpar)); |
| 132 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
140 | PetscCall(VecCreate(contpar,&contour->xdup)); |
| 133 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
140 | PetscCall(VecSetSizes(contour->xdup,mloc_sub,PETSC_DECIDE)); |
| 134 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
140 | PetscCall(VecSetType(contour->xdup,((PetscObject)v)->type_name)); |
| 135 | |||
| 136 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
140 | PetscCall(VecScatterDestroy(&contour->scatterin)); |
| 137 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
140 | PetscCall(VecGetSize(v,&m)); |
| 138 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
140 | PetscCall(VecGetOwnershipRange(v,&mstart,&mend)); |
| 139 | 140 | mlocal = mend - mstart; | |
| 140 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
140 | PetscCall(PetscMalloc2(contour->subcomm->n*mlocal,&idx1,contour->subcomm->n*mlocal,&idx2)); |
| 141 | j = 0; | ||
| 142 |
2/2✓ Branch 0 taken 8 times.
✓ Branch 1 taken 8 times.
|
452 | for (k=0;k<contour->subcomm->n;k++) { |
| 143 |
2/2✓ Branch 0 taken 8 times.
✓ Branch 1 taken 8 times.
|
38888 | for (i=mstart;i<mend;i++) { |
| 144 | 38576 | idx1[j] = i; | |
| 145 | 38576 | idx2[j++] = i + m*k; | |
| 146 | } | ||
| 147 | } | ||
| 148 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
140 | PetscCall(PetscSubcommGetParent(contour->subcomm,&parent)); |
| 149 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
140 | PetscCall(ISCreateGeneral(parent,contour->subcomm->n*mlocal,idx1,PETSC_COPY_VALUES,&is1)); |
| 150 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
140 | PetscCall(ISCreateGeneral(parent,contour->subcomm->n*mlocal,idx2,PETSC_COPY_VALUES,&is2)); |
| 151 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
140 | PetscCall(VecScatterCreate(v,is1,contour->xdup,is2,&contour->scatterin)); |
| 152 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
140 | PetscCall(ISDestroy(&is1)); |
| 153 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
140 | PetscCall(ISDestroy(&is2)); |
| 154 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
140 | PetscCall(PetscFree2(idx1,idx2)); |
| 155 |
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.
|
38 | PetscFunctionReturn(PETSC_SUCCESS); |
| 156 | } | ||
| 157 | |||
| 158 | /* | ||
| 159 | SlepcCISS_isGhost - Determine if any of the computed eigenpairs are spurious. | ||
| 160 | |||
| 161 | Input Parameters: | ||
| 162 | X - the matrix of eigenvectors (MATSEQDENSE) | ||
| 163 | n - the number of columns to consider | ||
| 164 | sigma - the singular values | ||
| 165 | thresh - threshold to decide whether a value is spurious | ||
| 166 | |||
| 167 | Output Parameter: | ||
| 168 | fl - array of n booleans | ||
| 169 | */ | ||
| 170 | 440 | PetscErrorCode SlepcCISS_isGhost(Mat X,PetscInt n,PetscReal *sigma,PetscReal thresh,PetscBool *fl) | |
| 171 | { | ||
| 172 | 440 | const PetscScalar *pX; | |
| 173 | 440 | PetscInt i,j,m,ld; | |
| 174 | 440 | PetscReal *tau,s1,s2,tau_max=0.0; | |
| 175 | |||
| 176 |
1/2✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
|
440 | PetscFunctionBegin; |
| 177 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
440 | PetscCall(MatGetSize(X,&m,NULL)); |
| 178 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
440 | PetscCall(MatDenseGetLDA(X,&ld)); |
| 179 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
440 | PetscCall(PetscMalloc1(n,&tau)); |
| 180 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
440 | PetscCall(MatDenseGetArrayRead(X,&pX)); |
| 181 |
2/2✓ Branch 0 taken 8 times.
✓ Branch 1 taken 8 times.
|
12891 | for (j=0;j<n;j++) { |
| 182 | s1 = 0.0; | ||
| 183 | s2 = 0.0; | ||
| 184 |
2/2✓ Branch 0 taken 8 times.
✓ Branch 1 taken 8 times.
|
702844 | for (i=0;i<m;i++) { |
| 185 | 690393 | s1 += PetscAbsScalar(PetscPowScalarInt(pX[i+j*ld],2)); | |
| 186 | 1380786 | s2 += PetscPowRealInt(PetscAbsScalar(pX[i+j*ld]),2)/sigma[i]; | |
| 187 | } | ||
| 188 | 12451 | tau[j] = s1/s2; | |
| 189 |
2/2✓ Branch 0 taken 8 times.
✓ Branch 1 taken 8 times.
|
14209 | tau_max = PetscMax(tau_max,tau[j]); |
| 190 | } | ||
| 191 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
440 | PetscCall(MatDenseRestoreArrayRead(X,&pX)); |
| 192 |
2/2✓ Branch 0 taken 8 times.
✓ Branch 1 taken 8 times.
|
12891 | for (j=0;j<n;j++) fl[j] = (tau[j]>=thresh*tau_max)? PETSC_TRUE: PETSC_FALSE; |
| 193 |
5/8✓ Branch 0 taken 8 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 6 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
|
440 | PetscCall(PetscFree(tau)); |
| 194 |
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.
|
114 | PetscFunctionReturn(PETSC_SUCCESS); |
| 195 | } | ||
| 196 | |||
| 197 | /* | ||
| 198 | SlepcCISS_BH_SVD - Compute SVD of block Hankel matrix and its rank. | ||
| 199 | |||
| 200 | Input Parameters: | ||
| 201 | H - block Hankel matrix obtained via CISS_BlockHankel() | ||
| 202 | ml - dimension of rows and columns, equal to M*L | ||
| 203 | delta - the tolerance used to determine the rank | ||
| 204 | |||
| 205 | Output Parameters: | ||
| 206 | sigma - computed singular values | ||
| 207 | rank - the rank of H | ||
| 208 | */ | ||
| 209 | 128 | PetscErrorCode SlepcCISS_BH_SVD(PetscScalar *H,PetscInt ml,PetscReal delta,PetscReal *sigma,PetscInt *rank) | |
| 210 | { | ||
| 211 | 128 | PetscInt i; | |
| 212 | 128 | PetscBLASInt m,n,lda,ldu,ldvt,lwork,info; | |
| 213 | 128 | PetscScalar *work; | |
| 214 | #if defined(PETSC_USE_COMPLEX) | ||
| 215 | 103 | PetscReal *rwork; | |
| 216 | #endif | ||
| 217 | |||
| 218 |
1/2✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
|
128 | PetscFunctionBegin; |
| 219 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
128 | PetscCall(PetscMalloc1(5*ml,&work)); |
| 220 | #if defined(PETSC_USE_COMPLEX) | ||
| 221 |
4/6✓ Branch 0 taken 1 times.
✓ Branch 1 taken 3 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
|
103 | PetscCall(PetscMalloc1(5*ml,&rwork)); |
| 222 | #endif | ||
| 223 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
128 | PetscCall(PetscBLASIntCast(ml,&m)); |
| 224 | 128 | n = m; lda = m; ldu = m; ldvt = m; lwork = 5*m; | |
| 225 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
128 | PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF)); |
| 226 | #if defined(PETSC_USE_COMPLEX) | ||
| 227 |
10/20✓ Branch 0 taken 1 times.
✓ Branch 1 taken 3 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.
|
103 | PetscCallBLAS("LAPACKgesvd",LAPACKgesvd_("N","N",&m,&n,H,&lda,sigma,NULL,&ldu,NULL,&ldvt,work,&lwork,rwork,&info)); |
| 228 | #else | ||
| 229 |
10/20✓ Branch 0 taken 1 times.
✓ Branch 1 taken 3 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.
|
25 | PetscCallBLAS("LAPACKgesvd",LAPACKgesvd_("N","N",&m,&n,H,&lda,sigma,NULL,&ldu,NULL,&ldvt,work,&lwork,&info)); |
| 230 | #endif | ||
| 231 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
128 | SlepcCheckLapackInfo("gesvd",info); |
| 232 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
128 | PetscCall(PetscFPTrapPop()); |
| 233 | 128 | (*rank) = 0; | |
| 234 |
2/2✓ Branch 0 taken 8 times.
✓ Branch 1 taken 8 times.
|
7928 | for (i=0;i<ml;i++) { |
| 235 |
4/4✓ Branch 0 taken 8 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 8 times.
✓ Branch 3 taken 8 times.
|
15320 | if (sigma[i]/PetscMax(sigma[0],1.0)>delta) (*rank)++; |
| 236 | } | ||
| 237 |
5/8✓ Branch 0 taken 8 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 6 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
|
128 | PetscCall(PetscFree(work)); |
| 238 | #if defined(PETSC_USE_COMPLEX) | ||
| 239 |
5/8✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1 times.
✓ Branch 3 taken 3 times.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 1 times.
|
103 | PetscCall(PetscFree(rwork)); |
| 240 | #endif | ||
| 241 |
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.
|
32 | PetscFunctionReturn(PETSC_SUCCESS); |
| 242 | } | ||
| 243 |