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