Actual source code: ij.c


  2: #include <../src/mat/impls/aij/seq/aij.h>

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

  9:     Description:
 10:     Take the data in the row-oriented sparse storage and build the
 11:     IJ data for the Matrix.  Return 0 on success,row + 1 on failure
 12:     at that row. Produces the ij for a symmetric matrix by using
 13:     the lower triangular part of the matrix if lower_triangular is PETSC_TRUE;
 14:     it uses the upper triangular otherwise.

 16:     Input Parameters:
 17: .   Matrix - matrix to convert
 18: .   lower_triangular - symmetrize the lower triangular part
 19: .   shiftin - the shift for the original matrix (0 or 1)
 20: .   shiftout - the shift required for the ordering routine (0 or 1)

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

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

 36:   PetscFunctionBegin;
 37:   /* allocate space for row pointers */
 38:   PetscCall(PetscCalloc1(m + 1, &ia));
 39:   *iia = ia;
 40:   PetscCall(PetscMalloc1(m + 1, &work));

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

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

 66:   /* allocate space for column pointers */
 67:   nz = ia[m] + (!shiftin);
 68:   PetscCall(PetscMalloc1(nz, &ja));
 69:   *jja = ja;

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