Completed the port of refblas

level1/dcopy.c

@@ -16,6 +16,7 @@ dcopy_(const int    *_n,
 //  Local scalars
 //
     int    i, m;
+
 //
 //  Quick return if possible
 //

level1/drot.c

@@ -1,10 +1,56 @@
 void
-drot_(const int      *n,
-      const double   *x,
-      const int      *incX,
-      const double   *y,
-      const int      *incY,
-      const double   *c,
-      const double   *s)
+drot_(const int      *_n,
+      double         *x,
+      const int      *_incX,
+      double         *y,
+      const int      *_incY,
+      const double   *_c,
+      const double   *_s)
 {
+//
+//  Dereference scalar parameters
+//
+    int n    = *_n;
+    int incX = *_incX;
+    int incY = *_incY;
+    double c = *_c;
+    double s = *_s;
+
+//
+//  Local scalars
+//
+    int    i;
+    double tmp;
+
+//
+//  Quick return if possible
+//
+    if (n==0) {
+        return;
+    }
+    if (incX==1 && incY==1) {
+//
+//      Code for both increments equal to 1
+//
+        for (i=0; i<n; ++i) {
+            tmp  = c*x[i] + s*y[i];
+            y[i] = c*y[i] - s*x[i];
+            x[i] = tmp;
+        }
+    } else {
+//
+//      Code for unequal increments or equal increments not equal to 1
+//
+        if (incX<0) {
+            x -= incX*(n-1);
+        }
+        if (incY<0) {
+            y -= incY*(n-1);
+        }
+        for (i=0; i<n; ++i, x+=incX, y+=incY) {
+            tmp  = c*(*x) + s*(*y);
+            (*y) = c*(*y) - s*(*x);
+            (*x) = tmp;
+        }
+    }
 }

level1/drotg.c

@@ -1,7 +1,45 @@
+#include <math.h>
+
+static double
+sign(const double x, const double y)
+{
+    if (y>=0.0) {
+        return fabs(x);
+    } else {
+        return -fabs(x);
+    }
+}
+
 void
 drotg_(double *a,
        double *b,
        double *c,
        double *s)
 {
+    double r, roe, scale, z;
+
+    if (fabs(*a)>fabs(*b)) {
+        roe = *a;
+    }
+    scale = fabs(*a) + fabs(*b);
+    if (scale==0.0) {
+        *c = 1.0;
+        *s = 0.0;
+        r = 0.0;
+        z = 0.0;
+    } else {
+        r = scale*sqrt(pow(*a/scale, 2.0) + pow(*b/scale, 2.0));
+        r = sign(1.0, roe)*r;
+        *c = *a / r;
+        *s = *b / r;
+        z = 1.0;
+        if (fabs(*a)>fabs(*b)) {
+            z = *s;
+        }
+        if (fabs(*b)>=fabs(*a) && *c!=0.0) {
+            z = 1.0/(*c);
+        }
+    }
+    *a = r;
+    *b = z;
 }