AOCL LAPACK: dgesvd fix to make singular values positive (flame#97)

Nikku, Deepika · varajago · commit aaaf45b2c4a0 · 2025-11-13T04:48:22.000Z
* AOCL LAPACK: dgesvd fix to make singular values positive

Added code changes for making singular values non negative and
ctest fixes in potri, hetrf

AMD-Internal: CPUPL-7326, CPUPL-7316

Signed-off-by: dnikku &lt;Deepika.Nikku@amd.com&gt;
diff --git a/src/lapack/dec/svd/ext/flamec/lapack_dbdsqr.c b/src/lapack/dec/svd/ext/flamec/lapack_dbdsqr.c
@@ -831,7 +831,11 @@
  for (i__ = 1;
  i__ <= i__1;
  ++i__) {
- if (d__[i__] < 0.) {
+ if (d__[i__] == 0.) {
+ /* Avoid -ZERO */
+ d__[i__] = 0.;
+ }
+ else if (d__[i__] < 0.) {
  d__[i__] = -d__[i__];
  /* Change sign of singular vectors, if desired */
  if (*ncvt > 0) {
diff --git a/src/lapack/dec/svd/ext/flamec/lapack_dbdsqr_small.c b/src/lapack/dec/svd/ext/flamec/lapack_dbdsqr_small.c
@@ -806,7 +806,12 @@ int lapack_dbdsqr_small(char *uplo, integer *n, integer *ncvt, integer *nru,
     i__1 = *n;
     for(i__ = 1; i__ <= i__1; ++i__)
     {
-        if(d__[i__] < 0.)
+        if(d__[i__] == 0.)
+        {
+            /* Avoid -ZERO */
+            d__[i__] = 0.;
+        }
+        else if(d__[i__] < 0.)
         {
             d__[i__] = -d__[i__];
             /* Change sign of singular vectors, if desired */
diff --git a/src/lapack/x86/avx2/fla_dgesvd_nn_small10_avx2.c b/src/lapack/x86/avx2/fla_dgesvd_nn_small10_avx2.c
@@ -11,6 +11,7 @@
 
 #if FLA_ENABLE_AMD_OPT
 
+
 extern void dlartg_(doublereal *da, doublereal *db, doublereal *c__, doublereal *s, doublereal *r);
 extern void dlasq1_(integer *, doublereal *, doublereal *, doublereal *, integer *);
 
@@ -163,13 +164,23 @@ void fla_dgesvd_xx_small10_avx2(integer wntu, integer wntv, integer *m, integer
             u[1 + 2 * *ldu] = -(p3 + p2);
             u[2 + 2 * *ldu] = p0 - p1;
         }
+
+        /* Normalize singular values and scale corresponding vectors for 2x2 case */
+        FLA_NORMALIZE_SINGULAR_VALUE_AND_VECTORS_2X2(1, wntu);
+        FLA_NORMALIZE_SINGULAR_VALUE_AND_VECTORS_2X2(2, wntu);
     }
     else
     {
         if(ncvt == 0 && nru == 0)
         {
             /* Compute Singular Values excluding computation of Singular Vectors */
             dlasq1_(n, &s[1], &e[1], &work[itauq - 1], info);
+
+            /* Ensure singular values are positive */
+            if(*info == 0)
+            {
+                FLA_ENSURE_POSITIVE_SINGULAR_VALUES(*n);
+            }
         }
         else
         {
diff --git a/src/lapack/x86/avx2/fla_dgesvd_small6_avx2.c b/src/lapack/x86/avx2/fla_dgesvd_small6_avx2.c
@@ -122,13 +122,23 @@ void fla_dgesvd_small6_avx2(integer wntus, integer wntvs, integer *m, integer *n
         {
             fla_drot_avx2(&nru, &u[1 + *ldu], &c__1, &u[1 + 2 * *ldu], &c__1, &cosl, &sinl);
         }
+
+        /* Normalize singular values and scale corresponding vectors for 2x2 case */
+        FLA_NORMALIZE_SINGULAR_VALUE_AND_VECTORS_2X2(1, wntus);
+        FLA_NORMALIZE_SINGULAR_VALUE_AND_VECTORS_2X2(2, wntus);
     }
     else
     {
         if(ncvt == 0 && nru == 0)
         {
             /* Compute Singular Values excluding computation of Singular Vectors */
             dlasq1_(n, &s[1], &e[1], &work[itau - 1], info);
+
+            /* Ensure singular values are positive */
+            if(*info == 0)
+            {
+                FLA_ENSURE_POSITIVE_SINGULAR_VALUES(*n);
+            }
         }
         else
         {
diff --git a/src/lapack/x86/avx2/fla_dgesvd_small_avx2.h b/src/lapack/x86/avx2/fla_dgesvd_small_avx2.h
@@ -79,5 +79,38 @@ doublereal d_sign(doublereal *, doublereal *);
     vt[1 + 2 * *ldvt] = scl1 * sr;                   \
     vt[2 + 2 * *ldvt] = scl2 * cr;
 
+/* Macro to normalize singular value sign and
+   scale corresponding singular vectors for 2x2 matrices */
+#define FLA_NORMALIZE_SINGULAR_VALUE_AND_VECTORS_2X2(idx, wntu_var)  \
+    if(s[idx] == 0.0)                                                \
+    {                                                                \
+        s[idx] = 0.0; /* Avoid -ZERO */                              \
+    }                                                                \
+    else if(s[idx] < 0.0)                                            \
+    {                                                                \
+        s[idx] = -s[idx]; /* Make singular value positive */         \
+        if(wntu_var && u != NULL)                                    \
+        {                                                            \
+            /* Negate corresponding left singular vector column */   \
+            u[1 + (idx) * *ldu] = -u[1 + (idx) * *ldu];              \
+            u[2 + (idx) * *ldu] = -u[2 + (idx) * *ldu];              \
+        }                                                            \
+    }
+
+/* Macro to ensure all singular values are positive
+   (values only, no vector adjustments) */
+#define FLA_ENSURE_POSITIVE_SINGULAR_VALUES(n_vals)   \
+    for(integer i__ = 1; i__ <= (n_vals); i__++)      \
+    {                                                 \
+        if(s[i__] == 0.0)                             \
+        {                                             \
+            s[i__] = 0.0; /* Avoid -ZERO */           \
+        }                                             \
+        else if(s[i__] < 0.0)                         \
+        {                                             \
+            s[i__] = -s[i__]; /* Make positive */     \
+        }                                             \
+    }
+
 #endif /* FLA_ENABLE_AMD_OPT */
 #endif /* FLA_DGESVD_SMALL_AVX2_DEFS_H */
diff --git a/test/main/src/test_gesdd.c b/test/main/src/test_gesdd.c
@@ -227,7 +227,6 @@ void fla_test_gesdd_experiment(char *tst_api, test_params_t *params, integer dat
     }
     create_realtype_vector(datatype, &s, fla_min(m, n));
     create_matrix(datatype, LAPACK_COL_MAJOR, m, n, &A_test, lda);
-    create_realtype_vector(datatype, &s_in, fla_min(m, n));
     create_vector(get_realtype(datatype), &scal, 1);
 
     if(!FLA_BRT_VERIFICATION_RUN)
@@ -238,6 +237,7 @@ void fla_test_gesdd_experiment(char *tst_api, test_params_t *params, integer dat
         }
         else
         {
+            create_realtype_vector(datatype, &s_in, fla_min(m, n));
             /* Generate matrix A from known singular values */
             create_svd_matrix(datatype, 'A', m, n, A, lda, s_in, s_one, s_one, i_one, i_one, info);
             if(FLA_OVERFLOW_UNDERFLOW_TEST)
diff --git a/test/main/validate_src/test_common.c b/test/main/validate_src/test_common.c
@@ -3037,7 +3037,7 @@ void print_matrix(char *desc, char *order, integer datatype, integer M, integer
     integer i, j, row_max = M, col_max = N, ldc = lda, ldr = 1;
 
     /* early return */
-    if(M <= 0 || N <= 0 || lda < M || A == NULL || desc == NULL || order == NULL)
+    if(M <= 0 || N <= 0 || lda <= 0 || A == NULL || desc == NULL || order == NULL)
     {
         return;
     }
diff --git a/test/main/validate_src/test_overflow_underflow.c b/test/main/validate_src/test_overflow_underflow.c
@@ -653,7 +653,7 @@ void scale_matrix_overflow_underflow_hetrf(integer datatype, integer m, void *A,
         }
         else if(m < 200)
         {
-            tuning_val = 37.0;
+            tuning_val = 45.0;
         }
         else
         {
@@ -2526,6 +2526,10 @@ void scale_matrix_overflow_underflow_potri(integer datatype, integer n, void *A,
     else if(same_char(imatrix_char, 'O'))
     {
         get_max_from_matrix(datatype, A, max_min, n, n, lda);
+        if (n < 100 && datatype == COMPLEX)
+        {
+            tuning_val = 2.0;
+        }
     }
     calculate_scale_value(datatype, scal, max_min, tuning_val, imatrix_char);
 

Original file line number	Diff line number	Diff line change
`@@ -806,7 +806,12 @@ int lapack_dbdsqr_small(char uplo, integer n, integer ncvt, integer nru,`
`806`	`806`	`i__1 = *n;`
`807`	`807`	`for(i__ = 1; i__ <= i__1; ++i__)`
`808`	`808`	`{`
`809`		`- if(d__[i__] < 0.)`
	`809`	`+ if(d__[i__] == 0.)`
	`810`	`+ {`
	`811`	`+ /* Avoid -ZERO */`
	`812`	`+ d__[i__] = 0.;`
	`813`	`+ }`
	`814`	`+ else if(d__[i__] < 0.)`
`810`	`815`	`{`
`811`	`816`	`d__[i__] = -d__[i__];`
`812`	`817`	`/* Change sign of singular vectors, if desired */`
Original file line number	Diff line number	Diff line change
`@@ -11,6 +11,7 @@`
`11`	`11`
`12`	`12`	`#if FLA_ENABLE_AMD_OPT`
`13`	`13`
	`14`	`+`
`14`	`15`	`extern void dlartg_(doublereal da, doublereal db, doublereal c__, doublereal s, doublereal *r);`
`15`	`16`	`extern void dlasq1_(integer , doublereal , doublereal , doublereal , integer *);`
`16`	`17`
`@@ -163,13 +164,23 @@ void fla_dgesvd_xx_small10_avx2(integer wntu, integer wntv, integer *m, integer`
`163`	`164`	`u[1 + 2 * *ldu] = -(p3 + p2);`
`164`	`165`	`u[2 + 2 * *ldu] = p0 - p1;`
`165`	`166`	`}`
	`167`	`+`
	`168`	`+ /* Normalize singular values and scale corresponding vectors for 2x2 case */`
	`169`	`+ FLA_NORMALIZE_SINGULAR_VALUE_AND_VECTORS_2X2(1, wntu);`
	`170`	`+ FLA_NORMALIZE_SINGULAR_VALUE_AND_VECTORS_2X2(2, wntu);`
`166`	`171`	`}`
`167`	`172`	`else`
`168`	`173`	`{`
`169`	`174`	`if(ncvt == 0 && nru == 0)`
`170`	`175`	`{`
`171`	`176`	`/* Compute Singular Values excluding computation of Singular Vectors */`
`172`	`177`	`dlasq1_(n, &s[1], &e[1], &work[itauq - 1], info);`
	`178`	`+`
	`179`	`+ /* Ensure singular values are positive */`
	`180`	`+ if(*info == 0)`
	`181`	`+ {`
	`182`	`+ FLA_ENSURE_POSITIVE_SINGULAR_VALUES(*n);`
	`183`	`+ }`
`173`	`184`	`}`
`174`	`185`	`else`
`175`	`186`	`{`
Original file line number	Diff line number	Diff line change
`@@ -122,13 +122,23 @@ void fla_dgesvd_small6_avx2(integer wntus, integer wntvs, integer m, integer n`
`122`	`122`	`{`
`123`	`123`	`fla_drot_avx2(&nru, &u[1 + ldu], &c__1, &u[1 + 2 *ldu], &c__1, &cosl, &sinl);`
`124`	`124`	`}`
	`125`	`+`
	`126`	`+ /* Normalize singular values and scale corresponding vectors for 2x2 case */`
	`127`	`+ FLA_NORMALIZE_SINGULAR_VALUE_AND_VECTORS_2X2(1, wntus);`
	`128`	`+ FLA_NORMALIZE_SINGULAR_VALUE_AND_VECTORS_2X2(2, wntus);`
`125`	`129`	`}`
`126`	`130`	`else`
`127`	`131`	`{`
`128`	`132`	`if(ncvt == 0 && nru == 0)`
`129`	`133`	`{`
`130`	`134`	`/* Compute Singular Values excluding computation of Singular Vectors */`
`131`	`135`	`dlasq1_(n, &s[1], &e[1], &work[itau - 1], info);`
	`136`	`+`
	`137`	`+ /* Ensure singular values are positive */`
	`138`	`+ if(*info == 0)`
	`139`	`+ {`
	`140`	`+ FLA_ENSURE_POSITIVE_SINGULAR_VALUES(*n);`
	`141`	`+ }`
`132`	`142`	`}`
`133`	`143`	`else`
`134`	`144`	`{`
Original file line number	Diff line number	Diff line change
`@@ -227,7 +227,6 @@ void fla_test_gesdd_experiment(char tst_api, test_params_t params, integer dat`
`227`	`227`	`}`
`228`	`228`	`create_realtype_vector(datatype, &s, fla_min(m, n));`
`229`	`229`	`create_matrix(datatype, LAPACK_COL_MAJOR, m, n, &A_test, lda);`
`230`		`- create_realtype_vector(datatype, &s_in, fla_min(m, n));`
`231`	`230`	`create_vector(get_realtype(datatype), &scal, 1);`
`232`	`231`
`233`	`232`	`if(!FLA_BRT_VERIFICATION_RUN)`
`@@ -238,6 +237,7 @@ void fla_test_gesdd_experiment(char tst_api, test_params_t params, integer dat`
`238`	`237`	`}`
`239`	`238`	`else`
`240`	`239`	`{`
	`240`	`+ create_realtype_vector(datatype, &s_in, fla_min(m, n));`
`241`	`241`	`/* Generate matrix A from known singular values */`
`242`	`242`	`create_svd_matrix(datatype, 'A', m, n, A, lda, s_in, s_one, s_one, i_one, i_one, info);`
`243`	`243`	`if(FLA_OVERFLOW_UNDERFLOW_TEST)`
Original file line number	Diff line number	Diff line change
`@@ -3037,7 +3037,7 @@ void print_matrix(char desc, char order, integer datatype, integer M, integer`
`3037`	`3037`	`integer i, j, row_max = M, col_max = N, ldc = lda, ldr = 1;`
`3038`	`3038`
`3039`	`3039`	`/* early return */`
`3040`		`- if(M <= 0 \|\| N <= 0 \|\| lda < M \|\| A == NULL \|\| desc == NULL \|\| order == NULL)`
	`3040`	`+ if(M <= 0 \|\| N <= 0 \|\| lda <= 0 \|\| A == NULL \|\| desc == NULL \|\| order == NULL)`
`3041`	`3041`	`{`
`3042`	`3042`	`return;`
`3043`	`3043`	`}`
Original file line number	Diff line number	Diff line change
`@@ -653,7 +653,7 @@ void scale_matrix_overflow_underflow_hetrf(integer datatype, integer m, void *A,`
`653`	`653`	`}`
`654`	`654`	`else if(m < 200)`
`655`	`655`	`{`
`656`		`- tuning_val = 37.0;`
	`656`	`+ tuning_val = 45.0;`
`657`	`657`	`}`
`658`	`658`	`else`
`659`	`659`	`{`
`@@ -2526,6 +2526,10 @@ void scale_matrix_overflow_underflow_potri(integer datatype, integer n, void *A,`
`2526`	`2526`	`else if(same_char(imatrix_char, 'O'))`
`2527`	`2527`	`{`
`2528`	`2528`	`get_max_from_matrix(datatype, A, max_min, n, n, lda);`
	`2529`	`+ if (n < 100 && datatype == COMPLEX)`
	`2530`	`+ {`
	`2531`	`+ tuning_val = 2.0;`
	`2532`	`+ }`
`2529`	`2533`	`}`
`2530`	`2534`	`calculate_scale_value(datatype, scal, max_min, tuning_val, imatrix_char);`
`2531`	`2535`