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