Actual source code: svdsetup.c

slepc-3.11.2 2019-07-30
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2019, 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>      /*I "slepcsvd.h" I*/

 16: /*@
 17:    SVDSetOperator - Set the matrix associated with the singular value problem.

 19:    Collective on SVD and Mat

 21:    Input Parameters:
 22: +  svd - the singular value solver context
 23: -  A  - the matrix associated with the singular value problem

 25:    Level: beginner

 27: .seealso: SVDSolve(), SVDGetOperator()
 28: @*/
 29: PetscErrorCode SVDSetOperator(SVD svd,Mat mat)
 30: {

 37:   PetscObjectReference((PetscObject)mat);
 38:   if (svd->state) { SVDReset(svd); }
 39:   else { MatDestroy(&svd->OP); }
 40:   svd->OP = mat;
 41:   svd->state = SVD_STATE_INITIAL;
 42:   return(0);
 43: }

 45: /*@
 46:    SVDGetOperator - Get the matrix associated with the singular value problem.

 48:    Not collective, though parallel Mats are returned if the SVD is parallel

 50:    Input Parameter:
 51: .  svd - the singular value solver context

 53:    Output Parameters:
 54: .  A    - the matrix associated with the singular value problem

 56:    Level: advanced

 58: .seealso: SVDSolve(), SVDSetOperator()
 59: @*/
 60: PetscErrorCode SVDGetOperator(SVD svd,Mat *A)
 61: {
 65:   *A = svd->OP;
 66:   return(0);
 67: }

 69: /*@
 70:    SVDSetUp - Sets up all the internal data structures necessary for the
 71:    execution of the singular value solver.

 73:    Collective on SVD

 75:    Input Parameter:
 76: .  svd   - singular value solver context

 78:    Notes:
 79:    This function need not be called explicitly in most cases, since SVDSolve()
 80:    calls it. It can be useful when one wants to measure the set-up time
 81:    separately from the solve time.

 83:    Level: developer

 85: .seealso: SVDCreate(), SVDSolve(), SVDDestroy()
 86: @*/
 87: PetscErrorCode SVDSetUp(SVD svd)
 88: {
 90:   PetscBool      expltrans,flg;
 91:   PetscInt       M,N,k;
 92:   SlepcSC        sc;
 93:   Vec            *T;

 97:   if (svd->state) return(0);
 98:   PetscLogEventBegin(SVD_SetUp,svd,0,0,0);

100:   /* reset the convergence flag from the previous solves */
101:   svd->reason = SVD_CONVERGED_ITERATING;

103:   /* Set default solver type (SVDSetFromOptions was not called) */
104:   if (!((PetscObject)svd)->type_name) {
105:     SVDSetType(svd,SVDCROSS);
106:   }
107:   if (!svd->ds) { SVDGetDS(svd,&svd->ds); }

109:   /* check matrix */
110:   if (!svd->OP) SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_WRONGSTATE,"SVDSetOperator must be called first");

112:   /* determine how to handle the transpose */
113:   expltrans = PETSC_TRUE;
114:   if (svd->impltrans) expltrans = PETSC_FALSE;
115:   else {
116:     MatHasOperation(svd->OP,MATOP_TRANSPOSE,&flg);
117:     if (!flg) expltrans = PETSC_FALSE;
118:     else {
119:       PetscObjectTypeCompare((PetscObject)svd,SVDLAPACK,&flg);
120:       if (flg) expltrans = PETSC_FALSE;
121:     }
122:   }

124:   /* build transpose matrix */
125:   MatDestroy(&svd->A);
126:   MatDestroy(&svd->AT);
127:   MatGetSize(svd->OP,&M,&N);
128:   PetscObjectReference((PetscObject)svd->OP);
129:   if (expltrans) {
130:     if (M>=N) {
131:       svd->A = svd->OP;
132:       MatTranspose(svd->OP,MAT_INITIAL_MATRIX,&svd->AT);
133:       MatConjugate(svd->AT);
134:     } else {
135:       MatTranspose(svd->OP,MAT_INITIAL_MATRIX,&svd->A);
136:       MatConjugate(svd->A);
137:       svd->AT = svd->OP;
138:     }
139:   } else {
140:     if (M>=N) {
141:       svd->A = svd->OP;
142:       svd->AT = NULL;
143:     } else {
144:       svd->A = NULL;
145:       svd->AT = svd->OP;
146:     }
147:   }

149:   /* swap initial vectors if necessary */
150:   if (M<N) {
151:     T=svd->ISL; svd->ISL=svd->IS; svd->IS=T;
152:     k=svd->ninil; svd->ninil=svd->nini; svd->nini=k;
153:   }

155:   if (svd->ncv > PetscMin(M,N)) svd->ncv = PetscMin(M,N);
156:   if (svd->nsv > PetscMin(M,N)) svd->nsv = PetscMin(M,N);
157:   if (svd->ncv && svd->nsv > svd->ncv) SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_OUTOFRANGE,"nsv bigger than ncv");

159:   /* call specific solver setup */
160:   (*svd->ops->setup)(svd);

162:   /* set tolerance if not yet set */
163:   if (svd->tol==PETSC_DEFAULT) svd->tol = SLEPC_DEFAULT_TOL;

165:   /* fill sorting criterion context */
166:   DSGetSlepcSC(svd->ds,&sc);
167:   sc->comparison    = (svd->which==SVD_LARGEST)? SlepcCompareLargestReal: SlepcCompareSmallestReal;
168:   sc->comparisonctx = NULL;
169:   sc->map           = NULL;
170:   sc->mapobj        = NULL;

172:   /* process initial vectors */
173:   if (svd->nini<0) {
174:     k = -svd->nini;
175:     if (k>svd->ncv) SETERRQ(PetscObjectComm((PetscObject)svd),1,"The number of initial vectors is larger than ncv");
176:     BVInsertVecs(svd->V,0,&k,svd->IS,PETSC_TRUE);
177:     SlepcBasisDestroy_Private(&svd->nini,&svd->IS);
178:     svd->nini = k;
179:   }
180:   if (svd->ninil<0) {
181:     k = 0;
182:     if (svd->leftbasis) {
183:       k = -svd->ninil;
184:       if (k>svd->ncv) SETERRQ(PetscObjectComm((PetscObject)svd),1,"The number of left initial vectors is larger than ncv");
185:       BVInsertVecs(svd->U,0,&k,svd->ISL,PETSC_TRUE);
186:     } else {
187:       PetscInfo(svd,"Ignoring initial left vectors\n");
188:     }
189:     SlepcBasisDestroy_Private(&svd->ninil,&svd->ISL);
190:     svd->ninil = k;
191:   }

193:   PetscLogEventEnd(SVD_SetUp,svd,0,0,0);
194:   svd->state = SVD_STATE_SETUP;
195:   return(0);
196: }

198: /*@C
199:    SVDSetInitialSpaces - Specify two basis of vectors that constitute the initial
200:    right and/or left spaces, that is, a rough approximation to the right and/or
201:    left singular subspaces from which the solver starts to iterate.

203:    Collective on SVD and Vec

205:    Input Parameter:
206: +  svd   - the singular value solver context
207: .  nr    - number of right vectors
208: .  isr   - set of basis vectors of the right initial space
209: .  nl    - number of left vectors
210: -  isl   - set of basis vectors of the left initial space

212:    Notes:
213:    It is not necessary to provide both sets of vectors.

215:    Some solvers start to iterate on a single vector (initial vector). In that case,
216:    the other vectors are ignored.

218:    These vectors do not persist from one SVDSolve() call to the other, so the
219:    initial space should be set every time.

221:    The vectors do not need to be mutually orthonormal, since they are explicitly
222:    orthonormalized internally.

224:    Common usage of this function is when the user can provide a rough approximation
225:    of the wanted singular space. Then, convergence may be faster.

227:    Level: intermediate
228: @*/
229: PetscErrorCode SVDSetInitialSpaces(SVD svd,PetscInt nr,Vec *isr,PetscInt nl,Vec *isl)
230: {

237:   if (nr<0) SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_OUTOFRANGE,"Argument nr cannot be negative");
238:   if (nl<0) SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_OUTOFRANGE,"Argument nl cannot be negative");
239:   SlepcBasisReference_Private(nr,isr,&svd->nini,&svd->IS);
240:   SlepcBasisReference_Private(nl,isl,&svd->ninil,&svd->ISL);
241:   if (nr>0 || nl>0) svd->state = SVD_STATE_INITIAL;
242:   return(0);
243: }

245: /*
246:   SVDSetDimensions_Default - Set reasonable values for ncv, mpd if not set
247:   by the user. This is called at setup.
248:  */
249: PetscErrorCode SVDSetDimensions_Default(SVD svd)
250: {
252:   PetscInt       N;

255:   SVDMatGetSize(svd,NULL,&N);
256:   if (svd->ncv) { /* ncv set */
257:     if (svd->ncv<svd->nsv) SETERRQ(PetscObjectComm((PetscObject)svd),1,"The value of ncv must be at least nsv");
258:   } else if (svd->mpd) { /* mpd set */
259:     svd->ncv = PetscMin(N,svd->nsv+svd->mpd);
260:   } else { /* neither set: defaults depend on nsv being small or large */
261:     if (svd->nsv<500) svd->ncv = PetscMin(N,PetscMax(2*svd->nsv,10));
262:     else {
263:       svd->mpd = 500;
264:       svd->ncv = PetscMin(N,svd->nsv+svd->mpd);
265:     }
266:   }
267:   if (!svd->mpd) svd->mpd = svd->ncv;
268:   return(0);
269: }

271: /*@
272:    SVDAllocateSolution - Allocate memory storage for common variables such
273:    as the singular values and the basis vectors.

275:    Collective on SVD

277:    Input Parameters:
278: +  svd   - eigensolver context
279: -  extra - number of additional positions, used for methods that require a
280:            working basis slightly larger than ncv

282:    Developers Notes:
283:    This is SLEPC_EXTERN because it may be required by user plugin SVD
284:    implementations.

286:    This is called at setup after setting the value of ncv and the flag leftbasis.

288:    Level: developer
289: @*/
290: PetscErrorCode SVDAllocateSolution(SVD svd,PetscInt extra)
291: {
293:   PetscInt       oldsize,requested;
294:   Vec            tr,tl;

297:   requested = svd->ncv + extra;

299:   /* oldsize is zero if this is the first time setup is called */
300:   BVGetSizes(svd->V,NULL,NULL,&oldsize);

302:   /* allocate sigma */
303:   if (requested != oldsize || !svd->sigma) {
304:     PetscFree3(svd->sigma,svd->perm,svd->errest);
305:     PetscMalloc3(requested,&svd->sigma,requested,&svd->perm,requested,&svd->errest);
306:     PetscLogObjectMemory((PetscObject)svd,PetscMax(0,requested-oldsize)*(2*sizeof(PetscReal)+sizeof(PetscInt)));
307:   }
308:   /* allocate V */
309:   if (!svd->V) { SVDGetBV(svd,&svd->V,NULL); }
310:   if (!oldsize) {
311:     if (!((PetscObject)(svd->V))->type_name) {
312:       BVSetType(svd->V,BVSVEC);
313:     }
314:     SVDMatCreateVecsEmpty(svd,&tr,NULL);
315:     BVSetSizesFromVec(svd->V,tr,requested);
316:     VecDestroy(&tr);
317:   } else {
318:     BVResize(svd->V,requested,PETSC_FALSE);
319:   }
320:   /* allocate U */
321:   if (svd->leftbasis) {
322:     if (!svd->U) { SVDGetBV(svd,NULL,&svd->U); }
323:     if (!oldsize) {
324:       if (!((PetscObject)(svd->U))->type_name) {
325:         BVSetType(svd->U,BVSVEC);
326:       }
327:       SVDMatCreateVecsEmpty(svd,NULL,&tl);
328:       BVSetSizesFromVec(svd->U,tl,requested);
329:       VecDestroy(&tl);
330:     } else {
331:       BVResize(svd->U,requested,PETSC_FALSE);
332:     }
333:   }
334:   return(0);
335: }