1: #include <../src/mat/impls/baij/seq/baij.h>
  2: #include <../src/mat/impls/sbaij/seq/sbaij.h>

  4: /*@
  5:   MatReorderingSeqSBAIJ - Prepare an updated index structure for a symmetric reordering of a `MATSEQSBAIJ` matrix.

  7:   Not Collective

  9:   Input Parameters:
 10: + A    - the `MATSEQSBAIJ` matrix
 11: - perm - the (assumed symmetric) permutation to be applied

 13:   Level: developer

 15:   Notes:
 16:   This routine currently raises `PETSC_ERR_SUP`; matrix reordering is not supported for `MATSEQSBAIJ` matrices,
 17:   so callers should convert to `MATSEQAIJ` first.

 19:   The intent (unimplemented) is to compute a new index set `(inew, jnew)` for `A` and a value map so that
 20:   all nonzero entries `A(perm(i), perm(k))` are stored in the upper triangle; the matrix itself is not permuted.

 22: .seealso: `Mat`, `MATSEQSBAIJ`, `MatGetOrdering()`, `MatPermute()`
 23: @*/
 24: PetscErrorCode MatReorderingSeqSBAIJ(Mat A, IS perm)
 25: {
 26:   Mat_SeqSBAIJ  *a   = (Mat_SeqSBAIJ *)A->data;
 27:   const PetscInt mbs = a->mbs;

 29:   PetscFunctionBegin;
 30:   if (!mbs) PetscFunctionReturn(PETSC_SUCCESS);
 31:   SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Matrix reordering is not supported for sbaij matrix. Use aij format");
 32: #if 0
 33:   const PetscInt *rip,*riip;
 34:   PetscInt       *ai,*aj,*r;
 35:   PetscInt       *nzr,nz,jmin,jmax,j,k,ajk,i;
 36:   IS             iperm;  /* inverse of perm */
 37:   PetscCall(ISGetIndices(perm,&rip));

 39:   PetscCall(ISInvertPermutation(perm,PETSC_DECIDE,&iperm));
 40:   PetscCall(ISGetIndices(iperm,&riip));

 42:   for (i=0; i<mbs; i++) PetscCheck(rip[i] == riip[i],PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Non-symmetric permutation, use symmetric permutation for symmetric matrices");
 43:   PetscCall(ISRestoreIndices(iperm,&riip));
 44:   PetscCall(ISDestroy(&iperm));

 46:   if (!a->inew) PetscCall(PetscMalloc2(mbs+1,&ai, 2*a->i[mbs],&aj));
 47:   else {
 48:     ai = a->inew;
 49:     aj = a->jnew;
 50:   }
 51:   PetscCall(PetscArraycpy(ai,a->i,mbs+1));
 52:   PetscCall(PetscArraycpy(aj,a->j,a->i[mbs]));

 54:   /*
 55:      Phase 1: Find row index r in which to store each nonzero.
 56:               Initialize count of nonzeros to be stored in each row (nzr).
 57:               At the end of this phase, a nonzero a(*,*)=a(r(),aj())
 58:               s.t. a(perm(r),perm(aj)) will fall into upper triangle part.
 59:   */

 61:   PetscCall(PetscMalloc1(mbs,&nzr));
 62:   PetscCall(PetscMalloc1(ai[mbs],&r));
 63:   for (i=0; i<mbs; i++) nzr[i] = 0;
 64:   for (i=0; i<ai[mbs]; i++) r[i] = 0;

 66:   /*  for each nonzero element */
 67:   for (i=0; i<mbs; i++) {
 68:     nz = ai[i+1] - ai[i];
 69:     j  = ai[i];
 70:     /* printf("nz = %d, j=%d\n",nz,j); */
 71:     while (nz--) {
 72:       /*  --- find row (=r[j]) and column (=aj[j]) in which to store a[j] ...*/
 73:       k = aj[j];                          /* col. index */
 74:       /* printf("nz = %d, k=%d\n", nz,k); */
 75:       /* for entry that will be permuted into lower triangle, swap row and col. index */
 76:       if (rip[k] < rip[i]) aj[j] = i;
 77:       else k = i;

 79:       r[j] = k; j++;
 80:       nzr[k]++;  /* increment count of nonzeros in that row */
 81:     }
 82:   }

 84:   /* Phase 2: Find new ai and permutation to apply to (aj,a).
 85:               Determine pointers (r) to delimit rows in permuted (aj,a).
 86:               Note: r is different from r used in phase 1.
 87:               At the end of this phase, (aj[j],a[j]) will be stored in
 88:               (aj[r(j)],a[r(j)]).
 89:   */
 90:   for (i=0; i<mbs; i++) {
 91:     ai[i+1] = ai[i] + nzr[i];
 92:     nzr[i]  = ai[i+1];
 93:   }

 95:   /* determine where each (aj[j], a[j]) is stored in new (aj,a)
 96:      for each nonzero element (in reverse order) */
 97:   jmin = ai[0]; jmax = ai[mbs];
 98:   nz   = jmax - jmin;
 99:   j    = jmax-1;
100:   while (nz--) {
101:     i = r[j];  /* row value */
102:     if (aj[j] == i) r[j] = ai[i]; /* put diagonal nonzero at beginning of row */
103:     else { /* put off-diagonal nonzero in last unused location in row */
104:       nzr[i]--; r[j] = nzr[i];
105:     }
106:     j--;
107:   }

109:   a->a2anew = aj + ai[mbs];
110:   PetscCall(PetscArraycpy(a->a2anew,r,ai[mbs]));

112:   /* Phase 3: permute (aj,a) to upper triangular form (wrt new ordering) */
113:   for (j=jmin; j<jmax; j++) {
114:     while (r[j] != j) {
115:       k   = r[j]; r[j] = r[k]; r[k] = k;
116:       ajk = aj[k]; aj[k] = aj[j]; aj[j] = ajk;
117:       /* ak = aa[k]; aa[k] = aa[j]; aa[j] = ak; */
118:     }
119:   }
120:   PetscCall(ISRestoreIndices(perm,&rip));

122:   a->inew = ai;
123:   a->jnew = aj;

125:   PetscCall(ISDestroy(&a->row));
126:   PetscCall(ISDestroy(&a->icol));
127:   PetscCall(PetscObjectReference((PetscObject)perm));
128:   PetscCall(ISDestroy(&a->row));
129:   a->row  = perm;
130:   PetscCall(PetscObjectReference((PetscObject)perm));
131:   PetscCall(ISDestroy(&a->icol));
132:   a->icol = perm;

134:   PetscCall(PetscFree(nzr));
135:   PetscCall(PetscFree(r));
136:   PetscFunctionReturn(PETSC_SUCCESS);
137: #endif
138: }