GCC Code Coverage Report

Directory: ./
File: src/sys/slepccontour.c
Date: 2025-10-04 04:19:13
Exec Total Coverage
Lines: 119 122 97.5%
Functions: 7 7 100.0%
Branches: 328 526 62.4%
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 557 PetscErrorCode SlepcContourDataCreate(PetscInt n,PetscInt npart,PetscObject parent,SlepcContourData *contour)
23 {
24
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 557 PetscFunctionBegin;
25
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 557 PetscCall(PetscNew(contour));
26 557 (*contour)->parent = parent;
27
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 557 PetscCall(PetscSubcommCreate(PetscObjectComm(parent),&(*contour)->subcomm));
28
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 557 PetscCall(PetscSubcommSetNumber((*contour)->subcomm,npart));
29
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 557 PetscCall(PetscSubcommSetType((*contour)->subcomm,PETSC_SUBCOMM_INTERLACED));
30 557 (*contour)->npoints = n / npart;
31
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
 557 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 419 PetscErrorCode SlepcContourDataReset(SlepcContourData contour)
40 {
41 419 PetscInt i;
42
43
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 419 PetscFunctionBegin;
44
1/2
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
 419 if (contour->ksp) {
45
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.
 9627 for (i=0;i<contour->npoints;i++) PetscCall(KSPReset(contour->ksp[i]));
46 }
47
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 419 if (contour->pA) {
48
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 174 PetscCall(MatDestroyMatrices(contour->nmat,&contour->pA));
49
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 174 PetscCall(MatDestroyMatrices(contour->nmat,&contour->pP));
50 174 contour->nmat = 0;
51 }
52
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 419 PetscCall(VecScatterDestroy(&contour->scatterin));
53
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 419 PetscCall(VecDestroy(&contour->xsub));
54
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 419 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 818 PetscErrorCode SlepcContourDataDestroy(SlepcContourData *contour)
62 {
63 818 PetscInt i;
64
65
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 818 PetscFunctionBegin;
66
8/14
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 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.
 818 if (!(*contour)) PetscFunctionReturn(PETSC_SUCCESS);
67
1/2
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
 557 if ((*contour)->ksp) {
68
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.
 13461 for (i=0;i<(*contour)->npoints;i++) PetscCall(KSPDestroy(&(*contour)->ksp[i]));
69
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.
 557 PetscCall(PetscFree((*contour)->ksp));
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.
 557 PetscCall(PetscSubcommDestroy(&(*contour)->subcomm));
72
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.
 557 PetscCall(PetscFree((*contour)));
73 557 *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.
 557 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 537 PetscErrorCode SlepcContourRedundantMat(SlepcContourData contour,PetscInt nmat,Mat *A,Mat *P)
86 {
87 537 PetscInt i;
88 537 MPI_Comm child;
89
90
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 537 PetscFunctionBegin;
91
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
 537 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 10 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 10 times.
✓ Branch 3 taken 10 times.
 537 if (contour->subcomm && contour->subcomm->n != 1) {
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.
 174 PetscCall(PetscSubcommGetChild(contour->subcomm,&child));
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.
 174 PetscCall(PetscCalloc1(nmat,&contour->pA));
99
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.
 548 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 10 times.
✓ Branch 1 taken 10 times.
 174 if (P) {
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.
 40 PetscCall(PetscCalloc1(nmat,&contour->pP));
102
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.
 140 for (i=0;i<nmat;i++) PetscCall(MatCreateRedundantMatrix(P[i],contour->subcomm->n,child,MAT_INITIAL_MATRIX,&contour->pP[i]));
103 }
104 174 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 174 PetscErrorCode SlepcContourScatterCreate(SlepcContourData contour,Vec v)
119 {
120 174 IS is1,is2;
121 174 PetscInt i,j,k,m,mstart,mend,mlocal;
122 174 PetscInt *idx1,*idx2,mloc_sub;
123 174 MPI_Comm contpar,parent;
124
125
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 174 PetscFunctionBegin;
126
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 174 PetscCall(VecDestroy(&contour->xsub));
127
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 174 PetscCall(MatCreateVecsEmpty(contour->pA[0],&contour->xsub,NULL));
128
129
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 174 PetscCall(VecDestroy(&contour->xdup));
130
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 174 PetscCall(MatGetLocalSize(contour->pA[0],&mloc_sub,NULL));
131
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 174 PetscCall(PetscSubcommGetContiguousParent(contour->subcomm,&contpar));
132
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 174 PetscCall(VecCreate(contpar,&contour->xdup));
133
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 174 PetscCall(VecSetSizes(contour->xdup,mloc_sub,PETSC_DECIDE));
134
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 174 PetscCall(VecSetType(contour->xdup,((PetscObject)v)->type_name));
135
136
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 174 PetscCall(VecScatterDestroy(&contour->scatterin));
137
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 174 PetscCall(VecGetSize(v,&m));
138
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 174 PetscCall(VecGetOwnershipRange(v,&mstart,&mend));
139 174 mlocal = mend - mstart;
140
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 174 PetscCall(PetscMalloc2(contour->subcomm->n*mlocal,&idx1,contour->subcomm->n*mlocal,&idx2));
141 j = 0;
142
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 562 for (k=0;k<contour->subcomm->n;k++) {
143
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 47708 for (i=mstart;i<mend;i++) {
144 47320 idx1[j] = i;
145 47320 idx2[j++] = i + m*k;
146 }
147 }
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.
 174 PetscCall(PetscSubcommGetParent(contour->subcomm,&parent));
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.
 174 PetscCall(ISCreateGeneral(parent,contour->subcomm->n*mlocal,idx1,PETSC_COPY_VALUES,&is1));
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.
 174 PetscCall(ISCreateGeneral(parent,contour->subcomm->n*mlocal,idx2,PETSC_COPY_VALUES,&is2));
151
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 174 PetscCall(VecScatterCreate(v,is1,contour->xdup,is2,&contour->scatterin));
152
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 174 PetscCall(ISDestroy(&is1));
153
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 174 PetscCall(ISDestroy(&is2));
154
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 174 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 548 PetscErrorCode SlepcCISS_isGhost(Mat X,PetscInt n,PetscReal *sigma,PetscReal thresh,PetscBool *fl)
171 {
172 548 const PetscScalar *pX;
173 548 PetscInt i,j,m,ld;
174 548 PetscReal *tau,s1,s2,tau_max=0.0;
175
176
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 548 PetscFunctionBegin;
177
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 548 PetscCall(MatGetSize(X,&m,NULL));
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.
 548 PetscCall(MatDenseGetLDA(X,&ld));
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.
 548 PetscCall(PetscMalloc1(n,&tau));
180
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 548 PetscCall(MatDenseGetArrayRead(X,&pX));
181
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 15893 for (j=0;j<n;j++) {
182 s1 = 0.0;
183 s2 = 0.0;
184
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 859546 for (i=0;i<m;i++) {
185 844201 s1 += PetscAbsScalar(PetscPowScalarInt(pX[i+j*ld],2));
186 1688402 s2 += PetscPowRealInt(PetscAbsScalar(pX[i+j*ld]),2)/sigma[i];
187 }
188 15345 tau[j] = s1/s2;
189
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 17536 tau_max = PetscMax(tau_max,tau[j]);
190 }
191
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 548 PetscCall(MatDenseRestoreArrayRead(X,&pX));
192
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 15893 for (j=0;j<n;j++) fl[j] = (tau[j]>=thresh*tau_max)? PETSC_TRUE: PETSC_FALSE;
193
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.
 548 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 160 PetscErrorCode SlepcCISS_BH_SVD(PetscScalar *H,PetscInt ml,PetscReal delta,PetscReal *sigma,PetscInt *rank)
210 {
211 160 PetscInt i;
212 160 PetscBLASInt m,n,lda,ldu,ldvt,lwork,info;
213 160 PetscScalar *work;
214 #if defined(PETSC_USE_COMPLEX)
215 128 PetscReal *rwork;
216 #endif
217
218
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 160 PetscFunctionBegin;
219
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 160 PetscCall(PetscMalloc1(5*ml,&work));
220 #if defined(PETSC_USE_COMPLEX)
221
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.
 128 PetscCall(PetscMalloc1(5*ml,&rwork));
222 #endif
223
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 160 PetscCall(PetscBLASIntCast(ml,&m));
224 160 n = m; lda = m; ldu = m; ldvt = m; lwork = 5*m;
225
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 160 PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
226 #if defined(PETSC_USE_COMPLEX)
227
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.
 128 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 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.
 32 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 10 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 160 SlepcCheckLapackInfo("gesvd",info);
232
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 160 PetscCall(PetscFPTrapPop());
233 160 (*rank) = 0;
234
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 9856 for (i=0;i<ml;i++) {
235
4/4
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 5 times.
✓ Branch 2 taken 10 times.
✓ Branch 3 taken 10 times.
 19000 if (sigma[i]/PetscMax(sigma[0],1.0)>delta) (*rank)++;
236 }
237
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.
 160 PetscCall(PetscFree(work));
238 #if defined(PETSC_USE_COMPLEX)
239
5/8
✓ 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 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 1 times.
 128 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