Line data Source code
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 : } SVDFeatureType;
56 :
57 : /*
58 : Defines the SVD data structure.
59 : */
60 : struct _p_SVD {
61 : PETSCHEADER(struct _SVDOps);
62 : /*------------------------- User parameters ---------------------------*/
63 : Mat OP,OPb; /* problem matrices */
64 : Vec omega; /* signature for hyperbolic problems */
65 : PetscInt max_it; /* max iterations */
66 : PetscInt nsv; /* number of requested values */
67 : PetscInt ncv; /* basis size */
68 : PetscInt mpd; /* maximum dimension of projected problem */
69 : PetscInt nini,ninil; /* number of initial vecs (negative means not copied yet) */
70 : PetscReal tol; /* tolerance */
71 : SVDConv conv; /* convergence test */
72 : SVDStop stop; /* stopping test */
73 : SVDWhich which; /* which singular values are computed */
74 : SVDProblemType problem_type; /* which kind of problem to be solved */
75 : PetscBool impltrans; /* implicit transpose mode */
76 : PetscBool trackall; /* whether all the residuals must be computed */
77 :
78 : /*-------------- User-provided functions and contexts -----------------*/
79 : PetscErrorCode (*converged)(SVD,PetscReal,PetscReal,PetscReal*,void*);
80 : PetscErrorCode (*convergeduser)(SVD,PetscReal,PetscReal,PetscReal*,void*);
81 : PetscErrorCode (*convergeddestroy)(void*);
82 : PetscErrorCode (*stopping)(SVD,PetscInt,PetscInt,PetscInt,PetscInt,SVDConvergedReason*,void*);
83 : PetscErrorCode (*stoppinguser)(SVD,PetscInt,PetscInt,PetscInt,PetscInt,SVDConvergedReason*,void*);
84 : PetscErrorCode (*stoppingdestroy)(void*);
85 : void *convergedctx;
86 : void *stoppingctx;
87 : PetscErrorCode (*monitor[MAXSVDMONITORS])(SVD,PetscInt,PetscInt,PetscReal*,PetscReal*,PetscInt,void*);
88 : PetscErrorCode (*monitordestroy[MAXSVDMONITORS])(void**);
89 : void *monitorcontext[MAXSVDMONITORS];
90 : PetscInt numbermonitors;
91 :
92 : /*----------------- Child objects and working data -------------------*/
93 : DS ds; /* direct solver object */
94 : BV U,V; /* left and right singular vectors */
95 : SlepcSC sc; /* sorting criterion data */
96 : Mat A,B; /* problem matrices */
97 : Mat AT,BT; /* transposed matrices */
98 : Vec *IS,*ISL; /* placeholder for references to user initial space */
99 : PetscReal *sigma; /* singular values */
100 : PetscReal *errest; /* error estimates */
101 : PetscReal *sign; /* +-1 for each singular value in hyperbolic problems=U'*Omega*U */
102 : PetscInt *perm; /* permutation for singular value ordering */
103 : PetscInt nworkl,nworkr; /* number of work vectors */
104 : Vec *workl,*workr; /* work vectors */
105 : void *data; /* placeholder for solver-specific stuff */
106 :
107 : /* ----------------------- Status variables -------------------------- */
108 : SVDStateType state; /* initial -> setup -> solved -> vectors */
109 : PetscInt nconv; /* number of converged values */
110 : PetscInt its; /* iteration counter */
111 : PetscBool leftbasis; /* if U is filled by the solver */
112 : PetscBool swapped; /* the U and V bases have been swapped (M<N) */
113 : PetscBool expltrans; /* explicit transpose created */
114 : PetscReal nrma,nrmb; /* computed matrix norms */
115 : PetscBool isgeneralized;
116 : PetscBool ishyperbolic;
117 : SVDConvergedReason reason;
118 : };
119 :
120 : /*
121 : Macros to test valid SVD arguments
122 : */
123 : #if !defined(PETSC_USE_DEBUG)
124 :
125 : #define SVDCheckSolved(h,arg) do {(void)(h);} while (0)
126 :
127 : #else
128 :
129 : #define SVDCheckSolved(h,arg) \
130 : do { \
131 : PetscCheck((h)->state>=SVD_STATE_SOLVED,PetscObjectComm((PetscObject)(h)),PETSC_ERR_ARG_WRONGSTATE,"Must call SVDSolve() first: Parameter #%d",arg); \
132 : } while (0)
133 :
134 : #endif
135 :
136 : /*
137 : Macros to check settings at SVDSetUp()
138 : */
139 :
140 : /* SVDCheckStandard: the problem is not GSVD */
141 : #define SVDCheckStandardCondition(svd,condition,msg) \
142 : do { \
143 : if (condition) { \
144 : PetscCheck(!(svd)->isgeneralized,PetscObjectComm((PetscObject)(svd)),PETSC_ERR_SUP,"The solver '%s'%s cannot be used for generalized problems",((PetscObject)(svd))->type_name,(msg)); \
145 : } \
146 : } while (0)
147 : #define SVDCheckStandard(svd) SVDCheckStandardCondition(svd,PETSC_TRUE,"")
148 :
149 : /* SVDCheckDefinite: the problem is not hyperbolic */
150 : #define SVDCheckDefiniteCondition(svd,condition,msg) \
151 : do { \
152 : if (condition) { \
153 : PetscCheck(!(svd)->ishyperbolic,PetscObjectComm((PetscObject)(svd)),PETSC_ERR_SUP,"The solver '%s'%s cannot be used for hyperbolic problems",((PetscObject)(svd))->type_name,(msg)); \
154 : } \
155 : } while (0)
156 : #define SVDCheckDefinite(svd) SVDCheckDefiniteCondition(svd,PETSC_TRUE,"")
157 :
158 : /* Check for unsupported features */
159 : #define SVDCheckUnsupportedCondition(svd,mask,condition,msg) \
160 : do { \
161 : if (condition) { \
162 : 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)); \
163 : 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)); \
164 : } \
165 : } while (0)
166 : #define SVDCheckUnsupported(svd,mask) SVDCheckUnsupportedCondition(svd,mask,PETSC_TRUE,"")
167 :
168 : /* Check for ignored features */
169 : #define SVDCheckIgnoredCondition(svd,mask,condition,msg) \
170 : do { \
171 : if (condition) { \
172 : 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))); \
173 : 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))); \
174 : } \
175 : } while (0)
176 : #define SVDCheckIgnored(svd,mask) SVDCheckIgnoredCondition(svd,mask,PETSC_TRUE,"")
177 :
178 : /*
179 : SVD_KSPSetOperators - Sets the KSP matrices
180 : */
181 78 : static inline PetscErrorCode SVD_KSPSetOperators(KSP ksp,Mat A,Mat B)
182 : {
183 78 : const char *prefix;
184 :
185 78 : PetscFunctionBegin;
186 78 : PetscCall(KSPSetOperators(ksp,A,B));
187 78 : PetscCall(MatGetOptionsPrefix(B,&prefix));
188 78 : if (!prefix) {
189 : /* set Mat prefix to be the same as KSP to enable setting command-line options (e.g. MUMPS)
190 : only applies if the Mat has no user-defined prefix */
191 78 : PetscCall(KSPGetOptionsPrefix(ksp,&prefix));
192 78 : PetscCall(MatSetOptionsPrefix(B,prefix));
193 : }
194 78 : PetscFunctionReturn(PETSC_SUCCESS);
195 : }
196 :
197 : /*
198 : Create the template vector for the left basis in GSVD, as in
199 : MatCreateVecsEmpty(Z,NULL,&t) for Z=[A;B] without forming Z.
200 : */
201 174 : static inline PetscErrorCode SVDCreateLeftTemplate(SVD svd,Vec *t)
202 : {
203 174 : PetscInt M,P,m,p;
204 174 : Vec v1,v2;
205 174 : VecType vec_type;
206 :
207 174 : PetscFunctionBegin;
208 174 : PetscCall(MatCreateVecsEmpty(svd->OP,NULL,&v1));
209 174 : PetscCall(VecGetSize(v1,&M));
210 174 : PetscCall(VecGetLocalSize(v1,&m));
211 174 : PetscCall(VecGetType(v1,&vec_type));
212 174 : PetscCall(MatCreateVecsEmpty(svd->OPb,NULL,&v2));
213 174 : PetscCall(VecGetSize(v2,&P));
214 174 : PetscCall(VecGetLocalSize(v2,&p));
215 174 : PetscCall(VecCreate(PetscObjectComm((PetscObject)(v1)),t));
216 174 : PetscCall(VecSetType(*t,vec_type));
217 174 : PetscCall(VecSetSizes(*t,m+p,M+P));
218 174 : PetscCall(VecSetUp(*t));
219 174 : PetscCall(VecDestroy(&v1));
220 174 : PetscCall(VecDestroy(&v2));
221 174 : PetscFunctionReturn(PETSC_SUCCESS);
222 : }
223 :
224 : SLEPC_INTERN PetscErrorCode SVDKrylovConvergence(SVD,PetscBool,PetscInt,PetscInt,PetscReal,PetscInt*);
225 : SLEPC_INTERN PetscErrorCode SVDTwoSideLanczos(SVD,PetscReal*,PetscReal*,BV,BV,PetscInt,PetscInt*,PetscBool*);
226 : SLEPC_INTERN PetscErrorCode SVDSetDimensions_Default(SVD);
227 : SLEPC_INTERN PetscErrorCode SVDComputeVectors(SVD);
228 : SLEPC_INTERN PetscErrorCode SVDComputeVectors_Left(SVD);
|
