110 lines
1.6 KiB
C
110 lines
1.6 KiB
C
#include <math.h>
|
|
#include "common.h"
|
|
#ifdef FUNCTION_PROFILE
|
|
#include "functable.h"
|
|
#endif
|
|
|
|
#ifndef CBLAS
|
|
|
|
void NAME(FLOAT *DA, FLOAT *DB, FLOAT *C, FLOAT *S){
|
|
|
|
#else
|
|
|
|
void CNAME(FLOAT *DA, FLOAT *DB, FLOAT *C, FLOAT *S){
|
|
|
|
#endif
|
|
|
|
|
|
#if defined(__i386__) || defined(__x86_64__) || defined(__ia64__)
|
|
|
|
long double da = *DA;
|
|
long double db = *DB;
|
|
long double c;
|
|
long double s;
|
|
long double r, roe, z;
|
|
|
|
long double ada = fabs(da);
|
|
long double adb = fabs(db);
|
|
long double scale = ada + adb;
|
|
|
|
#ifndef CBLAS
|
|
PRINT_DEBUG_NAME;
|
|
#else
|
|
PRINT_DEBUG_CNAME;
|
|
#endif
|
|
|
|
roe = db;
|
|
if (ada > adb) roe = da;
|
|
|
|
if (scale == ZERO) {
|
|
*C = ONE;
|
|
*S = ZERO;
|
|
*DA = ZERO;
|
|
*DB = ZERO;
|
|
} else {
|
|
r = sqrt(da * da + db * db);
|
|
if (roe < 0) r = -r;
|
|
c = da / r;
|
|
s = db / r;
|
|
z = ONE;
|
|
if (da != ZERO) {
|
|
if (ada > adb){
|
|
z = s;
|
|
} else {
|
|
z = ONE / c;
|
|
}
|
|
}
|
|
|
|
*C = c;
|
|
*S = s;
|
|
*DA = r;
|
|
*DB = z;
|
|
}
|
|
|
|
#else
|
|
FLOAT da = *DA;
|
|
FLOAT db = *DB;
|
|
FLOAT c = *C;
|
|
FLOAT s = *S;
|
|
FLOAT r, roe, z;
|
|
|
|
FLOAT ada = fabs(da);
|
|
FLOAT adb = fabs(db);
|
|
FLOAT scale = ada + adb;
|
|
|
|
#ifndef CBLAS
|
|
PRINT_DEBUG_NAME;
|
|
#else
|
|
PRINT_DEBUG_CNAME;
|
|
#endif
|
|
|
|
roe = db;
|
|
if (ada > adb) roe = da;
|
|
|
|
if (scale == ZERO) {
|
|
*C = ONE;
|
|
*S = ZERO;
|
|
*DA = ZERO;
|
|
*DB = ZERO;
|
|
} else {
|
|
FLOAT aa = da / scale;
|
|
FLOAT bb = db / scale;
|
|
|
|
r = scale * sqrt(aa * aa + bb * bb);
|
|
if (roe < 0) r = -r;
|
|
c = da / r;
|
|
s = db / r;
|
|
z = ONE;
|
|
if (ada > adb) z = s;
|
|
if ((ada < adb) && (c != ZERO)) z = ONE / c;
|
|
|
|
*C = c;
|
|
*S = s;
|
|
*DA = r;
|
|
*DB = z;
|
|
}
|
|
#endif
|
|
|
|
return;
|
|
}
|