z_csc_utils.c

/**
 *  PaStiX CSC management routines.
 *
 *  PaStiX is a software package provided by Inria Bordeaux - Sud-Ouest,
 *  LaBRI, University of Bordeaux 1 and IPB.
 *
 * @version 1.0.0
 * @author Mathieu Faverge
 * @author Pierre Ramet
 * @author Xavier Lacoste
 * @date 2011-11-11
 * @precisions normal z -> c d s
 *
 **/
#include "common.h"
#include "z_csc_utils.h"

/* /\* */
/*    Function: cmp_colrow */

/*    Used for qsort to sort arrays of pastix_int_t following their first element. */

/*    Returns the difference between the first element of *p1* and */
/*    the first element of *p2* */

/*    Parameters: */
/*      p1 - the first array to compare */
/*      p2 - the second array to compare */
/* *\/ */
/* int */
/* cmp_colrow(const void *p1, const void *p2) */
/* { */
/*   return ((* (pastix_int_t * const *) p1) - (*(pastix_int_t * const *) p2)); */
/* } */

/*
  Function: z_csc_symgraph


  Modify the CSC to a symetric graph one.
  Don't use it on a lower symetric CSC
  it would give you all the CSC upper + lower.

  External function.

  Parameters:
    n     - Number of columns/vertices
    ia    - Starting index of each column in *ja* and *a*
    ja    - Row index of each element
    a     - Value of each element,can be NULL
    newn  - New number of column
    newia - Starting index of each column in *ja* and *a*
    newja - Row index of each element
    newa  - Value of each element,can be NULL

 */
int z_csc_symgraph( pastix_int_t               n,
                    const pastix_int_t        *ia,
                    const pastix_int_t        *ja,
                    const pastix_complex64_t  *a,
                    pastix_int_t              *newn,
                    pastix_int_t             **newia,
                    pastix_int_t             **newja,
                    pastix_complex64_t       **newa)
{
  return z_csc_symgraph_int(n, ia, ja, a, newn, newia, newja, newa, API_NO);
}

/*
  Function: csc_symgraph_int


  Modify the CSC to a symetric graph one.
  Don't use it on a lower symetric CSC
  it would give you all the CSC upper + lower.

  Parameters:
    n           - Number of columns/vertices
    ia          - Starting index of each column in *ja* and *a*
    ja          - Row index of each element
    a           - Value of each element,can be NULL
    newn        - New number of column
    newia       - Starting index of each column in *ja* and *a*
    newja       - Row index of each element
    newa        - Value of each element,can be NULL
    malloc_flag - flag to indicate if function call is intern to pastix or extern.
 */
int z_csc_symgraph_int (pastix_int_t               n,
                        const pastix_int_t        *ia,
                        const pastix_int_t        *ja,
                        const pastix_complex64_t  *a,
                        pastix_int_t              *newn,
                        pastix_int_t             **newia,
                        pastix_int_t             **newja,
                        pastix_complex64_t       **newa,
                        int                        malloc_flag)
{
  pastix_int_t * nbrEltCol = NULL; /* nbrEltCol[i] = Number of elt to add in column i */
  pastix_int_t * cia       = NULL; /* ia of diff between good CSC and bad CSC */
  pastix_int_t * cja       = NULL; /* ja of diff between good CSC and bad CSC */
  pastix_int_t   nbr2add;          /* Number of elt to add */
  pastix_int_t   itercol, iterrow, iterrow2; /* iterators */
  pastix_int_t   l = ia[n] -1;
  pastix_int_t   newl;

  /* Ncol=Nrow don't need change */
  *newn = n;

  MALLOC_INTERN(nbrEltCol, n, pastix_int_t);
  /* !! Need check for malloc */

  /* Init nbrEltCol */
  for (itercol=0; itercol<n; itercol++)
    {
      nbrEltCol[itercol]=0;
    }

  /* Compute number of element by col to add for correcting the CSC */
  for (itercol=0; itercol<n; itercol++)
    {
      for (iterrow=ia[itercol]-1; iterrow<ia[itercol+1]-1; iterrow++)
        {
          if (ja[iterrow] != (itercol+1))
            {
              /* Not diagonal elt */
              /* So we have a (i,j) and we are looking for a (j,i) elt */
              /* i = itercol+1, j=ja[iterrow] */
              int rowidx=ja[iterrow]-1;
              int flag=0;

              for (iterrow2=ia[rowidx]-1; iterrow2<ia[rowidx+1]-1; iterrow2++)
                {
                  if (ja[iterrow2] == itercol+1)
                    {
                      /* Ok we found (j,i) so stop this madness */
                      flag = 1;
                      break;
                    }
                }

              if (flag==0)
                {
                  /* We never find (j,i) so increase nbrEltCol[j] */
                  (nbrEltCol[ja[iterrow]-1])++;
                }
            }
        }
    }

  /* Compute number of element to add */
  /* And cia=ia part of csc of element to add */
  /* kind of a diff between the corrected one and the original CSC */
  MALLOC_INTERN(cia, n+1, pastix_int_t);
  /* !! Need checking good alloc) */
  nbr2add=0;
  for (itercol=0;itercol<n;itercol++)
    {
      cia[itercol]=nbr2add;
      nbr2add += nbrEltCol[itercol];
    }
  cia[n]=nbr2add;
  /*fprintf(stderr, "nbr of elt to add %ld\n", nbr2add);*/

  if (nbr2add != 0)
    {
      /* Build cja */
      /* like cia, cja is ja part of diff CSC */
      MALLOC_INTERN(cja, nbr2add, pastix_int_t);
      /* !! again we need check of memAlloc */

      /* We walkthrough again the csc */
      for (itercol=0;itercol<n;itercol++)
        {
          for (iterrow=ia[itercol]-1;iterrow<ia[itercol+1]-1;iterrow++)
            {
              if (ja[iterrow] != itercol+1)
                {
                  /* we find (i,j) need to find (j,i) */
                  int rowidx=ja[iterrow]-1;
                  int flag=0;

                  for (iterrow2=ia[rowidx]-1;iterrow2<ia[rowidx+1]-1;iterrow2++)
                    {
                      if (ja[iterrow2] == itercol+1)
                        {
                          /* find (j,i) */
                          flag=1;
                          break;
                        }
                    }

                  if (flag==0)
                    {
                      /* We don't find (j,i) so put in diff CSC (cia,cja,0) */
                      pastix_int_t index=ja[iterrow]-1;
                      /* cia[index] = index to put in cja the elt */
                      cja[cia[index]] = itercol+1;
                      (cia[index])++;
                    }
                }
            }
        }

      /* Restore cia */
      cia[0]=0;
      for (itercol=0;itercol<n;itercol++)
        {
          cia[itercol+1]=cia[itercol]+nbrEltCol[itercol];
        }

      memFree_null(nbrEltCol);

      /* Build corrected csc */
      newl = l+nbr2add;

      if (malloc_flag == API_NO)
        {
          /* Ici on a des malloc car le free est externe */
          MALLOC_EXTERN(*newia, n+1, pastix_int_t);
          MALLOC_EXTERN(*newja, newl, pastix_int_t);
          if (a != NULL)
            MALLOC_EXTERN(*newa, newl, pastix_complex64_t);
        }
      else
        {
          /* Ici on a des memAlloc car le free est interne */
          MALLOC_INTERN(*newia, n+1, pastix_int_t);
          MALLOC_INTERN(*newja, newl, pastix_int_t);
          if (a != NULL)
            MALLOC_INTERN(*newa, newl, pastix_complex64_t);
        }
      iterrow2 = 0; /* iterator of the CSC diff */
      for (itercol=0; itercol<n; itercol++)
        {
          (*newia)[itercol] = ia[itercol]+iterrow2;
          for (iterrow=ia[itercol]-1;iterrow<ia[itercol+1]-1;iterrow++)
            {
              /* we add the new elt with respect of order in row */
              while ((iterrow2<cia[itercol+1]) &&
                     (ja[iterrow] > cja[iterrow2]))
                {
                  /* we have elt(s) to add with a row lower than ja[iterrow] */
                  (*newja)[iterrow+iterrow2]=cja[iterrow2];
                  if (a != NULL)
                    (*newa)[iterrow+iterrow2]=0.;
                  iterrow2++;
                }

              /* Put the elt from the origin CSC */
              (*newja)[iterrow+iterrow2] = ja[iterrow];
              if (a != NULL)
                (*newa)[iterrow+iterrow2] = a[iterrow];
            }

          /* Since we put elt with a row lower than elt in the origin CSC */
          /* We could have some elt to add after the last elt in the column */
          while(iterrow2<cia[itercol+1])
            {
              (*newja)[iterrow+iterrow2]=cja[iterrow2];
              if (a != NULL)
                (*newa)[iterrow+iterrow2]=0.;
              iterrow2++;
            }
        }

      (*newia)[n]=ia[n]+iterrow2;
      memFree_null(cja);

    }
  else
    {
      /* No correction to do */
      memFree_null(nbrEltCol);
      newl = l;
      if (malloc_flag == API_NO)
        {
          /* ici on a des mallocs car le free est externe */
          MALLOC_EXTERN(*newia, n+1, pastix_int_t);
          MALLOC_EXTERN(*newja, l, pastix_int_t);
          if (a != NULL)
            MALLOC_EXTERN(*newa, l, pastix_complex64_t);
        }
      else
        {
          /* ici on a des memAllocs car le free est interne */
          MALLOC_INTERN(*newia, n+1, pastix_int_t);
          MALLOC_INTERN(*newja, l,   pastix_int_t);
          if (a != NULL)
            MALLOC_INTERN(*newa, l, pastix_complex64_t);
        }
      memcpy((*newia), ia, (n+1)*sizeof(pastix_int_t));
      memcpy((*newja), ja, l * sizeof(pastix_int_t));
      if (a != NULL)
        memcpy((*newa) , a , l * sizeof(pastix_complex64_t));
    }
  memFree_null(cia);

  return EXIT_SUCCESS;
}


/**
    Function: z_csc_noDiag

    Supress diagonal term.
    After this call, *ja* can be reallocated to *ia[n] -1*.

    Parameters:
      n  - size of the matrix.
      ia - Index in *ja* and *a* of the first element of each column
      ja - row of each element
      a  - value of each element, can be set to NULL

    Returns:
      ia and ja tabulars modified.
*/
void z_csc_noDiag(pastix_int_t        baseval,
                  pastix_int_t        n,
                  pastix_int_t       *ia,
                  pastix_int_t       *ja,
                  pastix_complex64_t *a)
{
  pastix_int_t i, j;
  pastix_int_t indj;
  pastix_int_t *old_ia = NULL;

  MALLOC_INTERN(old_ia, n+1, pastix_int_t);
  memcpy(old_ia, ia, sizeof(pastix_int_t)*(n+1));

  ASSERT(ia[0]==baseval,MOD_SOPALIN);

  indj = 0;
  /*fprintf(stdout, "NNZ with diag = %ld \n", ia[n]);*/

  for(i=0;i<n;i++)
    {
      /* ia[i] = number of column already counted */
      ia[i] = indj+baseval;
      /* for each row number in each column i */
      for(j=old_ia[i];j<old_ia[i+1];j++)
        /* if element is not diagonal
           we add it in ja and we count it */
        if(ja[j-baseval] != i+baseval)
          {
            ja[indj] = ja[j-baseval];
            if (a != NULL)
              a[indj] = a [j -baseval];
            indj++;
          }
    }
  ia[n] = indj+baseval;

  assert( ia[n] <= old_ia[n] );

  /*fprintf(stdout, "NNZ without diag = %ld \n", ia[n]);*/
  memFree_null(old_ia);
}

/*
  Function: z_csc_check_doubles

  Check if the csc contains doubles and if correct if asked

  Assumes that the CSC is sorted.

  Assumes that the CSC is Fortran numeroted (base 1)

  Parameters:
    n         - Size of the matrix.
    colptr    - Index in *rows* and *values* of the first element of each column
    rows      - row of each element
    values    - value of each element (Can be NULL)
    dof       - Number of degrees of freedom
    flag      - Indicate if user wants correction (<API_BOOLEAN>)
    flagalloc - indicate if allocation on CSC uses internal malloc.

  Returns:
    API_YES - If the matrix contained no double or was successfully corrected.
    API_NO  - Otherwise.
*/
int z_csc_check_doubles(pastix_int_t         n,
                        pastix_int_t        *colptr,
                        pastix_int_t       **rowptr,
                        pastix_complex64_t **values,
                        int                  dof,
                        int                  flag,
                        int                  flagalloc)
{
  pastix_int_t     i,j,k,d;
  int     doubles = 0;
  pastix_int_t   * tmprows = NULL;
  pastix_complex64_t * tmpvals = NULL;
  pastix_int_t     index = 0;
  pastix_int_t     lastindex = 0;

  ASSERT(values == NULL || dof > 0, MOD_SOPALIN);
  ASSERT(flag == API_NO || flag == API_YES, MOD_SOPALIN);
  ASSERT(colptr[0] == 1, MOD_SOPALIN);
  ASSERT(n >= 0, MOD_SOPALIN);

  for (i = 0; i < n; i++)
    {
      for (j = colptr[i]-1; j < colptr[i+1]-1; j = k)
        {
          (*rowptr)[index]   = (*rowptr)[j];
          if (values != NULL)
            for (d = 0; d < dof*dof; d++)
              (*values)[index*dof*dof+d] = (*values)[j*dof*dof+d];

          k = j+1;
          while (k < colptr[i+1]-1 && (*rowptr)[j] == (*rowptr)[k])
            {
              if (flag == API_NO)
                return API_NO;
              if (values != NULL)
                for (d = 0; d < dof*dof; d++)
                  (*values)[index*dof*dof+d] += (*values)[k*dof*dof+d];
              doubles++;
              k++;
            }
          index++;
        }

      colptr[i] = lastindex+1;
      lastindex = index;
    }
  if (flag == API_NO)
    return API_YES;
  ASSERT(index == colptr[n]-1-doubles, MOD_SOPALIN);
  if (doubles > 0)
    {
      colptr[n] = lastindex+1;
      if (flagalloc == API_NO)
        {
          MALLOC_EXTERN(tmprows, lastindex, pastix_int_t);
          if (values != NULL)
            MALLOC_EXTERN(tmpvals, lastindex*dof*dof, pastix_complex64_t);
        }
      else
        {
          MALLOC_INTERN(tmprows, lastindex, pastix_int_t);
          if (values != NULL)
            MALLOC_INTERN(tmpvals, lastindex*dof*dof, pastix_complex64_t);
        }

      memcpy(tmprows, *rowptr,   lastindex*sizeof(pastix_int_t));
      if (values != NULL)
        memcpy(tmpvals, *values, lastindex*dof*dof*sizeof(pastix_complex64_t));
      if (flagalloc == API_NO)
        {
          free(*rowptr);
          if (values != NULL)
            free(*values);
        }
      else
        {
          memFree_null(*rowptr);
          if (values != NULL)
            memFree_null(*values);
        }
      *rowptr   = tmprows;
      if (values != NULL)
        *values = tmpvals;
    }
  return API_YES;

}

/*
  Function: z_csc_checksym

    Check if the CSC graph is symetric.

    For all local column C,

    For all row R in the column C,

    We look in column R if we have the row number C.

    If we can correct we had missing non zeros.

    Assumes that the CSC is Fortran numbered (1 based).

    Assumes that the matrix is sorted.

  Parameters:
    n        - Number of local columns
    colptr   - Starting index of each columns in *ja*
    rows     - Row of each element.
    values   - Value of each element.
    correct  - Flag indicating if we can correct the symmetry.
    alloc    - indicate if allocation on CSC uses internal malloc.
    dof      - Number of degrees of freedom.
*/
int z_csc_checksym(pastix_int_t         n,
                   pastix_int_t        *colptr,
                   pastix_int_t       **rows,
                   pastix_complex64_t **values,
                   int                  correct,
                   int                  alloc,
                   int                  dof)
{
  pastix_int_t            i,j,k,l,d;
  pastix_int_t            index1;
  pastix_int_t            index2;
  int            found;
  pastix_int_t            toaddsize;
  pastix_int_t            toaddcnt;
  pastix_int_t         *  toadd      = NULL;
  pastix_int_t         *  tmpcolptr  = NULL;
  pastix_int_t         *  tmprows    = NULL;
  pastix_complex64_t       *  tmpvals    = NULL;

  /* For all local column C,
     For all row R in the column C,

     If the row number R correspond to a local column,
     We look in column R if we have the row number C.

     Else,
   */
  toaddcnt  = 0;
  toaddsize = 0;
  for (i = 0; i < n; i++)
    {
      for (j = (colptr)[i]-1; j < (colptr)[i+1]-1; j++)
        {
          if ((*rows)[j] != i+1)
            {
              /* not in diagonal */
              k = (*rows)[j];
              found = 0;
              for (l = (colptr)[k-1]-1; l < (colptr)[k-1+1]-1; l++)
                {
                  if (i+1 == (*rows)[l])
                    {
                      found = 1;
                      break;
                    }
                  if (i+1 < (*rows)[l])
                    {
                      /* The CSC is sorted */
                      found = 0;
                      break;
                    }
                }
              if (found == 0)
                {
                  if (correct == API_NO)
                    return EXIT_FAILURE;
                  else
                    {
                      if (toaddsize == 0)
                        {
                          toaddsize = n/2;
                          MALLOC_INTERN(toadd, 2*toaddsize, pastix_int_t);
                        }
                      if (toaddcnt >= toaddsize)
                        {
                          toaddsize += toaddsize/2 + 1;
                          if (NULL ==
                              (toadd =
                               (pastix_int_t*)memRealloc(toadd,
                                                2*toaddsize*sizeof(pastix_int_t))))
                              MALLOC_ERROR("toadd");
                        }
                      toadd[2*toaddcnt]     = (*rows)[j];
                      toadd[2*toaddcnt + 1] = i+1;
                      /* fprintf(stdout, "Adding %ld, %ld\n", (long)(i+1), (long)(*rows)[j]); */
                      toaddcnt++;
                    }
                }
            }
        }
    }

  if (toaddcnt > 0)
    {

      intSort2asc1(toadd, toaddcnt);
      /* Correct here is API_YES, otherwise we would have return EXIT_FAILURE
         Or toaddcnt == 0*/
      MALLOC_INTERN(tmpcolptr, n + 1, pastix_int_t);
      if (alloc == API_NO)
        {
          MALLOC_EXTERN(tmprows, colptr[n]-1 + toaddcnt, pastix_int_t);
          if (values != NULL)
            {
              MALLOC_EXTERN(tmpvals, colptr[n]-1 + toaddcnt, pastix_complex64_t);
            }
        }
      else
        {
          MALLOC_INTERN(tmprows, colptr[n]-1 + toaddcnt, pastix_int_t);
          if (values != NULL)
            {
              MALLOC_INTERN(tmpvals, colptr[n]-1 + toaddcnt, pastix_complex64_t);
            }
        }
      /* Build tmpcolptr

         tmpcolptr[i+1] will contain the number of element of
         the column i
       */
      index1 = 0;
      index2 = 0;
      for (i = 0; i <  n; i++)
        {
          tmpcolptr[i] = index2+1;
          for (j = colptr[i]-1; j < colptr[i+1]-1; j++)
            {
              if (index1 < toaddcnt &&
                  (toadd[2*index1] == i+1) &&
                  (toadd[2*index1+1] < (*rows)[j]))
                {
                  tmprows[index2] = toadd[2*index1+1];
                  if (values != NULL)
                    {
                      for (d = 0; d < dof*dof ; d++)
                        tmpvals[index2*dof*dof+d] = 0.0;
                    }
                  index1++;
                  j--; /* hack do not increment j this step of the loop */
                }
              else
                {
                  tmprows[index2] = (*rows)[j];
                  if (values != NULL)
                    {
                      for (d = 0; d < dof*dof ; d++)
                        tmpvals[index2*dof*dof+d] = (*values)[j*dof*dof+d];
                    }
                }
              index2++;
            }

          while(index1 < toaddcnt && toadd[2*index1] == i+1)
            {
              tmprows[index2] = toadd[2*index1+1];
              if (values != NULL)
                {
                  for (d = 0; d < dof*dof ; d++)
                    tmpvals[index2*dof*dof+d] = 0.0;
                }
              index1++;
              index2++;
            }
        }
      tmpcolptr[n] = index2+1;
      ASSERT((tmpcolptr[n] - 1) == (colptr[n] - 1 + toaddcnt), MOD_SOPALIN);

      memcpy(colptr, tmpcolptr, (n+1)*sizeof(pastix_int_t));
      memFree_null(tmpcolptr);
      memFree_null(toadd);
      if (alloc == API_NO)
        {
          free(*rows);
          if (values != NULL)
            {
              free(*values);
            }
        }
      else
        {
          memFree_null(*rows);
          if (values != NULL)
            {
              memFree_null(*values);
            }
        }
      *rows   = tmprows;
      if (values != NULL)
        {
          *values = tmpvals;
        }
    }
  return EXIT_SUCCESS;
}


/*
  Function: z_csc_colPerm

  Performs column permutation on a CSC

  Parameters:
    n     - Size of the matrix.
    ia    - Index of first element of each column in *ia* and *a*
    ja    - Rows of non zeros of the matrix.
    a     - Values of non zeros of the matrix.
    cperm - Permutation to perform
*/
void z_csc_colPerm(pastix_int_t n, pastix_int_t *ia, pastix_int_t *ja, pastix_complex64_t *a, pastix_int_t *cperm)
{
  pastix_int_t i, k;
  pastix_int_t   *newja = NULL;
  pastix_int_t   *newia = NULL;
  pastix_complex64_t *newa  = NULL;
  int numflag, numflag2;

  numflag = ia[0];
  numflag2 = 1;
  for(i=0;i<n;i++)
    if(cperm[i] == 0)
      {
        numflag2 = 0;
        break;
      }

  if(numflag2 != numflag)
    {
      errorPrint("CSC_colPerm: rperm not in same numbering than the CSC.");
      exit(-1);
    }


  if(numflag == 1)
    {
      z_csc_Fnum2Cnum(ja, ia, n);
      for(i=0;i<n;i++)
        cperm[i]--;
    }

  MALLOC_INTERN(newia, n+1,   pastix_int_t);
  MALLOC_INTERN(newja, ia[n], pastix_int_t);
  MALLOC_INTERN(newa,  ia[n], pastix_complex64_t);


  newia[0] = 0;
  for(i=0;i<n;i++)
    {
#ifdef DEBUG_KASS
      ASSERT(cperm[i]>=0 && cperm[i] < n, MOD_KASS);
#endif
      newia[cperm[i]+1] = ia[i+1]-ia[i];
    }

#ifdef DEBUG_KASS
  for(i=1;i<=n;i++)
    ASSERT(newia[i] >0, MOD_KASS);
#endif

  for(i=1;i<=n;i++)
    newia[i] += newia[i-1];

#ifdef DEBUG_KASS
  ASSERT(newia[n] == ia[n], MOD_KASS);
#endif


  for(i=0;i<n;i++)
    {
      k = cperm[i];
#ifdef DEBUG_KASS
      ASSERT(newia[k+1]-newia[k] == ia[i+1]-ia[i], MOD_KASS);
#endif
      memcpy(newja + newia[k], ja + ia[i], sizeof(pastix_int_t)*(ia[i+1]-ia[i]));
      memcpy(newa + newia[k], a + ia[i], sizeof(pastix_complex64_t)*(ia[i+1]-ia[i]));

    }

  memcpy(ia, newia, sizeof(pastix_int_t)*(n+1));
  memcpy(ja, newja, sizeof(pastix_int_t)*ia[n]);
  memcpy(a, newa, sizeof(pastix_complex64_t)*ia[n]);

  memFree(newia);
  memFree(newja);
  memFree(newa);
  if(numflag == 1)
    {
      z_csc_Cnum2Fnum(ja, ia, n);
      for(i=0;i<n;i++)
        cperm[i]++;
    }
}


/*
  Function: z_csc_colScale

  Moved from kass, only used in MC64
*/
void z_csc_colScale(pastix_int_t        n,
                    pastix_int_t       *ia,
                    pastix_int_t       *ja,
                    pastix_complex64_t *a,
                    pastix_complex64_t *dcol)
{
  pastix_int_t i, j;
  int numflag;
  pastix_complex64_t d;
  numflag = ia[0];

  if(numflag == 1)
    z_csc_Fnum2Cnum(ja, ia, n);

  for(i=0;i<n;i++)
    {
      d = dcol[i];
      for(j=ia[i];j<ia[i+1];j++)
        {
          /***@@@ OIMBE DSCAL **/
          a[j] *= d;
        }
    }

  if(numflag == 1)
    z_csc_Cnum2Fnum(ja, ia, n);
}

/*
  Function: z_csc_rowScale

  Moved from kass, only used in MC64
*/
void z_csc_rowScale(pastix_int_t        n,
                    pastix_int_t       *ia,
                    pastix_int_t       *ja,
                    pastix_complex64_t *a,
                    pastix_complex64_t *drow)
{
  pastix_int_t i, j;
  int numflag;
  numflag = ia[0];

  if(numflag == 1)
    z_csc_Fnum2Cnum(ja, ia, n);

  for(i=0;i<n;i++)
    {
      for(j=ia[i];j<ia[i+1];j++)
        {
#ifdef DEBUG_KASS
          ASSERT(ja[j]>0 && ja[j] <n, MOD_KASS);
#endif
          a[j] *= drow[ja[j]];
        }
    }

  if(numflag == 1)
    z_csc_Cnum2Fnum(ja, ia, n);
}

/*
 * z_csc_sort:
 *
 * Sort CSC columns
 *
 * Parameters:
 *   n  - Number of columns
 *   ia - Index of first element of each column in *ia*.
 *   ja - Rows of each non zeros.
 *   a  - Values of each non zeros.
*/
void z_csc_sort(pastix_int_t        n,
                pastix_int_t       *ia,
                pastix_int_t       *ja,
                pastix_complex64_t *a,
                pastix_int_t        ndof)
{
  pastix_int_t i;
  int numflag;
  pastix_int_t ndof2 = ndof * ndof;
  void * sortptr[3];
  numflag = ia[0];
  if(numflag == 1)
    z_csc_Fnum2Cnum(ja, ia, n);
  if (a != NULL)
    {

      for(i=0;i<n;i++)
        {
          sortptr[0] = &ja[ia[i]];
          sortptr[1] = &a[ia[i]*ndof2];
          sortptr[2] = &ndof2;
          z_qsortIntFloatAsc(sortptr, ia[i+1] - ia[i]);
        }

    }
  else
    {
      for(i=0;i<n;i++)
        intSort1asc1(&ja[ia[i]], ia[i+1] - ia[i]);

    }
  if(numflag == 1)
    z_csc_Cnum2Fnum(ja, ia, n);
}

/*
  Function: z_csc_Fnum2Cnum

  Convert CSC numbering from fortran numbering to C numbering.

  Parameters:
    ja - Rows of each element.
    ia - First index of each column in *ja*
    n  - Number of columns
*/
void z_csc_Fnum2Cnum(pastix_int_t *ja,
                     pastix_int_t *ia,
                     pastix_int_t  n)
{
  pastix_int_t i, j;
  for(i=0;i<=n;i++)
    ia[i]--;

  for(i=0;i<n;i++)
    for(j=ia[i];j<ia[i+1];j++)
      ja[j]--;

}

/*
  Function: z_csc_Cnum2Fnum

  Convert CSC numbering from C numbering to Fortran numbering.

  Parameters:
    ja - Rows of each element.
    ia - First index of each column in *ja*
    n  - Number of columns
*/
void z_csc_Cnum2Fnum(pastix_int_t *ja,
                     pastix_int_t *ia,
                     pastix_int_t  n)
{
  pastix_int_t i, j;

  for(i=0;i<n;i++)
    for(j=ia[i];j<ia[i+1];j++)
      ja[j]++;

  for(i=0;i<=n;i++)
    ia[i]++;
}
/*
  Function: z_csc_buildZerosAndNonZerosGraphs

  Separate a graph in two graphs, following
  wether the diagonal term of a column is null or not.

  Parameters:
    n, colptr, rows, values  - The initial CSC
    n_nz, colptr_nz, rows_nz - The graph of the non-null diagonal part.
    n_z, colptr_z, rows_z    - The graph of the null diagonal part.
    perm                     - Permutation to go from the first graph to
                               the one composed of the two graph concatenated.
    revperm                  - Reverse permutation tabular.
    criteria                 - Value beside which a number is said null.
*/
int z_csc_buildZerosAndNonZerosGraphs(pastix_int_t         n,
                                      pastix_int_t        *colptr,
                                      pastix_int_t        *rows,
                                      pastix_complex64_t  *values,
                                      pastix_int_t        *n_nz,
                                      pastix_int_t       **colptr_nz,
                                      pastix_int_t       **rows_nz,
                                      pastix_int_t        *n_z,
                                      pastix_int_t       **colptr_z,
                                      pastix_int_t       **rows_z,
                                      pastix_int_t        *perm,
                                      pastix_int_t        *revperm,
                                      double               criteria)
{
  pastix_int_t  itercol;
  pastix_int_t  iterrow;

  pastix_int_t  ncoefszeros  = 0;
  pastix_int_t  ncoefsnzeros = 0;
  pastix_int_t  itercol_nz   = 0;
  pastix_int_t  itercol_z    = 0;
  int  seen;
  pastix_int_t  cntrows;

  for (itercol = 0; itercol <n; itercol++)
    {
      seen = 0;
      for (iterrow = colptr[itercol]-1; iterrow < colptr[itercol+1]-1; iterrow++)
        {
          if (itercol == rows[iterrow] -1 )
            {
              if (ABS_FLOAT(values[iterrow]) < criteria)
                {
                  (*n_z) ++;
                  ncoefszeros += colptr[itercol+1] - colptr[itercol];
                  seen = 1;
                }
              else
                {
                  (*n_nz) ++;
                  ncoefsnzeros += colptr[itercol+1] - colptr[itercol];
                  seen = 1;
                }
              break;
            }
        }
      if (colptr[itercol+1] == colptr[itercol]) /* empty column */
        {
          (*n_z)++;