floraachy
/
OpenBLAS
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Rewrite ROTMG based on the new implementation in GONUM 

based on the algorithm proposed by Tim Hopkins, see issue 1452 for the reference
This commit is contained in:
Martin Kroeker 2018-03-03 14:42:09 +01:00 committed by GitHub
parent 0ab5bf1746
commit 81f81cbc63
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
1 changed files with 58 additions and 68 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

126
interface/rotmg.c
View File

@ -64,6 +64,13 @@ void CNAME(FLOAT *dd1, FLOAT *dd2, FLOAT *dx1, FLOAT dy1, FLOAT *dparam){
FLOAT du, dp1, dp2, dq2, dq1, dh11=ZERO, dh21=ZERO, dh12=ZERO, dh22=ZERO, dflag=-ONE, dtemp; FLOAT du, dp1, dp2, dq2, dq1, dh11=ZERO, dh21=ZERO, dh12=ZERO, dh22=ZERO, dflag=-ONE, dtemp;
if (*dd2 == ZERO || dy1 == ZERO)
{
dflag = -TWO;
dparam[0] = dflag;
return;
}
if(*dd1 < ZERO) if(*dd1 < ZERO)
{ {
dflag = -ONE; dflag = -ONE;
@ -76,6 +83,16 @@ void CNAME(FLOAT *dd1, FLOAT *dd2, FLOAT *dx1, FLOAT dy1, FLOAT *dparam){
*dd2 = ZERO; *dd2 = ZERO;
*dx1 = ZERO; *dx1 = ZERO;
} }
else if ((*dd1 == ZERO || *dx1 == ZERO) && *dd2 > ZERO)
{
dflag = ONE;
dh12 = 1;
dh21 = -1;
*dx1 = dy1;
dtemp = *dd1;
*dd1 = *dd2;
*dd2 = dtemp;
}
else else
{ {
dp2 = *dd2 * dy1; dp2 = *dd2 * dy1;
@ -90,6 +107,9 @@ void CNAME(FLOAT *dd1, FLOAT *dd2, FLOAT *dx1, FLOAT dy1, FLOAT *dparam){
dq1 = dp1 * *dx1; dq1 = dp1 * *dx1;
if(ABS(dq1) > ABS(dq2)) if(ABS(dq1) > ABS(dq2))
{ {
dflag = ZERO;
dh11 = ONE;
dh22 = ONE;
dh21 = - dy1 / *dx1; dh21 = - dy1 / *dx1;
dh12 = dp2 / dp1; dh12 = dp2 / dp1;
@ -100,8 +120,19 @@ void CNAME(FLOAT *dd1, FLOAT *dd2, FLOAT *dx1, FLOAT dy1, FLOAT *dparam){
*dd1 = *dd1 / du; *dd1 = *dd1 / du;
*dd2 = *dd2 / du; *dd2 = *dd2 / du;
*dx1 = *dx1 * du; *dx1 = *dx1 * du;
} else {
dflag = -ONE;
dh11 = ZERO;
dh12 = ZERO;
dh21 = ZERO;
dh22 = ZERO;
*dd1 = ZERO;
*dd2 = ZERO;
*dx1 = ZERO;
} }
} }
else else
{ {
@ -120,7 +151,9 @@ void CNAME(FLOAT *dd1, FLOAT *dd2, FLOAT *dx1, FLOAT dy1, FLOAT *dparam){
} }
else else
{ {
dflag = ONE; dflag = ONE;
dh21 = -ONE;
dh12 = ONE;
dh11 = dp1 / dp2; dh11 = dp1 / dp2;
dh22 = *dx1 / dy1; dh22 = *dx1 / dy1;
@ -134,76 +167,33 @@ void CNAME(FLOAT *dd1, FLOAT *dd2, FLOAT *dx1, FLOAT dy1, FLOAT *dparam){
} }
if(*dd1 != ZERO) while ( *dd1 <= RGAMSQ && *dd1 != ZERO)
{ {
if( (*dd1 <= RGAMSQ) || (*dd1 >= GAMSQ) ) dflag = -ONE;
{ *dd1 = *dd1 * (GAM * GAM);
if(dflag == ZERO) *dx1 = *dx1 / GAM;
{ dh11 = dh11 / GAM;
dh11 = ONE; dh12 = dh12 / GAM;
dh22 = ONE; }
dflag = -ONE; while (ABS(*dd1) > GAMSQ) {
} dflag = -ONE;
else *dd1 = *dd1 / (GAM * GAM);
{ *dx1 = *dx1 * GAM;
dh21 = -ONE; dh11 = dh11 * GAM;
dh12 = ONE; dh12 = dh12 * GAM;
dflag = -ONE;
}
if( *dd1 <= RGAMSQ )
{
while (ABS(*dd1) <= RGAMSQ) {
*dd1 = *dd1 * (GAM * GAM);
*dx1 = *dx1 / GAM;
dh11 = dh11 / GAM;
dh12 = dh12 / GAM;
}
}
else
{
while (ABS(*dd1) >= GAMSQ) {
*dd1 = *dd1 / (GAM * GAM);
*dx1 = *dx1 * GAM;
dh11 = dh11 * GAM;
dh12 = dh12 * GAM;
}
}
}
} }
if(*dd2 != ZERO) while (ABS(*dd2) <= RGAMSQ && *dd2 != ZERO) {
{ dflag = -ONE;
if( (ABS(*dd2) <= RGAMSQ) || (ABS(*dd2) >= GAMSQ) ) *dd2 = *dd2 * (GAM * GAM);
{ dh21 = dh21 / GAM;
if(dflag == ZERO) dh22 = dh22 / GAM;
{ }
dh11 = ONE; while (ABS(*dd2) > GAMSQ) {
dh22 = ONE; dflag = -ONE;
dflag = -ONE; *dd2 = *dd2 / (GAM * GAM);
} dh21 = dh21 * GAM;
else dh22 = dh22 * GAM;
{
dh21 = -ONE;
dh12 = ONE;
dflag = -ONE;
}
if( ABS(*dd2) <= RGAMSQ )
{
while (ABS(*dd2) <= RGAMSQ) {
*dd2 = *dd2 * (GAM * GAM);
dh21 = dh21 / GAM;
dh22 = dh22 / GAM;
}
}
else
{
while (ABS(*dd2) >= GAMSQ) {
*dd2 = *dd2 / (GAM * GAM);
dh21 = dh21 * GAM;
dh22 = dh22 * GAM;
}
}
}
} }
} }