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