interface/rotg.c

   1 #include <math.h>
   2 #include "common.h"
   3 #ifdef FUNCTION_PROFILE
   4 #include "functable.h"
   5 #endif
   6
   7 #ifndef CBLAS
   8
   9 void NAME(FLOAT *DA, FLOAT *DB, FLOAT *C, FLOAT *S){
  10
  11 #else
  12
  13 void CNAME(FLOAT *DA, FLOAT *DB, FLOAT *C, FLOAT *S){
  14
  15 #endif
  16
  17 #if defined(__i386__) || defined(__x86_64__) || defined(__ia64__) || defined(_M_X64) || defined(_M_IX86)
  18
  19   long double da = *DA;
  20   long double db = *DB;
  21   long double c;
  22   long double s;
  23   long double r, roe, z;
  24
  25   long double ada = fabsl(da);
  26   long double adb = fabsl(db);
  27   long double scale = ada + adb;
  28
  29 #ifndef CBLAS
  30   PRINT_DEBUG_NAME;
  31 #else
  32   PRINT_DEBUG_CNAME;
  33 #endif
  34
  35   roe = db;
  36   if (ada > adb) roe = da;
  37
  38   if (scale == ZERO) {
  39     *C = ONE;
  40     *S = ZERO;
  41     *DA = ZERO;
  42     *DB = ZERO;
  43   } else {
  44     r = sqrt(da * da + db * db);
  45     if (roe < 0) r = -r;
  46     c = da / r;
  47     s = db / r;
  48     z = ONE;
  49     if (da != ZERO) {
  50       if (ada > adb){
  51         z = s;
  52       } else {
  53         z = ONE / c;
  54       }
  55     }
  56
  57     *C = c;
  58     *S = s;
  59     *DA = r;
  60     *DB = z;
  61   }
  62
  63 #else
  64   FLOAT da = *DA;
  65   FLOAT db = *DB;
  66   FLOAT c  = *C;
  67   FLOAT s  = *S;
  68   FLOAT r, roe, z;
  69
  70   FLOAT ada = fabs(da);
  71   FLOAT adb = fabs(db);
  72   FLOAT scale = ada + adb;
  73
  74 #ifndef CBLAS
  75   PRINT_DEBUG_NAME;
  76 #else
  77   PRINT_DEBUG_CNAME;
  78 #endif
  79
  80   roe = db;
  81   if (ada > adb) roe = da;
  82
  83   if (scale == ZERO) {
  84     *C = ONE;
  85     *S = ZERO;
  86     *DA = ZERO;
  87     *DB = ZERO;
  88   } else {
  89     FLOAT aa = da / scale;
  90     FLOAT bb = db / scale;
  91
  92     r = scale * sqrt(aa * aa + bb * bb);
  93     if (roe < 0) r = -r;
  94     c = da / r;
  95     s = db / r;
  96     z = ONE;
  97     if (ada > adb) z = s;
  98     if ((ada <= adb) && (c != ZERO)) z = ONE / c;
  99
 100     *C = c;
 101     *S = s;
 102     *DA = r;
 103     *DB = z;
 104   }
 105 #endif
 106
 107   return;
 108 }