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