Merge pull request #5116 from martin-frbg/issue5110

Handle INCX=0 in ?NRM2
OpenMathLib · Feb 9, 2025 · f42ce70 · f42ce70
2 parents c54f541 + 7478c10
commit f42ce70
Showing 1 changed file with 63 additions and 1 deletion.
diff --git 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,44 @@ FLOAT CNAME(blasint n, FLOAT *x, blasint incx){
 
   if (n <= 0) return 0.;
 
-  #ifndef COMPLEX
+#ifndef COMPLEX
   if (n == 1)
 #ifdef DOUBLE
     return fabs(x[0]);
 #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)