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 258 : PetscErrorCode SVDSetOperators(SVD svd,Mat A,Mat B)
31 : {
32 258 : PetscInt Ma,Na,Mb,Nb,ma,na,mb,nb,M0,N0,m0,n0;
33 258 : PetscBool samesize=PETSC_TRUE;
34 :
35 258 : PetscFunctionBegin;
36 258 : PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
37 258 : PetscValidHeaderSpecific(A,MAT_CLASSID,2);
38 258 : if (B) PetscValidHeaderSpecific(B,MAT_CLASSID,3);
39 258 : PetscCheckSameComm(svd,1,A,2);
40 258 : if (B) PetscCheckSameComm(svd,1,B,3);
41 :
42 : /* Check matrix sizes */
43 258 : PetscCall(MatGetSize(A,&Ma,&Na));
44 258 : PetscCall(MatGetLocalSize(A,&ma,&na));
45 258 : 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 258 : if (B) {
51 98 : PetscCall(MatGetSize(B,&Mb,&Nb));
52 98 : PetscCall(MatGetLocalSize(B,&mb,&nb));
53 98 : 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 98 : 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 98 : 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 258 : PetscCall(PetscObjectReference((PetscObject)A));
63 258 : if (B) PetscCall(PetscObjectReference((PetscObject)B));
64 258 : if (svd->state && !samesize) PetscCall(SVDReset(svd));
65 : else {
66 258 : PetscCall(MatDestroy(&svd->OP));
67 258 : PetscCall(MatDestroy(&svd->OPb));
68 258 : PetscCall(MatDestroy(&svd->A));
69 258 : PetscCall(MatDestroy(&svd->B));
70 258 : PetscCall(MatDestroy(&svd->AT));
71 258 : PetscCall(MatDestroy(&svd->BT));
72 : }
73 258 : svd->nrma = 0.0;
74 258 : svd->nrmb = 0.0;
75 258 : svd->OP = A;
76 258 : svd->OPb = B;
77 258 : svd->state = SVD_STATE_INITIAL;
78 258 : 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 36 : PetscErrorCode SVDSetSignature(SVD svd,Vec omega)
124 : {
125 36 : PetscInt N,Ma,n,ma;
126 :
127 36 : PetscFunctionBegin;
128 36 : PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
129 36 : if (omega) {
130 36 : PetscValidHeaderSpecific(omega,VEC_CLASSID,2);
131 36 : PetscCheckSameComm(svd,1,omega,2);
132 : }
133 :
134 36 : if (omega && svd->OP) { /* Check sizes */
135 35 : PetscCall(VecGetSize(omega,&N));
136 35 : PetscCall(VecGetLocalSize(omega,&n));
137 35 : PetscCall(MatGetSize(svd->OP,&Ma,NULL));
138 35 : PetscCall(MatGetLocalSize(svd->OP,&ma,NULL));
139 35 : 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 35 : 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 36 : PetscCall(VecDestroy(&svd->omega));
144 36 : if (omega) {
145 36 : PetscCall(VecDuplicate(omega,&svd->omega));
146 36 : PetscCall(VecCopy(omega,svd->omega));
147 : }
148 36 : svd->state = SVD_STATE_INITIAL;
149 36 : 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 493 : PetscErrorCode SVDSetDSType(SVD svd)
201 : {
202 493 : PetscFunctionBegin;
203 493 : PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
204 493 : PetscTryTypeMethod(svd,setdstype);
205 493 : 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 274 : PetscErrorCode SVDSetUp(SVD svd)
227 : {
228 274 : PetscBool flg;
229 274 : PetscInt M,N,P=0,k,maxnsol,m,Nom,nom;
230 274 : SlepcSC sc;
231 274 : Vec *T;
232 274 : BV bv;
233 :
234 274 : PetscFunctionBegin;
235 274 : PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
236 274 : if (svd->state) PetscFunctionReturn(PETSC_SUCCESS);
237 274 : PetscCall(PetscLogEventBegin(SVD_SetUp,svd,0,0,0));
238 :
239 : /* reset the convergence flag from the previous solves */
240 274 : svd->reason = SVD_CONVERGED_ITERATING;
241 :
242 : /* set default solver type (SVDSetFromOptions was not called) */
243 274 : if (!((PetscObject)svd)->type_name) PetscCall(SVDSetType(svd,SVDCROSS));
244 274 : if (!svd->ds) PetscCall(SVDGetDS(svd,&svd->ds));
245 :
246 : /* check matrices */
247 274 : PetscCheck(svd->OP,PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_WRONGSTATE,"SVDSetOperators() must be called first");
248 :
249 274 : if (!svd->problem_type) { /* set default problem type */
250 118 : if (svd->OPb) {
251 18 : PetscCheck(!svd->omega,PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"There is no support yet for generalized hyperbolic problems");
252 18 : PetscCall(SVDSetProblemType(svd,SVD_GENERALIZED));
253 : } else {
254 100 : if (svd->omega) PetscCall(SVDSetProblemType(svd,SVD_HYPERBOLIC));
255 99 : else PetscCall(SVDSetProblemType(svd,SVD_STANDARD));
256 : }
257 : } else { /* check consistency of problem type set by user */
258 156 : if (svd->OPb) {
259 80 : PetscCheck(svd->isgeneralized,PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_INCOMP,"Inconsistent SVD state: the problem type does not match the number of matrices");
260 80 : PetscCheck(!svd->omega,PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"There is no support yet for generalized hyperbolic problems");
261 : } else {
262 76 : PetscCheck(!svd->isgeneralized,PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_INCOMP,"Inconsistent SVD state: the problem type does not match the number of matrices");
263 76 : 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()");
264 41 : else PetscCheck(!svd->ishyperbolic,PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_INCOMP,"Inconsistent SVD state: a hyperbolic problem requires passing a signature with SVDSetSignature()");
265 : }
266 : }
267 :
268 : /* set DS type once the default problem type has been determined */
269 274 : PetscCall(SVDSetDSType(svd));
270 :
271 : /* determine how to handle the transpose */
272 274 : svd->expltrans = PETSC_TRUE;
273 274 : if (svd->impltrans) svd->expltrans = PETSC_FALSE;
274 : else {
275 255 : PetscCall(MatHasOperation(svd->OP,MATOP_TRANSPOSE,&flg));
276 255 : if (!flg) svd->expltrans = PETSC_FALSE;
277 : else {
278 254 : PetscCall(PetscObjectTypeCompareAny((PetscObject)svd,&flg,SVDLAPACK,SVDSCALAPACK,SVDKSVD,SVDELEMENTAL,""));
279 254 : if (flg) svd->expltrans = PETSC_FALSE;
280 : }
281 : }
282 :
283 : /* get matrix dimensions */
284 274 : PetscCall(MatGetSize(svd->OP,&M,&N));
285 274 : if (svd->isgeneralized) {
286 98 : PetscCall(MatGetSize(svd->OPb,&P,NULL));
287 98 : PetscCheck(M+P>=N,PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"The case when [A;B] has less rows than columns is not supported");
288 176 : } else if (svd->ishyperbolic) {
289 36 : PetscCall(VecGetSize(svd->omega,&Nom));
290 36 : PetscCall(VecGetLocalSize(svd->omega,&nom));
291 36 : PetscCall(MatGetLocalSize(svd->OP,&m,NULL));
292 36 : 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);
293 36 : 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);
294 : }
295 :
296 : /* build transpose matrix */
297 274 : PetscCall(MatDestroy(&svd->A));
298 274 : PetscCall(MatDestroy(&svd->AT));
299 274 : PetscCall(PetscObjectReference((PetscObject)svd->OP));
300 274 : if (svd->expltrans) {
301 230 : if (svd->isgeneralized || M>=N) {
302 195 : svd->A = svd->OP;
303 195 : PetscCall(MatHermitianTranspose(svd->OP,MAT_INITIAL_MATRIX,&svd->AT));
304 : } else {
305 35 : PetscCall(MatHermitianTranspose(svd->OP,MAT_INITIAL_MATRIX,&svd->A));
306 35 : svd->AT = svd->OP;
307 : }
308 : } else {
309 44 : if (svd->isgeneralized || M>=N) {
310 35 : svd->A = svd->OP;
311 35 : PetscCall(MatCreateHermitianTranspose(svd->OP,&svd->AT));
312 : } else {
313 9 : PetscCall(MatCreateHermitianTranspose(svd->OP,&svd->A));
314 9 : svd->AT = svd->OP;
315 : }
316 : }
317 :
318 : /* build transpose matrix B for GSVD */
319 274 : if (svd->isgeneralized) {
320 98 : PetscCall(MatDestroy(&svd->B));
321 98 : PetscCall(MatDestroy(&svd->BT));
322 98 : PetscCall(PetscObjectReference((PetscObject)svd->OPb));
323 98 : if (svd->expltrans) {
324 91 : svd->B = svd->OPb;
325 91 : PetscCall(MatHermitianTranspose(svd->OPb,MAT_INITIAL_MATRIX,&svd->BT));
326 : } else {
327 7 : svd->B = svd->OPb;
328 7 : PetscCall(MatCreateHermitianTranspose(svd->OPb,&svd->BT));
329 : }
330 : }
331 :
332 274 : if (!svd->isgeneralized && M<N) {
333 : /* swap initial vectors */
334 44 : if (svd->nini || svd->ninil) {
335 0 : T=svd->ISL; svd->ISL=svd->IS; svd->IS=T;
336 0 : k=svd->ninil; svd->ninil=svd->nini; svd->nini=k;
337 : }
338 : /* swap basis vectors */
339 44 : if (!svd->swapped) { /* only the first time in case of multiple calls */
340 37 : bv=svd->V; svd->V=svd->U; svd->U=bv;
341 37 : svd->swapped = PETSC_TRUE;
342 : }
343 : }
344 :
345 274 : maxnsol = svd->isgeneralized? PetscMin(PetscMin(M,N),P): PetscMin(M,N);
346 274 : svd->ncv = PetscMin(svd->ncv,maxnsol);
347 274 : svd->nsv = PetscMin(svd->nsv,maxnsol);
348 274 : PetscCheck(svd->ncv==PETSC_DEFAULT || svd->nsv<=svd->ncv,PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_OUTOFRANGE,"nsv bigger than ncv");
349 :
350 : /* relative convergence criterion is not allowed in GSVD */
351 399 : if (svd->conv==(SVDConv)-1) PetscCall(SVDSetConvergenceTest(svd,svd->isgeneralized?SVD_CONV_NORM:SVD_CONV_REL));
352 274 : PetscCheck(!svd->isgeneralized || svd->conv!=SVD_CONV_REL,PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"Relative convergence criterion is not allowed in GSVD");
353 :
354 : /* initialization of matrix norm (stardard case only, for GSVD it is done inside setup()) */
355 274 : if (!svd->isgeneralized && svd->conv==SVD_CONV_NORM && !svd->nrma) PetscCall(MatNorm(svd->OP,NORM_INFINITY,&svd->nrma));
356 :
357 : /* call specific solver setup */
358 274 : PetscUseTypeMethod(svd,setup);
359 :
360 : /* set tolerance if not yet set */
361 274 : if (svd->tol==(PetscReal)PETSC_DEFAULT) svd->tol = SLEPC_DEFAULT_TOL;
362 :
363 : /* fill sorting criterion context */
364 274 : PetscCall(DSGetSlepcSC(svd->ds,&sc));
365 274 : sc->comparison = (svd->which==SVD_LARGEST)? SlepcCompareLargestReal: SlepcCompareSmallestReal;
366 274 : sc->comparisonctx = NULL;
367 274 : sc->map = NULL;
368 274 : sc->mapobj = NULL;
369 :
370 : /* process initial vectors */
371 274 : if (svd->nini<0) {
372 25 : k = -svd->nini;
373 25 : PetscCheck(k<=svd->ncv,PetscObjectComm((PetscObject)svd),PETSC_ERR_USER_INPUT,"The number of initial vectors is larger than ncv");
374 25 : PetscCall(BVInsertVecs(svd->V,0,&k,svd->IS,PETSC_TRUE));
375 25 : PetscCall(SlepcBasisDestroy_Private(&svd->nini,&svd->IS));
376 25 : svd->nini = k;
377 : }
378 274 : if (svd->ninil<0) {
379 30 : k = 0;
380 30 : if (svd->leftbasis) {
381 23 : k = -svd->ninil;
382 23 : PetscCheck(k<=svd->ncv,PetscObjectComm((PetscObject)svd),PETSC_ERR_USER_INPUT,"The number of left initial vectors is larger than ncv");
383 23 : PetscCall(BVInsertVecs(svd->U,0,&k,svd->ISL,PETSC_TRUE));
384 7 : } else PetscCall(PetscInfo(svd,"Ignoring initial left vectors\n"));
385 30 : PetscCall(SlepcBasisDestroy_Private(&svd->ninil,&svd->ISL));
386 30 : svd->ninil = k;
387 : }
388 :
389 274 : PetscCall(PetscLogEventEnd(SVD_SetUp,svd,0,0,0));
390 274 : svd->state = SVD_STATE_SETUP;
391 274 : PetscFunctionReturn(PETSC_SUCCESS);
392 : }
393 :
394 : /*@C
395 : SVDSetInitialSpaces - Specify two basis of vectors that constitute the initial
396 : right and/or left spaces.
397 :
398 : Collective
399 :
400 : Input Parameters:
401 : + svd - the singular value solver context
402 : . nr - number of right vectors
403 : . isr - set of basis vectors of the right initial space
404 : . nl - number of left vectors
405 : - isl - set of basis vectors of the left initial space
406 :
407 : Notes:
408 : The initial right and left spaces are rough approximations to the right and/or
409 : left singular subspaces from which the solver starts to iterate.
410 : It is not necessary to provide both sets of vectors.
411 :
412 : Some solvers start to iterate on a single vector (initial vector). In that case,
413 : the other vectors are ignored.
414 :
415 : These vectors do not persist from one SVDSolve() call to the other, so the
416 : initial space should be set every time.
417 :
418 : The vectors do not need to be mutually orthonormal, since they are explicitly
419 : orthonormalized internally.
420 :
421 : Common usage of this function is when the user can provide a rough approximation
422 : of the wanted singular space. Then, convergence may be faster.
423 :
424 : Level: intermediate
425 :
426 : .seealso: SVDSetUp()
427 : @*/
428 34 : PetscErrorCode SVDSetInitialSpaces(SVD svd,PetscInt nr,Vec isr[],PetscInt nl,Vec isl[])
429 : {
430 34 : PetscFunctionBegin;
431 34 : PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
432 136 : PetscValidLogicalCollectiveInt(svd,nr,2);
433 136 : PetscValidLogicalCollectiveInt(svd,nl,4);
434 34 : PetscCheck(nr>=0,PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_OUTOFRANGE,"Argument nr cannot be negative");
435 34 : if (nr>0) {
436 34 : PetscAssertPointer(isr,3);
437 34 : PetscValidHeaderSpecific(*isr,VEC_CLASSID,3);
438 : }
439 34 : PetscCheck(nl>=0,PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_OUTOFRANGE,"Argument nl cannot be negative");
440 34 : if (nl>0) {
441 34 : PetscAssertPointer(isl,5);
442 34 : PetscValidHeaderSpecific(*isl,VEC_CLASSID,5);
443 : }
444 34 : PetscCall(SlepcBasisReference_Private(nr,isr,&svd->nini,&svd->IS));
445 34 : PetscCall(SlepcBasisReference_Private(nl,isl,&svd->ninil,&svd->ISL));
446 34 : if (nr>0 || nl>0) svd->state = SVD_STATE_INITIAL;
447 34 : PetscFunctionReturn(PETSC_SUCCESS);
448 : }
449 :
450 : /*
451 : SVDSetDimensions_Default - Set reasonable values for ncv, mpd if not set
452 : by the user. This is called at setup.
453 : */
454 145 : PetscErrorCode SVDSetDimensions_Default(SVD svd)
455 : {
456 145 : PetscInt N,M,P,maxnsol;
457 :
458 145 : PetscFunctionBegin;
459 145 : PetscCall(MatGetSize(svd->OP,&M,&N));
460 145 : maxnsol = PetscMin(M,N);
461 145 : if (svd->isgeneralized) {
462 64 : PetscCall(MatGetSize(svd->OPb,&P,NULL));
463 64 : maxnsol = PetscMin(maxnsol,P);
464 : }
465 145 : if (svd->ncv!=PETSC_DEFAULT) { /* ncv set */
466 61 : PetscCheck(svd->ncv>=svd->nsv,PetscObjectComm((PetscObject)svd),PETSC_ERR_USER_INPUT,"The value of ncv must be at least nsv");
467 84 : } else if (svd->mpd!=PETSC_DEFAULT) { /* mpd set */
468 2 : svd->ncv = PetscMin(maxnsol,svd->nsv+svd->mpd);
469 : } else { /* neither set: defaults depend on nsv being small or large */
470 82 : if (svd->nsv<500) svd->ncv = PetscMin(maxnsol,PetscMax(2*svd->nsv,10));
471 : else {
472 0 : svd->mpd = 500;
473 0 : svd->ncv = PetscMin(maxnsol,svd->nsv+svd->mpd);
474 : }
475 : }
476 145 : if (svd->mpd==PETSC_DEFAULT) svd->mpd = svd->ncv;
477 145 : PetscFunctionReturn(PETSC_SUCCESS);
478 : }
479 :
480 : /*@
481 : SVDAllocateSolution - Allocate memory storage for common variables such
482 : as the singular values and the basis vectors.
483 :
484 : Collective
485 :
486 : Input Parameters:
487 : + svd - eigensolver context
488 : - extra - number of additional positions, used for methods that require a
489 : working basis slightly larger than ncv
490 :
491 : Developer Notes:
492 : This is SLEPC_EXTERN because it may be required by user plugin SVD
493 : implementations.
494 :
495 : This is called at setup after setting the value of ncv and the flag leftbasis.
496 :
497 : Level: developer
498 :
499 : .seealso: SVDSetUp()
500 : @*/
501 274 : PetscErrorCode SVDAllocateSolution(SVD svd,PetscInt extra)
502 : {
503 274 : PetscInt oldsize,requested;
504 274 : Vec tr,tl;
505 :
506 274 : PetscFunctionBegin;
507 274 : requested = svd->ncv + extra;
508 :
509 : /* oldsize is zero if this is the first time setup is called */
510 274 : PetscCall(BVGetSizes(svd->V,NULL,NULL,&oldsize));
511 :
512 : /* allocate sigma */
513 274 : if (requested != oldsize || !svd->sigma) {
514 239 : PetscCall(PetscFree3(svd->sigma,svd->perm,svd->errest));
515 239 : if (svd->sign) PetscCall(PetscFree(svd->sign));
516 239 : PetscCall(PetscMalloc3(requested,&svd->sigma,requested,&svd->perm,requested,&svd->errest));
517 239 : if (svd->ishyperbolic) PetscCall(PetscMalloc1(requested,&svd->sign));
518 : }
519 : /* allocate V */
520 274 : if (!svd->V) PetscCall(SVDGetBV(svd,&svd->V,NULL));
521 274 : if (!oldsize) {
522 230 : if (!((PetscObject)svd->V)->type_name) PetscCall(BVSetType(svd->V,BVMAT));
523 230 : PetscCall(MatCreateVecsEmpty(svd->A,&tr,NULL));
524 230 : PetscCall(BVSetSizesFromVec(svd->V,tr,requested));
525 230 : PetscCall(VecDestroy(&tr));
526 44 : } else PetscCall(BVResize(svd->V,requested,PETSC_FALSE));
527 : /* allocate U */
528 274 : if (svd->leftbasis && !svd->isgeneralized) {
529 130 : if (!svd->U) PetscCall(SVDGetBV(svd,NULL,&svd->U));
530 130 : if (!oldsize) {
531 105 : if (!((PetscObject)svd->U)->type_name) PetscCall(BVSetType(svd->U,((PetscObject)svd->V)->type_name));
532 105 : PetscCall(MatCreateVecsEmpty(svd->A,NULL,&tl));
533 105 : PetscCall(BVSetSizesFromVec(svd->U,tl,requested));
534 105 : PetscCall(VecDestroy(&tl));
535 25 : } else PetscCall(BVResize(svd->U,requested,PETSC_FALSE));
536 144 : } else if (svd->isgeneralized) { /* left basis for the GSVD */
537 98 : if (!svd->U) PetscCall(SVDGetBV(svd,NULL,&svd->U));
538 98 : if (!oldsize) {
539 87 : if (!((PetscObject)svd->U)->type_name) PetscCall(BVSetType(svd->U,((PetscObject)svd->V)->type_name));
540 87 : PetscCall(SVDCreateLeftTemplate(svd,&tl));
541 87 : PetscCall(BVSetSizesFromVec(svd->U,tl,requested));
542 87 : PetscCall(VecDestroy(&tl));
543 11 : } else PetscCall(BVResize(svd->U,requested,PETSC_FALSE));
544 : }
545 274 : PetscFunctionReturn(PETSC_SUCCESS);
546 : }
|
