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 : SVD routines for setting up the solver
12 : */
13 :
14 : #include <slepc/private/svdimpl.h> /*I "slepcsvd.h" I*/
15 :
16 : /*@
17 : SVDSetOperators - Set the matrices associated with the singular value problem.
18 :
19 : Collective
20 :
21 : Input Parameters:
22 : + svd - the singular value solver context
23 : . A - the matrix associated with the singular value problem
24 : - B - the second matrix in the case of GSVD
25 :
26 : Level: beginner
27 :
28 : .seealso: SVDSolve(), SVDGetOperators()
29 : @*/
30 270 : PetscErrorCode SVDSetOperators(SVD svd,Mat A,Mat B)
31 : {
32 270 : PetscInt Ma,Na,Mb,Nb,ma,na,mb,nb,M0,N0,m0,n0;
33 270 : PetscBool samesize=PETSC_TRUE;
34 :
35 270 : PetscFunctionBegin;
36 270 : PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
37 270 : PetscValidHeaderSpecific(A,MAT_CLASSID,2);
38 270 : if (B) PetscValidHeaderSpecific(B,MAT_CLASSID,3);
39 270 : PetscCheckSameComm(svd,1,A,2);
40 270 : if (B) PetscCheckSameComm(svd,1,B,3);
41 :
42 : /* Check matrix sizes */
43 270 : PetscCall(MatGetSize(A,&Ma,&Na));
44 270 : PetscCall(MatGetLocalSize(A,&ma,&na));
45 270 : if (svd->OP) {
46 25 : PetscCall(MatGetSize(svd->OP,&M0,&N0));
47 25 : PetscCall(MatGetLocalSize(svd->OP,&m0,&n0));
48 25 : if (M0!=Ma || N0!=Na || m0!=ma || n0!=na) samesize = PETSC_FALSE;
49 : }
50 270 : if (B) {
51 102 : PetscCall(MatGetSize(B,&Mb,&Nb));
52 102 : PetscCall(MatGetLocalSize(B,&mb,&nb));
53 102 : PetscCheck(Na==Nb,PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_WRONG,"Different number of columns in A (%" PetscInt_FMT ") and B (%" PetscInt_FMT ")",Na,Nb);
54 102 : PetscCheck(na==nb,PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_WRONG,"Different local column size in A (%" PetscInt_FMT ") and B (%" PetscInt_FMT ")",na,nb);
55 102 : if (svd->OPb) {
56 11 : PetscCall(MatGetSize(svd->OPb,&M0,&N0));
57 11 : PetscCall(MatGetLocalSize(svd->OPb,&m0,&n0));
58 11 : if (M0!=Mb || N0!=Nb || m0!=mb || n0!=nb) samesize = PETSC_FALSE;
59 : }
60 : }
61 :
62 270 : PetscCall(PetscObjectReference((PetscObject)A));
63 270 : if (B) PetscCall(PetscObjectReference((PetscObject)B));
64 270 : if (svd->state && !samesize) PetscCall(SVDReset(svd));
65 : else {
66 270 : PetscCall(MatDestroy(&svd->OP));
67 270 : PetscCall(MatDestroy(&svd->OPb));
68 270 : PetscCall(MatDestroy(&svd->A));
69 270 : PetscCall(MatDestroy(&svd->B));
70 270 : PetscCall(MatDestroy(&svd->AT));
71 270 : PetscCall(MatDestroy(&svd->BT));
72 : }
73 270 : svd->nrma = 0.0;
74 270 : svd->nrmb = 0.0;
75 270 : svd->OP = A;
76 270 : svd->OPb = B;
77 270 : svd->state = SVD_STATE_INITIAL;
78 270 : PetscFunctionReturn(PETSC_SUCCESS);
79 : }
80 :
81 : /*@
82 : SVDGetOperators - Get the matrices associated with the singular value problem.
83 :
84 : Not Collective
85 :
86 : Input Parameter:
87 : . svd - the singular value solver context
88 :
89 : Output Parameters:
90 : + A - the matrix associated with the singular value problem
91 : - B - the second matrix in the case of GSVD
92 :
93 : Level: intermediate
94 :
95 : .seealso: SVDSolve(), SVDSetOperators()
96 : @*/
97 12 : PetscErrorCode SVDGetOperators(SVD svd,Mat *A,Mat *B)
98 : {
99 12 : PetscFunctionBegin;
100 12 : PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
101 12 : if (A) *A = svd->OP;
102 12 : if (B) *B = svd->OPb;
103 12 : PetscFunctionReturn(PETSC_SUCCESS);
104 : }
105 :
106 : /*@
107 : SVDSetSignature - Set the signature matrix defining a hyperbolic singular value problem.
108 :
109 : Collective
110 :
111 : Input Parameters:
112 : + svd - the singular value solver context
113 : - omega - a vector containing the diagonal elements of the signature matrix (or NULL)
114 :
115 : Notes:
116 : The signature matrix is relevant only for hyperbolic problems (HSVD).
117 : Use NULL to reset a previously set signature.
118 :
119 : Level: intermediate
120 :
121 : .seealso: SVDSetProblemType(), SVDSetOperators(), SVDGetSignature()
122 : @*/
123 38 : PetscErrorCode SVDSetSignature(SVD svd,Vec omega)
124 : {
125 38 : PetscInt N,Ma,n,ma;
126 :
127 38 : PetscFunctionBegin;
128 38 : PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
129 38 : if (omega) {
130 38 : PetscValidHeaderSpecific(omega,VEC_CLASSID,2);
131 38 : PetscCheckSameComm(svd,1,omega,2);
132 : }
133 :
134 38 : if (omega && svd->OP) { /* Check sizes */
135 37 : PetscCall(VecGetSize(omega,&N));
136 37 : PetscCall(VecGetLocalSize(omega,&n));
137 37 : PetscCall(MatGetSize(svd->OP,&Ma,NULL));
138 37 : PetscCall(MatGetLocalSize(svd->OP,&ma,NULL));
139 37 : PetscCheck(N==Ma,PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_SIZ,"Global size of signature (%" PetscInt_FMT ") does not match the row size of A (%" PetscInt_FMT ")",N,Ma);
140 37 : PetscCheck(n==ma,PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_SIZ,"Local size of signature (%" PetscInt_FMT ") does not match the local row size of A (%" PetscInt_FMT ")",n,ma);
141 : }
142 :
143 38 : PetscCall(VecDestroy(&svd->omega));
144 38 : if (omega) {
145 38 : PetscCall(VecDuplicate(omega,&svd->omega));
146 38 : PetscCall(VecCopy(omega,svd->omega));
147 : }
148 38 : svd->state = SVD_STATE_INITIAL;
149 38 : PetscFunctionReturn(PETSC_SUCCESS);
150 : }
151 :
152 : /*@
153 : SVDGetSignature - Get the signature matrix defining a hyperbolic singular value problem.
154 :
155 : Not Collective
156 :
157 : Input Parameter:
158 : . svd - the singular value solver context
159 :
160 : Output Parameter:
161 : . omega - a vector containing the diagonal elements of the signature matrix
162 :
163 : Notes:
164 : The signature matrix is relevant only for hyperbolic problems (HSVD).
165 : If no signature has been set, this function will return a vector of all ones.
166 :
167 : The user should pass a previously created Vec with the appropriate dimension.
168 :
169 : Level: intermediate
170 :
171 : .seealso: SVDSetSignature()
172 : @*/
173 1 : PetscErrorCode SVDGetSignature(SVD svd,Vec omega)
174 : {
175 1 : PetscFunctionBegin;
176 1 : PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
177 1 : PetscValidHeaderSpecific(omega,VEC_CLASSID,2);
178 1 : if (svd->omega) PetscCall(VecCopy(svd->omega,omega));
179 0 : else PetscCall(VecSet(omega,1.0));
180 1 : PetscFunctionReturn(PETSC_SUCCESS);
181 : }
182 :
183 : /*@
184 : SVDSetDSType - Sets the type of the internal DS object based on the current
185 : settings of the singular value solver.
186 :
187 : Collective
188 :
189 : Input Parameter:
190 : . svd - singular value solver context
191 :
192 : Note:
193 : This function need not be called explicitly, since it will be called at
194 : both SVDSetFromOptions() and SVDSetUp().
195 :
196 : Level: developer
197 :
198 : .seealso: SVDSetFromOptions(), SVDSetUp()
199 : @*/
200 517 : PetscErrorCode SVDSetDSType(SVD svd)
201 : {
202 517 : PetscFunctionBegin;
203 517 : PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
204 517 : PetscTryTypeMethod(svd,setdstype);
205 517 : PetscFunctionReturn(PETSC_SUCCESS);
206 : }
207 :
208 : /*@
209 : SVDSetUp - Sets up all the internal data structures necessary for the
210 : execution of the singular value solver.
211 :
212 : Collective
213 :
214 : Input Parameter:
215 : . svd - singular value solver context
216 :
217 : Notes:
218 : This function need not be called explicitly in most cases, since SVDSolve()
219 : calls it. It can be useful when one wants to measure the set-up time
220 : separately from the solve time.
221 :
222 : Level: developer
223 :
224 : .seealso: SVDCreate(), SVDSolve(), SVDDestroy()
225 : @*/
226 286 : PetscErrorCode SVDSetUp(SVD svd)
227 : {
228 286 : PetscBool flg;
229 286 : PetscInt M,N,P=0,k,maxnsol,m,Nom,nom;
230 286 : SlepcSC sc;
231 286 : Vec *T;
232 286 : BV bv;
233 286 : SVDStoppingCtx ctx;
234 :
235 286 : PetscFunctionBegin;
236 286 : PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
237 286 : if (svd->state) PetscFunctionReturn(PETSC_SUCCESS);
238 286 : PetscCall(PetscLogEventBegin(SVD_SetUp,svd,0,0,0));
239 :
240 : /* reset the convergence flag from the previous solves */
241 286 : svd->reason = SVD_CONVERGED_ITERATING;
242 :
243 : /* set default solver type (SVDSetFromOptions was not called) */
244 286 : if (!((PetscObject)svd)->type_name) PetscCall(SVDSetType(svd,SVDCROSS));
245 286 : if (!svd->ds) PetscCall(SVDGetDS(svd,&svd->ds));
246 :
247 : /* check matrices */
248 286 : PetscCheck(svd->OP,PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_WRONGSTATE,"SVDSetOperators() must be called first");
249 :
250 286 : if (!svd->problem_type) { /* set default problem type */
251 124 : if (svd->OPb) {
252 18 : PetscCheck(!svd->omega,PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"There is no support yet for generalized hyperbolic problems");
253 18 : PetscCall(SVDSetProblemType(svd,SVD_GENERALIZED));
254 : } else {
255 106 : if (svd->omega) PetscCall(SVDSetProblemType(svd,SVD_HYPERBOLIC));
256 105 : else PetscCall(SVDSetProblemType(svd,SVD_STANDARD));
257 : }
258 : } else { /* check consistency of problem type set by user */
259 162 : if (svd->OPb) {
260 84 : PetscCheck(svd->isgeneralized,PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_INCOMP,"Inconsistent SVD state: the problem type does not match the number of matrices");
261 84 : PetscCheck(!svd->omega,PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"There is no support yet for generalized hyperbolic problems");
262 : } else {
263 78 : PetscCheck(!svd->isgeneralized,PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_INCOMP,"Inconsistent SVD state: the problem type does not match the number of matrices");
264 78 : if (svd->omega) PetscCheck(svd->ishyperbolic,PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_INCOMP,"Inconsistent SVD state: the problem type must be set to hyperbolic when passing a signature with SVDSetSignature()");
265 41 : else PetscCheck(!svd->ishyperbolic,PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_INCOMP,"Inconsistent SVD state: a hyperbolic problem requires passing a signature with SVDSetSignature()");
266 : }
267 : }
268 :
269 : /* set DS type once the default problem type has been determined */
270 286 : PetscCall(SVDSetDSType(svd));
271 :
272 : /* determine how to handle the transpose */
273 286 : svd->expltrans = PETSC_TRUE;
274 286 : if (svd->impltrans) svd->expltrans = PETSC_FALSE;
275 : else {
276 267 : PetscCall(MatHasOperation(svd->OP,MATOP_TRANSPOSE,&flg));
277 267 : if (!flg) svd->expltrans = PETSC_FALSE;
278 : else {
279 266 : PetscCall(PetscObjectTypeCompareAny((PetscObject)svd,&flg,SVDLAPACK,SVDSCALAPACK,SVDKSVD,SVDELEMENTAL,""));
280 266 : if (flg) svd->expltrans = PETSC_FALSE;
281 : }
282 : }
283 :
284 : /* get matrix dimensions */
285 286 : PetscCall(MatGetSize(svd->OP,&M,&N));
286 286 : if (svd->isgeneralized) {
287 102 : PetscCall(MatGetSize(svd->OPb,&P,NULL));
288 102 : PetscCheck(M+P>=N,PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"The case when [A;B] has less rows than columns is not supported");
289 184 : } else if (svd->ishyperbolic) {
290 38 : PetscCall(VecGetSize(svd->omega,&Nom));
291 38 : PetscCall(VecGetLocalSize(svd->omega,&nom));
292 38 : PetscCall(MatGetLocalSize(svd->OP,&m,NULL));
293 38 : PetscCheck(Nom==M,PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_SIZ,"Global size of signature (%" PetscInt_FMT ") does not match the row size of A (%" PetscInt_FMT ")",Nom,M);
294 38 : PetscCheck(nom==m,PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_SIZ,"Local size of signature (%" PetscInt_FMT ") does not match the local row size of A (%" PetscInt_FMT ")",nom,m);
295 : }
296 :
297 : /* build transpose matrix */
298 286 : PetscCall(MatDestroy(&svd->A));
299 286 : PetscCall(MatDestroy(&svd->AT));
300 286 : PetscCall(PetscObjectReference((PetscObject)svd->OP));
301 286 : if (svd->expltrans) {
302 242 : if (svd->isgeneralized || M>=N) {
303 207 : svd->A = svd->OP;
304 207 : PetscCall(MatHermitianTranspose(svd->OP,MAT_INITIAL_MATRIX,&svd->AT));
305 : } else {
306 35 : PetscCall(MatHermitianTranspose(svd->OP,MAT_INITIAL_MATRIX,&svd->A));
307 35 : svd->AT = svd->OP;
308 : }
309 : } else {
310 44 : if (svd->isgeneralized || M>=N) {
311 35 : svd->A = svd->OP;
312 35 : PetscCall(MatCreateHermitianTranspose(svd->OP,&svd->AT));
313 : } else {
314 9 : PetscCall(MatCreateHermitianTranspose(svd->OP,&svd->A));
315 9 : svd->AT = svd->OP;
316 : }
317 : }
318 :
319 : /* build transpose matrix B for GSVD */
320 286 : if (svd->isgeneralized) {
321 102 : PetscCall(MatDestroy(&svd->B));
322 102 : PetscCall(MatDestroy(&svd->BT));
323 102 : PetscCall(PetscObjectReference((PetscObject)svd->OPb));
324 102 : if (svd->expltrans) {
325 95 : svd->B = svd->OPb;
326 95 : PetscCall(MatHermitianTranspose(svd->OPb,MAT_INITIAL_MATRIX,&svd->BT));
327 : } else {
328 7 : svd->B = svd->OPb;
329 7 : PetscCall(MatCreateHermitianTranspose(svd->OPb,&svd->BT));
330 : }
331 : }
332 :
333 286 : if (!svd->isgeneralized && M<N) {
334 : /* swap initial vectors */
335 44 : if (svd->nini || svd->ninil) {
336 0 : T=svd->ISL; svd->ISL=svd->IS; svd->IS=T;
337 0 : k=svd->ninil; svd->ninil=svd->nini; svd->nini=k;
338 : }
339 : /* swap basis vectors */
340 44 : if (!svd->swapped) { /* only the first time in case of multiple calls */
341 37 : bv=svd->V; svd->V=svd->U; svd->U=bv;
342 37 : svd->swapped = PETSC_TRUE;
343 : }
344 : }
345 :
346 286 : maxnsol = svd->isgeneralized? PetscMin(PetscMin(M,N),P): PetscMin(M,N);
347 286 : svd->ncv = PetscMin(svd->ncv,maxnsol);
348 286 : svd->nsv = PetscMin(svd->nsv,maxnsol);
349 286 : PetscCheck(svd->ncv==PETSC_DETERMINE || svd->nsv<=svd->ncv,PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_OUTOFRANGE,"nsv bigger than ncv");
350 :
351 : /* relative convergence criterion is not allowed in GSVD */
352 419 : if (svd->conv==(SVDConv)-1) PetscCall(SVDSetConvergenceTest(svd,svd->isgeneralized?SVD_CONV_NORM:SVD_CONV_REL));
353 286 : PetscCheck(!svd->isgeneralized || svd->conv!=SVD_CONV_REL,PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"Relative convergence criterion is not allowed in GSVD");
354 :
355 : /* initialization of matrix norm (standard case only, for GSVD it is done inside setup()) */
356 286 : if (!svd->isgeneralized && svd->conv==SVD_CONV_NORM && !svd->nrma) PetscCall(MatNorm(svd->OP,NORM_INFINITY,&svd->nrma));
357 :
358 : /* threshold stopping test */
359 286 : if (svd->stop==SVD_STOP_THRESHOLD) {
360 12 : PetscCheck(svd->thres>0.0,PetscObjectComm((PetscObject)svd),PETSC_ERR_USER_INPUT,"Must give a threshold value with SVDSetThreshold()");
361 12 : PetscCheck(svd->which==SVD_LARGEST || !svd->threlative,PetscObjectComm((PetscObject)svd),PETSC_ERR_USER_INPUT,"Can only use a relative threshold when which=SVD_LARGEST");
362 12 : PetscCall(PetscNew(&ctx));
363 12 : ctx->thres = svd->thres;
364 12 : ctx->threlative = svd->threlative;
365 12 : ctx->which = svd->which;
366 12 : PetscCall(SVDSetStoppingTestFunction(svd,SVDStoppingThreshold,ctx,PetscCtxDestroyDefault));
367 : }
368 :
369 : /* call specific solver setup */
370 286 : PetscUseTypeMethod(svd,setup);
371 :
372 : /* set tolerance if not yet set */
373 286 : if (svd->tol==(PetscReal)PETSC_DETERMINE) svd->tol = SLEPC_DEFAULT_TOL;
374 :
375 : /* fill sorting criterion context */
376 286 : PetscCall(DSGetSlepcSC(svd->ds,&sc));
377 286 : sc->comparison = (svd->which==SVD_LARGEST)? SlepcCompareLargestReal: SlepcCompareSmallestReal;
378 286 : sc->comparisonctx = NULL;
379 286 : sc->map = NULL;
380 286 : sc->mapobj = NULL;
381 :
382 : /* process initial vectors */
383 286 : if (svd->nini<0) {
384 25 : k = -svd->nini;
385 25 : PetscCheck(k<=svd->ncv,PetscObjectComm((PetscObject)svd),PETSC_ERR_USER_INPUT,"The number of initial vectors is larger than ncv");
386 25 : PetscCall(BVInsertVecs(svd->V,0,&k,svd->IS,PETSC_TRUE));
387 25 : PetscCall(SlepcBasisDestroy_Private(&svd->nini,&svd->IS));
388 25 : svd->nini = k;
389 : }
390 286 : if (svd->ninil<0) {
391 30 : k = 0;
392 30 : if (svd->leftbasis) {
393 23 : k = -svd->ninil;
394 23 : PetscCheck(k<=svd->ncv,PetscObjectComm((PetscObject)svd),PETSC_ERR_USER_INPUT,"The number of left initial vectors is larger than ncv");
395 23 : PetscCall(BVInsertVecs(svd->U,0,&k,svd->ISL,PETSC_TRUE));
396 7 : } else PetscCall(PetscInfo(svd,"Ignoring initial left vectors\n"));
397 30 : PetscCall(SlepcBasisDestroy_Private(&svd->ninil,&svd->ISL));
398 30 : svd->ninil = k;
399 : }
400 :
401 286 : PetscCall(PetscLogEventEnd(SVD_SetUp,svd,0,0,0));
402 286 : svd->state = SVD_STATE_SETUP;
403 286 : PetscFunctionReturn(PETSC_SUCCESS);
404 : }
405 :
406 : /*@
407 : SVDSetInitialSpaces - Specify two basis of vectors that constitute the initial
408 : right and/or left spaces.
409 :
410 : Collective
411 :
412 : Input Parameters:
413 : + svd - the singular value solver context
414 : . nr - number of right vectors
415 : . isr - set of basis vectors of the right initial space
416 : . nl - number of left vectors
417 : - isl - set of basis vectors of the left initial space
418 :
419 : Notes:
420 : The initial right and left spaces are rough approximations to the right and/or
421 : left singular subspaces from which the solver starts to iterate.
422 : It is not necessary to provide both sets of vectors.
423 :
424 : Some solvers start to iterate on a single vector (initial vector). In that case,
425 : the other vectors are ignored.
426 :
427 : These vectors do not persist from one SVDSolve() call to the other, so the
428 : initial space should be set every time.
429 :
430 : The vectors do not need to be mutually orthonormal, since they are explicitly
431 : orthonormalized internally.
432 :
433 : Common usage of this function is when the user can provide a rough approximation
434 : of the wanted singular space. Then, convergence may be faster.
435 :
436 : Level: intermediate
437 :
438 : .seealso: SVDSetUp()
439 : @*/
440 34 : PetscErrorCode SVDSetInitialSpaces(SVD svd,PetscInt nr,Vec isr[],PetscInt nl,Vec isl[])
441 : {
442 34 : PetscFunctionBegin;
443 34 : PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
444 102 : PetscValidLogicalCollectiveInt(svd,nr,2);
445 102 : PetscValidLogicalCollectiveInt(svd,nl,4);
446 34 : PetscCheck(nr>=0,PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_OUTOFRANGE,"Argument nr cannot be negative");
447 34 : if (nr>0) {
448 34 : PetscAssertPointer(isr,3);
449 34 : PetscValidHeaderSpecific(*isr,VEC_CLASSID,3);
450 : }
451 34 : PetscCheck(nl>=0,PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_OUTOFRANGE,"Argument nl cannot be negative");
452 34 : if (nl>0) {
453 34 : PetscAssertPointer(isl,5);
454 34 : PetscValidHeaderSpecific(*isl,VEC_CLASSID,5);
455 : }
456 34 : PetscCall(SlepcBasisReference_Private(nr,isr,&svd->nini,&svd->IS));
457 34 : PetscCall(SlepcBasisReference_Private(nl,isl,&svd->ninil,&svd->ISL));
458 34 : if (nr>0 || nl>0) svd->state = SVD_STATE_INITIAL;
459 34 : PetscFunctionReturn(PETSC_SUCCESS);
460 : }
461 :
462 : /*
463 : SVDSetDimensions_Default - Set reasonable values for ncv, mpd if not set
464 : by the user. This is called at setup.
465 : */
466 153 : PetscErrorCode SVDSetDimensions_Default(SVD svd)
467 : {
468 153 : PetscInt N,M,P,maxnsol;
469 :
470 153 : PetscFunctionBegin;
471 153 : PetscCall(MatGetSize(svd->OP,&M,&N));
472 153 : maxnsol = PetscMin(M,N);
473 153 : if (svd->isgeneralized) {
474 67 : PetscCall(MatGetSize(svd->OPb,&P,NULL));
475 67 : maxnsol = PetscMin(maxnsol,P);
476 : }
477 153 : if (svd->nsv==0 && svd->stop!=SVD_STOP_THRESHOLD) svd->nsv = 1;
478 153 : if (svd->ncv!=PETSC_DETERMINE) { /* ncv set */
479 65 : PetscCheck(svd->ncv>=svd->nsv,PetscObjectComm((PetscObject)svd),PETSC_ERR_USER_INPUT,"The value of ncv must be at least nsv");
480 88 : } else if (svd->mpd!=PETSC_DETERMINE) { /* mpd set */
481 2 : svd->ncv = PetscMin(maxnsol,svd->nsv+svd->mpd);
482 : } else { /* neither set: defaults depend on nsv being small or large */
483 86 : if (svd->nsv<500) svd->ncv = PetscMin(maxnsol,PetscMax(2*svd->nsv,10));
484 : else {
485 0 : svd->mpd = 500;
486 0 : svd->ncv = PetscMin(maxnsol,svd->nsv+svd->mpd);
487 : }
488 : }
489 153 : if (svd->mpd==PETSC_DETERMINE) svd->mpd = svd->ncv;
490 153 : PetscFunctionReturn(PETSC_SUCCESS);
491 : }
492 :
493 : /*@
494 : SVDAllocateSolution - Allocate memory storage for common variables such
495 : as the singular values and the basis vectors.
496 :
497 : Collective
498 :
499 : Input Parameters:
500 : + svd - singular value solver context
501 : - extra - number of additional positions, used for methods that require a
502 : working basis slightly larger than ncv
503 :
504 : Developer Notes:
505 : This is SLEPC_EXTERN because it may be required by user plugin SVD
506 : implementations.
507 :
508 : This is called at setup after setting the value of ncv and the flag leftbasis.
509 :
510 : Level: developer
511 :
512 : .seealso: SVDSetUp()
513 : @*/
514 286 : PetscErrorCode SVDAllocateSolution(SVD svd,PetscInt extra)
515 : {
516 286 : PetscInt oldsize,requested;
517 286 : Vec tr,tl;
518 :
519 286 : PetscFunctionBegin;
520 286 : requested = svd->ncv + extra;
521 :
522 : /* oldsize is zero if this is the first time setup is called */
523 286 : PetscCall(BVGetSizes(svd->V,NULL,NULL,&oldsize));
524 :
525 : /* allocate sigma */
526 286 : if (requested != oldsize || !svd->sigma) {
527 251 : PetscCall(PetscFree3(svd->sigma,svd->perm,svd->errest));
528 251 : if (svd->sign) PetscCall(PetscFree(svd->sign));
529 251 : PetscCall(PetscMalloc3(requested,&svd->sigma,requested,&svd->perm,requested,&svd->errest));
530 251 : if (svd->ishyperbolic) PetscCall(PetscMalloc1(requested,&svd->sign));
531 : }
532 : /* allocate V */
533 286 : if (!svd->V) PetscCall(SVDGetBV(svd,&svd->V,NULL));
534 286 : if (!oldsize) {
535 242 : if (!((PetscObject)svd->V)->type_name) PetscCall(BVSetType(svd->V,BVMAT));
536 242 : PetscCall(MatCreateVecsEmpty(svd->A,&tr,NULL));
537 242 : PetscCall(BVSetSizesFromVec(svd->V,tr,requested));
538 242 : PetscCall(VecDestroy(&tr));
539 44 : } else PetscCall(BVResize(svd->V,requested,PETSC_FALSE));
540 : /* allocate U */
541 286 : if (svd->leftbasis && !svd->isgeneralized) {
542 135 : if (!svd->U) PetscCall(SVDGetBV(svd,NULL,&svd->U));
543 135 : if (!oldsize) {
544 110 : if (!((PetscObject)svd->U)->type_name) PetscCall(BVSetType(svd->U,((PetscObject)svd->V)->type_name));
545 110 : PetscCall(MatCreateVecsEmpty(svd->A,NULL,&tl));
546 110 : PetscCall(BVSetSizesFromVec(svd->U,tl,requested));
547 110 : PetscCall(VecDestroy(&tl));
548 25 : } else PetscCall(BVResize(svd->U,requested,PETSC_FALSE));
549 151 : } else if (svd->isgeneralized) { /* left basis for the GSVD */
550 102 : if (!svd->U) PetscCall(SVDGetBV(svd,NULL,&svd->U));
551 102 : if (!oldsize) {
552 91 : if (!((PetscObject)svd->U)->type_name) PetscCall(BVSetType(svd->U,((PetscObject)svd->V)->type_name));
553 91 : PetscCall(SVDCreateLeftTemplate(svd,&tl));
554 91 : PetscCall(BVSetSizesFromVec(svd->U,tl,requested));
555 91 : PetscCall(VecDestroy(&tl));
556 11 : } else PetscCall(BVResize(svd->U,requested,PETSC_FALSE));
557 : }
558 286 : PetscFunctionReturn(PETSC_SUCCESS);
559 : }
560 :
561 : /*@
562 : SVDReallocateSolution - Reallocate memory storage for common variables such
563 : as the singular values and the basis vectors.
564 :
565 : Collective
566 :
567 : Input Parameters:
568 : + svd - singular value solver context
569 : - newsize - new size
570 :
571 : Developer Notes:
572 : This is SLEPC_EXTERN because it may be required by user plugin SVD
573 : implementations.
574 :
575 : This is called during the iteration in case the threshold stopping test has
576 : been selected.
577 :
578 : Level: developer
579 :
580 : .seealso: SVDAllocateSolution(), SVDSetThreshold()
581 : @*/
582 7 : PetscErrorCode SVDReallocateSolution(SVD svd,PetscInt newsize)
583 : {
584 7 : PetscInt oldsize,*nperm;
585 7 : PetscReal *nsigma,*nerrest,*nsign;
586 :
587 7 : PetscFunctionBegin;
588 7 : PetscCall(BVGetSizes(svd->V,NULL,NULL,&oldsize));
589 7 : if (oldsize>=newsize) PetscFunctionReturn(PETSC_SUCCESS);
590 7 : PetscCall(PetscInfo(svd,"Reallocating basis vectors to %" PetscInt_FMT " columns\n",newsize));
591 :
592 : /* reallocate sigma */
593 7 : PetscCall(PetscMalloc3(newsize,&nsigma,newsize,&nperm,newsize,&nerrest));
594 7 : PetscCall(PetscArraycpy(nsigma,svd->sigma,oldsize));
595 7 : PetscCall(PetscArraycpy(nperm,svd->perm,oldsize));
596 7 : PetscCall(PetscArraycpy(nerrest,svd->errest,oldsize));
597 7 : PetscCall(PetscFree3(svd->sigma,svd->perm,svd->errest));
598 7 : svd->sigma = nsigma;
599 7 : svd->perm = nperm;
600 7 : svd->errest = nerrest;
601 7 : if (svd->ishyperbolic) {
602 3 : PetscCall(PetscMalloc1(newsize,&nsign));
603 3 : PetscCall(PetscArraycpy(nsign,svd->sign,oldsize));
604 3 : PetscCall(PetscFree(svd->sign));
605 3 : svd->sign = nsign;
606 : }
607 : /* reallocate V,U */
608 7 : PetscCall(BVResize(svd->V,newsize,PETSC_TRUE));
609 7 : if (svd->leftbasis || svd->isgeneralized) PetscCall(BVResize(svd->U,newsize,PETSC_TRUE));
610 7 : PetscFunctionReturn(PETSC_SUCCESS);
611 : }
|
