From 0fd5448b2c3a3e16410bff49489f968c44e51d91 Mon Sep 17 00:00:00 2001 From: Martin Kroeker Date: Sat, 8 Feb 2025 19:33:05 +0100 Subject: [PATCH 1/2] Handle INCX=0 --- interface/nrm2.c | 70 +++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 69 insertions(+), 1 deletion(-) diff --git a/interface/nrm2.c b/interface/nrm2.c index 331ebc3d05..22b24396f0 100644 --- a/interface/nrm2.c +++ b/interface/nrm2.c @@ -61,6 +61,37 @@ FLOATRET NAME(blasint *N, FLOAT *x, blasint *INCX){ #else return fabsf(x[0]); #endif +#endif + + if (incx == 0) +#ifndef COMPLEX +#ifdef DOUBLE + return (sqrt((double)n)*fabs(x[0])); +#else + return (sqrt((float)n)*fabsf(x[0])); +#endif +#else +#ifdef DOUBLE + { + double fr=fabs(x[0]); + double fi=fabs(x[1]); + double fmin=MIN(fr,fi); + double fmax=MAX(fr,fi); + if (fmax==0.) return(fmax); + if (fmax==fmin) return(sqrt((double)n)*sqrt(2.)*fmax); + return (sqrt((double)n) * fmax * sqrt (1. + (fmin/fmax)*(fmin/fmax))); + } +#else + { + float fr=fabs(x[0]); + float fi=fabs(x[1]); + float fmin=MIN(fr,fi); + float fmax=MAX(fr,fi); + if (fmax==0.) return(fmax); + if (fmax==fmin) return(sqrt((float)n)*sqrt(2.)*fmax); + return (sqrt((float)n) * fmax * sqrt (1. + (fmin/fmax)*(fmin/fmax))); + } +#endif #endif if (incx < 0) @@ -97,13 +128,50 @@ FLOAT CNAME(blasint n, FLOAT *x, blasint incx){ if (n <= 0) return 0.; - #ifndef COMPLEX if (n == 1) +#ifndef COMPLEX #ifdef DOUBLE return fabs(x[0]); #else return fabsf(x[0]); #endif +#else +#ifdef DOUBLE + return fabs(x[0]+fabs(x[1])); +#else + return fabsf(x[0]+fabsf(x[1])); +#endif +#endif + + if (incx == 0) +#ifndef COMPLEX +#ifdef DOUBLE + return (sqrt((double)n)*fabs(x[0])); +#else + return (sqrt((float)n)*fabsf(x[0])); +#endif +#else +#ifdef DOUBLE + { + double fr=fabs(x[0]); + double fi=fabs(x[1]); + double fmin=MIN(fr,fi); + double fmax=MAX(fr,fi); + if (fmax==0.) return(fmax); + if (fmax==fmin) return(sqrt((double)n)*sqrt(2.)*fmax); + return (sqrt((double)n) * fmax * sqrt (1. + (fmin/fmax)*(fmin/fmax))); + } +#else + { + float fr=fabs(x[0]); + float fi=fabs(x[1]); + float fmin=MIN(fr,fi); + float fmax=MAX(fr,fi); + if (fmax==0.) return(fmax); + if (fmax==fmin) return(sqrt((float)n)*sqrt(2.)*fmax); + return (sqrt((float)n) * fmax * sqrt (1. + (fmin/fmax)*(fmin/fmax))); + } +#endif #endif if (incx < 0) From 60d0be0e971c064ea7a2bee0cb1247c387e18b3c Mon Sep 17 00:00:00 2001 From: Martin Kroeker Date: Sat, 8 Feb 2025 23:42:21 +0100 Subject: [PATCH 2/2] Update nrm2.c --- interface/nrm2.c | 8 +------- 1 file changed, 1 insertion(+), 7 deletions(-) diff --git a/interface/nrm2.c b/interface/nrm2.c index 22b24396f0..cfeb13df8f 100644 --- a/interface/nrm2.c +++ b/interface/nrm2.c @@ -128,19 +128,13 @@ FLOAT CNAME(blasint n, FLOAT *x, blasint incx){ if (n <= 0) return 0.; - if (n == 1) #ifndef COMPLEX + if (n == 1) #ifdef DOUBLE return fabs(x[0]); #else return fabsf(x[0]); #endif -#else -#ifdef DOUBLE - return fabs(x[0]+fabs(x[1])); -#else - return fabsf(x[0]+fabsf(x[1])); -#endif #endif if (incx == 0)