1 #include "petscsys.h" 2 #include "src/mat/impls/baij/seq/baij.h" 3 #include "src/vec/vecimpl.h" 4 #include "src/inline/spops.h" 5 #include "sbaij.h" 6 7 /* 8 Symmetric reordering the index of sparse symmetric matrix A 9 so that P*A*P^T is also stored symmetrically. 10 The permutation needs to be symmetric, i.e., P = P^T = inv(P). 11 -- modified from sor.f of YSMP 12 -- idea: reorder (ai, aj, a) (not A!) to ensure that all nonzero A_(p(i),isp(k)) 13 are stored in the upper triangle of P*A*P^T 14 */ 15 #undef __FUNC__ 16 #define __FUNC__ "MatReordering_SeqSBAIJ" 17 int MatReordering_SeqSBAIJ(Mat A,IS isp) 18 { 19 Mat_SeqSBAIJ *a=(Mat_SeqSBAIJ *)A->data; 20 int *r,ierr,i,mbs=a->mbs,*ai=a->i,*aj=a->j,*rip,*riip; 21 MatScalar *aa=a->a; 22 Scalar ak; 23 int *nzr,nz,jmin,jmax,j,k,ajk; 24 IS isip; /* inverse of isp */ 25 26 PetscFunctionBegin; 27 if (!mbs) PetscFunctionReturn(0); 28 29 ierr = ISGetIndices(isp,&rip);CHKERRQ(ierr); 30 ierr = ISInvertPermutation(isp,PETSC_DECIDE,&isip);CHKERRQ(ierr); 31 ierr = ISGetIndices(isip,&riip);CHKERRQ(ierr); 32 33 for (i=0; i<mbs; i++) { 34 if ( rip[i] - riip[i] != 0 ) { 35 printf("Non-symm. permutation, use symm. permutation or general matrix format\n"); 36 break; 37 } 38 } 39 40 /* Phase 1: find row in which to store each nonzero (r) 41 initialize count of nonzeros to be stored in each row (nzr) */ 42 43 nzr = (int*)PetscMalloc(mbs*sizeof(int));CHKPTRQ(nzr); 44 r = (int*)PetscMalloc(ai[mbs]*sizeof(int));CHKPTRQ(r); 45 for (i=0; i<mbs; i++) nzr[i] = 0; 46 for (i=0; i<ai[mbs]; i++) r[i] = 0; 47 48 /* for each nonzero element */ 49 for (i=0; i<mbs; i++){ 50 nz = ai[i+1] - ai[i]; 51 j = ai[i]; 52 while (nz--){ 53 /* --- find row (=r[j]) and column (=aj[j]) in which to store a[j] ...*/ 54 k = aj[j]; 55 if (rip[k] < rip[i]) aj[j] = i; 56 else k = i; 57 r[j] = k; j++; 58 nzr[k] ++; /* increment count of nonzeros in that row */ 59 } 60 } 61 62 /* Phase 2: find new ai and permutation to apply to (aj,a) 63 determine pointers (r) to delimit rows in permuted (aj,a) */ 64 for (i=0; i<mbs; i++){ 65 ai[i+1] = ai[i] + nzr[i]; 66 nzr[i] = ai[i+1]; 67 } 68 69 /* determine where each (aj[j], a[j]) is stored in permuted (aj,a) 70 for each nonzero element (in reverse order) */ 71 jmin = ai[0]; jmax = ai[mbs]; 72 nz = jmax - jmin; 73 j = jmax-1; 74 while (nz--){ 75 i = r[j]; /* row value */ 76 if (aj[j] == i) r[j] = ai[i]; /* put diagonal nonzero at beginning of row */ 77 else { /* put off-diagonal nonzero in last unused location in row */ 78 nzr[i]--; r[j] = nzr[i]; 79 } 80 j--; 81 } 82 83 /* Phase 3: permute (aj,a) to upper triangular form (wrt new ordering) */ 84 for (j=jmin; j<jmax; j++){ 85 while (r[j] != j){ 86 k = r[j]; r[j] = r[k]; r[k] = k; 87 ajk = aj[k]; aj[k] = aj[j]; aj[j] = ajk; 88 ak = aa[k]; aa[k] = aa[j]; aa[j] = ak; 89 } 90 } 91 92 ierr = ISRestoreIndices(isp,&rip);CHKERRQ(ierr); 93 94 a->row = isp; 95 a->col = isp; 96 a->icol = isp; 97 ierr = PetscObjectReference((PetscObject)isp); CHKERRQ(ierr); 98 PetscFunctionReturn(0); 99 } 100 101