Add subroutines to introduce diagonal elements in graph structure for fake factorizations

spm_gen_fake_values.c

@@ -34,8 +34,12 @@
  *          Array of size spm->n allocated on entry. On exit, contains the
  *          degree of each vertex in the spm matrix.
  *
+ *******************************************************************************
+ *
+ * @return the number of diagonal elements found during the computation.
+ *
  *******************************************************************************/
-static inline void
+static inline pastix_int_t
 spm_compute_degrees( const pastix_spm_t *spm,
                      pastix_int_t *degrees )
 {
@@ -43,6 +47,7 @@ spm_compute_degrees( const pastix_spm_t *spm,
     pastix_int_t *colptr = spm->colptr;
     pastix_int_t *rowptr = spm->rowptr;
     pastix_int_t baseval;
+    pastix_int_t diagval = 0;
 
     baseval = spmFindBase( spm );
     memset( degrees, 0, spm->n * sizeof(pastix_int_t) );
@@ -65,6 +70,9 @@ spm_compute_degrees( const pastix_spm_t *spm,
                         degrees[i] += 1;
                     }
                 }
+                else {
+                    diagval++;
+                }
             }
         }
         break;
@@ -81,8 +89,125 @@ spm_compute_degrees( const pastix_spm_t *spm,
                     degrees[i] += 1;
                 }
             }
+            else {
+                diagval++;
+            }
         }
     }
+
+    return diagval;
+}
+
+/**
+ *******************************************************************************
+ *
+ * @ingroup spm_dev_driver
+ *
+ * @brief Insert diagonal elements to the graph to have a full Laplacian
+ * generated
+ *
+ *******************************************************************************
+ *
+ * @param[in] spm
+ *          At start, the initial spm structure with missing diagonal elements.
+ *          At exit, contains the same sparse matrix with diagonal elements added.
+ *
+ * @param[in] diagval
+ *          The number of diagonal elements already present in the matrix.
+ *
+ *******************************************************************************/
+static inline void
+spm_add_diag( pastix_spm_t *spm,
+              pastix_int_t  diagval )
+{
+    pastix_spm_t oldspm;
+    pastix_int_t i, j, k;
+    pastix_int_t *oldcol = spm->colptr;
+    pastix_int_t *oldrow = spm->rowptr;
+    pastix_int_t *newrow, *newcol;
+    pastix_int_t baseval;
+
+    baseval = spmFindBase( spm );
+
+    memcpy( &oldspm, spm, sizeof(pastix_spm_t));
+
+    spm->nnz = oldspm.nnz + (spm->n - diagval);
+    newrow = malloc( spm->nnz * sizeof(pastix_int_t) );
+
+    switch(spm->fmttype)
+    {
+    case PastixCSR:
+        /* Swap pointers to call CSC */
+        oldcol = spm->rowptr;
+        oldrow = spm->colptr;
+        spm->colptr = newrow;
+
+    case PastixCSC:
+        newcol = oldcol;
+        if ( spm->fmttype == PastixCSC ) {
+            spm->rowptr = newrow;
+        }
+        diagval = 0;
+        for(j=0; j<spm->n; j++, newcol++) {
+            pastix_int_t nbelt = newcol[1] - newcol[0];
+            int diag = 0;
+
+            memcpy( newrow, oldrow, nbelt * sizeof(pastix_int_t) );
+            newrow += nbelt;
+
+            for(k=0; k<nbelt; k++, oldrow++) {
+                i = *oldrow - baseval;
+
+                if ( i == j ) {
+                    diag = 1;
+                }
+            }
+            newcol[0] += diagval;
+            if ( !diag ) {
+                *newrow = j + baseval;
+                newrow++;
+                diagval++;
+            }
+        }
+        newcol[0] += diagval;
+
+        if ( spm->fmttype == PastixCSC ) {
+            free( oldspm.rowptr );
+        }
+        else {
+            free( oldspm.colptr );
+        }
+        assert( diagval == spm->n );
+        break;
+
+    case PastixIJV:
+        newcol = malloc( spm->nnz * sizeof(pastix_int_t) );
+        spm->colptr = newcol;
+        spm->rowptr = newrow;
+
+        for(k=0; k<spm->n; k++, newrow++, newcol++)
+        {
+            *newrow = k + baseval;
+            *newcol = k + baseval;
+        }
+
+        for(k=0; k<spm->nnz; k++, oldrow++, oldcol++)
+        {
+            i = *oldrow - baseval;
+            j = *oldcol - baseval;
+
+            if ( i == j ) {
+                continue;
+            }
+
+            *newrow = i + baseval;
+            *newcol = j + baseval;
+            newrow++;
+            newcol++;
+        }
+        free( oldspm.colptr );
+        free( oldspm.rowptr );
+    }
 }
 
 /**
@@ -178,9 +303,9 @@ spm_generate_fake_values( pastix_spm_t *spm,
 void
 spmGenFakeValues( pastix_spm_t *spm )
 {
-    pastix_int_t *degrees;
+    pastix_int_t *degrees, diagval;
     double alpha = 10.;
-    double beta = 1.;
+    double beta  = 1.;
 
     assert( spm->flttype == PastixPattern );
     assert( spm->values == NULL );
@@ -216,7 +341,12 @@ spmGenFakeValues( pastix_spm_t *spm )
     }
 
     degrees = malloc( spm->n * sizeof(pastix_int_t));
-    spm_compute_degrees( spm, degrees );
+    diagval = spm_compute_degrees( spm, degrees );
+    if ( diagval != spm->n ) {
+        /* Diagonal elements must be added to the sparse matrix */
+        spm_add_diag( spm, diagval );
+        spmUpdateComputedFields( spm );
+    }
     spm_generate_fake_values( spm, degrees, alpha, beta );
     free( degrees );