Fix {s|c|z}pivotL.c for singular matrices.

xiaoyeli · Nov 26, 2023 · 32832c6 · 32832c6
1 parent 8f4f993
commit 32832c6
Show file tree

Hide file tree

Showing 4 changed files with 25 additions and 26 deletions.
diff --git a/SRC/cpivotL.c b/SRC/cpivotL.c
@@ -134,11 +134,12 @@ if ( jcol == MIN_COL ) {
 
     /* Test for singularity */
     if ( pivmax == 0.0 ) {
-#if 1
+#if 0
+        // There is no valid pivot.
+        // jcol represents the rank of U, 
+        // report the rank, let dgstrf handle the pivot
 	*pivrow = lsub_ptr[pivptr];
 	perm_r[*pivrow] = jcol;
-#else
-	perm_r[diagind] = jcol;
 #endif
 	*usepr = 0;
 	return (jcol+1);

diff --git a/SRC/dpivotL.c b/SRC/dpivotL.c
@@ -122,30 +122,26 @@ if ( jcol == MIN_COL ) {
     diag = EMPTY;
     old_pivptr = nsupc;
     for (isub = nsupc; isub < nsupr; ++isub) {
-	    rtemp = fabs (lu_col_ptr[isub]);
+	rtemp = fabs (lu_col_ptr[isub]);
 	if ( rtemp > pivmax ) {
-	        pivmax = rtemp;
-	        pivptr = isub;
-	    }
-    	if ( *usepr && lsub_ptr[isub] == *pivrow ) old_pivptr = isub;
-	    if ( lsub_ptr[isub] == diagind ) diag = isub;
+	    pivmax = rtemp;
+	    pivptr = isub;
+	}
+	if ( *usepr && lsub_ptr[isub] == *pivrow ) old_pivptr = isub;
+	if ( lsub_ptr[isub] == diagind ) diag = isub;
     }
 
     /* Test for singularity */
     if ( pivmax == 0.0 ) {
 #if 0
-        // There is no valid pivot.  
-        // jcol represents the rank of U
-        // report the rank let dgstrf handle the pivot
-#if 1
-    *pivrow = lsub_ptr[pivptr];
-    perm_r[*pivrow] = jcol;
-#else
-    perm_r[diagind] = jcol;
-#endif
+        // There is no valid pivot.
+        // jcol represents the rank of U, 
+        // report the rank, let dgstrf handle the pivot
+	*pivrow = lsub_ptr[pivptr];
+	perm_r[*pivrow] = jcol;
 #endif
-	    *usepr = 0;
-	    return (jcol+1);
+	*usepr = 0;
+	return (jcol+1);
     }
 
     thresh = u * pivmax;

diff --git a/SRC/spivotL.c b/SRC/spivotL.c
@@ -133,11 +133,12 @@ if ( jcol == MIN_COL ) {
 
     /* Test for singularity */
     if ( pivmax == 0.0 ) {
-#if 1
+#if 0
+        // There is no valid pivot.
+        // jcol represents the rank of U, 
+        // report the rank, let dgstrf handle the pivot
 	*pivrow = lsub_ptr[pivptr];
 	perm_r[*pivrow] = jcol;
-#else
-	perm_r[diagind] = jcol;
 #endif
 	*usepr = 0;
 	return (jcol+1);

diff --git a/SRC/zpivotL.c b/SRC/zpivotL.c
@@ -134,11 +134,12 @@ if ( jcol == MIN_COL ) {
 
     /* Test for singularity */
     if ( pivmax == 0.0 ) {
-#if 1
+#if 0
+        // There is no valid pivot.
+        // jcol represents the rank of U, 
+        // report the rank, let dgstrf handle the pivot
 	*pivrow = lsub_ptr[pivptr];
 	perm_r[*pivrow] = jcol;
-#else
-	perm_r[diagind] = jcol;
 #endif
 	*usepr = 0;
 	return (jcol+1);