1 #ifndef lint 2 static char vcid[] = "$Id: baij2.c,v 1.2 1996/04/30 18:59:24 balay Exp balay $"; 3 #endif 4 5 #include "baij.h" 6 #include "petsc.h" 7 #include "src/inline/bitarray.h" 8 9 int MatIncreaseOverlap_SeqBAIJ(Mat A,int is_max,IS *is,int ov) 10 { 11 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *) A->data; 12 int row, i,j,k,l,m,n, *idx,ierr, *nidx, isz, val; 13 int start, end, *ai, *aj; 14 char *table; 15 16 m = a->mbs; 17 ai = a->i; 18 aj = a->j; 19 20 if (ov < 0) SETERRQ(1,"MatIncreaseOverlap_SeqBAIJ: illegal overlap value used"); 21 22 table = (char *) PetscMalloc((m/BITSPERBYTE +1)*sizeof(char)); CHKPTRQ(table); 23 nidx = (int *) PetscMalloc((m+1)*sizeof(int)); CHKPTRQ(nidx); 24 25 for ( i=0; i<is_max; i++ ) { 26 /* Initialise the two local arrays */ 27 isz = 0; 28 PetscMemzero(table,(m/BITSPERBYTE +1)*sizeof(char)); 29 30 /* Extract the indices, assume there can be duplicate entries */ 31 ierr = ISGetIndices(is[i],&idx); CHKERRQ(ierr); 32 ierr = ISGetSize(is[i],&n); CHKERRQ(ierr); 33 34 /* Enter these into the temp arrays i.e mark table[row], enter row into new index */ 35 for ( j=0; j<n ; ++j){ 36 if (idx[j]>m) SETERRQ(1,"MatIncreaseOverlap_SeqBAIJ: index greater than mat-dim"); 37 if(!BT_LOOKUP(table, idx[j])) { nidx[isz++] = idx[j];} 38 } 39 ierr = ISRestoreIndices(is[i],&idx); CHKERRQ(ierr); 40 ierr = ISDestroy(is[i]); CHKERRQ(ierr); 41 42 k = 0; 43 for ( j=0; j<ov; j++){ /* for each overlap*/ 44 n = isz; 45 for ( ; k<n ; k++){ /* do only those rows in nidx[k], which are not done yet */ 46 row = nidx[k]; 47 start = ai[row]; 48 end = ai[row+1]; 49 for ( l = start; l<end ; l++){ 50 val = aj[l]; 51 if (!BT_LOOKUP(table,val)) {nidx[isz++] = val;} 52 } 53 } 54 } 55 ierr = ISCreateSeq(MPI_COMM_SELF, isz, nidx, (is+i)); CHKERRQ(ierr); 56 } 57 PetscFree(table); 58 PetscFree(nidx); 59 return 0; 60 } 61 int MatGetSubMatrix_SeqBAIJ(Mat A,IS isrow,IS iscol,MatGetSubMatrixCall scall,Mat *B) 62 { 63 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *) A->data,*c; 64 int nznew, *smap, i, k, kstart, kend, ierr, oldcols = a->nbs,*lens; 65 int row,mat_i,*mat_j,tcol,first,step,*mat_ilen; 66 register int sum,lensi; 67 int *irow, *icol, nrows, ncols,*ssmap,bs=a->bs, bs2=a->bs2; 68 int *starts,*j_new,*i_new,*aj = a->j, *ai = a->i,ii,*ailen = a->ilen; 69 Scalar *a_new,*mat_a; 70 Mat C; 71 72 ierr = ISSorted(iscol,(PetscTruth*)&i); 73 if (!i) SETERRQ(1,"MatGetSubmatrices_SeqBAIJ:IS is not sorted"); 74 75 ierr = ISGetIndices(isrow,&irow); CHKERRQ(ierr); 76 ierr = ISGetSize(isrow,&nrows); CHKERRQ(ierr); 77 ierr = ISGetSize(iscol,&ncols); CHKERRQ(ierr); 78 79 if (ISStrideGetInfo(iscol,&first,&step) && step == 1) { /* no need to sort */ 80 /* special case of contiguous rows */ 81 lens = (int *) PetscMalloc((ncols+nrows+1)*sizeof(int)); CHKPTRQ(lens); 82 starts = lens + ncols; 83 /* loop over new rows determining lens and starting points */ 84 for (i=0; i<nrows; i++) { 85 kstart = ai[irow[i]]; 86 kend = kstart + ailen[irow[i]]; 87 for ( k=kstart; k<kend; k++ ) { 88 if (aj[k] >= first) { 89 starts[i] = k; 90 break; 91 } 92 } 93 sum = 0; 94 while (k < kend) { 95 if (aj[k++] >= first+ncols) break; 96 sum++; 97 } 98 lens[i] = sum; 99 } 100 /* create submatrix */ 101 if (scall == MAT_REUSE_MATRIX) { 102 c = (Mat_SeqBAIJ *)((*B)->data); 103 104 if (c->mbs!=nrows || c->nbs!=ncols || c->bs!=bs) 105 SETERRQ(1,"MatGetSubMatrix_SeqBAIJ:"); 106 ierr = MatZeroEntries(*B); CHKERRQ(ierr); 107 C = *B; 108 } 109 else { 110 ierr = MatCreateSeqBAIJ(A->comm,bs,nrows*bs,ncols*bs,0,lens,&C);CHKERRQ(ierr); 111 } 112 c = (Mat_SeqBAIJ*) C->data; 113 114 /* loop over rows inserting into submatrix */ 115 a_new = c->a; 116 j_new = c->j; 117 i_new = c->i; 118 i_new[0] = 0; 119 for (i=0; i<nrows; i++) { 120 ii = starts[i]; 121 lensi = lens[i]; 122 for ( k=0; k<lensi; k++ ) { 123 *j_new++ = aj[ii+k] - first; 124 } 125 PetscMemcpy(a_new,a->a + starts[i],lensi*bs2*sizeof(Scalar)); 126 a_new += lensi*bs2; 127 i_new[i+1] = i_new[i] + lensi; 128 c->ilen[i] = lensi; 129 } 130 PetscFree(lens); 131 } 132 else { 133 ierr = ISGetIndices(iscol,&icol); CHKERRQ(ierr); 134 smap = (int *) PetscMalloc((1+oldcols)*sizeof(int)); CHKPTRQ(smap); 135 ssmap = smap; 136 lens = (int *) PetscMalloc((1+nrows)*sizeof(int)); CHKPTRQ(lens); 137 PetscMemzero(smap,oldcols*sizeof(int)); 138 for ( i=0; i<ncols; i++ ) smap[icol[i]] = i+1; 139 /* determine lens of each row */ 140 for (i=0; i<nrows; i++) { 141 kstart = ai[irow[i]]; 142 kend = kstart + a->ilen[irow[i]]; 143 lens[i] = 0; 144 for ( k=kstart; k<kend; k++ ) { 145 if (ssmap[aj[k]]) { 146 lens[i]++; 147 } 148 } 149 } 150 /* Create and fill new matrix */ 151 if (scall == MAT_REUSE_MATRIX) { 152 c = (Mat_SeqBAIJ *)((*B)->data); 153 154 if (c->mbs!=nrows || c->nbs!=ncols || c->bs!=bs) 155 SETERRQ(1,"MatGetSubMatrix_SeqBAIJ:"); 156 if (PetscMemcmp(c->ilen,lens, c->mbs *sizeof(int))) { 157 SETERRQ(1,"MatGetSubmatrix_SeqBAIJ:Cannot reuse matrix. wrong no of nonzeros"); 158 } 159 PetscMemzero(c->ilen,c->mbs*sizeof(int)); 160 C = *B; 161 } 162 else { 163 ierr = MatCreateSeqBAIJ(A->comm,bs,nrows*bs,ncols*bs,0,lens,&C);CHKERRQ(ierr); 164 } 165 c = (Mat_SeqBAIJ *)(C->data); 166 for (i=0; i<nrows; i++) { 167 row = irow[i]; 168 nznew = 0; 169 kstart = ai[row]; 170 kend = kstart + a->ilen[row]; 171 mat_i = c->i[i]; 172 mat_j = c->j + mat_i; 173 mat_a = c->a + mat_i; 174 mat_ilen = c->ilen + i; 175 for ( k=kstart; k<kend; k++ ) { 176 if ((tcol=ssmap[a->j[k]])) { 177 *mat_j++ = tcol - 1; 178 PetscMemcpy(mat_a,a->a+k*bs2,bs2*sizeof(Scalar)); mat_a+=bs2; 179 (*mat_ilen)++; 180 181 } 182 } 183 } 184 /* Free work space */ 185 ierr = ISRestoreIndices(iscol,&icol); CHKERRQ(ierr); 186 PetscFree(smap); PetscFree(lens); 187 } 188 ierr = MatAssemblyBegin(C,FINAL_ASSEMBLY); CHKERRQ(ierr); 189 ierr = MatAssemblyEnd(C,FINAL_ASSEMBLY); CHKERRQ(ierr); 190 191 ierr = ISRestoreIndices(isrow,&irow); CHKERRQ(ierr); 192 *B = C; 193 return 0; 194 } 195 196 int MatGetSubMatrices_SeqBAIJ(Mat A,int n, IS *irow,IS *icol,MatGetSubMatrixCall scall, 197 Mat **B) 198 { 199 int ierr,i; 200 201 if (scall == MAT_INITIAL_MATRIX) { 202 *B = (Mat *) PetscMalloc( (n+1)*sizeof(Mat) ); CHKPTRQ(*B); 203 } 204 205 for ( i=0; i<n; i++ ) { 206 ierr = MatGetSubMatrix(A,irow[i],icol[i],scall,&(*B)[i]); CHKERRQ(ierr); 207 } 208 return 0; 209 } 210