1 #include <../src/mat/impls/baij/seq/baij.h> 2 #include <../src/mat/impls/sbaij/seq/sbaij.h> 3 4 /* 5 This function is used before applying a 6 symmetric reordering to matrix A that is 7 in SBAIJ format. 8 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 13 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). 25 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; 35 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)); 45 46 PetscCall(ISInvertPermutation(perm,PETSC_DECIDE,&iperm)); 47 PetscCall(ISGetIndices(iperm,&riip)); 48 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)); 52 53 if (!a->inew) { 54 PetscCall(PetscMalloc2(mbs+1,&ai, 2*a->i[mbs],&aj)); 55 } else { 56 ai = a->inew; aj = a->jnew; 57 } 58 PetscCall(PetscArraycpy(ai,a->i,mbs+1)); 59 PetscCall(PetscArraycpy(aj,a->j,a->i[mbs])); 60 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 */ 67 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; 72 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; 85 86 r[j] = k; j++; 87 nzr[k]++; /* increment count of nonzeros in that row */ 88 } 89 } 90 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 } 101 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 } 115 116 a->a2anew = aj + ai[mbs]; 117 PetscCall(PetscArraycpy(a->a2anew,r,ai[mbs])); 118 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)); 128 129 a->inew = ai; 130 a->jnew = aj; 131 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; 140 141 PetscCall(PetscFree(nzr)); 142 PetscCall(PetscFree(r)); 143 PetscFunctionReturn(PETSC_SUCCESS); 144 #endif 145 } 146