Actual source code: sro.c
2: #include <../src/mat/impls/baij/seq/baij.h>
3: #include <../src/mat/impls/sbaij/seq/sbaij.h>
5: /*
6: This function is used before applying a
7: symmetric reordering to matrix A that is
8: in SBAIJ format.
10: The permutation is assumed to be symmetric, i.e.,
11: P = P^T (= inv(P)),
12: so the permuted matrix P*A*inv(P)=P*A*P^T is ensured to be symmetric.
13: - a wrong assumption! This code needs rework! -- Hong
15: The function is modified from sro.f of YSMP. The description from YSMP:
16: C THE NONZERO ENTRIES OF THE MATRIX M ARE ASSUMED TO BE STORED
17: C SYMMETRICALLY IN (IA,JA,A) FORMAT (I.E., NOT BOTH M(I,J) AND M(J,I)
18: C ARE STORED IF I NE J).
19: C
20: C SRO DOES NOT REARRANGE THE ORDER OF THE ROWS, BUT DOES MOVE
21: C NONZEROES FROM ONE ROW TO ANOTHER TO ENSURE THAT IF M(I,J) WILL BE
22: C IN THE UPPER TRIANGLE OF M WITH RESPECT TO THE NEW ORDERING, THEN
23: C M(I,J) IS STORED IN ROW I (AND THUS M(J,I) IS NOT STORED); WHEREAS
24: C IF M(I,J) WILL BE IN THE STRICT LOWER TRIANGLE OF M, THEN M(J,I) IS
25: C STORED IN ROW J (AND THUS M(I,J) IS NOT STORED).
27: -- output: new index set (inew, jnew) for A and a map a2anew that maps
28: values a to anew, such that all
29: nonzero A_(perm(i),iperm(k)) will be stored in the upper triangle.
30: Note: matrix A is not permuted by this function!
31: */
32: PetscErrorCode MatReorderingSeqSBAIJ(Mat A,IS perm)
33: {
34: Mat_SeqSBAIJ *a=(Mat_SeqSBAIJ*)A->data;
35: const PetscInt mbs=a->mbs;
38: if (!mbs) return(0);
39: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Matrix reordering is not supported for sbaij matrix. Use aij format");
40: #if 0
42: const PetscInt *rip,*riip;
43: PetscInt *ai,*aj,*r;
44: PetscInt *nzr,nz,jmin,jmax,j,k,ajk,i;
45: IS iperm; /* inverse of perm */
46: ISGetIndices(perm,&rip);
48: ISInvertPermutation(perm,PETSC_DECIDE,&iperm);
49: ISGetIndices(iperm,&riip);
51: for (i=0; i<mbs; i++) {
52: if (rip[i] != riip[i]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Non-symmetric permutation, use symmetric permutation for symmetric matrices");
53: }
54: ISRestoreIndices(iperm,&riip);
55: ISDestroy(&iperm);
57: if (!a->inew) {
58: PetscMalloc2(mbs+1,&ai, 2*a->i[mbs],&aj);
59: } else {
60: ai = a->inew; aj = a->jnew;
61: }
62: PetscArraycpy(ai,a->i,mbs+1);
63: PetscArraycpy(aj,a->j,a->i[mbs]);
65: /*
66: Phase 1: Find row index r in which to store each nonzero.
67: Initialize count of nonzeros to be stored in each row (nzr).
68: At the end of this phase, a nonzero a(*,*)=a(r(),aj())
69: s.t. a(perm(r),perm(aj)) will fall into upper triangle part.
70: */
72: PetscMalloc1(mbs,&nzr);
73: PetscMalloc1(ai[mbs],&r);
74: for (i=0; i<mbs; i++) nzr[i] = 0;
75: for (i=0; i<ai[mbs]; i++) r[i] = 0;
77: /* for each nonzero element */
78: for (i=0; i<mbs; i++) {
79: nz = ai[i+1] - ai[i];
80: j = ai[i];
81: /* printf("nz = %d, j=%d\n",nz,j); */
82: while (nz--) {
83: /* --- find row (=r[j]) and column (=aj[j]) in which to store a[j] ...*/
84: k = aj[j]; /* col. index */
85: /* printf("nz = %d, k=%d\n", nz,k); */
86: /* for entry that will be permuted into lower triangle, swap row and col. index */
87: if (rip[k] < rip[i]) aj[j] = i;
88: else k = i;
90: r[j] = k; j++;
91: nzr[k]++; /* increment count of nonzeros in that row */
92: }
93: }
95: /* Phase 2: Find new ai and permutation to apply to (aj,a).
96: Determine pointers (r) to delimit rows in permuted (aj,a).
97: Note: r is different from r used in phase 1.
98: At the end of this phase, (aj[j],a[j]) will be stored in
99: (aj[r(j)],a[r(j)]).
100: */
101: for (i=0; i<mbs; i++) {
102: ai[i+1] = ai[i] + nzr[i];
103: nzr[i] = ai[i+1];
104: }
106: /* determine where each (aj[j], a[j]) is stored in new (aj,a)
107: for each nonzero element (in reverse order) */
108: jmin = ai[0]; jmax = ai[mbs];
109: nz = jmax - jmin;
110: j = jmax-1;
111: while (nz--) {
112: i = r[j]; /* row value */
113: if (aj[j] == i) r[j] = ai[i]; /* put diagonal nonzero at beginning of row */
114: else { /* put off-diagonal nonzero in last unused location in row */
115: nzr[i]--; r[j] = nzr[i];
116: }
117: j--;
118: }
120: a->a2anew = aj + ai[mbs];
121: PetscArraycpy(a->a2anew,r,ai[mbs]);
123: /* Phase 3: permute (aj,a) to upper triangular form (wrt new ordering) */
124: for (j=jmin; j<jmax; j++) {
125: while (r[j] != j) {
126: k = r[j]; r[j] = r[k]; r[k] = k;
127: ajk = aj[k]; aj[k] = aj[j]; aj[j] = ajk;
128: /* ak = aa[k]; aa[k] = aa[j]; aa[j] = ak; */
129: }
130: }
131: ISRestoreIndices(perm,&rip);
133: a->inew = ai;
134: a->jnew = aj;
136: ISDestroy(&a->row);
137: ISDestroy(&a->icol);
138: PetscObjectReference((PetscObject)perm);
139: ISDestroy(&a->row);
140: a->row = perm;
141: PetscObjectReference((PetscObject)perm);
142: ISDestroy(&a->icol);
143: a->icol = perm;
145: PetscFree(nzr);
146: PetscFree(r);
147: return(0);
148: #endif
149: }