Add function to generate compresses matrices by duplicating the values

CMakeLists.txt

@@ -23,6 +23,7 @@ set(SOURCES
   z_spm_convert_to_ijv.c
   z_spm_dofs2flat.c
   z_spm_expand.c
+  z_spm_dof_extend.c
   z_spm_genrhs.c
   z_spm_integer.c
   z_spm_laplacian.c
@@ -47,7 +48,7 @@ set(spm_sources
   spm.c
   spm_io.c
   spm_integers.c
-  spm_expand.c
+  spm_dofs.c
   spm_read_driver.c
   drivers/skitf.f
   drivers/iohb.c

spm.h

@@ -79,12 +79,13 @@ struct pastix_spm_s {
     pastix_int_t      dof;     /**< Number of degrees of freedom per unknown,
                                     if > 0, constant degree of freedom
                                     otherwise, irregular degree of freedom (refer to dofs)      */
-    pastix_int_t     *dofs;    /**< Number of degrees of freedom per unknown (NULL, if dof > 0) */
+    pastix_int_t     *dofs;    /**< Array of the first column of each element in the
+                                    expanded matrix [+baseval]                                  */
     pastix_order_t    layout;  /**< PastixColMajor, or PastixRowMajor                           */
 
-    pastix_int_t     *colptr;  /**< List of indirections to rows for each vertex                */
-    pastix_int_t     *rowptr;  /**< List of edges for each vertex                               */
-    pastix_int_t     *loc2glob;/**< Corresponding numbering from local to global                */
+    pastix_int_t     *colptr;  /**< List of indirections to rows for each vertex [+baseval]     */
+    pastix_int_t     *rowptr;  /**< List of edges for each vertex [+baseval]                    */
+    pastix_int_t     *loc2glob;/**< Corresponding numbering from local to global [+baseval]     */
     void             *values;  /**< Values stored in the matrix                                 */
 };
 

spm_dofs.c

+/**
+ *
+ * @file spm_dofs.c
+ *
+ *  PaStiX spm routines
+ *  PaStiX is a software package provided by Inria Bordeaux - Sud-Ouest,
+ *  LaBRI, University of Bordeaux 1 and IPB.
+ *
+ * @version 5.1.0
+ * @author Mathieu Faverge
+ * @author Theophile Terraz
+ * @author Alban Bellot
+ * @date 2015-01-01
+ *
+ **/
+#include "common.h"
+#include "spm.h"
+
+#include "z_spm.h"
+#include "c_spm.h"
+#include "d_spm.h"
+#include "s_spm.h"
+#include "p_spm.h"
+
+
+pastix_spm_t *
+spmDofExtend( const int type,
+              const int dof,
+              const pastix_spm_t *spm )
+{
+    pastix_spm_t *newspm;
+
+    /* Quick return */
+    if ( dof == 1 )
+        return (pastix_spm_t *)spm;
+
+    if ( spm->dof != 1 ) {
+        pastix_error_print( "Cannot extend spm including dofs already\n" );
+        return (pastix_spm_t *)spm;
+    }
+
+    newspm = spmCopy( spm );
+
+    /**
+     * Generate constant dof
+     */
+    if (type == 0) {
+        newspm->dof     = dof;
+        newspm->nexp    = spm->n  * dof;
+        newspm->gNexp   = spm->gN * dof;
+        newspm->nnzexp  = spm->nnz  * dof * dof;
+        newspm->gnnzexp = spm->gnnz * dof * dof;
+    }
+    else {
+        pastix_int_t i, j, k, dofi, dofj, baseval;
+        pastix_int_t *dofptr, *colptr, *rowptr;
+
+        baseval = spmFindBase( spm );
+
+        newspm->dof  = -1;
+        newspm->dofs = malloc( (spm->n+1) * sizeof(pastix_int_t) );
+        dofptr = newspm->dofs;
+
+        /**
+         * Initialize the dofs array where the degree of freedom of vertex i is
+         * dof[i+1] - dof[i]
+         */
+        *dofptr = 0;
+        for(i=0; i<spm->n; i++, dofptr++) {
+            dofi = 1 + ( rand() % dof );
+            dofptr[1] = dofptr[0] + dofi;
+        }
+
+        newspm->nexp  = *dofptr;
+        newspm->gNexp = *dofptr;
+
+        /**
+         * Count the number of non zeroes
+         */
+        newspm->nnzexp = 0;
+        colptr = newspm->colptr;
+        rowptr = newspm->rowptr;
+        dofptr = newspm->dofs;
+
+        switch(spm->fmttype)
+        {
+        case PastixCSC:
+            for(j=0; j<newspm->n; j++, colptr++) {
+                dofj = dofptr[j+1] - dofptr[j];
+
+                for(k=colptr[0]; k<colptr[1]; k++, rowptr++) {
+                    i = *rowptr - baseval;
+                    dofi = dofptr[i+1] - dofptr[i];
+
+                    newspm->nnzexp += dofi * dofj;
+                }
+            }
+            break;
+        case PastixCSR:
+            for(i=0; i<newspm->n; i++, rowptr++) {
+                dofi = dofptr[i+1] - dofptr[i];
+
+                for(k=rowptr[0]; k<rowptr[1]; k++, rowptr++) {
+                    j = *colptr - baseval;
+                    dofj = dofptr[j+1] - dofptr[j];
+
+                    newspm->nnzexp += dofi * dofj;
+                }
+            }
+            break;
+        case PastixIJV:
+            for(k=0; k<newspm->nnz; k++, rowptr++, colptr++)
+            {
+                i = *rowptr - baseval;
+                j = *colptr - baseval;
+                dofi = dofptr[i+1] - dofptr[i];
+                dofj = dofptr[j+1] - dofptr[j];
+
+                newspm->nnzexp += dofi * dofj;
+            }
+        }
+        newspm->gnnzexp = newspm->nnzexp;
+    }
+
+    switch (spm->flttype) {
+    case PastixFloat:
+        s_spmDofExtend( newspm );
+        break;
+
+    case PastixDouble:
+        d_spmDofExtend( newspm );
+        break;
+
+    case PastixComplex32:
+        c_spmDofExtend( newspm );
+        break;
+
+    case PastixComplex64:
+        z_spmDofExtend( newspm );
+        break;
+
+    case PastixPattern:
+        ;
+    }
+
+    return newspm;
+}

spm_expand.c

-/**
- *
- * @file spm_expand.c
- *
- *  PaStiX spm routines
- *  PaStiX is a software package provided by Inria Bordeaux - Sud-Ouest,
- *  LaBRI, University of Bordeaux 1 and IPB.
- *
- * @version 5.1.0
- * @author Mathieu Faverge
- * @author Theophile Terraz
- * @author Alban Bellot
- * @date 2015-01-01
- *
- **/
-#include "common.h"
-#include "spm.h"
-
-#include "z_spm.h"
-#include "c_spm.h"
-#include "d_spm.h"
-#include "s_spm.h"
-#include "p_spm.h"
-
-void print_tab_int(pastix_int_t* tab, pastix_int_t size)
-{
-    int i;
-    for(i=0;i<size;i++)
-        printf("%d ",tab[i]);
-    printf("\n");
-}
-
-void print_tab_cpx(double* tab, pastix_int_t size)
-{
-    int i;
-    for(i=0;i<size;i++)
-        printf("%f ",tab[i]);
-    printf("\n");
-}
-
-
-void init_rand()
-{
-  srand(time(NULL));
-}
-
-int randminmax(pastix_int_t min, pastix_int_t max)
-{
-  return min + rand()%(max-min+1);
-}
-
-pastix_int_t* computeCorrespondIndex(const pastix_int_t *dofs, pastix_int_t n)
-{
-    pastix_int_t* dofs_coef=malloc(sizeof(pastix_int_t)*(n+1));
-    dofs_coef[0]=0;
-    int i;
-    for(i=1;i<n+1;i++)
-    {
-        dofs_coef[i]=dofs_coef[i-1]+dofs[i-1];
-    }
-    return dofs_coef;
-}
-
-void genDof(pastix_int_t step,pastix_int_t* dofs,pastix_int_t size)
-{
-  int i;
-  for(i=0;i<size;i++)
-      dofs[i]=i%step+1;
-}
-
-void genRandDof(pastix_int_t min,pastix_int_t max,pastix_int_t* dofs,pastix_int_t size)
-{
-  int i;
-  for(i=0;i<size;i++)
-    dofs[i]=randminmax(min,max);
-}
-
-pastix_int_t* computeEltPerCol(const pastix_spm_t *spm)
-{
-  pastix_int_t col,row,baseval;
-
-  pastix_int_t* coef=malloc(sizeof(pastix_int_t)*spm->n);
-  int i;
-  for(i=0;i<spm->n;i++)
-      coef[i]=0;
-  baseval=spmFindBase( spm );
-  row=0;
-  for( col=0; col < spm->n; col++ )
-    {
-      for( row=spm->colptr[col]-baseval; row<spm->colptr[col+1]-baseval; row++)
-	{
-            coef[col] += spm->dofs[col] * spm->dofs[spm->rowptr[row]-baseval];
-	}
-    }
-  return coef;
-}
-
-/**
- *******************************************************************************
- *
- * @ingroup spm_internal
- *
- * dofCst - Expand a CSC matrix without dofs into a CSC matrix with constant dofs
- *
- *
- *******************************************************************************
- *
- * @param[in] spm
- *           The spm to expand
- *
- *******************************************************************************/
-void
-dofCst(pastix_spm_t* spm,pastix_int_t dof)
-{
-    int i;
-    pastix_int_t* dofs = malloc((spm->n+1)*sizeof(pastix_int_t));
-
-    spm->dof      = dof;
-    spm->gNexp    = spm->n*dof;
-    spm->nexp     = spm->n*dof;
-    spm->nnzexp   = spm->n*dof*dof;
-    spm->gnnzexp  = spm->n*dof*dof;
-    spm->dofs     = dofs;
-
-    for(i=0;i<spm->n+1;i++)
-        dofs[i]=dof*i;
-
-    spmExpand(spm);
-}
-
-
-/**
- *******************************************************************************
- *
- * @ingroup spm_internal
- *
- * dofVar - Expand a CSC matrix without dofs into a CSC matrix with variable dofs
- *
- *
- *******************************************************************************
- *
- * @param[in] spm
- *           The spm to expand
- *
- *******************************************************************************/
-void
-dofVar(pastix_spm_t* spm)
-{
-    pastix_int_t *dofs=malloc(sizeof(pastix_int_t)*spm->n);
-    pastix_int_t i,nnzexp,nexp;
-
-    spm->gNexp    = spm->gN;
-    spm->nexp     = spm->n;
-    spm->gnnzexp  = spm->gnnz;
-    spm->nnzexp   = spm->nnz;
-    spm->dof      = -1;
-
-    genDof(3,dofs,spm->n); //cyclique, pas de 3
-    //genRandDof(/*min*/1,/*max*/5 ,spm->dofs,spm->n); //random dofs, entre 1 et 5
-
-    spm->dofs = dofs;
-    pastix_int_t *coef = computeEltPerCol(spm); // free coef
-
-    nnzexp = 0;
-    nexp   = 0;
-
-    for(i=0;i<spm->n;i++)
-    {
-        nnzexp += coef[i];
-        nexp   += spm->dofs[i];
-    }
-
-    spm->nnzexp  = nnzexp;
-    spm->gnnzexp = nnzexp;
-    spm->gNexp   = nexp;
-    spm->nexp    = nexp;
-
-    pastix_int_t* tab = computeCorrespondIndex(dofs,spm->n);
-
-    spm->dofs   = tab;
-    spmExpand(spm);
-
-    /*
-    print_tab_int(spm->colptr,spm->n);
-    print_tab_int(spm->rowptr,spm->nnz);
-    print_tab_cpx(spm->values,spm->nnzexp);
-    print_tab_int(spm->dofs,spm->n+1);
-     */
-}
-
-

z_spm.h

@@ -66,6 +66,6 @@ void z_spmPrint( const pastix_spm_t *spm );
 
 
 pastix_spm_t *z_spmExpand(const pastix_spm_t *spm);
-void          z_spmDofs2Flat(pastix_spm_t *spm);
+void          z_spmDofExtend(pastix_spm_t *spm);
 
 #endif /* _z_spm_H_ */

z_spm_dof_extend.c

+/**
+ *
+ * @file z_spm_dof_extend.c
+ *
+ *  PaStiX spm routines
+ *  PaStiX is a software package provided by Inria Bordeaux - Sud-Ouest,
+ *  LaBRI, University of Bordeaux 1 and IPB.
+ *
+ * @version 5.1.0
+ * @author Mathieu Faverge
+ * @author Alban Bellot
+ * @date 2015-01-01
+ *
+ * @precisions normal z -> c d s
+ **/
+#include "common.h"
+#include "spm.h"
+#include "z_spm.h"
+
+void
+z_spmDofExtend(pastix_spm_t *spm)
+{
+    pastix_int_t        i, j, k, ii, jj, dofi, dofj, baseval;
+    pastix_int_t       *colptr, *rowptr, *dofs;
+    pastix_complex64_t *newval, *oldval;
+
+    oldval = (pastix_complex64_t*)(spm->values);
+    newval = spm->values = malloc( spm->nnzexp * sizeof(pastix_complex64_t) );
+
+    baseval = spmFindBase( spm );
+    colptr = spm->colptr;
+    rowptr = spm->rowptr;
+    dofs   = spm->dofs;
+
+    switch(spm->fmttype)
+    {
+    case PastixCSC:
+        /**
+         * Loop on col
+         */
+        for(j=0; j<spm->n; j++, colptr++)
+        {
+            dofj = ( spm->dof > 0 ) ? spm->dof : dofs[j+1] - dofs[j];
+
+            /**
+             * Loop on rows
+             */
+            for(k=colptr[0]; k<colptr[1]; k++, rowptr++, oldval++)
+            {
+                i = *rowptr - baseval;
+                dofi = ( spm->dof > 0 ) ? spm->dof : dofs[i+1] - dofs[i];
+
+                for(jj=0; jj<dofj; jj++)
+                {
+                    for(ii=0; ii<dofi; ii++, newval++)
+                    {
+                        *newval = *oldval;
+                    }
+                }
+            }
+        }
+        break;
+    case PastixCSR:
+        /**
+         * Loop on row
+         */
+        for(i=0; i<spm->n; i++, rowptr++)
+        {
+            dofi = ( spm->dof > 0 ) ? spm->dof : dofs[i+1] - dofs[i];
+
+            /**
+             * Loop on cols
+             */
+            for(k=rowptr[0]; k<rowptr[1]; k++, colptr++, oldval++)
+            {
+                j = *colptr - baseval;
+                dofj = ( spm->dof > 0 ) ? spm->dof : dofs[j+1] - dofs[j];
+
+                for(jj=0; jj<dofj; jj++)
+                {
+                    for(ii=0; ii<dofi; ii++, newval++)
+                    {
+                        *newval = *oldval;
+                    }
+                }
+            }
+        }
+        break;
+    case PastixIJV:
+        /**
+         * Loop on coordinates
+         */
+        for(k=0; k<spm->nnz; k++, rowptr++, colptr++, oldval++)
+        {
+            i = *rowptr - baseval;
+            j = *colptr - baseval;
+            dofi = ( spm->dof > 0 ) ? spm->dof : dofs[i+1] - dofs[i];
+            dofj = ( spm->dof > 0 ) ? spm->dof : dofs[j+1] - dofs[j];
+
+            for(jj=0; jj<dofj; jj++)
+            {
+                for(ii=0; ii<dofi; ii++, newval++)
+                {
+                    *newval = *oldval;
+                }
+            }
+        }
+        break;
+    }
+    return;
+}