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 /* this function is sasme as MatCompressIndicesGeneral_MPIBAIJ -- should be removed! */ 14 #undef __FUNCT__ 15 #define __FUNCT__ "MatCompressIndicesGeneral_MPISBAIJ" 16 static int MatCompressIndicesGeneral_MPISBAIJ(Mat C,int imax,const IS is_in[],IS is_out[]) 17 { 18 Mat_MPISBAIJ *baij = (Mat_MPISBAIJ*)C->data; 19 int ierr,isz,bs = baij->bs,n,i,j,*idx,ival; 20 #if defined (PETSC_USE_CTABLE) 21 PetscTable gid1_lid1; 22 int tt, gid1, *nidx; 23 PetscTablePosition tpos; 24 #else 25 int Nbs,*nidx; 26 PetscBT table; 27 #endif 28 29 PetscFunctionBegin; 30 /* printf(" ...MatCompressIndicesGeneral_MPISBAIJ is called ...\n"); */ 31 #if defined (PETSC_USE_CTABLE) 32 ierr = PetscTableCreate(baij->mbs,&gid1_lid1);CHKERRQ(ierr); 33 #else 34 Nbs = baij->Nbs; 35 ierr = PetscMalloc((Nbs+1)*sizeof(int),&nidx);CHKERRQ(ierr); 36 ierr = PetscBTCreate(Nbs,table);CHKERRQ(ierr); 37 #endif 38 for (i=0; i<imax; i++) { 39 isz = 0; 40 #if defined (PETSC_USE_CTABLE) 41 ierr = PetscTableRemoveAll(gid1_lid1);CHKERRQ(ierr); 42 #else 43 ierr = PetscBTMemzero(Nbs,table);CHKERRQ(ierr); 44 #endif 45 ierr = ISGetIndices(is_in[i],&idx);CHKERRQ(ierr); 46 ierr = ISGetLocalSize(is_in[i],&n);CHKERRQ(ierr); 47 for (j=0; j<n ; j++) { 48 ival = idx[j]/bs; /* convert the indices into block indices */ 49 #if defined (PETSC_USE_CTABLE) 50 ierr = PetscTableFind(gid1_lid1,ival+1,&tt);CHKERRQ(ierr); 51 if (!tt) { 52 ierr = PetscTableAdd(gid1_lid1,ival+1,isz+1);CHKERRQ(ierr); 53 isz++; 54 } 55 #else 56 if (ival>Nbs) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"index greater than mat-dim"); 57 if(!PetscBTLookupSet(table,ival)) { nidx[isz++] = ival;} 58 #endif 59 } 60 ierr = ISRestoreIndices(is_in[i],&idx);CHKERRQ(ierr); 61 #if defined (PETSC_USE_CTABLE) 62 ierr = PetscMalloc((isz+1)*sizeof(int),&nidx);CHKERRQ(ierr); 63 ierr = PetscTableGetHeadPosition(gid1_lid1,&tpos);CHKERRQ(ierr); 64 j = 0; 65 while (tpos) { 66 ierr = PetscTableGetNext(gid1_lid1,&tpos,&gid1,&tt);CHKERRQ(ierr); 67 if (tt-- > isz) { SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"index greater than array-dim"); } 68 nidx[tt] = gid1 - 1; 69 j++; 70 } 71 if (j != isz) { SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"table error: jj != isz"); } 72 ierr = ISCreateGeneral(PETSC_COMM_SELF,isz,nidx,(is_out+i));CHKERRQ(ierr); 73 ierr = PetscFree(nidx);CHKERRQ(ierr); 74 #else 75 ierr = ISCreateGeneral(PETSC_COMM_SELF,isz,nidx,(is_out+i));CHKERRQ(ierr); 76 #endif 77 } 78 #if defined (PETSC_USE_CTABLE) 79 ierr = PetscTableDelete(gid1_lid1);CHKERRQ(ierr); 80 #else 81 ierr = PetscBTDestroy(table);CHKERRQ(ierr); 82 ierr = PetscFree(nidx);CHKERRQ(ierr); 83 #endif 84 PetscFunctionReturn(0); 85 } 86 87 #undef __FUNCT__ 88 #define __FUNCT__ "MatExpandIndices_MPISBAIJ" 89 static int MatExpandIndices_MPISBAIJ(Mat C,int imax,const IS is_in[],IS is_out[]) 90 { 91 Mat_MPISBAIJ *baij = (Mat_MPISBAIJ*)C->data; 92 int ierr,bs = baij->bs,n,i,j,k,*idx,*nidx; 93 #if defined (PETSC_USE_CTABLE) 94 int maxsz; 95 #else 96 int Nbs = baij->Nbs; 97 #endif 98 99 PetscFunctionBegin; 100 /* printf(" ... MatExpandIndices_MPISBAIJ is called ...\n"); */ 101 #if defined (PETSC_USE_CTABLE) 102 /* Now check max size */ 103 for (i=0,maxsz=0; i<imax; i++) { 104 ierr = ISGetIndices(is_in[i],&idx);CHKERRQ(ierr); 105 ierr = ISGetLocalSize(is_in[i],&n);CHKERRQ(ierr); 106 if (n*bs > maxsz) maxsz = n*bs; 107 } 108 ierr = PetscMalloc((maxsz+1)*sizeof(int),&nidx);CHKERRQ(ierr); 109 #else 110 ierr = PetscMalloc((Nbs*bs+1)*sizeof(int),&nidx);CHKERRQ(ierr); 111 #endif 112 113 for (i=0; i<imax; i++) { 114 ierr = ISGetIndices(is_in[i],&idx);CHKERRQ(ierr); 115 ierr = ISGetLocalSize(is_in[i],&n);CHKERRQ(ierr); 116 for (j=0; j<n ; ++j){ 117 for (k=0; k<bs; k++) 118 nidx[j*bs+k] = idx[j]*bs+k; 119 } 120 ierr = ISRestoreIndices(is_in[i],&idx);CHKERRQ(ierr); 121 ierr = ISCreateGeneral(PETSC_COMM_SELF,n*bs,nidx,is_out+i);CHKERRQ(ierr); 122 } 123 ierr = PetscFree(nidx);CHKERRQ(ierr); 124 PetscFunctionReturn(0); 125 } 126 127 #undef __FUNCT__ 128 #define __FUNCT__ "MatIncreaseOverlap_MPISBAIJ" 129 int MatIncreaseOverlap_MPISBAIJ(Mat C,int is_max,IS is[],int ov) 130 { 131 int i,ierr; 132 IS *is_new; 133 134 PetscFunctionBegin; 135 ierr = PetscMalloc(is_max*sizeof(IS),&is_new);CHKERRQ(ierr); 136 /* Convert the indices into block format */ 137 ierr = MatCompressIndicesGeneral_MPISBAIJ(C,is_max,is,is_new);CHKERRQ(ierr); 138 if (ov < 0){ SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Negative overlap specified\n");} 139 for (i=0; i<ov; ++i) { 140 ierr = MatIncreaseOverlap_MPISBAIJ_Once(C,is_max,is_new);CHKERRQ(ierr); 141 } 142 for (i=0; i<is_max; i++) {ierr = ISDestroy(is[i]);CHKERRQ(ierr);} 143 ierr = MatExpandIndices_MPISBAIJ(C,is_max,is_new,is);CHKERRQ(ierr); 144 for (i=0; i<is_max; i++) {ierr = ISDestroy(is_new[i]);CHKERRQ(ierr);} 145 ierr = PetscFree(is_new);CHKERRQ(ierr); 146 PetscFunctionReturn(0); 147 } 148 149 typedef enum {MINE,OTHER} WhoseOwner; 150 /* data1, odata1 and odata2 are packed in the format (for communication): 151 data[0] = is_max, no of is 152 data[1] = size of is[0] 153 ... 154 data[is_max] = size of is[is_max-1] 155 data[is_max + 1] = data(is[0]) 156 ... 157 data[is_max+1+sum(size of is[k]), k=0,...,i-1] = data(is[i]) 158 ... 159 data2 is packed in the format (for creating output is[]): 160 data[0] = is_max, no of is 161 data[1] = size of is[0] 162 ... 163 data[is_max] = size of is[is_max-1] 164 data[is_max + 1] = data(is[0]) 165 ... 166 data[is_max + 1 + Mbs*i) = data(is[i]) 167 ... 168 */ 169 #undef __FUNCT__ 170 #define __FUNCT__ "MatIncreaseOverlap_MPISBAIJ_Once" 171 static int MatIncreaseOverlap_MPISBAIJ_Once(Mat C,int is_max,IS is[]) 172 { 173 Mat_MPISBAIJ *c = (Mat_MPISBAIJ*)C->data; 174 int len,*idx_i,isz,col,*n,*data1,*data2,*data2_i,*data,*data_i, 175 size,rank,Mbs,i,j,k,ierr,nrqs,*odata1,*odata2, 176 tag1,tag2,flag,proc_id,**odata2_ptr; 177 char *t_p; 178 MPI_Comm comm; 179 MPI_Request *s_waits,r_req; 180 MPI_Status *s_status,r_status; 181 PetscBT *table; /* mark indices of this processor's is[] */ 182 PetscBT table_i; 183 PetscBT otable; /* mark indices of other processors' is[] */ 184 PetscFunctionBegin; 185 186 comm = C->comm; 187 size = c->size; 188 rank = c->rank; 189 Mbs = c->Mbs; 190 191 ierr = PetscObjectGetNewTag((PetscObject)C,&tag1);CHKERRQ(ierr); 192 ierr = PetscObjectGetNewTag((PetscObject)C,&tag2);CHKERRQ(ierr); 193 194 /* create tables for marking indices */ 195 len = is_max*sizeof(PetscBT) + (Mbs/PETSC_BITS_PER_BYTE+1)*is_max*sizeof(char) + 1; 196 ierr = PetscMalloc(len,&table);CHKERRQ(ierr); 197 t_p = (char *)(table + is_max); 198 for (i=0; i<is_max; i++) { 199 table[i] = t_p + (Mbs/PETSC_BITS_PER_BYTE+1)*i; 200 } 201 ierr = PetscBTCreate(Mbs,otable);CHKERRQ(ierr); 202 203 /* 1. Send this processor's is[] to all other processors */ 204 /*-------------------------------------------------------*/ 205 /* Allocate Memory for outgoing messages */ 206 len = is_max*sizeof(int); 207 ierr = PetscMalloc(len,&n);CHKERRQ(ierr); 208 209 len = 1 + is_max; 210 for (i=0; i<is_max; i++) { 211 ierr = ISGetLocalSize(is[i],&n[i]);CHKERRQ(ierr); 212 len += n[i]; 213 } 214 ierr = PetscMalloc(len*sizeof(int),&data1);CHKERRQ(ierr); 215 216 /* Form the outgoing messages */ 217 data1[0] = is_max; 218 k = is_max + 1; 219 for (i=0; i<is_max; i++) { 220 data1[1+i] = n[i]; 221 ierr = ISGetIndices(is[i],&idx_i);CHKERRQ(ierr); 222 for (j=0; j<data1[i+1]; j++){ 223 data1[k++] = idx_i[j]; 224 } 225 ierr = ISRestoreIndices(is[i],&idx_i);CHKERRQ(ierr); 226 ierr = ISDestroy(is[i]);CHKERRQ(ierr); 227 } 228 if (k != len) SETERRQ2(1,"Error on forming the outgoing messages: k %d != len %d",k,len); 229 ierr = PetscFree(n);CHKERRQ(ierr); 230 231 /* Now post the sends */ 232 ierr = PetscMalloc(size*sizeof(MPI_Request),&s_waits);CHKERRQ(ierr); 233 k = 0; 234 for (proc_id=0; proc_id<size; ++proc_id) { /* send data1 to processor [proc_id] */ 235 if (proc_id != rank){ 236 ierr = MPI_Isend(data1,len,MPI_INT,proc_id,tag1,comm,s_waits+k);CHKERRQ(ierr); 237 /* printf(" [%d] send %d msg to [%d], data1: \n",rank,len,proc_id); */ 238 k++; 239 } 240 } 241 242 /* 2. Do local work on this processor's is[] */ 243 /*-------------------------------------------*/ 244 ierr = MatIncreaseOverlap_MPISBAIJ_Local(C,data1,MINE,&data,table);CHKERRQ(ierr); 245 246 /* 3. Receive other's is[] and process. Then send back */ 247 /*-----------------------------------------------------*/ 248 /* Sending this processor's is[] is done */ 249 nrqs = size-1; 250 ierr = PetscMalloc(size*sizeof(MPI_Status),&s_status);CHKERRQ(ierr); 251 ierr = MPI_Waitall(nrqs,s_waits,s_status);CHKERRQ(ierr); 252 253 ierr = PetscMalloc(size*sizeof(int**),&odata2_ptr);CHKERRQ(ierr); 254 k = 0; 255 do { 256 /* Receive messages */ 257 ierr = MPI_Iprobe(MPI_ANY_SOURCE,tag1,comm,&flag,&r_status);CHKERRQ(ierr); 258 if (flag){ 259 ierr = MPI_Get_count(&r_status,MPI_INT,&len);CHKERRQ(ierr); 260 proc_id = r_status.MPI_SOURCE; 261 ierr = PetscMalloc(len*sizeof(int),&odata1);CHKERRQ(ierr); 262 ierr = MPI_Irecv(odata1,len,MPI_INT,proc_id,r_status.MPI_TAG,comm,&r_req);CHKERRQ(ierr); 263 /* printf(" [%d] recv %d msg from [%d]\n",rank,len,proc_id); */ 264 265 /* Process messages */ 266 ierr = MatIncreaseOverlap_MPISBAIJ_Local(C,odata1,OTHER,&odata2_ptr[k],&otable);CHKERRQ(ierr); 267 odata2 = odata2_ptr[k]; 268 len = 1 + odata2[0]; 269 for (i=0; i<odata2[0]; i++){ 270 len += odata2[1 + i]; 271 } 272 273 /* Send messages back */ 274 ierr = MPI_Isend(odata2,len,MPI_INT,proc_id,tag2,comm,s_waits+k);CHKERRQ(ierr); 275 /* printf(" [%d] send %d msg back to [%d] \n",rank,len,proc_id); */ 276 277 ierr = PetscFree(odata1);CHKERRQ(ierr); 278 k++; 279 } 280 } while (k < nrqs); 281 282 /* 4. Receive work done on other processors, then merge */ 283 /*--------------------------------------------------------*/ 284 /* Allocate memory for incoming data */ 285 len = (1+is_max*(Mbs+1)); 286 ierr = PetscMalloc(len*sizeof(int),&data2);CHKERRQ(ierr); 287 288 /* Sending others' is[] is done */ 289 ierr = MPI_Waitall(nrqs,s_waits,s_status);CHKERRQ(ierr); 290 for (k=0; k<nrqs; k++){ 291 ierr = PetscFree(odata2_ptr[k]);CHKERRQ(ierr); 292 } 293 ierr = PetscFree(odata2_ptr);CHKERRQ(ierr); 294 295 k = 0; 296 do { 297 /* Receive messages */ 298 ierr = MPI_Iprobe(MPI_ANY_SOURCE,tag2,comm,&flag,&r_status); 299 if (flag){ 300 ierr = MPI_Get_count(&r_status,MPI_INT,&len);CHKERRQ(ierr); 301 proc_id = r_status.MPI_SOURCE; 302 ierr = MPI_Irecv(data2,len,MPI_INT,proc_id,r_status.MPI_TAG,comm,&r_req);CHKERRQ(ierr); 303 /* printf(" [%d] recv %d msg from [%d], data2:\n",rank,len,proc_id); */ 304 305 /* Add data2 into data */ 306 data2_i = data2 + 1 + is_max; 307 for (i=0; i<is_max; i++){ 308 table_i = table[i]; 309 data_i = data + 1 + is_max + Mbs*i; 310 isz = data[1+i]; 311 for (j=0; j<data2[1+i]; j++){ 312 col = data2_i[j]; 313 if (!PetscBTLookupSet(table_i,col)) {data_i[isz++] = col;} 314 } 315 data[1+i] = isz; 316 if (i < is_max - 1) data2_i += data2[1+i]; 317 } 318 k++; 319 } 320 } while (k < nrqs); 321 322 /* 5. Create new is[] */ 323 /*--------------------*/ 324 for (i=0; i<is_max; i++) { 325 data_i = data + 1 + is_max + Mbs*i; 326 ierr = ISCreateGeneral(PETSC_COMM_SELF,data[1+i],data_i,is+i);CHKERRQ(ierr); 327 } 328 ierr = PetscFree(data1);CHKERRQ(ierr); 329 ierr = PetscFree(data2);CHKERRQ(ierr); 330 ierr = PetscFree(data);CHKERRQ(ierr); 331 ierr = PetscFree(s_waits);CHKERRQ(ierr); 332 ierr = PetscFree(s_status);CHKERRQ(ierr); 333 ierr = PetscFree(table);CHKERRQ(ierr); 334 ierr = PetscBTDestroy(otable);CHKERRQ(ierr); 335 PetscFunctionReturn(0); 336 } 337 338 #undef __FUNCT__ 339 #define __FUNCT__ "MatIncreaseOverlap_MPISBAIJ_Local" 340 /* 341 MatIncreaseOverlap_MPISBAIJ_Local - Called by MatIncreaseOverlap, to do 342 the work on the local processor. 343 344 Inputs: 345 C - MAT_MPISBAIJ; 346 data - holds is[]. See MatIncreaseOverlap_MPISBAIJ_Once() for the format. 347 whose - whose is[] to be processed, 348 MINE: this processor's is[] 349 OTHER: other processor's is[] 350 Output: 351 data_new - whose = MINE: 352 holds input and newly found indices in the same format as data 353 whose = OTHER: 354 only holds the newly found indices 355 table - table[i]: mark the indices of is[i], i=0,...,is_max. Used only in the case 'whose=MINE'. 356 */ 357 static int MatIncreaseOverlap_MPISBAIJ_Local(Mat C,int *data,int whose,int **data_new,PetscBT *table) 358 { 359 Mat_MPISBAIJ *c = (Mat_MPISBAIJ*)C->data; 360 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)(c->A)->data; 361 Mat_SeqBAIJ *b = (Mat_SeqBAIJ*)(c->B)->data; 362 int ierr,row,mbs,Mbs,*nidx,*nidx_i,col,isz,isz0,*ai,*aj,*bi,*bj,*garray,rstart,l; 363 int a_start,a_end,b_start,b_end,i,j,k,is_max,*idx_i,n; 364 PetscBT table0; /* mark the indices of input is[] for look up */ 365 PetscBT table_i; /* poits to i-th table. When whose=OTHER, a single table is used for all is[] */ 366 367 PetscFunctionBegin; 368 Mbs = c->Mbs; mbs = a->mbs; 369 ai = a->i; aj = a->j; 370 bi = b->i; bj = b->j; 371 garray = c->garray; 372 rstart = c->rstart; 373 is_max = data[0]; 374 375 ierr = PetscBTCreate(Mbs,table0);CHKERRQ(ierr); 376 377 ierr = PetscMalloc((1+is_max*(Mbs+1))*sizeof(int),&nidx);CHKERRQ(ierr); 378 nidx[0] = is_max; 379 380 idx_i = data + is_max + 1; /* ptr to input is[0] array */ 381 nidx_i = nidx + is_max + 1; /* ptr to output is[0] array */ 382 for (i=0; i<is_max; i++) { /* for each is */ 383 isz = 0; 384 n = data[1+i]; /* size of input is[i] */ 385 386 /* initialize table_i, set table0 */ 387 if (whose == MINE){ /* process this processor's is[] */ 388 table_i = table[i]; 389 nidx_i = nidx + 1+ is_max + Mbs*i; 390 } else { /* process other processor's is[] - only use one temp table */ 391 table_i = *table; 392 } 393 ierr = PetscBTMemzero(Mbs,table_i);CHKERRQ(ierr); 394 ierr = PetscBTMemzero(Mbs,table0);CHKERRQ(ierr); 395 if (n > 0) { 396 isz0 = 0; 397 for (j=0; j<n; j++){ 398 col = idx_i[j]; 399 if (col >= Mbs) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"index col %d >= Mbs %d",col,Mbs); 400 if(!PetscBTLookupSet(table_i,col)) { 401 ierr = PetscBTSet(table0,col);CHKERRQ(ierr); 402 if (whose == MINE) {nidx_i[isz0] = col;} 403 isz0++; 404 } 405 } 406 407 if (whose == MINE) {isz = isz0;} 408 k = 0; /* no. of indices from input is[i] that have been examined */ 409 for (row=0; row<mbs; row++){ 410 a_start = ai[row]; a_end = ai[row+1]; 411 b_start = bi[row]; b_end = bi[row+1]; 412 if (PetscBTLookup(table0,row+rstart)){ /* row is on input is[i]: 413 do row search: collect all col in this row */ 414 for (l = a_start; l<a_end ; l++){ /* Amat */ 415 col = aj[l] + rstart; 416 if (!PetscBTLookupSet(table_i,col)) {nidx_i[isz++] = col;} 417 } 418 for (l = b_start; l<b_end ; l++){ /* Bmat */ 419 col = garray[bj[l]]; 420 if (!PetscBTLookupSet(table_i,col)) {nidx_i[isz++] = col;} 421 } 422 k++; 423 if (k >= isz0) break; /* for (row=0; row<mbs; row++) */ 424 } else { /* row is not on input is[i]: 425 do col serach: add row onto nidx_i if there is a col in nidx_i */ 426 for (l = a_start; l<a_end ; l++){ /* Amat */ 427 col = aj[l] + rstart; 428 if (PetscBTLookup(table0,col)){ 429 if (!PetscBTLookupSet(table_i,row+rstart)) {nidx_i[isz++] = row+rstart;} 430 break; /* for l = start; l<end ; l++) */ 431 } 432 } 433 for (l = b_start; l<b_end ; l++){ /* Bmat */ 434 col = garray[bj[l]]; 435 if (PetscBTLookup(table0,col)){ 436 if (!PetscBTLookupSet(table_i,row+rstart)) {nidx_i[isz++] = row+rstart;} 437 break; /* for l = start; l<end ; l++) */ 438 } 439 } 440 } 441 } 442 } /* if (n > 0) */ 443 444 if (i < is_max - 1){ 445 idx_i += n; /* ptr to input is[i+1] array */ 446 nidx_i += isz; /* ptr to output is[i+1] array */ 447 } 448 nidx[1+i] = isz; /* size of new is[i] */ 449 } /* for each is */ 450 *data_new = nidx; 451 ierr = PetscBTDestroy(table0);CHKERRQ(ierr); 452 453 PetscFunctionReturn(0); 454 } 455 456 457