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 229 : PetscErrorCode SVDSetOperators(SVD svd,Mat A,Mat B)
31 : {
32 229 : PetscInt Ma,Na,Mb,Nb,ma,na,mb,nb,M0,N0,m0,n0;
33 229 : PetscBool samesize=PETSC_TRUE;
34 :
35 229 : PetscFunctionBegin;
36 229 : PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
37 229 : PetscValidHeaderSpecific(A,MAT_CLASSID,2);
38 229 : if (B) PetscValidHeaderSpecific(B,MAT_CLASSID,3);
39 229 : PetscCheckSameComm(svd,1,A,2);
40 229 : if (B) PetscCheckSameComm(svd,1,B,3);
41 :
42 : /* Check matrix sizes */
43 229 : PetscCall(MatGetSize(A,&Ma,&Na));
44 229 : PetscCall(MatGetLocalSize(A,&ma,&na));
45 229 : 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 229 : if (B) {
51 85 : PetscCall(MatGetSize(B,&Mb,&Nb));
52 85 : PetscCall(MatGetLocalSize(B,&mb,&nb));
53 85 : 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 85 : 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 85 : 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 229 : PetscCall(PetscObjectReference((PetscObject)A));
63 229 : if (B) PetscCall(PetscObjectReference((PetscObject)B));
64 229 : if (svd->state && !samesize) PetscCall(SVDReset(svd));
65 : else {
66 229 : PetscCall(MatDestroy(&svd->OP));
67 229 : PetscCall(MatDestroy(&svd->OPb));
68 229 : PetscCall(MatDestroy(&svd->A));
69 229 : PetscCall(MatDestroy(&svd->B));
70 229 : PetscCall(MatDestroy(&svd->AT));
71 229 : PetscCall(MatDestroy(&svd->BT));
72 : }
73 229 : svd->nrma = 0.0;
74 229 : svd->nrmb = 0.0;
75 229 : svd->OP = A;
76 229 : svd->OPb = B;
77 229 : svd->state = SVD_STATE_INITIAL;
78 229 : 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 11 : PetscErrorCode SVDGetOperators(SVD svd,Mat *A,Mat *B)
98 : {
99 11 : PetscFunctionBegin;
100 11 : PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
101 11 : if (A) *A = svd->OP;
102 11 : if (B) *B = svd->OPb;
103 11 : 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 26 : PetscErrorCode SVDSetSignature(SVD svd,Vec omega)
124 : {
125 26 : PetscInt N,Ma,n,ma;
126 :
127 26 : PetscFunctionBegin;
128 26 : PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
129 26 : if (omega) {
130 26 : PetscValidHeaderSpecific(omega,VEC_CLASSID,2);
131 26 : PetscCheckSameComm(svd,1,omega,2);
132 : }
133 :
134 26 : if (omega && svd->OP) { /* Check sizes */
135 25 : PetscCall(VecGetSize(omega,&N));
136 25 : PetscCall(VecGetLocalSize(omega,&n));
137 25 : PetscCall(MatGetSize(svd->OP,&Ma,NULL));
138 25 : PetscCall(MatGetLocalSize(svd->OP,&ma,NULL));
139 25 : 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 25 : 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 26 : PetscCall(VecDestroy(&svd->omega));
144 26 : if (omega) {
145 26 : PetscCall(VecDuplicate(omega,&svd->omega));
146 26 : PetscCall(VecCopy(omega,svd->omega));
147 : }
148 26 : svd->state = SVD_STATE_INITIAL;
149 26 : 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 435 : PetscErrorCode SVDSetDSType(SVD svd)
201 : {
202 435 : PetscFunctionBegin;
203 435 : PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
204 435 : PetscTryTypeMethod(svd,setdstype);
205 435 : 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 245 : PetscErrorCode SVDSetUp(SVD svd)
227 : {
228 245 : PetscBool flg;
229 245 : PetscInt M,N,P=0,k,maxnsol,m,Nom,nom;
230 245 : SlepcSC sc;
231 245 : Vec *T;
232 245 : BV bv;
233 :
234 245 : PetscFunctionBegin;
235 245 : PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
236 245 : if (svd->state) PetscFunctionReturn(PETSC_SUCCESS);
237 245 : PetscCall(PetscLogEventBegin(SVD_SetUp,svd,0,0,0));
238 :
239 : /* reset the convergence flag from the previous solves */
240 245 : svd->reason = SVD_CONVERGED_ITERATING;
241 :
242 : /* set default solver type (SVDSetFromOptions was not called) */
243 245 : if (!((PetscObject)svd)->type_name) PetscCall(SVDSetType(svd,SVDCROSS));
244 245 : if (!svd->ds) PetscCall(SVDGetDS(svd,&svd->ds));
245 :
246 : /* check matrices */
247 245 : PetscCheck(svd->OP,PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_WRONGSTATE,"SVDSetOperators() must be called first");
248 :
249 245 : if (!svd->problem_type) { /* set default problem type */
250 116 : 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 98 : if (svd->omega) PetscCall(SVDSetProblemType(svd,SVD_HYPERBOLIC));
255 97 : else PetscCall(SVDSetProblemType(svd,SVD_STANDARD));
256 : }
257 : } else { /* check consistency of problem type set by user */
258 129 : if (svd->OPb) {
259 67 : PetscCheck(svd->isgeneralized,PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_INCOMP,"Inconsistent SVD state: the problem type does not match the number of matrices");
260 67 : PetscCheck(!svd->omega,PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"There is no support yet for generalized hyperbolic problems");
261 : } else {
262 62 : PetscCheck(!svd->isgeneralized,PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_INCOMP,"Inconsistent SVD state: the problem type does not match the number of matrices");
263 62 : 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 37 : 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 245 : PetscCall(SVDSetDSType(svd));
270 :
271 : /* determine how to handle the transpose */
272 245 : svd->expltrans = PETSC_TRUE;
273 245 : if (svd->impltrans) svd->expltrans = PETSC_FALSE;
274 : else {
275 234 : PetscCall(MatHasOperation(svd->OP,MATOP_TRANSPOSE,&flg));
276 234 : if (!flg) svd->expltrans = PETSC_FALSE;
277 : else {
278 233 : PetscCall(PetscObjectTypeCompareAny((PetscObject)svd,&flg,SVDLAPACK,SVDSCALAPACK,SVDKSVD,SVDELEMENTAL,""));
279 233 : if (flg) svd->expltrans = PETSC_FALSE;
280 : }
281 : }
282 :
283 : /* get matrix dimensions */
284 245 : PetscCall(MatGetSize(svd->OP,&M,&N));
285 245 : if (svd->isgeneralized) {
286 85 : PetscCall(MatGetSize(svd->OPb,&P,NULL));
287 85 : PetscCheck(M+P>=N,PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"The case when [A;B] has less rows than columns is not supported");
288 160 : } else if (svd->ishyperbolic) {
289 26 : PetscCall(VecGetSize(svd->omega,&Nom));
290 26 : PetscCall(VecGetLocalSize(svd->omega,&nom));
291 26 : PetscCall(MatGetLocalSize(svd->OP,&m,NULL));
292 26 : 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 26 : 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 245 : PetscCall(MatDestroy(&svd->A));
298 245 : PetscCall(MatDestroy(&svd->AT));
299 245 : PetscCall(PetscObjectReference((PetscObject)svd->OP));
300 245 : if (svd->expltrans) {
301 210 : if (svd->isgeneralized || M>=N) {
302 175 : svd->A = svd->OP;
303 175 : 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 35 : if (svd->isgeneralized || M>=N) {
310 31 : svd->A = svd->OP;
311 31 : PetscCall(MatCreateHermitianTranspose(svd->OP,&svd->AT));
312 : } else {
313 4 : PetscCall(MatCreateHermitianTranspose(svd->OP,&svd->A));
314 4 : svd->AT = svd->OP;
315 : }
316 : }
317 :
318 : /* build transpose matrix B for GSVD */
319 245 : if (svd->isgeneralized) {
320 85 : PetscCall(MatDestroy(&svd->B));
321 85 : PetscCall(MatDestroy(&svd->BT));
322 85 : PetscCall(PetscObjectReference((PetscObject)svd->OPb));
323 85 : if (svd->expltrans) {
324 78 : svd->B = svd->OPb;
325 78 : 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 245 : if (!svd->isgeneralized && M<N) {
333 : /* swap initial vectors */
334 39 : 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 39 : if (!svd->swapped) { /* only the first time in case of multiple calls */
340 32 : bv=svd->V; svd->V=svd->U; svd->U=bv;
341 32 : svd->swapped = PETSC_TRUE;
342 : }
343 : }
344 :
345 245 : maxnsol = svd->isgeneralized? PetscMin(PetscMin(M,N),P): PetscMin(M,N);
346 245 : svd->ncv = PetscMin(svd->ncv,maxnsol);
347 245 : svd->nsv = PetscMin(svd->nsv,maxnsol);
348 245 : PetscCheck(svd->ncv==PETSC_DETERMINE || 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 355 : if (svd->conv==(SVDConv)-1) PetscCall(SVDSetConvergenceTest(svd,svd->isgeneralized?SVD_CONV_NORM:SVD_CONV_REL));
352 245 : 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 245 : 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 245 : PetscUseTypeMethod(svd,setup);
359 :
360 : /* set tolerance if not yet set */
361 245 : if (svd->tol==(PetscReal)PETSC_DETERMINE) svd->tol = SLEPC_DEFAULT_TOL;
362 :
363 : /* fill sorting criterion context */
364 245 : PetscCall(DSGetSlepcSC(svd->ds,&sc));
365 245 : sc->comparison = (svd->which==SVD_LARGEST)? SlepcCompareLargestReal: SlepcCompareSmallestReal;
366 245 : sc->comparisonctx = NULL;
367 245 : sc->map = NULL;
368 245 : sc->mapobj = NULL;
369 :
370 : /* process initial vectors */
371 245 : 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 245 : if (svd->ninil<0) {
379 29 : k = 0;
380 29 : 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 6 : } else PetscCall(PetscInfo(svd,"Ignoring initial left vectors\n"));
385 29 : PetscCall(SlepcBasisDestroy_Private(&svd->ninil,&svd->ISL));
386 29 : svd->ninil = k;
387 : }
388 :
389 245 : PetscCall(PetscLogEventEnd(SVD_SetUp,svd,0,0,0));
390 245 : svd->state = SVD_STATE_SETUP;
391 245 : PetscFunctionReturn(PETSC_SUCCESS);
392 : }
393 :
394 : /*@
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 33 : PetscErrorCode SVDSetInitialSpaces(SVD svd,PetscInt nr,Vec isr[],PetscInt nl,Vec isl[])
429 : {
430 33 : PetscFunctionBegin;
431 33 : PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
432 99 : PetscValidLogicalCollectiveInt(svd,nr,2);
433 99 : PetscValidLogicalCollectiveInt(svd,nl,4);
434 33 : PetscCheck(nr>=0,PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_OUTOFRANGE,"Argument nr cannot be negative");
435 33 : if (nr>0) {
436 33 : PetscAssertPointer(isr,3);
437 33 : PetscValidHeaderSpecific(*isr,VEC_CLASSID,3);
438 : }
439 33 : PetscCheck(nl>=0,PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_OUTOFRANGE,"Argument nl cannot be negative");
440 33 : if (nl>0) {
441 33 : PetscAssertPointer(isl,5);
442 33 : PetscValidHeaderSpecific(*isl,VEC_CLASSID,5);
443 : }
444 33 : PetscCall(SlepcBasisReference_Private(nr,isr,&svd->nini,&svd->IS));
445 33 : PetscCall(SlepcBasisReference_Private(nl,isl,&svd->ninil,&svd->ISL));
446 33 : if (nr>0 || nl>0) svd->state = SVD_STATE_INITIAL;
447 33 : 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 131 : PetscErrorCode SVDSetDimensions_Default(SVD svd)
455 : {
456 131 : PetscInt N,M,P,maxnsol;
457 :
458 131 : PetscFunctionBegin;
459 131 : PetscCall(MatGetSize(svd->OP,&M,&N));
460 131 : maxnsol = PetscMin(M,N);
461 131 : if (svd->isgeneralized) {
462 56 : PetscCall(MatGetSize(svd->OPb,&P,NULL));
463 56 : maxnsol = PetscMin(maxnsol,P);
464 : }
465 131 : if (svd->ncv!=PETSC_DETERMINE) { /* ncv set */
466 52 : PetscCheck(svd->ncv>=svd->nsv,PetscObjectComm((PetscObject)svd),PETSC_ERR_USER_INPUT,"The value of ncv must be at least nsv");
467 79 : } else if (svd->mpd!=PETSC_DETERMINE) { /* 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 77 : 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 131 : if (svd->mpd==PETSC_DETERMINE) svd->mpd = svd->ncv;
477 131 : 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 245 : PetscErrorCode SVDAllocateSolution(SVD svd,PetscInt extra)
502 : {
503 245 : PetscInt oldsize,requested;
504 245 : Vec tr,tl;
505 :
506 245 : PetscFunctionBegin;
507 245 : requested = svd->ncv + extra;
508 :
509 : /* oldsize is zero if this is the first time setup is called */
510 245 : PetscCall(BVGetSizes(svd->V,NULL,NULL,&oldsize));
511 :
512 : /* allocate sigma */
513 245 : if (requested != oldsize || !svd->sigma) {
514 210 : PetscCall(PetscFree3(svd->sigma,svd->perm,svd->errest));
515 210 : if (svd->sign) PetscCall(PetscFree(svd->sign));
516 210 : PetscCall(PetscMalloc3(requested,&svd->sigma,requested,&svd->perm,requested,&svd->errest));
517 210 : if (svd->ishyperbolic) PetscCall(PetscMalloc1(requested,&svd->sign));
518 : }
519 : /* allocate V */
520 245 : if (!svd->V) PetscCall(SVDGetBV(svd,&svd->V,NULL));
521 245 : if (!oldsize) {
522 201 : if (!((PetscObject)svd->V)->type_name) PetscCall(BVSetType(svd->V,BVMAT));
523 201 : PetscCall(MatCreateVecsEmpty(svd->A,&tr,NULL));
524 201 : PetscCall(BVSetSizesFromVec(svd->V,tr,requested));
525 201 : PetscCall(VecDestroy(&tr));
526 44 : } else PetscCall(BVResize(svd->V,requested,PETSC_FALSE));
527 : /* allocate U */
528 245 : if (svd->leftbasis && !svd->isgeneralized) {
529 122 : if (!svd->U) PetscCall(SVDGetBV(svd,NULL,&svd->U));
530 122 : if (!oldsize) {
531 97 : if (!((PetscObject)svd->U)->type_name) PetscCall(BVSetType(svd->U,((PetscObject)svd->V)->type_name));
532 97 : PetscCall(MatCreateVecsEmpty(svd->A,NULL,&tl));
533 97 : PetscCall(BVSetSizesFromVec(svd->U,tl,requested));
534 97 : PetscCall(VecDestroy(&tl));
535 25 : } else PetscCall(BVResize(svd->U,requested,PETSC_FALSE));
536 123 : } else if (svd->isgeneralized) { /* left basis for the GSVD */
537 85 : if (!svd->U) PetscCall(SVDGetBV(svd,NULL,&svd->U));
538 85 : if (!oldsize) {
539 74 : if (!((PetscObject)svd->U)->type_name) PetscCall(BVSetType(svd->U,((PetscObject)svd->V)->type_name));
540 74 : PetscCall(SVDCreateLeftTemplate(svd,&tl));
541 74 : PetscCall(BVSetSizesFromVec(svd->U,tl,requested));
542 74 : PetscCall(VecDestroy(&tl));
543 11 : } else PetscCall(BVResize(svd->U,requested,PETSC_FALSE));
544 : }
545 245 : PetscFunctionReturn(PETSC_SUCCESS);
546 : }
|
