Line data Source code
1 : /*
2 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3 : SLEPc - Scalable Library for Eigenvalue Problem Computations
4 : Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
5 :
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 : */
13 :
14 : #include <slepc/private/svdimpl.h> /*I "slepcsvd.h" I*/
15 :
16 : /* Logging support */
17 : PetscClassId SVD_CLASSID = 0;
18 : PetscLogEvent SVD_SetUp = 0,SVD_Solve = 0;
19 :
20 : /* List of registered SVD routines */
21 : PetscFunctionList SVDList = NULL;
22 : PetscBool SVDRegisterAllCalled = PETSC_FALSE;
23 :
24 : /* List of registered SVD monitors */
25 : PetscFunctionList SVDMonitorList = NULL;
26 : PetscFunctionList SVDMonitorCreateList = NULL;
27 : PetscFunctionList SVDMonitorDestroyList = NULL;
28 : PetscBool SVDMonitorRegisterAllCalled = PETSC_FALSE;
29 :
30 : /*@
31 : SVDCreate - Creates the default SVD context.
32 :
33 : Collective
34 :
35 : Input Parameter:
36 : . comm - MPI communicator
37 :
38 : Output Parameter:
39 : . outsvd - location to put the SVD context
40 :
41 : Note:
42 : The default SVD type is SVDCROSS
43 :
44 : Level: beginner
45 :
46 : .seealso: SVDSetUp(), SVDSolve(), SVDDestroy(), SVD
47 : @*/
48 200 : PetscErrorCode SVDCreate(MPI_Comm comm,SVD *outsvd)
49 : {
50 200 : SVD svd;
51 :
52 200 : PetscFunctionBegin;
53 200 : PetscAssertPointer(outsvd,2);
54 200 : *outsvd = NULL;
55 200 : PetscCall(SVDInitializePackage());
56 200 : PetscCall(SlepcHeaderCreate(svd,SVD_CLASSID,"SVD","Singular Value Decomposition","SVD",comm,SVDDestroy,SVDView));
57 :
58 200 : svd->OP = NULL;
59 200 : svd->OPb = NULL;
60 200 : svd->omega = NULL;
61 200 : svd->max_it = PETSC_DEFAULT;
62 200 : svd->nsv = 1;
63 200 : svd->ncv = PETSC_DEFAULT;
64 200 : svd->mpd = PETSC_DEFAULT;
65 200 : svd->nini = 0;
66 200 : svd->ninil = 0;
67 200 : svd->tol = PETSC_DEFAULT;
68 200 : svd->conv = (SVDConv)-1;
69 200 : svd->stop = SVD_STOP_BASIC;
70 200 : svd->which = SVD_LARGEST;
71 200 : svd->problem_type = (SVDProblemType)0;
72 200 : svd->impltrans = PETSC_FALSE;
73 200 : svd->trackall = PETSC_FALSE;
74 :
75 200 : svd->converged = NULL;
76 200 : svd->convergeduser = NULL;
77 200 : svd->convergeddestroy = NULL;
78 200 : svd->stopping = SVDStoppingBasic;
79 200 : svd->stoppinguser = NULL;
80 200 : svd->stoppingdestroy = NULL;
81 200 : svd->convergedctx = NULL;
82 200 : svd->stoppingctx = NULL;
83 200 : svd->numbermonitors = 0;
84 :
85 200 : svd->ds = NULL;
86 200 : svd->U = NULL;
87 200 : svd->V = NULL;
88 200 : svd->A = NULL;
89 200 : svd->B = NULL;
90 200 : svd->AT = NULL;
91 200 : svd->BT = NULL;
92 200 : svd->IS = NULL;
93 200 : svd->ISL = NULL;
94 200 : svd->sigma = NULL;
95 200 : svd->errest = NULL;
96 200 : svd->sign = NULL;
97 200 : svd->perm = NULL;
98 200 : svd->nworkl = 0;
99 200 : svd->nworkr = 0;
100 200 : svd->workl = NULL;
101 200 : svd->workr = NULL;
102 200 : svd->data = NULL;
103 :
104 200 : svd->state = SVD_STATE_INITIAL;
105 200 : svd->nconv = 0;
106 200 : svd->its = 0;
107 200 : svd->leftbasis = PETSC_FALSE;
108 200 : svd->swapped = PETSC_FALSE;
109 200 : svd->expltrans = PETSC_FALSE;
110 200 : svd->nrma = 0.0;
111 200 : svd->nrmb = 0.0;
112 200 : svd->isgeneralized = PETSC_FALSE;
113 200 : svd->reason = SVD_CONVERGED_ITERATING;
114 :
115 200 : PetscCall(PetscNew(&svd->sc));
116 200 : *outsvd = svd;
117 200 : PetscFunctionReturn(PETSC_SUCCESS);
118 : }
119 :
120 : /*@
121 : SVDReset - Resets the SVD context to the initial state (prior to setup)
122 : and destroys any allocated Vecs and Mats.
123 :
124 : Collective
125 :
126 : Input Parameter:
127 : . svd - singular value solver context obtained from SVDCreate()
128 :
129 : Level: advanced
130 :
131 : .seealso: SVDDestroy()
132 : @*/
133 206 : PetscErrorCode SVDReset(SVD svd)
134 : {
135 206 : PetscFunctionBegin;
136 206 : if (svd) PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
137 0 : if (!svd) PetscFunctionReturn(PETSC_SUCCESS);
138 206 : PetscTryTypeMethod(svd,reset);
139 206 : PetscCall(MatDestroy(&svd->OP));
140 206 : PetscCall(MatDestroy(&svd->OPb));
141 206 : PetscCall(VecDestroy(&svd->omega));
142 206 : PetscCall(MatDestroy(&svd->A));
143 206 : PetscCall(MatDestroy(&svd->B));
144 206 : PetscCall(MatDestroy(&svd->AT));
145 206 : PetscCall(MatDestroy(&svd->BT));
146 206 : PetscCall(BVDestroy(&svd->U));
147 206 : PetscCall(BVDestroy(&svd->V));
148 206 : PetscCall(VecDestroyVecs(svd->nworkl,&svd->workl));
149 206 : svd->nworkl = 0;
150 206 : PetscCall(VecDestroyVecs(svd->nworkr,&svd->workr));
151 206 : svd->nworkr = 0;
152 206 : svd->swapped = PETSC_FALSE;
153 206 : svd->state = SVD_STATE_INITIAL;
154 206 : PetscFunctionReturn(PETSC_SUCCESS);
155 : }
156 :
157 : /*@C
158 : SVDDestroy - Destroys the SVD context.
159 :
160 : Collective
161 :
162 : Input Parameter:
163 : . svd - singular value solver context obtained from SVDCreate()
164 :
165 : Level: beginner
166 :
167 : .seealso: SVDCreate(), SVDSetUp(), SVDSolve()
168 : @*/
169 200 : PetscErrorCode SVDDestroy(SVD *svd)
170 : {
171 200 : PetscFunctionBegin;
172 200 : if (!*svd) PetscFunctionReturn(PETSC_SUCCESS);
173 200 : PetscValidHeaderSpecific(*svd,SVD_CLASSID,1);
174 200 : if (--((PetscObject)*svd)->refct > 0) { *svd = NULL; PetscFunctionReturn(PETSC_SUCCESS); }
175 200 : PetscCall(SVDReset(*svd));
176 200 : PetscTryTypeMethod(*svd,destroy);
177 200 : if ((*svd)->sigma) PetscCall(PetscFree3((*svd)->sigma,(*svd)->perm,(*svd)->errest));
178 200 : if ((*svd)->sign) PetscCall(PetscFree((*svd)->sign));
179 200 : PetscCall(DSDestroy(&(*svd)->ds));
180 200 : PetscCall(PetscFree((*svd)->sc));
181 : /* just in case the initial vectors have not been used */
182 200 : PetscCall(SlepcBasisDestroy_Private(&(*svd)->nini,&(*svd)->IS));
183 200 : PetscCall(SlepcBasisDestroy_Private(&(*svd)->ninil,&(*svd)->ISL));
184 200 : PetscCall(SVDMonitorCancel(*svd));
185 200 : PetscCall(PetscHeaderDestroy(svd));
186 200 : PetscFunctionReturn(PETSC_SUCCESS);
187 : }
188 :
189 : /*@C
190 : SVDSetType - Selects the particular solver to be used in the SVD object.
191 :
192 : Logically Collective
193 :
194 : Input Parameters:
195 : + svd - the singular value solver context
196 : - type - a known method
197 :
198 : Options Database Key:
199 : . -svd_type <method> - Sets the method; use -help for a list
200 : of available methods
201 :
202 : Notes:
203 : See "slepc/include/slepcsvd.h" for available methods. The default
204 : is SVDCROSS.
205 :
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.
213 :
214 : Level: intermediate
215 :
216 : .seealso: SVDType
217 : @*/
218 201 : PetscErrorCode SVDSetType(SVD svd,SVDType type)
219 : {
220 201 : PetscErrorCode (*r)(SVD);
221 201 : PetscBool match;
222 :
223 201 : PetscFunctionBegin;
224 201 : PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
225 201 : PetscAssertPointer(type,2);
226 :
227 201 : PetscCall(PetscObjectTypeCompare((PetscObject)svd,type,&match));
228 201 : if (match) PetscFunctionReturn(PETSC_SUCCESS);
229 :
230 201 : PetscCall(PetscFunctionListFind(SVDList,type,&r));
231 201 : PetscCheck(r,PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown SVD type given: %s",type);
232 :
233 201 : PetscTryTypeMethod(svd,destroy);
234 201 : PetscCall(PetscMemzero(svd->ops,sizeof(struct _SVDOps)));
235 :
236 201 : svd->state = SVD_STATE_INITIAL;
237 201 : PetscCall(PetscObjectChangeTypeName((PetscObject)svd,type));
238 201 : PetscCall((*r)(svd));
239 201 : PetscFunctionReturn(PETSC_SUCCESS);
240 : }
241 :
242 : /*@C
243 : SVDGetType - Gets the SVD type as a string from the SVD object.
244 :
245 : Not Collective
246 :
247 : Input Parameter:
248 : . svd - the singular value solver context
249 :
250 : Output Parameter:
251 : . type - name of SVD method
252 :
253 : Level: intermediate
254 :
255 : .seealso: SVDSetType()
256 : @*/
257 44 : PetscErrorCode SVDGetType(SVD svd,SVDType *type)
258 : {
259 44 : PetscFunctionBegin;
260 44 : PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
261 44 : PetscAssertPointer(type,2);
262 44 : *type = ((PetscObject)svd)->type_name;
263 44 : PetscFunctionReturn(PETSC_SUCCESS);
264 : }
265 :
266 : /*@C
267 : SVDRegister - Adds a method to the singular value solver package.
268 :
269 : Not Collective
270 :
271 : Input Parameters:
272 : + name - name of a new user-defined solver
273 : - function - routine to create the solver context
274 :
275 : Notes:
276 : SVDRegister() may be called multiple times to add several user-defined solvers.
277 :
278 : Example Usage:
279 : .vb
280 : SVDRegister("my_solver",MySolverCreate);
281 : .ve
282 :
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
287 :
288 : Level: advanced
289 :
290 : .seealso: SVDRegisterAll()
291 : @*/
292 1407 : PetscErrorCode SVDRegister(const char *name,PetscErrorCode (*function)(SVD))
293 : {
294 1407 : PetscFunctionBegin;
295 1407 : PetscCall(SVDInitializePackage());
296 1407 : PetscCall(PetscFunctionListAdd(&SVDList,name,function));
297 1407 : PetscFunctionReturn(PETSC_SUCCESS);
298 : }
299 :
300 : /*@C
301 : SVDMonitorRegister - Adds SVD monitor routine.
302 :
303 : Not Collective
304 :
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
312 :
313 : Notes:
314 : SVDMonitorRegister() may be called multiple times to add several user-defined monitors.
315 :
316 : Example Usage:
317 : .vb
318 : SVDMonitorRegister("my_monitor",PETSCVIEWERASCII,PETSC_VIEWER_ASCII_INFO_DETAIL,MyMonitor,NULL,NULL);
319 : .ve
320 :
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
325 :
326 : Level: advanced
327 :
328 : .seealso: SVDMonitorRegisterAll()
329 : @*/
330 1407 : 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 1407 : char key[PETSC_MAX_PATH_LEN];
333 :
334 1407 : PetscFunctionBegin;
335 1407 : PetscCall(SVDInitializePackage());
336 1407 : PetscCall(SlepcMonitorMakeKey_Internal(name,vtype,format,key));
337 1407 : PetscCall(PetscFunctionListAdd(&SVDMonitorList,key,monitor));
338 1407 : if (create) PetscCall(PetscFunctionListAdd(&SVDMonitorCreateList,key,create));
339 1407 : if (destroy) PetscCall(PetscFunctionListAdd(&SVDMonitorDestroyList,key,destroy));
340 1407 : PetscFunctionReturn(PETSC_SUCCESS);
341 : }
342 :
343 : /*@
344 : SVDSetBV - Associates basis vectors objects to the singular value solver.
345 :
346 : Collective
347 :
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
352 :
353 : Note:
354 : Use SVDGetBV() to retrieve the basis vectors contexts (for example,
355 : to free them at the end of the computations).
356 :
357 : Level: advanced
358 :
359 : .seealso: SVDGetBV()
360 : @*/
361 3 : PetscErrorCode SVDSetBV(SVD svd,BV V,BV U)
362 : {
363 3 : PetscFunctionBegin;
364 3 : PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
365 3 : if (V) {
366 3 : PetscValidHeaderSpecific(V,BV_CLASSID,2);
367 3 : PetscCheckSameComm(svd,1,V,2);
368 3 : PetscCall(PetscObjectReference((PetscObject)V));
369 3 : PetscCall(BVDestroy(&svd->V));
370 3 : svd->V = V;
371 : }
372 3 : if (U) {
373 3 : PetscValidHeaderSpecific(U,BV_CLASSID,3);
374 3 : PetscCheckSameComm(svd,1,U,3);
375 3 : PetscCall(PetscObjectReference((PetscObject)U));
376 3 : PetscCall(BVDestroy(&svd->U));
377 3 : svd->U = U;
378 : }
379 3 : PetscFunctionReturn(PETSC_SUCCESS);
380 : }
381 :
382 : /*@
383 : SVDGetBV - Obtain the basis vectors objects associated to the singular
384 : value solver object.
385 :
386 : Not Collective
387 :
388 : Input Parameter:
389 : . svd - singular value solver context obtained from SVDCreate()
390 :
391 : Output Parameters:
392 : + V - basis vectors context for right singular vectors
393 : - U - basis vectors context for left singular vectors
394 :
395 : Level: advanced
396 :
397 : .seealso: SVDSetBV()
398 : @*/
399 404 : PetscErrorCode SVDGetBV(SVD svd,BV *V,BV *U)
400 : {
401 404 : PetscFunctionBegin;
402 404 : PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
403 404 : if (V) {
404 203 : if (!svd->V) {
405 203 : PetscCall(BVCreate(PetscObjectComm((PetscObject)svd),&svd->V));
406 203 : PetscCall(PetscObjectIncrementTabLevel((PetscObject)svd->V,(PetscObject)svd,0));
407 203 : PetscCall(PetscObjectSetOptions((PetscObject)svd->V,((PetscObject)svd)->options));
408 : }
409 203 : *V = svd->V;
410 : }
411 404 : if (U) {
412 201 : if (!svd->U) {
413 201 : PetscCall(BVCreate(PetscObjectComm((PetscObject)svd),&svd->U));
414 201 : PetscCall(PetscObjectIncrementTabLevel((PetscObject)svd->U,(PetscObject)svd,0));
415 201 : PetscCall(PetscObjectSetOptions((PetscObject)svd->U,((PetscObject)svd)->options));
416 : }
417 201 : *U = svd->U;
418 : }
419 404 : PetscFunctionReturn(PETSC_SUCCESS);
420 : }
421 :
422 : /*@
423 : SVDSetDS - Associates a direct solver object to the singular value solver.
424 :
425 : Collective
426 :
427 : Input Parameters:
428 : + svd - singular value solver context obtained from SVDCreate()
429 : - ds - the direct solver object
430 :
431 : Note:
432 : Use SVDGetDS() to retrieve the direct solver context (for example,
433 : to free it at the end of the computations).
434 :
435 : Level: advanced
436 :
437 : .seealso: SVDGetDS()
438 : @*/
439 0 : PetscErrorCode SVDSetDS(SVD svd,DS ds)
440 : {
441 0 : PetscFunctionBegin;
442 0 : PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
443 0 : PetscValidHeaderSpecific(ds,DS_CLASSID,2);
444 0 : PetscCheckSameComm(svd,1,ds,2);
445 0 : PetscCall(PetscObjectReference((PetscObject)ds));
446 0 : PetscCall(DSDestroy(&svd->ds));
447 0 : svd->ds = ds;
448 0 : PetscFunctionReturn(PETSC_SUCCESS);
449 : }
450 :
451 : /*@
452 : SVDGetDS - Obtain the direct solver object associated to the singular value
453 : solver object.
454 :
455 : Not Collective
456 :
457 : Input Parameters:
458 : . svd - singular value solver context obtained from SVDCreate()
459 :
460 : Output Parameter:
461 : . ds - direct solver context
462 :
463 : Level: advanced
464 :
465 : .seealso: SVDSetDS()
466 : @*/
467 200 : PetscErrorCode SVDGetDS(SVD svd,DS *ds)
468 : {
469 200 : PetscFunctionBegin;
470 200 : PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
471 200 : PetscAssertPointer(ds,2);
472 200 : if (!svd->ds) {
473 200 : PetscCall(DSCreate(PetscObjectComm((PetscObject)svd),&svd->ds));
474 200 : PetscCall(PetscObjectIncrementTabLevel((PetscObject)svd->ds,(PetscObject)svd,0));
475 200 : PetscCall(PetscObjectSetOptions((PetscObject)svd->ds,((PetscObject)svd)->options));
476 : }
477 200 : *ds = svd->ds;
478 200 : PetscFunctionReturn(PETSC_SUCCESS);
479 : }
|