Actual source code: svdbasic.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:    Basic SVD routines
 12: */

 14: #include <slepc/private/svdimpl.h>      /*I "slepcsvd.h" I*/

 16: PetscFunctionList SVDList = 0;
 17: PetscBool         SVDRegisterAllCalled = PETSC_FALSE;
 18: PetscClassId      SVD_CLASSID = 0;
 19: PetscLogEvent     SVD_SetUp = 0,SVD_Solve = 0;

 21: /*@
 22:    SVDCreate - Creates the default SVD context.

 24:    Collective on MPI_Comm

 26:    Input Parameter:
 27: .  comm - MPI communicator

 29:    Output Parameter:
 30: .  svd - location to put the SVD context

 32:    Note:
 33:    The default SVD type is SVDCROSS

 35:    Level: beginner

 37: .seealso: SVDSetUp(), SVDSolve(), SVDDestroy(), SVD
 38: @*/
 39: PetscErrorCode SVDCreate(MPI_Comm comm,SVD *outsvd)
 40: {
 42:   SVD            svd;

 46:   *outsvd = 0;
 47:   SVDInitializePackage();
 48:   SlepcHeaderCreate(svd,SVD_CLASSID,"SVD","Singular Value Decomposition","SVD",comm,SVDDestroy,SVDView);

 50:   svd->OP               = NULL;
 51:   svd->max_it           = 0;
 52:   svd->nsv              = 1;
 53:   svd->ncv              = 0;
 54:   svd->mpd              = 0;
 55:   svd->nini             = 0;
 56:   svd->ninil            = 0;
 57:   svd->tol              = PETSC_DEFAULT;
 58:   svd->conv             = SVD_CONV_REL;
 59:   svd->stop             = SVD_STOP_BASIC;
 60:   svd->which            = SVD_LARGEST;
 61:   svd->impltrans        = PETSC_FALSE;
 62:   svd->trackall         = PETSC_FALSE;

 64:   svd->converged        = SVDConvergedRelative;
 65:   svd->convergeduser    = NULL;
 66:   svd->convergeddestroy = NULL;
 67:   svd->stopping         = SVDStoppingBasic;
 68:   svd->stoppinguser     = NULL;
 69:   svd->stoppingdestroy  = NULL;
 70:   svd->convergedctx     = NULL;
 71:   svd->stoppingctx      = NULL;
 72:   svd->numbermonitors   = 0;

 74:   svd->ds               = NULL;
 75:   svd->U                = NULL;
 76:   svd->V                = NULL;
 77:   svd->A                = NULL;
 78:   svd->AT               = NULL;
 79:   svd->IS               = NULL;
 80:   svd->ISL              = NULL;
 81:   svd->sigma            = NULL;
 82:   svd->perm             = NULL;
 83:   svd->errest           = NULL;
 84:   svd->data             = NULL;

 86:   svd->state            = SVD_STATE_INITIAL;
 87:   svd->nconv            = 0;
 88:   svd->its              = 0;
 89:   svd->leftbasis        = PETSC_FALSE;
 90:   svd->reason           = SVD_CONVERGED_ITERATING;

 92:   PetscNewLog(svd,&svd->sc);
 93:   *outsvd = svd;
 94:   return(0);
 95: }

 97: /*@
 98:    SVDReset - Resets the SVD context to the initial state (prior to setup)
 99:    and destroys any allocated Vecs and Mats.

101:    Collective on SVD

103:    Input Parameter:
104: .  svd - singular value solver context obtained from SVDCreate()

106:    Level: advanced

108: .seealso: SVDDestroy()
109: @*/
110: PetscErrorCode SVDReset(SVD svd)
111: {

116:   if (!svd) return(0);
117:   if (svd->ops->reset) { (svd->ops->reset)(svd); }
118:   MatDestroy(&svd->OP);
119:   MatDestroy(&svd->A);
120:   MatDestroy(&svd->AT);
121:   BVDestroy(&svd->U);
122:   BVDestroy(&svd->V);
123:   svd->state = SVD_STATE_INITIAL;
124:   return(0);
125: }

127: /*@
128:    SVDDestroy - Destroys the SVD context.

130:    Collective on SVD

132:    Input Parameter:
133: .  svd - singular value solver context obtained from SVDCreate()

135:    Level: beginner

137: .seealso: SVDCreate(), SVDSetUp(), SVDSolve()
138: @*/
139: PetscErrorCode SVDDestroy(SVD *svd)
140: {

144:   if (!*svd) return(0);
146:   if (--((PetscObject)(*svd))->refct > 0) { *svd = 0; return(0); }
147:   SVDReset(*svd);
148:   if ((*svd)->ops->destroy) { (*(*svd)->ops->destroy)(*svd); }
149:   if ((*svd)->sigma) {
150:     PetscFree3((*svd)->sigma,(*svd)->perm,(*svd)->errest);
151:   }
152:   DSDestroy(&(*svd)->ds);
153:   PetscFree((*svd)->sc);
154:   /* just in case the initial vectors have not been used */
155:   SlepcBasisDestroy_Private(&(*svd)->nini,&(*svd)->IS);
156:   SlepcBasisDestroy_Private(&(*svd)->ninil,&(*svd)->ISL);
157:   SVDMonitorCancel(*svd);
158:   PetscHeaderDestroy(svd);
159:   return(0);
160: }

162: /*@C
163:    SVDSetType - Selects the particular solver to be used in the SVD object.

165:    Logically Collective on SVD

167:    Input Parameters:
168: +  svd      - the singular value solver context
169: -  type     - a known method

171:    Options Database Key:
172: .  -svd_type <method> - Sets the method; use -help for a list
173:     of available methods

175:    Notes:
176:    See "slepc/include/slepcsvd.h" for available methods. The default
177:    is SVDCROSS.

179:    Normally, it is best to use the SVDSetFromOptions() command and
180:    then set the SVD type from the options database rather than by using
181:    this routine.  Using the options database provides the user with
182:    maximum flexibility in evaluating the different available methods.
183:    The SVDSetType() routine is provided for those situations where it
184:    is necessary to set the iterative solver independently of the command
185:    line or options database.

187:    Level: intermediate

189: .seealso: SVDType
190: @*/
191: PetscErrorCode SVDSetType(SVD svd,SVDType type)
192: {
193:   PetscErrorCode ierr,(*r)(SVD);
194:   PetscBool      match;


200:   PetscObjectTypeCompare((PetscObject)svd,type,&match);
201:   if (match) return(0);

203:   PetscFunctionListFind(SVDList,type,&r);
204:   if (!r) SETERRQ1(PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown SVD type given: %s",type);

206:   if (svd->ops->destroy) { (*svd->ops->destroy)(svd); }
207:   PetscMemzero(svd->ops,sizeof(struct _SVDOps));

209:   svd->state = SVD_STATE_INITIAL;
210:   PetscObjectChangeTypeName((PetscObject)svd,type);
211:   (*r)(svd);
212:   return(0);
213: }

215: /*@C
216:    SVDGetType - Gets the SVD type as a string from the SVD object.

218:    Not Collective

220:    Input Parameter:
221: .  svd - the singular value solver context

223:    Output Parameter:
224: .  name - name of SVD method

226:    Level: intermediate

228: .seealso: SVDSetType()
229: @*/
230: PetscErrorCode SVDGetType(SVD svd,SVDType *type)
231: {
235:   *type = ((PetscObject)svd)->type_name;
236:   return(0);
237: }

239: /*@C
240:    SVDRegister - Adds a method to the singular value solver package.

242:    Not Collective

244:    Input Parameters:
245: +  name - name of a new user-defined solver
246: -  function - routine to create the solver context

248:    Notes:
249:    SVDRegister() may be called multiple times to add several user-defined solvers.

251:    Sample usage:
252: .vb
253:     SVDRegister("my_solver",MySolverCreate);
254: .ve

256:    Then, your solver can be chosen with the procedural interface via
257: $     SVDSetType(svd,"my_solver")
258:    or at runtime via the option
259: $     -svd_type my_solver

261:    Level: advanced

263: .seealso: SVDRegisterAll()
264: @*/
265: PetscErrorCode SVDRegister(const char *name,PetscErrorCode (*function)(SVD))
266: {

270:   SVDInitializePackage();
271:   PetscFunctionListAdd(&SVDList,name,function);
272:   return(0);
273: }

275: /*@
276:    SVDSetBV - Associates basis vectors objects to the singular value solver.

278:    Collective on SVD

280:    Input Parameters:
281: +  svd - singular value solver context obtained from SVDCreate()
282: .  V   - the basis vectors object for right singular vectors
283: -  U   - the basis vectors object for left singular vectors

285:    Note:
286:    Use SVDGetBV() to retrieve the basis vectors contexts (for example,
287:    to free them at the end of the computations).

289:    Level: advanced

291: .seealso: SVDGetBV()
292: @*/
293: PetscErrorCode SVDSetBV(SVD svd,BV V,BV U)
294: {

299:   if (V) {
302:     PetscObjectReference((PetscObject)V);
303:     BVDestroy(&svd->V);
304:     svd->V = V;
305:     PetscLogObjectParent((PetscObject)svd,(PetscObject)svd->V);
306:   }
307:   if (U) {
310:     PetscObjectReference((PetscObject)U);
311:     BVDestroy(&svd->U);
312:     svd->U = U;
313:     PetscLogObjectParent((PetscObject)svd,(PetscObject)svd->U);
314:   }
315:   return(0);
316: }

318: /*@
319:    SVDGetBV - Obtain the basis vectors objects associated to the singular
320:    value solver object.

322:    Not Collective

324:    Input Parameters:
325: .  svd - singular value solver context obtained from SVDCreate()

327:    Output Parameter:
328: +  V - basis vectors context for right singular vectors
329: -  U - basis vectors context for left singular vectors

331:    Level: advanced

333: .seealso: SVDSetBV()
334: @*/
335: PetscErrorCode SVDGetBV(SVD svd,BV *V,BV *U)
336: {

341:   if (V) {
342:     if (!svd->V) {
343:       BVCreate(PetscObjectComm((PetscObject)svd),&svd->V);
344:       PetscObjectIncrementTabLevel((PetscObject)svd->V,(PetscObject)svd,0);
345:       PetscLogObjectParent((PetscObject)svd,(PetscObject)svd->V);
346:       PetscObjectSetOptions((PetscObject)svd->V,((PetscObject)svd)->options);
347:     }
348:     *V = svd->V;
349:   }
350:   if (U) {
351:     if (!svd->U) {
352:       BVCreate(PetscObjectComm((PetscObject)svd),&svd->U);
353:       PetscObjectIncrementTabLevel((PetscObject)svd->U,(PetscObject)svd,0);
354:       PetscLogObjectParent((PetscObject)svd,(PetscObject)svd->U);
355:       PetscObjectSetOptions((PetscObject)svd->U,((PetscObject)svd)->options);
356:     }
357:     *U = svd->U;
358:   }
359:   return(0);
360: }

362: /*@
363:    SVDSetDS - Associates a direct solver object to the singular value solver.

365:    Collective on SVD

367:    Input Parameters:
368: +  svd - singular value solver context obtained from SVDCreate()
369: -  ds  - the direct solver object

371:    Note:
372:    Use SVDGetDS() to retrieve the direct solver context (for example,
373:    to free it at the end of the computations).

375:    Level: advanced

377: .seealso: SVDGetDS()
378: @*/
379: PetscErrorCode SVDSetDS(SVD svd,DS ds)
380: {

387:   PetscObjectReference((PetscObject)ds);
388:   DSDestroy(&svd->ds);
389:   svd->ds = ds;
390:   PetscLogObjectParent((PetscObject)svd,(PetscObject)svd->ds);
391:   return(0);
392: }

394: /*@
395:    SVDGetDS - Obtain the direct solver object associated to the singular value
396:    solver object.

398:    Not Collective

400:    Input Parameters:
401: .  svd - singular value solver context obtained from SVDCreate()

403:    Output Parameter:
404: .  ds - direct solver context

406:    Level: advanced

408: .seealso: SVDSetDS()
409: @*/
410: PetscErrorCode SVDGetDS(SVD svd,DS *ds)
411: {

417:   if (!svd->ds) {
418:     DSCreate(PetscObjectComm((PetscObject)svd),&svd->ds);
419:     PetscObjectIncrementTabLevel((PetscObject)svd->ds,(PetscObject)svd,0);
420:     PetscLogObjectParent((PetscObject)svd,(PetscObject)svd->ds);
421:     PetscObjectSetOptions((PetscObject)svd->ds,((PetscObject)svd)->options);
422:   }
423:   *ds = svd->ds;
424:   return(0);
425: }