| 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 | A polynomial filter based on a truncated Chebyshev series. | ||
| 12 | */ | ||
| 13 | |||
| 14 | #include <slepc/private/stimpl.h> | ||
| 15 | #include "filter.h" | ||
| 16 | |||
| 17 | #define DAMPING_LANCZOS_POW 2 | ||
| 18 | |||
| 19 | /* | ||
| 20 | Gateway to Chebyshev for evaluating y=p(A)*x | ||
| 21 | (Clenshaw with a vector) | ||
| 22 | */ | ||
| 23 | 6681 | static PetscErrorCode MatMult_Chebyshev(Mat A,Vec x,Vec y) | |
| 24 | { | ||
| 25 | 6681 | ST st; | |
| 26 | 6681 | ST_FILTER *ctx; | |
| 27 | 6681 | CHEBYSHEV_CTX cheby; | |
| 28 | 6681 | PetscInt degree,i,w0=0,w1=1; | |
| 29 | 6681 | PetscReal s1,s2; | |
| 30 | |||
| 31 |
1/2✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
|
6681 | PetscFunctionBegin; |
| 32 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
6681 | PetscCall(MatShellGetContext(A,&st)); |
| 33 | 6681 | ctx = (ST_FILTER*)st->data; | |
| 34 | 6681 | cheby = (CHEBYSHEV_CTX)ctx->data; | |
| 35 | 6681 | degree = ctx->polyDegree; | |
| 36 | 6681 | s1 = 4/(ctx->right - ctx->left); | |
| 37 | 6681 | s2 = -2*(ctx->right+ctx->left)/(ctx->right-ctx->left); | |
| 38 | |||
| 39 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
6681 | PetscCall(VecSet(st->work[w0],0.0)); |
| 40 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
6681 | PetscCall(VecCopy(x,y)); |
| 41 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
6681 | PetscCall(VecScale(y,cheby->damping_coeffs[degree]*cheby->coeffs[degree])); |
| 42 | |||
| 43 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
832881 | for (i=0;i<degree;i++) { |
| 44 | /* swap work vectors, save a copy of current y */ | ||
| 45 | 826200 | w0 = 1-w0; | |
| 46 | 826200 | w1 = 1-w1; | |
| 47 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
826200 | PetscCall(VecCopy(y,st->work[w0])); |
| 48 | /* compute next y */ | ||
| 49 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
826200 | if (i == degree-1) { |
| 50 | 6681 | s1 = s1/2; | |
| 51 | 6681 | s2 = s2/2; | |
| 52 | } | ||
| 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.
|
826200 | PetscCall(VecAXPBYPCZ(y,cheby->damping_coeffs[degree-i-1]*cheby->coeffs[degree-i-1],-1,s2,x,st->work[w1])); |
| 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.
|
826200 | PetscCall(MatMult(ctx->T,st->work[w0],st->work[w1])); |
| 55 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
826200 | PetscCall(VecAXPY(y,s1,st->work[w1])); |
| 56 | } | ||
| 57 |
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.
|
1263 | PetscFunctionReturn(PETSC_SUCCESS); |
| 58 | } | ||
| 59 | |||
| 60 | /* | ||
| 61 | Gateway to Chebyshev for evaluating C=p(A)*B | ||
| 62 | (Clenshaw with a matrix) | ||
| 63 | */ | ||
| 64 | 512 | static PetscErrorCode MatMatMult_Chebyshev(Mat A,Mat B,Mat C,void *pctx) | |
| 65 | { | ||
| 66 | 512 | ST st; | |
| 67 | 512 | ST_FILTER *ctx; | |
| 68 | 512 | CHEBYSHEV_CTX cheby; | |
| 69 | 512 | PetscInt degree,i,m1,m2,w0=0,w1=1; | |
| 70 | 512 | PetscReal s1,s2; | |
| 71 | |||
| 72 |
1/2✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
|
512 | PetscFunctionBegin; |
| 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.
|
512 | PetscCall(MatShellGetContext(A,&st)); |
| 74 | 512 | ctx = (ST_FILTER*)st->data; | |
| 75 | 512 | cheby = (CHEBYSHEV_CTX)ctx->data; | |
| 76 | 512 | degree = ctx->polyDegree; | |
| 77 | 512 | s1 = 4/(ctx->right - ctx->left); | |
| 78 | 512 | s2 = -2*(ctx->right+ctx->left)/(ctx->right-ctx->left); | |
| 79 | |||
| 80 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
512 | if (ctx->nW) { /* check if work matrices must be resized */ |
| 81 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
498 | PetscCall(MatGetSize(B,NULL,&m1)); |
| 82 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
498 | PetscCall(MatGetSize(ctx->W[0],NULL,&m2)); |
| 83 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
498 | if (m1!=m2) { |
| 84 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
57 | PetscCall(MatDestroyMatrices(ctx->nW,&ctx->W)); |
| 85 | 57 | ctx->nW = 0; | |
| 86 | } | ||
| 87 | } | ||
| 88 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
512 | if (!ctx->nW) { /* allocate work matrices */ |
| 89 | 71 | ctx->nW = 2; | |
| 90 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
71 | PetscCall(PetscMalloc1(ctx->nW,&ctx->W)); |
| 91 |
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.
|
213 | for (i=0;i<ctx->nW;i++) PetscCall(MatDuplicate(B,MAT_DO_NOT_COPY_VALUES,&ctx->W[i])); |
| 92 | } | ||
| 93 | |||
| 94 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
512 | PetscCall(MatScale(ctx->W[w0],0)); |
| 95 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
512 | PetscCall(MatCopy(B,C,SAME_NONZERO_PATTERN)); |
| 96 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
512 | PetscCall(MatScale(C,cheby->damping_coeffs[degree]*cheby->coeffs[degree])); |
| 97 | |||
| 98 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
102912 | for (i=0;i<degree;i++) { |
| 99 | /* swap work matrices, save a copy of current C */ | ||
| 100 | 102400 | w0 = 1-w0; | |
| 101 | 102400 | w1 = 1-w1; | |
| 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.
|
102400 | PetscCall(MatCopy(C,ctx->W[w0],SAME_NONZERO_PATTERN)); |
| 103 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
|
102400 | if (i == degree-1) { |
| 104 | 512 | s1 = s1/2; | |
| 105 | 512 | s2 = s2/2; | |
| 106 | } | ||
| 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.
|
102400 | PetscCall(MatScale(C,s2)); |
| 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.
|
102400 | PetscCall(MatAXPY(C,-1,ctx->W[w1],SAME_NONZERO_PATTERN)); |
| 109 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
102400 | PetscCall(MatAXPY(C,cheby->damping_coeffs[degree-i-1]*cheby->coeffs[degree-i-1],B,SAME_NONZERO_PATTERN)); |
| 110 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
102400 | PetscCall(MatMatMult(ctx->T,ctx->W[w0],MAT_REUSE_MATRIX,PETSC_CURRENT,&ctx->W[w1])); |
| 111 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
102400 | PetscCall(MatAXPY(C,s1,ctx->W[w1],SAME_NONZERO_PATTERN)); |
| 112 | } | ||
| 113 |
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.
|
73 | PetscFunctionReturn(PETSC_SUCCESS); |
| 114 | } | ||
| 115 | |||
| 116 | 98 | static PetscErrorCode Chebyshev_Compute_Coefficients(ST st) | |
| 117 | { | ||
| 118 | 98 | ST_FILTER *ctx = (ST_FILTER*)st->data; | |
| 119 | 98 | CHEBYSHEV_CTX cheby = (CHEBYSHEV_CTX)ctx->data; | |
| 120 | 98 | PetscInt i,degree; | |
| 121 | 98 | PetscReal s_inta,s_intb,acosa,acosb,t1,t2,s; | |
| 122 | |||
| 123 |
1/2✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
|
98 | PetscFunctionBegin; |
| 124 | 98 | degree = ctx->polyDegree; | |
| 125 |
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.
|
98 | PetscCall(PetscFree(cheby->coeffs)); |
| 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.
|
98 | PetscCall(PetscMalloc1(degree+1,&cheby->coeffs)); |
| 127 |
1/2✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
|
98 | s_inta = PetscMax((2.0*ctx->inta - (ctx->right + ctx->left))/(ctx->right - ctx->left),-1.0); |
| 128 |
1/2✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
|
98 | s_intb = PetscMin((2.0*ctx->intb - (ctx->right + ctx->left))/(ctx->right - ctx->left),1.0); |
| 129 | 98 | acosa = PetscAcosReal(s_inta); | |
| 130 | 98 | acosb = PetscAcosReal(s_intb); | |
| 131 | 98 | t1 = (acosa - acosb)/PETSC_PI; | |
| 132 | 98 | t2 = t1; | |
| 133 | 98 | cheby->coeffs[0] = t1; | |
| 134 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
13198 | for (i=0;i<degree;i++) { |
| 135 | 13100 | cheby->coeffs[i+1] = 2.0*(PetscSinReal((i+1)*acosa) - PetscSinReal((i+1)*acosb)) / ((i+1)*PETSC_PI); | |
| 136 | 13100 | t1 = t1 + cheby->coeffs[i+1]*PetscCosReal((i+1)*acosa); | |
| 137 | 13100 | t2 = t2 + cheby->coeffs[i+1]*PetscCosReal((i+1)*acosb); | |
| 138 | } | ||
| 139 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
98 | if (s_inta == -1.0) s = 0.5/PetscAbs(t1); |
| 140 |
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.
|
98 | else if (s_intb == 1.0) s = 0.5/PetscMin(PetscAbs(t1),PetscAbs(t2)); |
| 141 |
4/6✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 10 times.
✓ Branch 4 taken 10 times.
✓ Branch 5 taken 10 times.
|
168 | else s = 0.5/PetscMin(PetscAbs(t1),PetscAbs(t2)); |
| 142 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
13296 | for (i=0;i<degree+1;i++) cheby->coeffs[i] = s*cheby->coeffs[i]; |
| 143 |
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.
|
18 | PetscFunctionReturn(PETSC_SUCCESS); |
| 144 | } | ||
| 145 | |||
| 146 | 98 | static PetscErrorCode Chebyshev_Compute_Damping_Coefficients(ST st) | |
| 147 | { | ||
| 148 | 98 | ST_FILTER *ctx = (ST_FILTER*)st->data; | |
| 149 | 98 | CHEBYSHEV_CTX cheby = (CHEBYSHEV_CTX)ctx->data; | |
| 150 | 98 | PetscInt i,degree; | |
| 151 | 98 | PetscReal invdeg; | |
| 152 | |||
| 153 |
1/2✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
|
98 | PetscFunctionBegin; |
| 154 | 98 | degree = ctx->polyDegree; | |
| 155 |
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.
|
98 | PetscCall(PetscFree(cheby->damping_coeffs)); |
| 156 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
98 | PetscCall(PetscMalloc1(degree+1,&cheby->damping_coeffs)); |
| 157 | 98 | cheby->damping_coeffs[0] = 1; | |
| 158 |
4/5✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
✓ Branch 2 taken 10 times.
✓ Branch 3 taken 10 times.
✗ Branch 4 not taken.
|
98 | switch (ctx->damping) { |
| 159 | case ST_FILTER_DAMPING_NONE: | ||
| 160 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
8148 | for (i=1;i<=degree;i++) cheby->damping_coeffs[i] = 1; |
| 161 | break; | ||
| 162 | case ST_FILTER_DAMPING_LANCZOS: | ||
| 163 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
2020 | for (i=1;i<=degree;i++) cheby->damping_coeffs[i] = PetscPowRealInt(PetscSinReal(i*PETSC_PI/(degree+1))/(i*PETSC_PI/(degree+1)),DAMPING_LANCZOS_POW); |
| 164 | break; | ||
| 165 | 10 | case ST_FILTER_DAMPING_JACKSON: | |
| 166 | 10 | invdeg = 1/((PetscReal)degree+2); | |
| 167 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
1010 | for (i=1;i<=degree;i++) cheby->damping_coeffs[i] = invdeg*((degree+2-i)*PetscCosReal(PETSC_PI*i*invdeg)+PetscSinReal(PETSC_PI*i*invdeg)/PetscTanReal(PETSC_PI*invdeg)); |
| 168 | break; | ||
| 169 | case ST_FILTER_DAMPING_FEJER: | ||
| 170 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
2020 | for (i=1;i<=degree;i++) cheby->damping_coeffs[i] = 1-i/(PetscReal)(degree+1); |
| 171 | break; | ||
| 172 | ✗ | default: | |
| 173 | ✗ | SETERRQ(PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_WRONG,"Invalid Chebyshev filter damping type"); | |
| 174 | } | ||
| 175 |
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.
|
18 | PetscFunctionReturn(PETSC_SUCCESS); |
| 176 | } | ||
| 177 | |||
| 178 | /* | ||
| 179 | Computes the coefficients and damping coefficients for the Chebyshev filter | ||
| 180 | |||
| 181 | Creates the shifted (and scaled) matrix and the filter P(z). | ||
| 182 | G is a shell matrix whose MatMult() and MatMatMult() applies the filter. | ||
| 183 | */ | ||
| 184 | 98 | static PetscErrorCode STComputeOperator_Filter_Chebyshev(ST st,Mat *G) | |
| 185 | { | ||
| 186 | 98 | ST_FILTER *ctx = (ST_FILTER*)st->data; | |
| 187 | 98 | PetscInt n,m,N,M; | |
| 188 | |||
| 189 |
1/2✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
|
98 | PetscFunctionBegin; |
| 190 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
98 | PetscCall(STSetWorkVecs(st,2)); |
| 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.
|
98 | PetscCall(Chebyshev_Compute_Coefficients(st)); |
| 192 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
98 | PetscCall(Chebyshev_Compute_Damping_Coefficients(st)); |
| 193 | /* create shell matrix*/ | ||
| 194 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
98 | PetscCall(MatDestroy(&ctx->T)); |
| 195 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
98 | PetscCall(MatDuplicate(st->A[0],MAT_COPY_VALUES,&ctx->T)); |
| 196 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
98 | if (!*G) { |
| 197 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
78 | PetscCall(MatGetSize(ctx->T,&N,&M)); |
| 198 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
78 | PetscCall(MatGetLocalSize(ctx->T,&n,&m)); |
| 199 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
78 | PetscCall(MatCreateShell(PetscObjectComm((PetscObject)st),n,m,N,M,st,G)); |
| 200 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
78 | PetscCall(MatShellSetOperation(*G,MATOP_MULT,(PetscErrorCodeFn*)MatMult_Chebyshev)); |
| 201 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
78 | PetscCall(MatShellSetMatProductOperation(*G,MATPRODUCT_AB,NULL,MatMatMult_Chebyshev,NULL,MATDENSE,MATDENSE)); |
| 202 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
78 | PetscCall(MatShellSetMatProductOperation(*G,MATPRODUCT_AB,NULL,MatMatMult_Chebyshev,NULL,MATDENSECUDA,MATDENSECUDA)); |
| 203 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
78 | PetscCall(MatShellSetMatProductOperation(*G,MATPRODUCT_AB,NULL,MatMatMult_Chebyshev,NULL,MATDENSEHIP,MATDENSEHIP)); |
| 204 | } | ||
| 205 |
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.
|
18 | PetscFunctionReturn(PETSC_SUCCESS); |
| 206 | } | ||
| 207 | |||
| 208 | 334 | static PetscErrorCode Chebyshev_Clenshaw_Damping_Real(ST st,PetscReal x,PetscReal *px) | |
| 209 | { | ||
| 210 | 334 | ST_FILTER *ctx = (ST_FILTER*)st->data; | |
| 211 | 334 | CHEBYSHEV_CTX cheby = (CHEBYSHEV_CTX)ctx->data; | |
| 212 | 334 | PetscReal t1=0,t2=0; | |
| 213 | 334 | PetscInt i,degree; | |
| 214 | |||
| 215 |
1/2✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
|
334 | PetscFunctionBegin; |
| 216 | 334 | degree = ctx->polyDegree; | |
| 217 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
334 | PetscCheck(cheby->coeffs,PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_NULL,"Chebyshev coefficients are not computed"); |
| 218 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
334 | PetscCheck(cheby->damping_coeffs,PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_NULL,"Chebyshev damping coefficients are not computed"); |
| 219 | 334 | x = 2*x; | |
| 220 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
30234 | for (i=degree;i>1;i=i-2) { |
| 221 | 29900 | t2 = cheby->damping_coeffs[i]*cheby->coeffs[i] + x*t1 - t2; | |
| 222 | 29900 | t1 = cheby->damping_coeffs[i-1]*cheby->coeffs[i-1] + x*t2 - t1; | |
| 223 | } | ||
| 224 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
|
334 | if (degree % 2) { |
| 225 | ✗ | t2 = cheby->damping_coeffs[1]*cheby->coeffs[1] + x*t1 - t2; | |
| 226 | ✗ | t1 = t2; | |
| 227 | } | ||
| 228 | 334 | *px = cheby->damping_coeffs[0]*cheby->coeffs[0] + .5*x*t1 - t2; | |
| 229 |
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.
|
334 | PetscFunctionReturn(PETSC_SUCCESS); |
| 230 | } | ||
| 231 | |||
| 232 | /* | ||
| 233 | Computes the threshold for the Chebyshev filter | ||
| 234 | |||
| 235 | The threshold is P(scale(ctx->inta)), where ctx->inta is the lower bound of the interval. | ||
| 236 | It should be the same value if the upper bound was used instead. | ||
| 237 | */ | ||
| 238 | 334 | static PetscErrorCode STFilterGetThreshold_Filter_Chebyshev(ST st,PetscReal *gamma) | |
| 239 | { | ||
| 240 | 334 | ST_FILTER *ctx = (ST_FILTER*)st->data; | |
| 241 | 334 | PetscReal c,e; | |
| 242 | |||
| 243 |
1/2✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
|
334 | PetscFunctionBegin; |
| 244 | /* scale ctx->inta */ | ||
| 245 | 334 | c = (ctx->right + ctx->left) / 2.0; | |
| 246 | 334 | e = (ctx->right - ctx->left) / 2.0; | |
| 247 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
334 | PetscCall(Chebyshev_Clenshaw_Damping_Real(st,(ctx->inta-c)/e,gamma)); |
| 248 |
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.
|
52 | PetscFunctionReturn(PETSC_SUCCESS); |
| 249 | } | ||
| 250 | |||
| 251 | 78 | static PetscErrorCode STDestroy_Filter_Chebyshev(ST st) | |
| 252 | { | ||
| 253 | 78 | ST_FILTER *ctx = (ST_FILTER*)st->data; | |
| 254 | 78 | CHEBYSHEV_CTX cheby = (CHEBYSHEV_CTX)ctx->data; | |
| 255 | |||
| 256 |
1/2✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
|
78 | PetscFunctionBegin; |
| 257 |
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.
|
78 | PetscCall(PetscFree(cheby->coeffs)); |
| 258 |
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.
|
78 | PetscCall(PetscFree(cheby->damping_coeffs)); |
| 259 |
6/8✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 8 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
|
78 | PetscCall(PetscFree(cheby)); |
| 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.
|
14 | PetscFunctionReturn(PETSC_SUCCESS); |
| 261 | } | ||
| 262 | |||
| 263 | 78 | PetscErrorCode STCreate_Filter_Chebyshev(ST st) | |
| 264 | { | ||
| 265 | 78 | ST_FILTER *ctx = (ST_FILTER*)st->data; | |
| 266 | 78 | CHEBYSHEV_CTX cheby; | |
| 267 | |||
| 268 |
1/2✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
|
78 | PetscFunctionBegin; |
| 269 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
78 | PetscCall(PetscNew(&cheby)); |
| 270 | 78 | ctx->data = (void*)cheby; | |
| 271 | |||
| 272 | 78 | ctx->computeoperator = STComputeOperator_Filter_Chebyshev; | |
| 273 | 78 | ctx->getthreshold = STFilterGetThreshold_Filter_Chebyshev; | |
| 274 | 78 | ctx->destroy = STDestroy_Filter_Chebyshev; | |
| 275 |
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.
|
78 | PetscFunctionReturn(PETSC_SUCCESS); |
| 276 | } | ||
| 277 |