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