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