htswanalysis/MACS/lib/gsl/gsl-1.11/eigen/schur.c

   1 /* eigen/schur.c
   2  *
   3  * Copyright (C) 2006, 2007 Patrick Alken
   4  *
   5  * This program is free software; you can redistribute it and/or modify
   6  * it under the terms of the GNU General Public License as published by
   7  * the Free Software Foundation; either version 3 of the License, or (at
   8  * your option) any later version.
   9  *
  10  * This program is distributed in the hope that it will be useful, but
  11  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13  * General Public License for more details.
  14  *
  15  * You should have received a copy of the GNU General Public License
  16  * along with this program; if not, write to the Free Software
  17  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  18  */
  19
  20 #include <config.h>
  21 #include <gsl/gsl_eigen.h>
  22 #include <gsl/gsl_math.h>
  23 #include <gsl/gsl_matrix.h>
  24 #include <gsl/gsl_vector.h>
  25 #include <gsl/gsl_vector_complex.h>
  26 #include <gsl/gsl_blas.h>
  27 #include <gsl/gsl_complex.h>
  28 #include <gsl/gsl_complex_math.h>
  29
  30 /*
  31  * This module contains some routines related to manipulating the
  32  * Schur form of a matrix which are needed by the eigenvalue solvers
  33  *
  34  * This file contains routines based on original code from LAPACK
  35  * which is distributed under the modified BSD license.
  36  */
  37
  38 #define GSL_SCHUR_SMLNUM (2.0 * GSL_DBL_MIN)
  39 #define GSL_SCHUR_BIGNUM ((1.0 - GSL_DBL_EPSILON) / GSL_SCHUR_SMLNUM)
  40
  41 /*
  42 gsl_schur_gen_eigvals()
  43   Compute the eigenvalues of a 2-by-2 generalized block.
  44
  45 Inputs: A      - 2-by-2 matrix
  46         B      - 2-by-2 upper triangular matrix
  47         wr1    - (output) see notes
  48         wr2    - (output) see notes
  49         wi     - (output) see notes
  50         scale1 - (output) see notes
  51         scale2 - (output) see notes
  52
  53 Return: success
  54
  55 Notes:
  56
  57 1)
  58
  59 If the block contains real eigenvalues, then wi is set to 0,
  60 and wr1, wr2, scale1, and scale2 are set such that:
  61
  62 eval1 = wr1 * scale1
  63 eval2 = wr2 * scale2
  64
  65 If the block contains complex eigenvalues, then wr1, wr2, scale1,
  66 scale2, and wi are set such that:
  67
  68 wr1 = wr2 = scale1 * Re(eval)
  69 wi = scale1 * Im(eval)
  70
  71 wi is always non-negative
  72
  73 2) This routine is based on LAPACK DLAG2
  74 */
  75
  76 int
  77 gsl_schur_gen_eigvals(const gsl_matrix *A, const gsl_matrix *B, double *wr1,
  78                       double *wr2, double *wi, double *scale1,
  79                       double *scale2)
  80 {
  81   const double safemin = GSL_DBL_MIN * 1.0e2;
  82   const double safemax = 1.0 / safemin;
  83   const double rtmin = sqrt(safemin);
  84   const double rtmax = 1.0 / rtmin;
  85   double anorm, bnorm;
  86   double ascale, bscale, bsize;
  87   double s1, s2;
  88   double A11, A12, A21, A22;
  89   double B11, B12, B22;
  90   double binv11, binv22;
  91   double bmin;
  92   double as11, as12, as22, abi22;
  93   double pp, qq, shift, ss, discr, r;
  94
  95   /* scale A */
  96   anorm = GSL_MAX(GSL_MAX(fabs(gsl_matrix_get(A, 0, 0)) +
  97                           fabs(gsl_matrix_get(A, 1, 0)),
  98                           fabs(gsl_matrix_get(A, 0, 1)) +
  99                           fabs(gsl_matrix_get(A, 1, 1))),
 100                   safemin);
 101   ascale = 1.0 / anorm;
 102   A11 = ascale * gsl_matrix_get(A, 0, 0);
 103   A12 = ascale * gsl_matrix_get(A, 0, 1);
 104   A21 = ascale * gsl_matrix_get(A, 1, 0);
 105   A22 = ascale * gsl_matrix_get(A, 1, 1);
 106
 107   /* perturb B if necessary to ensure non-singularity */
 108   B11 = gsl_matrix_get(B, 0, 0);
 109   B12 = gsl_matrix_get(B, 0, 1);
 110   B22 = gsl_matrix_get(B, 1, 1);
 111   bmin = rtmin * GSL_MAX(fabs(B11),
 112                          GSL_MAX(fabs(B12), GSL_MAX(fabs(B22), rtmin)));
 113   if (fabs(B11) < bmin)
 114     B11 = GSL_SIGN(B11) * bmin;
 115   if (fabs(B22) < bmin)
 116     B22 = GSL_SIGN(B22) * bmin;
 117
 118   /* scale B */
 119   bnorm = GSL_MAX(fabs(B11), GSL_MAX(fabs(B12) + fabs(B22), safemin));
 120   bsize = GSL_MAX(fabs(B11), fabs(B22));
 121   bscale = 1.0 / bsize;
 122   B11 *= bscale;
 123   B12 *= bscale;
 124   B22 *= bscale;
 125
 126   /* compute larger eigenvalue */
 127
 128   binv11 = 1.0 / B11;
 129   binv22 = 1.0 / B22;
 130   s1 = A11 * binv11;
 131   s2 = A22 * binv22;
 132   if (fabs(s1) <= fabs(s2))
 133     {
 134       as12 = A12 - s1 * B12;
 135       as22 = A22 - s1 * B22;
 136       ss = A21 * (binv11 * binv22);
 137       abi22 = as22 * binv22 - ss * B12;
 138       pp = 0.5 * abi22;
 139       shift = s1;
 140     }
 141   else
 142     {
 143       as12 = A12 - s2 * B12;
 144       as11 = A11 - s2 * B11;
 145       ss = A21 * (binv11 * binv22);
 146       abi22 = -ss * B12;
 147       pp = 0.5 * (as11 * binv11 + abi22);
 148       shift = s2;
 149     }
 150
 151   qq = ss * as12;
 152   if (fabs(pp * rtmin) >= 1.0)
 153     {
 154       discr = (rtmin * pp) * (rtmin * pp) + qq * safemin;
 155       r = sqrt(fabs(discr)) * rtmax;
 156     }
 157   else if (pp * pp + fabs(qq) <= safemin)
 158     {
 159       discr = (rtmax * pp) * (rtmax * pp) + qq * safemax;
 160       r = sqrt(fabs(discr)) * rtmin;
 161     }
 162   else
 163     {
 164       discr = pp * pp + qq;
 165       r = sqrt(fabs(discr));
 166     }
 167
 168   if (discr >= 0.0 || r == 0.0)
 169     {
 170       double sum = pp + GSL_SIGN(pp) * r;
 171       double diff = pp - GSL_SIGN(pp) * r;
 172       double wbig = shift + sum;
 173       double wsmall = shift + diff;
 174
 175       /* compute smaller eigenvalue */
 176
 177       if (0.5 * fabs(wbig) > GSL_MAX(fabs(wsmall), safemin))
 178         {
 179           double wdet = (A11*A22 - A12*A21) * (binv11 * binv22);
 180           wsmall = wdet / wbig;
 181         }
 182
 183       /* choose (real) eigenvalue closest to 2,2 element of AB^{-1} for wr1 */
 184       if (pp > abi22)
 185         {
 186           *wr1 = GSL_MIN(wbig, wsmall);
 187           *wr2 = GSL_MAX(wbig, wsmall);
 188         }
 189       else
 190         {
 191           *wr1 = GSL_MAX(wbig, wsmall);
 192           *wr2 = GSL_MIN(wbig, wsmall);
 193         }
 194       *wi = 0.0;
 195     }
 196   else
 197     {
 198       /* complex eigenvalues */
 199       *wr1 = shift + pp;
 200       *wr2 = *wr1;
 201       *wi = r;
 202     }
 203
 204   /* compute scaling */
 205   {
 206     const double fuzzy1 = 1.0 + 1.0e-5;
 207     double c1, c2, c3, c4, c5;
 208     double wabs, wsize, wscale;
 209
 210     c1 = bsize * (safemin * GSL_MAX(1.0, ascale));
 211     c2 = safemin * GSL_MAX(1.0, bnorm);
 212     c3 = bsize * safemin;
 213     if (ascale <= 1.0 && bsize <= 1.0)
 214       c4 = GSL_MIN(1.0, (ascale / safemin) * bsize);
 215     else
 216       c4 = 1.0;
 217
 218     if (ascale <= 1.0 || bsize <= 1.0)
 219       c5 = GSL_MIN(1.0, ascale * bsize);
 220     else
 221       c5 = 1.0;
 222
 223     /* scale first eigenvalue */
 224     wabs = fabs(*wr1) + fabs(*wi);
 225     wsize = GSL_MAX(safemin,
 226               GSL_MAX(c1,
 227                 GSL_MAX(fuzzy1 * (wabs*c2 + c3),
 228                   GSL_MIN(c4, 0.5 * GSL_MAX(wabs, c5)))));
 229     if (wsize != 1.0)
 230       {
 231         wscale = 1.0 / wsize;
 232         if (wsize > 1.0)
 233           {
 234             *scale1 = (GSL_MAX(ascale, bsize) * wscale) *
 235                       GSL_MIN(ascale, bsize);
 236           }
 237         else
 238           {
 239             *scale1 = (GSL_MIN(ascale, bsize) * wscale) *
 240                       GSL_MAX(ascale, bsize);
 241           }
 242
 243         *wr1 *= wscale;
 244         if (*wi != 0.0)
 245           {
 246             *wi *= wscale;
 247             *wr2 = *wr1;
 248             *scale2 = *scale1;
 249           }
 250       }
 251     else
 252       {
 253         *scale1 = ascale * bsize;
 254         *scale2 = *scale1;
 255       }
 256
 257     /* scale second eigenvalue if real */
 258     if (*wi == 0.0)
 259       {
 260         wsize = GSL_MAX(safemin,
 261                   GSL_MAX(c1,
 262                     GSL_MAX(fuzzy1 * (fabs(*wr2) * c2 + c3),
 263                       GSL_MIN(c4, 0.5 * GSL_MAX(fabs(*wr2), c5)))));
 264         if (wsize != 1.0)
 265           {
 266             wscale = 1.0 / wsize;
 267             if (wsize > 1.0)
 268               {
 269                 *scale2 = (GSL_MAX(ascale, bsize) * wscale) *
 270                           GSL_MIN(ascale, bsize);
 271               }
 272             else
 273               {
 274                 *scale2 = (GSL_MIN(ascale, bsize) * wscale) *
 275                           GSL_MAX(ascale, bsize);
 276               }
 277
 278             *wr2 *= wscale;
 279           }
 280         else
 281           {
 282             *scale2 = ascale * bsize;
 283           }
 284       }
 285   }
 286
 287   return GSL_SUCCESS;
 288 } /* gsl_schur_gen_eigvals() */
 289
 290 /*
 291 gsl_schur_solve_equation()
 292
 293   Solve the equation which comes up in the back substitution
 294 when computing eigenvectors corresponding to real eigenvalues.
 295 The equation that is solved is:
 296
 297 (ca*A - z*D)*x = s*b
 298
 299 where
 300
 301 A is n-by-n with n = 1 or 2
 302 D is a n-by-n diagonal matrix
 303 b and x are n-by-1 real vectors
 304 s is a scaling factor set by this function to prevent overflow in x
 305
 306 Inputs: ca    - coefficient multiplying A
 307         A     - square matrix (n-by-n)
 308         z     - real scalar (eigenvalue)
 309         d1    - (1,1) element in diagonal matrix D
 310         d2    - (2,2) element in diagonal matrix D
 311         b     - right hand side vector
 312         x     - (output) where to store solution
 313         s     - (output) scale factor
 314         xnorm - (output) infinity norm of X
 315         smin  - lower bound on singular values of A - if ca*A - z*D
 316                 is less than this value, we'll use smin*I instead.
 317                 This value should be a safe distance above underflow.
 318
 319 Return: success
 320
 321 Notes: 1) A and b are not changed on output
 322        2) Based on lapack routine DLALN2
 323 */
 324
 325 int
 326 gsl_schur_solve_equation(double ca, const gsl_matrix *A, double z,
 327                          double d1, double d2, const gsl_vector *b,
 328                          gsl_vector *x, double *s, double *xnorm,
 329                          double smin)
 330 {
 331   size_t N = A->size1;
 332   double bnorm;
 333   double scale = 1.0;
 334
 335   if (N == 1)
 336     {
 337       double c,     /* denominator */
 338              cnorm; /* |c| */
 339
 340       /* we have a 1-by-1 (real) scalar system to solve */
 341
 342       c = ca * gsl_matrix_get(A, 0, 0) - z * d1;
 343       cnorm = fabs(c);
 344
 345       if (cnorm < smin)
 346         {
 347           /* set c = smin*I */
 348           c = smin;
 349           cnorm = smin;
 350         }
 351
 352       /* check scaling for x = b / c */
 353       bnorm = fabs(gsl_vector_get(b, 0));
 354       if (cnorm < 1.0 && bnorm > 1.0)
 355         {
 356           if (bnorm > GSL_SCHUR_BIGNUM*cnorm)
 357             scale = 1.0 / bnorm;
 358         }
 359
 360       /* compute x */
 361       gsl_vector_set(x, 0, gsl_vector_get(b, 0) * scale / c);
 362       *xnorm = fabs(gsl_vector_get(x, 0));
 363     } /* if (N == 1) */
 364   else
 365     {
 366       double cr[2][2];
 367       double *crv;
 368       double cmax;
 369       size_t icmax, j;
 370       double bval1, bval2;
 371       double ur11, ur12, ur22, ur11r;
 372       double cr21, cr22;
 373       double lr21;
 374       double b1, b2, bbnd;
 375       double x1, x2;
 376       double temp;
 377       size_t ipivot[4][4] = { { 0, 1, 2, 3 },
 378                               { 1, 0, 3, 2 },
 379                               { 2, 3, 0, 1 },
 380                               { 3, 2, 1, 0 } };
 381       int rswap[4] = { 0, 1, 0, 1 };
 382       int zswap[4] = { 0, 0, 1, 1 };
 383
 384       /*
 385        * we have a 2-by-2 real system to solve:
 386        *
 387        * ( ca * [ A11  A12 ] - z * [ D1 0  ] ) [ x1 ] = [ b1 ]
 388        * (      [ A21  A22 ]       [ 0  D2 ] ) [ x2 ]   [ b2 ]
 389        *
 390        * (z real)
 391        */
 392
 393       crv = (double *) cr;
 394
 395       /*
 396        * compute the real part of C = ca*A - z*D - use column ordering
 397        * here since porting from lapack
 398        */
 399       cr[0][0] = ca * gsl_matrix_get(A, 0, 0) - z * d1;
 400       cr[1][1] = ca * gsl_matrix_get(A, 1, 1) - z * d2;
 401       cr[0][1] = ca * gsl_matrix_get(A, 1, 0);
 402       cr[1][0] = ca * gsl_matrix_get(A, 0, 1);
 403
 404       /* find the largest element in C */
 405       cmax = 0.0;
 406       icmax = 0;
 407       for (j = 0; j < 4; ++j)
 408         {
 409           if (fabs(crv[j]) > cmax)
 410             {
 411               cmax = fabs(crv[j]);
 412               icmax = j;
 413             }
 414         }
 415
 416       bval1 = gsl_vector_get(b, 0);
 417       bval2 = gsl_vector_get(b, 1);
 418
 419       /* if norm(C) < smin, use smin*I */
 420
 421       if (cmax < smin)
 422         {
 423           bnorm = GSL_MAX(fabs(bval1), fabs(bval2));
 424           if (smin < 1.0 && bnorm > 1.0)
 425             {
 426               if (bnorm > GSL_SCHUR_BIGNUM*smin)
 427                 scale = 1.0 / bnorm;
 428             }
 429           temp = scale / smin;
 430           gsl_vector_set(x, 0, temp * bval1);
 431           gsl_vector_set(x, 1, temp * bval2);
 432           *xnorm = temp * bnorm;
 433           *s = scale;
 434           return GSL_SUCCESS;
 435         }
 436
 437       /* gaussian elimination with complete pivoting */
 438       ur11 = crv[icmax];
 439       cr21 = crv[ipivot[1][icmax]];
 440       ur12 = crv[ipivot[2][icmax]];
 441       cr22 = crv[ipivot[3][icmax]];
 442       ur11r = 1.0 / ur11;
 443       lr21 = ur11r * cr21;
 444       ur22 = cr22 - ur12 * lr21;
 445
 446       /* if smaller pivot < smin, use smin */
 447       if (fabs(ur22) < smin)
 448         ur22 = smin;
 449
 450       if (rswap[icmax])
 451         {
 452           b1 = bval2;
 453           b2 = bval1;
 454         }
 455       else
 456         {
 457           b1 = bval1;
 458           b2 = bval2;
 459         }
 460
 461       b2 -= lr21 * b1;
 462       bbnd = GSL_MAX(fabs(b1 * (ur22 * ur11r)), fabs(b2));
 463       if (bbnd > 1.0 && fabs(ur22) < 1.0)
 464         {
 465           if (bbnd >= GSL_SCHUR_BIGNUM * fabs(ur22))
 466             scale = 1.0 / bbnd;
 467         }
 468
 469       x2 = (b2 * scale) / ur22;
 470       x1 = (scale * b1) * ur11r - x2 * (ur11r * ur12);
 471       if (zswap[icmax])
 472         {
 473           gsl_vector_set(x, 0, x2);
 474           gsl_vector_set(x, 1, x1);
 475         }
 476       else
 477         {
 478           gsl_vector_set(x, 0, x1);
 479           gsl_vector_set(x, 1, x2);
 480         }
 481
 482       *xnorm = GSL_MAX(fabs(x1), fabs(x2));
 483
 484       /* further scaling if norm(A) norm(X) > overflow */
 485       if (*xnorm > 1.0 && cmax > 1.0)
 486         {
 487           if (*xnorm > GSL_SCHUR_BIGNUM / cmax)
 488             {
 489               temp = cmax / GSL_SCHUR_BIGNUM;
 490               gsl_blas_dscal(temp, x);
 491               *xnorm *= temp;
 492               scale *= temp;
 493             }
 494         }
 495     } /* if (N == 2) */
 496
 497   *s = scale;
 498   return GSL_SUCCESS;
 499 } /* gsl_schur_solve_equation() */
 500
 501 /*
 502 gsl_schur_solve_equation_z()
 503
 504   Solve the equation which comes up in the back substitution
 505 when computing eigenvectors corresponding to complex eigenvalues.
 506 The equation that is solved is:
 507
 508 (ca*A - z*D)*x = s*b
 509
 510 where
 511
 512 A is n-by-n with n = 1 or 2
 513 D is a n-by-n diagonal matrix
 514 b and x are n-by-1 complex vectors
 515 s is a scaling factor set by this function to prevent overflow in x
 516
 517 Inputs: ca    - coefficient multiplying A
 518         A     - square matrix (n-by-n)
 519         z     - complex scalar (eigenvalue)
 520         d1    - (1,1) element in diagonal matrix D
 521         d2    - (2,2) element in diagonal matrix D
 522         b     - right hand side vector
 523         x     - (output) where to store solution
 524         s     - (output) scale factor
 525         xnorm - (output) infinity norm of X
 526         smin  - lower bound on singular values of A - if ca*A - z*D
 527                 is less than this value, we'll use smin*I instead.
 528                 This value should be a safe distance above underflow.
 529
 530 Notes: 1) A and b are not changed on output
 531        2) Based on lapack routine DLALN2
 532 */
 533
 534 int
 535 gsl_schur_solve_equation_z(double ca, const gsl_matrix *A, gsl_complex *z,
 536                            double d1, double d2,
 537                            const gsl_vector_complex *b,
 538                            gsl_vector_complex *x, double *s, double *xnorm,
 539                            double smin)
 540 {
 541   size_t N = A->size1;
 542   double scale = 1.0;
 543   double bnorm;
 544
 545   if (N == 1)
 546     {
 547       double cr,    /* denominator */
 548              ci,
 549              cnorm; /* |c| */
 550       gsl_complex bval, c, xval, tmp;
 551
 552       /* we have a 1-by-1 (complex) scalar system to solve */
 553
 554       /* c = ca*a - z*d1 */
 555       cr = ca * gsl_matrix_get(A, 0, 0) - GSL_REAL(*z) * d1;
 556       ci = -GSL_IMAG(*z) * d1;
 557       cnorm = fabs(cr) + fabs(ci);
 558
 559       if (cnorm < smin)
 560         {
 561           /* set c = smin*I */
 562           cr = smin;
 563           ci = 0.0;
 564           cnorm = smin;
 565         }
 566
 567       /* check scaling for x = b / c */
 568       bval = gsl_vector_complex_get(b, 0);
 569       bnorm = fabs(GSL_REAL(bval)) + fabs(GSL_IMAG(bval));
 570       if (cnorm < 1.0 && bnorm > 1.0)
 571         {
 572           if (bnorm > GSL_SCHUR_BIGNUM*cnorm)
 573             scale = 1.0 / bnorm;
 574         }
 575
 576       /* compute x */
 577       GSL_SET_COMPLEX(&tmp, scale*GSL_REAL(bval), scale*GSL_IMAG(bval));
 578       GSL_SET_COMPLEX(&c, cr, ci);
 579       xval = gsl_complex_div(tmp, c);
 580
 581       gsl_vector_complex_set(x, 0, xval);
 582
 583       *xnorm = fabs(GSL_REAL(xval)) + fabs(GSL_IMAG(xval));
 584     } /* if (N == 1) */
 585   else
 586     {
 587       double cr[2][2], ci[2][2];
 588       double *civ, *crv;
 589       double cmax;
 590       gsl_complex bval1, bval2;
 591       gsl_complex xval1, xval2;
 592       double xr1, xi1;
 593       size_t icmax;
 594       size_t j;
 595       double temp;
 596       double ur11, ur12, ur22, ui11, ui12, ui22, ur11r, ui11r;
 597       double ur12s, ui12s;
 598       double u22abs;
 599       double lr21, li21;
 600       double cr21, cr22, ci21, ci22;
 601       double br1, bi1, br2, bi2, bbnd;
 602       gsl_complex b1, b2;
 603       size_t ipivot[4][4] = { { 0, 1, 2, 3 },
 604                               { 1, 0, 3, 2 },
 605                               { 2, 3, 0, 1 },
 606                               { 3, 2, 1, 0 } };
 607       int rswap[4] = { 0, 1, 0, 1 };
 608       int zswap[4] = { 0, 0, 1, 1 };
 609
 610       /*
 611        * complex 2-by-2 system:
 612        *
 613        * ( ca * [ A11 A12 ] - z * [ D1 0 ] ) [ X1 ] = [ B1 ]
 614        * (      [ A21 A22 ]       [ 0  D2] ) [ X2 ]   [ B2 ]
 615        *
 616        * (z complex)
 617        *
 618        * where the X and B values are complex.
 619        */
 620
 621       civ = (double *) ci;
 622       crv = (double *) cr;
 623
 624       /*
 625        * compute the real part of C = ca*A - z*D - use column ordering
 626        * here since porting from lapack
 627        */
 628       cr[0][0] = ca*gsl_matrix_get(A, 0, 0) - GSL_REAL(*z)*d1;
 629       cr[1][1] = ca*gsl_matrix_get(A, 1, 1) - GSL_REAL(*z)*d2;
 630       cr[0][1] = ca*gsl_matrix_get(A, 1, 0);
 631       cr[1][0] = ca*gsl_matrix_get(A, 0, 1);
 632
 633       /* compute the imaginary part */
 634       ci[0][0] = -GSL_IMAG(*z) * d1;
 635       ci[0][1] = 0.0;
 636       ci[1][0] = 0.0;
 637       ci[1][1] = -GSL_IMAG(*z) * d2;
 638
 639       cmax = 0.0;
 640       icmax = 0;
 641
 642       for (j = 0; j < 4; ++j)
 643         {
 644           if (fabs(crv[j]) + fabs(civ[j]) > cmax)
 645             {
 646               cmax = fabs(crv[j]) + fabs(civ[j]);
 647               icmax = j;
 648             }
 649         }
 650
 651       bval1 = gsl_vector_complex_get(b, 0);
 652       bval2 = gsl_vector_complex_get(b, 1);
 653
 654       /* if norm(C) < smin, use smin*I */
 655       if (cmax < smin)
 656         {
 657           bnorm = GSL_MAX(fabs(GSL_REAL(bval1)) + fabs(GSL_IMAG(bval1)),
 658                           fabs(GSL_REAL(bval2)) + fabs(GSL_IMAG(bval2)));
 659           if (smin < 1.0 && bnorm > 1.0)
 660             {
 661               if (bnorm > GSL_SCHUR_BIGNUM*smin)
 662                 scale = 1.0 / bnorm;
 663             }
 664
 665           temp = scale / smin;
 666           xval1 = gsl_complex_mul_real(bval1, temp);
 667           xval2 = gsl_complex_mul_real(bval2, temp);
 668           gsl_vector_complex_set(x, 0, xval1);
 669           gsl_vector_complex_set(x, 1, xval2);
 670           *xnorm = temp * bnorm;
 671           *s = scale;
 672           return GSL_SUCCESS;
 673         }
 674
 675       /* gaussian elimination with complete pivoting */
 676       ur11 = crv[icmax];
 677       ui11 = civ[icmax];
 678       cr21 = crv[ipivot[1][icmax]];
 679       ci21 = civ[ipivot[1][icmax]];
 680       ur12 = crv[ipivot[2][icmax]];
 681       ui12 = civ[ipivot[2][icmax]];
 682       cr22 = crv[ipivot[3][icmax]];
 683       ci22 = civ[ipivot[3][icmax]];
 684
 685       if (icmax == 0 || icmax == 3)
 686         {
 687           /* off diagonals of pivoted C are real */
 688           if (fabs(ur11) > fabs(ui11))
 689             {
 690               temp = ui11 / ur11;
 691               ur11r = 1.0 / (ur11 * (1.0 + temp*temp));
 692               ui11r = -temp * ur11r;
 693             }
 694           else
 695             {
 696               temp = ur11 / ui11;
 697               ui11r = -1.0 / (ui11 * (1.0 + temp*temp));
 698               ur11r = -temp*ui11r;
 699             }
 700           lr21 = cr21 * ur11r;
 701           li21 = cr21 * ui11r;
 702           ur12s = ur12 * ur11r;
 703           ui12s = ur12 * ui11r;
 704           ur22 = cr22 - ur12 * lr21;
 705           ui22 = ci22 - ur12 * li21;
 706         }
 707       else
 708         {
 709           /* diagonals of pivoted C are real */
 710           ur11r = 1.0 / ur11;
 711           ui11r = 0.0;
 712           lr21 = cr21 * ur11r;
 713           li21 = ci21 * ur11r;
 714           ur12s = ur12 * ur11r;
 715           ui12s = ui12 * ur11r;
 716           ur22 = cr22 - ur12 * lr21 + ui12 * li21;
 717           ui22 = -ur12 * li21 - ui12 * lr21;
 718         }
 719
 720       u22abs = fabs(ur22) + fabs(ui22);
 721
 722       /* if smaller pivot < smin, use smin */
 723       if (u22abs < smin)
 724         {
 725           ur22 = smin;
 726           ui22 = 0.0;
 727         }
 728
 729       if (rswap[icmax])
 730         {
 731           br2 = GSL_REAL(bval1);
 732           bi2 = GSL_IMAG(bval1);
 733           br1 = GSL_REAL(bval2);
 734           bi1 = GSL_IMAG(bval2);
 735         }
 736       else
 737         {
 738           br1 = GSL_REAL(bval1);
 739           bi1 = GSL_IMAG(bval1);
 740           br2 = GSL_REAL(bval2);
 741           bi2 = GSL_IMAG(bval2);
 742         }
 743
 744       br2 += li21*bi1 - lr21*br1;
 745       bi2 -= li21*br1 + lr21*bi1;
 746       bbnd = GSL_MAX((fabs(br1) + fabs(bi1)) *
 747                      (u22abs * (fabs(ur11r) + fabs(ui11r))),
 748                      fabs(br2) + fabs(bi2));
 749       if (bbnd > 1.0 && u22abs < 1.0)
 750         {
 751           if (bbnd >= GSL_SCHUR_BIGNUM*u22abs)
 752             {
 753               scale = 1.0 / bbnd;
 754               br1 *= scale;
 755               bi1 *= scale;
 756               br2 *= scale;
 757               bi2 *= scale;
 758             }
 759         }
 760
 761       GSL_SET_COMPLEX(&b1, br2, bi2);
 762       GSL_SET_COMPLEX(&b2, ur22, ui22);
 763       xval2 = gsl_complex_div(b1, b2);
 764
 765       xr1 = ur11r*br1 - ui11r*bi1 - ur12s*GSL_REAL(xval2) + ui12s*GSL_IMAG(xval2);
 766       xi1 = ui11r*br1 + ur11r*bi1 - ui12s*GSL_REAL(xval2) - ur12s*GSL_IMAG(xval2);
 767       GSL_SET_COMPLEX(&xval1, xr1, xi1);
 768
 769       if (zswap[icmax])
 770         {
 771           gsl_vector_complex_set(x, 0, xval2);
 772           gsl_vector_complex_set(x, 1, xval1);
 773         }
 774       else
 775         {
 776           gsl_vector_complex_set(x, 0, xval1);
 777           gsl_vector_complex_set(x, 1, xval2);
 778         }
 779
 780       *xnorm = GSL_MAX(fabs(GSL_REAL(xval1)) + fabs(GSL_IMAG(xval1)),
 781                        fabs(GSL_REAL(xval2)) + fabs(GSL_IMAG(xval2)));
 782
 783       /* further scaling if norm(A) norm(X) > overflow */
 784       if (*xnorm > 1.0 && cmax > 1.0)
 785         {
 786           if (*xnorm > GSL_SCHUR_BIGNUM / cmax)
 787             {
 788               temp = cmax / GSL_SCHUR_BIGNUM;
 789               gsl_blas_zdscal(temp, x);
 790               *xnorm *= temp;
 791               scale *= temp;
 792             }
 793         }
 794     } /* if (N == 2) */
 795
 796   *s = scale;
 797   return GSL_SUCCESS;
 798 } /* gsl_schur_solve_equation_z() */