3 - SPM/Dist_basic_routines

Make sure some of the basic routines are working with distributed matrices, or prevent issue by returning an error:

• SpmPrint
• Spm2Dense
• SpmConvert
• SpmGenFakeValues

Some of these routines have also been largely factorized to try to simplify them.

src/z_spm_2dense.c

@@ -19,6 +19,229 @@
 #include "common.h"
 #include "z_spm.h"
 
+typedef spm_complex64_t (*__conj_fct_t)( spm_complex64_t );
+
+static inline spm_complex64_t
+__fct_id( spm_complex64_t val ) {
+    return val;
+}
+
+#if defined(PRECISION_c) || defined(PRECISION_z)
+static inline spm_complex64_t
+__fct_conj( spm_complex64_t val ) {
+    return conj( val );
+}
+#endif
+
+static inline void
+z_spm_2dense_elt_sym_diag_col( spm_complex64_t       *A,
+                               const spm_int_t        lda,
+                               const spm_int_t        row,
+                               const spm_int_t        dofi,
+                               const spm_int_t        col,
+                               const spm_int_t        dofj,
+                               const __conj_fct_t     conjfct,
+                               const spm_complex64_t *valptr )
+{
+    spm_int_t ii, jj;
+
+    for(jj=0; jj<dofj; jj++)
+    {
+        /* Skip unused upper triangular part */
+        for(ii=0; ii<jj; ii++) {
+            valptr++;
+        }
+
+        /* Diagonal element */
+        A[ lda * (col + jj) + (row + ii) ] = *valptr;
+        valptr++;
+
+        for(ii=jj+1; ii<dofi; ii++, valptr++)
+        {
+            /* Lower part */
+            A[ lda * (col + jj) + (row + ii) ] = *valptr;
+            /* Upper part */
+            A[ lda * (row + ii) + (col + jj) ] = conjfct(*valptr);
+        }
+    }
+    (void)conjfct;
+}
+
+static inline void
+z_spm_2dense_elt_sym_diag_row( spm_complex64_t       *A,
+                               const spm_int_t        lda,
+                               const spm_int_t        row,
+                               const spm_int_t        dofi,
+                               const spm_int_t        col,
+                               const spm_int_t        dofj,
+                               const __conj_fct_t     conjfct,
+                               const spm_complex64_t *valptr )
+{
+    spm_int_t ii, jj;
+
+    for(ii=0; ii<dofi; ii++)
+    {
+        for(jj=0; jj<ii; jj++, valptr++)
+        {
+            /* Lower part */
+            A[ lda * (col + jj) + (row + ii) ] = *valptr;
+
+            /* Upper part */
+            A[ lda * (row + ii) + (col + jj) ] = conjfct(*valptr);
+        }
+
+        /* Diagonal element */
+        A[ lda * (col + jj) + (row + ii) ] = *valptr;
+        valptr++;
+
+        /* Skip unused upper triangular part */
+        for(jj=ii+1; jj<dofj; jj++) {
+            valptr++;
+        }
+    }
+    (void)conjfct;
+}
+
+static inline void
+z_spm_2dense_elt_sym_diag( spm_complex64_t       *A,
+                           const spm_int_t        lda,
+                           const spm_layout_t     layout,
+                           const spm_int_t        row,
+                           const spm_int_t        dofi,
+                           const spm_int_t        col,
+                           const spm_int_t        dofj,
+                           const __conj_fct_t     conjfct,
+                           const spm_complex64_t *valptr )
+{
+    if ( layout == SpmColMajor ) {
+        z_spm_2dense_elt_sym_diag_col( A, lda, row, dofi, col, dofj, conjfct, valptr );
+    }
+    else {
+        z_spm_2dense_elt_sym_diag_row( A, lda, row, dofi, col, dofj, conjfct, valptr );
+    }
+    (void)conjfct;
+}
+
+static inline void
+z_spm_2dense_elt_gen_col( spm_complex64_t       *A,
+                          const spm_int_t        lda,
+                          const spm_int_t        row,
+                          const spm_int_t        dofi,
+                          const spm_int_t        col,
+                          const spm_int_t        dofj,
+                          const __conj_fct_t     conjfct,
+                          const spm_complex64_t *valptr )
+{
+    spm_int_t ii, jj;
+
+    for(jj=0; jj<dofj; jj++)
+    {
+        for(ii=0; ii<dofi; ii++, valptr++)
+        {
+            A[ lda * (col + jj) + (row + ii) ] = conjfct(*valptr);
+        }
+    }
+    (void)conjfct;
+}
+
+static inline void
+z_spm_2dense_elt_gen_row( spm_complex64_t       *A,
+                          const spm_int_t        lda,
+                          const spm_int_t        row,
+                          const spm_int_t        dofi,
+                          const spm_int_t        col,
+                          const spm_int_t        dofj,
+                          const __conj_fct_t     conjfct,
+                          const spm_complex64_t *valptr )
+{
+    spm_int_t ii, jj;
+
+    for(ii=0; ii<dofi; ii++)
+    {
+        for(jj=0; jj<dofj; jj++, valptr++)
+        {
+            A[ lda * (col + jj) + (row + ii) ] = conjfct(*valptr);
+        }
+    }
+    (void)conjfct;
+}
+
+static inline void
+z_spm_2dense_elt_gen( spm_complex64_t       *A,
+                      const spm_int_t        lda,
+                      const spm_layout_t     layout,
+                      const spm_int_t        row,
+                      const spm_int_t        dofi,
+                      const spm_int_t        col,
+                      const spm_int_t        dofj,
+                      const __conj_fct_t     conjfct,
+                      const spm_complex64_t *valptr )
+{
+    if ( layout == SpmColMajor ) {
+        z_spm_2dense_elt_gen_col( A, lda, row, dofi, col, dofj, conjfct, valptr );
+    }
+    else {
+        z_spm_2dense_elt_gen_row( A, lda, row, dofi, col, dofj, conjfct, valptr );
+    }
+}
+
+static inline void
+z_spm_2dense_elt_sym_offd( spm_complex64_t       *A,
+                           const spm_int_t        lda,
+                           const spm_layout_t     layout,
+                           const spm_int_t        row,
+                           const spm_int_t        dofi,
+                           const spm_int_t        col,
+                           const spm_int_t        dofj,
+                           const __conj_fct_t     conjfct,
+                           const spm_complex64_t *valptr )
+{
+    if ( layout == SpmColMajor ) {
+        z_spm_2dense_elt_gen_col( A, lda, row, dofi, col, dofj, __fct_id, valptr );
+        z_spm_2dense_elt_gen_row( A, lda, col, dofj, row, dofi, conjfct,  valptr );
+    }
+    else {
+        z_spm_2dense_elt_gen_row( A, lda, row, dofi, col, dofj, __fct_id, valptr );
+        z_spm_2dense_elt_gen_col( A, lda, col, dofj, row, dofi, conjfct,  valptr );
+    }
+}
+
+static inline void
+z_spm_2dense_elt( spm_complex64_t       *A,
+                  const spm_int_t        lda,
+                  const spm_mtxtype_t    mtxtype,
+                  const spm_layout_t     layout,
+                  const spm_int_t        row,
+                  const spm_int_t        dofi,
+                  const spm_int_t        col,
+                  const spm_int_t        dofj,
+                  const spm_complex64_t *valptr )
+{
+    if ( mtxtype == SpmGeneral ) {
+        z_spm_2dense_elt_gen( A, lda, layout, row, dofi, col, dofj, __fct_id, valptr );
+    }
+    else {
+        __conj_fct_t conjfct;
+
+#if defined(PRECISION_c) || defined(PRECISION_z)
+        if ( mtxtype == SpmHermitian ) {
+            conjfct = __fct_conj;
+        }
+        else
+#endif
+        {
+            conjfct = __fct_id;
+        }
+
+        if ( row == col ) {
+            z_spm_2dense_elt_sym_diag( A, lda, layout, row, dofi, col, dofj, conjfct, valptr );
+        }
+        else {
+            z_spm_2dense_elt_sym_offd( A, lda, layout, row, dofi, col, dofj, conjfct, valptr );
+        }
+    }
+}
+
 /**
  *******************************************************************************
  *
@@ -40,12 +263,18 @@
  * @return A dense matrix in Lapack layout format
  *
  *******************************************************************************/
-spm_complex64_t *
+static inline spm_complex64_t *
 z_spmCSC2dense( const spmatrix_t *spm )
 {
-    spm_int_t i, j, k, lda, baseval;
-    spm_complex64_t *A, *valptr;
-    spm_int_t *colptr, *rowptr;
+    spm_int_t              j, k, lda, baseval;
+    spm_int_t              ig, dofi, row;
+    spm_int_t              jg, dofj, col;
+    const spm_int_t       *colptr;
+    const spm_int_t       *rowptr;
+    const spm_int_t       *dofs;
+    const spm_int_t       *loc2glob;
+    const spm_complex64_t *valptr;
+    spm_complex64_t       *A;
 
     assert( spm->fmttype == SpmCSC );
     assert( spm->flttype == SpmComplex64 );
@@ -55,146 +284,43 @@ z_spmCSC2dense( const spmatrix_t *spm )
     memset( A, 0, lda * lda * sizeof(spm_complex64_t));
 
     baseval = spmFindBase( spm );
-    i = 0;
-    j = 0;
 
-    colptr = spm->colptr;
-    rowptr = spm->rowptr;
-    valptr = (spm_complex64_t*)(spm->values);
+    colptr   = spm->colptr;
+    rowptr   = spm->rowptr;
+    valptr   = (spm_complex64_t*)(spm->values);
+    dofs     = spm->dofs;
+    loc2glob = spm->loc2glob;
 
-    /**
-     * Constant degree of fredom of 1
-     */
-    if ( spm->dof == 1 ) {
-        switch( spm->mtxtype ){
-#if defined(PRECISION_z) || defined(PRECISION_c)
-        case SpmHermitian:
-            for(j=0; j<spm->n; j++, colptr++)
-            {
-                for(k=colptr[0]; k<colptr[1]; k++, rowptr++, valptr++)
-                {
-                    i = (*rowptr-baseval);
-                    if( i == j ) {
-                        /* Make sure the matrix is hermitian */
-                        A[ j * lda + i ] = creal(*valptr) + I * 0.;
-                    }
-                    else {
-                        A[ j * lda + i ] = *valptr;
-                        A[ i * lda + j ] = conj(*valptr);
-                    }
-                }
-            }
-            break;
-#endif
-        case SpmSymmetric:
-            for(j=0; j<spm->n; j++, colptr++)
-            {
-                for(k=colptr[0]; k<colptr[1]; k++, rowptr++, valptr++)
-                {
-                    i = (*rowptr-baseval);
-                    A[ j * lda + i ] = *valptr;
-                    A[ i * lda + j ] = *valptr;
-                }
-            }
-            break;
-        case SpmGeneral:
-        default:
-            for(j=0; j<spm->n; j++, colptr++)
-            {
-                for(k=colptr[0]; k<colptr[1]; k++, rowptr++, valptr++)
-                {
-                    i = (*rowptr-baseval);
-                    A[ j * lda + i ] = *valptr;
-                }
-            }
+    for(j=0; j<spm->n; j++, colptr++, loc2glob++)
+    {
+        jg = (spm->loc2glob == NULL) ? j : (*loc2glob) - baseval;
+        if ( spm->dof > 0 ) {
+            dofj = spm->dof;
+            col  = spm->dof * jg;
         }
-    }
-    /**
-     * General degree of freedom (constant or variable)
-     */
-    else {
-        spm_int_t  k, ii, jj, dofi, dofj, col, row;
-        spm_int_t *dofs = spm->dofs;
-
-        switch( spm->mtxtype ){
-#if defined(PRECISION_z) || defined(PRECISION_c)
-        case SpmHermitian:
-            for(j=0; j<spm->n; j++, colptr++)
-            {
-                dofj = ( spm->dof > 1 ) ?  spm->dof      : dofs[j+1] - dofs[j];
-                col  = ( spm->dof > 1 ) ? (spm->dof * j) : dofs[j] - baseval;
-
-                for(k=colptr[0]; k<colptr[1]; k++, rowptr++)
-                {
-                    i = (*rowptr - baseval);
-                    dofi = ( spm->dof > 1 ) ?  spm->dof      : dofs[i+1] - dofs[i];
-                    row  = ( spm->dof > 1 ) ? (spm->dof * i) : dofs[i] - baseval;
-
-                    for(jj=0; jj<dofj; jj++)
-                    {
-                        for(ii=0; ii<dofi; ii++, valptr++)
-                        {
-                            if( col+jj == row+ii ) {
-                                /* Make sure the matrix is hermitian */
-                                A[ (col + jj) * lda + (row + ii) ] = creal(*valptr) + I * 0.;
-                            }
-                            else {
-                                A[ (col + jj) * lda + (row + ii) ] = *valptr;
-                                A[ (row + ii) * lda + (col + jj) ] = conj(*valptr);
-                            }
-                        }
-                    }
-                }
-            }
-            break;
-#endif
-        case SpmSymmetric:
-            for(j=0; j<spm->n; j++, colptr++)
-            {
-                dofj = ( spm->dof > 1 ) ? spm->dof     : dofs[j+1] - dofs[j];
-                col  = ( spm->dof > 1 ) ? spm->dof * j : dofs[j] - baseval;
-
-                for(k=colptr[0]; k<colptr[1]; k++, rowptr++)
-                {
-                    i = (*rowptr - baseval);
-                    dofi = ( spm->dof > 1 ) ? spm->dof     : dofs[i+1] - dofs[i];
-                    row  = ( spm->dof > 1 ) ? spm->dof * i : dofs[i] - baseval;
-
-                    for(jj=0; jj<dofj; jj++)
-                    {
-                        for(ii=0; ii<dofi; ii++, valptr++)
-                        {
-                            A[ (col + jj) * lda + (row + ii) ] = *valptr;
-                            A[ (row + ii) * lda + (col + jj) ] = *valptr;
-                        }
-                    }
-                }
+        else {
+            dofj = dofs[jg+1] - dofs[jg];
+            col  = dofs[jg] - baseval;
+        }
+
+        for(k=colptr[0]; k<colptr[1]; k++, rowptr++)
+        {
+            ig = (*rowptr - baseval);
+            if ( spm->dof > 0 ) {
+                dofi = spm->dof;
+                row  = spm->dof * ig;
             }
-            break;
-        case SpmGeneral:
-        default:
-            for(j=0; j<spm->n; j++, colptr++)
-            {
-                dofj = ( spm->dof > 1 ) ? spm->dof     : dofs[j+1] - dofs[j];
-                col  = ( spm->dof > 1 ) ? spm->dof * j : dofs[j] - baseval;
-
-                for(k=colptr[0]; k<colptr[1]; k++, rowptr++)
-                {
-                    i = (*rowptr - baseval);
-                    dofi = ( spm->dof > 1 ) ? spm->dof     : dofs[i+1] - dofs[i];
-                    row  = ( spm->dof > 1 ) ? spm->dof * i : dofs[i] - baseval;
-
-                    for(jj=0; jj<dofj; jj++)
-                    {
-                        for(ii=0; ii<dofi; ii++, valptr++)
-                        {
-                            A[ (col + jj) * lda + (row + ii) ] = *valptr;
-                        }
-                    }
-                }
+            else {
+                dofi = dofs[ig+1] - dofs[ig];
+                row  = dofs[ig] - baseval;
             }
+
+            z_spm_2dense_elt( A, lda, spm->mtxtype, spm->layout,
+                              row, dofi, col, dofj, valptr );
+            valptr += dofi * dofj;
         }
     }
+
     return A;
 }
 
@@ -219,12 +345,18 @@ z_spmCSC2dense( const spmatrix_t *spm )
  * @return A dense matrix in Lapack layout format
  *
  *******************************************************************************/
-spm_complex64_t *
+static inline spm_complex64_t *
 z_spmCSR2dense( const spmatrix_t *spm )
 {
-    spm_int_t i, j, k, lda, baseval;
-    spm_complex64_t *A, *valptr;
-    spm_int_t *colptr, *rowptr;
+    spm_int_t              i, k, lda, baseval;
+    spm_int_t              ig, dofi, row;
+    spm_int_t              jg, dofj, col;
+    const spm_int_t       *colptr;
+    const spm_int_t       *rowptr;
+    const spm_int_t       *dofs;
+    const spm_int_t       *loc2glob;
+    const spm_complex64_t *valptr;
+    spm_complex64_t       *A;
 
     assert( spm->fmttype == SpmCSR );
     assert( spm->flttype == SpmComplex64 );
@@ -234,146 +366,43 @@ z_spmCSR2dense( const spmatrix_t *spm )
     memset( A, 0, lda * lda * sizeof(spm_complex64_t));
 
     baseval = spmFindBase( spm );
-    i = 0;
-    j = 0;
 
-    colptr = spm->colptr;
-    rowptr = spm->rowptr;
-    valptr = (spm_complex64_t*)(spm->values);
+    colptr   = spm->colptr;
+    rowptr   = spm->rowptr;
+    valptr   = (spm_complex64_t*)(spm->values);
+    dofs     = spm->dofs;
+    loc2glob = spm->loc2glob;
 
-    /**
-     * Constant degree of fredom of 1
-     */
-    if ( spm->dof == 1 ) {
-        switch( spm->mtxtype ){
-#if defined(PRECISION_z) || defined(PRECISION_c)
-        case SpmHermitian:
-            for(i=0; i<spm->n; i++, rowptr++)
-            {
-                for(k=rowptr[0]; k<rowptr[1]; k++, colptr++, valptr++)
-                {
-                    j = (*colptr-baseval);
-                    if( i == j ) {
-                        /* Make sure the matrix is hermitian */
-                        A[ j * lda + i ] = creal(*valptr) + I * 0.;
-                    }
-                    else {
-                        A[ j * lda + i ] = *valptr;
-                        A[ i * lda + j ] = conj(*valptr);
-                    }
-                }
-            }
-            break;
-#endif
-        case SpmSymmetric:
-            for(i=0; i<spm->n; i++, rowptr++)
-            {
-                for(k=rowptr[0]; k<rowptr[1]; k++, colptr++, valptr++)
-                {
-                    j = (*colptr-baseval);
-                    A[ j * lda + i ] = *valptr;
-                    A[ i * lda + j ] = *valptr;
-                }
-            }
-            break;
-        case SpmGeneral:
-        default:
-            for(i=0; i<spm->n; i++, rowptr++)
-            {
-                for(k=rowptr[0]; k<rowptr[1]; k++, colptr++, valptr++)
-                {
-                    j = (*colptr-baseval);
-                    A[ j * lda + i ] = *valptr;
-                }
-            }
+    for(i=0; i<spm->n; i++, rowptr++, loc2glob++)
+    {
+        ig = (spm->loc2glob == NULL) ? i : (*loc2glob) - baseval;
+        if ( spm->dof > 0 ) {
+            dofi = spm->dof;
+            row  = spm->dof * ig;
         }
-    }
-    /**
-     * General degree of freedom (constant or variable)
-     */
-    else {
-        spm_int_t  k, ii, jj, dofi, dofj, col, row;
-        spm_int_t *dofs = spm->dofs;
-
-        switch( spm->mtxtype ){
-#if defined(PRECISION_z) || defined(PRECISION_c)
-        case SpmHermitian:
-            for(i=0; i<spm->n; i++, rowptr++)
-            {
-                dofi = ( spm->dof > 1 ) ? spm->dof     : dofs[i+1] - dofs[i];
-                row  = ( spm->dof > 1 ) ? spm->dof * i : dofs[i] - baseval;
-
-                for(k=rowptr[0]; k<rowptr[1]; k++, colptr++)
-                {
-                    j = (*colptr - baseval);
-                    dofj = ( spm->dof > 1 ) ? spm->dof     : dofs[j+1] - dofs[j];
-                    col  = ( spm->dof > 1 ) ? spm->dof * j : dofs[j] - baseval;
-
-                    for(jj=0; jj<dofj; jj++)
-                    {
-                        for(ii=0; ii<dofi; ii++, valptr++)
-                        {
-                            if( col+jj == row+ii ) {
-                                /* Make sure the matrix is hermitian */
-                                A[ (col + jj) * lda + (row + ii) ] = creal(*valptr) + I * 0.;
-                            }
-                            else {
-                                A[ (col + jj) * lda + (row + ii) ] = *valptr;
-                                A[ (row + ii) * lda + (col + jj) ] = conj(*valptr);
-                            }
-                        }
-                    }
-                }
-            }
-            break;
-#endif
-        case SpmSymmetric:
-            for(i=0; i<spm->n; i++, rowptr++)
-            {
-                dofi = ( spm->dof > 1 ) ? spm->dof     : dofs[i+1] - dofs[i];
-                row  = ( spm->dof > 1 ) ? spm->dof * i : dofs[i] - baseval;
-
-                for(k=rowptr[0]; k<rowptr[1]; k++, colptr++)
-                {
-                    j = (*colptr - baseval);
-                    dofj = ( spm->dof > 1 ) ? spm->dof     : dofs[j+1] - dofs[j];
-                    col  = ( spm->dof > 1 ) ? spm->dof * j : dofs[j] - baseval;
-
-                    for(jj=0; jj<dofj; jj++)
-                    {
-                        for(ii=0; ii<dofi; ii++, valptr++)
-                        {
-                            A[ (col + jj) * lda + (row + ii) ] = *valptr;
-                            A[ (row + ii) * lda + (col + jj) ] = *valptr;
-                        }
-                    }
-                }
+        else {
+            dofi = dofs[ig+1] - dofs[ig];
+            row  = dofs[ig] - baseval;
+        }
+
+        for(k=rowptr[0]; k<rowptr[1]; k++, colptr++)
+        {
+            jg = (*colptr - baseval);
+            if ( spm->dof > 0 ) {
+                dofj = spm->dof;
+                col  = spm->dof * jg;
             }
-            break;
-        case SpmGeneral:
-        default:
-            for(i=0; i<spm->n; i++, rowptr++)
-            {
-                dofi = ( spm->dof > 1 ) ? spm->dof     : dofs[i+1] - dofs[i];
-                row  = ( spm->dof > 1 ) ? spm->dof * i : dofs[i] - baseval;
-
-                for(k=rowptr[0]; k<rowptr[1]; k++, colptr++)
-                {
-                    j = (*colptr - baseval);
-                    dofj = ( spm->dof > 1 ) ? spm->dof     : dofs[j+1] - dofs[j];
-                    col  = ( spm->dof > 1 ) ? spm->dof * j : dofs[j] - baseval;
-
-                    for(jj=0; jj<dofj; jj++)
-                    {
-                        for(ii=0; ii<dofi; ii++, valptr++)
-                        {
-                            A[ (col + jj) * lda + (row + ii) ] = *valptr;
-                        }
-                    }
-                }
+            else {
+                dofj = dofs[jg+1] - dofs[jg];
+                col  = dofs[jg] - baseval;
             }
+
+            z_spm_2dense_elt( A, lda, spm->mtxtype, spm->layout,
+                              row, dofi, col, dofj, valptr );
+            valptr += dofi * dofj;
         }
     }
+
     return A;
 }
 
@@ -398,12 +427,17 @@ z_spmCSR2dense( const spmatrix_t *spm )
  * @return A dense matrix in Lapack layout format
  *
  *******************************************************************************/
-spm_complex64_t *
+static inline spm_complex64_t *
 z_spmIJV2dense( const spmatrix_t *spm )
 {
-    spm_int_t i, j, k, lda, baseval;
-    spm_complex64_t *A, *valptr;
-    spm_int_t *colptr, *rowptr;
+    spm_int_t              k, lda, baseval;
+    spm_int_t              i, dofi, row;
+    spm_int_t              j, dofj, col;
+    const spm_int_t       *colptr;
+    const spm_int_t       *rowptr;
+    const spm_int_t       *dofs;
+    const spm_complex64_t *valptr;
+    spm_complex64_t       *A;
 
     assert( spm->fmttype == SpmIJV );
     assert( spm->flttype == SpmComplex64 );
@@ -413,162 +447,35 @@ z_spmIJV2dense( const spmatrix_t *spm )
     memset( A, 0, lda * lda * sizeof(spm_complex64_t));
 
     baseval = spmFindBase( spm );
-    i = 0;
-    j = 0;
 
     colptr = spm->colptr;
     rowptr = spm->rowptr;
     valptr = (spm_complex64_t*)(spm->values);
+    dofs   = spm->dofs;
 
-    /**
-     * Constant degree of fredom of 1
-     */
-    if ( spm->dof == 1 ) {
-        switch( spm->mtxtype ){
-#if defined(PRECISION_z) || defined(PRECISION_c)
-        case SpmHermitian:
-            for(k=0; k<spm->nnz; k++, rowptr++, colptr++, valptr++)
-            {
-                i = *rowptr - baseval;
-                j = *colptr - baseval;
-
-                if( i == j ) {
-                    /* Make sure the matrix is hermitian */
-                    A[ i * lda + i ] = creal(*valptr) + I * 0.;
-                }
-                else {
-                    A[ j * lda + i ] = *valptr;
-                    A[ i * lda + j ] = conj(*valptr);
-                }
-            }
-            break;
-#endif
-        case SpmSymmetric:
-            for(k=0; k<spm->nnz; k++, rowptr++, colptr++, valptr++)
-            {
-                i = *rowptr - baseval;
-                j = *colptr - baseval;
-
-                A[ j * lda + i ] = *valptr;
-                A[ i * lda + j ] = *valptr;
-            }
-            break;
-        case SpmGeneral:
-        default:
-            for(k=0; k<spm->nnz; k++, rowptr++, colptr++, valptr++)
-            {
-                i = *rowptr - baseval;
-                j = *colptr - baseval;
-
-                A[ j * lda + i ] = *valptr;
-            }
+    for(k=0; k<spm->nnz; k++, rowptr++, colptr++)
+    {
+        i = *rowptr - baseval;
+        j = *colptr - baseval;
+
+        if ( spm->dof > 0 ) {
+            dofi = spm->dof;
+            row  = spm->dof * i;
+            dofj = spm->dof;
+            col  = spm->dof * j;
         }
-    }
-    /**
-     * General degree of freedom (constant or variable)
-     */
-    else {
-        spm_int_t  k, ii, jj, dofi, dofj, col, row;
-        spm_int_t *dofs = spm->dofs;
-
-        switch( spm->mtxtype ){
-#if defined(PRECISION_z) || defined(PRECISION_c)
-        case SpmHermitian:
-            for(k=0; k<spm->nnz; k++, rowptr++, colptr++)
-            {
-                i = *rowptr - baseval;
-                j = *colptr - baseval;
-
-                if ( spm->dof > 1 ) {
-                    dofi = spm->dof;
-                    row  = spm->dof * i;
-                    dofj = spm->dof;
-                    col  = spm->dof * j;
-                }
-                else {
-                    dofi = dofs[i+1] - dofs[i];
-                    row  = dofs[i] - baseval;
-                    dofj = dofs[j+1] - dofs[j];
-                    col  = dofs[j] - baseval;
-                }
-
-                for(jj=0; jj<dofj; jj++)
-                {
-                    for(ii=0; ii<dofi; ii++, valptr++)
-                    {
-                        if( col+jj == row+ii ) {
-                            /* Make sure the matrix is hermitian */
-                            A[ (col+jj) * lda + (row+ii) ] = creal(*valptr) + I * 0.;
-                        }
-                        else {
-                            A[ (col+jj) * lda + (row+ii) ] = *valptr;
-                            A[ (row+ii) * lda + (col+jj) ] = conj(*valptr);
-                        }
-                    }
-                }
-            }
-            break;
-#endif
-        case SpmSymmetric:
-            for(k=0; k<spm->nnz; k++, rowptr++, colptr++)
-            {
-                i = *rowptr - baseval;
-                j = *colptr - baseval;
-
-                if ( spm->dof > 1 ) {
-                    dofi = spm->dof;
-                    row  = spm->dof * i;
-                    dofj = spm->dof;
-                    col  = spm->dof * j;
-                }
-                else {
-                    dofi = dofs[i+1] - dofs[i];
-                    row  = dofs[i] - baseval;
-                    dofj = dofs[j+1] - dofs[j];
-                    col  = dofs[j] - baseval;
-                }
-
-                for(jj=0; jj<dofj; jj++)
-                {
-                    for(ii=0; ii<dofi; ii++, valptr++)
-                    {
-                        A[ (col+jj) * lda + (row+ii) ] = *valptr;
-                        A[ (row+ii) * lda + (col+jj) ] = *valptr;
-                    }
-                }
-            }
-            break;
-
-        case SpmGeneral:
-        default:
-            for(k=0; k<spm->nnz; k++, rowptr++, colptr++)
-            {
-                i = *rowptr - baseval;
-                j = *colptr - baseval;
-
-                if ( spm->dof > 1 ) {
-                    dofi = spm->dof;
-                    row  = spm->dof * i;
-                    dofj = spm->dof;
-                    col  = spm->dof * j;
-                }
-                else {
-                    dofi = dofs[i+1] - dofs[i];
-                    row  = dofs[i] - baseval;
-                    dofj = dofs[j+1] - dofs[j];
-                    col  = dofs[j] - baseval;
-                }
-
-                for(jj=0; jj<dofj; jj++)
-                {
-                    for(ii=0; ii<dofi; ii++, valptr++)
-                    {
-                        A[ (col+jj) * lda + (row+ii) ] = *valptr;
-                    }
-                }
-            }
+        else {
+            dofi = dofs[i+1] - dofs[i];
+            row  = dofs[i] - baseval;
+            dofj = dofs[j+1] - dofs[j];
+            col  = dofs[j] - baseval;
         }
+
+        z_spm_2dense_elt( A, lda, spm->mtxtype, spm->layout,
+                          row, dofi, col, dofj, valptr );
+        valptr += dofi * dofj;
     }
+
     return A;
 }
 
@@ -596,15 +503,26 @@ z_spmIJV2dense( const spmatrix_t *spm )
 spm_complex64_t *
 z_spm2dense( const spmatrix_t *spm )
 {
+    spm_complex64_t *A;
+
+    if ( spm->loc2glob != NULL ) {
+        fprintf( stderr, "spm2dense: Conversion to dense matrix with distributed spm is not available\n");
+        return NULL;
+    }
+
     switch (spm->fmttype) {
     case SpmCSC:
-        return z_spmCSC2dense( spm );
+        A = z_spmCSC2dense( spm );
+        break;
     case SpmCSR:
-        return z_spmCSR2dense( spm );
+        A = z_spmCSR2dense( spm );
+        break;
     case SpmIJV:
-        return z_spmIJV2dense( spm );
+        A = z_spmIJV2dense( spm );
+        break;
     }
-    return NULL;
+
+    return A;
 }
 
 /**