LoongArch64: Fixed snrm2_lsx.S

Author: XiWeiGu

kernel/arm/nrm2.c

@@ -52,18 +52,32 @@ USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 FLOAT CNAME(BLASLONG n, FLOAT *x, BLASLONG inc_x)
 {
 	BLASLONG i=0;
+    BLASLONG j=0;
 	FLOAT scale = 0.0;
-	FLOAT ssq   = 1.0;
+	FLOAT ssq   = 0.0;
 	FLOAT absxi = 0.0;
+    FLOAT max = 0.0;
 
 
 	if (n <= 0 || inc_x == 0) return(0.0);
 	if ( n == 1 ) return( ABS(x[0]) );
 
 	n *= inc_x;
+    // Find Max value
+	while(abs(j) < abs(n))
+	{
+        if (max < ABS(x[j])) max = ABS(x[j]);
+
+		j += inc_x;
+	}
+
+    if (max == 0.0) return (0.0);
+    
+
 	while(abs(i) < abs(n))
 	{
 
+#if 0
 		if ( x[i] != 0.0 )
 		{
 			absxi = ABS( x[i] );
@@ -78,9 +92,11 @@ FLOAT CNAME(BLASLONG n, FLOAT *x, BLASLONG inc_x)
 			}
 
 		}
+#endif
+        ssq += (x[i] / max) * (x[i] / max);
 		i += inc_x;
 	}
-	scale = scale * sqrt( ssq );
+	scale =  max * sqrt( ssq );
 	return(scale);
 
 }

kernel/loongarch64/KERNEL.LA264

@@ -130,12 +130,12 @@ ZROTKERNEL  =  crot_lsx.S
 #CROTKERNEL   = ../arm/zrot.c
 #ZROTKERNEL   = ../arm/zrot.c
 
-#SNRM2KERNEL =  snrm2_lsx.S
-#DNRM2KERNEL =  dnrm2_lsx.S
+SNRM2KERNEL =  snrm2_lsx.S
+DNRM2KERNEL =  dnrm2_lsx.S
 #CNRM2KERNEL =  cnrm2_lsx.S
 #ZNRM2KERNEL =  znrm2_lsx.S
-SNRM2KERNEL  = ../arm/nrm2.c
-DNRM2KERNEL  = ../arm/nrm2.c
+#SNRM2KERNEL  = ../arm/nrm2.c
+#DNRM2KERNEL  = ../arm/nrm2.c
 CNRM2KERNEL  = ../arm/znrm2.c
 ZNRM2KERNEL  = ../arm/znrm2.c
 

kernel/loongarch64/snrm2_lsx.S

@@ -52,17 +52,52 @@ USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 /* Don't change following FR unless you know the effects. */
 #define res1   $vr19
 #define res2   $vr20
+#define RCP    $f2
+#define VALPHA $vr3
+
+// The optimization for snrm2 cannot simply involve
+// extending the data type from float to double and
+// then summing the squares of the data. LAPACK tests
+// have shown that this approach can still lead to data overflow.
+// Instead, we need to find the maximum absolute value in the entire
+// array and divide each data element by this maximum value before
+// performing the calculation. This approach can avoid overflow (and does not require extending the data type).
 
    PROLOGUE
 
 #ifdef F_INTERFACE
    LDINT   N,     0(N)
    LDINT   INCX,  0(INCX)
 #endif
-   vxor.v res1, res1, res1
-   vxor.v res2, res2, res2
    bge $r0,    N, .L999
    beq $r0, INCX, .L999
+
+   addi.d  $sp,     $sp,     -64
+   st.d    $ra,     $sp,     0
+   st.d    $s0,     $sp,     8
+   st.d    $s1,     $sp,     16
+   st.d    $s2,     $sp,     24
+   st.d    $s3,     $sp,     32
+   st.d    N,       $sp,     40
+   st.d    X,       $sp,     48
+   st.d    INCX,    $sp,     56
+   bl samax_k
+   ld.d    $ra,     $sp,     0
+   ld.d    $s0,     $sp,     8
+   ld.d    $s1,     $sp,     16
+   ld.d    $s2,     $sp,     24
+   ld.d    $s3,     $sp,     32
+   ld.d    N,       $sp,     40
+   ld.d    X,       $sp,     48
+   ld.d    INCX,    $sp,     56
+   addi.d  $sp,     $sp,     64
+
+   frecip.s RCP, $f0
+   vreplvei.w   VALPHA, $vr2, 0
+   vxor.v res1, res1, res1
+   vxor.v res2, res2, res2
+   fcmp.ceq.s  $fcc0, $f0,  $f19
+   bcnez  $fcc0, .L999
    li.d  TEMP, SIZE
    slli.d INCX, INCX, BASE_SHIFT
    srai.d I, N, 3
@@ -75,14 +110,12 @@ USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
    vld  VX5, X, 4 * SIZE
    addi.d I, I, -1
    addi.d X, X, 8 * SIZE
-   vfcvtl.d.s VX1, VX0
-   vfcvth.d.s VX2, VX0
-   vfcvtl.d.s VX3, VX5
-   vfcvth.d.s VX4, VX5
-   vfmadd.d  res1, VX1, VX1, res1
-   vfmadd.d  res2, VX2, VX2, res2
-   vfmadd.d  res1, VX3, VX3, res1
-   vfmadd.d  res2, VX4, VX4, res2
+
+   vfmul.s  VX0,  VX0,  VALPHA
+   vfmul.s  VX5,  VX5,  VALPHA
+
+   vfmadd.s res1, VX0, VX0, res1
+   vfmadd.s res2, VX5, VX5, res2
    blt $r0, I, .L10
    b .L996
    .align 3
@@ -104,10 +137,9 @@ USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
    vinsgr2vr.w VX0, t2, 1
    vinsgr2vr.w VX0, t3, 2
    vinsgr2vr.w VX0, t4, 3
-   vfcvtl.d.s VX1, VX0
-   vfcvth.d.s VX2, VX0
-   vfmadd.d res1, VX1, VX1, res1
-   vfmadd.d res2, VX2, VX2, res2
+   vfmul.s  VX0,  VX0,  VALPHA
+   vfmadd.s res1, VX0, VX0, res1
+
    ld.w t1, X, 0
    add.d X, X, INCX
    ld.w t2, X, 0
@@ -120,19 +152,20 @@ USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
    vinsgr2vr.w VX0, t2, 1
    vinsgr2vr.w VX0, t3, 2
    vinsgr2vr.w VX0, t4, 3
-   vfcvtl.d.s VX3, VX0
-   vfcvth.d.s VX4, VX0
-   vfmadd.d res1, VX3, VX3, res1
-   vfmadd.d res2, VX4, VX4, res2
+   vfmul.s  VX0,  VX0,  VALPHA
+   vfmadd.s res2, VX0, VX0, res2
    addi.d  I, I, -1
    blt $r0, I, .L21
-   b .L996
    .align 3
 
 .L996:
-   vfadd.d res1, res1, res2
-   vreplvei.d VX1, res1, 1
-   vfadd.d res1, VX1, res1
+   vfadd.s res1, res1, res2
+   vreplvei.w VX1, res1, 1
+   vreplvei.w VX2, res1, 2
+   vreplvei.w VX3, res1, 3
+   vfadd.s res1, VX1, res1
+   vfadd.s res1, VX2, res1
+   vfadd.s res1, VX3, res1
    .align 3
 
 .L997:
@@ -143,16 +176,16 @@ USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 .L998:
    fld.s $f15, X, 0
    addi.d I, I, -1
-   fcvt.d.s $f15, $f15
-   fmadd.d $f19, $f15, $f15, $f19
+   fmul.s  $f15, $f15, RCP
+   fmadd.s $f19, $f15, $f15, $f19
    add.d X, X, INCX
    blt $r0, I, .L998
    .align 3
 
 .L999:
-   fsqrt.d $f19, $f19
+   fsqrt.s $f19, $f19
+   fmul.s  $f0, $f19, $f0
    move $r4, $r17
-   fcvt.s.d $f0, $f19
    jirl $r0, $r1, 0x0
    .align 3