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 | #pragma once | ||
12 | |||
13 | #include <slepcsvd.h> | ||
14 | #include <slepc/private/slepcimpl.h> | ||
15 | |||
16 | /* SUBMANSEC = SVD */ | ||
17 | |||
18 | SLEPC_EXTERN PetscBool SVDRegisterAllCalled; | ||
19 | SLEPC_EXTERN PetscBool SVDMonitorRegisterAllCalled; | ||
20 | SLEPC_EXTERN PetscErrorCode SVDRegisterAll(void); | ||
21 | SLEPC_EXTERN PetscErrorCode SVDMonitorRegisterAll(void); | ||
22 | SLEPC_EXTERN PetscLogEvent SVD_SetUp,SVD_Solve; | ||
23 | |||
24 | typedef struct _SVDOps *SVDOps; | ||
25 | |||
26 | struct _SVDOps { | ||
27 | PetscErrorCode (*solve)(SVD); | ||
28 | PetscErrorCode (*solveg)(SVD); | ||
29 | PetscErrorCode (*solveh)(SVD); | ||
30 | PetscErrorCode (*setup)(SVD); | ||
31 | PetscErrorCode (*setfromoptions)(SVD,PetscOptionItems); | ||
32 | PetscErrorCode (*publishoptions)(SVD); | ||
33 | PetscErrorCode (*destroy)(SVD); | ||
34 | PetscErrorCode (*reset)(SVD); | ||
35 | PetscErrorCode (*view)(SVD,PetscViewer); | ||
36 | PetscErrorCode (*computevectors)(SVD); | ||
37 | PetscErrorCode (*setdstype)(SVD); | ||
38 | }; | ||
39 | |||
40 | /* | ||
41 | Maximum number of monitors you can run with a single SVD | ||
42 | */ | ||
43 | #define MAXSVDMONITORS 5 | ||
44 | |||
45 | typedef enum { SVD_STATE_INITIAL, | ||
46 | SVD_STATE_SETUP, | ||
47 | SVD_STATE_SOLVED, | ||
48 | SVD_STATE_VECTORS } SVDStateType; | ||
49 | |||
50 | /* | ||
51 | To check for unsupported features at SVDSetUp_XXX() | ||
52 | */ | ||
53 | typedef enum { SVD_FEATURE_CONVERGENCE=16, /* convergence test selected by user */ | ||
54 | SVD_FEATURE_STOPPING=32, /* stopping test */ | ||
55 | SVD_FEATURE_THRESHOLD=64 /* threshold stopping test */ | ||
56 | } SVDFeatureType; | ||
57 | |||
58 | /* | ||
59 | Defines the SVD data structure. | ||
60 | */ | ||
61 | struct _p_SVD { | ||
62 | PETSCHEADER(struct _SVDOps); | ||
63 | /*------------------------- User parameters ---------------------------*/ | ||
64 | Mat OP,OPb; /* problem matrices */ | ||
65 | Vec omega; /* signature for hyperbolic problems */ | ||
66 | PetscInt max_it; /* max iterations */ | ||
67 | PetscInt nsv; /* number of requested values */ | ||
68 | PetscInt ncv; /* basis size */ | ||
69 | PetscInt mpd; /* maximum dimension of projected problem */ | ||
70 | PetscInt nini,ninil; /* number of initial vecs (negative means not copied yet) */ | ||
71 | PetscReal tol; /* tolerance */ | ||
72 | PetscReal thres; /* threshold value */ | ||
73 | PetscBool threlative; /* threshold is relative */ | ||
74 | SVDConv conv; /* convergence test */ | ||
75 | SVDStop stop; /* stopping test */ | ||
76 | SVDWhich which; /* which singular values are computed */ | ||
77 | SVDProblemType problem_type; /* which kind of problem to be solved */ | ||
78 | PetscBool impltrans; /* implicit transpose mode */ | ||
79 | PetscBool trackall; /* whether all the residuals must be computed */ | ||
80 | |||
81 | /*-------------- User-provided functions and contexts -----------------*/ | ||
82 | SVDConvergenceTestFn *converged; | ||
83 | SVDConvergenceTestFn *convergeduser; | ||
84 | PetscCtxDestroyFn *convergeddestroy; | ||
85 | SVDStoppingTestFn *stopping; | ||
86 | SVDStoppingTestFn *stoppinguser; | ||
87 | PetscCtxDestroyFn *stoppingdestroy; | ||
88 | void *convergedctx; | ||
89 | void *stoppingctx; | ||
90 | SVDMonitorFn *monitor[MAXSVDMONITORS]; | ||
91 | PetscCtxDestroyFn *monitordestroy[MAXSVDMONITORS]; | ||
92 | void *monitorcontext[MAXSVDMONITORS]; | ||
93 | PetscInt numbermonitors; | ||
94 | |||
95 | /*----------------- Child objects and working data -------------------*/ | ||
96 | DS ds; /* direct solver object */ | ||
97 | BV U,V; /* left and right singular vectors */ | ||
98 | SlepcSC sc; /* sorting criterion data */ | ||
99 | Mat A,B; /* problem matrices */ | ||
100 | Mat AT,BT; /* transposed matrices */ | ||
101 | Vec *IS,*ISL; /* placeholder for references to user initial space */ | ||
102 | PetscReal *sigma; /* singular values */ | ||
103 | PetscReal *errest; /* error estimates */ | ||
104 | PetscReal *sign; /* +-1 for each singular value in hyperbolic problems=U'*Omega*U */ | ||
105 | PetscInt *perm; /* permutation for singular value ordering */ | ||
106 | PetscInt nworkl,nworkr; /* number of work vectors */ | ||
107 | Vec *workl,*workr; /* work vectors */ | ||
108 | void *data; /* placeholder for solver-specific stuff */ | ||
109 | |||
110 | /* ----------------------- Status variables -------------------------- */ | ||
111 | SVDStateType state; /* initial -> setup -> solved -> vectors */ | ||
112 | PetscInt nconv; /* number of converged values */ | ||
113 | PetscInt its; /* iteration counter */ | ||
114 | PetscBool leftbasis; /* if U is filled by the solver */ | ||
115 | PetscBool swapped; /* the U and V bases have been swapped (M<N) */ | ||
116 | PetscBool expltrans; /* explicit transpose created */ | ||
117 | PetscReal nrma,nrmb; /* computed matrix norms */ | ||
118 | PetscBool isgeneralized; | ||
119 | PetscBool ishyperbolic; | ||
120 | SVDConvergedReason reason; | ||
121 | }; | ||
122 | |||
123 | /* | ||
124 | Macros to test valid SVD arguments | ||
125 | */ | ||
126 | #if !defined(PETSC_USE_DEBUG) | ||
127 | |||
128 | #define SVDCheckSolved(h,arg) do {(void)(h);} while (0) | ||
129 | |||
130 | #else | ||
131 | |||
132 | #define SVDCheckSolved(h,arg) \ | ||
133 | do { \ | ||
134 | PetscCheck((h)->state>=SVD_STATE_SOLVED,PetscObjectComm((PetscObject)(h)),PETSC_ERR_ARG_WRONGSTATE,"Must call SVDSolve() first: Parameter #%d",arg); \ | ||
135 | } while (0) | ||
136 | |||
137 | #endif | ||
138 | |||
139 | /* | ||
140 | Macros to check settings at SVDSetUp() | ||
141 | */ | ||
142 | |||
143 | /* SVDCheckStandard: the problem is not GSVD */ | ||
144 | #define SVDCheckStandardCondition(svd,condition,msg) \ | ||
145 | do { \ | ||
146 | if (condition) { \ | ||
147 | PetscCheck(!(svd)->isgeneralized,PetscObjectComm((PetscObject)(svd)),PETSC_ERR_SUP,"The solver '%s'%s cannot be used for generalized problems",((PetscObject)(svd))->type_name,(msg)); \ | ||
148 | } \ | ||
149 | } while (0) | ||
150 | #define SVDCheckStandard(svd) SVDCheckStandardCondition(svd,PETSC_TRUE,"") | ||
151 | |||
152 | /* SVDCheckDefinite: the problem is not hyperbolic */ | ||
153 | #define SVDCheckDefiniteCondition(svd,condition,msg) \ | ||
154 | do { \ | ||
155 | if (condition) { \ | ||
156 | PetscCheck(!(svd)->ishyperbolic,PetscObjectComm((PetscObject)(svd)),PETSC_ERR_SUP,"The solver '%s'%s cannot be used for hyperbolic problems",((PetscObject)(svd))->type_name,(msg)); \ | ||
157 | } \ | ||
158 | } while (0) | ||
159 | #define SVDCheckDefinite(svd) SVDCheckDefiniteCondition(svd,PETSC_TRUE,"") | ||
160 | |||
161 | /* Check for unsupported features */ | ||
162 | #define SVDCheckUnsupportedCondition(svd,mask,condition,msg) \ | ||
163 | do { \ | ||
164 | if (condition) { \ | ||
165 | PetscCheck(!((mask) & SVD_FEATURE_CONVERGENCE) || (svd)->converged==SVDConvergedRelative,PetscObjectComm((PetscObject)(svd)),PETSC_ERR_SUP,"The solver '%s'%s only supports the default convergence test",((PetscObject)(svd))->type_name,(msg)); \ | ||
166 | PetscCheck(!((mask) & SVD_FEATURE_STOPPING) || (svd)->stopping==SVDStoppingBasic,PetscObjectComm((PetscObject)(svd)),PETSC_ERR_SUP,"The solver '%s'%s only supports the default stopping test",((PetscObject)(svd))->type_name,(msg)); \ | ||
167 | PetscCheck(!((mask) & SVD_FEATURE_THRESHOLD) || (svd)->stopping!=SVDStoppingThreshold,PetscObjectComm((PetscObject)(svd)),PETSC_ERR_SUP,"The solver '%s'%s does not support the threshold stopping test",((PetscObject)(svd))->type_name,(msg)); \ | ||
168 | } \ | ||
169 | } while (0) | ||
170 | #define SVDCheckUnsupported(svd,mask) SVDCheckUnsupportedCondition(svd,mask,PETSC_TRUE,"") | ||
171 | |||
172 | /* Check for ignored features */ | ||
173 | #define SVDCheckIgnoredCondition(svd,mask,condition,msg) \ | ||
174 | do { \ | ||
175 | if (condition) { \ | ||
176 | if (((mask) & SVD_FEATURE_CONVERGENCE) && (svd)->converged!=SVDConvergedRelative) PetscCall(PetscInfo((svd),"The solver '%s'%s ignores the convergence test settings\n",((PetscObject)(svd))->type_name,(msg))); \ | ||
177 | if (((mask) & SVD_FEATURE_STOPPING) && (svd)->stopping!=SVDStoppingBasic) PetscCall(PetscInfo((svd),"The solver '%s'%s ignores the stopping test settings\n",((PetscObject)(svd))->type_name,(msg))); \ | ||
178 | } \ | ||
179 | } while (0) | ||
180 | #define SVDCheckIgnored(svd,mask) SVDCheckIgnoredCondition(svd,mask,PETSC_TRUE,"") | ||
181 | |||
182 | /* | ||
183 | SVDSetCtxThreshold - Fills SVDStoppingCtx with data needed for the threshold stopping test | ||
184 | */ | ||
185 | #define SVDSetCtxThreshold(svd,sigma,k) \ | ||
186 | do { \ | ||
187 | if (svd->stop==SVD_STOP_THRESHOLD && k) { \ | ||
188 | ((SVDStoppingCtx)svd->stoppingctx)->firstsv = sigma[0]; \ | ||
189 | ((SVDStoppingCtx)svd->stoppingctx)->lastsv = sigma[k-1]; \ | ||
190 | } \ | ||
191 | } while (0) | ||
192 | |||
193 | /* | ||
194 | SVD_KSPSetOperators - Sets the KSP matrices | ||
195 | */ | ||
196 | 729 | static inline PetscErrorCode SVD_KSPSetOperators(KSP ksp,Mat A,Mat B) | |
197 | { | ||
198 | 729 | const char *prefix; | |
199 | |||
200 |
1/2✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
|
729 | PetscFunctionBegin; |
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.
|
729 | PetscCall(KSPSetOperators(ksp,A,B)); |
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.
|
729 | PetscCall(MatGetOptionsPrefix(B,&prefix)); |
203 |
1/2✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
|
729 | if (!prefix) { |
204 | /* set Mat prefix to be the same as KSP to enable setting command-line options (e.g. MUMPS) | ||
205 | only applies if the Mat has no user-defined prefix */ | ||
206 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
729 | PetscCall(KSPGetOptionsPrefix(ksp,&prefix)); |
207 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
729 | PetscCall(MatSetOptionsPrefix(B,prefix)); |
208 | } | ||
209 |
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.
|
147 | PetscFunctionReturn(PETSC_SUCCESS); |
210 | } | ||
211 | |||
212 | /* | ||
213 | Create the template vector for the left basis in GSVD, as in | ||
214 | MatCreateVecsEmpty(Z,NULL,&t) for Z=[A;B] without forming Z. | ||
215 | */ | ||
216 | 1702 | static inline PetscErrorCode SVDCreateLeftTemplate(SVD svd,Vec *t) | |
217 | { | ||
218 | 1702 | PetscInt M,P,m,p; | |
219 | 1702 | Vec v1,v2; | |
220 | 1702 | VecType vec_type; | |
221 | |||
222 |
1/2✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
|
1702 | PetscFunctionBegin; |
223 |
4/6✓ Branch 0 taken 4 times.
✓ Branch 1 taken 16 times.
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 4 times.
|
1702 | PetscCall(MatCreateVecsEmpty(svd->OP,NULL,&v1)); |
224 |
4/6✓ Branch 0 taken 4 times.
✓ Branch 1 taken 16 times.
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 4 times.
|
1702 | PetscCall(VecGetSize(v1,&M)); |
225 |
4/6✓ Branch 0 taken 4 times.
✓ Branch 1 taken 16 times.
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 4 times.
|
1702 | PetscCall(VecGetLocalSize(v1,&m)); |
226 |
4/6✓ Branch 0 taken 4 times.
✓ Branch 1 taken 16 times.
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 4 times.
|
1702 | PetscCall(VecGetType(v1,&vec_type)); |
227 |
4/6✓ Branch 0 taken 4 times.
✓ Branch 1 taken 16 times.
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 4 times.
|
1702 | PetscCall(MatCreateVecsEmpty(svd->OPb,NULL,&v2)); |
228 |
4/6✓ Branch 0 taken 4 times.
✓ Branch 1 taken 16 times.
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 4 times.
|
1702 | PetscCall(VecGetSize(v2,&P)); |
229 |
4/6✓ Branch 0 taken 4 times.
✓ Branch 1 taken 16 times.
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 4 times.
|
1702 | PetscCall(VecGetLocalSize(v2,&p)); |
230 |
4/6✓ Branch 0 taken 4 times.
✓ Branch 1 taken 16 times.
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 4 times.
|
1702 | PetscCall(VecCreate(PetscObjectComm((PetscObject)(v1)),t)); |
231 |
4/6✓ Branch 0 taken 4 times.
✓ Branch 1 taken 16 times.
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 4 times.
|
1702 | PetscCall(VecSetType(*t,vec_type)); |
232 |
4/6✓ Branch 0 taken 4 times.
✓ Branch 1 taken 16 times.
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 4 times.
|
1702 | PetscCall(VecSetSizes(*t,m+p,M+P)); |
233 |
4/6✓ Branch 0 taken 4 times.
✓ Branch 1 taken 16 times.
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 4 times.
|
1702 | PetscCall(VecSetUp(*t)); |
234 |
4/6✓ Branch 0 taken 4 times.
✓ Branch 1 taken 16 times.
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 4 times.
|
1702 | PetscCall(VecDestroy(&v1)); |
235 |
4/6✓ Branch 0 taken 4 times.
✓ Branch 1 taken 16 times.
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 4 times.
|
1702 | PetscCall(VecDestroy(&v2)); |
236 |
6/12✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 4 times.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 4 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 4 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 4 times.
|
342 | PetscFunctionReturn(PETSC_SUCCESS); |
237 | } | ||
238 | |||
239 | SLEPC_INTERN PetscErrorCode SVDKrylovConvergence(SVD,PetscBool,PetscInt,PetscInt,PetscReal,PetscInt*); | ||
240 | SLEPC_INTERN PetscErrorCode SVDTwoSideLanczos(SVD,PetscReal*,PetscReal*,BV,BV,PetscInt,PetscInt*,PetscBool*); | ||
241 | SLEPC_INTERN PetscErrorCode SVDSetDimensions_Default(SVD); | ||
242 | SLEPC_INTERN PetscErrorCode SVDComputeVectors(SVD); | ||
243 | SLEPC_INTERN PetscErrorCode SVDComputeVectors_Left(SVD); | ||
244 |