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 Mat_MPISBAIJ *c = (Mat_MPISBAIJ*)C->data; 132 int i,ierr,size; 133 IS *is_new; 134 135 PetscFunctionBegin; 136 ierr = PetscMalloc(is_max*sizeof(IS),&is_new);CHKERRQ(ierr); 137 /* Convert the indices into block format */ 138 ierr = MatCompressIndicesGeneral_MPISBAIJ(C,is_max,is,is_new);CHKERRQ(ierr); 139 if (ov < 0){ SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Negative overlap specified\n");} 140 for (i=0; i<ov; ++i) { 141 ierr = MatIncreaseOverlap_MPISBAIJ_Once(C,is_max,is_new);CHKERRQ(ierr); 142 } 143 for (i=0; i<is_max; i++) {ierr = ISDestroy(is[i]);CHKERRQ(ierr);} 144 ierr = MatExpandIndices_MPISBAIJ(C,is_max,is_new,is);CHKERRQ(ierr); 145 for (i=0; i<is_max; i++) {ierr = ISDestroy(is_new[i]);CHKERRQ(ierr);} 146 ierr = PetscFree(is_new);CHKERRQ(ierr); 147 PetscFunctionReturn(0); 148 } 149 150 typedef enum {MINE,OTHER} whose; 151 /* data1, odata1 and odata2 are packed in the format (for communication): 152 data[0] = is_max, no of is 153 data[1] = size of is[0] 154 ... 155 data[is_max] = size of is[is_max-1] 156 data[is_max + 1] = data(is[0]) 157 ... 158 data[is_max+1+sum(size of is[k]), k=0,...,i-1] = data(is[i]) 159 ... 160 data2 is packed in the format (for creating output is[]): 161 data[0] = is_max, no of is 162 data[1] = size of is[0] 163 ... 164 data[is_max] = size of is[is_max-1] 165 data[is_max + 1] = data(is[0]) 166 ... 167 data[is_max + 1 + Mbs*i) = data(is[i]) 168 ... 169 */ 170 #undef __FUNCT__ 171 #define __FUNCT__ "MatIncreaseOverlap_MPISBAIJ_Once" 172 static int MatIncreaseOverlap_MPISBAIJ_Once(Mat C,int is_max,IS is[]) 173 { 174 Mat_MPISBAIJ *c = (Mat_MPISBAIJ*)C->data; 175 int len,*idx_i,isz,col,*n,*data1,*data2,*data2_i,*data,*data_i, 176 size,rank,Mbs,i,j,k,ierr,nrqs,*odata1,*odata2, 177 tag1,tag2,flag,proc_id,**odata2_ptr; 178 char *t_p; 179 MPI_Comm comm; 180 MPI_Request *s_waits,r_req; 181 MPI_Status *s_status,r_status; 182 PetscBT *table; /* mark indices of this processor's is[] */ 183 PetscBT table_i; 184 PetscBT otable; /* mark indices of other processors' is[] */ 185 PetscFunctionBegin; 186 187 comm = C->comm; 188 size = c->size; 189 rank = c->rank; 190 Mbs = c->Mbs; 191 192 ierr = PetscObjectGetNewTag((PetscObject)C,&tag1);CHKERRQ(ierr); 193 ierr = PetscObjectGetNewTag((PetscObject)C,&tag2);CHKERRQ(ierr); 194 195 /* create tables for marking indices */ 196 len = is_max*sizeof(PetscBT) + (Mbs/PETSC_BITS_PER_BYTE+1)*is_max*sizeof(char) + 1; 197 ierr = PetscMalloc(len,&table);CHKERRQ(ierr); 198 t_p = (char *)(table + is_max); 199 for (i=0; i<is_max; i++) { 200 table[i] = t_p + (Mbs/PETSC_BITS_PER_BYTE+1)*i; 201 } 202 ierr = PetscBTCreate(Mbs,otable);CHKERRQ(ierr); 203 204 /* 1. Send this processor's is[] to all other processors */ 205 /*-------------------------------------------------------*/ 206 /* Allocate Memory for outgoing messages */ 207 len = is_max*sizeof(int); 208 ierr = PetscMalloc(len,&n);CHKERRQ(ierr); 209 210 len = 1 + is_max; 211 for (i=0; i<is_max; i++) { 212 ierr = ISGetLocalSize(is[i],&n[i]);CHKERRQ(ierr); 213 len += n[i]; 214 } 215 ierr = PetscMalloc(len*sizeof(int),&data1);CHKERRQ(ierr); 216 217 /* Form the outgoing messages */ 218 data1[0] = is_max; 219 k = is_max + 1; 220 for (i=0; i<is_max; i++) { 221 data1[1+i] = n[i]; 222 ierr = ISGetIndices(is[i],&idx_i);CHKERRQ(ierr); 223 for (j=0; j<data1[i+1]; j++){ 224 data1[k++] = idx_i[j]; 225 } 226 ierr = ISRestoreIndices(is[i],&idx_i);CHKERRQ(ierr); 227 ierr = ISDestroy(is[i]);CHKERRQ(ierr); 228 } 229 if (k != len) SETERRQ2(1,"Error on forming the outgoing messages: k %d != len %d",k,len); 230 ierr = PetscFree(n);CHKERRQ(ierr); 231 232 /* Now post the sends */ 233 ierr = PetscMalloc(size*sizeof(MPI_Request),&s_waits);CHKERRQ(ierr); 234 k = 0; 235 for (proc_id=0; proc_id<size; ++proc_id) { /* send data1 to processor [proc_id] */ 236 if (proc_id != rank){ 237 ierr = MPI_Isend(data1,len,MPI_INT,proc_id,tag1,comm,s_waits+k);CHKERRQ(ierr); 238 /* printf(" [%d] send %d msg to [%d], data1: \n",rank,len,proc_id); */ 239 k++; 240 } 241 } 242 243 /* 2. Do local work on this processor's is[] */ 244 /*-------------------------------------------*/ 245 ierr = MatIncreaseOverlap_MPISBAIJ_Local(C,data1,MINE,&data,table);CHKERRQ(ierr); 246 247 /* 3. Receive other's is[] and process. Then send back */ 248 /*-----------------------------------------------------*/ 249 /* Sending this processor's is[] is done */ 250 nrqs = size-1; 251 ierr = PetscMalloc(size*sizeof(MPI_Status),&s_status);CHKERRQ(ierr); 252 ierr = MPI_Waitall(nrqs,s_waits,s_status);CHKERRQ(ierr); 253 254 ierr = PetscMalloc(size*sizeof(int**),&odata2_ptr);CHKERRQ(ierr); 255 k = 0; 256 do { 257 /* Receive messages */ 258 ierr = MPI_Iprobe(MPI_ANY_SOURCE,tag1,comm,&flag,&r_status);CHKERRQ(ierr); 259 if (flag){ 260 ierr = MPI_Get_count(&r_status,MPI_INT,&len);CHKERRQ(ierr); 261 proc_id = r_status.MPI_SOURCE; 262 ierr = PetscMalloc(len*sizeof(int),&odata1);CHKERRQ(ierr); 263 ierr = MPI_Irecv(odata1,len,MPI_INT,proc_id,r_status.MPI_TAG,comm,&r_req);CHKERRQ(ierr); 264 /* printf(" [%d] recv %d msg from [%d]\n",rank,len,proc_id); */ 265 266 /* Process messages */ 267 ierr = MatIncreaseOverlap_MPISBAIJ_Local(C,odata1,OTHER,&odata2_ptr[k],&otable);CHKERRQ(ierr); 268 odata2 = odata2_ptr[k]; 269 len = 1 + odata2[0]; 270 for (i=0; i<odata2[0]; i++){ 271 len += odata2[1 + i]; 272 } 273 274 /* Send messages back */ 275 ierr = MPI_Isend(odata2,len,MPI_INT,proc_id,tag2,comm,s_waits+k);CHKERRQ(ierr); 276 /* printf(" [%d] send %d msg back to [%d] \n",rank,len,proc_id); */ 277 278 ierr = PetscFree(odata1);CHKERRQ(ierr); 279 k++; 280 } 281 } while (k < nrqs); 282 283 /* 4. Receive work done on other processors, then merge */ 284 /*--------------------------------------------------------*/ 285 /* Allocate memory for incoming data */ 286 len = (1+is_max*(Mbs+1)); 287 ierr = PetscMalloc(len*sizeof(int),&data2);CHKERRQ(ierr); 288 289 /* Sending others' is[] is done */ 290 ierr = MPI_Waitall(nrqs,s_waits,s_status);CHKERRQ(ierr); 291 for (k=0; k<nrqs; k++){ 292 ierr = PetscFree(odata2_ptr[k]);CHKERRQ(ierr); 293 } 294 ierr = PetscFree(odata2_ptr);CHKERRQ(ierr); 295 296 k = 0; 297 do { 298 /* Receive messages */ 299 ierr = MPI_Iprobe(MPI_ANY_SOURCE,tag2,comm,&flag,&r_status); 300 if (flag){ 301 ierr = MPI_Get_count(&r_status,MPI_INT,&len);CHKERRQ(ierr); 302 proc_id = r_status.MPI_SOURCE; 303 ierr = MPI_Irecv(data2,len,MPI_INT,proc_id,r_status.MPI_TAG,comm,&r_req);CHKERRQ(ierr); 304 /* printf(" [%d] recv %d msg from [%d], data2:\n",rank,len,proc_id); */ 305 306 /* Add data2 into data */ 307 data2_i = data2 + 1 + is_max; 308 for (i=0; i<is_max; i++){ 309 table_i = table[i]; 310 data_i = data + 1 + is_max + Mbs*i; 311 isz = data[1+i]; 312 for (j=0; j<data2[1+i]; j++){ 313 col = data2_i[j]; 314 if (!PetscBTLookupSet(table_i,col)) {data_i[isz++] = col;} 315 } 316 data[1+i] = isz; 317 if (i < is_max - 1) data2_i += data2[1+i]; 318 } 319 k++; 320 } 321 } while (k < nrqs); 322 323 /* 5. Create new is[] */ 324 /*--------------------*/ 325 for (i=0; i<is_max; i++) { 326 data_i = data + 1 + is_max + Mbs*i; 327 ierr = ISCreateGeneral(PETSC_COMM_SELF,data[1+i],data_i,is+i);CHKERRQ(ierr); 328 } 329 ierr = PetscFree(data1);CHKERRQ(ierr); 330 ierr = PetscFree(data2);CHKERRQ(ierr); 331 ierr = PetscFree(data);CHKERRQ(ierr); 332 ierr = PetscFree(s_waits);CHKERRQ(ierr); 333 ierr = PetscFree(s_status);CHKERRQ(ierr); 334 ierr = PetscFree(table);CHKERRQ(ierr); 335 ierr = PetscBTDestroy(otable);CHKERRQ(ierr); 336 PetscFunctionReturn(0); 337 } 338 339 #undef __FUNCT__ 340 #define __FUNCT__ "MatIncreaseOverlap_MPISBAIJ_Local" 341 /* 342 MatIncreaseOverlap_MPISBAIJ_Local - Called by MatIncreaseOverlap, to do 343 the work on the local processor. 344 345 Inputs: 346 C - MAT_MPISBAIJ; 347 data - holds is[]. See MatIncreaseOverlap_MPISBAIJ_Once() for the format. 348 whose - whose is[] to be processed, 349 MINE: this processor's is[] 350 OTHER: other processor's is[] 351 Output: 352 data_new - whose = MINE: 353 holds input and newly found indices in the same format as data 354 whose = OTHER: 355 only holds the newly found indices 356 table - table[i]: mark the indices of is[i], i=0,...,is_max. Used only in the case 'whose=MINE'. 357 */ 358 static int MatIncreaseOverlap_MPISBAIJ_Local(Mat C,int *data,int whose,int **data_new,PetscBT *table) 359 { 360 Mat_MPISBAIJ *c = (Mat_MPISBAIJ*)C->data; 361 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)(c->A)->data; 362 Mat_SeqBAIJ *b = (Mat_SeqBAIJ*)(c->B)->data; 363 int ierr,row,mbs,Mbs,*nidx,*nidx_i,col,isz,isz0,*ai,*aj,*bi,*bj,*garray,rstart,l; 364 int a_start,a_end,b_start,b_end,i,j,k,is_max,*idx_i,n; 365 PetscBT table0; /* mark the indices of input is[] for look up */ 366 PetscBT table_i; /* poits to i-th table. When whose=OTHER, a single table is used for all is[] */ 367 368 PetscFunctionBegin; 369 Mbs = c->Mbs; mbs = a->mbs; 370 ai = a->i; aj = a->j; 371 bi = b->i; bj = b->j; 372 garray = c->garray; 373 rstart = c->rstart; 374 is_max = data[0]; 375 376 ierr = PetscBTCreate(Mbs,table0);CHKERRQ(ierr); 377 378 ierr = PetscMalloc((1+is_max*(Mbs+1))*sizeof(int),&nidx);CHKERRQ(ierr); 379 nidx[0] = is_max; 380 381 idx_i = data + is_max + 1; /* ptr to input is[0] array */ 382 nidx_i = nidx + is_max + 1; /* ptr to output is[0] array */ 383 for (i=0; i<is_max; i++) { /* for each is */ 384 isz = 0; 385 n = data[1+i]; /* size of input is[i] */ 386 387 /* initialize table_i, set table0 */ 388 if (whose == MINE){ /* process this processor's is[] */ 389 table_i = table[i]; 390 nidx_i = nidx + 1+ is_max + Mbs*i; 391 } else { /* process other processor's is[] - only use one temp table */ 392 table_i = *table; 393 } 394 ierr = PetscBTMemzero(Mbs,table_i);CHKERRQ(ierr); 395 ierr = PetscBTMemzero(Mbs,table0);CHKERRQ(ierr); 396 if (n > 0) { 397 isz0 = 0; 398 for (j=0; j<n; j++){ 399 col = idx_i[j]; 400 if (col >= Mbs) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"index col %d >= Mbs %d",col,Mbs); 401 if(!PetscBTLookupSet(table_i,col)) { 402 ierr = PetscBTSet(table0,col);CHKERRQ(ierr); 403 if (whose == MINE) {nidx_i[isz0] = col;} 404 isz0++; 405 } 406 } 407 408 if (whose == MINE) {isz = isz0;} 409 k = 0; /* no. of indices from input is[i] that have been examined */ 410 for (row=0; row<mbs; row++){ 411 a_start = ai[row]; a_end = ai[row+1]; 412 b_start = bi[row]; b_end = bi[row+1]; 413 if (PetscBTLookup(table0,row+rstart)){ /* row is on input is[i]: 414 do row search: collect all col in this row */ 415 for (l = a_start; l<a_end ; l++){ /* Amat */ 416 col = aj[l] + rstart; 417 if (!PetscBTLookupSet(table_i,col)) {nidx_i[isz++] = col;} 418 } 419 for (l = b_start; l<b_end ; l++){ /* Bmat */ 420 col = garray[bj[l]]; 421 if (!PetscBTLookupSet(table_i,col)) {nidx_i[isz++] = col;} 422 } 423 k++; 424 if (k >= isz0) break; /* for (row=0; row<mbs; row++) */ 425 } else { /* row is not on input is[i]: 426 do col serach: add row onto nidx_i if there is a col in nidx_i */ 427 for (l = a_start; l<a_end ; l++){ /* Amat */ 428 col = aj[l] + rstart; 429 if (PetscBTLookup(table0,col)){ 430 if (!PetscBTLookupSet(table_i,row+rstart)) {nidx_i[isz++] = row+rstart;} 431 break; /* for l = start; l<end ; l++) */ 432 } 433 } 434 for (l = b_start; l<b_end ; l++){ /* Bmat */ 435 col = garray[bj[l]]; 436 if (PetscBTLookup(table0,col)){ 437 if (!PetscBTLookupSet(table_i,row+rstart)) {nidx_i[isz++] = row+rstart;} 438 break; /* for l = start; l<end ; l++) */ 439 } 440 } 441 } 442 } 443 } /* if (n > 0) */ 444 445 if (i < is_max - 1){ 446 idx_i += n; /* ptr to input is[i+1] array */ 447 nidx_i += isz; /* ptr to output is[i+1] array */ 448 } 449 nidx[1+i] = isz; /* size of new is[i] */ 450 } /* /* for each is */ 451 *data_new = nidx; 452 ierr = PetscBTDestroy(table0);CHKERRQ(ierr); 453 454 PetscFunctionReturn(0); 455 } 456 457 458