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