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: `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: 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 (M == 0 || N == 0) PetscFunctionReturn(PETSC_SUCCESS);
287: if (svd->isgeneralized) {
288: PetscCall(MatGetSize(svd->OPb,&P,NULL));
289: PetscCheck(M+P>=N,PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"The case when [A;B] has less rows than columns is not supported");
290: } else if (svd->ishyperbolic) {
291: PetscCall(VecGetSize(svd->omega,&Nom));
292: PetscCall(VecGetLocalSize(svd->omega,&nom));
293: PetscCall(MatGetLocalSize(svd->OP,&m,NULL));
294: 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);
295: 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);
296: }
298: /* build transpose matrix */
299: PetscCall(MatDestroy(&svd->A));
300: PetscCall(MatDestroy(&svd->AT));
301: PetscCall(PetscObjectReference((PetscObject)svd->OP));
302: if (svd->expltrans) {
303: if (svd->isgeneralized || M>=N) {
304: svd->A = svd->OP;
305: PetscCall(MatHermitianTranspose(svd->OP,MAT_INITIAL_MATRIX,&svd->AT));
306: } else {
307: PetscCall(MatHermitianTranspose(svd->OP,MAT_INITIAL_MATRIX,&svd->A));
308: svd->AT = svd->OP;
309: }
310: } else {
311: if (svd->isgeneralized || M>=N) {
312: svd->A = svd->OP;
313: PetscCall(MatCreateHermitianTranspose(svd->OP,&svd->AT));
314: } else {
315: PetscCall(MatCreateHermitianTranspose(svd->OP,&svd->A));
316: svd->AT = svd->OP;
317: }
318: }
320: /* build transpose matrix B for GSVD */
321: if (svd->isgeneralized) {
322: PetscCall(MatDestroy(&svd->B));
323: PetscCall(MatDestroy(&svd->BT));
324: PetscCall(PetscObjectReference((PetscObject)svd->OPb));
325: if (svd->expltrans) {
326: svd->B = svd->OPb;
327: PetscCall(MatHermitianTranspose(svd->OPb,MAT_INITIAL_MATRIX,&svd->BT));
328: } else {
329: svd->B = svd->OPb;
330: PetscCall(MatCreateHermitianTranspose(svd->OPb,&svd->BT));
331: }
332: }
334: if (!svd->isgeneralized && M<N) {
335: /* swap initial vectors */
336: if (svd->nini || svd->ninil) {
337: T=svd->ISL; svd->ISL=svd->IS; svd->IS=T;
338: k=svd->ninil; svd->ninil=svd->nini; svd->nini=k;
339: }
340: /* swap basis vectors */
341: if (!svd->swapped) { /* only the first time in case of multiple calls */
342: bv=svd->V; svd->V=svd->U; svd->U=bv;
343: svd->swapped = PETSC_TRUE;
344: }
345: }
347: maxnsol = svd->isgeneralized? PetscMin(PetscMin(M,N),P): PetscMin(M,N);
348: svd->ncv = PetscMin(svd->ncv,maxnsol);
349: svd->nsv = PetscMin(svd->nsv,maxnsol);
350: PetscCheck(svd->ncv==PETSC_DETERMINE || svd->nsv<=svd->ncv,PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_OUTOFRANGE,"nsv bigger than ncv");
352: /* relative convergence criterion is not allowed in GSVD */
353: if (svd->conv==(SVDConv)-1) PetscCall(SVDSetConvergenceTest(svd,svd->isgeneralized?SVD_CONV_NORM:SVD_CONV_REL));
354: PetscCheck(!svd->isgeneralized || svd->conv!=SVD_CONV_REL,PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"Relative convergence criterion is not allowed in GSVD");
356: /* initialization of matrix norm (standard case only, for GSVD it is done inside setup()) */
357: if (!svd->isgeneralized && svd->conv==SVD_CONV_NORM && !svd->nrma) PetscCall(MatNorm(svd->OP,NORM_INFINITY,&svd->nrma));
359: /* threshold stopping test */
360: if (svd->stop==SVD_STOP_THRESHOLD) {
361: PetscCheck(svd->thres>0.0,PetscObjectComm((PetscObject)svd),PETSC_ERR_USER_INPUT,"Must give a threshold value with SVDSetThreshold()");
362: PetscCheck(svd->which==SVD_LARGEST || !svd->threlative,PetscObjectComm((PetscObject)svd),PETSC_ERR_USER_INPUT,"Can only use a relative threshold when which=SVD_LARGEST");
363: PetscCall(PetscNew(&ctx));
364: ctx->thres = svd->thres;
365: ctx->threlative = svd->threlative;
366: ctx->which = svd->which;
367: PetscCall(SVDSetStoppingTestFunction(svd,SVDStoppingThreshold,ctx,PetscCtxDestroyDefault));
368: }
370: /* call specific solver setup */
371: PetscUseTypeMethod(svd,setup);
373: /* set tolerance if not yet set */
374: if (svd->tol==(PetscReal)PETSC_DETERMINE) svd->tol = SLEPC_DEFAULT_TOL;
376: /* fill sorting criterion context */
377: PetscCall(DSGetSlepcSC(svd->ds,&sc));
378: sc->comparison = (svd->which==SVD_LARGEST)? SlepcCompareLargestReal: SlepcCompareSmallestReal;
379: sc->comparisonctx = NULL;
380: sc->map = NULL;
381: sc->mapobj = NULL;
383: /* process initial vectors */
384: if (svd->nini<0) {
385: k = -svd->nini;
386: PetscCheck(k<=svd->ncv,PetscObjectComm((PetscObject)svd),PETSC_ERR_USER_INPUT,"The number of initial vectors is larger than ncv");
387: PetscCall(BVInsertVecs(svd->V,0,&k,svd->IS,PETSC_TRUE));
388: PetscCall(SlepcBasisDestroy_Private(&svd->nini,&svd->IS));
389: svd->nini = k;
390: }
391: if (svd->ninil<0) {
392: k = 0;
393: if (svd->leftbasis) {
394: k = -svd->ninil;
395: PetscCheck(k<=svd->ncv,PetscObjectComm((PetscObject)svd),PETSC_ERR_USER_INPUT,"The number of left initial vectors is larger than ncv");
396: PetscCall(BVInsertVecs(svd->U,0,&k,svd->ISL,PETSC_TRUE));
397: } else PetscCall(PetscInfo(svd,"Ignoring initial left vectors\n"));
398: PetscCall(SlepcBasisDestroy_Private(&svd->ninil,&svd->ISL));
399: svd->ninil = k;
400: }
402: PetscCall(PetscLogEventEnd(SVD_SetUp,svd,0,0,0));
403: svd->state = SVD_STATE_SETUP;
404: PetscFunctionReturn(PETSC_SUCCESS);
405: }
407: /*@
408: SVDSetInitialSpaces - Specify two basis of vectors that constitute the initial
409: right and/or left spaces.
411: Collective
413: Input Parameters:
414: + svd - the singular value solver context
415: . nr - number of right vectors
416: . isr - set of basis vectors of the right initial space
417: . nl - number of left vectors
418: - isl - set of basis vectors of the left initial space
420: Notes:
421: The initial right and left spaces are rough approximations to the right and/or
422: left singular subspaces from which the solver starts to iterate.
423: It is not necessary to provide both sets of vectors.
425: Some solvers start to iterate on a single vector (initial vector). In that case,
426: the other vectors are ignored.
428: These vectors do not persist from one SVDSolve() call to the other, so the
429: initial space should be set every time.
431: The vectors do not need to be mutually orthonormal, since they are explicitly
432: orthonormalized internally.
434: Common usage of this function is when the user can provide a rough approximation
435: of the wanted singular space. Then, convergence may be faster.
437: Level: intermediate
439: .seealso: `SVDSetUp()`
440: @*/
441: PetscErrorCode SVDSetInitialSpaces(SVD svd,PetscInt nr,Vec isr[],PetscInt nl,Vec isl[])
442: {
443: PetscFunctionBegin;
447: PetscCheck(nr>=0,PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_OUTOFRANGE,"Argument nr cannot be negative");
448: if (nr>0) {
449: PetscAssertPointer(isr,3);
451: }
452: PetscCheck(nl>=0,PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_OUTOFRANGE,"Argument nl cannot be negative");
453: if (nl>0) {
454: PetscAssertPointer(isl,5);
456: }
457: PetscCall(SlepcBasisReference_Private(nr,isr,&svd->nini,&svd->IS));
458: PetscCall(SlepcBasisReference_Private(nl,isl,&svd->ninil,&svd->ISL));
459: if (nr>0 || nl>0) svd->state = SVD_STATE_INITIAL;
460: PetscFunctionReturn(PETSC_SUCCESS);
461: }
463: /*
464: SVDSetDimensions_Default - Set reasonable values for ncv, mpd if not set
465: by the user. This is called at setup.
466: */
467: PetscErrorCode SVDSetDimensions_Default(SVD svd)
468: {
469: PetscInt N,M,P,maxnsol;
471: PetscFunctionBegin;
472: PetscCall(MatGetSize(svd->OP,&M,&N));
473: maxnsol = PetscMin(M,N);
474: if (svd->isgeneralized) {
475: PetscCall(MatGetSize(svd->OPb,&P,NULL));
476: maxnsol = PetscMin(maxnsol,P);
477: }
478: if (svd->nsv==0 && svd->stop!=SVD_STOP_THRESHOLD) svd->nsv = 1;
479: if (svd->ncv!=PETSC_DETERMINE) { /* ncv set */
480: PetscCheck(svd->ncv>=svd->nsv,PetscObjectComm((PetscObject)svd),PETSC_ERR_USER_INPUT,"The value of ncv must be at least nsv");
481: } else if (svd->mpd!=PETSC_DETERMINE) { /* mpd set */
482: svd->ncv = PetscMin(maxnsol,svd->nsv+svd->mpd);
483: } else { /* neither set: defaults depend on nsv being small or large */
484: if (svd->nsv<500) svd->ncv = PetscMin(maxnsol,PetscMax(2*svd->nsv,10));
485: else {
486: svd->mpd = 500;
487: svd->ncv = PetscMin(maxnsol,svd->nsv+svd->mpd);
488: }
489: }
490: if (svd->mpd==PETSC_DETERMINE) svd->mpd = svd->ncv;
491: PetscFunctionReturn(PETSC_SUCCESS);
492: }
494: /*@
495: SVDAllocateSolution - Allocate memory storage for common variables such
496: as the singular values and the basis vectors.
498: Collective
500: Input Parameters:
501: + svd - singular value solver context
502: - extra - number of additional positions, used for methods that require a
503: working basis slightly larger than ncv
505: Developer Notes:
506: This is SLEPC_EXTERN because it may be required by user plugin SVD
507: implementations.
509: This is called at setup after setting the value of ncv and the flag leftbasis.
511: Level: developer
513: .seealso: `SVDSetUp()`
514: @*/
515: PetscErrorCode SVDAllocateSolution(SVD svd,PetscInt extra)
516: {
517: PetscInt oldsize,requested;
518: Vec tr,tl;
520: PetscFunctionBegin;
521: requested = svd->ncv + extra;
523: /* oldsize is zero if this is the first time setup is called */
524: PetscCall(BVGetSizes(svd->V,NULL,NULL,&oldsize));
526: /* allocate sigma */
527: if (requested != oldsize || !svd->sigma) {
528: PetscCall(PetscFree3(svd->sigma,svd->perm,svd->errest));
529: if (svd->sign) PetscCall(PetscFree(svd->sign));
530: PetscCall(PetscMalloc3(requested,&svd->sigma,requested,&svd->perm,requested,&svd->errest));
531: if (svd->ishyperbolic) PetscCall(PetscMalloc1(requested,&svd->sign));
532: }
533: /* allocate V */
534: if (!svd->V) PetscCall(SVDGetBV(svd,&svd->V,NULL));
535: if (!oldsize) {
536: if (!((PetscObject)svd->V)->type_name) PetscCall(BVSetType(svd->V,BVMAT));
537: PetscCall(MatCreateVecsEmpty(svd->A,&tr,NULL));
538: PetscCall(BVSetSizesFromVec(svd->V,tr,requested));
539: PetscCall(VecDestroy(&tr));
540: } else PetscCall(BVResize(svd->V,requested,PETSC_FALSE));
541: /* allocate U */
542: if (svd->leftbasis && !svd->isgeneralized) {
543: if (!svd->U) PetscCall(SVDGetBV(svd,NULL,&svd->U));
544: if (!oldsize) {
545: if (!((PetscObject)svd->U)->type_name) PetscCall(BVSetType(svd->U,((PetscObject)svd->V)->type_name));
546: PetscCall(MatCreateVecsEmpty(svd->A,NULL,&tl));
547: PetscCall(BVSetSizesFromVec(svd->U,tl,requested));
548: PetscCall(VecDestroy(&tl));
549: } else PetscCall(BVResize(svd->U,requested,PETSC_FALSE));
550: } else if (svd->isgeneralized) { /* left basis for the GSVD */
551: if (!svd->U) PetscCall(SVDGetBV(svd,NULL,&svd->U));
552: if (!oldsize) {
553: if (!((PetscObject)svd->U)->type_name) PetscCall(BVSetType(svd->U,((PetscObject)svd->V)->type_name));
554: PetscCall(SVDCreateLeftTemplate(svd,&tl));
555: PetscCall(BVSetSizesFromVec(svd->U,tl,requested));
556: PetscCall(VecDestroy(&tl));
557: } else PetscCall(BVResize(svd->U,requested,PETSC_FALSE));
558: }
559: PetscFunctionReturn(PETSC_SUCCESS);
560: }
562: /*@
563: SVDReallocateSolution - Reallocate memory storage for common variables such
564: as the singular values and the basis vectors.
566: Collective
568: Input Parameters:
569: + svd - singular value solver context
570: - newsize - new size
572: Developer Notes:
573: This is SLEPC_EXTERN because it may be required by user plugin SVD
574: implementations.
576: This is called during the iteration in case the threshold stopping test has
577: been selected.
579: Level: developer
581: .seealso: `SVDAllocateSolution()`, `SVDSetThreshold()`
582: @*/
583: PetscErrorCode SVDReallocateSolution(SVD svd,PetscInt newsize)
584: {
585: PetscInt oldsize,*nperm;
586: PetscReal *nsigma,*nerrest,*nsign;
588: PetscFunctionBegin;
589: PetscCall(BVGetSizes(svd->V,NULL,NULL,&oldsize));
590: if (oldsize>=newsize) PetscFunctionReturn(PETSC_SUCCESS);
591: PetscCall(PetscInfo(svd,"Reallocating basis vectors to %" PetscInt_FMT " columns\n",newsize));
593: /* reallocate sigma */
594: PetscCall(PetscMalloc3(newsize,&nsigma,newsize,&nperm,newsize,&nerrest));
595: PetscCall(PetscArraycpy(nsigma,svd->sigma,oldsize));
596: PetscCall(PetscArraycpy(nperm,svd->perm,oldsize));
597: PetscCall(PetscArraycpy(nerrest,svd->errest,oldsize));
598: PetscCall(PetscFree3(svd->sigma,svd->perm,svd->errest));
599: svd->sigma = nsigma;
600: svd->perm = nperm;
601: svd->errest = nerrest;
602: if (svd->ishyperbolic) {
603: PetscCall(PetscMalloc1(newsize,&nsign));
604: PetscCall(PetscArraycpy(nsign,svd->sign,oldsize));
605: PetscCall(PetscFree(svd->sign));
606: svd->sign = nsign;
607: }
608: /* reallocate V,U */
609: PetscCall(BVResize(svd->V,newsize,PETSC_TRUE));
610: if (svd->leftbasis || svd->isgeneralized) PetscCall(BVResize(svd->U,newsize,PETSC_TRUE));
611: PetscFunctionReturn(PETSC_SUCCESS);
612: }