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 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: {
32: PetscInt Ma,Na,Mb,Nb,ma,na,mb,nb,M0,N0,m0,n0;
33: PetscBool samesize=PETSC_TRUE;
41: /* Check matrix sizes */
42: MatGetSize(A,&Ma,&Na);
43: MatGetLocalSize(A,&ma,&na);
44: if (svd->OP) {
45: MatGetSize(svd->OP,&M0,&N0);
46: MatGetLocalSize(svd->OP,&m0,&n0);
47: if (M0!=Ma || N0!=Na || m0!=ma || n0!=na) samesize = PETSC_FALSE;
48: }
49: if (B) {
50: MatGetSize(B,&Mb,&Nb);
51: MatGetLocalSize(B,&mb,&nb);
54: if (svd->OPb) {
55: MatGetSize(svd->OPb,&M0,&N0);
56: MatGetLocalSize(svd->OPb,&m0,&n0);
57: if (M0!=Mb || N0!=Nb || m0!=mb || n0!=nb) samesize = PETSC_FALSE;
58: }
59: }
61: PetscObjectReference((PetscObject)A);
62: if (B) PetscObjectReference((PetscObject)B);
63: if (svd->state && !samesize) SVDReset(svd);
64: else {
65: MatDestroy(&svd->OP);
66: MatDestroy(&svd->OPb);
67: MatDestroy(&svd->A);
68: MatDestroy(&svd->B);
69: MatDestroy(&svd->AT);
70: MatDestroy(&svd->BT);
71: }
72: svd->nrma = 0.0;
73: svd->nrmb = 0.0;
74: svd->OP = A;
75: svd->OPb = B;
76: svd->state = SVD_STATE_INITIAL;
77: PetscFunctionReturn(0);
78: }
80: /*@
81: SVDGetOperators - Get the matrices associated with the singular value problem.
83: Collective on svd
85: Input Parameter:
86: . svd - the singular value solver context
88: Output Parameters:
89: + A - the matrix associated with the singular value problem
90: - B - the second matrix in the case of GSVD
92: Level: intermediate
94: .seealso: SVDSolve(), SVDSetOperators()
95: @*/
96: PetscErrorCode SVDGetOperators(SVD svd,Mat *A,Mat *B) 97: {
99: if (A) *A = svd->OP;
100: if (B) *B = svd->OPb;
101: PetscFunctionReturn(0);
102: }
104: /*@
105: SVDSetUp - Sets up all the internal data structures necessary for the
106: execution of the singular value solver.
108: Collective on svd
110: Input Parameter:
111: . svd - singular value solver context
113: Notes:
114: This function need not be called explicitly in most cases, since SVDSolve()
115: calls it. It can be useful when one wants to measure the set-up time
116: separately from the solve time.
118: Level: developer
120: .seealso: SVDCreate(), SVDSolve(), SVDDestroy()
121: @*/
122: PetscErrorCode SVDSetUp(SVD svd)123: {
124: PetscBool flg;
125: PetscInt M,N,P=0,k,maxnsol;
126: SlepcSC sc;
127: Vec *T;
128: BV bv;
131: if (svd->state) PetscFunctionReturn(0);
132: PetscLogEventBegin(SVD_SetUp,svd,0,0,0);
134: /* reset the convergence flag from the previous solves */
135: svd->reason = SVD_CONVERGED_ITERATING;
137: /* Set default solver type (SVDSetFromOptions was not called) */
138: if (!((PetscObject)svd)->type_name) SVDSetType(svd,SVDCROSS);
139: if (!svd->ds) SVDGetDS(svd,&svd->ds);
141: /* check matrices */
144: /* Set default problem type */
145: if (!svd->problem_type) {
146: if (svd->OPb) SVDSetProblemType(svd,SVD_GENERALIZED);
147: else SVDSetProblemType(svd,SVD_STANDARD);
148: } else if (!svd->OPb && svd->isgeneralized) {
149: PetscInfo(svd,"Problem type set as generalized but no matrix B was provided; reverting to a standard singular value problem\n");
150: svd->isgeneralized = PETSC_FALSE;
151: svd->problem_type = SVD_STANDARD;
154: /* determine how to handle the transpose */
155: svd->expltrans = PETSC_TRUE;
156: if (svd->impltrans) svd->expltrans = PETSC_FALSE;
157: else {
158: MatHasOperation(svd->OP,MATOP_TRANSPOSE,&flg);
159: if (!flg) svd->expltrans = PETSC_FALSE;
160: else {
161: PetscObjectTypeCompareAny((PetscObject)svd,&flg,SVDLAPACK,SVDSCALAPACK,SVDELEMENTAL,"");
162: if (flg) svd->expltrans = PETSC_FALSE;
163: }
164: }
166: /* get matrix dimensions */
167: MatGetSize(svd->OP,&M,&N);
168: if (svd->isgeneralized) {
169: MatGetSize(svd->OPb,&P,NULL);
171: }
173: /* build transpose matrix */
174: MatDestroy(&svd->A);
175: MatDestroy(&svd->AT);
176: PetscObjectReference((PetscObject)svd->OP);
177: if (svd->expltrans) {
178: if (svd->isgeneralized || M>=N) {
179: svd->A = svd->OP;
180: MatHermitianTranspose(svd->OP,MAT_INITIAL_MATRIX,&svd->AT);
181: } else {
182: MatHermitianTranspose(svd->OP,MAT_INITIAL_MATRIX,&svd->A);
183: svd->AT = svd->OP;
184: }
185: } else {
186: if (svd->isgeneralized || M>=N) {
187: svd->A = svd->OP;
188: MatCreateHermitianTranspose(svd->OP,&svd->AT);
189: } else {
190: MatCreateHermitianTranspose(svd->OP,&svd->A);
191: svd->AT = svd->OP;
192: }
193: }
195: /* build transpose matrix B for GSVD */
196: if (svd->isgeneralized) {
197: MatDestroy(&svd->B);
198: MatDestroy(&svd->BT);
199: PetscObjectReference((PetscObject)svd->OPb);
200: if (svd->expltrans) {
201: svd->B = svd->OPb;
202: MatHermitianTranspose(svd->OPb,MAT_INITIAL_MATRIX,&svd->BT);
203: } else {
204: svd->B = svd->OPb;
205: MatCreateHermitianTranspose(svd->OPb,&svd->BT);
206: }
207: }
209: if (!svd->isgeneralized && M<N) {
210: /* swap initial vectors */
211: if (svd->nini || svd->ninil) {
212: T=svd->ISL; svd->ISL=svd->IS; svd->IS=T;
213: k=svd->ninil; svd->ninil=svd->nini; svd->nini=k;
214: }
215: /* swap basis vectors */
216: if (!svd->swapped) { /* only the first time in case of multiple calls */
217: bv=svd->V; svd->V=svd->U; svd->U=bv;
218: svd->swapped = PETSC_TRUE;
219: }
220: }
222: maxnsol = svd->isgeneralized? PetscMin(PetscMin(M,N),P): PetscMin(M,N);
223: svd->ncv = PetscMin(svd->ncv,maxnsol);
224: svd->nsv = PetscMin(svd->nsv,maxnsol);
227: /* relative convergence criterion is not allowed in GSVD */
228: if (svd->conv==(SVDConv)-1) SVDSetConvergenceTest(svd,svd->isgeneralized?SVD_CONV_NORM:SVD_CONV_REL);
231: /* initialization of matrix norm (stardard case only, for GSVD it is done inside setup()) */
232: if (!svd->isgeneralized && svd->conv==SVD_CONV_NORM && !svd->nrma) MatNorm(svd->OP,NORM_INFINITY,&svd->nrma);
234: /* call specific solver setup */
235: (*svd->ops->setup)(svd);
237: /* set tolerance if not yet set */
238: if (svd->tol==PETSC_DEFAULT) svd->tol = SLEPC_DEFAULT_TOL;
240: /* fill sorting criterion context */
241: DSGetSlepcSC(svd->ds,&sc);
242: sc->comparison = (svd->which==SVD_LARGEST)? SlepcCompareLargestReal: SlepcCompareSmallestReal;
243: sc->comparisonctx = NULL;
244: sc->map = NULL;
245: sc->mapobj = NULL;
247: /* process initial vectors */
248: if (svd->nini<0) {
249: k = -svd->nini;
251: BVInsertVecs(svd->V,0,&k,svd->IS,PETSC_TRUE);
252: SlepcBasisDestroy_Private(&svd->nini,&svd->IS);
253: svd->nini = k;
254: }
255: if (svd->ninil<0) {
256: k = 0;
257: if (svd->leftbasis) {
258: k = -svd->ninil;
260: BVInsertVecs(svd->U,0,&k,svd->ISL,PETSC_TRUE);
261: } else PetscInfo(svd,"Ignoring initial left vectors\n");
262: SlepcBasisDestroy_Private(&svd->ninil,&svd->ISL);
263: svd->ninil = k;
264: }
266: PetscLogEventEnd(SVD_SetUp,svd,0,0,0);
267: svd->state = SVD_STATE_SETUP;
268: PetscFunctionReturn(0);
269: }
271: /*@C
272: SVDSetInitialSpaces - Specify two basis of vectors that constitute the initial
273: right and/or left spaces.
275: Collective on svd
277: Input Parameters:
278: + svd - the singular value solver context
279: . nr - number of right vectors
280: . isr - set of basis vectors of the right initial space
281: . nl - number of left vectors
282: - isl - set of basis vectors of the left initial space
284: Notes:
285: The initial right and left spaces are rough approximations to the right and/or
286: left singular subspaces from which the solver starts to iterate.
287: It is not necessary to provide both sets of vectors.
289: Some solvers start to iterate on a single vector (initial vector). In that case,
290: the other vectors are ignored.
292: These vectors do not persist from one SVDSolve() call to the other, so the
293: initial space should be set every time.
295: The vectors do not need to be mutually orthonormal, since they are explicitly
296: orthonormalized internally.
298: Common usage of this function is when the user can provide a rough approximation
299: of the wanted singular space. Then, convergence may be faster.
301: Level: intermediate
303: .seealso: SVDSetUp()
304: @*/
305: PetscErrorCode SVDSetInitialSpaces(SVD svd,PetscInt nr,Vec isr[],PetscInt nl,Vec isl[])306: {
311: if (nr>0) {
314: }
316: if (nl>0) {
319: }
320: SlepcBasisReference_Private(nr,isr,&svd->nini,&svd->IS);
321: SlepcBasisReference_Private(nl,isl,&svd->ninil,&svd->ISL);
322: if (nr>0 || nl>0) svd->state = SVD_STATE_INITIAL;
323: PetscFunctionReturn(0);
324: }
326: /*
327: SVDSetDimensions_Default - Set reasonable values for ncv, mpd if not set
328: by the user. This is called at setup.
329: */
330: PetscErrorCode SVDSetDimensions_Default(SVD svd)331: {
332: PetscInt N,M,P,maxnsol;
334: MatGetSize(svd->OP,&M,&N);
335: maxnsol = PetscMin(M,N);
336: if (svd->isgeneralized) {
337: MatGetSize(svd->OPb,&P,NULL);
338: maxnsol = PetscMin(maxnsol,P);
339: }
340: if (svd->ncv!=PETSC_DEFAULT) { /* ncv set */
342: } else if (svd->mpd!=PETSC_DEFAULT) { /* mpd set */
343: svd->ncv = PetscMin(maxnsol,svd->nsv+svd->mpd);
344: } else { /* neither set: defaults depend on nsv being small or large */
345: if (svd->nsv<500) svd->ncv = PetscMin(maxnsol,PetscMax(2*svd->nsv,10));
346: else {
347: svd->mpd = 500;
348: svd->ncv = PetscMin(maxnsol,svd->nsv+svd->mpd);
349: }
350: }
351: if (svd->mpd==PETSC_DEFAULT) svd->mpd = svd->ncv;
352: PetscFunctionReturn(0);
353: }
355: /*@
356: SVDAllocateSolution - Allocate memory storage for common variables such
357: as the singular values and the basis vectors.
359: Collective on svd
361: Input Parameters:
362: + svd - eigensolver context
363: - extra - number of additional positions, used for methods that require a
364: working basis slightly larger than ncv
366: Developer Notes:
367: This is SLEPC_EXTERN because it may be required by user plugin SVD368: implementations.
370: This is called at setup after setting the value of ncv and the flag leftbasis.
372: Level: developer
374: .seealso: SVDSetUp()
375: @*/
376: PetscErrorCode SVDAllocateSolution(SVD svd,PetscInt extra)377: {
378: PetscInt oldsize,requested;
379: Vec tr,tl;
381: requested = svd->ncv + extra;
383: /* oldsize is zero if this is the first time setup is called */
384: BVGetSizes(svd->V,NULL,NULL,&oldsize);
386: /* allocate sigma */
387: if (requested != oldsize || !svd->sigma) {
388: PetscFree3(svd->sigma,svd->perm,svd->errest);
389: PetscMalloc3(requested,&svd->sigma,requested,&svd->perm,requested,&svd->errest);
390: PetscLogObjectMemory((PetscObject)svd,PetscMax(0,requested-oldsize)*(2*sizeof(PetscReal)+sizeof(PetscInt)));
391: }
392: /* allocate V */
393: if (!svd->V) SVDGetBV(svd,&svd->V,NULL);
394: if (!oldsize) {
395: if (!((PetscObject)(svd->V))->type_name) BVSetType(svd->V,BVSVEC);
396: MatCreateVecsEmpty(svd->A,&tr,NULL);
397: BVSetSizesFromVec(svd->V,tr,requested);
398: VecDestroy(&tr);
399: } else BVResize(svd->V,requested,PETSC_FALSE);
400: /* allocate U */
401: if (svd->leftbasis && !svd->isgeneralized) {
402: if (!svd->U) SVDGetBV(svd,NULL,&svd->U);
403: if (!oldsize) {
404: if (!((PetscObject)(svd->U))->type_name) BVSetType(svd->U,((PetscObject)(svd->V))->type_name);
405: MatCreateVecsEmpty(svd->A,NULL,&tl);
406: BVSetSizesFromVec(svd->U,tl,requested);
407: VecDestroy(&tl);
408: } else BVResize(svd->U,requested,PETSC_FALSE);
409: } else if (svd->isgeneralized) { /* left basis for the GSVD */
410: if (!svd->U) SVDGetBV(svd,NULL,&svd->U);
411: if (!oldsize) {
412: if (!((PetscObject)(svd->U))->type_name) BVSetType(svd->U,((PetscObject)(svd->V))->type_name);
413: SVDCreateLeftTemplate(svd,&tl);
414: BVSetSizesFromVec(svd->U,tl,requested);
415: VecDestroy(&tl);
416: } else BVResize(svd->U,requested,PETSC_FALSE);
417: }
418: PetscFunctionReturn(0);
419: }