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 : ST interface routines related to the KSP object associated with it
12 : */
13 :
14 : #include <slepc/private/stimpl.h> /*I "slepcst.h" I*/
15 :
16 : /*
17 : This is used to set a default type for the KSP and PC objects.
18 : It is called at STSetFromOptions (before KSPSetFromOptions)
19 : and also at STSetUp (in case STSetFromOptions was not called).
20 : */
21 1141 : PetscErrorCode STSetDefaultKSP(ST st)
22 : {
23 1141 : PetscFunctionBegin;
24 1141 : PetscValidHeaderSpecific(st,ST_CLASSID,1);
25 1141 : PetscValidType(st,1);
26 1141 : if (!st->ksp) PetscCall(STGetKSP(st,&st->ksp));
27 1141 : PetscTryTypeMethod(st,setdefaultksp);
28 1141 : PetscFunctionReturn(PETSC_SUCCESS);
29 : }
30 :
31 : /*
32 : This is done by all ST types except PRECOND.
33 : The default is an LU direct solver, or GMRES+Jacobi if matmode=shell.
34 : */
35 850 : PetscErrorCode STSetDefaultKSP_Default(ST st)
36 : {
37 850 : PC pc;
38 850 : PCType pctype;
39 850 : KSPType ksptype;
40 :
41 850 : PetscFunctionBegin;
42 850 : PetscCall(KSPGetPC(st->ksp,&pc));
43 850 : PetscCall(KSPGetType(st->ksp,&ksptype));
44 850 : PetscCall(PCGetType(pc,&pctype));
45 850 : if (!pctype && !ksptype) {
46 592 : if (st->Pmat || st->Psplit) {
47 21 : PetscCall(KSPSetType(st->ksp,KSPBCGS));
48 21 : PetscCall(PCSetType(pc,PCBJACOBI));
49 571 : } else if (st->matmode == ST_MATMODE_SHELL) {
50 49 : PetscCall(KSPSetType(st->ksp,KSPGMRES));
51 49 : PetscCall(PCSetType(pc,PCJACOBI));
52 : } else {
53 522 : PetscCall(KSPSetType(st->ksp,KSPPREONLY));
54 1014 : PetscCall(PCSetType(pc,st->asymm?PCCHOLESKY:PCLU));
55 : }
56 : }
57 850 : PetscCall(KSPSetErrorIfNotConverged(st->ksp,PETSC_TRUE));
58 850 : PetscFunctionReturn(PETSC_SUCCESS);
59 : }
60 :
61 : /*@
62 : STMatMult - Computes the matrix-vector product y = T[k] x, where T[k] is
63 : the k-th matrix of the spectral transformation.
64 :
65 : Neighbor-wise Collective
66 :
67 : Input Parameters:
68 : + st - the spectral transformation context
69 : . k - index of matrix to use
70 : - x - the vector to be multiplied
71 :
72 : Output Parameter:
73 : . y - the result
74 :
75 : Level: developer
76 :
77 : .seealso: STMatMultTranspose(), STMatMultHermitianTranspose()
78 : @*/
79 64912 : PetscErrorCode STMatMult(ST st,PetscInt k,Vec x,Vec y)
80 : {
81 64912 : PetscFunctionBegin;
82 64912 : PetscValidHeaderSpecific(st,ST_CLASSID,1);
83 194736 : PetscValidLogicalCollectiveInt(st,k,2);
84 64912 : PetscValidHeaderSpecific(x,VEC_CLASSID,3);
85 64912 : PetscValidHeaderSpecific(y,VEC_CLASSID,4);
86 64912 : STCheckMatrices(st,1);
87 64912 : PetscCheck(k>=0 && k<PetscMax(2,st->nmat),PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_OUTOFRANGE,"k must be between 0 and %" PetscInt_FMT,st->nmat);
88 64912 : PetscCheck(x!=y,PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
89 64912 : PetscCall(VecSetErrorIfLocked(y,3));
90 :
91 64912 : if (st->state!=ST_STATE_SETUP) PetscCall(STSetUp(st));
92 64912 : PetscCall(VecLockReadPush(x));
93 64912 : PetscCall(PetscLogEventBegin(ST_MatMult,st,x,y,0));
94 64912 : if (!st->T[k]) PetscCall(VecCopy(x,y)); /* T[k]=NULL means identity matrix */
95 64912 : else PetscCall(MatMult(st->T[k],x,y));
96 64912 : PetscCall(PetscLogEventEnd(ST_MatMult,st,x,y,0));
97 64912 : PetscCall(VecLockReadPop(x));
98 64912 : PetscFunctionReturn(PETSC_SUCCESS);
99 : }
100 :
101 : /*@
102 : STMatMultTranspose - Computes the matrix-vector product y = T[k]^T x, where T[k] is
103 : the k-th matrix of the spectral transformation.
104 :
105 : Neighbor-wise Collective
106 :
107 : Input Parameters:
108 : + st - the spectral transformation context
109 : . k - index of matrix to use
110 : - x - the vector to be multiplied
111 :
112 : Output Parameter:
113 : . y - the result
114 :
115 : Level: developer
116 :
117 : .seealso: STMatMult(), STMatMultHermitianTranspose()
118 : @*/
119 6 : PetscErrorCode STMatMultTranspose(ST st,PetscInt k,Vec x,Vec y)
120 : {
121 6 : PetscFunctionBegin;
122 6 : PetscValidHeaderSpecific(st,ST_CLASSID,1);
123 18 : PetscValidLogicalCollectiveInt(st,k,2);
124 6 : PetscValidHeaderSpecific(x,VEC_CLASSID,3);
125 6 : PetscValidHeaderSpecific(y,VEC_CLASSID,4);
126 6 : STCheckMatrices(st,1);
127 6 : PetscCheck(k>=0 && k<PetscMax(2,st->nmat),PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_OUTOFRANGE,"k must be between 0 and %" PetscInt_FMT,st->nmat);
128 6 : PetscCheck(x!=y,PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
129 6 : PetscCall(VecSetErrorIfLocked(y,3));
130 :
131 6 : if (st->state!=ST_STATE_SETUP) PetscCall(STSetUp(st));
132 6 : PetscCall(VecLockReadPush(x));
133 6 : PetscCall(PetscLogEventBegin(ST_MatMultTranspose,st,x,y,0));
134 6 : if (!st->T[k]) PetscCall(VecCopy(x,y)); /* T[k]=NULL means identity matrix */
135 6 : else PetscCall(MatMultTranspose(st->T[k],x,y));
136 6 : PetscCall(PetscLogEventEnd(ST_MatMultTranspose,st,x,y,0));
137 6 : PetscCall(VecLockReadPop(x));
138 6 : PetscFunctionReturn(PETSC_SUCCESS);
139 : }
140 :
141 : /*@
142 : STMatMultHermitianTranspose - Computes the matrix-vector product y = T[k]^H x, where T[k] is
143 : the k-th matrix of the spectral transformation.
144 :
145 : Neighbor-wise Collective
146 :
147 : Input Parameters:
148 : + st - the spectral transformation context
149 : . k - index of matrix to use
150 : - x - the vector to be multiplied
151 :
152 : Output Parameter:
153 : . y - the result
154 :
155 : Level: developer
156 :
157 : .seealso: STMatMult(), STMatMultTranspose()
158 : @*/
159 3 : PetscErrorCode STMatMultHermitianTranspose(ST st,PetscInt k,Vec x,Vec y)
160 : {
161 3 : PetscFunctionBegin;
162 3 : PetscValidHeaderSpecific(st,ST_CLASSID,1);
163 9 : PetscValidLogicalCollectiveInt(st,k,2);
164 3 : PetscValidHeaderSpecific(x,VEC_CLASSID,3);
165 3 : PetscValidHeaderSpecific(y,VEC_CLASSID,4);
166 3 : STCheckMatrices(st,1);
167 3 : PetscCheck(k>=0 && k<PetscMax(2,st->nmat),PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_OUTOFRANGE,"k must be between 0 and %" PetscInt_FMT,st->nmat);
168 3 : PetscCheck(x!=y,PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
169 3 : PetscCall(VecSetErrorIfLocked(y,3));
170 :
171 3 : if (st->state!=ST_STATE_SETUP) PetscCall(STSetUp(st));
172 3 : PetscCall(VecLockReadPush(x));
173 3 : PetscCall(PetscLogEventBegin(ST_MatMultTranspose,st,x,y,0));
174 3 : if (!st->T[k]) PetscCall(VecCopy(x,y)); /* T[k]=NULL means identity matrix */
175 3 : else PetscCall(MatMultHermitianTranspose(st->T[k],x,y));
176 3 : PetscCall(PetscLogEventEnd(ST_MatMultTranspose,st,x,y,0));
177 3 : PetscCall(VecLockReadPop(x));
178 3 : PetscFunctionReturn(PETSC_SUCCESS);
179 : }
180 :
181 : /*@
182 : STMatSolve - Solves P x = b, where P is the preconditioner matrix of
183 : the spectral transformation, using a KSP object stored internally.
184 :
185 : Collective
186 :
187 : Input Parameters:
188 : + st - the spectral transformation context
189 : - b - right hand side vector
190 :
191 : Output Parameter:
192 : . x - computed solution
193 :
194 : Level: developer
195 :
196 : .seealso: STMatSolveTranspose(), STMatSolveHermitianTranspose(), STMatMatSolve()
197 : @*/
198 53183 : PetscErrorCode STMatSolve(ST st,Vec b,Vec x)
199 : {
200 53183 : PetscFunctionBegin;
201 53183 : PetscValidHeaderSpecific(st,ST_CLASSID,1);
202 53183 : PetscValidHeaderSpecific(b,VEC_CLASSID,2);
203 53183 : PetscValidHeaderSpecific(x,VEC_CLASSID,3);
204 53183 : STCheckMatrices(st,1);
205 53183 : PetscCheck(x!=b,PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_IDN,"x and b must be different vectors");
206 53183 : PetscCall(VecSetErrorIfLocked(x,3));
207 :
208 53183 : if (st->state!=ST_STATE_SETUP) PetscCall(STSetUp(st));
209 53183 : PetscCall(VecLockReadPush(b));
210 53183 : PetscCall(PetscLogEventBegin(ST_MatSolve,st,b,x,0));
211 53183 : if (!st->P) PetscCall(VecCopy(b,x)); /* P=NULL means identity matrix */
212 53175 : else PetscCall(KSPSolve(st->ksp,b,x));
213 53183 : PetscCall(PetscLogEventEnd(ST_MatSolve,st,b,x,0));
214 53183 : PetscCall(VecLockReadPop(b));
215 53183 : PetscFunctionReturn(PETSC_SUCCESS);
216 : }
217 :
218 : /*@
219 : STMatMatSolve - Solves P X = B, where P is the preconditioner matrix of
220 : the spectral transformation, using a KSP object stored internally.
221 :
222 : Collective
223 :
224 : Input Parameters:
225 : + st - the spectral transformation context
226 : - B - right hand side vectors
227 :
228 : Output Parameter:
229 : . X - computed solutions
230 :
231 : Level: developer
232 :
233 : .seealso: STMatSolve()
234 : @*/
235 3362 : PetscErrorCode STMatMatSolve(ST st,Mat B,Mat X)
236 : {
237 3362 : PetscFunctionBegin;
238 3362 : PetscValidHeaderSpecific(st,ST_CLASSID,1);
239 3362 : PetscValidHeaderSpecific(B,MAT_CLASSID,2);
240 3362 : PetscValidHeaderSpecific(X,MAT_CLASSID,3);
241 3362 : STCheckMatrices(st,1);
242 :
243 3362 : if (st->state!=ST_STATE_SETUP) PetscCall(STSetUp(st));
244 3362 : PetscCall(PetscLogEventBegin(ST_MatSolve,st,B,X,0));
245 3362 : if (!st->P) PetscCall(MatCopy(B,X,SAME_NONZERO_PATTERN)); /* P=NULL means identity matrix */
246 3362 : else PetscCall(KSPMatSolve(st->ksp,B,X));
247 3362 : PetscCall(PetscLogEventEnd(ST_MatSolve,st,B,X,0));
248 3362 : PetscFunctionReturn(PETSC_SUCCESS);
249 : }
250 :
251 : /*@
252 : STMatSolveTranspose - Solves P^T x = b, where P is the preconditioner matrix of
253 : the spectral transformation, using a KSP object stored internally.
254 :
255 : Collective
256 :
257 : Input Parameters:
258 : + st - the spectral transformation context
259 : - b - right hand side vector
260 :
261 : Output Parameter:
262 : . x - computed solution
263 :
264 : Level: developer
265 :
266 : .seealso: STMatSolve()
267 : @*/
268 18 : PetscErrorCode STMatSolveTranspose(ST st,Vec b,Vec x)
269 : {
270 18 : PetscFunctionBegin;
271 18 : PetscValidHeaderSpecific(st,ST_CLASSID,1);
272 18 : PetscValidHeaderSpecific(b,VEC_CLASSID,2);
273 18 : PetscValidHeaderSpecific(x,VEC_CLASSID,3);
274 18 : STCheckMatrices(st,1);
275 18 : PetscCheck(x!=b,PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_IDN,"x and b must be different vectors");
276 18 : PetscCall(VecSetErrorIfLocked(x,3));
277 :
278 18 : if (st->state!=ST_STATE_SETUP) PetscCall(STSetUp(st));
279 18 : PetscCall(VecLockReadPush(b));
280 18 : PetscCall(PetscLogEventBegin(ST_MatSolveTranspose,st,b,x,0));
281 18 : if (!st->P) PetscCall(VecCopy(b,x)); /* P=NULL means identity matrix */
282 18 : else PetscCall(KSPSolveTranspose(st->ksp,b,x));
283 18 : PetscCall(PetscLogEventEnd(ST_MatSolveTranspose,st,b,x,0));
284 18 : PetscCall(VecLockReadPop(b));
285 18 : PetscFunctionReturn(PETSC_SUCCESS);
286 : }
287 :
288 : /*@
289 : STMatSolveHermitianTranspose - Solves P^H x = b, where P is the preconditioner matrix of
290 : the spectral transformation, using a KSP object stored internally.
291 :
292 : Collective
293 :
294 : Input Parameters:
295 : + st - the spectral transformation context
296 : - b - right hand side vector
297 :
298 : Output Parameter:
299 : . x - computed solution
300 :
301 : Level: developer
302 :
303 : .seealso: STMatSolve()
304 : @*/
305 564 : PetscErrorCode STMatSolveHermitianTranspose(ST st,Vec b,Vec x)
306 : {
307 564 : Vec w;
308 :
309 564 : PetscFunctionBegin;
310 564 : PetscValidHeaderSpecific(st,ST_CLASSID,1);
311 564 : PetscValidHeaderSpecific(b,VEC_CLASSID,2);
312 564 : PetscValidHeaderSpecific(x,VEC_CLASSID,3);
313 564 : STCheckMatrices(st,1);
314 564 : PetscCheck(x!=b,PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_IDN,"x and b must be different vectors");
315 564 : PetscCall(VecSetErrorIfLocked(x,3));
316 :
317 564 : if (st->state!=ST_STATE_SETUP) PetscCall(STSetUp(st));
318 564 : PetscCall(VecLockReadPush(b));
319 564 : PetscCall(PetscLogEventBegin(ST_MatSolveTranspose,st,b,x,0));
320 564 : if (!st->P) PetscCall(VecCopy(b,x)); /* P=NULL means identity matrix */
321 : else {
322 564 : PetscCall(VecDuplicate(b,&w));
323 564 : PetscCall(VecCopy(b,w));
324 564 : PetscCall(VecConjugate(w));
325 564 : PetscCall(KSPSolveTranspose(st->ksp,w,x));
326 564 : PetscCall(VecDestroy(&w));
327 564 : PetscCall(VecConjugate(x));
328 : }
329 564 : PetscCall(PetscLogEventEnd(ST_MatSolveTranspose,st,b,x,0));
330 564 : PetscCall(VecLockReadPop(b));
331 564 : PetscFunctionReturn(PETSC_SUCCESS);
332 : }
333 :
334 1333 : PetscErrorCode STCheckFactorPackage(ST st)
335 : {
336 1333 : PC pc;
337 1333 : PetscMPIInt size;
338 1333 : PetscBool flg;
339 1333 : MatSolverType stype;
340 :
341 1333 : PetscFunctionBegin;
342 1333 : PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)st),&size));
343 1333 : if (size==1) PetscFunctionReturn(PETSC_SUCCESS);
344 82 : PetscCall(KSPGetPC(st->ksp,&pc));
345 82 : PetscCall(PCFactorGetMatSolverType(pc,&stype));
346 82 : if (stype) { /* currently selected PC is a factorization */
347 0 : PetscCall(PetscStrcmp(stype,MATSOLVERPETSC,&flg));
348 0 : PetscCheck(!flg,PetscObjectComm((PetscObject)st),PETSC_ERR_SUP,"You chose to solve linear systems with a factorization, but in parallel runs you need to select an external package; see the users guide for details");
349 : }
350 82 : PetscFunctionReturn(PETSC_SUCCESS);
351 : }
352 :
353 : /*@
354 : STSetKSP - Sets the KSP object associated with the spectral
355 : transformation.
356 :
357 : Collective
358 :
359 : Input Parameters:
360 : + st - the spectral transformation context
361 : - ksp - the linear system context
362 :
363 : Level: advanced
364 :
365 : .seealso: STGetKSP()
366 : @*/
367 6 : PetscErrorCode STSetKSP(ST st,KSP ksp)
368 : {
369 6 : PetscFunctionBegin;
370 6 : PetscValidHeaderSpecific(st,ST_CLASSID,1);
371 6 : PetscValidHeaderSpecific(ksp,KSP_CLASSID,2);
372 6 : PetscCheckSameComm(st,1,ksp,2);
373 6 : STCheckNotSeized(st,1);
374 6 : PetscCall(PetscObjectReference((PetscObject)ksp));
375 6 : PetscCall(KSPDestroy(&st->ksp));
376 6 : st->ksp = ksp;
377 6 : PetscFunctionReturn(PETSC_SUCCESS);
378 : }
379 :
380 : /*@
381 : STGetKSP - Gets the KSP object associated with the spectral
382 : transformation.
383 :
384 : Collective
385 :
386 : Input Parameter:
387 : . st - the spectral transformation context
388 :
389 : Output Parameter:
390 : . ksp - the linear system context
391 :
392 : Level: intermediate
393 :
394 : .seealso: STSetKSP()
395 : @*/
396 2708 : PetscErrorCode STGetKSP(ST st,KSP* ksp)
397 : {
398 2708 : PetscFunctionBegin;
399 2708 : PetscValidHeaderSpecific(st,ST_CLASSID,1);
400 2708 : PetscAssertPointer(ksp,2);
401 2708 : if (!st->ksp) {
402 783 : PetscCall(KSPCreate(PetscObjectComm((PetscObject)st),&st->ksp));
403 783 : PetscCall(PetscObjectIncrementTabLevel((PetscObject)st->ksp,(PetscObject)st,1));
404 783 : PetscCall(KSPSetOptionsPrefix(st->ksp,((PetscObject)st)->prefix));
405 783 : PetscCall(KSPAppendOptionsPrefix(st->ksp,"st_"));
406 783 : PetscCall(PetscObjectSetOptions((PetscObject)st->ksp,((PetscObject)st)->options));
407 783 : PetscCall(KSPSetTolerances(st->ksp,SLEPC_DEFAULT_TOL,PETSC_CURRENT,PETSC_CURRENT,PETSC_CURRENT));
408 : }
409 2708 : *ksp = st->ksp;
410 2708 : PetscFunctionReturn(PETSC_SUCCESS);
411 : }
412 :
413 1 : PetscErrorCode STCheckNullSpace_Default(ST st,BV V)
414 : {
415 1 : PetscInt nc,i,c;
416 1 : PetscReal norm;
417 1 : Vec *T,w,vi;
418 1 : Mat A;
419 1 : PC pc;
420 1 : MatNullSpace nullsp;
421 :
422 1 : PetscFunctionBegin;
423 1 : PetscCall(BVGetNumConstraints(V,&nc));
424 1 : PetscCall(PetscMalloc1(nc,&T));
425 1 : if (!st->ksp) PetscCall(STGetKSP(st,&st->ksp));
426 1 : PetscCall(KSPGetPC(st->ksp,&pc));
427 1 : PetscCall(PCGetOperators(pc,&A,NULL));
428 1 : PetscCall(MatCreateVecs(A,NULL,&w));
429 : c = 0;
430 2 : for (i=0;i<nc;i++) {
431 1 : PetscCall(BVGetColumn(V,-nc+i,&vi));
432 1 : PetscCall(MatMult(A,vi,w));
433 1 : PetscCall(VecNorm(w,NORM_2,&norm));
434 1 : if (norm < 10.0*PETSC_SQRT_MACHINE_EPSILON) {
435 1 : PetscCall(PetscInfo(st,"Vector %" PetscInt_FMT " included in the nullspace of OP, norm=%g\n",i,(double)norm));
436 1 : PetscCall(BVCreateVec(V,T+c));
437 1 : PetscCall(VecCopy(vi,T[c]));
438 1 : PetscCall(VecNormalize(T[c],NULL));
439 1 : c++;
440 0 : } else PetscCall(PetscInfo(st,"Vector %" PetscInt_FMT " discarded as possible nullspace of OP, norm=%g\n",i,(double)norm));
441 1 : PetscCall(BVRestoreColumn(V,-nc+i,&vi));
442 : }
443 1 : PetscCall(VecDestroy(&w));
444 1 : if (c>0) {
445 1 : PetscCall(MatNullSpaceCreate(PetscObjectComm((PetscObject)st),PETSC_FALSE,c,T,&nullsp));
446 1 : PetscCall(MatSetNullSpace(A,nullsp));
447 1 : PetscCall(MatNullSpaceDestroy(&nullsp));
448 1 : PetscCall(VecDestroyVecs(c,&T));
449 0 : } else PetscCall(PetscFree(T));
450 1 : PetscFunctionReturn(PETSC_SUCCESS);
451 : }
452 :
453 : /*@
454 : STCheckNullSpace - Tests if constraint vectors are nullspace vectors.
455 :
456 : Collective
457 :
458 : Input Parameters:
459 : + st - the spectral transformation context
460 : - V - basis vectors to be checked
461 :
462 : Notes:
463 : Given a basis vectors object, this function tests each of its constraint
464 : vectors to be a nullspace vector of the coefficient matrix of the
465 : associated KSP object. All these nullspace vectors are passed to the KSP
466 : object.
467 :
468 : This function allows handling singular pencils and solving some problems
469 : in which the nullspace is important (see the users guide for details).
470 :
471 : Level: developer
472 :
473 : .seealso: EPSSetDeflationSpace()
474 : @*/
475 14 : PetscErrorCode STCheckNullSpace(ST st,BV V)
476 : {
477 14 : PetscInt nc;
478 :
479 14 : PetscFunctionBegin;
480 14 : PetscValidHeaderSpecific(st,ST_CLASSID,1);
481 14 : PetscValidHeaderSpecific(V,BV_CLASSID,2);
482 14 : PetscValidType(st,1);
483 14 : PetscCheckSameComm(st,1,V,2);
484 14 : PetscCheck(st->state,PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_WRONGSTATE,"Must call STSetUp() first");
485 :
486 14 : PetscCall(BVGetNumConstraints(V,&nc));
487 14 : if (nc) PetscTryTypeMethod(st,checknullspace,V);
488 14 : PetscFunctionReturn(PETSC_SUCCESS);
489 : }
|
