1
2 #include <../src/mat/impls/baij/seq/baij.h>
3 #include <../src/mat/impls/sbaij/seq/sbaij.h>
4
5 /*
6 This function is used before applying a
7 symmetric reordering to matrix A that is
8 in SBAIJ format.
9
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
14
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).
26
27
28 -- output: new index set (inew, jnew) for A and a map a2anew that maps
29 values a to anew, such that all
30 nonzero A_(perm(i),iperm(k)) will be stored in the upper triangle.
31 Note: matrix A is not permuted by this function!
32 */
MatReorderingSeqSBAIJ(Mat A,IS perm)33 PetscErrorCode MatReorderingSeqSBAIJ(Mat A,IS perm)
34 {
35 Mat_SeqSBAIJ *a=(Mat_SeqSBAIJ*)A->data;
36 const PetscInt mbs=a->mbs;
37
38 PetscFunctionBegin;
39 if (!mbs) PetscFunctionReturn(0);
40 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Matrix reordering is not supported for sbaij matrix. Use aij format");
41 #if 0
42 PetscErrorCode ierr;
43 const PetscInt *rip,*riip;
44 PetscInt *ai,*aj,*r;
45 PetscInt *nzr,nz,jmin,jmax,j,k,ajk,i;
46 IS iperm; /* inverse of perm */
47 ierr = ISGetIndices(perm,&rip);CHKERRQ(ierr);
48
49 ierr = ISInvertPermutation(perm,PETSC_DECIDE,&iperm);CHKERRQ(ierr);
50 ierr = ISGetIndices(iperm,&riip);CHKERRQ(ierr);
51
52 for (i=0; i<mbs; i++) {
53 if (rip[i] != riip[i]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Non-symmetric permutation, use symmetric permutation for symmetric matrices");
54 }
55 ierr = ISRestoreIndices(iperm,&riip);CHKERRQ(ierr);
56 ierr = ISDestroy(&iperm);CHKERRQ(ierr);
57
58 if (!a->inew) {
59 ierr = PetscMalloc2(mbs+1,&ai, 2*a->i[mbs],&aj);CHKERRQ(ierr);
60 } else {
61 ai = a->inew; aj = a->jnew;
62 }
63 ierr = PetscArraycpy(ai,a->i,mbs+1);CHKERRQ(ierr);
64 ierr = PetscArraycpy(aj,a->j,a->i[mbs]);CHKERRQ(ierr);
65
66 /*
67 Phase 1: Find row index r in which to store each nonzero.
68 Initialize count of nonzeros to be stored in each row (nzr).
69 At the end of this phase, a nonzero a(*,*)=a(r(),aj())
70 s.t. a(perm(r),perm(aj)) will fall into upper triangle part.
71 */
72
73 ierr = PetscMalloc1(mbs,&nzr);CHKERRQ(ierr);
74 ierr = PetscMalloc1(ai[mbs],&r);CHKERRQ(ierr);
75 for (i=0; i<mbs; i++) nzr[i] = 0;
76 for (i=0; i<ai[mbs]; i++) r[i] = 0;
77
78 /* for each nonzero element */
79 for (i=0; i<mbs; i++) {
80 nz = ai[i+1] - ai[i];
81 j = ai[i];
82 /* printf("nz = %d, j=%d\n",nz,j); */
83 while (nz--) {
84 /* --- find row (=r[j]) and column (=aj[j]) in which to store a[j] ...*/
85 k = aj[j]; /* col. index */
86 /* printf("nz = %d, k=%d\n", nz,k); */
87 /* for entry that will be permuted into lower triangle, swap row and col. index */
88 if (rip[k] < rip[i]) aj[j] = i;
89 else k = i;
90
91 r[j] = k; j++;
92 nzr[k]++; /* increment count of nonzeros in that row */
93 }
94 }
95
96 /* Phase 2: Find new ai and permutation to apply to (aj,a).
97 Determine pointers (r) to delimit rows in permuted (aj,a).
98 Note: r is different from r used in phase 1.
99 At the end of this phase, (aj[j],a[j]) will be stored in
100 (aj[r(j)],a[r(j)]).
101 */
102 for (i=0; i<mbs; i++) {
103 ai[i+1] = ai[i] + nzr[i];
104 nzr[i] = ai[i+1];
105 }
106
107 /* determine where each (aj[j], a[j]) is stored in new (aj,a)
108 for each nonzero element (in reverse order) */
109 jmin = ai[0]; jmax = ai[mbs];
110 nz = jmax - jmin;
111 j = jmax-1;
112 while (nz--) {
113 i = r[j]; /* row value */
114 if (aj[j] == i) r[j] = ai[i]; /* put diagonal nonzero at beginning of row */
115 else { /* put off-diagonal nonzero in last unused location in row */
116 nzr[i]--; r[j] = nzr[i];
117 }
118 j--;
119 }
120
121 a->a2anew = aj + ai[mbs];
122 ierr = PetscArraycpy(a->a2anew,r,ai[mbs]);CHKERRQ(ierr);
123
124 /* Phase 3: permute (aj,a) to upper triangular form (wrt new ordering) */
125 for (j=jmin; j<jmax; j++) {
126 while (r[j] != j) {
127 k = r[j]; r[j] = r[k]; r[k] = k;
128 ajk = aj[k]; aj[k] = aj[j]; aj[j] = ajk;
129 /* ak = aa[k]; aa[k] = aa[j]; aa[j] = ak; */
130 }
131 }
132 ierr= ISRestoreIndices(perm,&rip);CHKERRQ(ierr);
133
134 a->inew = ai;
135 a->jnew = aj;
136
137 ierr = ISDestroy(&a->row);CHKERRQ(ierr);
138 ierr = ISDestroy(&a->icol);CHKERRQ(ierr);
139 ierr = PetscObjectReference((PetscObject)perm);CHKERRQ(ierr);
140 ierr = ISDestroy(&a->row);CHKERRQ(ierr);
141 a->row = perm;
142 ierr = PetscObjectReference((PetscObject)perm);CHKERRQ(ierr);
143 ierr = ISDestroy(&a->icol);CHKERRQ(ierr);
144 a->icol = perm;
145
146 ierr = PetscFree(nzr);CHKERRQ(ierr);
147 ierr = PetscFree(r);CHKERRQ(ierr);
148 PetscFunctionReturn(0);
149 #endif
150 }
151
152
153