1 /*$Id: sro.c,v 1.7 2000/09/07 16:11:21 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 This function is used before applying a 11 symmetric reordering to matrix A that is 12 in SBAIJ format. 13 14 The permutation is assumed to be symmetric, i.e., 15 P = P^T (= inv(P)), 16 so the permuted matrix P*A*inv(P)=P*A*P^T is ensured to be symmetric. 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 (ai, aj, a) for A such that all 32 nonzero A_(p(i),isp(k)) will be stored in the upper triangle. 33 Note: matrix A is not permuted by this function! 34 */ 35 #undef __FUNC__ 36 #define __FUNC__ "MatReorderingSeqSBAIJ" 37 int MatReIndexingSeqSBAIJ(Mat A,IS isp) 38 { 39 Mat_SeqSBAIJ *a=(Mat_SeqSBAIJ *)A->data; 40 int *r,ierr,i,mbs=a->mbs,*ai=a->i,*aj=a->j,*rip,*riip; 41 MatScalar *aa=a->a; 42 Scalar ak; 43 int *nzr,nz,jmin,jmax,j,k,ajk; 44 IS isip; /* inverse of isp */ 45 46 PetscFunctionBegin; 47 if (!mbs) PetscFunctionReturn(0); 48 49 ierr = ISGetIndices(isp,&rip);CHKERRQ(ierr); 50 ierr = ISInvertPermutation(isp,PETSC_DECIDE,&isip);CHKERRQ(ierr); 51 ierr = ISGetIndices(isip,&riip);CHKERRQ(ierr); 52 53 for (i=0; i<mbs; i++) { 54 if (rip[i] - riip[i] != 0) SETERRQ(1,1,"Non-symm. permutation, use symm. permutation or general matrix format"); 55 } 56 57 /* Phase 1: find row in which to store each nonzero (r) 58 initialize count of nonzeros to be stored in each row (nzr) */ 59 60 nzr = (int*)PetscMalloc(mbs*sizeof(int));CHKPTRQ(nzr); 61 r = (int*)PetscMalloc(ai[mbs]*sizeof(int));CHKPTRQ(r); 62 for (i=0; i<mbs; i++) nzr[i] = 0; 63 for (i=0; i<ai[mbs]; i++) r[i] = 0; 64 65 /* for each nonzero element */ 66 for (i=0; i<mbs; i++){ 67 nz = ai[i+1] - ai[i]; 68 j = ai[i]; 69 while (nz--){ 70 /* --- find row (=r[j]) and column (=aj[j]) in which to store a[j] ...*/ 71 k = aj[j]; 72 if (rip[k] < rip[i]) aj[j] = i; 73 else k = i; 74 r[j] = k; j++; 75 nzr[k] ++; /* increment count of nonzeros in that row */ 76 } 77 } 78 79 /* Phase 2: find new ai and permutation to apply to (aj,a) 80 determine pointers (r) to delimit rows in permuted (aj,a) */ 81 for (i=0; i<mbs; i++){ 82 ai[i+1] = ai[i] + nzr[i]; 83 nzr[i] = ai[i+1]; 84 } 85 86 /* determine where each (aj[j], a[j]) is stored in permuted (aj,a) 87 for each nonzero element (in reverse order) */ 88 jmin = ai[0]; jmax = ai[mbs]; 89 nz = jmax - jmin; 90 j = jmax-1; 91 while (nz--){ 92 i = r[j]; /* row value */ 93 if (aj[j] == i) r[j] = ai[i]; /* put diagonal nonzero at beginning of row */ 94 else { /* put off-diagonal nonzero in last unused location in row */ 95 nzr[i]--; r[j] = nzr[i]; 96 } 97 j--; 98 } 99 100 /* Phase 3: permute (aj,a) to upper triangular form (wrt new ordering) */ 101 for (j=jmin; j<jmax; j++){ 102 while (r[j] != j){ 103 k = r[j]; r[j] = r[k]; r[k] = k; 104 ajk = aj[k]; aj[k] = aj[j]; aj[j] = ajk; 105 ak = aa[k]; aa[k] = aa[j]; aa[j] = ak; 106 } 107 } 108 109 ierr = ISRestoreIndices(isp,&rip);CHKERRQ(ierr); 110 111 a->row = isp; 112 a->icol = isp; 113 ierr = PetscObjectReference((PetscObject)isp);CHKERRQ(ierr); 114 ierr = PetscObjectReference((PetscObject)isp);CHKERRQ(ierr); 115 PetscFunctionReturn(0); 116 } 117 118