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: }