genv.c

来自「math library from gnu」· C语言 代码 · 共 924 行 · 第 1/2 页

              scale = GSL_MIN(scale,
                        1.0 / (GSL_DBL_MIN *
                               GSL_MAX(1.0,
                                 GSL_MAX(fabs(acoef), fabs(bcoefr)))));
              if (lsa)
                acoef = ascale * (scale * sbeta);
              else
                acoef *= scale;

              if (lsb)
                bcoefr = bscale * (scale * salfar);
              else
                bcoefr *= scale;
            }

          acoefa = fabs(acoef);
          bcoefa = fabs(bcoefr);

          /* first component is 1 */
          gsl_vector_set(w->work3, je, 1.0);
          xmax = 1.0;

          /* compute contribution from column je of A and B to sum */

          for (i = 0; i < je; ++i)
            {
              gsl_vector_set(w->work3, i,
                bcoefr*gsl_matrix_get(T, i, je) -
                acoef * gsl_matrix_get(S, i, je));
            }
        }
      else
        {
          gsl_matrix_const_view vs =
            gsl_matrix_const_submatrix(S, je - 1, je - 1, 2, 2);
          gsl_matrix_const_view vt =
            gsl_matrix_const_submatrix(T, je - 1, je - 1, 2, 2);

          /* complex eigenvalue */

          gsl_schur_gen_eigvals(&vs.matrix,
                                &vt.matrix,
                                &bcoefr,
                                &temp2,
                                &bcoefi,
                                &acoef,
                                &temp);
          if (bcoefi == 0.0)
            {
              GSL_ERROR("gsl_schur_gen_eigvals failed on complex block", GSL_FAILURE);
            }

          /* scale to avoid over/underflow */
          acoefa = fabs(acoef);
          bcoefa = fabs(bcoefr) + fabs(bcoefi);
          scale = 1.0;

          if (acoefa*GSL_DBL_EPSILON < GSL_DBL_MIN && acoefa >= GSL_DBL_MIN)
            scale = (GSL_DBL_MIN / GSL_DBL_EPSILON) / acoefa;
          if (bcoefa*GSL_DBL_EPSILON < GSL_DBL_MIN && bcoefa >= GSL_DBL_MIN)
            scale = GSL_MAX(scale, (GSL_DBL_MIN/GSL_DBL_EPSILON) / bcoefa);
          if (GSL_DBL_MIN*acoefa > ascale)
            scale = ascale / (GSL_DBL_MIN * acoefa);
          if (GSL_DBL_MIN*bcoefa > bscale)
            scale = GSL_MIN(scale, bscale / (GSL_DBL_MIN*bcoefa));
          if (scale != 1.0)
            {
              acoef *= scale;
              acoefa = fabs(acoef);
              bcoefr *= scale;
              bcoefi *= scale;
              bcoefa = fabs(bcoefr) + fabs(bcoefi);
            }

          /* compute first two components of eigenvector */

          temp = acoef * gsl_matrix_get(S, je, je - 1);
          temp2r = acoef * gsl_matrix_get(S, je, je) -
                   bcoefr * gsl_matrix_get(T, je, je);
          temp2i = -bcoefi * gsl_matrix_get(T, je, je);

          if (fabs(temp) >= fabs(temp2r) + fabs(temp2i))
            {
              gsl_vector_set(w->work3, je, 1.0);
              gsl_vector_set(w->work4, je, 0.0);
              gsl_vector_set(w->work3, je - 1, -temp2r / temp);
              gsl_vector_set(w->work4, je - 1, -temp2i / temp);
            }
          else
            {
              gsl_vector_set(w->work3, je - 1, 1.0);
              gsl_vector_set(w->work4, je - 1, 0.0);
              temp = acoef * gsl_matrix_get(S, je - 1, je);
              gsl_vector_set(w->work3, je,
                (bcoefr*gsl_matrix_get(T, je - 1, je - 1) -
                 acoef*gsl_matrix_get(S, je - 1, je - 1)) / temp);
              gsl_vector_set(w->work4, je,
                bcoefi*gsl_matrix_get(T, je - 1, je - 1) / temp);
            }

          xmax = GSL_MAX(fabs(gsl_vector_get(w->work3, je)) +
                         fabs(gsl_vector_get(w->work4, je)),
                         fabs(gsl_vector_get(w->work3, je - 1)) +
                         fabs(gsl_vector_get(w->work4, je - 1)));

          /* compute contribution from column je and je - 1 */

          creala = acoef * gsl_vector_get(w->work3, je - 1);
          cimaga = acoef * gsl_vector_get(w->work4, je - 1);
          crealb = bcoefr * gsl_vector_get(w->work3, je - 1) -
                   bcoefi * gsl_vector_get(w->work4, je - 1);
          cimagb = bcoefi * gsl_vector_get(w->work3, je - 1) +
                   bcoefr * gsl_vector_get(w->work4, je - 1);
          cre2a = acoef * gsl_vector_get(w->work3, je);
          cim2a = acoef * gsl_vector_get(w->work4, je);
          cre2b = bcoefr * gsl_vector_get(w->work3, je) -
                  bcoefi * gsl_vector_get(w->work4, je);
          cim2b = bcoefi * gsl_vector_get(w->work3, je) +
                  bcoefr * gsl_vector_get(w->work4, je);

          for (i = 0; i < je - 1; ++i)
            {
              gsl_vector_set(w->work3, i,
                -creala * gsl_matrix_get(S, i, je - 1) +
                crealb * gsl_matrix_get(T, i, je - 1) -
                cre2a * gsl_matrix_get(S, i, je) +
                cre2b * gsl_matrix_get(T, i, je));
              gsl_vector_set(w->work4, i,
                -cimaga * gsl_matrix_get(S, i, je - 1) +
                cimagb * gsl_matrix_get(T, i, je - 1) -
                cim2a * gsl_matrix_get(S, i, je) +
                cim2b * gsl_matrix_get(T, i, je));
            }
        }

      dmin = GSL_MAX(GSL_DBL_MIN,
               GSL_MAX(GSL_DBL_EPSILON*acoefa*anorm,
                       GSL_DBL_EPSILON*bcoefa*bnorm));

      /* triangular solve of (a A - b B) x = 0 */

      il2by2 = 0;
      for (is = (int) je - (int) nw; is >= 0; --is)
        {
          j = (size_t) is;

          if (!il2by2 && j > 0)
            {
              if (gsl_matrix_get(S, j, j - 1) != 0.0)
                {
                  il2by2 = 1;
                  continue;
                }
            }

          bdiag[0] = gsl_matrix_get(T, j, j);
          if (il2by2)
            {
              na = 2;
              bdiag[1] = gsl_matrix_get(T, j + 1, j + 1);
            }
          else
            na = 1;


          if (nw == 1)
            {
              gsl_matrix_const_view sv =
                gsl_matrix_const_submatrix(S, j, j, na, na);
              gsl_vector_view xv, bv;

              bv = gsl_vector_subvector(w->work3, j, na);

              /*
               * the loop below expects the solution in the first column
               * of sum, so set stride to 2
               */
              xv = gsl_vector_view_array_with_stride(sum, 2, na);

              gsl_schur_solve_equation(acoef,
                                       &sv.matrix,
                                       bcoefr,
                                       bdiag[0],
                                       bdiag[1],
                                       &bv.vector,
                                       &xv.vector,
                                       &scale,
                                       &temp,
                                       dmin);
            }
          else
            {
              double bdat[4];
              gsl_matrix_const_view sv =
                gsl_matrix_const_submatrix(S, j, j, na, na);
              gsl_vector_complex_view xv =
                gsl_vector_complex_view_array(sum, na);
              gsl_vector_complex_view bv =
                gsl_vector_complex_view_array(bdat, na);
              gsl_complex z;

              bdat[0] = gsl_vector_get(w->work3, j);
              bdat[1] = gsl_vector_get(w->work4, j);
              if (na == 2)
                {
                  bdat[2] = gsl_vector_get(w->work3, j + 1);
                  bdat[3] = gsl_vector_get(w->work4, j + 1);
                }

              GSL_SET_COMPLEX(&z, bcoefr, bcoefi);

              gsl_schur_solve_equation_z(acoef,
                                         &sv.matrix,
                                         &z,
                                         bdiag[0],
                                         bdiag[1],
                                         &bv.vector,
                                         &xv.vector,
                                         &scale,
                                         &temp,
                                         dmin);
            }

          if (scale < 1.0)
            {
              for (jr = 0; jr <= je; ++jr)
                {
                  gsl_vector_set(w->work3, jr,
                    scale * gsl_vector_get(w->work3, jr));
                  if (nw == 2)
                    {
                      gsl_vector_set(w->work4, jr,
                        scale * gsl_vector_get(w->work4, jr));
                    }
                }
            }

          xmax = GSL_MAX(scale * xmax, temp);

          for (jr = 0; jr < na; ++jr)
            {
              gsl_vector_set(w->work3, j + jr, sum[jr*na]);
              if (nw == 2)
                gsl_vector_set(w->work4, j + jr, sum[jr*na + 1]);
            }

          if (j > 0)
            {
              xscale = 1.0 / GSL_MAX(1.0, xmax);
              temp = acoefa * gsl_vector_get(w->work1, j) +
                     bcoefa * gsl_vector_get(w->work2, j);
              if (il2by2)
                {
                  temp = GSL_MAX(temp,
                           acoefa * gsl_vector_get(w->work1, j + 1) +
                           bcoefa * gsl_vector_get(w->work2, j + 1));
                }

              temp = GSL_MAX(temp, GSL_MAX(acoefa, bcoefa));
              if (temp > bignum * xscale)
                {
                  for (jr = 0; jr <= je; ++jr)
                    {
                      gsl_vector_set(w->work3, jr,
                        xscale * gsl_vector_get(w->work3, jr));
                      if (nw == 2)
                        {
                          gsl_vector_set(w->work4, jr,
                            xscale * gsl_vector_get(w->work4, jr));
                        }
                    }
                  xmax *= xscale;
                }

              for (ja = 0; ja < na; ++ja)
                {
                  if (complex_pair)
                    {
                      creala = acoef * gsl_vector_get(w->work3, j + ja);
                      cimaga = acoef * gsl_vector_get(w->work4, j + ja);
                      crealb = bcoefr * gsl_vector_get(w->work3, j + ja) -
                               bcoefi * gsl_vector_get(w->work4, j + ja);
                      cimagb = bcoefi * gsl_vector_get(w->work3, j + ja) +
                               bcoefr * gsl_vector_get(w->work4, j + ja);
                      for (jr = 0; jr <= j - 1; ++jr)
                        {
                          gsl_vector_set(w->work3, jr,
                            gsl_vector_get(w->work3, jr) -
                            creala * gsl_matrix_get(S, jr, j + ja) +
                            crealb * gsl_matrix_get(T, jr, j + ja));
                          gsl_vector_set(w->work4, jr,
                            gsl_vector_get(w->work4, jr) -
                            cimaga * gsl_matrix_get(S, jr, j + ja) +
                            cimagb * gsl_matrix_get(T, jr, j + ja));
                        }
                    }
                  else
                    {
                      creala = acoef * gsl_vector_get(w->work3, j + ja);
                      crealb = bcoefr * gsl_vector_get(w->work3, j + ja);
                      for (jr = 0; jr <= j - 1; ++jr)
                        {
                          gsl_vector_set(w->work3, jr,
                            gsl_vector_get(w->work3, jr) -
                            creala * gsl_matrix_get(S, jr, j + ja) +
                            crealb * gsl_matrix_get(T, jr, j + ja));
                        }
                    } /* if (!complex_pair) */
                } /* for (ja = 0; ja < na; ++ja) */
            } /* if (j > 0) */

          il2by2 = 0;
        } /* for (i = 0; i < je - nw; ++i) */

      for (jr = 0; jr < N; ++jr)
        {
          gsl_vector_set(w->work5, jr,
            gsl_vector_get(w->work3, 0) * gsl_matrix_get(Z, jr, 0));
          if (nw == 2)
            {
              gsl_vector_set(w->work6, jr,
                gsl_vector_get(w->work4, 0) * gsl_matrix_get(Z, jr, 0));
            }
        }

      for (jc = 1; jc <= je; ++jc)
        {
          for (jr = 0; jr < N; ++jr)
            {
              gsl_vector_set(w->work5, jr,
                gsl_vector_get(w->work5, jr) +
                gsl_vector_get(w->work3, jc) * gsl_matrix_get(Z, jr, jc));
              if (nw == 2)
                {
                  gsl_vector_set(w->work6, jr,
                    gsl_vector_get(w->work6, jr) +
                    gsl_vector_get(w->work4, jc) * gsl_matrix_get(Z, jr, jc));
                }
            }
        }

      /* store the eigenvector */

      if (complex_pair)
        {
          ecol = gsl_matrix_complex_column(evec, je - 1);
          re = gsl_vector_complex_real(&ecol.vector);
          im = gsl_vector_complex_imag(&ecol.vector);

          ecol = gsl_matrix_complex_column(evec, je);
          re2 = gsl_vector_complex_real(&ecol.vector);
          im2 = gsl_vector_complex_imag(&ecol.vector);
        }
      else
        {
          ecol = gsl_matrix_complex_column(evec, je);
          re = gsl_vector_complex_real(&ecol.vector);
          im = gsl_vector_complex_imag(&ecol.vector);
        }

      for (jr = 0; jr < N; ++jr)
        {
          gsl_vector_set(&re.vector, jr, gsl_vector_get(w->work5, jr));
          if (complex_pair)
            {
              gsl_vector_set(&im.vector, jr, gsl_vector_get(w->work6, jr));
              gsl_vector_set(&re2.vector, jr, gsl_vector_get(w->work5, jr));
              gsl_vector_set(&im2.vector, jr, -gsl_vector_get(w->work6, jr));
            }
          else
            {
              gsl_vector_set(&im.vector, jr, 0.0);
            }
        }

      /* scale eigenvector */
      xmax = 0.0;
      if (complex_pair)
        {
          for (j = 0; j < N; ++j)
            {
              xmax = GSL_MAX(xmax,
                             fabs(gsl_vector_get(&re.vector, j)) +
                             fabs(gsl_vector_get(&im.vector, j)));
            }
        }
      else
        {
          for (j = 0; j < N; ++j)
            {
              xmax = GSL_MAX(xmax, fabs(gsl_vector_get(&re.vector, j)));
            }
        }

      if (xmax > GSL_DBL_MIN)
        {
          xscale = 1.0 / xmax;
          for (j = 0; j < N; ++j)
            {
              gsl_vector_set(&re.vector, j,
                             gsl_vector_get(&re.vector, j) * xscale);
              if (complex_pair)
                {
                  gsl_vector_set(&im.vector, j,
                                 gsl_vector_get(&im.vector, j) * xscale);
                  gsl_vector_set(&re2.vector, j,
                                 gsl_vector_get(&re2.vector, j) * xscale);
                  gsl_vector_set(&im2.vector, j,
                                 gsl_vector_get(&im2.vector, j) * xscale);
                }
            }
        }
    } /* for (k = 0; k < N; ++k) */

  return GSL_SUCCESS;
} /* genv_get_right_eigenvectors() */

/*
genv_normalize_eigenvectors()
  Normalize eigenvectors so that their Euclidean norm is 1

Inputs: alpha - eigenvalue numerators
        evec  - eigenvectors
*/

static void
genv_normalize_eigenvectors(gsl_vector_complex *alpha,
                            gsl_matrix_complex *evec)
{
  const size_t N = evec->size1;
  size_t i;     /* looping */
  gsl_complex ai;
  gsl_vector_complex_view vi;
  gsl_vector_view re, im;
  double scale; /* scaling factor */

  for (i = 0; i < N; ++i)
    {
      ai = gsl_vector_complex_get(alpha, i);
      vi = gsl_matrix_complex_column(evec, i);

      re = gsl_vector_complex_real(&vi.vector);

      if (GSL_IMAG(ai) == 0.0)
        {
          scale = 1.0 / gsl_blas_dnrm2(&re.vector);
          gsl_blas_dscal(scale, &re.vector);
        }
      else if (GSL_IMAG(ai) > 0.0)
        {
          im = gsl_vector_complex_imag(&vi.vector);

          scale = 1.0 / gsl_hypot(gsl_blas_dnrm2(&re.vector),
                                  gsl_blas_dnrm2(&im.vector));
          gsl_blas_zdscal(scale, &vi.vector);

          vi = gsl_matrix_complex_column(evec, i + 1);
          gsl_blas_zdscal(scale, &vi.vector);
        }
    }
} /* genv_normalize_eigenvectors() */