1 #define PETSCKSP_DLL
2 /*
3 Functions in this file reorder the Ritz values in the (modified) Leja order.
4
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>
8
KSPAGMRESLejafmaxarray(PetscScalar * re,PetscInt pt,PetscInt n,PetscInt * pos)9 static PetscErrorCode KSPAGMRESLejafmaxarray(PetscScalar *re, PetscInt pt, PetscInt n,PetscInt *pos)
10 {
11 PetscInt i;
12 PetscScalar mx;
13
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(0);
24 }
25
KSPAGMRESLejaCfpdMax(PetscScalar * rm,PetscScalar * im,PetscInt * spos,PetscInt nbre,PetscInt n,PetscInt * rpos)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;
30
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(0);
48 }
49
KSPAGMRESLejaOrdering(PetscScalar * re,PetscScalar * im,PetscScalar * rre,PetscScalar * rim,PetscInt m)50 PetscErrorCode KSPAGMRESLejaOrdering(PetscScalar *re, PetscScalar *im, PetscScalar *rre, PetscScalar *rim, PetscInt m)
51 {
52 PetscInt *spos;
53 PetscScalar *n_cmpl,temp;
54 PetscErrorCode ierr;
55 PetscInt i, pos, j;
56
57 PetscFunctionBegin;
58 ierr = PetscMalloc1(m, &n_cmpl);CHKERRQ(ierr);
59 ierr = PetscMalloc1(m, &spos);CHKERRQ(ierr);
60 /* Check the proper order of complex conjugate pairs */
61 j = 0;
62 while (j < m) {
63 if (im[j] != 0.0) { /* complex eigenvalue */
64 if (im[j] < 0.0) { /* change the order */
65 temp = im[j+1];
66 im[j+1] = im[j];
67 im[j] = temp;
68 }
69 j += 2;
70 } else j++;
71 }
72
73 for (i = 0; i < m; i++) n_cmpl[i] = PetscSqrtReal(re[i]*re[i]+im[i]*im[i]);
74 ierr = KSPAGMRESLejafmaxarray(n_cmpl, 0, m, &pos);CHKERRQ(ierr);
75 j = 0;
76 if (im[pos] >= 0.0) {
77 rre[0] = re[pos];
78 rim[0] = im[pos];
79 j++;
80 spos[0] = pos;
81 }
82 while (j < (m)) {
83 if (im[pos] > 0) {
84 rre[j] = re[pos+1];
85 rim[j] = im[pos+1];
86 spos[j] = pos + 1;
87 j++;
88 }
89 ierr = KSPAGMRESLejaCfpdMax(re, im, spos, j, m, &pos);CHKERRQ(ierr);
90 if (im[pos] < 0) pos--;
91
92 if ((im[pos] >= 0) && (j < m)) {
93 rre[j] = re[pos];
94 rim[j] = im[pos];
95 spos[j] = pos;
96 j++;
97 }
98 }
99 ierr = PetscFree(spos);CHKERRQ(ierr);
100 ierr = PetscFree(n_cmpl);CHKERRQ(ierr);
101 PetscFunctionReturn(0);
102 }
103