Actual source code: svdbasic.c

slepc-3.22.1 2024-10-28
Report Typos and Errors
  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:    Basic SVD routines
 12: */

 14: #include <slepc/private/svdimpl.h>

 16: /* Logging support */
 17: PetscClassId      SVD_CLASSID = 0;
 18: PetscLogEvent     SVD_SetUp = 0,SVD_Solve = 0;

 20: /* List of registered SVD routines */
 21: PetscFunctionList SVDList = NULL;
 22: PetscBool         SVDRegisterAllCalled = PETSC_FALSE;

 24: /* List of registered SVD monitors */
 25: PetscFunctionList SVDMonitorList              = NULL;
 26: PetscFunctionList SVDMonitorCreateList        = NULL;
 27: PetscFunctionList SVDMonitorDestroyList       = NULL;
 28: PetscBool         SVDMonitorRegisterAllCalled = PETSC_FALSE;

 30: /*@
 31:    SVDCreate - Creates the default SVD context.

 33:    Collective

 35:    Input Parameter:
 36: .  comm - MPI communicator

 38:    Output Parameter:
 39: .  outsvd - location to put the SVD context

 41:    Note:
 42:    The default SVD type is SVDCROSS

 44:    Level: beginner

 46: .seealso: SVDSetUp(), SVDSolve(), SVDDestroy(), SVD
 47: @*/
 48: PetscErrorCode SVDCreate(MPI_Comm comm,SVD *outsvd)
 49: {
 50:   SVD            svd;

 52:   PetscFunctionBegin;
 53:   PetscAssertPointer(outsvd,2);
 54:   PetscCall(SVDInitializePackage());
 55:   PetscCall(SlepcHeaderCreate(svd,SVD_CLASSID,"SVD","Singular Value Decomposition","SVD",comm,SVDDestroy,SVDView));

 57:   svd->OP               = NULL;
 58:   svd->OPb              = NULL;
 59:   svd->omega            = NULL;
 60:   svd->max_it           = PETSC_DETERMINE;
 61:   svd->nsv              = 1;
 62:   svd->ncv              = PETSC_DETERMINE;
 63:   svd->mpd              = PETSC_DETERMINE;
 64:   svd->nini             = 0;
 65:   svd->ninil            = 0;
 66:   svd->tol              = PETSC_DETERMINE;
 67:   svd->conv             = (SVDConv)-1;
 68:   svd->stop             = SVD_STOP_BASIC;
 69:   svd->which            = SVD_LARGEST;
 70:   svd->problem_type     = (SVDProblemType)0;
 71:   svd->impltrans        = PETSC_FALSE;
 72:   svd->trackall         = PETSC_FALSE;

 74:   svd->converged        = NULL;
 75:   svd->convergeduser    = NULL;
 76:   svd->convergeddestroy = NULL;
 77:   svd->stopping         = SVDStoppingBasic;
 78:   svd->stoppinguser     = NULL;
 79:   svd->stoppingdestroy  = NULL;
 80:   svd->convergedctx     = NULL;
 81:   svd->stoppingctx      = NULL;
 82:   svd->numbermonitors   = 0;

 84:   svd->ds               = NULL;
 85:   svd->U                = NULL;
 86:   svd->V                = NULL;
 87:   svd->A                = NULL;
 88:   svd->B                = NULL;
 89:   svd->AT               = NULL;
 90:   svd->BT               = NULL;
 91:   svd->IS               = NULL;
 92:   svd->ISL              = NULL;
 93:   svd->sigma            = NULL;
 94:   svd->errest           = NULL;
 95:   svd->sign             = NULL;
 96:   svd->perm             = NULL;
 97:   svd->nworkl           = 0;
 98:   svd->nworkr           = 0;
 99:   svd->workl            = NULL;
100:   svd->workr            = NULL;
101:   svd->data             = NULL;

103:   svd->state            = SVD_STATE_INITIAL;
104:   svd->nconv            = 0;
105:   svd->its              = 0;
106:   svd->leftbasis        = PETSC_FALSE;
107:   svd->swapped          = PETSC_FALSE;
108:   svd->expltrans        = PETSC_FALSE;
109:   svd->nrma             = 0.0;
110:   svd->nrmb             = 0.0;
111:   svd->isgeneralized    = PETSC_FALSE;
112:   svd->reason           = SVD_CONVERGED_ITERATING;

114:   PetscCall(PetscNew(&svd->sc));
115:   *outsvd = svd;
116:   PetscFunctionReturn(PETSC_SUCCESS);
117: }

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

123:    Collective

125:    Input Parameter:
126: .  svd - singular value solver context obtained from SVDCreate()

128:    Level: advanced

130: .seealso: SVDDestroy()
131: @*/
132: PetscErrorCode SVDReset(SVD svd)
133: {
134:   PetscFunctionBegin;
136:   if (!svd) PetscFunctionReturn(PETSC_SUCCESS);
137:   PetscTryTypeMethod(svd,reset);
138:   PetscCall(MatDestroy(&svd->OP));
139:   PetscCall(MatDestroy(&svd->OPb));
140:   PetscCall(VecDestroy(&svd->omega));
141:   PetscCall(MatDestroy(&svd->A));
142:   PetscCall(MatDestroy(&svd->B));
143:   PetscCall(MatDestroy(&svd->AT));
144:   PetscCall(MatDestroy(&svd->BT));
145:   PetscCall(BVDestroy(&svd->U));
146:   PetscCall(BVDestroy(&svd->V));
147:   PetscCall(VecDestroyVecs(svd->nworkl,&svd->workl));
148:   svd->nworkl = 0;
149:   PetscCall(VecDestroyVecs(svd->nworkr,&svd->workr));
150:   svd->nworkr = 0;
151:   svd->swapped = PETSC_FALSE;
152:   svd->state = SVD_STATE_INITIAL;
153:   PetscFunctionReturn(PETSC_SUCCESS);
154: }

156: /*@
157:    SVDDestroy - Destroys the SVD context.

159:    Collective

161:    Input Parameter:
162: .  svd - singular value solver context obtained from SVDCreate()

164:    Level: beginner

166: .seealso: SVDCreate(), SVDSetUp(), SVDSolve()
167: @*/
168: PetscErrorCode SVDDestroy(SVD *svd)
169: {
170:   PetscFunctionBegin;
171:   if (!*svd) PetscFunctionReturn(PETSC_SUCCESS);
173:   if (--((PetscObject)*svd)->refct > 0) { *svd = NULL; PetscFunctionReturn(PETSC_SUCCESS); }
174:   PetscCall(SVDReset(*svd));
175:   PetscTryTypeMethod(*svd,destroy);
176:   if ((*svd)->sigma) PetscCall(PetscFree3((*svd)->sigma,(*svd)->perm,(*svd)->errest));
177:   if ((*svd)->sign) PetscCall(PetscFree((*svd)->sign));
178:   PetscCall(DSDestroy(&(*svd)->ds));
179:   PetscCall(PetscFree((*svd)->sc));
180:   /* just in case the initial vectors have not been used */
181:   PetscCall(SlepcBasisDestroy_Private(&(*svd)->nini,&(*svd)->IS));
182:   PetscCall(SlepcBasisDestroy_Private(&(*svd)->ninil,&(*svd)->ISL));
183:   PetscCall(SVDMonitorCancel(*svd));
184:   PetscCall(PetscHeaderDestroy(svd));
185:   PetscFunctionReturn(PETSC_SUCCESS);
186: }

188: /*@
189:    SVDSetType - Selects the particular solver to be used in the SVD object.

191:    Logically Collective

193:    Input Parameters:
194: +  svd      - the singular value solver context
195: -  type     - a known method

197:    Options Database Key:
198: .  -svd_type <method> - Sets the method; use -help for a list
199:     of available methods

201:    Notes:
202:    See "slepc/include/slepcsvd.h" for available methods. The default
203:    is SVDCROSS.

205:    Normally, it is best to use the SVDSetFromOptions() command and
206:    then set the SVD type from the options database rather than by using
207:    this routine.  Using the options database provides the user with
208:    maximum flexibility in evaluating the different available methods.
209:    The SVDSetType() routine is provided for those situations where it
210:    is necessary to set the iterative solver independently of the command
211:    line or options database.

213:    Level: intermediate

215: .seealso: SVDType
216: @*/
217: PetscErrorCode SVDSetType(SVD svd,SVDType type)
218: {
219:   PetscErrorCode (*r)(SVD);
220:   PetscBool      match;

222:   PetscFunctionBegin;
224:   PetscAssertPointer(type,2);

226:   PetscCall(PetscObjectTypeCompare((PetscObject)svd,type,&match));
227:   if (match) PetscFunctionReturn(PETSC_SUCCESS);

229:   PetscCall(PetscFunctionListFind(SVDList,type,&r));
230:   PetscCheck(r,PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown SVD type given: %s",type);

232:   PetscTryTypeMethod(svd,destroy);
233:   PetscCall(PetscMemzero(svd->ops,sizeof(struct _SVDOps)));

235:   svd->state = SVD_STATE_INITIAL;
236:   PetscCall(PetscObjectChangeTypeName((PetscObject)svd,type));
237:   PetscCall((*r)(svd));
238:   PetscFunctionReturn(PETSC_SUCCESS);
239: }

241: /*@
242:    SVDGetType - Gets the SVD type as a string from the SVD object.

244:    Not Collective

246:    Input Parameter:
247: .  svd - the singular value solver context

249:    Output Parameter:
250: .  type - name of SVD method

252:    Level: intermediate

254: .seealso: SVDSetType()
255: @*/
256: PetscErrorCode SVDGetType(SVD svd,SVDType *type)
257: {
258:   PetscFunctionBegin;
260:   PetscAssertPointer(type,2);
261:   *type = ((PetscObject)svd)->type_name;
262:   PetscFunctionReturn(PETSC_SUCCESS);
263: }

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

268:    Not Collective

270:    Input Parameters:
271: +  name - name of a new user-defined solver
272: -  function - routine to create the solver context

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

277:    Example Usage:
278: .vb
279:     SVDRegister("my_solver",MySolverCreate);
280: .ve

282:    Then, your solver can be chosen with the procedural interface via
283: $     SVDSetType(svd,"my_solver")
284:    or at runtime via the option
285: $     -svd_type my_solver

287:    Level: advanced

289: .seealso: SVDRegisterAll()
290: @*/
291: PetscErrorCode SVDRegister(const char *name,PetscErrorCode (*function)(SVD))
292: {
293:   PetscFunctionBegin;
294:   PetscCall(SVDInitializePackage());
295:   PetscCall(PetscFunctionListAdd(&SVDList,name,function));
296:   PetscFunctionReturn(PETSC_SUCCESS);
297: }

299: /*@C
300:    SVDMonitorRegister - Adds SVD monitor routine.

302:    Not Collective

304:    Input Parameters:
305: +  name    - name of a new monitor routine
306: .  vtype   - a PetscViewerType for the output
307: .  format  - a PetscViewerFormat for the output
308: .  monitor - monitor routine
309: .  create  - creation routine, or NULL
310: -  destroy - destruction routine, or NULL

312:    Notes:
313:    SVDMonitorRegister() may be called multiple times to add several user-defined monitors.

315:    Example Usage:
316: .vb
317:    SVDMonitorRegister("my_monitor",PETSCVIEWERASCII,PETSC_VIEWER_ASCII_INFO_DETAIL,MyMonitor,NULL,NULL);
318: .ve

320:    Then, your monitor can be chosen with the procedural interface via
321: $      SVDMonitorSetFromOptions(svd,"-svd_monitor_my_monitor","my_monitor",NULL)
322:    or at runtime via the option
323: $      -svd_monitor_my_monitor

325:    Level: advanced

327: .seealso: SVDMonitorRegisterAll()
328: @*/
329: PetscErrorCode SVDMonitorRegister(const char name[],PetscViewerType vtype,PetscViewerFormat format,PetscErrorCode (*monitor)(SVD,PetscInt,PetscInt,PetscReal*,PetscReal*,PetscInt,PetscViewerAndFormat*),PetscErrorCode (*create)(PetscViewer,PetscViewerFormat,void*,PetscViewerAndFormat**),PetscErrorCode (*destroy)(PetscViewerAndFormat**))
330: {
331:   char           key[PETSC_MAX_PATH_LEN];

333:   PetscFunctionBegin;
334:   PetscCall(SVDInitializePackage());
335:   PetscCall(SlepcMonitorMakeKey_Internal(name,vtype,format,key));
336:   PetscCall(PetscFunctionListAdd(&SVDMonitorList,key,monitor));
337:   if (create)  PetscCall(PetscFunctionListAdd(&SVDMonitorCreateList,key,create));
338:   if (destroy) PetscCall(PetscFunctionListAdd(&SVDMonitorDestroyList,key,destroy));
339:   PetscFunctionReturn(PETSC_SUCCESS);
340: }

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

345:    Collective

347:    Input Parameters:
348: +  svd - singular value solver context obtained from SVDCreate()
349: .  V   - the basis vectors object for right singular vectors
350: -  U   - the basis vectors object for left singular vectors

352:    Note:
353:    Use SVDGetBV() to retrieve the basis vectors contexts (for example,
354:    to free them at the end of the computations).

356:    Level: advanced

358: .seealso: SVDGetBV()
359: @*/
360: PetscErrorCode SVDSetBV(SVD svd,BV V,BV U)
361: {
362:   PetscFunctionBegin;
364:   if (V) {
366:     PetscCheckSameComm(svd,1,V,2);
367:     PetscCall(PetscObjectReference((PetscObject)V));
368:     PetscCall(BVDestroy(&svd->V));
369:     svd->V = V;
370:   }
371:   if (U) {
373:     PetscCheckSameComm(svd,1,U,3);
374:     PetscCall(PetscObjectReference((PetscObject)U));
375:     PetscCall(BVDestroy(&svd->U));
376:     svd->U = U;
377:   }
378:   PetscFunctionReturn(PETSC_SUCCESS);
379: }

381: /*@
382:    SVDGetBV - Obtain the basis vectors objects associated to the singular
383:    value solver object.

385:    Not Collective

387:    Input Parameter:
388: .  svd - singular value solver context obtained from SVDCreate()

390:    Output Parameters:
391: +  V - basis vectors context for right singular vectors
392: -  U - basis vectors context for left singular vectors

394:    Level: advanced

396: .seealso: SVDSetBV()
397: @*/
398: PetscErrorCode SVDGetBV(SVD svd,BV *V,BV *U)
399: {
400:   PetscFunctionBegin;
402:   if (V) {
403:     if (!svd->V) {
404:       PetscCall(BVCreate(PetscObjectComm((PetscObject)svd),&svd->V));
405:       PetscCall(PetscObjectIncrementTabLevel((PetscObject)svd->V,(PetscObject)svd,0));
406:       PetscCall(PetscObjectSetOptions((PetscObject)svd->V,((PetscObject)svd)->options));
407:     }
408:     *V = svd->V;
409:   }
410:   if (U) {
411:     if (!svd->U) {
412:       PetscCall(BVCreate(PetscObjectComm((PetscObject)svd),&svd->U));
413:       PetscCall(PetscObjectIncrementTabLevel((PetscObject)svd->U,(PetscObject)svd,0));
414:       PetscCall(PetscObjectSetOptions((PetscObject)svd->U,((PetscObject)svd)->options));
415:     }
416:     *U = svd->U;
417:   }
418:   PetscFunctionReturn(PETSC_SUCCESS);
419: }

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

424:    Collective

426:    Input Parameters:
427: +  svd - singular value solver context obtained from SVDCreate()
428: -  ds  - the direct solver object

430:    Note:
431:    Use SVDGetDS() to retrieve the direct solver context (for example,
432:    to free it at the end of the computations).

434:    Level: advanced

436: .seealso: SVDGetDS()
437: @*/
438: PetscErrorCode SVDSetDS(SVD svd,DS ds)
439: {
440:   PetscFunctionBegin;
443:   PetscCheckSameComm(svd,1,ds,2);
444:   PetscCall(PetscObjectReference((PetscObject)ds));
445:   PetscCall(DSDestroy(&svd->ds));
446:   svd->ds = ds;
447:   PetscFunctionReturn(PETSC_SUCCESS);
448: }

450: /*@
451:    SVDGetDS - Obtain the direct solver object associated to the singular value
452:    solver object.

454:    Not Collective

456:    Input Parameters:
457: .  svd - singular value solver context obtained from SVDCreate()

459:    Output Parameter:
460: .  ds - direct solver context

462:    Level: advanced

464: .seealso: SVDSetDS()
465: @*/
466: PetscErrorCode SVDGetDS(SVD svd,DS *ds)
467: {
468:   PetscFunctionBegin;
470:   PetscAssertPointer(ds,2);
471:   if (!svd->ds) {
472:     PetscCall(DSCreate(PetscObjectComm((PetscObject)svd),&svd->ds));
473:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)svd->ds,(PetscObject)svd,0));
474:     PetscCall(PetscObjectSetOptions((PetscObject)svd->ds,((PetscObject)svd)->options));
475:   }
476:   *ds = svd->ds;
477:   PetscFunctionReturn(PETSC_SUCCESS);
478: }