#include "petscmat.h"
#include "petscis.h"
Go to the source code of this file.
Functions | |
int | MatSplitPoints (Mat A, int np, IS *splitpoints) |
int MatSplitPoints | ( | Mat | A, | |
int | np, | |||
IS * | splitpoints | |||
) |
Definition at line 308 of file icmk.c.
References ChainSplits(), CORNERUP, GetUpBiDiagSplits(), IF_TOO_FAR_RIGHT_BREAK, LEFTDOWN, LEFTUP, MatIsNull(), MatSubMatIsNull(), MPIAllGatherIntV(), OVERFLOW_DOWN, OVERFLOW_UP, RIGHTDOWN, RIGHTUP, SET_J, SET_LAST_COL, SET_STATUS, SET_TS, SET_TSS, SPLIT_BLOCK, Splits(), STATUS, VERY_FIRST_ROW, and VERY_LAST_ROW.
Referenced by compute_icm_splits().
00309 { 00310 MPI_Comm comm; 00311 int M,N,ntids,mytid, first,last, isFirst,isLast, 00312 nsplits,Nsplits,*split_array,*Split_array, 00313 *Splits,S_top, 00314 ierr; 00315 00316 PetscFunctionBegin; 00317 00318 ierr = PetscObjectGetComm((PetscObject)A,&comm); CHKERRQ(ierr); 00319 MPI_Comm_size(comm,&ntids); MPI_Comm_rank(comm,&mytid); 00320 isFirst = (mytid==0); isLast = (mytid==ntids-1); 00321 ierr = MatGetSize(A,&M,&N); CHKERRQ(ierr); 00322 ierr = MatGetOwnershipRange(A,&first,&last); CHKERRQ(ierr); 00323 00324 ierr = GetUpBiDiagSplits 00325 (A,first,last,isFirst,isLast, 00326 &split_array,&nsplits,&Split_array,&Nsplits); CHKERRQ(ierr); 00327 Split_array++; /* zero location is size */ 00328 00329 { 00330 IS *isses,*jsses; Mat *subs; 00331 int *splits,r_top,s_top, isplit,jsplit,*jsplits,*joffsets, 00332 ntest,nsave,*status; 00333 00334 #define SET_J j = Split_array[jsplit]; 00335 #define SET_TS ts = PetscMax(PetscMin(bs/10,i),1); 00336 #define SET_TSS \ 00337 tsr = PetscMin(ts,N-j); tsd = PetscMin(ts,M-i); tsc = PetscMin(ts,M-i); 00338 #define SET_LAST_COL \ 00339 ierr = MatGetRow(A,i,&ncols,&cols,PETSC_NULL); CHKERRQ(ierr); \ 00340 lastcol = cols[ncols-1]; \ 00341 ierr = MatRestoreRow(A,i,&ncols,&cols,PETSC_NULL); CHKERRQ(ierr); 00342 #define IF_TOO_FAR_RIGHT_BREAK if (j-ts>lastcol) break; 00343 /* 00344 * what is the number of splits for each of the local splits? 00345 * (we only need this for efficient alloation of the "status" array; 00346 * straightforward allocation of nsplits*Nsplits is far too much.) 00347 */ 00348 ierr = PetscMalloc(nsplits*sizeof(int),&joffsets); CHKERRQ(ierr); 00349 ierr = PetscMemzero(joffsets,nsplits*sizeof(int)); CHKERRQ(ierr); 00350 ierr = PetscMalloc((nsplits+1)*sizeof(int),&jsplits); CHKERRQ(ierr); 00351 ierr = PetscMemzero(jsplits,(nsplits+1)*sizeof(int)); CHKERRQ(ierr); 00352 jsplits[0] = 0; 00353 for (isplit=0; isplit<nsplits; isplit++) { 00354 int i = split_array[isplit],ncols,lastcol,njsplits=0; const int *cols; 00355 if (!VERY_LAST_ROW) { 00356 SET_LAST_COL; 00357 for (jsplit=0; jsplit<Nsplits; jsplit++) { 00358 int j,bs,ts; 00359 SET_J; bs = j-i; SET_TS; 00360 if (bs<=0) {joffsets[isplit]++; continue;} 00361 IF_TOO_FAR_RIGHT_BREAK; 00362 njsplits++; /* count how many splits we actually use */ 00363 } 00364 } 00365 jsplits[isplit+1] = jsplits[isplit]+njsplits; 00366 } 00367 ierr = PetscMalloc(jsplits[nsplits]*sizeof(int),&status); CHKERRQ(ierr); 00368 ierr = PetscMemzero(status,jsplits[nsplits]*sizeof(int)); CHKERRQ(ierr); 00369 #define LEFTDOWN 1 00370 #define RIGHTDOWN 2 00371 #define LEFTUP 4 00372 #define RIGHTUP 8 00373 #define CORNERUP 16 00374 #define STATUS(isp,jsp) status[jsplits[isp]+(jsp-joffsets[isplit])] 00375 #define SET_STATUS(isp,jsp,s) STATUS(isp,jsp) += s 00376 ntest = 3*jsplits[nsplits]; 00377 #if TRACE_MSGS 00378 printf("total # tests: %d\n",ntest); 00379 #endif 00380 ierr = PetscMalloc(ntest*sizeof(IS),&isses); CHKERRQ(ierr); 00381 ierr = PetscMalloc(ntest*sizeof(IS),&jsses); CHKERRQ(ierr); 00382 #define OVERFLOW_DOWN i+tsd-1>=last 00383 /* why did that say >=last-1 ? */ 00384 #define OVERFLOW_UP i-ts<first 00385 ntest = 0; nsave = 0; 00386 for (isplit=0; isplit<nsplits; isplit++) { 00387 int i=split_array[isplit],ncols,lastcol; 00388 const int *cols; PetscTruth flag; 00389 if (VERY_LAST_ROW) break; 00390 #if TRACE_MSGS > 1 00391 printf("from split %d=>%d: ",isplit,i); 00392 #endif 00393 SET_LAST_COL; 00394 for (jsplit=joffsets[isplit]; jsplit<Nsplits; jsplit++) { 00395 IS left,right,up,down,corner; int j,bs,ts,tsr,tsd,tsc; 00396 SET_J; bs = j-i; SET_TS; SET_TSS; 00397 IF_TOO_FAR_RIGHT_BREAK; 00398 #if TRACE_MSGS > 1 00399 printf("split %d=>%d, ",jsplit,j); 00400 #endif 00401 if ((OVERFLOW_DOWN) || (OVERFLOW_UP)) { 00402 ierr = ISCreateStride(MPI_COMM_SELF,tsr,j,1,&right); CHKERRQ(ierr); 00403 ierr = ISCreateStride(MPI_COMM_SELF,ts,j-ts,1,&left); CHKERRQ(ierr);} 00404 if (!VERY_LAST_ROW) { 00405 if (OVERFLOW_DOWN) { 00406 /*printf("saving down i=%d j=%d\n",i,j);*/ 00407 ierr = ISCreateStride(MPI_COMM_SELF,tsd,i,1,&down); CHKERRQ(ierr); 00408 isses[nsave] = left; jsses[nsave] = down; nsave++; 00409 isses[nsave] = right; jsses[nsave] = down; nsave++; 00410 } else { 00411 ierr = MatSubMatIsNull 00412 (A,i,i+tsd-1,j-ts,j-1,&flag); CHKERRQ(ierr); 00413 #if TRACE_MSGS > 1 00414 printf("submat (%d,%d) x (%d,%d) : n=%d, ", 00415 i,i+tsd-1,j-ts,j-1,flag); 00416 #endif 00417 SET_STATUS(isplit,jsplit,flag*LEFTDOWN); 00418 ierr = MatSubMatIsNull 00419 (A,i,i+tsd-1,j,j+tsr-1,&flag); CHKERRQ(ierr); 00420 #if TRACE_MSGS > 1 00421 printf("submat (%d,%d) x (%d,%d) : n=%d, ", 00422 i,i+tsd-1,j,j+tsr-1,flag); 00423 #endif 00424 SET_STATUS(isplit,jsplit,flag*RIGHTDOWN); 00425 } 00426 ntest += 2; 00427 } 00428 if (!VERY_FIRST_ROW) { 00429 if (OVERFLOW_UP) { 00430 ierr = ISCreateStride(MPI_COMM_SELF,ts,i-ts,1,&up); CHKERRQ(ierr); 00431 ierr = ISCreateStride(MPI_COMM_SELF,tsc,i,1,&corner); CHKERRQ(ierr); 00432 /*printf("saving up i=%d j=%d\n",i,j);*/ 00433 isses[nsave] = left; jsses[nsave] = up; nsave++; 00434 isses[nsave] = right; jsses[nsave] = up; nsave++; 00435 isses[nsave] = corner; jsses[nsave] = up; nsave++; 00436 } else { 00437 ierr = MatSubMatIsNull 00438 (A,i-ts,i-1,j-ts-1,j-1,&flag); CHKERRQ(ierr); 00439 #if TRACE_MSGS > 1 00440 printf("submat (%d,%d) x (%d,%d) : n=%d, ", 00441 i-ts,i-1,j-ts-1,j-1,flag); 00442 #endif 00443 SET_STATUS(isplit,jsplit,flag*LEFTUP); 00444 ierr = MatSubMatIsNull 00445 (A,i-ts,i-1,j,j+tsr-1,&flag); CHKERRQ(ierr); 00446 SET_STATUS(isplit,jsplit,flag*RIGHTUP); 00447 ierr = MatSubMatIsNull 00448 (A,i-ts,i-1,i,i+tsc-1,&flag); CHKERRQ(ierr); 00449 SET_STATUS(isplit,jsplit,flag*CORNERUP); 00450 } 00451 ntest += 3; 00452 } 00453 } 00454 #if TRACE_MSGS > 1 00455 printf("\n\n"); 00456 #endif 00457 } 00458 00459 #if TRACE_MSGS 00460 printf("saved %d straddling matrices out of %d tests\n",nsave,ntest); 00461 #endif 00462 ierr = MatGetSubMatrices 00463 (A,nsave,isses,jsses,MAT_INITIAL_MATRIX,&subs); CHKERRQ(ierr); 00464 00465 ierr = PetscMalloc(SPLIT_BLOCK*ntest*sizeof(int),&splits); CHKERRQ(ierr); 00466 ntest = 0; r_top = 0; s_top = 0; 00467 for (isplit=0; isplit<nsplits; isplit++) { 00468 int i=split_array[isplit],ncols,lastcol,restart=1,split=1,c=0; 00469 const int *cols; 00470 if (VERY_LAST_ROW) break; 00471 #if TRACE_MSGS > 1 00472 printf("from split %d=>%d: ",isplit,i); 00473 #endif 00474 SET_LAST_COL; 00475 for (jsplit=joffsets[isplit]; jsplit<Nsplits; jsplit++) { 00476 int j,bs,ts,tsr,tsd,tsc,rank=0; 00477 SET_J; bs = j-i; SET_TS; SET_TSS; 00478 #if TRACE_MSGS > 1 00479 if (j-ts>lastcol) printf("%d:a ",j); 00480 #endif 00481 IF_TOO_FAR_RIGHT_BREAK; 00482 if (!VERY_LAST_ROW) { 00483 Mat leftdown,rightdown; int flag_leftdown,flag_rightdown; 00484 if (OVERFLOW_DOWN) { 00485 leftdown = subs[ntest++]; rightdown = subs[ntest++]; 00486 ierr = MatIsNull(leftdown,&flag_leftdown); CHKERRQ(ierr); 00487 SET_STATUS(isplit,jsplit,flag_leftdown*LEFTDOWN); 00488 ierr = MatIsNull(rightdown,&flag_rightdown); CHKERRQ(ierr); 00489 SET_STATUS(isplit,jsplit,flag_rightdown*RIGHTDOWN); 00490 } else { 00491 flag_leftdown = STATUS(isplit,jsplit) & LEFTDOWN; 00492 flag_rightdown = STATUS(isplit,jsplit) & RIGHTDOWN; 00493 } 00494 restart = flag_leftdown && (!flag_rightdown || j==N); 00495 #if TRACE_MSGS > 1 00496 if (restart) printf("%d:dn:r, ",j); else printf("%d:dn:x, ",j); 00497 #endif 00498 if (restart) rank=1; 00499 } 00500 if (VERY_FIRST_ROW) { 00501 if (rank) rank=3; 00502 } else { 00503 Mat leftup,rightup,corner; int flag_leftup,flag_rightup,flag_corner; 00504 if (OVERFLOW_UP) { 00505 leftup = subs[ntest++]; rightup = subs[ntest++]; 00506 ierr = MatIsNull(leftup,&flag_leftup); CHKERRQ(ierr); 00507 ierr = MatIsNull(rightup,&flag_rightup); CHKERRQ(ierr); 00508 corner = subs[ntest++]; 00509 ierr = MatIsNull(corner,&flag_corner); CHKERRQ(ierr); 00510 } else { 00511 flag_leftup = STATUS(isplit,jsplit) & LEFTUP; 00512 flag_rightup = STATUS(isplit,jsplit) & RIGHTUP; 00513 flag_corner = STATUS(isplit,jsplit) & CORNERUP; 00514 } 00515 split = (flag_rightup || j==N) && !flag_leftup; 00516 #if TRACE_MSGS > 1 00517 if (split) printf("%d:up:s, ",j); else printf("%d:up:x, ",j); 00518 #endif 00519 if (rank && split) rank=2; 00520 if (rank==2 && flag_corner) rank=3; 00521 } 00522 if (rank>0) { 00523 c++; 00524 splits[s_top] = i; splits[s_top+1] = j; splits[s_top+2] = rank; 00525 s_top+=SPLIT_BLOCK;} 00526 } 00527 if (!c) { /* append next split as default case */ 00528 splits[s_top] = i; splits[s_top+1] = split_array[isplit+1]; 00529 splits[s_top+2] = 1; s_top+=SPLIT_BLOCK;} 00530 #if TRACE_MSGS > 1 00531 printf("\n"); 00532 #endif 00533 } 00534 splits[s_top] = splits[s_top-SPLIT_BLOCK+1]; s_top++; 00535 splits[s_top++] = N; splits[s_top++] = 1; 00536 #if TRACE_MSGS > 1 00537 { 00538 int i; 00539 printf("local final splits: "); 00540 for (i=0; i<s_top; i+=SPLIT_BLOCK) { 00541 printf("[%d,%d],",splits[i],splits[i+1]);} 00542 printf("\n");} 00543 #endif 00544 ierr = MPIAllGatherIntV(splits,s_top,&Splits,&S_top,comm); CHKERRQ(ierr); 00545 Splits++; 00546 #ifdef TRACE_MSGS 00547 if (mytid==0) { 00548 int i; 00549 printf("global final splits: "); 00550 for (i=0; i<S_top; i+=SPLIT_BLOCK) { 00551 printf("[%d,%d:%d],",Splits[i],Splits[i+1],Splits[i+2]); 00552 } 00553 printf("\n");} 00554 #endif 00555 ierr = PetscFree(splits); CHKERRQ(ierr); 00556 ierr = PetscFree(joffsets); CHKERRQ(ierr); 00557 ierr = PetscFree(jsplits); CHKERRQ(ierr); 00558 00559 for (isplit=0; isplit<nsave; isplit++) { 00560 ierr = ISDestroy(isses[isplit]); CHKERRQ(ierr); 00561 if (isplit>0 && (jsses[isplit]!=jsses[isplit-1])) { 00562 ierr = ISDestroy(jsses[isplit]); CHKERRQ(ierr); 00563 } 00564 } 00565 PetscFree(isses); PetscFree(jsses); 00566 00567 ierr = MatDestroyMatrices(ntest,&subs); CHKERRQ(ierr); /*also frees subs*/ 00568 } 00569 { 00570 IS isret; 00571 ierr = ChainSplits(N,Splits,S_top,np,&isret); CHKERRQ(ierr); 00572 *splitpoints = isret; 00573 } 00574 Splits--; 00575 ierr = PetscFree(Splits); CHKERRQ(ierr); 00576 00577 PetscFunctionReturn(0); 00578 }