schur.c

math library from gnu

      cr[0][1] = ca * gsl_matrix_get(A, 1, 0);
      cr[1][0] = ca * gsl_matrix_get(A, 0, 1);

      /* find the largest element in C */
      cmax = 0.0;
      icmax = 0;
      for (j = 0; j < 4; ++j)
        {
          if (fabs(crv[j]) > cmax)
            {
              cmax = fabs(crv[j]);
              icmax = j;
            }
        }

      bval1 = gsl_vector_get(b, 0);
      bval2 = gsl_vector_get(b, 1);

      /* if norm(C) < smin, use smin*I */

      if (cmax < smin)
        {
          bnorm = GSL_MAX(fabs(bval1), fabs(bval2));
          if (smin < 1.0 && bnorm > 1.0)
            {
              if (bnorm > GSL_SCHUR_BIGNUM*smin)
                scale = 1.0 / bnorm;
            }
          temp = scale / smin;
          gsl_vector_set(x, 0, temp * bval1);
          gsl_vector_set(x, 1, temp * bval2);
          *xnorm = temp * bnorm;
          *s = scale;
          return GSL_SUCCESS;
        }

      /* gaussian elimination with complete pivoting */
      ur11 = crv[icmax];
      cr21 = crv[ipivot[1][icmax]];
      ur12 = crv[ipivot[2][icmax]];
      cr22 = crv[ipivot[3][icmax]];
      ur11r = 1.0 / ur11;
      lr21 = ur11r * cr21;
      ur22 = cr22 - ur12 * lr21;

      /* if smaller pivot < smin, use smin */
      if (fabs(ur22) < smin)
        ur22 = smin;

      if (rswap[icmax])
        {
          b1 = bval2;
          b2 = bval1;
        }
      else
        {
          b1 = bval1;
          b2 = bval2;
        }

      b2 -= lr21 * b1;
      bbnd = GSL_MAX(fabs(b1 * (ur22 * ur11r)), fabs(b2));
      if (bbnd > 1.0 && fabs(ur22) < 1.0)
        {
          if (bbnd >= GSL_SCHUR_BIGNUM * fabs(ur22))
            scale = 1.0 / bbnd;
        }

      x2 = (b2 * scale) / ur22;
      x1 = (scale * b1) * ur11r - x2 * (ur11r * ur12);
      if (zswap[icmax])
        {
          gsl_vector_set(x, 0, x2);
          gsl_vector_set(x, 1, x1);
        }
      else
        {
          gsl_vector_set(x, 0, x1);
          gsl_vector_set(x, 1, x2);
        }

      *xnorm = GSL_MAX(fabs(x1), fabs(x2));

      /* further scaling if norm(A) norm(X) > overflow */
      if (*xnorm > 1.0 && cmax > 1.0)
        {
          if (*xnorm > GSL_SCHUR_BIGNUM / cmax)
            {
              temp = cmax / GSL_SCHUR_BIGNUM;
              gsl_blas_dscal(temp, x);
              *xnorm *= temp;
              scale *= temp;
            }
        }
    } /* if (N == 2) */

  *s = scale;
  return GSL_SUCCESS;
} /* gsl_schur_solve_equation() */

/*
gsl_schur_solve_equation_z()

  Solve the equation which comes up in the back substitution
when computing eigenvectors corresponding to complex eigenvalues.
The equation that is solved is:

(ca*A - z*D)*x = s*b

where

A is n-by-n with n = 1 or 2
D is a n-by-n diagonal matrix
b and x are n-by-1 complex vectors
s is a scaling factor set by this function to prevent overflow in x

Inputs: ca    - coefficient multiplying A
        A     - square matrix (n-by-n)
        z     - complex scalar (eigenvalue)
        d1    - (1,1) element in diagonal matrix D
        d2    - (2,2) element in diagonal matrix D
        b     - right hand side vector
        x     - (output) where to store solution
        s     - (output) scale factor
        xnorm - (output) infinity norm of X
        smin  - lower bound on singular values of A - if ca*A - z*D
                is less than this value, we'll use smin*I instead.
                This value should be a safe distance above underflow.

Notes: 1) A and b are not changed on output
       2) Based on lapack routine DLALN2
*/

int
gsl_schur_solve_equation_z(double ca, const gsl_matrix *A, gsl_complex *z,
                           double d1, double d2,
                           const gsl_vector_complex *b,
                           gsl_vector_complex *x, double *s, double *xnorm,
                           double smin)
{
  size_t N = A->size1;
  double scale = 1.0;
  double bnorm;

  if (N == 1)
    {
      double cr,    /* denominator */
             ci,
             cnorm; /* |c| */
      gsl_complex bval, c, xval, tmp;

      /* we have a 1-by-1 (complex) scalar system to solve */

      /* c = ca*a - z*d1 */
      cr = ca * gsl_matrix_get(A, 0, 0) - GSL_REAL(*z) * d1;
      ci = -GSL_IMAG(*z) * d1;
      cnorm = fabs(cr) + fabs(ci);

      if (cnorm < smin)
        {
          /* set c = smin*I */
          cr = smin;
          ci = 0.0;
          cnorm = smin;
        }

      /* check scaling for x = b / c */
      bval = gsl_vector_complex_get(b, 0);
      bnorm = fabs(GSL_REAL(bval)) + fabs(GSL_IMAG(bval));
      if (cnorm < 1.0 && bnorm > 1.0)
        {
          if (bnorm > GSL_SCHUR_BIGNUM*cnorm)
            scale = 1.0 / bnorm;
        }

      /* compute x */
      GSL_SET_COMPLEX(&tmp, scale*GSL_REAL(bval), scale*GSL_IMAG(bval));
      GSL_SET_COMPLEX(&c, cr, ci);
      xval = gsl_complex_div(tmp, c);

      gsl_vector_complex_set(x, 0, xval);

      *xnorm = fabs(GSL_REAL(xval)) + fabs(GSL_IMAG(xval));
    } /* if (N == 1) */
  else
    {
      double cr[2][2], ci[2][2];
      double *civ, *crv;
      double cmax;
      gsl_complex bval1, bval2;
      gsl_complex xval1, xval2;
      double xr1, xi1;
      size_t icmax;
      size_t j;
      double temp;
      double ur11, ur12, ur22, ui11, ui12, ui22, ur11r, ui11r;
      double ur12s, ui12s;
      double u22abs;
      double lr21, li21;
      double cr21, cr22, ci21, ci22;
      double br1, bi1, br2, bi2, bbnd;
      gsl_complex b1, b2;
      size_t ipivot[4][4] = { { 0, 1, 2, 3 },
                              { 1, 0, 3, 2 },
                              { 2, 3, 0, 1 },
                              { 3, 2, 1, 0 } };
      int rswap[4] = { 0, 1, 0, 1 };
      int zswap[4] = { 0, 0, 1, 1 };

      /*
       * complex 2-by-2 system:
       *
       * ( ca * [ A11 A12 ] - z * [ D1 0 ] ) [ X1 ] = [ B1 ]
       * (      [ A21 A22 ]       [ 0  D2] ) [ X2 ]   [ B2 ]
       *
       * (z complex)
       *
       * where the X and B values are complex.
       */

      civ = (double *) ci;
      crv = (double *) cr;

      /*
       * compute the real part of C = ca*A - z*D - use column ordering
       * here since porting from lapack
       */
      cr[0][0] = ca*gsl_matrix_get(A, 0, 0) - GSL_REAL(*z)*d1;
      cr[1][1] = ca*gsl_matrix_get(A, 1, 1) - GSL_REAL(*z)*d2;
      cr[0][1] = ca*gsl_matrix_get(A, 1, 0);
      cr[1][0] = ca*gsl_matrix_get(A, 0, 1);

      /* compute the imaginary part */
      ci[0][0] = -GSL_IMAG(*z) * d1;
      ci[0][1] = 0.0;
      ci[1][0] = 0.0;
      ci[1][1] = -GSL_IMAG(*z) * d2;

      cmax = 0.0;
      icmax = 0;

      for (j = 0; j < 4; ++j)
        {
          if (fabs(crv[j]) + fabs(civ[j]) > cmax)
            {
              cmax = fabs(crv[j]) + fabs(civ[j]);
              icmax = j;
            }
        }

      bval1 = gsl_vector_complex_get(b, 0);
      bval2 = gsl_vector_complex_get(b, 1);

      /* if norm(C) < smin, use smin*I */
      if (cmax < smin)
        {
          bnorm = GSL_MAX(fabs(GSL_REAL(bval1)) + fabs(GSL_IMAG(bval1)),
                          fabs(GSL_REAL(bval2)) + fabs(GSL_IMAG(bval2)));
          if (smin < 1.0 && bnorm > 1.0)
            {
              if (bnorm > GSL_SCHUR_BIGNUM*smin)
                scale = 1.0 / bnorm;
            }

          temp = scale / smin;
          xval1 = gsl_complex_mul_real(bval1, temp);
          xval2 = gsl_complex_mul_real(bval2, temp);
          gsl_vector_complex_set(x, 0, xval1);
          gsl_vector_complex_set(x, 1, xval2);
          *xnorm = temp * bnorm;
          *s = scale;
          return GSL_SUCCESS;
        }

      /* gaussian elimination with complete pivoting */
      ur11 = crv[icmax];
      ui11 = civ[icmax];
      cr21 = crv[ipivot[1][icmax]];
      ci21 = civ[ipivot[1][icmax]];
      ur12 = crv[ipivot[2][icmax]];
      ui12 = civ[ipivot[2][icmax]];
      cr22 = crv[ipivot[3][icmax]];
      ci22 = civ[ipivot[3][icmax]];

      if (icmax == 0 || icmax == 3)
        {
          /* off diagonals of pivoted C are real */
          if (fabs(ur11) > fabs(ui11))
            {
              temp = ui11 / ur11;
              ur11r = 1.0 / (ur11 * (1.0 + temp*temp));
              ui11r = -temp * ur11r;
            }
          else
            {
              temp = ur11 / ui11;
              ui11r = -1.0 / (ui11 * (1.0 + temp*temp));
              ur11r = -temp*ui11r;
            }
          lr21 = cr21 * ur11r;
          li21 = cr21 * ui11r;
          ur12s = ur12 * ur11r;
          ui12s = ur12 * ui11r;
          ur22 = cr22 - ur12 * lr21;
          ui22 = ci22 - ur12 * li21;
        }
      else
        {
          /* diagonals of pivoted C are real */
          ur11r = 1.0 / ur11;
          ui11r = 0.0;
          lr21 = cr21 * ur11r;
          li21 = ci21 * ur11r;
          ur12s = ur12 * ur11r;
          ui12s = ui12 * ur11r;
          ur22 = cr22 - ur12 * lr21 + ui12 * li21;
          ui22 = -ur12 * li21 - ui12 * lr21;
        }

      u22abs = fabs(ur22) + fabs(ui22);

      /* if smaller pivot < smin, use smin */
      if (u22abs < smin)
        {
          ur22 = smin;
          ui22 = 0.0;
        }

      if (rswap[icmax])
        {
          br2 = GSL_REAL(bval1);
          bi2 = GSL_IMAG(bval1);
          br1 = GSL_REAL(bval2);
          bi1 = GSL_IMAG(bval2);
        }
      else
        {
          br1 = GSL_REAL(bval1);
          bi1 = GSL_IMAG(bval1);
          br2 = GSL_REAL(bval2);
          bi2 = GSL_IMAG(bval2);
        }

      br2 += li21*bi1 - lr21*br1;
      bi2 -= li21*br1 + lr21*bi1;
      bbnd = GSL_MAX((fabs(br1) + fabs(bi1)) *
                     (u22abs * (fabs(ur11r) + fabs(ui11r))),
                     fabs(br2) + fabs(bi2));
      if (bbnd > 1.0 && u22abs < 1.0)
        {
          if (bbnd >= GSL_SCHUR_BIGNUM*u22abs)
            {
              scale = 1.0 / bbnd;
              br1 *= scale;
              bi1 *= scale;
              br2 *= scale;
              bi2 *= scale;
            }
        }

      GSL_SET_COMPLEX(&b1, br2, bi2);
      GSL_SET_COMPLEX(&b2, ur22, ui22);
      xval2 = gsl_complex_div(b1, b2);

      xr1 = ur11r*br1 - ui11r*bi1 - ur12s*GSL_REAL(xval2) + ui12s*GSL_IMAG(xval2);
      xi1 = ui11r*br1 + ur11r*bi1 - ui12s*GSL_REAL(xval2) - ur12s*GSL_IMAG(xval2);
      GSL_SET_COMPLEX(&xval1, xr1, xi1);

      if (zswap[icmax])
        {
          gsl_vector_complex_set(x, 0, xval2);
          gsl_vector_complex_set(x, 1, xval1);
        }
      else
        {
          gsl_vector_complex_set(x, 0, xval1);
          gsl_vector_complex_set(x, 1, xval2);
        }

      *xnorm = GSL_MAX(fabs(GSL_REAL(xval1)) + fabs(GSL_IMAG(xval1)),
                       fabs(GSL_REAL(xval2)) + fabs(GSL_IMAG(xval2)));

      /* further scaling if norm(A) norm(X) > overflow */
      if (*xnorm > 1.0 && cmax > 1.0)
        {
          if (*xnorm > GSL_SCHUR_BIGNUM / cmax)
            {
              temp = cmax / GSL_SCHUR_BIGNUM;
              gsl_blas_zdscal(temp, x);
              *xnorm *= temp;
              scale *= temp;
            }
        }
    } /* if (N == 2) */

  *s = scale;
  return GSL_SUCCESS;
} /* gsl_schur_solve_equation_z() */