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 -- 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; 36 37 PetscFunctionBegin; 38 if (!mbs) PetscFunctionReturn(0); 39 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Matrix reordering is not supported for sbaij matrix. Use aij format"); 40 #if 0 41 PetscErrorCode ierr; 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 ierr = ISGetIndices(perm,&rip);CHKERRQ(ierr); 47 48 ierr = ISInvertPermutation(perm,PETSC_DECIDE,&iperm);CHKERRQ(ierr); 49 ierr = ISGetIndices(iperm,&riip);CHKERRQ(ierr); 50 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 ierr = ISRestoreIndices(iperm,&riip);CHKERRQ(ierr); 55 ierr = ISDestroy(&iperm);CHKERRQ(ierr); 56 57 if (!a->inew) { 58 ierr = PetscMalloc2(mbs+1,&ai, 2*a->i[mbs],&aj);CHKERRQ(ierr); 59 } else { 60 ai = a->inew; aj = a->jnew; 61 } 62 ierr = PetscArraycpy(ai,a->i,mbs+1);CHKERRQ(ierr); 63 ierr = PetscArraycpy(aj,a->j,a->i[mbs]);CHKERRQ(ierr); 64 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 */ 71 72 ierr = PetscMalloc1(mbs,&nzr);CHKERRQ(ierr); 73 ierr = PetscMalloc1(ai[mbs],&r);CHKERRQ(ierr); 74 for (i=0; i<mbs; i++) nzr[i] = 0; 75 for (i=0; i<ai[mbs]; i++) r[i] = 0; 76 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; 89 90 r[j] = k; j++; 91 nzr[k]++; /* increment count of nonzeros in that row */ 92 } 93 } 94 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 } 105 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 } 119 120 a->a2anew = aj + ai[mbs]; 121 ierr = PetscArraycpy(a->a2anew,r,ai[mbs]);CHKERRQ(ierr); 122 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 ierr = ISRestoreIndices(perm,&rip);CHKERRQ(ierr); 132 133 a->inew = ai; 134 a->jnew = aj; 135 136 ierr = ISDestroy(&a->row);CHKERRQ(ierr); 137 ierr = ISDestroy(&a->icol);CHKERRQ(ierr); 138 ierr = PetscObjectReference((PetscObject)perm);CHKERRQ(ierr); 139 ierr = ISDestroy(&a->row);CHKERRQ(ierr); 140 a->row = perm; 141 ierr = PetscObjectReference((PetscObject)perm);CHKERRQ(ierr); 142 ierr = ISDestroy(&a->icol);CHKERRQ(ierr); 143 a->icol = perm; 144 145 ierr = PetscFree(nzr);CHKERRQ(ierr); 146 ierr = PetscFree(r);CHKERRQ(ierr); 147 PetscFunctionReturn(0); 148 #endif 149 } 150 151