From d70f2b089dac8b9e4c19295dfa6014c36afee2ec Mon Sep 17 00:00:00 2001 From: "Field G. Van Zee" Date: Sat, 2 Nov 2013 17:19:40 -0500 Subject: [PATCH] Added scaling to abval2s, sqrt2s macros. Details: - Re-defined abval2s and sqrt2s macros to use scaling to avoid underflow and overflow from squaring the real and imaginary components. (This is the same technique used to fix recent bugs in invscals/invscaljs and inverts.) --- frame/include/level0/bli_abval2s.h | 56 +++++++++++++++++++-------- frame/include/level0/bli_sqrt2s.h | 62 +++++++++++++++++++----------- 2 files changed, 79 insertions(+), 39 deletions(-) diff --git a/frame/include/level0/bli_abval2s.h b/frame/include/level0/bli_abval2s.h index 92c7eb253..bc3a44681 100644 --- a/frame/include/level0/bli_abval2s.h +++ b/frame/include/level0/bli_abval2s.h @@ -84,46 +84,70 @@ #define bli_csabval2s( x, a ) \ { \ - bli_csabsq2s( x, a ); \ - bli_sssqrt2s( a, a ); \ + float s = bli_fmaxabs( bli_creal(x), bli_cimag(x) ); \ + float mag = sqrtf( s ) * \ + sqrtf( ( bli_creal(x) / s ) * bli_creal(x) + \ + ( bli_cimag(x) / s ) * bli_cimag(x) ); \ + bli_sssetris( mag, 0.0, a ); \ } #define bli_zsabval2s( x, a ) \ { \ - bli_zsabsq2s( x, a ); \ - bli_sssqrt2s( a, a ); \ + double s = bli_fmaxabs( bli_zreal(x), bli_zimag(x) ); \ + double mag = sqrt( s ) * \ + sqrt( ( bli_zreal(x) / s ) * bli_zreal(x) + \ + ( bli_zimag(x) / s ) * bli_zimag(x) ); \ + bli_dssetris( mag, 0.0, a ); \ } #define bli_cdabval2s( x, a ) \ { \ - bli_cdabsq2s( x, a ); \ - bli_ddsqrt2s( a, a ); \ + double s = bli_fmaxabs( bli_creal(x), bli_cimag(x) ); \ + double mag = sqrt( s ) * \ + sqrt( ( bli_creal(x) / s ) * bli_creal(x) + \ + ( bli_cimag(x) / s ) * bli_cimag(x) ); \ + bli_ddsetris( mag, 0.0, a ); \ } #define bli_zdabval2s( x, a ) \ { \ - bli_zdabsq2s( x, a ); \ - bli_ddsqrt2s( a, a ); \ + double s = bli_fmaxabs( bli_zreal(x), bli_zimag(x) ); \ + double mag = sqrt( s ) * \ + sqrt( ( bli_zreal(x) / s ) * bli_zreal(x) + \ + ( bli_zimag(x) / s ) * bli_zimag(x) ); \ + bli_ddsetris( mag, 0.0, a ); \ } #define bli_ccabval2s( x, a ) \ { \ - bli_ccabsq2s( x, a ); \ - bli_ccsqrt2s( a, a ); \ + float s = bli_fmaxabs( bli_creal(x), bli_cimag(x) ); \ + float mag = sqrtf( s ) * \ + sqrtf( ( bli_creal(x) / s ) * bli_creal(x) + \ + ( bli_cimag(x) / s ) * bli_cimag(x) ); \ + bli_scsetris( mag, 0.0, a ); \ } #define bli_zcabval2s( x, a ) \ { \ - bli_zcabsq2s( x, a ); \ - bli_ccsqrt2s( a, a ); \ + double s = bli_fmaxabs( bli_zreal(x), bli_zimag(x) ); \ + double mag = sqrt( s ) * \ + sqrt( ( bli_zreal(x) / s ) * bli_zreal(x) + \ + ( bli_zimag(x) / s ) * bli_zimag(x) ); \ + bli_dcsetris( mag, 0.0, a ); \ } #define bli_czabval2s( x, a ) \ { \ - bli_czabsq2s( x, a ); \ - bli_zzsqrt2s( a, a ); \ + double s = bli_fmaxabs( bli_creal(x), bli_cimag(x) ); \ + double mag = sqrt( s ) * \ + sqrt( ( bli_creal(x) / s ) * bli_creal(x) + \ + ( bli_cimag(x) / s ) * bli_cimag(x) ); \ + bli_dzsetris( mag, 0.0, a ); \ } #define bli_zzabval2s( x, a ) \ { \ - bli_zzabsq2s( x, a ); \ - bli_zzsqrt2s( a, a ); \ + double s = bli_fmaxabs( bli_zreal(x), bli_zimag(x) ); \ + double mag = sqrt( s ) * \ + sqrt( ( bli_zreal(x) / s ) * bli_zreal(x) + \ + ( bli_zimag(x) / s ) * bli_zimag(x) ); \ + bli_dzsetris( mag, 0.0, a ); \ } diff --git a/frame/include/level0/bli_sqrt2s.h b/frame/include/level0/bli_sqrt2s.h index 396e4afc5..0ad60c1a5 100644 --- a/frame/include/level0/bli_sqrt2s.h +++ b/frame/include/level0/bli_sqrt2s.h @@ -55,15 +55,19 @@ } #define bli_cssqrt2s( x, a ) \ { \ - float mag = sqrtf( bli_creal(x) * bli_creal(x) + \ - bli_cimag(x) * bli_cimag(x) ); \ + float s = bli_fmaxabs( bli_creal(x), bli_cimag(x) ); \ + float mag = sqrtf( s ) * \ + sqrtf( ( bli_creal(x) / s ) * bli_creal(x) + \ + ( bli_cimag(x) / s ) * bli_cimag(x) ); \ \ - (a) = sqrtf( ( mag + bli_creal(x) ) / 2.0F ); \ + (a) = sqrtf( ( mag + bli_creal(x) ) / 2.0 ); \ } #define bli_zssqrt2s( x, a ) \ { \ - double mag = sqrt( bli_zreal(x) * bli_zreal(x) + \ - bli_zimag(x) * bli_zimag(x) ); \ + double s = bli_fmaxabs( bli_zreal(x), bli_zimag(x) ); \ + double mag = sqrt( s ) * \ + sqrt( ( bli_zreal(x) / s ) * bli_zreal(x) + \ + ( bli_zimag(x) / s ) * bli_zimag(x) ); \ \ (a) = sqrt( ( mag + bli_zreal(x) ) / 2.0 ); \ } @@ -79,15 +83,19 @@ } #define bli_cdsqrt2s( x, a ) \ { \ - float mag = sqrtf( bli_creal(x) * bli_creal(x) + \ - bli_cimag(x) * bli_cimag(x) ); \ + double s = bli_fmaxabs( bli_creal(x), bli_cimag(x) ); \ + double mag = sqrt( s ) * \ + sqrt( ( bli_creal(x) / s ) * bli_creal(x) + \ + ( bli_cimag(x) / s ) * bli_cimag(x) ); \ \ - (a) = sqrtf( ( mag + bli_creal(x) ) / 2.0F ); \ + (a) = sqrt( ( mag + bli_creal(x) ) / 2.0 ); \ } #define bli_zdsqrt2s( x, a ) \ { \ - double mag = sqrt( bli_zreal(x) * bli_zreal(x) + \ - bli_zimag(x) * bli_zimag(x) ); \ + double s = bli_fmaxabs( bli_zreal(x), bli_zimag(x) ); \ + double mag = sqrt( s ) * \ + sqrt( ( bli_zreal(x) / s ) * bli_zreal(x) + \ + ( bli_zimag(x) / s ) * bli_zimag(x) ); \ \ (a) = sqrt( ( mag + bli_zreal(x) ) / 2.0 ); \ } @@ -105,16 +113,20 @@ } #define bli_ccsqrt2s( x, a ) \ { \ - float mag = sqrtf( bli_creal(x) * bli_creal(x) + \ - bli_cimag(x) * bli_cimag(x) ); \ + float s = bli_fmaxabs( bli_creal(x), bli_cimag(x) ); \ + float mag = sqrtf( s ) * \ + sqrtf( ( bli_creal(x) / s ) * bli_creal(x) + \ + ( bli_cimag(x) / s ) * bli_cimag(x) ); \ \ - bli_creal(a) = sqrtf( ( mag + bli_creal(x) ) / 2.0F ); \ - bli_cimag(a) = sqrtf( ( mag - bli_cimag(x) ) / 2.0F ); \ + bli_creal(a) = sqrtf( ( mag + bli_creal(x) ) / 2.0 ); \ + bli_cimag(a) = sqrtf( ( mag - bli_cimag(x) ) / 2.0 ); \ } #define bli_zcsqrt2s( x, a ) \ { \ - double mag = sqrt( bli_zreal(x) * bli_zreal(x) + \ - bli_zimag(x) * bli_zimag(x) ); \ + double s = bli_fmaxabs( bli_zreal(x), bli_zimag(x) ); \ + double mag = sqrt( s ) * \ + sqrt( ( bli_zreal(x) / s ) * bli_zreal(x) + \ + ( bli_zimag(x) / s ) * bli_zimag(x) ); \ \ bli_creal(a) = sqrt( ( mag + bli_zreal(x) ) / 2.0 ); \ bli_cimag(a) = sqrt( ( mag - bli_zimag(x) ) / 2.0 ); \ @@ -129,20 +141,24 @@ #define bli_dzsqrt2s( x, a ) \ { \ bli_zreal(a) = sqrt( (x) ); \ - bli_zimag(a) = 0.0F; \ + bli_zimag(a) = 0.0; \ } #define bli_czsqrt2s( x, a ) \ { \ - float mag = sqrtf( bli_creal(x) * bli_creal(x) + \ - bli_cimag(x) * bli_cimag(x) ); \ + double s = bli_fmaxabs( bli_creal(x), bli_cimag(x) ); \ + double mag = sqrt( s ) * \ + sqrt( ( bli_creal(x) / s ) * bli_creal(x) + \ + ( bli_cimag(x) / s ) * bli_cimag(x) ); \ \ - bli_zreal(a) = sqrtf( ( mag + bli_creal(x) ) / 2.0F ); \ - bli_zimag(a) = sqrtf( ( mag - bli_cimag(x) ) / 2.0F ); \ + bli_zreal(a) = sqrt( ( mag + bli_creal(x) ) / 2.0 ); \ + bli_zimag(a) = sqrt( ( mag - bli_cimag(x) ) / 2.0 ); \ } #define bli_zzsqrt2s( x, a ) \ { \ - double mag = sqrt( bli_zreal(x) * bli_zreal(x) + \ - bli_zimag(x) * bli_zimag(x) ); \ + double s = bli_fmaxabs( bli_zreal(x), bli_zimag(x) ); \ + double mag = sqrt( s ) * \ + sqrt( ( bli_zreal(x) / s ) * bli_zreal(x) + \ + ( bli_zimag(x) / s ) * bli_zimag(x) ); \ \ bli_zreal(a) = sqrt( ( mag + bli_zreal(x) ) / 2.0 ); \ bli_zimag(a) = sqrt( ( mag - bli_zimag(x) ) / 2.0 ); \