Actual source code: ij.c

      1: #define PETSCMAT_DLL

     3:  #include src/mat/impls/aij/seq/aij.h

      7: /*
      8:   MatToSymmetricIJ_SeqAIJ - Convert a (generally nonsymmetric) sparse AIJ matrix
      9:            to IJ format (ignore the "A" part) Allocates the space needed. Uses only
     10:            the lower triangular part of the matrix.

     12:     Description:
     13:     Take the data in the row-oriented sparse storage and build the
     14:     IJ data for the Matrix.  Return 0 on success,row + 1 on failure
     15:     at that row. Produces the ij for a symmetric matrix by only using
     16:     the lower triangular part of the matrix.

     18:     Input Parameters:
     19: .   Matrix - matrix to convert
     20: .   shiftin - the shift for the original matrix (0 or 1)
     21: .   shiftout - the shift required for the ordering routine (0 or 1)

     23:     Output Parameters:
     24: .   ia     - ia part of IJ representation (row information)
     25: .   ja     - ja part (column indices)

     27:     Notes:
     28:     Both ia and ja may be freed with PetscFree();
     29:     This routine is provided for ordering routines that require a
     30:     symmetric structure.  It is required since those routines call
     31:     SparsePak routines that expect a symmetric  matrix.
     32: */
     33: PetscErrorCode MatToSymmetricIJ_SeqAIJ(PetscInt m,PetscInt *ai,PetscInt *aj,PetscInt shiftin,PetscInt shiftout,PetscInt **iia,PetscInt **jja)
     34: {
     36:   PetscInt       *work,*ia,*ja,*j,i,nz,row,col;

     39:   /* allocate space for row pointers */
     40:   PetscMalloc((m+1)*sizeof(PetscInt),&ia);
     41:   *iia = ia;
     42:   PetscMemzero(ia,(m+1)*sizeof(PetscInt));
     43:   PetscMalloc((m+1)*sizeof(PetscInt),&work);

     45:   /* determine the number of columns in each row */
     46:   ia[0] = shiftout;
     47:   for (row = 0; row < m; row++) {
     48:     nz = ai[row+1] - ai[row];
     49:     j  = aj + ai[row] + shiftin;
     50:     while (nz--) {
     51:        col = *j++ + shiftin;
     52:        if (col > row) { break;}
     53:        if (col != row) ia[row+1]++;
     54:        ia[col+1]++;
     55:     }
     56:   }

     58:   /* shiftin ia[i] to point to next row */
     59:   for (i=1; i<m+1; i++) {
     60:     row       = ia[i-1];
     61:     ia[i]     += row;
     62:     work[i-1] = row - shiftout;
     63:   }

     65:   /* allocate space for column pointers */
     66:   nz   = ia[m] + (!shiftin);
     67:   PetscMalloc(nz*sizeof(PetscInt),&ja);
     68:   *jja = ja;

     70:   /* loop over lower triangular part putting into ja */
     71:   for (row = 0; row < m; row++) {
     72:     nz = ai[row+1] - ai[row];
     73:     j  = aj + ai[row] + shiftin;
     74:     while (nz--) {
     75:       col = *j++ + shiftin;
     76:       if (col > row) { break;}
     77:       if (col != row) {ja[work[col]++] = row + shiftout; }
     78:       ja[work[row]++] = col + shiftout;
     79:     }
     80:   }
     81:   PetscFree(work);
     82:   return(0);
     83: }