Actual source code: svdbasic.c
slepc-3.17.2 2022-08-09
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 = 0;
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;
53: *outsvd = 0;
54: SVDInitializePackage();
55: SlepcHeaderCreate(svd,SVD_CLASSID,"SVD","Singular Value Decomposition","SVD",comm,SVDDestroy,SVDView);
57: svd->OP = NULL;
58: svd->OPb = NULL;
59: svd->max_it = PETSC_DEFAULT;
60: svd->nsv = 1;
61: svd->ncv = PETSC_DEFAULT;
62: svd->mpd = PETSC_DEFAULT;
63: svd->nini = 0;
64: svd->ninil = 0;
65: svd->tol = PETSC_DEFAULT;
66: svd->conv = (SVDConv)-1;
67: svd->stop = SVD_STOP_BASIC;
68: svd->which = SVD_LARGEST;
69: svd->problem_type = (SVDProblemType)0;
70: svd->impltrans = PETSC_FALSE;
71: svd->trackall = PETSC_FALSE;
73: svd->converged = NULL;
74: svd->convergeduser = NULL;
75: svd->convergeddestroy = NULL;
76: svd->stopping = SVDStoppingBasic;
77: svd->stoppinguser = NULL;
78: svd->stoppingdestroy = NULL;
79: svd->convergedctx = NULL;
80: svd->stoppingctx = NULL;
81: svd->numbermonitors = 0;
83: svd->ds = NULL;
84: svd->U = NULL;
85: svd->V = NULL;
86: svd->A = NULL;
87: svd->B = NULL;
88: svd->AT = NULL;
89: svd->BT = NULL;
90: svd->IS = NULL;
91: svd->ISL = NULL;
92: svd->sigma = NULL;
93: svd->errest = NULL;
94: svd->perm = NULL;
95: svd->nworkl = 0;
96: svd->nworkr = 0;
97: svd->workl = NULL;
98: svd->workr = NULL;
99: svd->data = NULL;
101: svd->state = SVD_STATE_INITIAL;
102: svd->nconv = 0;
103: svd->its = 0;
104: svd->leftbasis = PETSC_FALSE;
105: svd->swapped = PETSC_FALSE;
106: svd->expltrans = PETSC_FALSE;
107: svd->nrma = 0.0;
108: svd->nrmb = 0.0;
109: svd->isgeneralized = PETSC_FALSE;
110: svd->reason = SVD_CONVERGED_ITERATING;
112: PetscNewLog(svd,&svd->sc);
113: *outsvd = svd;
114: PetscFunctionReturn(0);
115: }
117: /*@
118: SVDReset - Resets the SVD context to the initial state (prior to setup)
119: and destroys any allocated Vecs and Mats.
121: Collective on svd
123: Input Parameter:
124: . svd - singular value solver context obtained from SVDCreate()
126: Level: advanced
128: .seealso: SVDDestroy()
129: @*/
130: PetscErrorCode SVDReset(SVD svd)
131: {
133: if (!svd) PetscFunctionReturn(0);
134: if (svd->ops->reset) (svd->ops->reset)(svd);
135: MatDestroy(&svd->OP);
136: MatDestroy(&svd->OPb);
137: MatDestroy(&svd->A);
138: MatDestroy(&svd->B);
139: MatDestroy(&svd->AT);
140: MatDestroy(&svd->BT);
141: BVDestroy(&svd->U);
142: BVDestroy(&svd->V);
143: VecDestroyVecs(svd->nworkl,&svd->workl);
144: svd->nworkl = 0;
145: VecDestroyVecs(svd->nworkr,&svd->workr);
146: svd->nworkr = 0;
147: svd->state = SVD_STATE_INITIAL;
148: PetscFunctionReturn(0);
149: }
151: /*@C
152: SVDDestroy - Destroys the SVD context.
154: Collective on svd
156: Input Parameter:
157: . svd - singular value solver context obtained from SVDCreate()
159: Level: beginner
161: .seealso: SVDCreate(), SVDSetUp(), SVDSolve()
162: @*/
163: PetscErrorCode SVDDestroy(SVD *svd)
164: {
165: if (!*svd) PetscFunctionReturn(0);
167: if (--((PetscObject)(*svd))->refct > 0) { *svd = 0; PetscFunctionReturn(0); }
168: SVDReset(*svd);
169: if ((*svd)->ops->destroy) (*(*svd)->ops->destroy)(*svd);
170: if ((*svd)->sigma) PetscFree3((*svd)->sigma,(*svd)->perm,(*svd)->errest);
171: DSDestroy(&(*svd)->ds);
172: PetscFree((*svd)->sc);
173: /* just in case the initial vectors have not been used */
174: SlepcBasisDestroy_Private(&(*svd)->nini,&(*svd)->IS);
175: SlepcBasisDestroy_Private(&(*svd)->ninil,&(*svd)->ISL);
176: SVDMonitorCancel(*svd);
177: PetscHeaderDestroy(svd);
178: PetscFunctionReturn(0);
179: }
181: /*@C
182: SVDSetType - Selects the particular solver to be used in the SVD object.
184: Logically Collective on svd
186: Input Parameters:
187: + svd - the singular value solver context
188: - type - a known method
190: Options Database Key:
191: . -svd_type <method> - Sets the method; use -help for a list
192: of available methods
194: Notes:
195: See "slepc/include/slepcsvd.h" for available methods. The default
196: is SVDCROSS.
198: Normally, it is best to use the SVDSetFromOptions() command and
199: then set the SVD type from the options database rather than by using
200: this routine. Using the options database provides the user with
201: maximum flexibility in evaluating the different available methods.
202: The SVDSetType() routine is provided for those situations where it
203: is necessary to set the iterative solver independently of the command
204: line or options database.
206: Level: intermediate
208: .seealso: SVDType
209: @*/
210: PetscErrorCode SVDSetType(SVD svd,SVDType type)
211: {
212: PetscErrorCode (*r)(SVD);
213: PetscBool match;
218: PetscObjectTypeCompare((PetscObject)svd,type,&match);
219: if (match) PetscFunctionReturn(0);
221: PetscFunctionListFind(SVDList,type,&r);
224: if (svd->ops->destroy) (*svd->ops->destroy)(svd);
225: PetscMemzero(svd->ops,sizeof(struct _SVDOps));
227: svd->state = SVD_STATE_INITIAL;
228: PetscObjectChangeTypeName((PetscObject)svd,type);
229: (*r)(svd);
230: PetscFunctionReturn(0);
231: }
233: /*@C
234: SVDGetType - Gets the SVD type as a string from the SVD object.
236: Not Collective
238: Input Parameter:
239: . svd - the singular value solver context
241: Output Parameter:
242: . type - name of SVD method
244: Level: intermediate
246: .seealso: SVDSetType()
247: @*/
248: PetscErrorCode SVDGetType(SVD svd,SVDType *type)
249: {
252: *type = ((PetscObject)svd)->type_name;
253: PetscFunctionReturn(0);
254: }
256: /*@C
257: SVDRegister - Adds a method to the singular value solver package.
259: Not Collective
261: Input Parameters:
262: + name - name of a new user-defined solver
263: - function - routine to create the solver context
265: Notes:
266: SVDRegister() may be called multiple times to add several user-defined solvers.
268: Sample usage:
269: .vb
270: SVDRegister("my_solver",MySolverCreate);
271: .ve
273: Then, your solver can be chosen with the procedural interface via
274: $ SVDSetType(svd,"my_solver")
275: or at runtime via the option
276: $ -svd_type my_solver
278: Level: advanced
280: .seealso: SVDRegisterAll()
281: @*/
282: PetscErrorCode SVDRegister(const char *name,PetscErrorCode (*function)(SVD))
283: {
284: SVDInitializePackage();
285: PetscFunctionListAdd(&SVDList,name,function);
286: PetscFunctionReturn(0);
287: }
289: /*@C
290: SVDMonitorRegister - Adds SVD monitor routine.
292: Not Collective
294: Input Parameters:
295: + name - name of a new monitor routine
296: . vtype - a PetscViewerType for the output
297: . format - a PetscViewerFormat for the output
298: . monitor - monitor routine
299: . create - creation routine, or NULL
300: - destroy - destruction routine, or NULL
302: Notes:
303: SVDMonitorRegister() may be called multiple times to add several user-defined monitors.
305: Sample usage:
306: .vb
307: SVDMonitorRegister("my_monitor",PETSCVIEWERASCII,PETSC_VIEWER_ASCII_INFO_DETAIL,MyMonitor,NULL,NULL);
308: .ve
310: Then, your monitor can be chosen with the procedural interface via
311: $ SVDMonitorSetFromOptions(svd,"-svd_monitor_my_monitor","my_monitor",NULL)
312: or at runtime via the option
313: $ -svd_monitor_my_monitor
315: Level: advanced
317: .seealso: SVDMonitorRegisterAll()
318: @*/
319: 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**))
320: {
321: char key[PETSC_MAX_PATH_LEN];
323: SVDInitializePackage();
324: SlepcMonitorMakeKey_Internal(name,vtype,format,key);
325: PetscFunctionListAdd(&SVDMonitorList,key,monitor);
326: if (create) PetscFunctionListAdd(&SVDMonitorCreateList,key,create);
327: if (destroy) PetscFunctionListAdd(&SVDMonitorDestroyList,key,destroy);
328: PetscFunctionReturn(0);
329: }
331: /*@
332: SVDSetBV - Associates basis vectors objects to the singular value solver.
334: Collective on svd
336: Input Parameters:
337: + svd - singular value solver context obtained from SVDCreate()
338: . V - the basis vectors object for right singular vectors
339: - U - the basis vectors object for left singular vectors
341: Note:
342: Use SVDGetBV() to retrieve the basis vectors contexts (for example,
343: to free them at the end of the computations).
345: Level: advanced
347: .seealso: SVDGetBV()
348: @*/
349: PetscErrorCode SVDSetBV(SVD svd,BV V,BV U)
350: {
352: if (V) {
355: PetscObjectReference((PetscObject)V);
356: BVDestroy(&svd->V);
357: svd->V = V;
358: PetscLogObjectParent((PetscObject)svd,(PetscObject)svd->V);
359: }
360: if (U) {
363: PetscObjectReference((PetscObject)U);
364: BVDestroy(&svd->U);
365: svd->U = U;
366: PetscLogObjectParent((PetscObject)svd,(PetscObject)svd->U);
367: }
368: PetscFunctionReturn(0);
369: }
371: /*@
372: SVDGetBV - Obtain the basis vectors objects associated to the singular
373: value solver object.
375: Not Collective
377: Input Parameter:
378: . svd - singular value solver context obtained from SVDCreate()
380: Output Parameters:
381: + V - basis vectors context for right singular vectors
382: - U - basis vectors context for left singular vectors
384: Level: advanced
386: .seealso: SVDSetBV()
387: @*/
388: PetscErrorCode SVDGetBV(SVD svd,BV *V,BV *U)
389: {
391: if (V) {
392: if (!svd->V) {
393: BVCreate(PetscObjectComm((PetscObject)svd),&svd->V);
394: PetscObjectIncrementTabLevel((PetscObject)svd->V,(PetscObject)svd,0);
395: PetscLogObjectParent((PetscObject)svd,(PetscObject)svd->V);
396: PetscObjectSetOptions((PetscObject)svd->V,((PetscObject)svd)->options);
397: }
398: *V = svd->V;
399: }
400: if (U) {
401: if (!svd->U) {
402: BVCreate(PetscObjectComm((PetscObject)svd),&svd->U);
403: PetscObjectIncrementTabLevel((PetscObject)svd->U,(PetscObject)svd,0);
404: PetscLogObjectParent((PetscObject)svd,(PetscObject)svd->U);
405: PetscObjectSetOptions((PetscObject)svd->U,((PetscObject)svd)->options);
406: }
407: *U = svd->U;
408: }
409: PetscFunctionReturn(0);
410: }
412: /*@
413: SVDSetDS - Associates a direct solver object to the singular value solver.
415: Collective on svd
417: Input Parameters:
418: + svd - singular value solver context obtained from SVDCreate()
419: - ds - the direct solver object
421: Note:
422: Use SVDGetDS() to retrieve the direct solver context (for example,
423: to free it at the end of the computations).
425: Level: advanced
427: .seealso: SVDGetDS()
428: @*/
429: PetscErrorCode SVDSetDS(SVD svd,DS ds)
430: {
434: PetscObjectReference((PetscObject)ds);
435: DSDestroy(&svd->ds);
436: svd->ds = ds;
437: PetscLogObjectParent((PetscObject)svd,(PetscObject)svd->ds);
438: PetscFunctionReturn(0);
439: }
441: /*@
442: SVDGetDS - Obtain the direct solver object associated to the singular value
443: solver object.
445: Not Collective
447: Input Parameters:
448: . svd - singular value solver context obtained from SVDCreate()
450: Output Parameter:
451: . ds - direct solver context
453: Level: advanced
455: .seealso: SVDSetDS()
456: @*/
457: PetscErrorCode SVDGetDS(SVD svd,DS *ds)
458: {
461: if (!svd->ds) {
462: DSCreate(PetscObjectComm((PetscObject)svd),&svd->ds);
463: PetscObjectIncrementTabLevel((PetscObject)svd->ds,(PetscObject)svd,0);
464: PetscLogObjectParent((PetscObject)svd,(PetscObject)svd->ds);
465: PetscObjectSetOptions((PetscObject)svd->ds,((PetscObject)svd)->options);
466: }
467: *ds = svd->ds;
468: PetscFunctionReturn(0);
469: }