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