1 /*$Id: sbaijov.c,v 1.65 2001/08/06 21:15:42 bsmith Exp $*/ 2 3 /* 4 Routines to compute overlapping regions of a parallel MPI matrix. 5 Used for finding submatrices that were shared across processors. 6 */ 7 #include "src/mat/impls/sbaij/mpi/mpisbaij.h" 8 #include "petscbt.h" 9 10 static int MatIncreaseOverlap_MPISBAIJ_Once(Mat,int,IS *); 11 static int MatIncreaseOverlap_MPISBAIJ_Local(Mat,int *,int,int **,PetscBT*); 12 13 #undef __FUNCT__ 14 #define __FUNCT__ "MatIncreaseOverlap_MPISBAIJ" 15 int MatIncreaseOverlap_MPISBAIJ(Mat C,int is_max,IS is[],int ov) 16 { 17 Mat_MPISBAIJ *c = (Mat_MPISBAIJ*)C->data; 18 int i,ierr,N=C->N, bs=c->bs; 19 IS *is_new; 20 21 PetscFunctionBegin; 22 ierr = PetscMalloc(is_max*sizeof(IS),&is_new);CHKERRQ(ierr); 23 /* Convert the indices into block format */ 24 ierr = ISCompressIndicesGeneral(N,bs,is_max,is,is_new);CHKERRQ(ierr); 25 if (ov < 0){ SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Negative overlap specified\n");} 26 for (i=0; i<ov; ++i) { 27 ierr = MatIncreaseOverlap_MPISBAIJ_Once(C,is_max,is_new);CHKERRQ(ierr); 28 } 29 for (i=0; i<is_max; i++) {ierr = ISDestroy(is[i]);CHKERRQ(ierr);} 30 ierr = ISExpandIndicesGeneral(N,bs,is_max,is_new,is);CHKERRQ(ierr); 31 for (i=0; i<is_max; i++) {ierr = ISDestroy(is_new[i]);CHKERRQ(ierr);} 32 ierr = PetscFree(is_new);CHKERRQ(ierr); 33 PetscFunctionReturn(0); 34 } 35 36 typedef enum {MINE,OTHER} WhoseOwner; 37 /* data1, odata1 and odata2 are packed in the format (for communication): 38 data[0] = is_max, no of is 39 data[1] = size of is[0] 40 ... 41 data[is_max] = size of is[is_max-1] 42 data[is_max + 1] = data(is[0]) 43 ... 44 data[is_max+1+sum(size of is[k]), k=0,...,i-1] = data(is[i]) 45 ... 46 data2 is packed in the format (for creating output is[]): 47 data[0] = is_max, no of is 48 data[1] = size of is[0] 49 ... 50 data[is_max] = size of is[is_max-1] 51 data[is_max + 1] = data(is[0]) 52 ... 53 data[is_max + 1 + Mbs*i) = data(is[i]) 54 ... 55 */ 56 #undef __FUNCT__ 57 #define __FUNCT__ "MatIncreaseOverlap_MPISBAIJ_Once" 58 static int MatIncreaseOverlap_MPISBAIJ_Once(Mat C,int is_max,IS is[]) 59 { 60 Mat_MPISBAIJ *c = (Mat_MPISBAIJ*)C->data; 61 int len,*idx_i,isz,col,*n,*data1,*data2,*data2_i,*data,*data_i, 62 size,rank,Mbs,i,j,k,ierr,nrqs,*odata1,*odata2, 63 tag1,tag2,flag,proc_id,**odata2_ptr; 64 char *t_p; 65 MPI_Comm comm; 66 MPI_Request *s_waits,r_req; 67 MPI_Status *s_status,r_status; 68 PetscBT *table; /* mark indices of this processor's is[] */ 69 PetscBT table_i; 70 PetscBT otable; /* mark indices of other processors' is[] */ 71 PetscFunctionBegin; 72 73 comm = C->comm; 74 size = c->size; 75 rank = c->rank; 76 Mbs = c->Mbs; 77 78 ierr = PetscObjectGetNewTag((PetscObject)C,&tag1);CHKERRQ(ierr); 79 ierr = PetscObjectGetNewTag((PetscObject)C,&tag2);CHKERRQ(ierr); 80 81 /* create tables for marking indices */ 82 len = is_max*sizeof(PetscBT) + (Mbs/PETSC_BITS_PER_BYTE+1)*is_max*sizeof(char) + 1; 83 ierr = PetscMalloc(len,&table);CHKERRQ(ierr); 84 t_p = (char *)(table + is_max); 85 for (i=0; i<is_max; i++) { 86 table[i] = t_p + (Mbs/PETSC_BITS_PER_BYTE+1)*i; 87 } 88 ierr = PetscBTCreate(Mbs,otable);CHKERRQ(ierr); 89 90 /* 1. Send this processor's is[] to all other processors */ 91 /*-------------------------------------------------------*/ 92 /* Allocate Memory for outgoing messages */ 93 len = is_max*sizeof(int); 94 ierr = PetscMalloc(len,&n);CHKERRQ(ierr); 95 96 len = 1 + is_max; 97 for (i=0; i<is_max; i++) { 98 ierr = ISGetLocalSize(is[i],&n[i]);CHKERRQ(ierr); 99 len += n[i]; 100 } 101 ierr = PetscMalloc(len*sizeof(int),&data1);CHKERRQ(ierr); 102 103 /* Form the outgoing messages */ 104 data1[0] = is_max; 105 k = is_max + 1; 106 for (i=0; i<is_max; i++) { 107 data1[1+i] = n[i]; 108 ierr = ISGetIndices(is[i],&idx_i);CHKERRQ(ierr); 109 for (j=0; j<data1[i+1]; j++){ 110 data1[k++] = idx_i[j]; 111 } 112 ierr = ISRestoreIndices(is[i],&idx_i);CHKERRQ(ierr); 113 ierr = ISDestroy(is[i]);CHKERRQ(ierr); 114 } 115 if (k != len) SETERRQ2(1,"Error on forming the outgoing messages: k %d != len %d",k,len); 116 ierr = PetscFree(n);CHKERRQ(ierr); 117 118 /* Now post the sends */ 119 ierr = PetscMalloc(size*sizeof(MPI_Request),&s_waits);CHKERRQ(ierr); 120 k = 0; 121 for (proc_id=0; proc_id<size; ++proc_id) { /* send data1 to processor [proc_id] */ 122 if (proc_id != rank){ 123 ierr = MPI_Isend(data1,len,MPI_INT,proc_id,tag1,comm,s_waits+k);CHKERRQ(ierr); 124 /* printf(" [%d] send %d msg to [%d], data1: \n",rank,len,proc_id); */ 125 k++; 126 } 127 } 128 129 /* 2. Do local work on this processor's is[] */ 130 /*-------------------------------------------*/ 131 ierr = MatIncreaseOverlap_MPISBAIJ_Local(C,data1,MINE,&data,table);CHKERRQ(ierr); 132 133 /* 3. Receive other's is[] and process. Then send back */ 134 /*-----------------------------------------------------*/ 135 /* Sending this processor's is[] is done */ 136 nrqs = size-1; 137 ierr = PetscMalloc(size*sizeof(MPI_Status),&s_status);CHKERRQ(ierr); 138 ierr = MPI_Waitall(nrqs,s_waits,s_status);CHKERRQ(ierr); 139 140 ierr = PetscMalloc(size*sizeof(int**),&odata2_ptr);CHKERRQ(ierr); 141 k = 0; 142 do { 143 /* Receive messages */ 144 ierr = MPI_Iprobe(MPI_ANY_SOURCE,tag1,comm,&flag,&r_status);CHKERRQ(ierr); 145 if (flag){ 146 ierr = MPI_Get_count(&r_status,MPI_INT,&len);CHKERRQ(ierr); 147 proc_id = r_status.MPI_SOURCE; 148 ierr = PetscMalloc(len*sizeof(int),&odata1);CHKERRQ(ierr); 149 ierr = MPI_Irecv(odata1,len,MPI_INT,proc_id,r_status.MPI_TAG,comm,&r_req);CHKERRQ(ierr); 150 /* printf(" [%d] recv %d msg from [%d]\n",rank,len,proc_id); */ 151 152 /* Process messages */ 153 ierr = MatIncreaseOverlap_MPISBAIJ_Local(C,odata1,OTHER,&odata2_ptr[k],&otable);CHKERRQ(ierr); 154 odata2 = odata2_ptr[k]; 155 len = 1 + odata2[0]; 156 for (i=0; i<odata2[0]; i++){ 157 len += odata2[1 + i]; 158 } 159 160 /* Send messages back */ 161 ierr = MPI_Isend(odata2,len,MPI_INT,proc_id,tag2,comm,s_waits+k);CHKERRQ(ierr); 162 /* printf(" [%d] send %d msg back to [%d] \n",rank,len,proc_id); */ 163 164 ierr = PetscFree(odata1);CHKERRQ(ierr); 165 k++; 166 } 167 } while (k < nrqs); 168 169 /* 4. Receive work done on other processors, then merge */ 170 /*--------------------------------------------------------*/ 171 /* Allocate memory for incoming data */ 172 len = (1+is_max*(Mbs+1)); 173 ierr = PetscMalloc(len*sizeof(int),&data2);CHKERRQ(ierr); 174 175 /* Sending others' is[] is done */ 176 ierr = MPI_Waitall(nrqs,s_waits,s_status);CHKERRQ(ierr); 177 for (k=0; k<nrqs; k++){ 178 ierr = PetscFree(odata2_ptr[k]);CHKERRQ(ierr); 179 } 180 ierr = PetscFree(odata2_ptr);CHKERRQ(ierr); 181 182 k = 0; 183 do { 184 /* Receive messages */ 185 ierr = MPI_Iprobe(MPI_ANY_SOURCE,tag2,comm,&flag,&r_status); 186 if (flag){ 187 ierr = MPI_Get_count(&r_status,MPI_INT,&len);CHKERRQ(ierr); 188 proc_id = r_status.MPI_SOURCE; 189 ierr = MPI_Irecv(data2,len,MPI_INT,proc_id,r_status.MPI_TAG,comm,&r_req);CHKERRQ(ierr); 190 /* printf(" [%d] recv %d msg from [%d], data2:\n",rank,len,proc_id); */ 191 192 /* Add data2 into data */ 193 data2_i = data2 + 1 + is_max; 194 for (i=0; i<is_max; i++){ 195 table_i = table[i]; 196 data_i = data + 1 + is_max + Mbs*i; 197 isz = data[1+i]; 198 for (j=0; j<data2[1+i]; j++){ 199 col = data2_i[j]; 200 if (!PetscBTLookupSet(table_i,col)) {data_i[isz++] = col;} 201 } 202 data[1+i] = isz; 203 if (i < is_max - 1) data2_i += data2[1+i]; 204 } 205 k++; 206 } 207 } while (k < nrqs); 208 209 /* 5. Create new is[] */ 210 /*--------------------*/ 211 for (i=0; i<is_max; i++) { 212 data_i = data + 1 + is_max + Mbs*i; 213 ierr = ISCreateGeneral(PETSC_COMM_SELF,data[1+i],data_i,is+i);CHKERRQ(ierr); 214 } 215 ierr = PetscFree(data1);CHKERRQ(ierr); 216 ierr = PetscFree(data2);CHKERRQ(ierr); 217 ierr = PetscFree(data);CHKERRQ(ierr); 218 ierr = PetscFree(s_waits);CHKERRQ(ierr); 219 ierr = PetscFree(s_status);CHKERRQ(ierr); 220 ierr = PetscFree(table);CHKERRQ(ierr); 221 ierr = PetscBTDestroy(otable);CHKERRQ(ierr); 222 PetscFunctionReturn(0); 223 } 224 225 #undef __FUNCT__ 226 #define __FUNCT__ "MatIncreaseOverlap_MPISBAIJ_Local" 227 /* 228 MatIncreaseOverlap_MPISBAIJ_Local - Called by MatIncreaseOverlap, to do 229 the work on the local processor. 230 231 Inputs: 232 C - MAT_MPISBAIJ; 233 data - holds is[]. See MatIncreaseOverlap_MPISBAIJ_Once() for the format. 234 whose - whose is[] to be processed, 235 MINE: this processor's is[] 236 OTHER: other processor's is[] 237 Output: 238 data_new - whose = MINE: 239 holds input and newly found indices in the same format as data 240 whose = OTHER: 241 only holds the newly found indices 242 table - table[i]: mark the indices of is[i], i=0,...,is_max. Used only in the case 'whose=MINE'. 243 */ 244 static int MatIncreaseOverlap_MPISBAIJ_Local(Mat C,int *data,int whose,int **data_new,PetscBT *table) 245 { 246 Mat_MPISBAIJ *c = (Mat_MPISBAIJ*)C->data; 247 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)(c->A)->data; 248 Mat_SeqBAIJ *b = (Mat_SeqBAIJ*)(c->B)->data; 249 int ierr,row,mbs,Mbs,*nidx,*nidx_i,col,isz,isz0,*ai,*aj,*bi,*bj,*garray,rstart,l; 250 int a_start,a_end,b_start,b_end,i,j,k,is_max,*idx_i,n; 251 PetscBT table0; /* mark the indices of input is[] for look up */ 252 PetscBT table_i; /* poits to i-th table. When whose=OTHER, a single table is used for all is[] */ 253 254 PetscFunctionBegin; 255 Mbs = c->Mbs; mbs = a->mbs; 256 ai = a->i; aj = a->j; 257 bi = b->i; bj = b->j; 258 garray = c->garray; 259 rstart = c->rstart; 260 is_max = data[0]; 261 262 ierr = PetscBTCreate(Mbs,table0);CHKERRQ(ierr); 263 264 ierr = PetscMalloc((1+is_max*(Mbs+1))*sizeof(int),&nidx);CHKERRQ(ierr); 265 nidx[0] = is_max; 266 267 idx_i = data + is_max + 1; /* ptr to input is[0] array */ 268 nidx_i = nidx + is_max + 1; /* ptr to output is[0] array */ 269 for (i=0; i<is_max; i++) { /* for each is */ 270 isz = 0; 271 n = data[1+i]; /* size of input is[i] */ 272 273 /* initialize table_i, set table0 */ 274 if (whose == MINE){ /* process this processor's is[] */ 275 table_i = table[i]; 276 nidx_i = nidx + 1+ is_max + Mbs*i; 277 } else { /* process other processor's is[] - only use one temp table */ 278 table_i = *table; 279 } 280 ierr = PetscBTMemzero(Mbs,table_i);CHKERRQ(ierr); 281 ierr = PetscBTMemzero(Mbs,table0);CHKERRQ(ierr); 282 if (n > 0) { 283 isz0 = 0; 284 for (j=0; j<n; j++){ 285 col = idx_i[j]; 286 if (col >= Mbs) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"index col %d >= Mbs %d",col,Mbs); 287 if(!PetscBTLookupSet(table_i,col)) { 288 ierr = PetscBTSet(table0,col);CHKERRQ(ierr); 289 if (whose == MINE) {nidx_i[isz0] = col;} 290 isz0++; 291 } 292 } 293 294 if (whose == MINE) {isz = isz0;} 295 k = 0; /* no. of indices from input is[i] that have been examined */ 296 for (row=0; row<mbs; row++){ 297 a_start = ai[row]; a_end = ai[row+1]; 298 b_start = bi[row]; b_end = bi[row+1]; 299 if (PetscBTLookup(table0,row+rstart)){ /* row is on input is[i]: 300 do row search: collect all col in this row */ 301 for (l = a_start; l<a_end ; l++){ /* Amat */ 302 col = aj[l] + rstart; 303 if (!PetscBTLookupSet(table_i,col)) {nidx_i[isz++] = col;} 304 } 305 for (l = b_start; l<b_end ; l++){ /* Bmat */ 306 col = garray[bj[l]]; 307 if (!PetscBTLookupSet(table_i,col)) {nidx_i[isz++] = col;} 308 } 309 k++; 310 if (k >= isz0) break; /* for (row=0; row<mbs; row++) */ 311 } else { /* row is not on input is[i]: 312 do col serach: add row onto nidx_i if there is a col in nidx_i */ 313 for (l = a_start; l<a_end ; l++){ /* Amat */ 314 col = aj[l] + rstart; 315 if (PetscBTLookup(table0,col)){ 316 if (!PetscBTLookupSet(table_i,row+rstart)) {nidx_i[isz++] = row+rstart;} 317 break; /* for l = start; l<end ; l++) */ 318 } 319 } 320 for (l = b_start; l<b_end ; l++){ /* Bmat */ 321 col = garray[bj[l]]; 322 if (PetscBTLookup(table0,col)){ 323 if (!PetscBTLookupSet(table_i,row+rstart)) {nidx_i[isz++] = row+rstart;} 324 break; /* for l = start; l<end ; l++) */ 325 } 326 } 327 } 328 } 329 } /* if (n > 0) */ 330 331 if (i < is_max - 1){ 332 idx_i += n; /* ptr to input is[i+1] array */ 333 nidx_i += isz; /* ptr to output is[i+1] array */ 334 } 335 nidx[1+i] = isz; /* size of new is[i] */ 336 } /* for each is */ 337 *data_new = nidx; 338 ierr = PetscBTDestroy(table0);CHKERRQ(ierr); 339 340 PetscFunctionReturn(0); 341 } 342 343 344