Actual source code: svdsetup.c
slepc-3.20.2 2024-03-15
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_WRONG,"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_WRONG,"Local size of signature (%" PetscInt_FMT ") does not match the local row size of A (%" PetscInt_FMT ")",n,ma);
141: }
143: if (omega) PetscCall(PetscObjectReference((PetscObject)omega));
144: PetscCall(VecDestroy(&svd->omega));
145: svd->omega = omega;
146: svd->state = SVD_STATE_INITIAL;
147: PetscFunctionReturn(PETSC_SUCCESS);
148: }
150: /*@
151: SVDGetSignature - Get the signature matrix defining a hyperbolic singular value problem.
153: Not Collective
155: Input Parameter:
156: . svd - the singular value solver context
158: Output Parameter:
159: . omega - a vector containing the diagonal elements of the signature matrix
161: Level: intermediate
163: .seealso: SVDSetSignature()
164: @*/
165: PetscErrorCode SVDGetSignature(SVD svd,Vec *omega)
166: {
167: PetscFunctionBegin;
169: PetscAssertPointer(omega,2);
170: *omega = svd->omega;
171: PetscFunctionReturn(PETSC_SUCCESS);
172: }
174: /*@
175: SVDSetDSType - Sets the type of the internal DS object based on the current
176: settings of the singular value solver.
178: Collective
180: Input Parameter:
181: . svd - singular value solver context
183: Note:
184: This function need not be called explicitly, since it will be called at
185: both SVDSetFromOptions() and SVDSetUp().
187: Level: developer
189: .seealso: SVDSetFromOptions(), SVDSetUp()
190: @*/
191: PetscErrorCode SVDSetDSType(SVD svd)
192: {
193: PetscFunctionBegin;
195: PetscTryTypeMethod(svd,setdstype);
196: PetscFunctionReturn(PETSC_SUCCESS);
197: }
199: /*@
200: SVDSetUp - Sets up all the internal data structures necessary for the
201: execution of the singular value solver.
203: Collective
205: Input Parameter:
206: . svd - singular value solver context
208: Notes:
209: This function need not be called explicitly in most cases, since SVDSolve()
210: calls it. It can be useful when one wants to measure the set-up time
211: separately from the solve time.
213: Level: developer
215: .seealso: SVDCreate(), SVDSolve(), SVDDestroy()
216: @*/
217: PetscErrorCode SVDSetUp(SVD svd)
218: {
219: PetscBool flg;
220: PetscInt M,N,P=0,k,maxnsol;
221: SlepcSC sc;
222: Vec *T;
223: BV bv;
225: PetscFunctionBegin;
227: if (svd->state) PetscFunctionReturn(PETSC_SUCCESS);
228: PetscCall(PetscLogEventBegin(SVD_SetUp,svd,0,0,0));
230: /* reset the convergence flag from the previous solves */
231: svd->reason = SVD_CONVERGED_ITERATING;
233: /* set default solver type (SVDSetFromOptions was not called) */
234: if (!((PetscObject)svd)->type_name) PetscCall(SVDSetType(svd,SVDCROSS));
235: if (!svd->ds) PetscCall(SVDGetDS(svd,&svd->ds));
236: PetscCall(SVDSetDSType(svd));
238: /* check matrices */
239: PetscCheck(svd->OP,PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_WRONGSTATE,"SVDSetOperators() must be called first");
241: if (!svd->problem_type) { /* set default problem type */
242: if (svd->OPb) {
243: PetscCheck(!svd->omega,PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"There is no support yet for generalized hyperbolic problems");
244: PetscCall(SVDSetProblemType(svd,SVD_GENERALIZED));
245: } else {
246: if (svd->omega) PetscCall(SVDSetProblemType(svd,SVD_HYPERBOLIC));
247: else PetscCall(SVDSetProblemType(svd,SVD_STANDARD));
248: }
249: } else { /* check consistency of problem type set by user */
250: if (svd->OPb) {
251: PetscCheck(svd->isgeneralized,PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_INCOMP,"Inconsistent SVD state: the problem type does not match the number of matrices");
252: PetscCheck(!svd->omega,PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"There is no support yet for generalized hyperbolic problems");
253: } else {
254: PetscCheck(!svd->isgeneralized,PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_INCOMP,"Inconsistent SVD state: the problem type does not match the number of matrices");
255: 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()");
256: else PetscCheck(!svd->ishyperbolic,PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_INCOMP,"Inconsistent SVD state: a hyperbolic problem requires passing a signature with SVDSetSignature()");
257: }
258: }
260: /* determine how to handle the transpose */
261: svd->expltrans = PETSC_TRUE;
262: if (svd->impltrans) svd->expltrans = PETSC_FALSE;
263: else {
264: PetscCall(MatHasOperation(svd->OP,MATOP_TRANSPOSE,&flg));
265: if (!flg) svd->expltrans = PETSC_FALSE;
266: else {
267: PetscCall(PetscObjectTypeCompareAny((PetscObject)svd,&flg,SVDLAPACK,SVDSCALAPACK,SVDKSVD,SVDELEMENTAL,""));
268: if (flg) svd->expltrans = PETSC_FALSE;
269: }
270: }
272: /* get matrix dimensions */
273: PetscCall(MatGetSize(svd->OP,&M,&N));
274: if (svd->isgeneralized) {
275: PetscCall(MatGetSize(svd->OPb,&P,NULL));
276: PetscCheck(M+P>=N,PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"The case when [A;B] has less rows than columns is not supported");
277: }
279: /* build transpose matrix */
280: PetscCall(MatDestroy(&svd->A));
281: PetscCall(MatDestroy(&svd->AT));
282: PetscCall(PetscObjectReference((PetscObject)svd->OP));
283: if (svd->expltrans) {
284: if (svd->isgeneralized || M>=N) {
285: svd->A = svd->OP;
286: PetscCall(MatHermitianTranspose(svd->OP,MAT_INITIAL_MATRIX,&svd->AT));
287: } else {
288: PetscCall(MatHermitianTranspose(svd->OP,MAT_INITIAL_MATRIX,&svd->A));
289: svd->AT = svd->OP;
290: }
291: } else {
292: if (svd->isgeneralized || M>=N) {
293: svd->A = svd->OP;
294: PetscCall(MatCreateHermitianTranspose(svd->OP,&svd->AT));
295: } else {
296: PetscCall(MatCreateHermitianTranspose(svd->OP,&svd->A));
297: svd->AT = svd->OP;
298: }
299: }
301: /* build transpose matrix B for GSVD */
302: if (svd->isgeneralized) {
303: PetscCall(MatDestroy(&svd->B));
304: PetscCall(MatDestroy(&svd->BT));
305: PetscCall(PetscObjectReference((PetscObject)svd->OPb));
306: if (svd->expltrans) {
307: svd->B = svd->OPb;
308: PetscCall(MatHermitianTranspose(svd->OPb,MAT_INITIAL_MATRIX,&svd->BT));
309: } else {
310: svd->B = svd->OPb;
311: PetscCall(MatCreateHermitianTranspose(svd->OPb,&svd->BT));
312: }
313: }
315: if (!svd->isgeneralized && M<N) {
316: /* swap initial vectors */
317: if (svd->nini || svd->ninil) {
318: T=svd->ISL; svd->ISL=svd->IS; svd->IS=T;
319: k=svd->ninil; svd->ninil=svd->nini; svd->nini=k;
320: }
321: /* swap basis vectors */
322: if (!svd->swapped) { /* only the first time in case of multiple calls */
323: bv=svd->V; svd->V=svd->U; svd->U=bv;
324: svd->swapped = PETSC_TRUE;
325: }
326: }
328: maxnsol = svd->isgeneralized? PetscMin(PetscMin(M,N),P): PetscMin(M,N);
329: svd->ncv = PetscMin(svd->ncv,maxnsol);
330: svd->nsv = PetscMin(svd->nsv,maxnsol);
331: PetscCheck(svd->ncv==PETSC_DEFAULT || svd->nsv<=svd->ncv,PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_OUTOFRANGE,"nsv bigger than ncv");
333: /* relative convergence criterion is not allowed in GSVD */
334: if (svd->conv==(SVDConv)-1) PetscCall(SVDSetConvergenceTest(svd,svd->isgeneralized?SVD_CONV_NORM:SVD_CONV_REL));
335: PetscCheck(!svd->isgeneralized || svd->conv!=SVD_CONV_REL,PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"Relative convergence criterion is not allowed in GSVD");
337: /* initialization of matrix norm (stardard case only, for GSVD it is done inside setup()) */
338: if (!svd->isgeneralized && svd->conv==SVD_CONV_NORM && !svd->nrma) PetscCall(MatNorm(svd->OP,NORM_INFINITY,&svd->nrma));
340: /* call specific solver setup */
341: PetscUseTypeMethod(svd,setup);
343: /* set tolerance if not yet set */
344: if (svd->tol==(PetscReal)PETSC_DEFAULT) svd->tol = SLEPC_DEFAULT_TOL;
346: /* fill sorting criterion context */
347: PetscCall(DSGetSlepcSC(svd->ds,&sc));
348: sc->comparison = (svd->which==SVD_LARGEST)? SlepcCompareLargestReal: SlepcCompareSmallestReal;
349: sc->comparisonctx = NULL;
350: sc->map = NULL;
351: sc->mapobj = NULL;
353: /* process initial vectors */
354: if (svd->nini<0) {
355: k = -svd->nini;
356: PetscCheck(k<=svd->ncv,PetscObjectComm((PetscObject)svd),PETSC_ERR_USER_INPUT,"The number of initial vectors is larger than ncv");
357: PetscCall(BVInsertVecs(svd->V,0,&k,svd->IS,PETSC_TRUE));
358: PetscCall(SlepcBasisDestroy_Private(&svd->nini,&svd->IS));
359: svd->nini = k;
360: }
361: if (svd->ninil<0) {
362: k = 0;
363: if (svd->leftbasis) {
364: k = -svd->ninil;
365: PetscCheck(k<=svd->ncv,PetscObjectComm((PetscObject)svd),PETSC_ERR_USER_INPUT,"The number of left initial vectors is larger than ncv");
366: PetscCall(BVInsertVecs(svd->U,0,&k,svd->ISL,PETSC_TRUE));
367: } else PetscCall(PetscInfo(svd,"Ignoring initial left vectors\n"));
368: PetscCall(SlepcBasisDestroy_Private(&svd->ninil,&svd->ISL));
369: svd->ninil = k;
370: }
372: PetscCall(PetscLogEventEnd(SVD_SetUp,svd,0,0,0));
373: svd->state = SVD_STATE_SETUP;
374: PetscFunctionReturn(PETSC_SUCCESS);
375: }
377: /*@C
378: SVDSetInitialSpaces - Specify two basis of vectors that constitute the initial
379: right and/or left spaces.
381: Collective
383: Input Parameters:
384: + svd - the singular value solver context
385: . nr - number of right vectors
386: . isr - set of basis vectors of the right initial space
387: . nl - number of left vectors
388: - isl - set of basis vectors of the left initial space
390: Notes:
391: The initial right and left spaces are rough approximations to the right and/or
392: left singular subspaces from which the solver starts to iterate.
393: It is not necessary to provide both sets of vectors.
395: Some solvers start to iterate on a single vector (initial vector). In that case,
396: the other vectors are ignored.
398: These vectors do not persist from one SVDSolve() call to the other, so the
399: initial space should be set every time.
401: The vectors do not need to be mutually orthonormal, since they are explicitly
402: orthonormalized internally.
404: Common usage of this function is when the user can provide a rough approximation
405: of the wanted singular space. Then, convergence may be faster.
407: Level: intermediate
409: .seealso: SVDSetUp()
410: @*/
411: PetscErrorCode SVDSetInitialSpaces(SVD svd,PetscInt nr,Vec isr[],PetscInt nl,Vec isl[])
412: {
413: PetscFunctionBegin;
417: PetscCheck(nr>=0,PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_OUTOFRANGE,"Argument nr cannot be negative");
418: if (nr>0) {
419: PetscAssertPointer(isr,3);
421: }
422: PetscCheck(nl>=0,PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_OUTOFRANGE,"Argument nl cannot be negative");
423: if (nl>0) {
424: PetscAssertPointer(isl,5);
426: }
427: PetscCall(SlepcBasisReference_Private(nr,isr,&svd->nini,&svd->IS));
428: PetscCall(SlepcBasisReference_Private(nl,isl,&svd->ninil,&svd->ISL));
429: if (nr>0 || nl>0) svd->state = SVD_STATE_INITIAL;
430: PetscFunctionReturn(PETSC_SUCCESS);
431: }
433: /*
434: SVDSetDimensions_Default - Set reasonable values for ncv, mpd if not set
435: by the user. This is called at setup.
436: */
437: PetscErrorCode SVDSetDimensions_Default(SVD svd)
438: {
439: PetscInt N,M,P,maxnsol;
441: PetscFunctionBegin;
442: PetscCall(MatGetSize(svd->OP,&M,&N));
443: maxnsol = PetscMin(M,N);
444: if (svd->isgeneralized) {
445: PetscCall(MatGetSize(svd->OPb,&P,NULL));
446: maxnsol = PetscMin(maxnsol,P);
447: }
448: if (svd->ncv!=PETSC_DEFAULT) { /* ncv set */
449: PetscCheck(svd->ncv>=svd->nsv,PetscObjectComm((PetscObject)svd),PETSC_ERR_USER_INPUT,"The value of ncv must be at least nsv");
450: } else if (svd->mpd!=PETSC_DEFAULT) { /* mpd set */
451: svd->ncv = PetscMin(maxnsol,svd->nsv+svd->mpd);
452: } else { /* neither set: defaults depend on nsv being small or large */
453: if (svd->nsv<500) svd->ncv = PetscMin(maxnsol,PetscMax(2*svd->nsv,10));
454: else {
455: svd->mpd = 500;
456: svd->ncv = PetscMin(maxnsol,svd->nsv+svd->mpd);
457: }
458: }
459: if (svd->mpd==PETSC_DEFAULT) svd->mpd = svd->ncv;
460: PetscFunctionReturn(PETSC_SUCCESS);
461: }
463: /*@
464: SVDAllocateSolution - Allocate memory storage for common variables such
465: as the singular values and the basis vectors.
467: Collective
469: Input Parameters:
470: + svd - eigensolver context
471: - extra - number of additional positions, used for methods that require a
472: working basis slightly larger than ncv
474: Developer Notes:
475: This is SLEPC_EXTERN because it may be required by user plugin SVD
476: implementations.
478: This is called at setup after setting the value of ncv and the flag leftbasis.
480: Level: developer
482: .seealso: SVDSetUp()
483: @*/
484: PetscErrorCode SVDAllocateSolution(SVD svd,PetscInt extra)
485: {
486: PetscInt oldsize,requested;
487: Vec tr,tl;
489: PetscFunctionBegin;
490: requested = svd->ncv + extra;
492: /* oldsize is zero if this is the first time setup is called */
493: PetscCall(BVGetSizes(svd->V,NULL,NULL,&oldsize));
495: /* allocate sigma */
496: if (requested != oldsize || !svd->sigma) {
497: PetscCall(PetscFree3(svd->sigma,svd->perm,svd->errest));
498: if (svd->sign) PetscCall(PetscFree(svd->sign));
499: PetscCall(PetscMalloc3(requested,&svd->sigma,requested,&svd->perm,requested,&svd->errest));
500: if (svd->ishyperbolic) PetscCall(PetscMalloc1(requested,&svd->sign));
501: }
502: /* allocate V */
503: if (!svd->V) PetscCall(SVDGetBV(svd,&svd->V,NULL));
504: if (!oldsize) {
505: if (!((PetscObject)(svd->V))->type_name) PetscCall(BVSetType(svd->V,BVMAT));
506: PetscCall(MatCreateVecsEmpty(svd->A,&tr,NULL));
507: PetscCall(BVSetSizesFromVec(svd->V,tr,requested));
508: PetscCall(VecDestroy(&tr));
509: } else PetscCall(BVResize(svd->V,requested,PETSC_FALSE));
510: /* allocate U */
511: if (svd->leftbasis && !svd->isgeneralized) {
512: if (!svd->U) PetscCall(SVDGetBV(svd,NULL,&svd->U));
513: if (!oldsize) {
514: if (!((PetscObject)(svd->U))->type_name) PetscCall(BVSetType(svd->U,((PetscObject)(svd->V))->type_name));
515: PetscCall(MatCreateVecsEmpty(svd->A,NULL,&tl));
516: PetscCall(BVSetSizesFromVec(svd->U,tl,requested));
517: PetscCall(VecDestroy(&tl));
518: } else PetscCall(BVResize(svd->U,requested,PETSC_FALSE));
519: } else if (svd->isgeneralized) { /* left basis for the GSVD */
520: if (!svd->U) PetscCall(SVDGetBV(svd,NULL,&svd->U));
521: if (!oldsize) {
522: if (!((PetscObject)(svd->U))->type_name) PetscCall(BVSetType(svd->U,((PetscObject)(svd->V))->type_name));
523: PetscCall(SVDCreateLeftTemplate(svd,&tl));
524: PetscCall(BVSetSizesFromVec(svd->U,tl,requested));
525: PetscCall(VecDestroy(&tl));
526: } else PetscCall(BVResize(svd->U,requested,PETSC_FALSE));
527: }
528: PetscFunctionReturn(PETSC_SUCCESS);
529: }