1 /*
2 Functions in this file reorder the Ritz values in the (modified) Leja order.
3 */
4 #include <../src/ksp/ksp/impls/gmres/agmres/agmresimpl.h>
5
KSPAGMRESLejafmaxarray(PetscScalar * re,PetscInt pt,PetscInt n,PetscInt * pos)6 static PetscErrorCode KSPAGMRESLejafmaxarray(PetscScalar *re, PetscInt pt, PetscInt n, PetscInt *pos)
7 {
8 PetscInt i;
9 PetscScalar mx;
10
11 PetscFunctionBegin;
12 mx = re[0];
13 *pos = 0;
14 for (i = pt; i < n; i++) {
15 if (mx < re[i]) {
16 mx = re[i];
17 *pos = i;
18 }
19 }
20 PetscFunctionReturn(PETSC_SUCCESS);
21 }
22
KSPAGMRESLejaCfpdMax(PetscScalar * rm,PetscScalar * im,PetscInt * spos,PetscInt nbre,PetscInt n,PetscInt * rpos)23 static PetscErrorCode KSPAGMRESLejaCfpdMax(PetscScalar *rm, PetscScalar *im, PetscInt *spos, PetscInt nbre, PetscInt n, PetscInt *rpos)
24 {
25 PetscScalar rd, id, pd, max;
26 PetscInt i, j;
27
28 PetscFunctionBegin;
29 pd = 1.0;
30 max = 0.0;
31 *rpos = 0;
32 for (i = 0; i < n; i++) {
33 for (j = 0; j < nbre; j++) {
34 rd = rm[i] - rm[spos[j]];
35 id = im[i] - im[spos[j]];
36 pd = pd * PetscSqrtReal(rd * rd + id * id);
37 }
38 if (max < pd) {
39 *rpos = i;
40 max = pd;
41 }
42 pd = 1.0;
43 }
44 PetscFunctionReturn(PETSC_SUCCESS);
45 }
46
KSPAGMRESLejaOrdering(PetscScalar * re,PetscScalar * im,PetscScalar * rre,PetscScalar * rim,PetscInt m)47 PetscErrorCode KSPAGMRESLejaOrdering(PetscScalar *re, PetscScalar *im, PetscScalar *rre, PetscScalar *rim, PetscInt m)
48 {
49 PetscInt *spos;
50 PetscScalar *n_cmpl, temp;
51 PetscInt i, pos, j;
52
53 PetscFunctionBegin;
54 PetscCall(PetscMalloc1(m, &n_cmpl));
55 PetscCall(PetscMalloc1(m, &spos));
56 /* Check the proper order of complex conjugate pairs */
57 j = 0;
58 while (j < m) {
59 if (im[j] != 0.0) { /* complex eigenvalue */
60 if (im[j] < 0.0) { /* change the order */
61 temp = im[j + 1];
62 im[j + 1] = im[j];
63 im[j] = temp;
64 }
65 j += 2;
66 } else j++;
67 }
68
69 for (i = 0; i < m; i++) n_cmpl[i] = PetscSqrtReal(re[i] * re[i] + im[i] * im[i]);
70 PetscCall(KSPAGMRESLejafmaxarray(n_cmpl, 0, m, &pos));
71 j = 0;
72 if (im[pos] >= 0.0) {
73 rre[0] = re[pos];
74 rim[0] = im[pos];
75 j++;
76 spos[0] = pos;
77 }
78 while (j < (m)) {
79 if (im[pos] > 0) {
80 rre[j] = re[pos + 1];
81 rim[j] = im[pos + 1];
82 spos[j] = pos + 1;
83 j++;
84 }
85 PetscCall(KSPAGMRESLejaCfpdMax(re, im, spos, j, m, &pos));
86 if (im[pos] < 0) pos--;
87
88 if ((im[pos] >= 0) && (j < m)) {
89 rre[j] = re[pos];
90 rim[j] = im[pos];
91 spos[j] = pos;
92 j++;
93 }
94 }
95 PetscCall(PetscFree(spos));
96 PetscCall(PetscFree(n_cmpl));
97 PetscFunctionReturn(PETSC_SUCCESS);
98 }
99