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