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 | SLEPc matrix function solver: "krylov" | ||
12 | |||
13 | Method: Arnoldi with Eiermann-Ernst restart | ||
14 | |||
15 | Algorithm: | ||
16 | |||
17 | Build Arnoldi approximations using f(H) for the Hessenberg matrix H, | ||
18 | restart by discarding the Krylov basis but keeping H. | ||
19 | |||
20 | References: | ||
21 | |||
22 | [1] M. Eiermann and O. Ernst, "A restarted Krylov subspace method | ||
23 | for the evaluation of matrix functions", SIAM J. Numer. Anal. | ||
24 | 44(6):2481-2504, 2006. | ||
25 | */ | ||
26 | |||
27 | #include <slepc/private/mfnimpl.h> | ||
28 | #include <slepcblaslapack.h> | ||
29 | |||
30 | 160 | static PetscErrorCode MFNSetUp_Krylov(MFN mfn) | |
31 | { | ||
32 | 160 | PetscInt N; | |
33 | |||
34 |
1/2✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
|
160 | PetscFunctionBegin; |
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.
|
160 | PetscCall(MatGetSize(mfn->A,&N,NULL)); |
36 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
160 | if (mfn->ncv==PETSC_DETERMINE) mfn->ncv = PetscMin(30,N); |
37 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
160 | if (mfn->max_it==PETSC_DETERMINE) mfn->max_it = 100; |
38 |
4/6✓ Branch 0 taken 2 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(MFNAllocateSolution(mfn,1)); |
39 |
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.
|
29 | PetscFunctionReturn(PETSC_SUCCESS); |
40 | } | ||
41 | |||
42 | 964 | static PetscErrorCode MFNSolve_Krylov(MFN mfn,Vec b,Vec x) | |
43 | { | ||
44 | 964 | PetscInt n=0,m,ld,ldh,j; | |
45 | 964 | PetscBLASInt m_,inc=1; | |
46 | 964 | Mat M,G=NULL,H=NULL; | |
47 | 964 | Vec F=NULL; | |
48 | 964 | PetscScalar *marray,*farray,*harray; | |
49 | 964 | const PetscScalar *garray; | |
50 | 964 | PetscReal beta,betaold=0.0,nrm=1.0; | |
51 | 964 | PetscBool breakdown; | |
52 | |||
53 |
1/2✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
|
964 | PetscFunctionBegin; |
54 | 964 | m = mfn->ncv; | |
55 | 964 | ld = m+1; | |
56 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
964 | PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,ld,m,NULL,&M)); |
57 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
964 | PetscCall(MatDenseGetArray(M,&marray)); |
58 | |||
59 | /* set initial vector to b/||b|| */ | ||
60 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
964 | PetscCall(BVInsertVec(mfn->V,0,b)); |
61 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
964 | PetscCall(BVScaleColumn(mfn->V,0,1.0/mfn->bnorm)); |
62 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
964 | PetscCall(VecSet(x,0.0)); |
63 | |||
64 | /* Restart loop */ | ||
65 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
3262 | while (mfn->reason == MFN_CONVERGED_ITERATING) { |
66 | 2298 | mfn->its++; | |
67 | |||
68 | /* compute Arnoldi factorization */ | ||
69 |
7/8✓ Branch 0 taken 6 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 8 times.
✓ Branch 4 taken 2 times.
✓ Branch 5 taken 2 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
|
2298 | PetscCall(BVMatArnoldi(mfn->V,mfn->transpose_solve?mfn->AT:mfn->A,M,0,&m,&beta,&breakdown)); |
70 | |||
71 | /* save previous Hessenberg matrix in G; allocate new storage for H and f(H) */ | ||
72 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
2298 | if (mfn->its>1) { G = H; H = NULL; } |
73 | 2298 | ldh = n+m; | |
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.
|
2298 | PetscCall(MFN_CreateVec(ldh,&F)); |
75 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
2298 | PetscCall(MFN_CreateDenseMat(ldh,&H)); |
76 | |||
77 | /* glue together the previous H and the new H obtained with Arnoldi */ | ||
78 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
2298 | PetscCall(MatDenseGetArray(H,&harray)); |
79 |
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.
|
51270 | for (j=0;j<m;j++) PetscCall(PetscArraycpy(harray+n+(j+n)*ldh,marray+j*ld,m)); |
80 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
2298 | if (mfn->its>1) { |
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.
|
1334 | PetscCall(MatDenseGetArrayRead(G,&garray)); |
82 |
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.
|
44270 | for (j=0;j<n;j++) PetscCall(PetscArraycpy(harray+j*ldh,garray+j*n,n)); |
83 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
1334 | PetscCall(MatDenseRestoreArrayRead(G,&garray)); |
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.
|
1334 | PetscCall(MatDestroy(&G)); |
85 | 1334 | harray[n+(n-1)*ldh] = betaold; | |
86 | } | ||
87 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
2298 | PetscCall(MatDenseRestoreArray(H,&harray)); |
88 | |||
89 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
2298 | if (mfn->its==1) { |
90 | /* set symmetry flag of H from A */ | ||
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.
|
964 | PetscCall(MatPropagateSymmetryOptions(mfn->A,H)); |
92 | } | ||
93 | |||
94 | /* evaluate f(H) */ | ||
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.
|
2298 | PetscCall(FNEvaluateFunctionMatVec(mfn->fn,H,F)); |
96 | |||
97 | /* x += ||b||*V*f(H)*e_1 */ | ||
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.
|
2298 | PetscCall(VecGetArray(F,&farray)); |
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.
|
2298 | PetscCall(PetscBLASIntCast(m,&m_)); |
100 | 2298 | nrm = BLASnrm2_(&m_,farray+n,&inc); /* relative norm of the update ||u||/||b|| */ | |
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.
|
2298 | PetscCall(MFNMonitor(mfn,mfn->its,nrm)); |
102 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
51270 | for (j=0;j<m;j++) farray[j+n] *= mfn->bnorm; |
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.
|
2298 | PetscCall(BVSetActiveColumns(mfn->V,0,m)); |
104 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
2298 | PetscCall(BVMultVec(mfn->V,1.0,1.0,x,farray+n)); |
105 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
2298 | PetscCall(VecRestoreArray(F,&farray)); |
106 | |||
107 | /* check convergence */ | ||
108 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
|
2298 | if (mfn->its >= mfn->max_it) mfn->reason = MFN_DIVERGED_ITS; |
109 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
2298 | if (mfn->its>1) { |
110 |
5/8✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 10 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 10 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 10 times.
✓ Branch 7 taken 10 times.
|
1334 | if (m<mfn->ncv || breakdown || beta==0.0 || nrm<mfn->tol) mfn->reason = MFN_CONVERGED_TOL; |
111 | } | ||
112 | |||
113 | /* restart with vector v_{m+1} */ | ||
114 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
2298 | if (mfn->reason == MFN_CONVERGED_ITERATING) { |
115 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
1334 | PetscCall(BVCopyColumn(mfn->V,m,0)); |
116 | 1334 | n += m; | |
117 | 1334 | betaold = beta; | |
118 | } | ||
119 | } | ||
120 | |||
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.
|
964 | PetscCall(MatDestroy(&H)); |
122 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
964 | PetscCall(MatDestroy(&G)); |
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.
|
964 | PetscCall(VecDestroy(&F)); |
124 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
964 | PetscCall(MatDenseRestoreArray(M,&marray)); |
125 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
964 | PetscCall(MatDestroy(&M)); |
126 |
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.
|
189 | PetscFunctionReturn(PETSC_SUCCESS); |
127 | } | ||
128 | |||
129 | 174 | SLEPC_EXTERN PetscErrorCode MFNCreate_Krylov(MFN mfn) | |
130 | { | ||
131 |
1/2✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
|
174 | PetscFunctionBegin; |
132 | 174 | mfn->ops->solve = MFNSolve_Krylov; | |
133 | 174 | mfn->ops->setup = MFNSetUp_Krylov; | |
134 |
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.
|
174 | PetscFunctionReturn(PETSC_SUCCESS); |
135 | } | ||
136 |