Actual source code: sro.c

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

      4: /*
      5: This function is used before applying a
      6: symmetric reordering to matrix A that is
      7: in SBAIJ format.

      9: The permutation is assumed to be symmetric, i.e.,
     10: P = P^T (= inv(P)),
     11: so the permuted matrix P*A*inv(P)=P*A*P^T is ensured to be symmetric.
     12:  - a wrong assumption! This code needs rework!  -- Hong

     14: The function is modified from sro.f of YSMP. The description from YSMP:
     15: C    THE NONZERO ENTRIES OF THE MATRIX M ARE ASSUMED TO BE STORED
     16: C    SYMMETRICALLY IN (IA,JA,A) FORMAT (I.E., NOT BOTH M(I,J) AND M(J,I)
     17: C    ARE STORED IF I NE J).
     18: C
     19: C    SRO DOES NOT REARRANGE THE ORDER OF THE ROWS, BUT DOES MOVE
     20: C    NONZEROES FROM ONE ROW TO ANOTHER TO ENSURE THAT IF M(I,J) WILL BE
     21: C    IN THE UPPER TRIANGLE OF M WITH RESPECT TO THE NEW ORDERING, THEN
     22: C    M(I,J) IS STORED IN ROW I (AND THUS M(J,I) IS NOT STORED);  WHEREAS
     23: C    IF M(I,J) WILL BE IN THE STRICT LOWER TRIANGLE OF M, THEN M(J,I) IS
     24: C    STORED IN ROW J (AND THUS M(I,J) IS NOT STORED).

     26:   -- output: new index set (inew, jnew) for A and a map a2anew that maps
     27:              values a to anew, such that all
     28:              nonzero A_(perm(i),iperm(k)) will be stored in the upper triangle.
     29:              Note: matrix A is not permuted by this function!
     30: */
     31: PetscErrorCode MatReorderingSeqSBAIJ(Mat A, IS perm)
     32: {
     33:   Mat_SeqSBAIJ  *a   = (Mat_SeqSBAIJ *)A->data;
     34:   const PetscInt mbs = a->mbs;

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

     46:   PetscCall(ISInvertPermutation(perm,PETSC_DECIDE,&iperm));
     47:   PetscCall(ISGetIndices(iperm,&riip));

     49:   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");
     50:   PetscCall(ISRestoreIndices(iperm,&riip));
     51:   PetscCall(ISDestroy(&iperm));

     53:   if (!a->inew) PetscCall(PetscMalloc2(mbs+1,&ai, 2*a->i[mbs],&aj));
     54:   else {
     55:     ai = a->inew;
     56:     aj = a->jnew;
     57:   }
     58:   PetscCall(PetscArraycpy(ai,a->i,mbs+1));
     59:   PetscCall(PetscArraycpy(aj,a->j,a->i[mbs]));

     61:   /*
     62:      Phase 1: Find row index r in which to store each nonzero.
     63:               Initialize count of nonzeros to be stored in each row (nzr).
     64:               At the end of this phase, a nonzero a(*,*)=a(r(),aj())
     65:               s.t. a(perm(r),perm(aj)) will fall into upper triangle part.
     66:   */

     68:   PetscCall(PetscMalloc1(mbs,&nzr));
     69:   PetscCall(PetscMalloc1(ai[mbs],&r));
     70:   for (i=0; i<mbs; i++) nzr[i] = 0;
     71:   for (i=0; i<ai[mbs]; i++) r[i] = 0;

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

     86:       r[j] = k; j++;
     87:       nzr[k]++;  /* increment count of nonzeros in that row */
     88:     }
     89:   }

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

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

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

    119:   /* Phase 3: permute (aj,a) to upper triangular form (wrt new ordering) */
    120:   for (j=jmin; j<jmax; j++) {
    121:     while (r[j] != j) {
    122:       k   = r[j]; r[j] = r[k]; r[k] = k;
    123:       ajk = aj[k]; aj[k] = aj[j]; aj[j] = ajk;
    124:       /* ak = aa[k]; aa[k] = aa[j]; aa[j] = ak; */
    125:     }
    126:   }
    127:   PetscCall(ISRestoreIndices(perm,&rip));

    129:   a->inew = ai;
    130:   a->jnew = aj;

    132:   PetscCall(ISDestroy(&a->row));
    133:   PetscCall(ISDestroy(&a->icol));
    134:   PetscCall(PetscObjectReference((PetscObject)perm));
    135:   PetscCall(ISDestroy(&a->row));
    136:   a->row  = perm;
    137:   PetscCall(PetscObjectReference((PetscObject)perm));
    138:   PetscCall(ISDestroy(&a->icol));
    139:   a->icol = perm;

    141:   PetscCall(PetscFree(nzr));
    142:   PetscCall(PetscFree(r));
    143:   PetscFunctionReturn(PETSC_SUCCESS);
    144: #endif
    145: }