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 : DS operations: DSSolve(), DSVectors(), etc
12 : */
13 :
14 : #include <slepc/private/dsimpl.h> /*I "slepcds.h" I*/
15 :
16 : /*@
17 : DSGetLeadingDimension - Returns the leading dimension of the allocated
18 : matrices.
19 :
20 : Not Collective
21 :
22 : Input Parameter:
23 : . ds - the direct solver context
24 :
25 : Output Parameter:
26 : . ld - leading dimension (maximum allowed dimension for the matrices)
27 :
28 : Level: advanced
29 :
30 : .seealso: DSAllocate(), DSSetDimensions()
31 : @*/
32 37175 : PetscErrorCode DSGetLeadingDimension(DS ds,PetscInt *ld)
33 : {
34 37175 : PetscFunctionBegin;
35 37175 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
36 37175 : PetscAssertPointer(ld,2);
37 37175 : *ld = ds->ld;
38 37175 : PetscFunctionReturn(PETSC_SUCCESS);
39 : }
40 :
41 : /*@
42 : DSSetState - Change the state of the DS object.
43 :
44 : Logically Collective
45 :
46 : Input Parameters:
47 : + ds - the direct solver context
48 : - state - the new state
49 :
50 : Notes:
51 : The state indicates that the dense system is in an initial state (raw),
52 : in an intermediate state (such as tridiagonal, Hessenberg or
53 : Hessenberg-triangular), in a condensed state (such as diagonal, Schur or
54 : generalized Schur), or in a truncated state.
55 :
56 : This function is normally used to return to the raw state when the
57 : condensed structure is destroyed.
58 :
59 : Level: advanced
60 :
61 : .seealso: DSGetState()
62 : @*/
63 22434 : PetscErrorCode DSSetState(DS ds,DSStateType state)
64 : {
65 22434 : PetscFunctionBegin;
66 22434 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
67 89736 : PetscValidLogicalCollectiveEnum(ds,state,2);
68 22434 : switch (state) {
69 22434 : case DS_STATE_RAW:
70 : case DS_STATE_INTERMEDIATE:
71 : case DS_STATE_CONDENSED:
72 : case DS_STATE_TRUNCATED:
73 22434 : if (ds->state!=state) PetscCall(PetscInfo(ds,"State has changed from %s to %s\n",DSStateTypes[ds->state],DSStateTypes[state]));
74 22434 : ds->state = state;
75 22434 : break;
76 0 : default:
77 0 : SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONG,"Wrong state");
78 : }
79 22434 : PetscFunctionReturn(PETSC_SUCCESS);
80 : }
81 :
82 : /*@
83 : DSGetState - Returns the current state.
84 :
85 : Not Collective
86 :
87 : Input Parameter:
88 : . ds - the direct solver context
89 :
90 : Output Parameter:
91 : . state - current state
92 :
93 : Level: advanced
94 :
95 : .seealso: DSSetState()
96 : @*/
97 2 : PetscErrorCode DSGetState(DS ds,DSStateType *state)
98 : {
99 2 : PetscFunctionBegin;
100 2 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
101 2 : PetscAssertPointer(state,2);
102 2 : *state = ds->state;
103 2 : PetscFunctionReturn(PETSC_SUCCESS);
104 : }
105 :
106 : /*@
107 : DSSetDimensions - Resize the matrices in the DS object.
108 :
109 : Logically Collective
110 :
111 : Input Parameters:
112 : + ds - the direct solver context
113 : . n - the new size
114 : . l - number of locked (inactive) leading columns
115 : - k - intermediate dimension (e.g., position of arrow)
116 :
117 : Notes:
118 : The internal arrays are not reallocated.
119 :
120 : Some DS types have additional dimensions, e.g. the number of columns
121 : in DSSVD. For these, you should call a specific interface function.
122 :
123 : Level: intermediate
124 :
125 : .seealso: DSGetDimensions(), DSAllocate(), DSSVDSetDimensions()
126 : @*/
127 31234 : PetscErrorCode DSSetDimensions(DS ds,PetscInt n,PetscInt l,PetscInt k)
128 : {
129 31234 : PetscInt on,ol,ok;
130 31234 : PetscBool issvd;
131 :
132 31234 : PetscFunctionBegin;
133 31234 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
134 31234 : DSCheckAlloc(ds,1);
135 124936 : PetscValidLogicalCollectiveInt(ds,n,2);
136 124936 : PetscValidLogicalCollectiveInt(ds,l,3);
137 124936 : PetscValidLogicalCollectiveInt(ds,k,4);
138 31234 : on = ds->n; ol = ds->l; ok = ds->k;
139 31234 : if (n==PETSC_DECIDE || n==PETSC_DEFAULT) {
140 131 : ds->n = ds->extrarow? ds->ld-1: ds->ld;
141 : } else {
142 31103 : PetscCheck(n>=0 && n<=ds->ld,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of n. Must be between 0 and ld");
143 31103 : PetscCall(PetscObjectTypeCompareAny((PetscObject)ds,&issvd,DSSVD,DSGSVD,"")); /* SVD and GSVD have extra column instead of extra row */
144 31103 : PetscCheck(!ds->extrarow || n<ds->ld || issvd,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"A value of n equal to ld leaves no room for extra row");
145 31103 : ds->n = n;
146 : }
147 31234 : ds->t = ds->n; /* truncated length equal to the new dimension */
148 31234 : if (l==PETSC_DECIDE || l==PETSC_DEFAULT) {
149 487 : ds->l = 0;
150 : } else {
151 30747 : PetscCheck(l>=0 && l<=ds->n,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of l. Must be between 0 and n");
152 30747 : ds->l = l;
153 : }
154 31234 : if (k==PETSC_DECIDE || k==PETSC_DEFAULT) {
155 1419 : ds->k = ds->n/2;
156 : } else {
157 29815 : PetscCheck(k>=0 || k<=ds->n,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of k. Must be between 0 and n");
158 29815 : ds->k = k;
159 : }
160 31234 : if (on!=ds->n || ol!=ds->l || ok!=ds->k) PetscCall(PetscInfo(ds,"New dimensions are: n=%" PetscInt_FMT ", l=%" PetscInt_FMT ", k=%" PetscInt_FMT "\n",ds->n,ds->l,ds->k));
161 31234 : PetscFunctionReturn(PETSC_SUCCESS);
162 : }
163 :
164 : /*@
165 : DSGetDimensions - Returns the current dimensions.
166 :
167 : Not Collective
168 :
169 : Input Parameter:
170 : . ds - the direct solver context
171 :
172 : Output Parameters:
173 : + n - the current size
174 : . l - number of locked (inactive) leading columns
175 : . k - intermediate dimension (e.g., position of arrow)
176 : - t - truncated length
177 :
178 : Note:
179 : The t parameter makes sense only if DSTruncate() has been called.
180 : Otherwise its value equals n.
181 :
182 : Some DS types have additional dimensions, e.g. the number of columns
183 : in DSSVD. For these, you should call a specific interface function.
184 :
185 : Level: intermediate
186 :
187 : .seealso: DSSetDimensions(), DSTruncate(), DSSVDGetDimensions()
188 : @*/
189 22251 : PetscErrorCode DSGetDimensions(DS ds,PetscInt *n,PetscInt *l,PetscInt *k,PetscInt *t)
190 : {
191 22251 : PetscFunctionBegin;
192 22251 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
193 22251 : DSCheckAlloc(ds,1);
194 22251 : if (n) *n = ds->n;
195 22251 : if (l) *l = ds->l;
196 22251 : if (k) *k = ds->k;
197 22251 : if (t) *t = ds->t;
198 22251 : PetscFunctionReturn(PETSC_SUCCESS);
199 : }
200 :
201 : /*@
202 : DSTruncate - Truncates the system represented in the DS object.
203 :
204 : Logically Collective
205 :
206 : Input Parameters:
207 : + ds - the direct solver context
208 : . n - the new size
209 : - trim - a flag to indicate if the factorization must be trimmed
210 :
211 : Note:
212 : The new size is set to n. Note that in some cases the new size could
213 : be n+1 or n-1 to avoid breaking a 2x2 diagonal block (e.g. in real
214 : Schur form). In cases where the extra row is meaningful, the first
215 : n elements are kept as the extra row for the new system.
216 :
217 : If the flag trim is turned on, it resets the locked and intermediate
218 : dimensions to zero, see DSSetDimensions(), and sets the state to RAW.
219 : It also cleans the extra row if being used.
220 :
221 : The typical usage of trim=true is to truncate the Schur decomposition
222 : at the end of a Krylov iteration. In this case, the state must be
223 : changed to RAW so that DSVectors() computes eigenvectors from scratch.
224 :
225 : Level: advanced
226 :
227 : .seealso: DSSetDimensions(), DSSetExtraRow(), DSStateType
228 : @*/
229 9026 : PetscErrorCode DSTruncate(DS ds,PetscInt n,PetscBool trim)
230 : {
231 9026 : DSStateType old;
232 :
233 9026 : PetscFunctionBegin;
234 9026 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
235 9026 : PetscValidType(ds,1);
236 9026 : DSCheckAlloc(ds,1);
237 36104 : PetscValidLogicalCollectiveInt(ds,n,2);
238 36104 : PetscValidLogicalCollectiveBool(ds,trim,3);
239 9026 : PetscCheck(n>=ds->l && n<=ds->n,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of n (%" PetscInt_FMT "). Must be between l (%" PetscInt_FMT ") and n (%" PetscInt_FMT ")",n,ds->l,ds->n);
240 9026 : PetscCall(PetscLogEventBegin(DS_Other,ds,0,0,0));
241 9026 : PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
242 9026 : PetscUseTypeMethod(ds,truncate,n,trim);
243 9026 : PetscCall(PetscFPTrapPop());
244 9026 : PetscCall(PetscLogEventEnd(DS_Other,ds,0,0,0));
245 17363 : PetscCall(PetscInfo(ds,"Decomposition %s to size n=%" PetscInt_FMT "\n",trim?"trimmed":"truncated",ds->n));
246 9026 : old = ds->state;
247 9026 : ds->state = trim? DS_STATE_RAW: DS_STATE_TRUNCATED;
248 9026 : if (old!=ds->state) PetscCall(PetscInfo(ds,"State has changed from %s to %s\n",DSStateTypes[old],DSStateTypes[ds->state]));
249 9026 : PetscCall(PetscObjectStateIncrease((PetscObject)ds));
250 9026 : PetscFunctionReturn(PETSC_SUCCESS);
251 : }
252 :
253 : /*@
254 : DSMatGetSize - Returns the numbers of rows and columns of one of the DS matrices.
255 :
256 : Not Collective
257 :
258 : Input Parameters:
259 : + ds - the direct solver context
260 : - t - the requested matrix
261 :
262 : Output Parameters:
263 : + n - the number of rows
264 : - m - the number of columns
265 :
266 : Note:
267 : This is equivalent to MatGetSize() on a matrix obtained with DSGetMat().
268 :
269 : Level: developer
270 :
271 : .seealso: DSSetDimensions(), DSGetMat()
272 : @*/
273 58634 : PetscErrorCode DSMatGetSize(DS ds,DSMatType t,PetscInt *m,PetscInt *n)
274 : {
275 58634 : PetscInt rows,cols;
276 :
277 58634 : PetscFunctionBegin;
278 58634 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
279 58634 : PetscValidType(ds,1);
280 58634 : DSCheckValidMat(ds,t,2);
281 58634 : if (ds->ops->matgetsize) PetscUseTypeMethod(ds,matgetsize,t,&rows,&cols);
282 : else {
283 46391 : if (ds->state==DS_STATE_TRUNCATED && t>=DS_MAT_Q) rows = ds->t;
284 39680 : else rows = (t==DS_MAT_A && ds->extrarow)? ds->n+1: ds->n;
285 46391 : if (t==DS_MAT_T) cols = PetscDefined(USE_COMPLEX)? 2: 3;
286 44026 : else if (t==DS_MAT_D) cols = 1;
287 40700 : else cols = ds->n;
288 : }
289 58634 : if (m) *m = rows;
290 58634 : if (n) *n = cols;
291 58634 : PetscFunctionReturn(PETSC_SUCCESS);
292 : }
293 :
294 : /*@
295 : DSMatIsHermitian - Checks if one of the DS matrices is known to be Hermitian.
296 :
297 : Not Collective
298 :
299 : Input Parameters:
300 : + ds - the direct solver context
301 : - t - the requested matrix
302 :
303 : Output Parameter:
304 : . flg - the Hermitian flag
305 :
306 : Note:
307 : Does not check the matrix values directly. The flag is set according to the
308 : problem structure. For instance, in DSHEP matrix A is Hermitian.
309 :
310 : Level: developer
311 :
312 : .seealso: DSGetMat()
313 : @*/
314 58606 : PetscErrorCode DSMatIsHermitian(DS ds,DSMatType t,PetscBool *flg)
315 : {
316 58606 : PetscFunctionBegin;
317 58606 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
318 58606 : PetscValidType(ds,1);
319 58606 : DSCheckValidMat(ds,t,2);
320 58606 : PetscAssertPointer(flg,3);
321 58606 : *flg = PETSC_FALSE;
322 58606 : PetscTryTypeMethod(ds,hermitian,t,flg);
323 58606 : PetscFunctionReturn(PETSC_SUCCESS);
324 : }
325 :
326 1632 : PetscErrorCode DSGetTruncateSize_Default(DS ds,PetscInt l,PetscInt n,PetscInt *k)
327 : {
328 : #if !defined(PETSC_USE_COMPLEX)
329 1632 : PetscScalar val;
330 : #endif
331 :
332 1632 : PetscFunctionBegin;
333 : #if !defined(PETSC_USE_COMPLEX)
334 1632 : PetscCall(MatGetValue(ds->omat[DS_MAT_A],l+(*k),l+(*k)-1,&val));
335 1632 : if (val != 0.0) {
336 433 : if (l+(*k)<n-1) (*k)++;
337 0 : else (*k)--;
338 : }
339 : #endif
340 1632 : PetscFunctionReturn(PETSC_SUCCESS);
341 : }
342 :
343 : /*@
344 : DSGetTruncateSize - Gets the correct size to be used in DSTruncate()
345 : to avoid breaking a 2x2 block.
346 :
347 : Not Collective
348 :
349 : Input Parameters:
350 : + ds - the direct solver context
351 : . l - the size of the locked part (set to 0 to use ds->l)
352 : - n - the total matrix size (set to 0 to use ds->n)
353 :
354 : Output Parameter:
355 : . k - the wanted truncation size (possibly modified)
356 :
357 : Notes:
358 : This should be called before DSTruncate() to make sure that the truncation
359 : does not break a 2x2 block corresponding to a complex conjugate eigenvalue.
360 :
361 : The total size is n (either user-provided or ds->n if 0 is passed). The
362 : size where the truncation is intended is equal to l+k (where l can be
363 : equal to the locked size ds->l if set to 0). Then if there is a 2x2 block
364 : at the l+k limit, the value of k is increased (or decreased) by 1.
365 :
366 : Level: advanced
367 :
368 : .seealso: DSTruncate(), DSSetDimensions()
369 : @*/
370 4939 : PetscErrorCode DSGetTruncateSize(DS ds,PetscInt l,PetscInt n,PetscInt *k)
371 : {
372 4939 : PetscFunctionBegin;
373 4939 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
374 4939 : DSCheckAlloc(ds,1);
375 19756 : PetscValidLogicalCollectiveInt(ds,l,2);
376 19756 : PetscValidLogicalCollectiveInt(ds,n,3);
377 4939 : PetscAssertPointer(k,4);
378 4939 : PetscUseTypeMethod(ds,gettruncatesize,l?l:ds->l,n?n:ds->n,k);
379 4939 : PetscFunctionReturn(PETSC_SUCCESS);
380 : }
381 :
382 : /*@
383 : DSGetMat - Returns a sequential dense Mat object containing the requested
384 : matrix.
385 :
386 : Not Collective
387 :
388 : Input Parameters:
389 : + ds - the direct solver context
390 : - m - the requested matrix
391 :
392 : Output Parameter:
393 : . A - Mat object
394 :
395 : Notes:
396 : The returned Mat has sizes equal to the current DS dimensions (nxm),
397 : and contains the values that would be obtained with DSGetArray()
398 : (not DSGetArrayReal()). If the DS was truncated, then the number of rows
399 : is equal to the dimension prior to truncation, see DSTruncate().
400 : The communicator is always PETSC_COMM_SELF.
401 :
402 : It is implemented with MatDenseGetSubMatrix(), and when no longer needed
403 : the user must call DSRestoreMat() which will invoke MatDenseRestoreSubMatrix().
404 :
405 : For matrices DS_MAT_T and DS_MAT_D, this function will return a Mat object
406 : that cannot be used directly for computations, since it uses compact storage
407 : (three and one diagonals for T and D, respectively). In complex scalars, the
408 : internal array stores real values, so it is sufficient with 2 columns for T.
409 :
410 : Level: advanced
411 :
412 : .seealso: DSRestoreMat(), DSSetDimensions(), DSGetArray(), DSGetArrayReal(), DSTruncate()
413 : @*/
414 58606 : PetscErrorCode DSGetMat(DS ds,DSMatType m,Mat *A)
415 : {
416 58606 : PetscInt rows,cols;
417 58606 : PetscBool flg;
418 :
419 58606 : PetscFunctionBegin;
420 58606 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
421 58606 : DSCheckAlloc(ds,1);
422 58606 : DSCheckValidMat(ds,m,2);
423 58606 : PetscAssertPointer(A,3);
424 :
425 58606 : PetscCall(DSMatGetSize(ds,m,&rows,&cols));
426 58606 : PetscCheck(rows && cols,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"Must call DSSetDimensions() first");
427 58606 : PetscCall(MatDenseGetSubMatrix(ds->omat[m],0,rows,0,cols,A));
428 :
429 : /* set Hermitian flag */
430 58606 : PetscCall(DSMatIsHermitian(ds,m,&flg));
431 58606 : PetscCall(MatSetOption(*A,MAT_SYMMETRIC,flg));
432 58606 : PetscCall(MatSetOption(*A,MAT_HERMITIAN,flg));
433 58606 : PetscFunctionReturn(PETSC_SUCCESS);
434 : }
435 :
436 : /*@
437 : DSRestoreMat - Restores the matrix after DSGetMat() was called.
438 :
439 : Not Collective
440 :
441 : Input Parameters:
442 : + ds - the direct solver context
443 : . m - the requested matrix
444 : - A - the fetched Mat object
445 :
446 : Note:
447 : A call to this function must match a previous call of DSGetMat().
448 :
449 : Level: advanced
450 :
451 : .seealso: DSGetMat(), DSRestoreArray(), DSRestoreArrayReal()
452 : @*/
453 58606 : PetscErrorCode DSRestoreMat(DS ds,DSMatType m,Mat *A)
454 : {
455 58606 : PetscFunctionBegin;
456 58606 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
457 58606 : DSCheckAlloc(ds,1);
458 58606 : DSCheckValidMat(ds,m,2);
459 58606 : PetscAssertPointer(A,3);
460 :
461 58606 : PetscCall(MatDenseRestoreSubMatrix(ds->omat[m],A));
462 58606 : PetscCall(PetscObjectStateIncrease((PetscObject)ds));
463 58606 : PetscFunctionReturn(PETSC_SUCCESS);
464 : }
465 :
466 : /*@
467 : DSGetMatAndColumn - Returns a sequential dense Mat object containing the requested
468 : matrix and one of its columns as a Vec.
469 :
470 : Not Collective
471 :
472 : Input Parameters:
473 : + ds - the direct solver context
474 : . m - the requested matrix
475 : - col - the requested column
476 :
477 : Output Parameters:
478 : + A - Mat object
479 : - v - Vec object (the column)
480 :
481 : Notes:
482 : This calls DSGetMat() and then it extracts the selected column.
483 : The user must call DSRestoreMatAndColumn() to recover the original state.
484 : For matrices DS_MAT_T and DS_MAT_D, in complex scalars this function implies
485 : copying from real values stored internally to scalar values in the Vec.
486 :
487 : Level: advanced
488 :
489 : .seealso: DSRestoreMatAndColumn(), DSGetMat()
490 : @*/
491 3690 : PetscErrorCode DSGetMatAndColumn(DS ds,DSMatType m,PetscInt col,Mat *A,Vec *v)
492 : {
493 3690 : PetscFunctionBegin;
494 3690 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
495 3690 : DSCheckAlloc(ds,1);
496 3690 : DSCheckValidMat(ds,m,2);
497 3690 : PetscAssertPointer(A,4);
498 3690 : PetscAssertPointer(v,5);
499 :
500 3690 : PetscCall(DSGetMat(ds,m,A));
501 3690 : if (PetscDefined(USE_COMPLEX) && (m==DS_MAT_T || m==DS_MAT_D)) {
502 : const PetscScalar *as;
503 : PetscScalar *vs;
504 : PetscReal *ar;
505 : PetscInt i,n,lda;
506 : PetscCall(MatCreateVecs(*A,NULL,v));
507 : PetscCall(VecGetSize(*v,&n));
508 : PetscCall(MatDenseGetLDA(*A,&lda));
509 : PetscCall(MatDenseGetArrayRead(*A,&as));
510 : PetscCall(VecGetArrayWrite(*v,&vs));
511 : ar = (PetscReal*)as;
512 : for (i=0;i<n;i++) vs[i] = ar[i+col*lda];
513 : PetscCall(VecRestoreArrayWrite(*v,&vs));
514 : PetscCall(MatDenseRestoreArrayRead(*A,&as));
515 3690 : } else PetscCall(MatDenseGetColumnVec(*A,col,v));
516 3690 : PetscFunctionReturn(PETSC_SUCCESS);
517 : }
518 :
519 : /*@
520 : DSRestoreMatAndColumn - Restores the matrix and vector after DSGetMatAndColumn()
521 : was called.
522 :
523 : Not Collective
524 :
525 : Input Parameters:
526 : + ds - the direct solver context
527 : . m - the requested matrix
528 : . col - the requested column
529 : . A - the fetched Mat object
530 : - v - the fetched Vec object
531 :
532 : Note:
533 : A call to this function must match a previous call of DSGetMatAndColumn().
534 :
535 : Level: advanced
536 :
537 : .seealso: DSGetMatAndColumn(), DSRestoreMat()
538 : @*/
539 3690 : PetscErrorCode DSRestoreMatAndColumn(DS ds,DSMatType m,PetscInt col,Mat *A,Vec *v)
540 : {
541 3690 : PetscFunctionBegin;
542 3690 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
543 3690 : DSCheckAlloc(ds,1);
544 3690 : DSCheckValidMat(ds,m,2);
545 3690 : PetscAssertPointer(A,4);
546 3690 : PetscAssertPointer(v,5);
547 :
548 3690 : if (PetscDefined(USE_COMPLEX) && (m==DS_MAT_T || m==DS_MAT_D)) {
549 : const PetscScalar *vs;
550 : PetscScalar *as;
551 : PetscReal *ar;
552 : PetscInt i,n,lda;
553 : PetscCall(VecGetSize(*v,&n));
554 : PetscCall(MatDenseGetLDA(*A,&lda));
555 : PetscCall(MatDenseGetArray(*A,&as));
556 : PetscCall(VecGetArrayRead(*v,&vs));
557 : ar = (PetscReal*)as;
558 : for (i=0;i<n;i++) ar[i+col*lda] = PetscRealPart(vs[i]);
559 : PetscCall(VecRestoreArrayRead(*v,&vs));
560 : PetscCall(MatDenseRestoreArray(*A,&as));
561 : PetscCall(VecDestroy(v));
562 3690 : } else PetscCall(MatDenseRestoreColumnVec(*A,col,v));
563 3690 : PetscCall(DSRestoreMat(ds,m,A));
564 3690 : PetscFunctionReturn(PETSC_SUCCESS);
565 : }
566 :
567 : /*@C
568 : DSGetArray - Returns a pointer to the internal array of one of the
569 : matrices. You MUST call DSRestoreArray() when you no longer
570 : need to access the array.
571 :
572 : Not Collective
573 :
574 : Input Parameters:
575 : + ds - the direct solver context
576 : - m - the requested matrix
577 :
578 : Output Parameter:
579 : . a - pointer to the values
580 :
581 : Note:
582 : To get read-only access, use DSGetMat() followed by MatDenseGetArrayRead().
583 :
584 : Level: advanced
585 :
586 : .seealso: DSRestoreArray(), DSGetArrayReal()
587 : @*/
588 62860 : PetscErrorCode DSGetArray(DS ds,DSMatType m,PetscScalar *a[])
589 : {
590 62860 : PetscFunctionBegin;
591 62860 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
592 62860 : DSCheckAlloc(ds,1);
593 62860 : DSCheckValidMat(ds,m,2);
594 62860 : PetscAssertPointer(a,3);
595 62860 : PetscCall(MatDenseGetArray(ds->omat[m],a));
596 62860 : PetscFunctionReturn(PETSC_SUCCESS);
597 : }
598 :
599 : /*@C
600 : DSRestoreArray - Restores the matrix after DSGetArray() was called.
601 :
602 : Not Collective
603 :
604 : Input Parameters:
605 : + ds - the direct solver context
606 : . m - the requested matrix
607 : - a - pointer to the values
608 :
609 : Level: advanced
610 :
611 : .seealso: DSGetArray(), DSGetArrayReal()
612 : @*/
613 62759 : PetscErrorCode DSRestoreArray(DS ds,DSMatType m,PetscScalar *a[])
614 : {
615 62759 : PetscFunctionBegin;
616 62759 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
617 62759 : DSCheckAlloc(ds,1);
618 62759 : DSCheckValidMat(ds,m,2);
619 62759 : PetscAssertPointer(a,3);
620 62759 : PetscCall(MatDenseRestoreArray(ds->omat[m],a));
621 62759 : PetscCall(PetscObjectStateIncrease((PetscObject)ds));
622 62759 : PetscFunctionReturn(PETSC_SUCCESS);
623 : }
624 :
625 : /*@C
626 : DSGetArrayReal - Returns a real pointer to the internal array of T or D.
627 : You MUST call DSRestoreArrayReal() when you no longer need to access the array.
628 :
629 : Not Collective
630 :
631 : Input Parameters:
632 : + ds - the direct solver context
633 : - m - the requested matrix
634 :
635 : Output Parameter:
636 : . a - pointer to the values
637 :
638 : Note:
639 : This function can be used only for DS_MAT_T and DS_MAT_D. These matrices always
640 : store real values, even in complex scalars, that is why the returned pointer is
641 : PetscReal.
642 :
643 : Level: advanced
644 :
645 : .seealso: DSRestoreArrayReal(), DSGetArray()
646 : @*/
647 137795 : PetscErrorCode DSGetArrayReal(DS ds,DSMatType m,PetscReal *a[])
648 : {
649 : #if defined(PETSC_USE_COMPLEX)
650 : PetscScalar *as;
651 : #endif
652 :
653 137795 : PetscFunctionBegin;
654 137795 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
655 137795 : DSCheckAlloc(ds,1);
656 137795 : DSCheckValidMatReal(ds,m,2);
657 137795 : PetscAssertPointer(a,3);
658 : #if defined(PETSC_USE_COMPLEX)
659 : PetscCall(MatDenseGetArray(ds->omat[m],&as));
660 : *a = (PetscReal*)as;
661 : #else
662 137795 : PetscCall(MatDenseGetArray(ds->omat[m],a));
663 : #endif
664 137795 : PetscFunctionReturn(PETSC_SUCCESS);
665 : }
666 :
667 : /*@C
668 : DSRestoreArrayReal - Restores the matrix after DSGetArrayReal() was called.
669 :
670 : Not Collective
671 :
672 : Input Parameters:
673 : + ds - the direct solver context
674 : . m - the requested matrix
675 : - a - pointer to the values
676 :
677 : Level: advanced
678 :
679 : .seealso: DSGetArrayReal(), DSGetArray()
680 : @*/
681 137795 : PetscErrorCode DSRestoreArrayReal(DS ds,DSMatType m,PetscReal *a[])
682 : {
683 : #if defined(PETSC_USE_COMPLEX)
684 : PetscScalar *as;
685 : #endif
686 :
687 137795 : PetscFunctionBegin;
688 137795 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
689 137795 : DSCheckAlloc(ds,1);
690 137795 : DSCheckValidMatReal(ds,m,2);
691 137795 : PetscAssertPointer(a,3);
692 : #if defined(PETSC_USE_COMPLEX)
693 : PetscCall(MatDenseRestoreArray(ds->omat[m],&as));
694 : *a = NULL;
695 : #else
696 137795 : PetscCall(MatDenseRestoreArray(ds->omat[m],a));
697 : #endif
698 137795 : PetscCall(PetscObjectStateIncrease((PetscObject)ds));
699 137795 : PetscFunctionReturn(PETSC_SUCCESS);
700 : }
701 :
702 : /*@
703 : DSSolve - Solves the problem.
704 :
705 : Logically Collective
706 :
707 : Input Parameters:
708 : + ds - the direct solver context
709 : . eigr - array to store the computed eigenvalues (real part)
710 : - eigi - array to store the computed eigenvalues (imaginary part)
711 :
712 : Note:
713 : This call brings the dense system to condensed form. No ordering
714 : of the eigenvalues is enforced (for this, call DSSort() afterwards).
715 :
716 : Level: intermediate
717 :
718 : .seealso: DSSort(), DSStateType
719 : @*/
720 22423 : PetscErrorCode DSSolve(DS ds,PetscScalar eigr[],PetscScalar eigi[])
721 : {
722 22423 : PetscFunctionBegin;
723 22423 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
724 22423 : PetscValidType(ds,1);
725 22423 : DSCheckAlloc(ds,1);
726 22423 : PetscAssertPointer(eigr,2);
727 22423 : if (ds->state>=DS_STATE_CONDENSED) PetscFunctionReturn(PETSC_SUCCESS);
728 22398 : PetscCheck(ds->ops->solve[ds->method],PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"The specified method number does not exist for this DS");
729 22398 : PetscCall(PetscInfo(ds,"Starting solve with problem sizes: n=%" PetscInt_FMT ", l=%" PetscInt_FMT ", k=%" PetscInt_FMT "\n",ds->n,ds->l,ds->k));
730 22398 : PetscCall(PetscLogEventBegin(DS_Solve,ds,0,0,0));
731 22398 : PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
732 22398 : PetscUseTypeMethod(ds,solve[ds->method],eigr,eigi);
733 22398 : PetscCall(PetscFPTrapPop());
734 22398 : PetscCall(PetscLogEventEnd(DS_Solve,ds,0,0,0));
735 22398 : PetscCall(PetscInfo(ds,"State has changed from %s to CONDENSED\n",DSStateTypes[ds->state]));
736 22398 : ds->state = DS_STATE_CONDENSED;
737 22398 : PetscCall(PetscObjectStateIncrease((PetscObject)ds));
738 22398 : PetscFunctionReturn(PETSC_SUCCESS);
739 : }
740 :
741 : /*@C
742 : DSSort - Sorts the result of DSSolve() according to a given sorting
743 : criterion.
744 :
745 : Logically Collective
746 :
747 : Input Parameters:
748 : + ds - the direct solver context
749 : . rr - (optional) array containing auxiliary values (real part)
750 : - ri - (optional) array containing auxiliary values (imaginary part)
751 :
752 : Input/Output Parameters:
753 : + eigr - array containing the computed eigenvalues (real part)
754 : . eigi - array containing the computed eigenvalues (imaginary part)
755 : - k - (optional) number of elements in the leading group
756 :
757 : Notes:
758 : This routine sorts the arrays provided in eigr and eigi, and also
759 : sorts the dense system stored inside ds (assumed to be in condensed form).
760 : The sorting criterion is specified with DSSetSlepcSC().
761 :
762 : If arrays rr and ri are provided, then a (partial) reordering based on these
763 : values rather than on the eigenvalues is performed. In symmetric problems
764 : a total order is obtained (parameter k is ignored), but otherwise the result
765 : is sorted only partially. In this latter case, it is only guaranteed that
766 : all the first k elements satisfy the comparison with any of the last n-k
767 : elements. The output value of parameter k is the final number of elements in
768 : the first set.
769 :
770 : Level: intermediate
771 :
772 : .seealso: DSSolve(), DSSetSlepcSC(), DSSortWithPermutation()
773 : @*/
774 29819 : PetscErrorCode DSSort(DS ds,PetscScalar *eigr,PetscScalar *eigi,PetscScalar *rr,PetscScalar *ri,PetscInt *k)
775 : {
776 29819 : PetscInt i;
777 :
778 29819 : PetscFunctionBegin;
779 29819 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
780 29819 : PetscValidType(ds,1);
781 29819 : DSCheckSolved(ds,1);
782 29819 : PetscAssertPointer(eigr,2);
783 29819 : if (rr) PetscAssertPointer(rr,4);
784 29819 : PetscCheck(ds->state<DS_STATE_TRUNCATED,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"Cannot sort a truncated DS");
785 29819 : PetscCheck(ds->sc,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"Must provide a sorting criterion first");
786 29819 : PetscCheck(!k || rr,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONG,"Argument k can only be used together with rr");
787 :
788 476754 : for (i=0;i<ds->n;i++) ds->perm[i] = i; /* initialize to trivial permutation */
789 29819 : PetscCall(PetscLogEventBegin(DS_Other,ds,0,0,0));
790 29819 : PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
791 29819 : PetscUseTypeMethod(ds,sort,eigr,eigi,rr,ri,k);
792 29819 : PetscCall(PetscFPTrapPop());
793 29819 : PetscCall(PetscLogEventEnd(DS_Other,ds,0,0,0));
794 29819 : PetscCall(PetscObjectStateIncrease((PetscObject)ds));
795 29819 : PetscCall(PetscInfo(ds,"Finished sorting\n"));
796 29819 : PetscFunctionReturn(PETSC_SUCCESS);
797 : }
798 :
799 : /*@C
800 : DSSortWithPermutation - Reorders the result of DSSolve() according to a given
801 : permutation.
802 :
803 : Logically Collective
804 :
805 : Input Parameters:
806 : + ds - the direct solver context
807 : - perm - permutation that indicates the new ordering
808 :
809 : Input/Output Parameters:
810 : + eigr - array with the reordered eigenvalues (real part)
811 : - eigi - array with the reordered eigenvalues (imaginary part)
812 :
813 : Notes:
814 : This routine reorders the arrays provided in eigr and eigi, and also the dense
815 : system stored inside ds (assumed to be in condensed form). There is no sorting
816 : criterion, as opposed to DSSort(). Instead, the new ordering is given in argument perm.
817 :
818 : Level: advanced
819 :
820 : .seealso: DSSolve(), DSSort()
821 : @*/
822 1 : PetscErrorCode DSSortWithPermutation(DS ds,PetscInt *perm,PetscScalar *eigr,PetscScalar *eigi)
823 : {
824 1 : PetscFunctionBegin;
825 1 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
826 1 : PetscValidType(ds,1);
827 1 : DSCheckSolved(ds,1);
828 1 : PetscAssertPointer(perm,2);
829 1 : PetscAssertPointer(eigr,3);
830 1 : PetscCheck(ds->state<DS_STATE_TRUNCATED,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"Cannot sort a truncated DS");
831 :
832 1 : PetscCall(PetscLogEventBegin(DS_Other,ds,0,0,0));
833 1 : PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
834 1 : PetscUseTypeMethod(ds,sortperm,perm,eigr,eigi);
835 1 : PetscCall(PetscFPTrapPop());
836 1 : PetscCall(PetscLogEventEnd(DS_Other,ds,0,0,0));
837 1 : PetscCall(PetscObjectStateIncrease((PetscObject)ds));
838 1 : PetscCall(PetscInfo(ds,"Finished sorting\n"));
839 1 : PetscFunctionReturn(PETSC_SUCCESS);
840 : }
841 :
842 : /*@
843 : DSSynchronize - Make sure that all processes have the same data, performing
844 : communication if necessary.
845 :
846 : Collective
847 :
848 : Input Parameter:
849 : . ds - the direct solver context
850 :
851 : Input/Output Parameters:
852 : + eigr - (optional) array with the computed eigenvalues (real part)
853 : - eigi - (optional) array with the computed eigenvalues (imaginary part)
854 :
855 : Notes:
856 : When the DS has been created with a communicator with more than one process,
857 : the internal data, especially the computed matrices, may diverge in the
858 : different processes. This happens when using multithreaded BLAS and may
859 : cause numerical issues in some ill-conditioned problems. This function
860 : performs the necessary communication among the processes so that the
861 : internal data is exactly equal in all of them.
862 :
863 : Depending on the parallel mode as set with DSSetParallel(), this function
864 : will either do nothing or synchronize the matrices computed by DSSolve()
865 : and DSSort(). The arguments eigr and eigi are typically those used in the
866 : calls to DSSolve() and DSSort().
867 :
868 : Level: developer
869 :
870 : .seealso: DSSetParallel(), DSSolve(), DSSort()
871 : @*/
872 22171 : PetscErrorCode DSSynchronize(DS ds,PetscScalar eigr[],PetscScalar eigi[])
873 : {
874 22171 : PetscMPIInt size;
875 :
876 22171 : PetscFunctionBegin;
877 22171 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
878 22171 : PetscValidType(ds,1);
879 22171 : DSCheckAlloc(ds,1);
880 22171 : PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)ds),&size));
881 22171 : if (size>1 && ds->pmode==DS_PARALLEL_SYNCHRONIZED) {
882 340 : PetscCall(PetscLogEventBegin(DS_Synchronize,ds,0,0,0));
883 340 : PetscUseTypeMethod(ds,synchronize,eigr,eigi);
884 340 : PetscCall(PetscLogEventEnd(DS_Synchronize,ds,0,0,0));
885 340 : PetscCall(PetscObjectStateIncrease((PetscObject)ds));
886 340 : PetscCall(PetscInfo(ds,"Synchronization completed (%s)\n",DSParallelTypes[ds->pmode]));
887 : }
888 22171 : PetscFunctionReturn(PETSC_SUCCESS);
889 : }
890 :
891 : /*@C
892 : DSVectors - Compute vectors associated to the dense system such
893 : as eigenvectors.
894 :
895 : Logically Collective
896 :
897 : Input Parameters:
898 : + ds - the direct solver context
899 : - mat - the matrix, used to indicate which vectors are required
900 :
901 : Input/Output Parameter:
902 : . j - (optional) index of vector to be computed
903 :
904 : Output Parameter:
905 : . rnorm - (optional) computed residual norm
906 :
907 : Notes:
908 : Allowed values for mat are DS_MAT_X, DS_MAT_Y, DS_MAT_U and DS_MAT_V, to
909 : compute right or left eigenvectors, or left or right singular vectors,
910 : respectively.
911 :
912 : If NULL is passed in argument j then all vectors are computed,
913 : otherwise j indicates which vector must be computed. In real non-symmetric
914 : problems, on exit the index j will be incremented when a complex conjugate
915 : pair is found.
916 :
917 : This function can be invoked after the dense problem has been solved,
918 : to get the residual norm estimate of the associated Ritz pair. In that
919 : case, the relevant information is returned in rnorm.
920 :
921 : For computing eigenvectors, LAPACK's _trevc is used so the matrix must
922 : be in (quasi-)triangular form, or call DSSolve() first.
923 :
924 : Level: intermediate
925 :
926 : .seealso: DSSolve()
927 : @*/
928 46428 : PetscErrorCode DSVectors(DS ds,DSMatType mat,PetscInt *j,PetscReal *rnorm)
929 : {
930 46428 : PetscFunctionBegin;
931 46428 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
932 46428 : PetscValidType(ds,1);
933 46428 : DSCheckAlloc(ds,1);
934 185712 : PetscValidLogicalCollectiveEnum(ds,mat,2);
935 46428 : PetscCheck(mat<DS_NUM_MAT,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONG,"Invalid matrix");
936 46428 : PetscCheck(!rnorm || j,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"Must give a value of j");
937 46428 : if (!ds->omat[mat]) PetscCall(DSAllocateMat_Private(ds,mat));
938 46428 : if (!j) PetscCall(PetscInfo(ds,"Computing all vectors on %s\n",DSMatName[mat]));
939 46428 : PetscCall(PetscLogEventBegin(DS_Vectors,ds,0,0,0));
940 46428 : PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
941 46428 : PetscUseTypeMethod(ds,vectors,mat,j,rnorm);
942 46428 : PetscCall(PetscFPTrapPop());
943 46428 : PetscCall(PetscLogEventEnd(DS_Vectors,ds,0,0,0));
944 46428 : PetscCall(PetscObjectStateIncrease((PetscObject)ds));
945 46428 : PetscFunctionReturn(PETSC_SUCCESS);
946 : }
947 :
948 : /*@
949 : DSUpdateExtraRow - Performs all necessary operations so that the extra
950 : row gets up-to-date after a call to DSSolve().
951 :
952 : Logically Collective
953 :
954 : Input Parameters:
955 : . ds - the direct solver context
956 :
957 : Level: advanced
958 :
959 : .seealso: DSSolve(), DSSetExtraRow()
960 : @*/
961 10488 : PetscErrorCode DSUpdateExtraRow(DS ds)
962 : {
963 10488 : PetscFunctionBegin;
964 10488 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
965 10488 : PetscValidType(ds,1);
966 10488 : DSCheckAlloc(ds,1);
967 10488 : PetscCheck(ds->extrarow,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONGSTATE,"Should have called DSSetExtraRow");
968 10488 : PetscCall(PetscInfo(ds,"Updating extra row\n"));
969 10488 : PetscCall(PetscLogEventBegin(DS_Other,ds,0,0,0));
970 10488 : PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
971 10488 : PetscUseTypeMethod(ds,update);
972 10488 : PetscCall(PetscFPTrapPop());
973 10488 : PetscCall(PetscLogEventEnd(DS_Other,ds,0,0,0));
974 10488 : PetscFunctionReturn(PETSC_SUCCESS);
975 : }
976 :
977 : /*@
978 : DSCond - Compute the condition number.
979 :
980 : Logically Collective
981 :
982 : Input Parameters:
983 : + ds - the direct solver context
984 : - cond - the computed condition number
985 :
986 : Notes:
987 : In standard eigenvalue problems, returns the inf-norm condition number of the first
988 : matrix, computed as cond(A) = norm(A)*norm(inv(A)).
989 :
990 : In GSVD problems, returns the maximum of cond(A) and cond(B), where cond(.) is
991 : computed as the ratio of the largest and smallest singular values.
992 :
993 : Does not take into account the extra row.
994 :
995 : Level: advanced
996 :
997 : .seealso: DSSolve(), DSSetExtraRow()
998 : @*/
999 164 : PetscErrorCode DSCond(DS ds,PetscReal *cond)
1000 : {
1001 164 : PetscFunctionBegin;
1002 164 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
1003 164 : PetscValidType(ds,1);
1004 164 : DSCheckAlloc(ds,1);
1005 164 : PetscAssertPointer(cond,2);
1006 164 : PetscCall(PetscLogEventBegin(DS_Other,ds,0,0,0));
1007 164 : PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
1008 164 : PetscUseTypeMethod(ds,cond,cond);
1009 164 : PetscCall(PetscFPTrapPop());
1010 164 : PetscCall(PetscLogEventEnd(DS_Other,ds,0,0,0));
1011 164 : PetscCall(PetscInfo(ds,"Computed condition number = %g\n",(double)*cond));
1012 164 : PetscFunctionReturn(PETSC_SUCCESS);
1013 : }
1014 :
1015 : /*@C
1016 : DSTranslateHarmonic - Computes a translation of the dense system.
1017 :
1018 : Logically Collective
1019 :
1020 : Input Parameters:
1021 : + ds - the direct solver context
1022 : . tau - the translation amount
1023 : . beta - last component of vector b
1024 : - recover - boolean flag to indicate whether to recover or not
1025 :
1026 : Output Parameters:
1027 : + g - the computed vector (optional)
1028 : - gamma - scale factor (optional)
1029 :
1030 : Notes:
1031 : This function is intended for use in the context of Krylov methods only.
1032 : It computes a translation of a Krylov decomposition in order to extract
1033 : eigenpair approximations by harmonic Rayleigh-Ritz.
1034 : The matrix is updated as A + g*b' where g = (A-tau*eye(n))'\b and
1035 : vector b is assumed to be beta*e_n^T.
1036 :
1037 : The gamma factor is defined as sqrt(1+g'*g) and can be interpreted as
1038 : the factor by which the residual of the Krylov decomposition is scaled.
1039 :
1040 : If the recover flag is activated, the computed translation undoes the
1041 : translation done previously. In that case, parameter tau is ignored.
1042 :
1043 : Level: developer
1044 :
1045 : .seealso: DSTranslateRKS()
1046 : @*/
1047 81 : PetscErrorCode DSTranslateHarmonic(DS ds,PetscScalar tau,PetscReal beta,PetscBool recover,PetscScalar *g,PetscReal *gamma)
1048 : {
1049 81 : PetscFunctionBegin;
1050 81 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
1051 81 : PetscValidType(ds,1);
1052 81 : DSCheckAlloc(ds,1);
1053 81 : if (recover) PetscCall(PetscInfo(ds,"Undoing the translation\n"));
1054 57 : else PetscCall(PetscInfo(ds,"Computing the translation\n"));
1055 81 : PetscCall(PetscLogEventBegin(DS_Other,ds,0,0,0));
1056 81 : PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
1057 81 : PetscUseTypeMethod(ds,transharm,tau,beta,recover,g,gamma);
1058 81 : PetscCall(PetscFPTrapPop());
1059 81 : PetscCall(PetscLogEventEnd(DS_Other,ds,0,0,0));
1060 81 : ds->state = DS_STATE_RAW;
1061 81 : PetscCall(PetscObjectStateIncrease((PetscObject)ds));
1062 81 : PetscFunctionReturn(PETSC_SUCCESS);
1063 : }
1064 :
1065 : /*@
1066 : DSTranslateRKS - Computes a modification of the dense system corresponding
1067 : to an update of the shift in a rational Krylov method.
1068 :
1069 : Logically Collective
1070 :
1071 : Input Parameters:
1072 : + ds - the direct solver context
1073 : - alpha - the translation amount
1074 :
1075 : Notes:
1076 : This function is intended for use in the context of Krylov methods only.
1077 : It takes the leading (k+1,k) submatrix of A, containing the truncated
1078 : Rayleigh quotient of a Krylov-Schur relation computed from a shift
1079 : sigma1 and transforms it to obtain a Krylov relation as if computed
1080 : from a different shift sigma2. The new matrix is computed as
1081 : 1.0/alpha*(eye(k)-Q*inv(R)), where [Q,R]=qr(eye(k)-alpha*A) and
1082 : alpha = sigma1-sigma2.
1083 :
1084 : Matrix Q is placed in DS_MAT_Q so that it can be used to update the
1085 : Krylov basis.
1086 :
1087 : Level: developer
1088 :
1089 : .seealso: DSTranslateHarmonic()
1090 : @*/
1091 12 : PetscErrorCode DSTranslateRKS(DS ds,PetscScalar alpha)
1092 : {
1093 12 : PetscFunctionBegin;
1094 12 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
1095 12 : PetscValidType(ds,1);
1096 12 : DSCheckAlloc(ds,1);
1097 12 : PetscCall(PetscInfo(ds,"Translating with alpha=%g\n",(double)PetscRealPart(alpha)));
1098 12 : PetscCall(PetscLogEventBegin(DS_Other,ds,0,0,0));
1099 12 : PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
1100 12 : PetscUseTypeMethod(ds,transrks,alpha);
1101 12 : PetscCall(PetscFPTrapPop());
1102 12 : PetscCall(PetscLogEventEnd(DS_Other,ds,0,0,0));
1103 12 : ds->state = DS_STATE_RAW;
1104 12 : ds->compact = PETSC_FALSE;
1105 12 : PetscCall(PetscObjectStateIncrease((PetscObject)ds));
1106 12 : PetscFunctionReturn(PETSC_SUCCESS);
1107 : }
|