1 #ifndef lint 2 static char vcid[] = "$Id: mpiov.c,v 1.9 1996/01/31 20:22:46 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 20 int MatIncreaseOverlap_MPIAIJ_private(Mat C, int is_max, IS *is) 21 { 22 Mat_MPIAIJ *c = (Mat_MPIAIJ *) C->data; 23 Mat A = c->A, B = c->B; 24 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*b = (Mat_SeqAIJ*)B->data; 25 int **idx, *n, *w1, *w2, *w3, *w4, *rtable,**data; 26 int size, rank, m,i,j,k, ierr, **rbuf, row, proc, mct, msz, **outdat, **ptr; 27 int *ctr, sum, *pa, tag, *tmp,bsz, nmsg , *isz, *isz1, **xdata,rstart; 28 int cstart, *ai, *aj, *bi, *bj, *garray, bsz1, **rbuf2, ashift, bshift; 29 char **table, *xtable; 30 MPI_Comm comm; 31 MPI_Request *send_waits,*recv_waits,*send_waits2,*recv_waits2 ; 32 MPI_Status *send_status ,*recv_status; 33 double space, fr, maxs; 34 /* assume overlap = 1 */ 35 comm = C->comm; 36 tag = C->tag; 37 size = c->size; 38 rank = c->rank; 39 m = c->M; 40 rstart = c->rstart; 41 cstart = c->cstart; 42 ashift = a->indexshift; 43 ai = a->i; 44 aj = a->j +ashift; 45 bshift = b->indexshift; 46 bi = b->i; 47 bj = b->j +bshift; 48 garray = c->garray; 49 50 51 TrSpace( &space, &fr, &maxs ); 52 /* MPIU_printf(MPI_COMM_SELF,"[%d] acclcated space = %f fragments = %f max ever allocated = %f\n", rank, space, fr, maxs); */ 53 54 idx = (int **)PetscMalloc((is_max+1)*sizeof(int *)); CHKPTRQ(idx); 55 n = (int *)PetscMalloc((is_max+1)*sizeof(int )); CHKPTRQ(n); 56 rtable = (int *)PetscMalloc((m+1)*sizeof(int )); CHKPTRQ(rtable); 57 /* Hash table for maping row ->proc */ 58 59 for ( i=0 ; i<is_max ; i++) { 60 ierr = ISGetIndices(is[i],&idx[i]); CHKERRQ(ierr); 61 ierr = ISGetLocalSize(is[i],&n[i]); CHKERRQ(ierr); 62 } 63 64 /* Create hash table for the mapping :row -> proc*/ 65 for( i=0, j=0; i< size; i++) { 66 for (; j <c->rowners[i+1]; j++) { 67 rtable[j] = i; 68 } 69 } 70 71 /* evaluate communication - mesg to who, length of mesg, and buffer space 72 required. Based on this, buffers are allocated, and data copied into them*/ 73 w1 = (int *)PetscMalloc((size)*4*sizeof(int )); CHKPTRQ(w1); /* mesg size */ 74 w2 = w1 + size; /* if w2[i] marked, then a message to proc i*/ 75 w3 = w2 + size; /* no of IS that needs to be sent to proc i */ 76 w4 = w3 + size; /* temp work space used in determining w1, w2, w3 */ 77 PetscMemzero(w1,(size)*3*sizeof(int)); /* initialise work vector*/ 78 for ( i=0; i<is_max ; i++) { 79 PetscMemzero(w4,(size)*sizeof(int)); /* initialise work vector*/ 80 for ( j =0 ; j < n[i] ; j++) { 81 row = idx[i][j]; 82 proc = rtable[row]; 83 w4[proc]++; 84 } 85 for( j = 0; j < size; j++){ 86 if( w4[j] ) { w1[j] += w4[j]; w3[j] += 1;} 87 } 88 } 89 90 mct = 0; /* no of outgoing messages */ 91 msz = 0; /* total mesg length (for all proc */ 92 w1[rank] = 0; /* no mesg sent to intself */ 93 w3[rank] = 0; 94 for (i =0; i < size ; i++) { 95 if (w1[i]) { w2[i] = 1; mct++;} /* there exists a message to proc i */ 96 } 97 pa = (int *)PetscMalloc((mct +1)*sizeof(int)); CHKPTRQ(pa); /* (proc -array) */ 98 for (i =0, j=0; i < size ; i++) { 99 if (w1[i]) { pa[j] = i; j++; } 100 } 101 102 /* Each message would have a header = 1 + 2*(no of IS) + data */ 103 for (i = 0; i<mct ; i++) { 104 j = pa[i]; 105 w1[j] += w2[j] + 2* w3[j]; 106 msz += w1[j]; 107 } 108 109 /* Allocate Memory for outgoing messages */ 110 outdat = (int **)PetscMalloc( 2*size*sizeof(int*)); CHKPTRQ(outdat); 111 PetscMemzero(outdat, 2*size*sizeof(int*)); 112 tmp = (int *)PetscMalloc((msz+1) *sizeof (int)); CHKPTRQ(tmp); /*mrsg arr */ 113 ptr = outdat +size; /* Pointers to the data in outgoing buffers */ 114 ctr = (int *)PetscMalloc( size*sizeof(int)); CHKPTRQ(ctr); 115 116 { 117 int *iptr = tmp; 118 int ict = 0; 119 for (i = 0; i < mct ; i++) { 120 j = pa[i]; 121 iptr += ict; 122 outdat[j] = iptr; 123 ict = w1[j]; 124 } 125 } 126 127 /* Form the outgoing messages */ 128 /*plug in the headers*/ 129 for ( i=0 ; i<mct ; i++) { 130 j = pa[i]; 131 outdat[j][0] = 0; 132 PetscMemzero(outdat[j]+1, 2 * w3[j]*sizeof(int)); 133 ptr[j] = outdat[j] + 2*w3[j] +1; 134 } 135 136 /* Memory for doing local proc's work*/ 137 table = (char **)PetscMalloc((is_max+1)*sizeof(int *)); CHKPTRQ(table); 138 data = (int **)PetscMalloc((is_max+1)*sizeof(int *)); CHKPTRQ(data); 139 table[0] = (char *)PetscMalloc((m/8 +1)*(is_max)); CHKPTRQ(table[0]); 140 data [0] = (int *)PetscMalloc((m+1)*(is_max)*sizeof(int)); CHKPTRQ(data[0]); 141 142 for(i = 1; i<is_max ; i++) { 143 table[i] = table[0] + (m/8+1)*i; 144 data[i] = data[0] + (m+1)*i; 145 } 146 147 PetscMemzero((void*)*table,(m/8+1)*(is_max)); 148 isz = (int *)PetscMalloc((is_max+1) *sizeof(int)); CHKPTRQ(isz); 149 PetscMemzero((void *)isz,(is_max+1) *sizeof(int)); 150 151 /* Parse the IS and update local tables and the outgoing buf with the data*/ 152 for ( i=0 ; i<is_max ; i++) { 153 PetscMemzero(ctr,size*sizeof(int)); 154 for( j=0; j<n[i]; j++) { /* parse the indices of each IS */ 155 row = idx[i][j]; 156 proc = rtable[row]; 157 if (proc != rank) { /* copy to the outgoing buf*/ 158 ctr[proc]++; 159 *ptr[proc] = row; 160 ptr[proc]++; 161 } 162 else { /* Update the table */ 163 if ( !BT_LOOKUP(table[i],row)) { data[i][isz[i]++] = row;} 164 /* if(!table[i][row]++) { data[i][isz[i]++] = row;}*/ 165 } 166 } 167 /* Update the headers for the current IS */ 168 for( j = 0; j<size; j++) { /* Can Optimise this loop too */ 169 if (ctr[j]) { 170 k= ++outdat[j][0]; 171 outdat[j][2*k] = ctr[j]; 172 outdat[j][2*k-1] = i; 173 } 174 } 175 } 176 177 /* Check Validity of the outgoing messages */ 178 for ( i=0 ; i<mct ; i++) { 179 j = pa[i]; 180 if (w3[j] != outdat[j][0]) {SETERRQ(1,"MatIncreaseOverlap_MPIAIJ: Blew it! Header[1] mismatch!\n"); } 181 } 182 183 for ( i=0 ; i<mct ; i++) { 184 j = pa[i]; 185 sum = 1; 186 for (k = 1; k <= w3[j]; k++) sum += outdat[j][2*k]+2; 187 if (sum != w1[j]) { SETERRQ(1,"MatIncreaseOverlap_MPIAIJ: Blew it! Header[2-n] mismatch! \n"); } 188 } 189 190 191 /* Do a global reduction to determine how many messages to expect*/ 192 { 193 int *rw1, *rw2; 194 rw1 = (int *)PetscMalloc(2*size*sizeof(int)); CHKPTRQ(rw1); 195 rw2 = rw1+size; 196 MPI_Allreduce((void *)w1, rw1, size, MPI_INT, MPI_MAX, comm); 197 bsz = rw1[rank]; 198 MPI_Allreduce((void *)w2, rw2, size, MPI_INT, MPI_SUM, comm); 199 nmsg = rw2[rank]; 200 PetscFree(rw1); 201 } 202 203 /* Allocate memory for recv buffers . Prob none if nmsg = 0 ???? */ 204 rbuf = (int**) PetscMalloc((nmsg+1) *sizeof(int*)); CHKPTRQ(rbuf); 205 rbuf[0] = (int *) PetscMalloc((nmsg *bsz+1) * sizeof(int)); CHKPTRQ(rbuf[0]); 206 for (i=1; i<nmsg ; ++i) rbuf[i] = rbuf[i-1] + bsz; 207 208 /* Now post the receives */ 209 recv_waits = (MPI_Request *) PetscMalloc((nmsg+1)*sizeof(MPI_Request)); 210 CHKPTRQ(recv_waits); 211 for ( i=0; i<nmsg; ++i){ 212 MPI_Irecv((void *)rbuf[i], bsz, MPI_INT, MPI_ANY_SOURCE, tag, comm, recv_waits+i); 213 } 214 215 /* Now post the sends */ 216 send_waits = (MPI_Request *) PetscMalloc((mct+1)*sizeof(MPI_Request)); 217 CHKPTRQ(send_waits); 218 for( i =0; i< mct; ++i){ 219 j = pa[i]; 220 MPI_Isend( (void *)outdat[j], w1[j], MPI_INT, j, tag, comm, send_waits+i); 221 } 222 223 /* Do Local work*/ 224 /* Extract the matrices */ 225 { 226 int start, end, val, max; 227 228 for( i=0; i<is_max; i++) { 229 for ( j=0, max =isz[i] ; j< max; j++) { 230 row = data[i][j] - rstart; 231 start = ai[row]; 232 end = ai[row+1]; 233 for ( k=start; k < end; k++) { /* Amat */ 234 val = aj[k] + ashift + cstart; 235 if(!BT_LOOKUP(table[i],val)) { data[i][isz[i]++] = val;} 236 /* if(!table[i][val]++) { data[i][isz[i]++] = val;} */ 237 } 238 start = bi[row]; 239 end = bi[row+1]; 240 for ( k=start; k < end; k++) { /* Bmat */ 241 val = garray[bj[k]+bshift] ; 242 if(! BT_LOOKUP(table[i],val)) { data[i][isz[i]++] = val;} 243 /* if(!table[i][val]++) { data[i][isz[i]++] = val;} */ 244 } 245 } 246 } 247 } 248 /* Receive messages*/ 249 { 250 int index; 251 252 recv_status = (MPI_Status *) PetscMalloc( (nmsg+1)*sizeof(MPI_Status) ); 253 CHKPTRQ(recv_status); 254 for ( i=0; i< nmsg; ++i) { 255 MPI_Waitany(nmsg, recv_waits, &index, recv_status+i); 256 } 257 258 send_status = (MPI_Status *) PetscMalloc( (mct+1)*sizeof(MPI_Status) ); 259 CHKPTRQ(send_status); 260 MPI_Waitall(mct,send_waits,send_status); 261 } 262 263 264 /* Do work on recieved messages*/ 265 { 266 int total_sz,ct, ct1, ct2, mem_estimate; 267 int oct2, l, start, end, val, max1, max2; 268 269 270 for (i =0, ct =0, total_sz =0; i< nmsg; ++i){ 271 ct+= rbuf[i][0]; 272 for ( j = 1; j <= rbuf[i][0] ; j++ ) { total_sz += rbuf[i][2*j]; } 273 } 274 max1 = ct*(a->nz +b->nz)/c->m; 275 mem_estimate = 3* (total_sz > max1?total_sz:max1) 276 ; 277 xdata = (int **)PetscMalloc((nmsg+1)*sizeof(int *)); CHKPTRQ(xdata); 278 xdata[0] = (int *)PetscMalloc(mem_estimate *sizeof(int)); CHKPTRQ(xdata[0]); 279 xtable = (char *)PetscMalloc((m/8+1)); CHKPTRQ(xtable); 280 isz1 = (int *)PetscMalloc((nmsg+1) *sizeof(int)); CHKPTRQ(isz1); 281 PetscMemzero((void *)isz1,(nmsg+1) *sizeof(int)); 282 283 284 for (i =0; i< nmsg; i++) { /* for easch mesg from proc i */ 285 ct1 = 2*rbuf[i][0]+1; 286 ct2 = 2*rbuf[i][0]+1; 287 for (j = 1, max1= rbuf[i][0]; j<=max1; j++) { /* for each IS from proc i*/ 288 PetscMemzero((void *)xtable,(m/8+1)); 289 oct2 = ct2; 290 for (k =0; k < rbuf[i][2*j]; k++, ct1++) { 291 row = rbuf[i][ct1]; 292 if(!BT_LOOKUP(xtable,row)) { xdata[i][ct2++] = row;} 293 /* if(!xtable[row]++) { xdata[i][ct2++] = row;} */ 294 } 295 for ( k=oct2, max2 =ct2 ; k< max2; k++) { 296 row = xdata[i][k] - rstart; 297 start = ai[row]; 298 end = ai[row+1]; 299 for ( l=start; l < end; l++) { 300 val = aj[l] +ashift + cstart; 301 if(!BT_LOOKUP(xtable,val)) { xdata[i][ct2++] = val;} 302 /* if(!xtable[val]++) { xdata[i][ct2++] = val;} */ 303 } 304 start = bi[row]; 305 end = bi[row+1]; 306 for ( l=start; l < end; l++) { 307 val = garray[bj[l]+bshift] ; 308 if(!BT_LOOKUP(xtable,val)) { xdata[i][ct2++] = val;} 309 /* if(!xtable[val]++) { xdata[i][ct2++] = val;} */ 310 } 311 } 312 /* Update the header*/ 313 xdata[i][2*j] = ct2-oct2; /* Undo the vector isz1 and use only a var*/ 314 xdata[i][2*j-1] = rbuf[i][2*j-1]; 315 } 316 xdata[i][0] = rbuf[i][0]; 317 xdata[i+1] = xdata[i] +ct2; 318 isz1[i] = ct2; /* size of each message */ 319 } 320 321 } 322 /* need isz, xdata;*/ 323 324 /* Send the data back*/ 325 /* Do a global reduction to know the buffer space req for incoming messages*/ 326 { 327 int *rw1, *rw2; 328 329 rw1 = (int *)PetscMalloc(2*size*sizeof(int)); CHKPTRQ(rw1); 330 PetscMemzero((void*)rw1,2*size*sizeof(int)); 331 rw2 = rw1+size; 332 for (i =0; i < nmsg ; ++i) { 333 proc = recv_status[i].MPI_SOURCE; 334 rw1[proc] = isz1[i]; 335 } 336 337 MPI_Allreduce((void *)rw1, (void *)rw2, size, MPI_INT, MPI_MAX, comm); 338 bsz1 = rw2[rank]; 339 PetscFree(rw1); 340 } 341 342 /* Allocate buffers*/ 343 344 /* Allocate memory for recv buffers . Prob none if nmsg = 0 ???? */ 345 rbuf2 = (int**) PetscMalloc((mct+1) *sizeof(int*)); CHKPTRQ(rbuf2); 346 rbuf2[0] = (int *) PetscMalloc((mct*bsz1+1) * sizeof(int)); CHKPTRQ(rbuf2[0]); 347 for (i=1; i<mct ; ++i) rbuf2[i] = rbuf2[i-1] + bsz1; 348 349 /* Now post the receives */ 350 recv_waits2 = (MPI_Request *)PetscMalloc((mct+1)*sizeof(MPI_Request)); CHKPTRQ(recv_waits2) 351 CHKPTRQ(recv_waits2); 352 for ( i=0; i<mct; ++i){ 353 MPI_Irecv((void *)rbuf2[i], bsz1, MPI_INT, MPI_ANY_SOURCE, tag, comm, recv_waits2+i); 354 } 355 356 /* Now post the sends */ 357 send_waits2 = (MPI_Request *) PetscMalloc((nmsg+1)*sizeof(MPI_Request)); 358 CHKPTRQ(send_waits2); 359 for( i =0; i< nmsg; ++i){ 360 j = recv_status[i].MPI_SOURCE; 361 MPI_Isend( (void *)xdata[i], isz1[i], MPI_INT, j, tag, comm, send_waits2+i); 362 } 363 364 /* recieve work done on other processors*/ 365 { 366 int index, is_no, ct1, max; 367 MPI_Status *send_status2 ,*recv_status2; 368 369 recv_status2 = (MPI_Status *) PetscMalloc( (mct+1)*sizeof(MPI_Status) ); 370 CHKPTRQ(recv_status2); 371 372 373 for ( i=0; i< mct; ++i) { 374 MPI_Waitany(mct, recv_waits2, &index, recv_status2+i); 375 /* Process the message*/ 376 ct1 = 2*rbuf2[index][0]+1; 377 for (j=1; j<=rbuf2[index][0]; j++) { 378 max = rbuf2[index][2*j]; 379 is_no = rbuf2[index][2*j-1]; 380 for (k=0; k < max ; k++, ct1++) { 381 row = rbuf2[index][ct1]; 382 if(!BT_LOOKUP(table[is_no],row)) { data[is_no][isz[is_no]++] = row;} 383 /* if(!table[is_no][row]++) { data[is_no][isz[is_no]++] = row;} */ 384 } 385 } 386 } 387 388 389 send_status2 = (MPI_Status *) PetscMalloc( (nmsg+1)*sizeof(MPI_Status) ); 390 CHKPTRQ(send_status2); 391 MPI_Waitall(nmsg,send_waits2,send_status2); 392 393 PetscFree(send_status2); PetscFree(recv_status2); 394 } 395 for( i=0; i< is_max; ++i) { 396 ierr = ISRestoreIndices(is[i], idx+i); CHKERRQ(ierr); 397 } 398 for( i=0; i< is_max; ++i) { 399 ierr = ISDestroy(is[i]); CHKERRQ(ierr); 400 } 401 for ( i=0; i<is_max; ++i) { 402 ierr = ISCreateSeq(MPI_COMM_SELF, isz[i], data[i], is+i); CHKERRQ(ierr); 403 } 404 405 TrSpace( &space, &fr, &maxs ); 406 /* MPIU_printf(MPI_COMM_SELF,"[%d] acclcated space = %f fragments = %f max ever allocated = %f\n", rank, space, fr, maxs);*/ 407 408 /* whew! already done? check again :) */ 409 PetscFree(idx); 410 PetscFree(n); 411 PetscFree(rtable); 412 PetscFree(w1); 413 PetscFree(tmp); 414 PetscFree(outdat); 415 PetscFree(ctr); 416 PetscFree(pa); 417 PetscFree(rbuf[0]); 418 PetscFree(rbuf); 419 PetscFree(rbuf2[0]); 420 PetscFree(rbuf2); 421 PetscFree(send_waits); 422 PetscFree(recv_waits); 423 PetscFree(send_waits2); 424 PetscFree(recv_waits2); 425 PetscFree(table[0]); 426 PetscFree(table); 427 PetscFree(data[0]); 428 PetscFree(data); 429 PetscFree(send_status); 430 PetscFree(recv_status); 431 PetscFree(isz); 432 PetscFree(isz1); 433 PetscFree(xtable); 434 PetscFree(xdata[0]); 435 PetscFree(xdata); 436 437 /* Dont forget to ISRestoreIndices */ 438 return 0; 439 } 440