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