htswanalysis/MACS/lib/gsl/gsl-1.11/cblas/source_rotmg.h

   1 /* blas/source_rotmg.h
   2  *
   3  * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Brian Gough
   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 {
  21   const BASE G = 4096.0, G2 = G * G;
  22   BASE D1 = *d1, D2 = *d2, x = *b1, y = b2;
  23   BASE h11, h12, h21, h22, u;
  24
  25   BASE c, s;
  26
  27   /* case of d1 < 0, appendix A, second to last paragraph */
  28
  29   if (D1 < 0.0) {
  30     P[0] = -1;
  31     P[1] = 0;
  32     P[2] = 0;
  33     P[3] = 0;
  34     P[4] = 0;
  35     *d1 = 0;
  36     *d2 = 0;
  37     *b1 = 0;
  38     return;
  39   }
  40
  41   if (D2 * y == 0.0) {
  42     P[0] = -2;                  /* case of H = I, page 315 */
  43     return;
  44   }
  45
  46   c = fabs(D1 * x * x);
  47   s = fabs(D2 * y * y);
  48
  49   if (c > s) {
  50     /* case of equation A6 */
  51
  52     P[0] = 0.0;
  53
  54     h11 = 1;
  55     h12 = (D2 * y) / (D1 * x);
  56     h21 = -y / x;
  57     h22 = 1;
  58
  59     u = 1 - h21 * h12;
  60
  61     if (u <= 0.0) {             /* the case u <= 0 is rejected */
  62       P[0] = -1;
  63       P[1] = 0;
  64       P[2] = 0;
  65       P[3] = 0;
  66       P[4] = 0;
  67       *d1 = 0;
  68       *d2 = 0;
  69       *b1 = 0;
  70       return;
  71     }
  72
  73     D1 /= u;
  74     D2 /= u;
  75     x *= u;
  76   } else {
  77     /* case of equation A7 */
  78
  79     if (D2 * y * y < 0.0) {
  80       P[0] = -1;
  81       P[1] = 0;
  82       P[2] = 0;
  83       P[3] = 0;
  84       P[4] = 0;
  85       *d1 = 0;
  86       *d2 = 0;
  87       *b1 = 0;
  88       return;
  89     }
  90
  91     P[0] = 1;
  92
  93     h11 = (D1 * x) / (D2 * y);
  94     h12 = 1;
  95     h21 = -1;
  96     h22 = x / y;
  97
  98     u = 1 + h11 * h22;
  99
 100     D1 /= u;
 101     D2 /= u;
 102
 103     {
 104       BASE tmp = D2;
 105       D2 = D1;
 106       D1 = tmp;
 107     }
 108
 109     x = y * u;
 110   }
 111
 112   /* rescale D1 to range [1/G2,G2] */
 113
 114   while (D1 <= 1.0 / G2 && D1 != 0.0) {
 115     P[0] = -1;
 116     D1 *= G2;
 117     x /= G;
 118     h11 /= G;
 119     h12 /= G;
 120   }
 121
 122   while (D1 >= G2) {
 123     P[0] = -1;
 124     D1 /= G2;
 125     x *= G;
 126     h11 *= G;
 127     h12 *= G;
 128   }
 129
 130   /* rescale D2 to range [1/G2,G2] */
 131
 132   while (fabs(D2) <= 1.0 / G2 && D2 != 0.0) {
 133     P[0] = -1;
 134     D2 *= G2;
 135     h21 /= G;
 136     h22 /= G;
 137   }
 138
 139   while (fabs(D2) >= G2) {
 140     P[0] = -1;
 141     D2 /= G2;
 142     h21 *= G;
 143     h22 *= G;
 144   }
 145
 146   *d1 = D1;
 147   *d2 = D2;
 148   *b1 = x;
 149
 150   if (P[0] == -1.0) {
 151     P[1] = h11;
 152     P[2] = h21;
 153     P[3] = h12;
 154     P[4] = h22;
 155   } else if (P[0] == 0.0) {
 156     P[2] = h21;
 157     P[3] = h12;
 158   } else if (P[0] == 1.0) {
 159     P[1] = h11;
 160     P[4] = h22;
 161   }
 162 }