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: mx = re[0];
15: *pos = 0;
16: for (i = pt; i < n; i++) {
17: if (mx < re[i]) {
18: mx = re[i];
19: *pos = i;
20: }
21: }
22: return 0;
23: }
25: static PetscErrorCode KSPAGMRESLejaCfpdMax(PetscScalar *rm, PetscScalar *im, PetscInt *spos, PetscInt nbre, PetscInt n, PetscInt *rpos)
26: {
27: PetscScalar rd, id, pd, max;
28: PetscInt i, j;
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: return 0;
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: PetscMalloc1(m, &n_cmpl);
55: 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: }
69: for (i = 0; i < m; i++) n_cmpl[i] = PetscSqrtReal(re[i] * re[i] + im[i] * im[i]);
70: 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: KSPAGMRESLejaCfpdMax(re, im, spos, j, m, &pos);
86: if (im[pos] < 0) pos--;
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: PetscFree(spos);
96: PetscFree(n_cmpl);
97: return 0;
98: }