1 files changed, 162 insertions, 0 deletions
diff --git a/gsl-1.9/cblas/source_rotmg.h b/gsl-1.9/cblas/source_rotmg.h
new file mode 100644
index 0000000..7a41538
--- /dev/null
+++ b/gsl-1.9/cblas/source_rotmg.h
@@ -0,0 +1,162 @@
+/* blas/source_rotmg.h
+ * 
+ * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough
+ * 
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 2 of the License, or (at
+ * your option) any later version.
+ * 
+ * This program is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * General Public License for more details.
+ * 
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ */
+
+{
+  const BASE G = 4096.0, G2 = G * G;
+  BASE D1 = *d1, D2 = *d2, x = *b1, y = b2;
+  BASE h11, h12, h21, h22, u;
+
+  BASE c, s;
+
+  /* case of d1 < 0, appendix A, second to last paragraph */
+
+  if (D1 < 0.0) {
+    P[0] = -1;
+    P[1] = 0;
+    P[2] = 0;
+    P[3] = 0;
+    P[4] = 0;
+    *d1 = 0;
+    *d2 = 0;
+    *b1 = 0;
+    return;
+  }
+
+  if (D2 * y == 0.0) {
+    P[0] = -2;                  /* case of H = I, page 315 */
+    return;
+  }
+
+  c = fabs(D1 * x * x);
+  s = fabs(D2 * y * y);
+
+  if (c > s) {
+    /* case of equation A6 */
+
+    P[0] = 0.0;
+
+    h11 = 1;
+    h12 = (D2 * y) / (D1 * x);
+    h21 = -y / x;
+    h22 = 1;
+
+    u = 1 - h21 * h12;
+
+    if (u <= 0.0) {             /* the case u <= 0 is rejected */
+      P[0] = -1;
+      P[1] = 0;
+      P[2] = 0;
+      P[3] = 0;
+      P[4] = 0;
+      *d1 = 0;
+      *d2 = 0;
+      *b1 = 0;
+      return;
+    }
+
+    D1 /= u;
+    D2 /= u;
+    x *= u;
+  } else {
+    /* case of equation A7 */
+
+    if (D2 * y * y < 0.0) {
+      P[0] = -1;
+      P[1] = 0;
+      P[2] = 0;
+      P[3] = 0;
+      P[4] = 0;
+      *d1 = 0;
+      *d2 = 0;
+      *b1 = 0;
+      return;
+    }
+
+    P[0] = 1;
+
+    h11 = (D1 * x) / (D2 * y);
+    h12 = 1;
+    h21 = -1;
+    h22 = x / y;
+
+    u = 1 + h11 * h22;
+
+    D1 /= u;
+    D2 /= u;
+
+    {
+      BASE tmp = D2;
+      D2 = D1;
+      D1 = tmp;
+    }
+
+    x = y * u;
+  }
+
+  /* rescale D1 to range [1/G2,G2] */
+
+  while (D1 <= 1.0 / G2 && D1 != 0.0) {
+    P[0] = -1;
+    D1 *= G2;
+    x /= G;
+    h11 /= G;
+    h12 /= G;
+  }
+
+  while (D1 >= G2) {
+    P[0] = -1;
+    D1 /= G2;
+    x *= G;
+    h11 *= G;
+    h12 *= G;
+  }
+
+  /* rescale D2 to range [1/G2,G2] */
+
+  while (fabs(D2) <= 1.0 / G2 && D2 != 0.0) {
+    P[0] = -1;
+    D2 *= G2;
+    h21 /= G;
+    h22 /= G;
+  }
+
+  while (fabs(D2) >= G2) {
+    P[0] = -1;
+    D2 /= G2;
+    h21 *= G;
+    h22 *= G;
+  }
+
+  *d1 = D1;
+  *d2 = D2;
+  *b1 = x;
+
+  if (P[0] == -1.0) {
+    P[1] = h11;
+    P[2] = h21;
+    P[3] = h12;
+    P[4] = h22;
+  } else if (P[0] == 0.0) {
+    P[2] = h21;
+    P[3] = h12;
+  } else if (P[0] == 1.0) {
+    P[1] = h11;
+    P[4] = h22;
+  }
+}