Home Browse Files
      1
      2
      3
      4
      5
      6
      7
      8
      9
     10
     11
     12
     13
     14
     15
     16
     17
     18
     19
     20
     21
     22
     23
     24
     25
     26
     27
     28
     29
     30
     31
     32
     33
     34
     35
     36
     37
     38
     39
     40
     41
     42
     43
     44
     45
     46
     47
     48
     49
     50
     51
     52
     53
     54
     55
#include <ulmblas.h>
#include <math.h>

static double
sign(const double x, const double y)
{
    if (y>=0.0) {
        return fabs(x);
    } else {
        return -fabs(x);
    }
}

void
ULMBLAS(drotg)(double *a,
               double *b,
               double *c,
               double *s)
{
    double r, roe = *b, 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;
}

void
F77BLAS(drotg)(double *a,
               double *b,
               double *c,
               double *s)
{
    ULMBLAS(drotg)(a, b, c, s);
}





































































































Home