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