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

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

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

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

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

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

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

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

 63:   /* shiftin ia[i] to point to next row */
 64:   for (i=1; i<m+1; i++) {
 65:     row       = ia[i-1];
 66:     ia[i]    += row;
 67:     work[i-1] = row - shiftout;
 68:   }

 70:   /* allocate space for column pointers */
 71:   nz   = ia[m] + (!shiftin);
 72:   PetscMalloc1(nz,&ja);
 73:   *jja = ja;

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