Actual source code: svdbasic.c
slepc-3.21.0 2024-03-30
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: *outsvd = NULL;
55: PetscCall(SVDInitializePackage());
56: PetscCall(SlepcHeaderCreate(svd,SVD_CLASSID,"SVD","Singular Value Decomposition","SVD",comm,SVDDestroy,SVDView));
58: svd->OP = NULL;
59: svd->OPb = NULL;
60: svd->omega = NULL;
61: svd->max_it = PETSC_DEFAULT;
62: svd->nsv = 1;
63: svd->ncv = PETSC_DEFAULT;
64: svd->mpd = PETSC_DEFAULT;
65: svd->nini = 0;
66: svd->ninil = 0;
67: svd->tol = PETSC_DEFAULT;
68: svd->conv = (SVDConv)-1;
69: svd->stop = SVD_STOP_BASIC;
70: svd->which = SVD_LARGEST;
71: svd->problem_type = (SVDProblemType)0;
72: svd->impltrans = PETSC_FALSE;
73: svd->trackall = PETSC_FALSE;
75: svd->converged = NULL;
76: svd->convergeduser = NULL;
77: svd->convergeddestroy = NULL;
78: svd->stopping = SVDStoppingBasic;
79: svd->stoppinguser = NULL;
80: svd->stoppingdestroy = NULL;
81: svd->convergedctx = NULL;
82: svd->stoppingctx = NULL;
83: svd->numbermonitors = 0;
85: svd->ds = NULL;
86: svd->U = NULL;
87: svd->V = NULL;
88: svd->A = NULL;
89: svd->B = NULL;
90: svd->AT = NULL;
91: svd->BT = NULL;
92: svd->IS = NULL;
93: svd->ISL = NULL;
94: svd->sigma = NULL;
95: svd->errest = NULL;
96: svd->sign = NULL;
97: svd->perm = NULL;
98: svd->nworkl = 0;
99: svd->nworkr = 0;
100: svd->workl = NULL;
101: svd->workr = NULL;
102: svd->data = NULL;
104: svd->state = SVD_STATE_INITIAL;
105: svd->nconv = 0;
106: svd->its = 0;
107: svd->leftbasis = PETSC_FALSE;
108: svd->swapped = PETSC_FALSE;
109: svd->expltrans = PETSC_FALSE;
110: svd->nrma = 0.0;
111: svd->nrmb = 0.0;
112: svd->isgeneralized = PETSC_FALSE;
113: svd->reason = SVD_CONVERGED_ITERATING;
115: PetscCall(PetscNew(&svd->sc));
116: *outsvd = svd;
117: PetscFunctionReturn(PETSC_SUCCESS);
118: }
120: /*@
121: SVDReset - Resets the SVD context to the initial state (prior to setup)
122: and destroys any allocated Vecs and Mats.
124: Collective
126: Input Parameter:
127: . svd - singular value solver context obtained from SVDCreate()
129: Level: advanced
131: .seealso: SVDDestroy()
132: @*/
133: PetscErrorCode SVDReset(SVD svd)
134: {
135: PetscFunctionBegin;
137: if (!svd) PetscFunctionReturn(PETSC_SUCCESS);
138: PetscTryTypeMethod(svd,reset);
139: PetscCall(MatDestroy(&svd->OP));
140: PetscCall(MatDestroy(&svd->OPb));
141: PetscCall(VecDestroy(&svd->omega));
142: PetscCall(MatDestroy(&svd->A));
143: PetscCall(MatDestroy(&svd->B));
144: PetscCall(MatDestroy(&svd->AT));
145: PetscCall(MatDestroy(&svd->BT));
146: PetscCall(BVDestroy(&svd->U));
147: PetscCall(BVDestroy(&svd->V));
148: PetscCall(VecDestroyVecs(svd->nworkl,&svd->workl));
149: svd->nworkl = 0;
150: PetscCall(VecDestroyVecs(svd->nworkr,&svd->workr));
151: svd->nworkr = 0;
152: svd->swapped = PETSC_FALSE;
153: svd->state = SVD_STATE_INITIAL;
154: PetscFunctionReturn(PETSC_SUCCESS);
155: }
157: /*@C
158: SVDDestroy - Destroys the SVD context.
160: Collective
162: Input Parameter:
163: . svd - singular value solver context obtained from SVDCreate()
165: Level: beginner
167: .seealso: SVDCreate(), SVDSetUp(), SVDSolve()
168: @*/
169: PetscErrorCode SVDDestroy(SVD *svd)
170: {
171: PetscFunctionBegin;
172: if (!*svd) PetscFunctionReturn(PETSC_SUCCESS);
174: if (--((PetscObject)*svd)->refct > 0) { *svd = NULL; PetscFunctionReturn(PETSC_SUCCESS); }
175: PetscCall(SVDReset(*svd));
176: PetscTryTypeMethod(*svd,destroy);
177: if ((*svd)->sigma) PetscCall(PetscFree3((*svd)->sigma,(*svd)->perm,(*svd)->errest));
178: if ((*svd)->sign) PetscCall(PetscFree((*svd)->sign));
179: PetscCall(DSDestroy(&(*svd)->ds));
180: PetscCall(PetscFree((*svd)->sc));
181: /* just in case the initial vectors have not been used */
182: PetscCall(SlepcBasisDestroy_Private(&(*svd)->nini,&(*svd)->IS));
183: PetscCall(SlepcBasisDestroy_Private(&(*svd)->ninil,&(*svd)->ISL));
184: PetscCall(SVDMonitorCancel(*svd));
185: PetscCall(PetscHeaderDestroy(svd));
186: PetscFunctionReturn(PETSC_SUCCESS);
187: }
189: /*@C
190: SVDSetType - Selects the particular solver to be used in the SVD object.
192: Logically Collective
194: Input Parameters:
195: + svd - the singular value solver context
196: - type - a known method
198: Options Database Key:
199: . -svd_type <method> - Sets the method; use -help for a list
200: of available methods
202: Notes:
203: See "slepc/include/slepcsvd.h" for available methods. The default
204: is SVDCROSS.
206: Normally, it is best to use the SVDSetFromOptions() command and
207: then set the SVD type from the options database rather than by using
208: this routine. Using the options database provides the user with
209: maximum flexibility in evaluating the different available methods.
210: The SVDSetType() routine is provided for those situations where it
211: is necessary to set the iterative solver independently of the command
212: line or options database.
214: Level: intermediate
216: .seealso: SVDType
217: @*/
218: PetscErrorCode SVDSetType(SVD svd,SVDType type)
219: {
220: PetscErrorCode (*r)(SVD);
221: PetscBool match;
223: PetscFunctionBegin;
225: PetscAssertPointer(type,2);
227: PetscCall(PetscObjectTypeCompare((PetscObject)svd,type,&match));
228: if (match) PetscFunctionReturn(PETSC_SUCCESS);
230: PetscCall(PetscFunctionListFind(SVDList,type,&r));
231: PetscCheck(r,PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown SVD type given: %s",type);
233: PetscTryTypeMethod(svd,destroy);
234: PetscCall(PetscMemzero(svd->ops,sizeof(struct _SVDOps)));
236: svd->state = SVD_STATE_INITIAL;
237: PetscCall(PetscObjectChangeTypeName((PetscObject)svd,type));
238: PetscCall((*r)(svd));
239: PetscFunctionReturn(PETSC_SUCCESS);
240: }
242: /*@C
243: SVDGetType - Gets the SVD type as a string from the SVD object.
245: Not Collective
247: Input Parameter:
248: . svd - the singular value solver context
250: Output Parameter:
251: . type - name of SVD method
253: Level: intermediate
255: .seealso: SVDSetType()
256: @*/
257: PetscErrorCode SVDGetType(SVD svd,SVDType *type)
258: {
259: PetscFunctionBegin;
261: PetscAssertPointer(type,2);
262: *type = ((PetscObject)svd)->type_name;
263: PetscFunctionReturn(PETSC_SUCCESS);
264: }
266: /*@C
267: SVDRegister - Adds a method to the singular value solver package.
269: Not Collective
271: Input Parameters:
272: + name - name of a new user-defined solver
273: - function - routine to create the solver context
275: Notes:
276: SVDRegister() may be called multiple times to add several user-defined solvers.
278: Example Usage:
279: .vb
280: SVDRegister("my_solver",MySolverCreate);
281: .ve
283: Then, your solver can be chosen with the procedural interface via
284: $ SVDSetType(svd,"my_solver")
285: or at runtime via the option
286: $ -svd_type my_solver
288: Level: advanced
290: .seealso: SVDRegisterAll()
291: @*/
292: PetscErrorCode SVDRegister(const char *name,PetscErrorCode (*function)(SVD))
293: {
294: PetscFunctionBegin;
295: PetscCall(SVDInitializePackage());
296: PetscCall(PetscFunctionListAdd(&SVDList,name,function));
297: PetscFunctionReturn(PETSC_SUCCESS);
298: }
300: /*@C
301: SVDMonitorRegister - Adds SVD monitor routine.
303: Not Collective
305: Input Parameters:
306: + name - name of a new monitor routine
307: . vtype - a PetscViewerType for the output
308: . format - a PetscViewerFormat for the output
309: . monitor - monitor routine
310: . create - creation routine, or NULL
311: - destroy - destruction routine, or NULL
313: Notes:
314: SVDMonitorRegister() may be called multiple times to add several user-defined monitors.
316: Example Usage:
317: .vb
318: SVDMonitorRegister("my_monitor",PETSCVIEWERASCII,PETSC_VIEWER_ASCII_INFO_DETAIL,MyMonitor,NULL,NULL);
319: .ve
321: Then, your monitor can be chosen with the procedural interface via
322: $ SVDMonitorSetFromOptions(svd,"-svd_monitor_my_monitor","my_monitor",NULL)
323: or at runtime via the option
324: $ -svd_monitor_my_monitor
326: Level: advanced
328: .seealso: SVDMonitorRegisterAll()
329: @*/
330: 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**))
331: {
332: char key[PETSC_MAX_PATH_LEN];
334: PetscFunctionBegin;
335: PetscCall(SVDInitializePackage());
336: PetscCall(SlepcMonitorMakeKey_Internal(name,vtype,format,key));
337: PetscCall(PetscFunctionListAdd(&SVDMonitorList,key,monitor));
338: if (create) PetscCall(PetscFunctionListAdd(&SVDMonitorCreateList,key,create));
339: if (destroy) PetscCall(PetscFunctionListAdd(&SVDMonitorDestroyList,key,destroy));
340: PetscFunctionReturn(PETSC_SUCCESS);
341: }
343: /*@
344: SVDSetBV - Associates basis vectors objects to the singular value solver.
346: Collective
348: Input Parameters:
349: + svd - singular value solver context obtained from SVDCreate()
350: . V - the basis vectors object for right singular vectors
351: - U - the basis vectors object for left singular vectors
353: Note:
354: Use SVDGetBV() to retrieve the basis vectors contexts (for example,
355: to free them at the end of the computations).
357: Level: advanced
359: .seealso: SVDGetBV()
360: @*/
361: PetscErrorCode SVDSetBV(SVD svd,BV V,BV U)
362: {
363: PetscFunctionBegin;
365: if (V) {
367: PetscCheckSameComm(svd,1,V,2);
368: PetscCall(PetscObjectReference((PetscObject)V));
369: PetscCall(BVDestroy(&svd->V));
370: svd->V = V;
371: }
372: if (U) {
374: PetscCheckSameComm(svd,1,U,3);
375: PetscCall(PetscObjectReference((PetscObject)U));
376: PetscCall(BVDestroy(&svd->U));
377: svd->U = U;
378: }
379: PetscFunctionReturn(PETSC_SUCCESS);
380: }
382: /*@
383: SVDGetBV - Obtain the basis vectors objects associated to the singular
384: value solver object.
386: Not Collective
388: Input Parameter:
389: . svd - singular value solver context obtained from SVDCreate()
391: Output Parameters:
392: + V - basis vectors context for right singular vectors
393: - U - basis vectors context for left singular vectors
395: Level: advanced
397: .seealso: SVDSetBV()
398: @*/
399: PetscErrorCode SVDGetBV(SVD svd,BV *V,BV *U)
400: {
401: PetscFunctionBegin;
403: if (V) {
404: if (!svd->V) {
405: PetscCall(BVCreate(PetscObjectComm((PetscObject)svd),&svd->V));
406: PetscCall(PetscObjectIncrementTabLevel((PetscObject)svd->V,(PetscObject)svd,0));
407: PetscCall(PetscObjectSetOptions((PetscObject)svd->V,((PetscObject)svd)->options));
408: }
409: *V = svd->V;
410: }
411: if (U) {
412: if (!svd->U) {
413: PetscCall(BVCreate(PetscObjectComm((PetscObject)svd),&svd->U));
414: PetscCall(PetscObjectIncrementTabLevel((PetscObject)svd->U,(PetscObject)svd,0));
415: PetscCall(PetscObjectSetOptions((PetscObject)svd->U,((PetscObject)svd)->options));
416: }
417: *U = svd->U;
418: }
419: PetscFunctionReturn(PETSC_SUCCESS);
420: }
422: /*@
423: SVDSetDS - Associates a direct solver object to the singular value solver.
425: Collective
427: Input Parameters:
428: + svd - singular value solver context obtained from SVDCreate()
429: - ds - the direct solver object
431: Note:
432: Use SVDGetDS() to retrieve the direct solver context (for example,
433: to free it at the end of the computations).
435: Level: advanced
437: .seealso: SVDGetDS()
438: @*/
439: PetscErrorCode SVDSetDS(SVD svd,DS ds)
440: {
441: PetscFunctionBegin;
444: PetscCheckSameComm(svd,1,ds,2);
445: PetscCall(PetscObjectReference((PetscObject)ds));
446: PetscCall(DSDestroy(&svd->ds));
447: svd->ds = ds;
448: PetscFunctionReturn(PETSC_SUCCESS);
449: }
451: /*@
452: SVDGetDS - Obtain the direct solver object associated to the singular value
453: solver object.
455: Not Collective
457: Input Parameters:
458: . svd - singular value solver context obtained from SVDCreate()
460: Output Parameter:
461: . ds - direct solver context
463: Level: advanced
465: .seealso: SVDSetDS()
466: @*/
467: PetscErrorCode SVDGetDS(SVD svd,DS *ds)
468: {
469: PetscFunctionBegin;
471: PetscAssertPointer(ds,2);
472: if (!svd->ds) {
473: PetscCall(DSCreate(PetscObjectComm((PetscObject)svd),&svd->ds));
474: PetscCall(PetscObjectIncrementTabLevel((PetscObject)svd->ds,(PetscObject)svd,0));
475: PetscCall(PetscObjectSetOptions((PetscObject)svd->ds,((PetscObject)svd)->options));
476: }
477: *ds = svd->ds;
478: PetscFunctionReturn(PETSC_SUCCESS);
479: }