Ensure random objects' 1-norms are non-zero.

Details: - Fixed an innocuous bug that manifested when running the testsuite on extremely small matrices with randomization via the "powers of 2 in narrow precision range" option enabled. When the randomization function emits a perfect 0.0 to fill a 1x1 matrix, the testsuite will then compute 0.0/0.0 during the normalization process, which leads to NaN residuals. The solution entails smarter implementaions of randv, randnv, randm, and randnm, each of which will compute the 1-norm of the vector or matrix in question. If the object has a 1-norm of 0.0, the object is re-randomized until the 1-norm is not 0.0. Thanks to Kiran Varaganti for reporting this issue (#413). - Updated the implementation of randm_unb_var1() so that it loops over a call to the randv_unb_var1() implementation directly rather than calling it indirectly via randv(). This was done to avoid the overhead of multiple calls to norm1v() when randomizing the rows/columns of a matrix. - Updated comments.
2026-04-20 07:38:53 +00:00 · 2020-06-17 16:42:24 -05:00
parent 35e38fb693
commit b5b604e106
3 changed files with 100 additions and 43 deletions
--- a/frame/include/level0/bli_randnp2s.h
+++ b/frame/include/level0/bli_randnp2s.h
@@ -42,6 +42,8 @@
 { \
 	bli_drandnp2s( a ); \
 }
+
+#if 0
 #define bli_drandnp2s_prev( a ) \
 { \
 	const double m_max  = 3.0; \
@@ -95,6 +97,8 @@
 	   down to float. */ \
 	a = r_val; \
 }
+#endif
+
 #define bli_drandnp2s( a ) \
 { \
 	const double m_max  = 6.0; \
@@ -108,15 +112,19 @@
 	   represents the largest power of two we will use to generate the
 	   random numbers. */ \
 \
-	/* Generate a random real number t on the interval: [0.0, 6.0]. */ \
-	t = ( ( double ) rand() / ( double ) RAND_MAX ) * m_max2; \
+	do \
+	{ \
+		/* Generate a random real number t on the interval: [0.0, 6.0]. */ \
+		t = ( ( double ) rand() / ( double ) RAND_MAX ) * m_max2; \
 \
-	/* Modify t to guarantee that is never equal to the upper bound of
-	   the interval (in this case, 6.0). */ \
-	if ( t == m_max2 ) t = t - 1.0; \
-\
-	/* Transform the interval into the set of integers, {0,1,2,3,4,5}. */ \
-	t = floor( t ); \
+		/* Transform the interval into the set of integers, {0,1,2,3,4,5}.
+		   Note that 6 is prohibited by the loop guard below. */ \
+		t = floor( t ); \
+	} \
+	/* If t is ever equal to m_max2, we re-randomize. The guard against
+	   m_max2 < t is for sanity and shouldn't happen, unless perhaps there
+	   is weirdness in the typecasting to double when computing t above. */ \
+	while ( m_max2 <= t ); \
 \
 	/* Map values of t == 0 to a final value of 0. */ \
 	if ( t == 0.0 ) r_val = 0.0; \
@@ -126,7 +134,7 @@
 \
 		double s_val; \
 \
-		/* Compute r_val = 2^s where s = +/-(t-1) = {-4,-3,-2,-1,0}. */ \
+		/* Compute r_val = 2^s where s = -(t-1) = {-4,-3,-2,-1,0}. */ \
 		r_val = pow( 2.0, -(t - 1.0) ); \
 \
 		/* Compute a random number to determine the sign of the final
--- a/frame/util/bli_util_tapi.c
+++ b/frame/util/bli_util_tapi.c
@@ -271,8 +271,8 @@ void PASTEMAC2(ch,opname,EX_SUF) \
 INSERT_GENTFUNC_BASIC_I( printm, fprintm )


-#undef  GENTFUNC
-#define GENTFUNC( ctype, ch, opname ) \
+#undef  GENTFUNCR
+#define GENTFUNCR( ctype, ctype_r, ch, chr, opname ) \
 \
 void PASTEMAC2(ch,opname,EX_SUF) \
     ( \
@@ -291,23 +291,44 @@ void PASTEMAC2(ch,opname,EX_SUF) \
 	/* Obtain a valid context from the gks if necessary. */ \
 	/*if ( cntx == NULL ) cntx = bli_gks_query_cntx();*/ \
 \
-	/* Invoke the helper variant, which loops over the appropriate kernel
-	   to implement the current operation. */ \
-	PASTEMAC2(ch,opname,_unb_var1) \
-	( \
-	  n, \
-	  x, incx, \
-	  cntx, \
-	  rntm  \
-	); \
+	ctype_r norm; \
+\
+	/* Set the norm to zero. */ \
+	PASTEMAC(chr,set0s)( norm ); \
+\
+	/* Iterate at least once, but continue iterating until the norm is not zero. */ \
+	while ( PASTEMAC(chr,eq0)( norm ) ) \
+	{ \
+		/* Invoke the helper variant, which loops over the appropriate kernel
+		   to implement the current operation. */ \
+		PASTEMAC2(ch,opname,_unb_var1) \
+		( \
+		  n, \
+		  x, incx, \
+		  cntx, \
+		  rntm  \
+		); \
+\
+		/* Check the 1-norm of the randomzied vector. In the unlikely event that
+		   the 1-norm is zero, it means that *all* elements are zero, in which
+		   case we want to re-randomize until the 1-norm is not zero. */ \
+		PASTEMAC2(ch,norm1v,BLIS_TAPI_EX_SUF) \
+		( \
+		  n, \
+		  x, incx, \
+		  &norm, \
+		  cntx, \
+		  rntm  \
+		); \
+	} \
 }

-INSERT_GENTFUNC_BASIC0( randv )
-INSERT_GENTFUNC_BASIC0( randnv )
+INSERT_GENTFUNCR_BASIC0( randv )
+INSERT_GENTFUNCR_BASIC0( randnv )


-#undef  GENTFUNC
-#define GENTFUNC( ctype, ch, opname ) \
+#undef  GENTFUNCR
+#define GENTFUNCR( ctype, ctype_r, ch, chr, opname ) \
 \
 void PASTEMAC2(ch,opname,EX_SUF) \
     ( \
@@ -329,22 +350,47 @@ void PASTEMAC2(ch,opname,EX_SUF) \
 	/* Obtain a valid context from the gks if necessary. */ \
 	/*if ( cntx == NULL ) cntx = bli_gks_query_cntx();*/ \
 \
-	/* Invoke the helper variant, which loops over the appropriate kernel
-	   to implement the current operation. */ \
-	PASTEMAC2(ch,opname,_unb_var1) \
-	( \
-	  diagoffx, \
-	  uplox, \
-	  m, \
-	  n, \
-	  x, rs_x, cs_x, \
-	  cntx, \
-	  rntm  \
-	); \
+	ctype_r norm; \
+\
+	/* Set the norm to zero. */ \
+	PASTEMAC(chr,set0s)( norm ); \
+\
+	/* Iterate at least once, but continue iterating until the norm is not zero. */ \
+	while ( PASTEMAC(chr,eq0)( norm ) ) \
+	{ \
+		/* Invoke the helper variant, which loops over the appropriate kernel
+		   to implement the current operation. */ \
+		PASTEMAC2(ch,opname,_unb_var1) \
+		( \
+		  diagoffx, \
+		  uplox, \
+		  m, \
+		  n, \
+		  x, rs_x, cs_x, \
+		  cntx, \
+		  rntm  \
+		); \
+\
+		/* Check the 1-norm of the randomzied matrix. In the unlikely event that
+		   the 1-norm is zero, it means that *all* elements are zero, in which
+		   case we want to re-randomize until the 1-norm is not zero. */ \
+		PASTEMAC2(ch,norm1m,BLIS_TAPI_EX_SUF) \
+		( \
+		  diagoffx, \
+		  BLIS_NONUNIT_DIAG, \
+		  uplox, \
+		  m, \
+		  n, \
+		  x, rs_x, cs_x, \
+		  &norm, \
+		  cntx, \
+		  rntm  \
+		); \
+	} \
 }

-INSERT_GENTFUNC_BASIC0( randm )
-INSERT_GENTFUNC_BASIC0( randnm )
+INSERT_GENTFUNCR_BASIC0( randm )
+INSERT_GENTFUNCR_BASIC0( randnm )


 #undef  GENTFUNCR
--- a/frame/util/bli_util_unb_var1.c
+++ b/frame/util/bli_util_unb_var1.c
@@ -1019,7 +1019,8 @@ void PASTEMAC(ch,varname) \
 \
 			x1     = x + (j  )*ldx + (0  )*incx; \
 \
-			PASTEMAC2(ch,kername,BLIS_TAPI_EX_SUF) \
+			/*PASTEMAC2(ch,kername,BLIS_TAPI_EX_SUF)*/ \
+			PASTEMAC(ch,kername) \
 			( \
 			  n_elem, \
 			  x1, incx, \
@@ -1046,7 +1047,8 @@ void PASTEMAC(ch,varname) \
 				x0     = x1; \
 				chi1   = x1 + (n_elem-1)*incx; \
 \
-				PASTEMAC2(ch,kername,BLIS_TAPI_EX_SUF) \
+				/*PASTEMAC2(ch,kername,BLIS_TAPI_EX_SUF)*/ \
+				PASTEMAC(ch,kername) \
 				( \
 				  n_elem, \
 				  x1, incx, \
@@ -1086,7 +1088,8 @@ void PASTEMAC(ch,varname) \
 				x2     = x1 + incx; \
 				chi1   = x1; \
 \
-				PASTEMAC2(ch,kername,BLIS_TAPI_EX_SUF) \
+				/*PASTEMAC2(ch,kername,BLIS_TAPI_EX_SUF)*/ \
+				PASTEMAC(ch,kername) \
 				( \
 				  n_elem, \
 				  x1, incx, \
@@ -1118,8 +1121,8 @@ void PASTEMAC(ch,varname) \
 	} \
 }

-INSERT_GENTFUNC_BASIC( randm_unb_var1,  randv )
-INSERT_GENTFUNC_BASIC( randnm_unb_var1, randnv )
+INSERT_GENTFUNC_BASIC( randm_unb_var1,  randv_unb_var1 )
+INSERT_GENTFUNC_BASIC( randnm_unb_var1, randnv_unb_var1 )


 #undef  GENTFUNCR