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