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