1 #ifndef lint 2 static char vcid[] = "$Id: mpiov.c,v 1.10 1996/02/01 00:24:06 balay Exp balay $"; 3 #endif 4 5 #include "mpiaij.h" 6 #include "inline/bitarray.h" 7 8 int MatIncreaseOverlap_MPIAIJ_private(Mat, int, IS *); 9 10 int MatIncreaseOverlap_MPIAIJ(Mat C, int is_max, IS *is, int ov) 11 { 12 int i, ierr; 13 if (ov < 0){ SETERRQ(1," MatIncreaseOverlap_MPIAIJ: negative overlap specified\n");} 14 for (i =0; i<ov; ++i) { 15 ierr = MatIncreaseOverlap_MPIAIJ_private(C, is_max, is); CHKERRQ(ierr); 16 } 17 return 0; 18 } 19 int FindOverlapLocal(Mat , int , char **,int*, int**); 20 int FindOverlapRecievedMesg(Mat , int, int **, int**, int* ); 21 22 /* 23 Sample message format: 24 If a processor A wants processor B to process some elements corresponding 25 to index sets 1s[1], is[5] 26 mesg [0] = 2 ( no of index sets in the mesg) 27 ----------- 28 mesg [1] = 1 => is[1] 29 mesg [2] = sizeof(is[1]); 30 ----------- 31 mesg [5] = 5 => is[5] 32 mesg [6] = sizeof(is[5]); 33 ----------- 34 mesg [7] 35 mesg [n] datas[1] 36 ----------- 37 mesg[n+1] 38 mesg[m] data(is[5]) 39 ----------- 40 */ 41 int MatIncreaseOverlap_MPIAIJ_private(Mat C, int is_max, IS *is) 42 { 43 Mat_MPIAIJ *c = (Mat_MPIAIJ *) C->data; 44 int **idx, *n, *w1, *w2, *w3, *w4, *rtable,**data; 45 int size, rank, m,i,j,k, ierr, **rbuf, row, proc, mct, msz, **outdat, **ptr; 46 int *ctr, *pa, tag, *tmp,bsz, nmsg , *isz, *isz1, **xdata; 47 int bsz1, **rbuf2; 48 char **table; 49 MPI_Comm comm; 50 MPI_Request *send_waits,*recv_waits,*send_waits2,*recv_waits2 ; 51 MPI_Status *send_status ,*recv_status; 52 double space, fr, maxs; 53 54 comm = C->comm; 55 tag = C->tag; 56 size = c->size; 57 rank = c->rank; 58 m = c->M; 59 60 61 TrSpace( &space, &fr, &maxs ); 62 /* MPIU_printf(MPI_COMM_SELF,"[%d] allocated space = %f fragments = %f max ever allocated = %f\n", rank, space, fr, maxs); */ 63 64 idx = (int **)PetscMalloc((is_max+1)*sizeof(int *)); CHKPTRQ(idx); 65 n = (int *)PetscMalloc((is_max+1)*sizeof(int )); CHKPTRQ(n); 66 rtable = (int *)PetscMalloc((m+1)*sizeof(int )); CHKPTRQ(rtable); 67 /* Hash table for maping row ->proc */ 68 69 for ( i=0 ; i<is_max ; i++) { 70 ierr = ISGetIndices(is[i],&idx[i]); CHKERRQ(ierr); 71 ierr = ISGetLocalSize(is[i],&n[i]); CHKERRQ(ierr); 72 } 73 74 /* Create hash table for the mapping :row -> proc*/ 75 for( i=0, j=0; i< size; i++) { 76 for (; j <c->rowners[i+1]; j++) { 77 rtable[j] = i; 78 } 79 } 80 81 /* evaluate communication - mesg to who, length of mesg, and buffer space 82 required. Based on this, buffers are allocated, and data copied into them*/ 83 w1 = (int *)PetscMalloc((size)*4*sizeof(int )); CHKPTRQ(w1); /* mesg size */ 84 w2 = w1 + size; /* if w2[i] marked, then a message to proc i*/ 85 w3 = w2 + size; /* no of IS that needs to be sent to proc i */ 86 w4 = w3 + size; /* temp work space used in determining w1, w2, w3 */ 87 PetscMemzero(w1,(size)*3*sizeof(int)); /* initialise work vector*/ 88 for ( i=0; i<is_max ; i++) { 89 PetscMemzero(w4,(size)*sizeof(int)); /* initialise work vector*/ 90 for ( j =0 ; j < n[i] ; j++) { 91 row = idx[i][j]; 92 proc = rtable[row]; 93 w4[proc]++; 94 } 95 for( j = 0; j < size; j++){ 96 if( w4[j] ) { w1[j] += w4[j]; w3[j] += 1;} 97 } 98 } 99 100 mct = 0; /* no of outgoing messages */ 101 msz = 0; /* total mesg length (for all proc */ 102 w1[rank] = 0; /* no mesg sent to intself */ 103 w3[rank] = 0; 104 for (i =0; i < size ; i++) { 105 if (w1[i]) { w2[i] = 1; mct++;} /* there exists a message to proc i */ 106 } 107 pa = (int *)PetscMalloc((mct +1)*sizeof(int)); CHKPTRQ(pa); /* (proc -array) */ 108 for (i =0, j=0; i < size ; i++) { 109 if (w1[i]) { pa[j] = i; j++; } 110 } 111 112 /* Each message would have a header = 1 + 2*(no of IS) + data */ 113 for (i = 0; i<mct ; i++) { 114 j = pa[i]; 115 w1[j] += w2[j] + 2* w3[j]; 116 msz += w1[j]; 117 } 118 119 /* Allocate Memory for outgoing messages */ 120 outdat = (int **)PetscMalloc( 2*size*sizeof(int*)); CHKPTRQ(outdat); 121 PetscMemzero(outdat, 2*size*sizeof(int*)); 122 tmp = (int *)PetscMalloc((msz+1) *sizeof (int)); CHKPTRQ(tmp); /*mrsg arr */ 123 ptr = outdat +size; /* Pointers to the data in outgoing buffers */ 124 ctr = (int *)PetscMalloc( size*sizeof(int)); CHKPTRQ(ctr); 125 126 { 127 int *iptr = tmp; 128 int ict = 0; 129 for (i = 0; i < mct ; i++) { 130 j = pa[i]; 131 iptr += ict; 132 outdat[j] = iptr; 133 ict = w1[j]; 134 } 135 } 136 137 /* Form the outgoing messages */ 138 /*plug in the headers*/ 139 for ( i=0 ; i<mct ; i++) { 140 j = pa[i]; 141 outdat[j][0] = 0; 142 PetscMemzero(outdat[j]+1, 2 * w3[j]*sizeof(int)); 143 ptr[j] = outdat[j] + 2*w3[j] +1; 144 } 145 146 /* Memory for doing local proc's work*/ 147 table = (char **)PetscMalloc((is_max+1)*sizeof(int *)); CHKPTRQ(table); 148 data = (int **)PetscMalloc((is_max+1)*sizeof(int *)); CHKPTRQ(data); 149 table[0] = (char *)PetscMalloc((m/BITSPERBYTE +1)*(is_max)); CHKPTRQ(table[0]); 150 data [0] = (int *)PetscMalloc((m+1)*(is_max)*sizeof(int)); CHKPTRQ(data[0]); 151 152 for(i = 1; i<is_max ; i++) { 153 table[i] = table[0] + (m/BITSPERBYTE+1)*i; 154 data[i] = data[0] + (m+1)*i; 155 } 156 157 PetscMemzero((void*)*table,(m/BITSPERBYTE+1)*(is_max)); 158 isz = (int *)PetscMalloc((is_max+1) *sizeof(int)); CHKPTRQ(isz); 159 PetscMemzero((void *)isz,(is_max+1) *sizeof(int)); 160 161 /* Parse the IS and update local tables and the outgoing buf with the data*/ 162 for ( i=0 ; i<is_max ; i++) { 163 PetscMemzero(ctr,size*sizeof(int)); 164 for( j=0; j<n[i]; j++) { /* parse the indices of each IS */ 165 row = idx[i][j]; 166 proc = rtable[row]; 167 if (proc != rank) { /* copy to the outgoing buf*/ 168 ctr[proc]++; 169 *ptr[proc] = row; 170 ptr[proc]++; 171 } 172 else { /* Update the table */ 173 if ( !BT_LOOKUP(table[i],row)) { data[i][isz[i]++] = row;} 174 } 175 } 176 /* Update the headers for the current IS */ 177 for( j = 0; j<size; j++) { /* Can Optimise this loop too */ 178 if (ctr[j]) { 179 k= ++outdat[j][0]; 180 outdat[j][2*k] = ctr[j]; 181 outdat[j][2*k-1] = i; 182 } 183 } 184 } 185 186 /* I nolonger need the original indices*/ 187 for( i=0; i< is_max; ++i) { 188 ierr = ISRestoreIndices(is[i], idx+i); CHKERRQ(ierr); 189 } 190 PetscFree(idx); 191 PetscFree(n); 192 PetscFree(rtable); 193 for( i=0; i< is_max; ++i) { 194 ierr = ISDestroy(is[i]); CHKERRQ(ierr); 195 } 196 197 198 /* Do a global reduction to determine how many messages to expect*/ 199 { 200 int *rw1, *rw2; 201 rw1 = (int *)PetscMalloc(2*size*sizeof(int)); CHKPTRQ(rw1); 202 rw2 = rw1+size; 203 MPI_Allreduce((void *)w1, rw1, size, MPI_INT, MPI_MAX, comm); 204 bsz = rw1[rank]; 205 MPI_Allreduce((void *)w2, rw2, size, MPI_INT, MPI_SUM, comm); 206 nmsg = rw2[rank]; 207 PetscFree(rw1); 208 } 209 210 /* Allocate memory for recv buffers . Prob none if nmsg = 0 ???? */ 211 rbuf = (int**) PetscMalloc((nmsg+1) *sizeof(int*)); CHKPTRQ(rbuf); 212 rbuf[0] = (int *) PetscMalloc((nmsg *bsz+1) * sizeof(int)); CHKPTRQ(rbuf[0]); 213 for (i=1; i<nmsg ; ++i) rbuf[i] = rbuf[i-1] + bsz; 214 215 /* Now post the receives */ 216 recv_waits = (MPI_Request *) PetscMalloc((nmsg+1)*sizeof(MPI_Request)); 217 CHKPTRQ(recv_waits); 218 for ( i=0; i<nmsg; ++i){ 219 MPI_Irecv((void *)rbuf[i], bsz, MPI_INT, MPI_ANY_SOURCE, tag, comm, recv_waits+i); 220 } 221 222 /* Now post the sends */ 223 send_waits = (MPI_Request *) PetscMalloc((mct+1)*sizeof(MPI_Request)); 224 CHKPTRQ(send_waits); 225 for( i =0; i< mct; ++i){ 226 j = pa[i]; 227 MPI_Isend( (void *)outdat[j], w1[j], MPI_INT, j, tag, comm, send_waits+i); 228 } 229 230 /* Do Local work*/ 231 ierr = FindOverlapLocal(C, is_max, table,isz, data); CHKERRQ(ierr); 232 /* Extract the matrices */ 233 234 /* Receive messages*/ 235 { 236 int index; 237 238 recv_status = (MPI_Status *) PetscMalloc( (nmsg+1)*sizeof(MPI_Status) ); 239 CHKPTRQ(recv_status); 240 for ( i=0; i< nmsg; ++i) { 241 MPI_Waitany(nmsg, recv_waits, &index, recv_status+i); 242 } 243 244 send_status = (MPI_Status *) PetscMalloc( (mct+1)*sizeof(MPI_Status) ); 245 CHKPTRQ(send_status); 246 MPI_Waitall(mct,send_waits,send_status); 247 } 248 /* Pahse 1 sends are complete - deallocate buffers */ 249 PetscFree(outdat); 250 PetscFree(w1); 251 PetscFree(tmp); 252 253 /* int FindOverlapRecievedMesg(Mat C, int is_max, int *isz, char **table, int **data)*/ 254 xdata = (int **)PetscMalloc((nmsg+1)*sizeof(int *)); CHKPTRQ(xdata); 255 isz1 = (int *)PetscMalloc((nmsg+1) *sizeof(int)); CHKPTRQ(isz1); 256 ierr = FindOverlapRecievedMesg(C, nmsg, rbuf,xdata,isz1); CHKERRQ(ierr); 257 258 /* Nolonger need rbuf. */ 259 PetscFree(rbuf[0]); 260 PetscFree(rbuf); 261 262 263 /* Send the data back*/ 264 /* Do a global reduction to know the buffer space req for incoming messages*/ 265 { 266 int *rw1, *rw2; 267 268 rw1 = (int *)PetscMalloc(2*size*sizeof(int)); CHKPTRQ(rw1); 269 PetscMemzero((void*)rw1,2*size*sizeof(int)); 270 rw2 = rw1+size; 271 for (i =0; i < nmsg ; ++i) { 272 proc = recv_status[i].MPI_SOURCE; 273 rw1[proc] = isz1[i]; 274 } 275 276 MPI_Allreduce((void *)rw1, (void *)rw2, size, MPI_INT, MPI_MAX, comm); 277 bsz1 = rw2[rank]; 278 PetscFree(rw1); 279 } 280 281 /* Allocate buffers*/ 282 283 /* Allocate memory for recv buffers . Prob none if nmsg = 0 ???? */ 284 rbuf2 = (int**) PetscMalloc((mct+1) *sizeof(int*)); CHKPTRQ(rbuf2); 285 rbuf2[0] = (int *) PetscMalloc((mct*bsz1+1) * sizeof(int)); CHKPTRQ(rbuf2[0]); 286 for (i=1; i<mct ; ++i) rbuf2[i] = rbuf2[i-1] + bsz1; 287 288 /* Now post the receives */ 289 recv_waits2 = (MPI_Request *)PetscMalloc((mct+1)*sizeof(MPI_Request)); CHKPTRQ(recv_waits2) 290 CHKPTRQ(recv_waits2); 291 for ( i=0; i<mct; ++i){ 292 MPI_Irecv((void *)rbuf2[i], bsz1, MPI_INT, MPI_ANY_SOURCE, tag, comm, recv_waits2+i); 293 } 294 295 /* Now post the sends */ 296 send_waits2 = (MPI_Request *) PetscMalloc((nmsg+1)*sizeof(MPI_Request)); 297 CHKPTRQ(send_waits2); 298 for( i =0; i< nmsg; ++i){ 299 j = recv_status[i].MPI_SOURCE; 300 MPI_Isend( (void *)xdata[i], isz1[i], MPI_INT, j, tag, comm, send_waits2+i); 301 } 302 303 /* recieve work done on other processors*/ 304 { 305 int index, is_no, ct1, max; 306 MPI_Status *send_status2 ,*recv_status2; 307 308 recv_status2 = (MPI_Status *) PetscMalloc( (mct+1)*sizeof(MPI_Status) ); 309 CHKPTRQ(recv_status2); 310 311 312 for ( i=0; i< mct; ++i) { 313 MPI_Waitany(mct, recv_waits2, &index, recv_status2+i); 314 /* Process the message*/ 315 ct1 = 2*rbuf2[index][0]+1; 316 for (j=1; j<=rbuf2[index][0]; j++) { 317 max = rbuf2[index][2*j]; 318 is_no = rbuf2[index][2*j-1]; 319 for (k=0; k < max ; k++, ct1++) { 320 row = rbuf2[index][ct1]; 321 if(!BT_LOOKUP(table[is_no],row)) { data[is_no][isz[is_no]++] = row;} 322 } 323 } 324 } 325 326 327 send_status2 = (MPI_Status *) PetscMalloc( (nmsg+1)*sizeof(MPI_Status) ); 328 CHKPTRQ(send_status2); 329 MPI_Waitall(nmsg,send_waits2,send_status2); 330 331 PetscFree(send_status2); PetscFree(recv_status2); 332 } 333 334 TrSpace( &space, &fr, &maxs ); 335 /* MPIU_printf(MPI_COMM_SELF,"[%d] allocated space = %f fragments = %f max ever allocated = %f\n", rank, space, fr, maxs);*/ 336 337 PetscFree(ctr); 338 PetscFree(pa); 339 PetscFree(rbuf2[0]); 340 PetscFree(rbuf2); 341 PetscFree(send_waits); 342 PetscFree(recv_waits); 343 PetscFree(send_waits2); 344 PetscFree(recv_waits2); 345 PetscFree(table[0]); 346 PetscFree(table); 347 PetscFree(send_status); 348 PetscFree(recv_status); 349 PetscFree(isz1); 350 PetscFree(xdata[0]); 351 PetscFree(xdata); 352 353 for ( i=0; i<is_max; ++i) { 354 ierr = ISCreateSeq(MPI_COMM_SELF, isz[i], data[i], is+i); CHKERRQ(ierr); 355 } 356 PetscFree(isz); 357 PetscFree(data[0]); 358 PetscFree(data); 359 360 return 0; 361 } 362 363 /* FindOverlapLocal() - Called by MatincreaseOverlap, to do the work on 364 the local processor. 365 366 Inputs: 367 C - MAT_MPIAIJ; 368 is_max - total no of index sets processed at a time; 369 table - an array of char - size = m bits. 370 371 Output: 372 isz - array containing the count of the solution elements correspondign 373 to each index set; 374 data - pointer to the solutions 375 */ 376 int FindOverlapLocal(Mat C, int is_max, char **table, int *isz,int **data) 377 { 378 Mat_MPIAIJ *c = (Mat_MPIAIJ *) C->data; 379 Mat A = c->A, B = c->B; 380 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*b = (Mat_SeqAIJ*)B->data; 381 int start, end, val, max, rstart,cstart, ashift, bshift,*ai, *aj; 382 int *bi, *bj, *garray, i, j, k, row; 383 384 rstart = c->rstart; 385 cstart = c->cstart; 386 ashift = a->indexshift; 387 ai = a->i; 388 aj = a->j +ashift; 389 bshift = b->indexshift; 390 bi = b->i; 391 bj = b->j +bshift; 392 garray = c->garray; 393 394 395 for( i=0; i<is_max; i++) { 396 for ( j=0, max =isz[i] ; j< max; j++) { 397 row = data[i][j] - rstart; 398 start = ai[row]; 399 end = ai[row+1]; 400 for ( k=start; k < end; k++) { /* Amat */ 401 val = aj[k] + ashift + cstart; 402 if(!BT_LOOKUP(table[i],val)) { data[i][isz[i]++] = val;} 403 } 404 start = bi[row]; 405 end = bi[row+1]; 406 for ( k=start; k < end; k++) { /* Bmat */ 407 val = garray[bj[k]+bshift] ; 408 if(! BT_LOOKUP(table[i],val)) { data[i][isz[i]++] = val;} 409 } 410 } 411 } 412 413 return 0; 414 } 415 /* FindOverlapRecievedMesg - Process the recieved messages, 416 and return the output 417 418 Input: 419 C - the matrix 420 nmsg - no of messages being processed. 421 rbuf - an array of pointers to the recieved requests 422 423 Output: 424 xdata - array of messages to be sent back 425 isz1 - size of each message 426 */ 427 int FindOverlapRecievedMesg(Mat C, int nmsg, int ** rbuf, int ** xdata, int * isz1 ) 428 { 429 Mat_MPIAIJ *c = (Mat_MPIAIJ *) C->data; 430 Mat A = c->A, B = c->B; 431 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*b = (Mat_SeqAIJ*)B->data; 432 int rstart,cstart, ashift, bshift,*ai, *aj, *bi, *bj, *garray, i, j, k; 433 int row,total_sz,ct, ct1, ct2, ct3,mem_estimate, oct2, l, start, end; 434 int val, max1, max2, rank, m, no_malloc =0, *tmp, new_estimate, ctr; 435 char *xtable; 436 437 rank = c->rank; 438 m = c->M; 439 rstart = c->rstart; 440 cstart = c->cstart; 441 ashift = a->indexshift; 442 ai = a->i; 443 aj = a->j +ashift; 444 bshift = b->indexshift; 445 bi = b->i; 446 bj = b->j +bshift; 447 garray = c->garray; 448 449 450 for (i =0, ct =0, total_sz =0; i< nmsg; ++i){ 451 ct+= rbuf[i][0]; 452 for ( j = 1; j <= rbuf[i][0] ; j++ ) { total_sz += rbuf[i][2*j]; } 453 } 454 455 max1 = ct*(a->nz +b->nz)/c->m; 456 mem_estimate = 3*((total_sz > max1?total_sz:max1)+1); 457 xdata[0] = (int *)PetscMalloc(mem_estimate *sizeof(int)); CHKPTRQ(xdata[0]); 458 ++no_malloc; 459 xtable = (char *)PetscMalloc((m/BITSPERBYTE+1)); CHKPTRQ(xtable); 460 PetscMemzero((void *)isz1,(nmsg+1) *sizeof(int)); 461 462 ct3 = 0; 463 for (i =0; i< nmsg; i++) { /* for easch mesg from proc i */ 464 ct1 = 2*rbuf[i][0]+1; 465 ct2 = ct1; 466 ct3+= ct1; 467 for (j = 1, max1= rbuf[i][0]; j<=max1; j++) { /* for each IS from proc i*/ 468 PetscMemzero((void *)xtable,(m/BITSPERBYTE+1)); 469 oct2 = ct2; 470 for (k =0; k < rbuf[i][2*j]; k++, ct1++) { 471 row = rbuf[i][ct1]; 472 if(!BT_LOOKUP(xtable,row)) { 473 if (!(ct3 < mem_estimate)) { 474 new_estimate = 1.5*mem_estimate; 475 tmp = (int*) PetscMalloc(new_estimate * sizeof(int)); CHKPTRQ(tmp); 476 PetscMemcpy((char *)tmp,(char *)xdata[0],mem_estimate*sizeof(int)); 477 PetscFree(xdata[0]); 478 xdata[0] = tmp; 479 mem_estimate = new_estimate; ++no_malloc; 480 for (ctr =1; ctr <=i; ctr++) { xdata[ctr] = xdata[ctr-1] + isz1[ctr-1];} 481 } 482 xdata[i][ct2++] = row;ct3++; 483 } 484 } 485 for ( k=oct2, max2 =ct2 ; k< max2; k++) { 486 row = xdata[i][k] - rstart; 487 start = ai[row]; 488 end = ai[row+1]; 489 for ( l=start; l < end; l++) { 490 val = aj[l] +ashift + cstart; 491 if(!BT_LOOKUP(xtable,val)) { 492 if (!(ct3 < mem_estimate)) { 493 new_estimate = 1.5*mem_estimate; 494 tmp = (int*) PetscMalloc(new_estimate * sizeof(int)); CHKPTRQ(tmp); 495 PetscMemcpy((char *)tmp,(char *)xdata[0],mem_estimate*sizeof(int)); 496 PetscFree(xdata[0]); 497 xdata[0] = tmp; 498 mem_estimate = new_estimate; ++no_malloc; 499 for (ctr =1; ctr <=i; ctr++) { xdata[ctr] = xdata[ctr-1] + isz1[ctr-1];} 500 } 501 xdata[i][ct2++] = val;ct3++; 502 } 503 } 504 start = bi[row]; 505 end = bi[row+1]; 506 for ( l=start; l < end; l++) { 507 val = garray[bj[l]+bshift] ; 508 if(!BT_LOOKUP(xtable,val)) { 509 if (!(ct3 < mem_estimate)) { 510 new_estimate = 1.5*mem_estimate; 511 tmp = (int*) PetscMalloc(new_estimate * sizeof(int)); CHKPTRQ(tmp); 512 PetscMemcpy((char *)tmp,(char *)xdata[0],mem_estimate*sizeof(int)); 513 PetscFree(xdata[0]); 514 xdata[0] = tmp; 515 mem_estimate = new_estimate; ++no_malloc; 516 for (ctr =1; ctr <=i; ctr++) { xdata[ctr] = xdata[ctr-1] + isz1[ctr-1];} 517 } 518 xdata[i][ct2++] = val;ct3++; 519 } 520 } 521 } 522 /* Update the header*/ 523 xdata[i][2*j] = ct2-oct2; /* Undo the vector isz1 and use only a var*/ 524 xdata[i][2*j-1] = rbuf[i][2*j-1]; 525 } 526 xdata[i][0] = rbuf[i][0]; 527 xdata[i+1] = xdata[i] +ct2; 528 isz1[i] = ct2; /* size of each message */ 529 } 530 PetscFree(xtable); 531 PLogInfo(0,"MatIncreaseOverlap_MPIAIJ:[%d] Allocated %d bytes, required %d bytes, no of mallocs = %d\n",rank,mem_estimate, ct3,no_malloc); 532 return 0; 533 } 534