1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2021, 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 on svd
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: {
33: PetscInt Ma,Na,Mb,Nb,ma,na,mb,nb,M0,N0,m0,n0;
34: PetscBool samesize=PETSC_TRUE;
43: /* Check matrix sizes */
44: MatGetSize(A,&Ma,&Na);
45: MatGetLocalSize(A,&ma,&na);
46: if (svd->OP) {
47: MatGetSize(svd->OP,&M0,&N0);
48: MatGetLocalSize(svd->OP,&m0,&n0);
49: if (M0!=Ma || N0!=Na || m0!=ma || n0!=na) samesize = PETSC_FALSE;
50: }
51: if (B) {
52: MatGetSize(B,&Mb,&Nb);
53: MatGetLocalSize(B,&mb,&nb);
54: if (Na!=Nb) SETERRQ2(PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_WRONG,"Different number of columns in A (%D) and B (%D)",Na,Nb);
55: if (na!=nb) SETERRQ2(PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_WRONG,"Different local column size in A (%D) and B (%D)",na,nb);
56: if (svd->OPb) {
57: MatGetSize(svd->OPb,&M0,&N0);
58: MatGetLocalSize(svd->OPb,&m0,&n0);
59: if (M0!=Mb || N0!=Nb || m0!=mb || n0!=nb) samesize = PETSC_FALSE;
60: }
61: }
63: PetscObjectReference((PetscObject)A);
64: if (B) { PetscObjectReference((PetscObject)B); }
65: if (svd->state && !samesize) {
66: SVDReset(svd);
67: } else {
68: MatDestroy(&svd->OP);
69: MatDestroy(&svd->OPb);
70: MatDestroy(&svd->A);
71: MatDestroy(&svd->B);
72: MatDestroy(&svd->AT);
73: MatDestroy(&svd->BT);
74: }
75: svd->nrma = 0.0;
76: svd->nrmb = 0.0;
77: svd->OP = A;
78: svd->OPb = B;
79: svd->state = SVD_STATE_INITIAL;
80: return(0);
81: }
83: /*@
84: SVDGetOperators - Get the matrices associated with the singular value problem.
86: Collective on svd
88: Input Parameter:
89: . svd - the singular value solver context
91: Output Parameters:
92: + A - the matrix associated with the singular value problem
93: - B - the second matrix in the case of GSVD
95: Level: intermediate
97: .seealso: SVDSolve(), SVDSetOperators()
98: @*/
99: PetscErrorCode SVDGetOperators(SVD svd,Mat *A,Mat *B)100: {
103: if (A) *A = svd->OP;
104: if (B) *B = svd->OPb;
105: return(0);
106: }
108: /*@
109: SVDSetUp - Sets up all the internal data structures necessary for the
110: execution of the singular value solver.
112: Collective on svd
114: Input Parameter:
115: . svd - singular value solver context
117: Notes:
118: This function need not be called explicitly in most cases, since SVDSolve()
119: calls it. It can be useful when one wants to measure the set-up time
120: separately from the solve time.
122: Level: developer
124: .seealso: SVDCreate(), SVDSolve(), SVDDestroy()
125: @*/
126: PetscErrorCode SVDSetUp(SVD svd)127: {
129: PetscBool flg;
130: PetscInt M,N,P=0,k,maxnsol;
131: SlepcSC sc;
132: Vec *T;
133: BV bv;
137: if (svd->state) return(0);
138: PetscLogEventBegin(SVD_SetUp,svd,0,0,0);
140: /* reset the convergence flag from the previous solves */
141: svd->reason = SVD_CONVERGED_ITERATING;
143: /* Set default solver type (SVDSetFromOptions was not called) */
144: if (!((PetscObject)svd)->type_name) {
145: SVDSetType(svd,SVDCROSS);
146: }
147: if (!svd->ds) { SVDGetDS(svd,&svd->ds); }
149: /* check matrices */
150: if (!svd->OP) SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_WRONGSTATE,"SVDSetOperators() must be called first");
152: /* Set default problem type */
153: if (!svd->problem_type) {
154: if (svd->OPb) {
155: SVDSetProblemType(svd,SVD_GENERALIZED);
156: } else {
157: SVDSetProblemType(svd,SVD_STANDARD);
158: }
159: } else if (!svd->OPb && svd->isgeneralized) {
160: PetscInfo(svd,"Problem type set as generalized but no matrix B was provided; reverting to a standard singular value problem\n");
161: svd->isgeneralized = PETSC_FALSE;
162: svd->problem_type = SVD_STANDARD;
163: } else if (svd->OPb && !svd->isgeneralized) SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_INCOMP,"Inconsistent SVD state: the problem type does not match the number of matrices");
165: /* determine how to handle the transpose */
166: svd->expltrans = PETSC_TRUE;
167: if (svd->impltrans) svd->expltrans = PETSC_FALSE;
168: else {
169: MatHasOperation(svd->OP,MATOP_TRANSPOSE,&flg);
170: if (!flg) svd->expltrans = PETSC_FALSE;
171: else {
172: PetscObjectTypeCompareAny((PetscObject)svd,&flg,SVDLAPACK,SVDSCALAPACK,SVDELEMENTAL,"");
173: if (flg) svd->expltrans = PETSC_FALSE;
174: }
175: }
177: /* get matrix dimensions */
178: MatGetSize(svd->OP,&M,&N);
179: if (svd->isgeneralized) {
180: MatGetSize(svd->OPb,&P,NULL);
181: if (M+P<N) SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"The case when [A;B] has less rows than columns is not supported");
182: }
184: /* build transpose matrix */
185: MatDestroy(&svd->A);
186: MatDestroy(&svd->AT);
187: PetscObjectReference((PetscObject)svd->OP);
188: if (svd->expltrans) {
189: if (svd->isgeneralized || M>=N) {
190: svd->A = svd->OP;
191: MatHermitianTranspose(svd->OP,MAT_INITIAL_MATRIX,&svd->AT);
192: } else {
193: MatHermitianTranspose(svd->OP,MAT_INITIAL_MATRIX,&svd->A);
194: svd->AT = svd->OP;
195: }
196: } else {
197: if (svd->isgeneralized || M>=N) {
198: svd->A = svd->OP;
199: MatCreateHermitianTranspose(svd->OP,&svd->AT);
200: } else {
201: MatCreateHermitianTranspose(svd->OP,&svd->A);
202: svd->AT = svd->OP;
203: }
204: }
206: /* build transpose matrix B for GSVD */
207: if (svd->isgeneralized) {
208: MatDestroy(&svd->B);
209: MatDestroy(&svd->BT);
210: PetscObjectReference((PetscObject)svd->OPb);
211: if (svd->expltrans) {
212: svd->B = svd->OPb;
213: MatHermitianTranspose(svd->OPb,MAT_INITIAL_MATRIX,&svd->BT);
214: } else {
215: svd->B = svd->OPb;
216: MatCreateHermitianTranspose(svd->OPb,&svd->BT);
217: }
218: }
220: if (!svd->isgeneralized && M<N) {
221: /* swap initial vectors */
222: if (svd->nini || svd->ninil) {
223: T=svd->ISL; svd->ISL=svd->IS; svd->IS=T;
224: k=svd->ninil; svd->ninil=svd->nini; svd->nini=k;
225: }
226: /* swap basis vectors */
227: if (!svd->swapped) { /* only the first time in case of multiple calls */
228: bv=svd->V; svd->V=svd->U; svd->U=bv;
229: svd->swapped = PETSC_TRUE;
230: }
231: }
233: maxnsol = svd->isgeneralized? PetscMin(PetscMin(M,N),P): PetscMin(M,N);
234: svd->ncv = PetscMin(svd->ncv,maxnsol);
235: svd->nsv = PetscMin(svd->nsv,maxnsol);
236: if (svd->ncv!=PETSC_DEFAULT && svd->nsv > svd->ncv) SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_OUTOFRANGE,"nsv bigger than ncv");
238: /* initialization of matrix norms */
239: if (svd->conv==SVD_CONV_NORM) {
240: if (!svd->nrma) {
241: MatNorm(svd->OP,NORM_INFINITY,&svd->nrma);
242: }
243: if (svd->isgeneralized && !svd->nrmb) {
244: MatNorm(svd->OPb,NORM_INFINITY,&svd->nrmb);
245: }
246: }
248: /* call specific solver setup */
249: (*svd->ops->setup)(svd);
251: /* set tolerance if not yet set */
252: if (svd->tol==PETSC_DEFAULT) svd->tol = SLEPC_DEFAULT_TOL;
254: /* fill sorting criterion context */
255: DSGetSlepcSC(svd->ds,&sc);
256: sc->comparison = (svd->which==SVD_LARGEST)? SlepcCompareLargestReal: SlepcCompareSmallestReal;
257: sc->comparisonctx = NULL;
258: sc->map = NULL;
259: sc->mapobj = NULL;
261: /* process initial vectors */
262: if (svd->nini<0) {
263: k = -svd->nini;
264: if (k>svd->ncv) SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_USER_INPUT,"The number of initial vectors is larger than ncv");
265: BVInsertVecs(svd->V,0,&k,svd->IS,PETSC_TRUE);
266: SlepcBasisDestroy_Private(&svd->nini,&svd->IS);
267: svd->nini = k;
268: }
269: if (svd->ninil<0) {
270: k = 0;
271: if (svd->leftbasis) {
272: k = -svd->ninil;
273: if (k>svd->ncv) SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_USER_INPUT,"The number of left initial vectors is larger than ncv");
274: BVInsertVecs(svd->U,0,&k,svd->ISL,PETSC_TRUE);
275: } else {
276: PetscInfo(svd,"Ignoring initial left vectors\n");
277: }
278: SlepcBasisDestroy_Private(&svd->ninil,&svd->ISL);
279: svd->ninil = k;
280: }
282: PetscLogEventEnd(SVD_SetUp,svd,0,0,0);
283: svd->state = SVD_STATE_SETUP;
284: return(0);
285: }
287: /*@C
288: SVDSetInitialSpaces - Specify two basis of vectors that constitute the initial
289: right and/or left spaces, that is, a rough approximation to the right and/or
290: left singular subspaces from which the solver starts to iterate.
292: Collective on svd
294: Input Parameters:
295: + svd - the singular value solver context
296: . nr - number of right vectors
297: . isr - set of basis vectors of the right initial space
298: . nl - number of left vectors
299: - isl - set of basis vectors of the left initial space
301: Notes:
302: It is not necessary to provide both sets of vectors.
304: Some solvers start to iterate on a single vector (initial vector). In that case,
305: the other vectors are ignored.
307: These vectors do not persist from one SVDSolve() call to the other, so the
308: initial space should be set every time.
310: The vectors do not need to be mutually orthonormal, since they are explicitly
311: orthonormalized internally.
313: Common usage of this function is when the user can provide a rough approximation
314: of the wanted singular space. Then, convergence may be faster.
316: Level: intermediate
317: @*/
318: PetscErrorCode SVDSetInitialSpaces(SVD svd,PetscInt nr,Vec isr[],PetscInt nl,Vec isl[])319: {
326: if (nr<0) SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_OUTOFRANGE,"Argument nr cannot be negative");
327: if (nr>0) {
330: }
331: if (nl<0) SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_OUTOFRANGE,"Argument nl cannot be negative");
332: if (nl>0) {
335: }
336: SlepcBasisReference_Private(nr,isr,&svd->nini,&svd->IS);
337: SlepcBasisReference_Private(nl,isl,&svd->ninil,&svd->ISL);
338: if (nr>0 || nl>0) svd->state = SVD_STATE_INITIAL;
339: return(0);
340: }
342: /*
343: SVDSetDimensions_Default - Set reasonable values for ncv, mpd if not set
344: by the user. This is called at setup.
345: */
346: PetscErrorCode SVDSetDimensions_Default(SVD svd)347: {
349: PetscInt N,M,P,maxnsol;
352: MatGetSize(svd->OP,&M,&N);
353: maxnsol = PetscMin(M,N);
354: if (svd->isgeneralized) {
355: MatGetSize(svd->OPb,&P,NULL);
356: maxnsol = PetscMin(maxnsol,P);
357: }
358: if (svd->ncv!=PETSC_DEFAULT) { /* ncv set */
359: if (svd->ncv<svd->nsv) SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_USER_INPUT,"The value of ncv must be at least nsv");
360: } else if (svd->mpd!=PETSC_DEFAULT) { /* mpd set */
361: svd->ncv = PetscMin(maxnsol,svd->nsv+svd->mpd);
362: } else { /* neither set: defaults depend on nsv being small or large */
363: if (svd->nsv<500) svd->ncv = PetscMin(maxnsol,PetscMax(2*svd->nsv,10));
364: else {
365: svd->mpd = 500;
366: svd->ncv = PetscMin(maxnsol,svd->nsv+svd->mpd);
367: }
368: }
369: if (svd->mpd==PETSC_DEFAULT) svd->mpd = svd->ncv;
370: return(0);
371: }
373: /*@
374: SVDAllocateSolution - Allocate memory storage for common variables such
375: as the singular values and the basis vectors.
377: Collective on svd
379: Input Parameters:
380: + svd - eigensolver context
381: - extra - number of additional positions, used for methods that require a
382: working basis slightly larger than ncv
384: Developers Notes:
385: This is SLEPC_EXTERN because it may be required by user plugin SVD386: implementations.
388: This is called at setup after setting the value of ncv and the flag leftbasis.
390: Level: developer
391: @*/
392: PetscErrorCode SVDAllocateSolution(SVD svd,PetscInt extra)393: {
395: PetscInt oldsize,requested;
396: Vec tr,tl;
399: requested = svd->ncv + extra;
401: /* oldsize is zero if this is the first time setup is called */
402: BVGetSizes(svd->V,NULL,NULL,&oldsize);
404: /* allocate sigma */
405: if (requested != oldsize || !svd->sigma) {
406: PetscFree3(svd->sigma,svd->perm,svd->errest);
407: PetscMalloc3(requested,&svd->sigma,requested,&svd->perm,requested,&svd->errest);
408: PetscLogObjectMemory((PetscObject)svd,PetscMax(0,requested-oldsize)*(2*sizeof(PetscReal)+sizeof(PetscInt)));
409: }
410: /* allocate V */
411: if (!svd->V) { SVDGetBV(svd,&svd->V,NULL); }
412: if (!oldsize) {
413: if (!((PetscObject)(svd->V))->type_name) {
414: BVSetType(svd->V,BVSVEC);
415: }
416: MatCreateVecsEmpty(svd->A,&tr,NULL);
417: BVSetSizesFromVec(svd->V,tr,requested);
418: VecDestroy(&tr);
419: } else {
420: BVResize(svd->V,requested,PETSC_FALSE);
421: }
422: /* allocate U */
423: if (svd->leftbasis && !svd->isgeneralized) {
424: if (!svd->U) { SVDGetBV(svd,NULL,&svd->U); }
425: if (!oldsize) {
426: if (!((PetscObject)(svd->U))->type_name) {
427: BVSetType(svd->U,BVSVEC);
428: }
429: MatCreateVecsEmpty(svd->A,NULL,&tl);
430: BVSetSizesFromVec(svd->U,tl,requested);
431: VecDestroy(&tl);
432: } else {
433: BVResize(svd->U,requested,PETSC_FALSE);
434: }
435: } else if (svd->isgeneralized) { /* left basis for the GSVD */
436: if (!svd->U) { SVDGetBV(svd,NULL,&svd->U); }
437: if (!oldsize) {
438: if (!((PetscObject)(svd->U))->type_name) {
439: BVSetType(svd->U,BVSVEC);
440: }
441: SVDCreateLeftTemplate(svd,&tl);
442: BVSetSizesFromVec(svd->U,tl,requested);
443: VecDestroy(&tl);
444: } else {
445: BVResize(svd->U,requested,PETSC_FALSE);
446: }
447: }
448: return(0);
449: }