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 : SVD routines for setting up the solver
12 : */
13 :
14 : #include <slepc/private/svdimpl.h> /*I "slepcsvd.h" I*/
15 :
16 : /*@
17 : SVDSetOperators - Set the matrices associated with the singular value problem.
18 :
19 : Collective
20 :
21 : Input Parameters:
22 : + svd - the singular value solver context
23 : . A - the matrix associated with the singular value problem
24 : - B - the second matrix in the case of GSVD
25 :
26 : Level: beginner
27 :
28 : .seealso: SVDSolve(), SVDGetOperators()
29 : @*/
30 241 : PetscErrorCode SVDSetOperators(SVD svd,Mat A,Mat B)
31 : {
32 241 : PetscInt Ma,Na,Mb,Nb,ma,na,mb,nb,M0,N0,m0,n0;
33 241 : PetscBool samesize=PETSC_TRUE;
34 :
35 241 : PetscFunctionBegin;
36 241 : PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
37 241 : PetscValidHeaderSpecific(A,MAT_CLASSID,2);
38 241 : if (B) PetscValidHeaderSpecific(B,MAT_CLASSID,3);
39 241 : PetscCheckSameComm(svd,1,A,2);
40 241 : if (B) PetscCheckSameComm(svd,1,B,3);
41 :
42 : /* Check matrix sizes */
43 241 : PetscCall(MatGetSize(A,&Ma,&Na));
44 241 : PetscCall(MatGetLocalSize(A,&ma,&na));
45 241 : if (svd->OP) {
46 25 : PetscCall(MatGetSize(svd->OP,&M0,&N0));
47 25 : PetscCall(MatGetLocalSize(svd->OP,&m0,&n0));
48 25 : if (M0!=Ma || N0!=Na || m0!=ma || n0!=na) samesize = PETSC_FALSE;
49 : }
50 241 : if (B) {
51 89 : PetscCall(MatGetSize(B,&Mb,&Nb));
52 89 : PetscCall(MatGetLocalSize(B,&mb,&nb));
53 89 : PetscCheck(Na==Nb,PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_WRONG,"Different number of columns in A (%" PetscInt_FMT ") and B (%" PetscInt_FMT ")",Na,Nb);
54 89 : PetscCheck(na==nb,PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_WRONG,"Different local column size in A (%" PetscInt_FMT ") and B (%" PetscInt_FMT ")",na,nb);
55 89 : if (svd->OPb) {
56 11 : PetscCall(MatGetSize(svd->OPb,&M0,&N0));
57 11 : PetscCall(MatGetLocalSize(svd->OPb,&m0,&n0));
58 11 : if (M0!=Mb || N0!=Nb || m0!=mb || n0!=nb) samesize = PETSC_FALSE;
59 : }
60 : }
61 :
62 241 : PetscCall(PetscObjectReference((PetscObject)A));
63 241 : if (B) PetscCall(PetscObjectReference((PetscObject)B));
64 241 : if (svd->state && !samesize) PetscCall(SVDReset(svd));
65 : else {
66 241 : PetscCall(MatDestroy(&svd->OP));
67 241 : PetscCall(MatDestroy(&svd->OPb));
68 241 : PetscCall(MatDestroy(&svd->A));
69 241 : PetscCall(MatDestroy(&svd->B));
70 241 : PetscCall(MatDestroy(&svd->AT));
71 241 : PetscCall(MatDestroy(&svd->BT));
72 : }
73 241 : svd->nrma = 0.0;
74 241 : svd->nrmb = 0.0;
75 241 : svd->OP = A;
76 241 : svd->OPb = B;
77 241 : svd->state = SVD_STATE_INITIAL;
78 241 : PetscFunctionReturn(PETSC_SUCCESS);
79 : }
80 :
81 : /*@
82 : SVDGetOperators - Get the matrices associated with the singular value problem.
83 :
84 : Not Collective
85 :
86 : Input Parameter:
87 : . svd - the singular value solver context
88 :
89 : Output Parameters:
90 : + A - the matrix associated with the singular value problem
91 : - B - the second matrix in the case of GSVD
92 :
93 : Level: intermediate
94 :
95 : .seealso: SVDSolve(), SVDSetOperators()
96 : @*/
97 11 : PetscErrorCode SVDGetOperators(SVD svd,Mat *A,Mat *B)
98 : {
99 11 : PetscFunctionBegin;
100 11 : PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
101 11 : if (A) *A = svd->OP;
102 11 : if (B) *B = svd->OPb;
103 11 : PetscFunctionReturn(PETSC_SUCCESS);
104 : }
105 :
106 : /*@
107 : SVDSetSignature - Set the signature matrix defining a hyperbolic singular value problem.
108 :
109 : Collective
110 :
111 : Input Parameters:
112 : + svd - the singular value solver context
113 : - omega - a vector containing the diagonal elements of the signature matrix (or NULL)
114 :
115 : Notes:
116 : The signature matrix is relevant only for hyperbolic problems (HSVD).
117 : Use NULL to reset a previously set signature.
118 :
119 : Level: intermediate
120 :
121 : .seealso: SVDSetProblemType(), SVDSetOperators(), SVDGetSignature()
122 : @*/
123 28 : PetscErrorCode SVDSetSignature(SVD svd,Vec omega)
124 : {
125 28 : PetscInt N,Ma,n,ma;
126 :
127 28 : PetscFunctionBegin;
128 28 : PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
129 28 : if (omega) {
130 28 : PetscValidHeaderSpecific(omega,VEC_CLASSID,2);
131 28 : PetscCheckSameComm(svd,1,omega,2);
132 : }
133 :
134 28 : if (omega && svd->OP) { /* Check sizes */
135 27 : PetscCall(VecGetSize(omega,&N));
136 27 : PetscCall(VecGetLocalSize(omega,&n));
137 27 : PetscCall(MatGetSize(svd->OP,&Ma,NULL));
138 27 : PetscCall(MatGetLocalSize(svd->OP,&ma,NULL));
139 27 : PetscCheck(N==Ma,PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_SIZ,"Global size of signature (%" PetscInt_FMT ") does not match the row size of A (%" PetscInt_FMT ")",N,Ma);
140 27 : PetscCheck(n==ma,PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_SIZ,"Local size of signature (%" PetscInt_FMT ") does not match the local row size of A (%" PetscInt_FMT ")",n,ma);
141 : }
142 :
143 28 : PetscCall(VecDestroy(&svd->omega));
144 28 : if (omega) {
145 28 : PetscCall(VecDuplicate(omega,&svd->omega));
146 28 : PetscCall(VecCopy(omega,svd->omega));
147 : }
148 28 : svd->state = SVD_STATE_INITIAL;
149 28 : PetscFunctionReturn(PETSC_SUCCESS);
150 : }
151 :
152 : /*@
153 : SVDGetSignature - Get the signature matrix defining a hyperbolic singular value problem.
154 :
155 : Not Collective
156 :
157 : Input Parameter:
158 : . svd - the singular value solver context
159 :
160 : Output Parameter:
161 : . omega - a vector containing the diagonal elements of the signature matrix
162 :
163 : Notes:
164 : The signature matrix is relevant only for hyperbolic problems (HSVD).
165 : If no signature has been set, this function will return a vector of all ones.
166 :
167 : The user should pass a previously created Vec with the appropriate dimension.
168 :
169 : Level: intermediate
170 :
171 : .seealso: SVDSetSignature()
172 : @*/
173 1 : PetscErrorCode SVDGetSignature(SVD svd,Vec omega)
174 : {
175 1 : PetscFunctionBegin;
176 1 : PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
177 1 : PetscValidHeaderSpecific(omega,VEC_CLASSID,2);
178 1 : if (svd->omega) PetscCall(VecCopy(svd->omega,omega));
179 0 : else PetscCall(VecSet(omega,1.0));
180 1 : PetscFunctionReturn(PETSC_SUCCESS);
181 : }
182 :
183 : /*@
184 : SVDSetDSType - Sets the type of the internal DS object based on the current
185 : settings of the singular value solver.
186 :
187 : Collective
188 :
189 : Input Parameter:
190 : . svd - singular value solver context
191 :
192 : Note:
193 : This function need not be called explicitly, since it will be called at
194 : both SVDSetFromOptions() and SVDSetUp().
195 :
196 : Level: developer
197 :
198 : .seealso: SVDSetFromOptions(), SVDSetUp()
199 : @*/
200 459 : PetscErrorCode SVDSetDSType(SVD svd)
201 : {
202 459 : PetscFunctionBegin;
203 459 : PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
204 459 : PetscTryTypeMethod(svd,setdstype);
205 459 : PetscFunctionReturn(PETSC_SUCCESS);
206 : }
207 :
208 : /*@
209 : SVDSetUp - Sets up all the internal data structures necessary for the
210 : execution of the singular value solver.
211 :
212 : Collective
213 :
214 : Input Parameter:
215 : . svd - singular value solver context
216 :
217 : Notes:
218 : This function need not be called explicitly in most cases, since SVDSolve()
219 : calls it. It can be useful when one wants to measure the set-up time
220 : separately from the solve time.
221 :
222 : Level: developer
223 :
224 : .seealso: SVDCreate(), SVDSolve(), SVDDestroy()
225 : @*/
226 257 : PetscErrorCode SVDSetUp(SVD svd)
227 : {
228 257 : PetscBool flg;
229 257 : PetscInt M,N,P=0,k,maxnsol,m,Nom,nom;
230 257 : SlepcSC sc;
231 257 : Vec *T;
232 257 : BV bv;
233 257 : SVDStoppingCtx ctx;
234 :
235 257 : PetscFunctionBegin;
236 257 : PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
237 257 : if (svd->state) PetscFunctionReturn(PETSC_SUCCESS);
238 257 : PetscCall(PetscLogEventBegin(SVD_SetUp,svd,0,0,0));
239 :
240 : /* reset the convergence flag from the previous solves */
241 257 : svd->reason = SVD_CONVERGED_ITERATING;
242 :
243 : /* set default solver type (SVDSetFromOptions was not called) */
244 257 : if (!((PetscObject)svd)->type_name) PetscCall(SVDSetType(svd,SVDCROSS));
245 257 : if (!svd->ds) PetscCall(SVDGetDS(svd,&svd->ds));
246 :
247 : /* check matrices */
248 257 : PetscCheck(svd->OP,PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_WRONGSTATE,"SVDSetOperators() must be called first");
249 :
250 257 : if (!svd->problem_type) { /* set default problem type */
251 122 : if (svd->OPb) {
252 18 : PetscCheck(!svd->omega,PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"There is no support yet for generalized hyperbolic problems");
253 18 : PetscCall(SVDSetProblemType(svd,SVD_GENERALIZED));
254 : } else {
255 104 : if (svd->omega) PetscCall(SVDSetProblemType(svd,SVD_HYPERBOLIC));
256 103 : else PetscCall(SVDSetProblemType(svd,SVD_STANDARD));
257 : }
258 : } else { /* check consistency of problem type set by user */
259 135 : if (svd->OPb) {
260 71 : PetscCheck(svd->isgeneralized,PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_INCOMP,"Inconsistent SVD state: the problem type does not match the number of matrices");
261 71 : PetscCheck(!svd->omega,PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"There is no support yet for generalized hyperbolic problems");
262 : } else {
263 64 : PetscCheck(!svd->isgeneralized,PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_INCOMP,"Inconsistent SVD state: the problem type does not match the number of matrices");
264 64 : if (svd->omega) PetscCheck(svd->ishyperbolic,PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_INCOMP,"Inconsistent SVD state: the problem type must be set to hyperbolic when passing a signature with SVDSetSignature()");
265 37 : else PetscCheck(!svd->ishyperbolic,PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_INCOMP,"Inconsistent SVD state: a hyperbolic problem requires passing a signature with SVDSetSignature()");
266 : }
267 : }
268 :
269 : /* set DS type once the default problem type has been determined */
270 257 : PetscCall(SVDSetDSType(svd));
271 :
272 : /* determine how to handle the transpose */
273 257 : svd->expltrans = PETSC_TRUE;
274 257 : if (svd->impltrans) svd->expltrans = PETSC_FALSE;
275 : else {
276 246 : PetscCall(MatHasOperation(svd->OP,MATOP_TRANSPOSE,&flg));
277 246 : if (!flg) svd->expltrans = PETSC_FALSE;
278 : else {
279 245 : PetscCall(PetscObjectTypeCompareAny((PetscObject)svd,&flg,SVDLAPACK,SVDSCALAPACK,SVDKSVD,SVDELEMENTAL,""));
280 245 : if (flg) svd->expltrans = PETSC_FALSE;
281 : }
282 : }
283 :
284 : /* get matrix dimensions */
285 257 : PetscCall(MatGetSize(svd->OP,&M,&N));
286 257 : if (svd->isgeneralized) {
287 89 : PetscCall(MatGetSize(svd->OPb,&P,NULL));
288 89 : PetscCheck(M+P>=N,PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"The case when [A;B] has less rows than columns is not supported");
289 168 : } else if (svd->ishyperbolic) {
290 28 : PetscCall(VecGetSize(svd->omega,&Nom));
291 28 : PetscCall(VecGetLocalSize(svd->omega,&nom));
292 28 : PetscCall(MatGetLocalSize(svd->OP,&m,NULL));
293 28 : PetscCheck(Nom==M,PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_SIZ,"Global size of signature (%" PetscInt_FMT ") does not match the row size of A (%" PetscInt_FMT ")",Nom,M);
294 28 : PetscCheck(nom==m,PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_SIZ,"Local size of signature (%" PetscInt_FMT ") does not match the local row size of A (%" PetscInt_FMT ")",nom,m);
295 : }
296 :
297 : /* build transpose matrix */
298 257 : PetscCall(MatDestroy(&svd->A));
299 257 : PetscCall(MatDestroy(&svd->AT));
300 257 : PetscCall(PetscObjectReference((PetscObject)svd->OP));
301 257 : if (svd->expltrans) {
302 222 : if (svd->isgeneralized || M>=N) {
303 187 : svd->A = svd->OP;
304 187 : PetscCall(MatHermitianTranspose(svd->OP,MAT_INITIAL_MATRIX,&svd->AT));
305 : } else {
306 35 : PetscCall(MatHermitianTranspose(svd->OP,MAT_INITIAL_MATRIX,&svd->A));
307 35 : svd->AT = svd->OP;
308 : }
309 : } else {
310 35 : if (svd->isgeneralized || M>=N) {
311 31 : svd->A = svd->OP;
312 31 : PetscCall(MatCreateHermitianTranspose(svd->OP,&svd->AT));
313 : } else {
314 4 : PetscCall(MatCreateHermitianTranspose(svd->OP,&svd->A));
315 4 : svd->AT = svd->OP;
316 : }
317 : }
318 :
319 : /* build transpose matrix B for GSVD */
320 257 : if (svd->isgeneralized) {
321 89 : PetscCall(MatDestroy(&svd->B));
322 89 : PetscCall(MatDestroy(&svd->BT));
323 89 : PetscCall(PetscObjectReference((PetscObject)svd->OPb));
324 89 : if (svd->expltrans) {
325 82 : svd->B = svd->OPb;
326 82 : PetscCall(MatHermitianTranspose(svd->OPb,MAT_INITIAL_MATRIX,&svd->BT));
327 : } else {
328 7 : svd->B = svd->OPb;
329 7 : PetscCall(MatCreateHermitianTranspose(svd->OPb,&svd->BT));
330 : }
331 : }
332 :
333 257 : if (!svd->isgeneralized && M<N) {
334 : /* swap initial vectors */
335 39 : if (svd->nini || svd->ninil) {
336 0 : T=svd->ISL; svd->ISL=svd->IS; svd->IS=T;
337 0 : k=svd->ninil; svd->ninil=svd->nini; svd->nini=k;
338 : }
339 : /* swap basis vectors */
340 39 : if (!svd->swapped) { /* only the first time in case of multiple calls */
341 32 : bv=svd->V; svd->V=svd->U; svd->U=bv;
342 32 : svd->swapped = PETSC_TRUE;
343 : }
344 : }
345 :
346 257 : maxnsol = svd->isgeneralized? PetscMin(PetscMin(M,N),P): PetscMin(M,N);
347 257 : svd->ncv = PetscMin(svd->ncv,maxnsol);
348 257 : svd->nsv = PetscMin(svd->nsv,maxnsol);
349 257 : PetscCheck(svd->ncv==PETSC_DETERMINE || svd->nsv<=svd->ncv,PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_OUTOFRANGE,"nsv bigger than ncv");
350 :
351 : /* relative convergence criterion is not allowed in GSVD */
352 375 : if (svd->conv==(SVDConv)-1) PetscCall(SVDSetConvergenceTest(svd,svd->isgeneralized?SVD_CONV_NORM:SVD_CONV_REL));
353 257 : PetscCheck(!svd->isgeneralized || svd->conv!=SVD_CONV_REL,PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"Relative convergence criterion is not allowed in GSVD");
354 :
355 : /* initialization of matrix norm (standard case only, for GSVD it is done inside setup()) */
356 257 : if (!svd->isgeneralized && svd->conv==SVD_CONV_NORM && !svd->nrma) PetscCall(MatNorm(svd->OP,NORM_INFINITY,&svd->nrma));
357 :
358 : /* threshold stopping test */
359 257 : if (svd->stop==SVD_STOP_THRESHOLD) {
360 12 : PetscCheck(svd->thres>0.0,PetscObjectComm((PetscObject)svd),PETSC_ERR_USER_INPUT,"Must give a threshold value with SVDSetThreshold()");
361 12 : PetscCheck(svd->which==SVD_LARGEST || !svd->threlative,PetscObjectComm((PetscObject)svd),PETSC_ERR_USER_INPUT,"Can only use a relative threshold when which=SVD_LARGEST");
362 12 : PetscCall(PetscNew(&ctx));
363 12 : ctx->thres = svd->thres;
364 12 : ctx->threlative = svd->threlative;
365 12 : ctx->which = svd->which;
366 12 : PetscCall(SVDSetStoppingTestFunction(svd,SVDStoppingThreshold,ctx,PetscCtxDestroyDefault));
367 : }
368 :
369 : /* call specific solver setup */
370 257 : PetscUseTypeMethod(svd,setup);
371 :
372 : /* set tolerance if not yet set */
373 257 : if (svd->tol==(PetscReal)PETSC_DETERMINE) svd->tol = SLEPC_DEFAULT_TOL;
374 :
375 : /* fill sorting criterion context */
376 257 : PetscCall(DSGetSlepcSC(svd->ds,&sc));
377 257 : sc->comparison = (svd->which==SVD_LARGEST)? SlepcCompareLargestReal: SlepcCompareSmallestReal;
378 257 : sc->comparisonctx = NULL;
379 257 : sc->map = NULL;
380 257 : sc->mapobj = NULL;
381 :
382 : /* process initial vectors */
383 257 : if (svd->nini<0) {
384 25 : k = -svd->nini;
385 25 : PetscCheck(k<=svd->ncv,PetscObjectComm((PetscObject)svd),PETSC_ERR_USER_INPUT,"The number of initial vectors is larger than ncv");
386 25 : PetscCall(BVInsertVecs(svd->V,0,&k,svd->IS,PETSC_TRUE));
387 25 : PetscCall(SlepcBasisDestroy_Private(&svd->nini,&svd->IS));
388 25 : svd->nini = k;
389 : }
390 257 : if (svd->ninil<0) {
391 29 : k = 0;
392 29 : if (svd->leftbasis) {
393 23 : k = -svd->ninil;
394 23 : PetscCheck(k<=svd->ncv,PetscObjectComm((PetscObject)svd),PETSC_ERR_USER_INPUT,"The number of left initial vectors is larger than ncv");
395 23 : PetscCall(BVInsertVecs(svd->U,0,&k,svd->ISL,PETSC_TRUE));
396 6 : } else PetscCall(PetscInfo(svd,"Ignoring initial left vectors\n"));
397 29 : PetscCall(SlepcBasisDestroy_Private(&svd->ninil,&svd->ISL));
398 29 : svd->ninil = k;
399 : }
400 :
401 257 : PetscCall(PetscLogEventEnd(SVD_SetUp,svd,0,0,0));
402 257 : svd->state = SVD_STATE_SETUP;
403 257 : PetscFunctionReturn(PETSC_SUCCESS);
404 : }
405 :
406 : /*@
407 : SVDSetInitialSpaces - Specify two basis of vectors that constitute the initial
408 : right and/or left spaces.
409 :
410 : Collective
411 :
412 : Input Parameters:
413 : + svd - the singular value solver context
414 : . nr - number of right vectors
415 : . isr - set of basis vectors of the right initial space
416 : . nl - number of left vectors
417 : - isl - set of basis vectors of the left initial space
418 :
419 : Notes:
420 : The initial right and left spaces are rough approximations to the right and/or
421 : left singular subspaces from which the solver starts to iterate.
422 : It is not necessary to provide both sets of vectors.
423 :
424 : Some solvers start to iterate on a single vector (initial vector). In that case,
425 : the other vectors are ignored.
426 :
427 : These vectors do not persist from one SVDSolve() call to the other, so the
428 : initial space should be set every time.
429 :
430 : The vectors do not need to be mutually orthonormal, since they are explicitly
431 : orthonormalized internally.
432 :
433 : Common usage of this function is when the user can provide a rough approximation
434 : of the wanted singular space. Then, convergence may be faster.
435 :
436 : Level: intermediate
437 :
438 : .seealso: SVDSetUp()
439 : @*/
440 33 : PetscErrorCode SVDSetInitialSpaces(SVD svd,PetscInt nr,Vec isr[],PetscInt nl,Vec isl[])
441 : {
442 33 : PetscFunctionBegin;
443 33 : PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
444 99 : PetscValidLogicalCollectiveInt(svd,nr,2);
445 99 : PetscValidLogicalCollectiveInt(svd,nl,4);
446 33 : PetscCheck(nr>=0,PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_OUTOFRANGE,"Argument nr cannot be negative");
447 33 : if (nr>0) {
448 33 : PetscAssertPointer(isr,3);
449 33 : PetscValidHeaderSpecific(*isr,VEC_CLASSID,3);
450 : }
451 33 : PetscCheck(nl>=0,PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_OUTOFRANGE,"Argument nl cannot be negative");
452 33 : if (nl>0) {
453 33 : PetscAssertPointer(isl,5);
454 33 : PetscValidHeaderSpecific(*isl,VEC_CLASSID,5);
455 : }
456 33 : PetscCall(SlepcBasisReference_Private(nr,isr,&svd->nini,&svd->IS));
457 33 : PetscCall(SlepcBasisReference_Private(nl,isl,&svd->ninil,&svd->ISL));
458 33 : if (nr>0 || nl>0) svd->state = SVD_STATE_INITIAL;
459 33 : PetscFunctionReturn(PETSC_SUCCESS);
460 : }
461 :
462 : /*
463 : SVDSetDimensions_Default - Set reasonable values for ncv, mpd if not set
464 : by the user. This is called at setup.
465 : */
466 139 : PetscErrorCode SVDSetDimensions_Default(SVD svd)
467 : {
468 139 : PetscInt N,M,P,maxnsol;
469 :
470 139 : PetscFunctionBegin;
471 139 : PetscCall(MatGetSize(svd->OP,&M,&N));
472 139 : maxnsol = PetscMin(M,N);
473 139 : if (svd->isgeneralized) {
474 59 : PetscCall(MatGetSize(svd->OPb,&P,NULL));
475 59 : maxnsol = PetscMin(maxnsol,P);
476 : }
477 139 : if (svd->nsv==0 && svd->stop!=SVD_STOP_THRESHOLD) svd->nsv = 1;
478 139 : if (svd->ncv!=PETSC_DETERMINE) { /* ncv set */
479 59 : PetscCheck(svd->ncv>=svd->nsv,PetscObjectComm((PetscObject)svd),PETSC_ERR_USER_INPUT,"The value of ncv must be at least nsv");
480 80 : } else if (svd->mpd!=PETSC_DETERMINE) { /* mpd set */
481 2 : svd->ncv = PetscMin(maxnsol,svd->nsv+svd->mpd);
482 : } else { /* neither set: defaults depend on nsv being small or large */
483 78 : if (svd->nsv<500) svd->ncv = PetscMin(maxnsol,PetscMax(2*svd->nsv,10));
484 : else {
485 0 : svd->mpd = 500;
486 0 : svd->ncv = PetscMin(maxnsol,svd->nsv+svd->mpd);
487 : }
488 : }
489 139 : if (svd->mpd==PETSC_DETERMINE) svd->mpd = svd->ncv;
490 139 : PetscFunctionReturn(PETSC_SUCCESS);
491 : }
492 :
493 : /*@
494 : SVDAllocateSolution - Allocate memory storage for common variables such
495 : as the singular values and the basis vectors.
496 :
497 : Collective
498 :
499 : Input Parameters:
500 : + svd - singular value solver context
501 : - extra - number of additional positions, used for methods that require a
502 : working basis slightly larger than ncv
503 :
504 : Developer Notes:
505 : This is SLEPC_EXTERN because it may be required by user plugin SVD
506 : implementations.
507 :
508 : This is called at setup after setting the value of ncv and the flag leftbasis.
509 :
510 : Level: developer
511 :
512 : .seealso: SVDSetUp()
513 : @*/
514 257 : PetscErrorCode SVDAllocateSolution(SVD svd,PetscInt extra)
515 : {
516 257 : PetscInt oldsize,requested;
517 257 : Vec tr,tl;
518 :
519 257 : PetscFunctionBegin;
520 257 : requested = svd->ncv + extra;
521 :
522 : /* oldsize is zero if this is the first time setup is called */
523 257 : PetscCall(BVGetSizes(svd->V,NULL,NULL,&oldsize));
524 :
525 : /* allocate sigma */
526 257 : if (requested != oldsize || !svd->sigma) {
527 222 : PetscCall(PetscFree3(svd->sigma,svd->perm,svd->errest));
528 222 : if (svd->sign) PetscCall(PetscFree(svd->sign));
529 222 : PetscCall(PetscMalloc3(requested,&svd->sigma,requested,&svd->perm,requested,&svd->errest));
530 222 : if (svd->ishyperbolic) PetscCall(PetscMalloc1(requested,&svd->sign));
531 : }
532 : /* allocate V */
533 257 : if (!svd->V) PetscCall(SVDGetBV(svd,&svd->V,NULL));
534 257 : if (!oldsize) {
535 213 : if (!((PetscObject)svd->V)->type_name) PetscCall(BVSetType(svd->V,BVMAT));
536 213 : PetscCall(MatCreateVecsEmpty(svd->A,&tr,NULL));
537 213 : PetscCall(BVSetSizesFromVec(svd->V,tr,requested));
538 213 : PetscCall(VecDestroy(&tr));
539 44 : } else PetscCall(BVResize(svd->V,requested,PETSC_FALSE));
540 : /* allocate U */
541 257 : if (svd->leftbasis && !svd->isgeneralized) {
542 127 : if (!svd->U) PetscCall(SVDGetBV(svd,NULL,&svd->U));
543 127 : if (!oldsize) {
544 102 : if (!((PetscObject)svd->U)->type_name) PetscCall(BVSetType(svd->U,((PetscObject)svd->V)->type_name));
545 102 : PetscCall(MatCreateVecsEmpty(svd->A,NULL,&tl));
546 102 : PetscCall(BVSetSizesFromVec(svd->U,tl,requested));
547 102 : PetscCall(VecDestroy(&tl));
548 25 : } else PetscCall(BVResize(svd->U,requested,PETSC_FALSE));
549 130 : } else if (svd->isgeneralized) { /* left basis for the GSVD */
550 89 : if (!svd->U) PetscCall(SVDGetBV(svd,NULL,&svd->U));
551 89 : if (!oldsize) {
552 78 : if (!((PetscObject)svd->U)->type_name) PetscCall(BVSetType(svd->U,((PetscObject)svd->V)->type_name));
553 78 : PetscCall(SVDCreateLeftTemplate(svd,&tl));
554 78 : PetscCall(BVSetSizesFromVec(svd->U,tl,requested));
555 78 : PetscCall(VecDestroy(&tl));
556 11 : } else PetscCall(BVResize(svd->U,requested,PETSC_FALSE));
557 : }
558 257 : PetscFunctionReturn(PETSC_SUCCESS);
559 : }
560 :
561 : /*@
562 : SVDReallocateSolution - Reallocate memory storage for common variables such
563 : as the singular values and the basis vectors.
564 :
565 : Collective
566 :
567 : Input Parameters:
568 : + svd - singular value solver context
569 : - newsize - new size
570 :
571 : Developer Notes:
572 : This is SLEPC_EXTERN because it may be required by user plugin SVD
573 : implementations.
574 :
575 : This is called during the iteration in case the threshold stopping test has
576 : been selected.
577 :
578 : Level: developer
579 :
580 : .seealso: SVDAllocateSolution(), SVDSetThreshold()
581 : @*/
582 4 : PetscErrorCode SVDReallocateSolution(SVD svd,PetscInt newsize)
583 : {
584 4 : PetscInt oldsize,*nperm;
585 4 : PetscReal *nsigma,*nerrest,*nsign;
586 :
587 4 : PetscFunctionBegin;
588 4 : PetscCall(BVGetSizes(svd->V,NULL,NULL,&oldsize));
589 4 : if (oldsize>=newsize) PetscFunctionReturn(PETSC_SUCCESS);
590 4 : PetscCall(PetscInfo(svd,"Reallocating basis vectors to %" PetscInt_FMT " columns\n",newsize));
591 :
592 : /* reallocate sigma */
593 4 : PetscCall(PetscMalloc3(newsize,&nsigma,newsize,&nperm,newsize,&nerrest));
594 4 : PetscCall(PetscArraycpy(nsigma,svd->sigma,oldsize));
595 4 : PetscCall(PetscArraycpy(nperm,svd->perm,oldsize));
596 4 : PetscCall(PetscArraycpy(nerrest,svd->errest,oldsize));
597 4 : PetscCall(PetscFree3(svd->sigma,svd->perm,svd->errest));
598 4 : svd->sigma = nsigma;
599 4 : svd->perm = nperm;
600 4 : svd->errest = nerrest;
601 4 : if (svd->ishyperbolic) {
602 1 : PetscCall(PetscMalloc1(newsize,&nsign));
603 1 : PetscCall(PetscArraycpy(nsign,svd->sign,oldsize));
604 1 : PetscCall(PetscFree(svd->sign));
605 1 : svd->sign = nsign;
606 : }
607 : /* reallocate V,U */
608 4 : PetscCall(BVResize(svd->V,newsize,PETSC_TRUE));
609 4 : if (svd->leftbasis || svd->isgeneralized) PetscCall(BVResize(svd->U,newsize,PETSC_TRUE));
610 4 : PetscFunctionReturn(PETSC_SUCCESS);
611 : }
|
