SALSA Analysis Modules
|
#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().
{ MPI_Comm comm; int M,N,ntids,mytid, first,last, isFirst,isLast, nsplits,Nsplits,*split_array,*Split_array, *Splits,S_top, ierr; PetscFunctionBegin; ierr = PetscObjectGetComm((PetscObject)A,&comm); CHKERRQ(ierr); MPI_Comm_size(comm,&ntids); MPI_Comm_rank(comm,&mytid); isFirst = (mytid==0); isLast = (mytid==ntids-1); ierr = MatGetSize(A,&M,&N); CHKERRQ(ierr); ierr = MatGetOwnershipRange(A,&first,&last); CHKERRQ(ierr); ierr = GetUpBiDiagSplits (A,first,last,isFirst,isLast, &split_array,&nsplits,&Split_array,&Nsplits); CHKERRQ(ierr); Split_array++; /* zero location is size */ { IS *isses,*jsses; Mat *subs; int *splits,r_top,s_top, isplit,jsplit,*jsplits,*joffsets, ntest,nsave,*status; #define SET_J j = Split_array[jsplit]; #define SET_TS ts = PetscMax(PetscMin(bs/10,i),1); #define SET_TSS \ tsr = PetscMin(ts,N-j); tsd = PetscMin(ts,M-i); tsc = PetscMin(ts,M-i); #define SET_LAST_COL \ ierr = MatGetRow(A,i,&ncols,&cols,PETSC_NULL); CHKERRQ(ierr); \ lastcol = cols[ncols-1]; \ ierr = MatRestoreRow(A,i,&ncols,&cols,PETSC_NULL); CHKERRQ(ierr); #define IF_TOO_FAR_RIGHT_BREAK if (j-ts>lastcol) break; /* * what is the number of splits for each of the local splits? * (we only need this for efficient alloation of the "status" array; * straightforward allocation of nsplits*Nsplits is far too much.) */ ierr = PetscMalloc(nsplits*sizeof(int),&joffsets); CHKERRQ(ierr); ierr = PetscMemzero(joffsets,nsplits*sizeof(int)); CHKERRQ(ierr); ierr = PetscMalloc((nsplits+1)*sizeof(int),&jsplits); CHKERRQ(ierr); ierr = PetscMemzero(jsplits,(nsplits+1)*sizeof(int)); CHKERRQ(ierr); jsplits[0] = 0; for (isplit=0; isplit<nsplits; isplit++) { int i = split_array[isplit],ncols,lastcol,njsplits=0; const int *cols; if (!VERY_LAST_ROW) { SET_LAST_COL; for (jsplit=0; jsplit<Nsplits; jsplit++) { int j,bs,ts; SET_J; bs = j-i; SET_TS; if (bs<=0) {joffsets[isplit]++; continue;} IF_TOO_FAR_RIGHT_BREAK; njsplits++; /* count how many splits we actually use */ } } jsplits[isplit+1] = jsplits[isplit]+njsplits; } ierr = PetscMalloc(jsplits[nsplits]*sizeof(int),&status); CHKERRQ(ierr); ierr = PetscMemzero(status,jsplits[nsplits]*sizeof(int)); CHKERRQ(ierr); #define LEFTDOWN 1 #define RIGHTDOWN 2 #define LEFTUP 4 #define RIGHTUP 8 #define CORNERUP 16 #define STATUS(isp,jsp) status[jsplits[isp]+(jsp-joffsets[isplit])] #define SET_STATUS(isp,jsp,s) STATUS(isp,jsp) += s ntest = 3*jsplits[nsplits]; #if TRACE_MSGS printf("total # tests: %d\n",ntest); #endif ierr = PetscMalloc(ntest*sizeof(IS),&isses); CHKERRQ(ierr); ierr = PetscMalloc(ntest*sizeof(IS),&jsses); CHKERRQ(ierr); #define OVERFLOW_DOWN i+tsd-1>=last /* why did that say >=last-1 ? */ #define OVERFLOW_UP i-ts<first ntest = 0; nsave = 0; for (isplit=0; isplit<nsplits; isplit++) { int i=split_array[isplit],ncols,lastcol; const int *cols; PetscTruth flag; if (VERY_LAST_ROW) break; #if TRACE_MSGS > 1 printf("from split %d=>%d: ",isplit,i); #endif SET_LAST_COL; for (jsplit=joffsets[isplit]; jsplit<Nsplits; jsplit++) { IS left,right,up,down,corner; int j,bs,ts,tsr,tsd,tsc; SET_J; bs = j-i; SET_TS; SET_TSS; IF_TOO_FAR_RIGHT_BREAK; #if TRACE_MSGS > 1 printf("split %d=>%d, ",jsplit,j); #endif if ((OVERFLOW_DOWN) || (OVERFLOW_UP)) { ierr = ISCreateStride(MPI_COMM_SELF,tsr,j,1,&right); CHKERRQ(ierr); ierr = ISCreateStride(MPI_COMM_SELF,ts,j-ts,1,&left); CHKERRQ(ierr);} if (!VERY_LAST_ROW) { if (OVERFLOW_DOWN) { /*printf("saving down i=%d j=%d\n",i,j);*/ ierr = ISCreateStride(MPI_COMM_SELF,tsd,i,1,&down); CHKERRQ(ierr); isses[nsave] = left; jsses[nsave] = down; nsave++; isses[nsave] = right; jsses[nsave] = down; nsave++; } else { ierr = MatSubMatIsNull (A,i,i+tsd-1,j-ts,j-1,&flag); CHKERRQ(ierr); #if TRACE_MSGS > 1 printf("submat (%d,%d) x (%d,%d) : n=%d, ", i,i+tsd-1,j-ts,j-1,flag); #endif SET_STATUS(isplit,jsplit,flag*LEFTDOWN); ierr = MatSubMatIsNull (A,i,i+tsd-1,j,j+tsr-1,&flag); CHKERRQ(ierr); #if TRACE_MSGS > 1 printf("submat (%d,%d) x (%d,%d) : n=%d, ", i,i+tsd-1,j,j+tsr-1,flag); #endif SET_STATUS(isplit,jsplit,flag*RIGHTDOWN); } ntest += 2; } if (!VERY_FIRST_ROW) { if (OVERFLOW_UP) { ierr = ISCreateStride(MPI_COMM_SELF,ts,i-ts,1,&up); CHKERRQ(ierr); ierr = ISCreateStride(MPI_COMM_SELF,tsc,i,1,&corner); CHKERRQ(ierr); /*printf("saving up i=%d j=%d\n",i,j);*/ isses[nsave] = left; jsses[nsave] = up; nsave++; isses[nsave] = right; jsses[nsave] = up; nsave++; isses[nsave] = corner; jsses[nsave] = up; nsave++; } else { ierr = MatSubMatIsNull (A,i-ts,i-1,j-ts-1,j-1,&flag); CHKERRQ(ierr); #if TRACE_MSGS > 1 printf("submat (%d,%d) x (%d,%d) : n=%d, ", i-ts,i-1,j-ts-1,j-1,flag); #endif SET_STATUS(isplit,jsplit,flag*LEFTUP); ierr = MatSubMatIsNull (A,i-ts,i-1,j,j+tsr-1,&flag); CHKERRQ(ierr); SET_STATUS(isplit,jsplit,flag*RIGHTUP); ierr = MatSubMatIsNull (A,i-ts,i-1,i,i+tsc-1,&flag); CHKERRQ(ierr); SET_STATUS(isplit,jsplit,flag*CORNERUP); } ntest += 3; } } #if TRACE_MSGS > 1 printf("\n\n"); #endif } #if TRACE_MSGS printf("saved %d straddling matrices out of %d tests\n",nsave,ntest); #endif ierr = MatGetSubMatrices (A,nsave,isses,jsses,MAT_INITIAL_MATRIX,&subs); CHKERRQ(ierr); ierr = PetscMalloc(SPLIT_BLOCK*ntest*sizeof(int),&splits); CHKERRQ(ierr); ntest = 0; r_top = 0; s_top = 0; for (isplit=0; isplit<nsplits; isplit++) { int i=split_array[isplit],ncols,lastcol,restart=1,split=1,c=0; const int *cols; if (VERY_LAST_ROW) break; #if TRACE_MSGS > 1 printf("from split %d=>%d: ",isplit,i); #endif SET_LAST_COL; for (jsplit=joffsets[isplit]; jsplit<Nsplits; jsplit++) { int j,bs,ts,tsr,tsd,tsc,rank=0; SET_J; bs = j-i; SET_TS; SET_TSS; #if TRACE_MSGS > 1 if (j-ts>lastcol) printf("%d:a ",j); #endif IF_TOO_FAR_RIGHT_BREAK; if (!VERY_LAST_ROW) { Mat leftdown,rightdown; int flag_leftdown,flag_rightdown; if (OVERFLOW_DOWN) { leftdown = subs[ntest++]; rightdown = subs[ntest++]; ierr = MatIsNull(leftdown,&flag_leftdown); CHKERRQ(ierr); SET_STATUS(isplit,jsplit,flag_leftdown*LEFTDOWN); ierr = MatIsNull(rightdown,&flag_rightdown); CHKERRQ(ierr); SET_STATUS(isplit,jsplit,flag_rightdown*RIGHTDOWN); } else { flag_leftdown = STATUS(isplit,jsplit) & LEFTDOWN; flag_rightdown = STATUS(isplit,jsplit) & RIGHTDOWN; } restart = flag_leftdown && (!flag_rightdown || j==N); #if TRACE_MSGS > 1 if (restart) printf("%d:dn:r, ",j); else printf("%d:dn:x, ",j); #endif if (restart) rank=1; } if (VERY_FIRST_ROW) { if (rank) rank=3; } else { Mat leftup,rightup,corner; int flag_leftup,flag_rightup,flag_corner; if (OVERFLOW_UP) { leftup = subs[ntest++]; rightup = subs[ntest++]; ierr = MatIsNull(leftup,&flag_leftup); CHKERRQ(ierr); ierr = MatIsNull(rightup,&flag_rightup); CHKERRQ(ierr); corner = subs[ntest++]; ierr = MatIsNull(corner,&flag_corner); CHKERRQ(ierr); } else { flag_leftup = STATUS(isplit,jsplit) & LEFTUP; flag_rightup = STATUS(isplit,jsplit) & RIGHTUP; flag_corner = STATUS(isplit,jsplit) & CORNERUP; } split = (flag_rightup || j==N) && !flag_leftup; #if TRACE_MSGS > 1 if (split) printf("%d:up:s, ",j); else printf("%d:up:x, ",j); #endif if (rank && split) rank=2; if (rank==2 && flag_corner) rank=3; } if (rank>0) { c++; splits[s_top] = i; splits[s_top+1] = j; splits[s_top+2] = rank; s_top+=SPLIT_BLOCK;} } if (!c) { /* append next split as default case */ splits[s_top] = i; splits[s_top+1] = split_array[isplit+1]; splits[s_top+2] = 1; s_top+=SPLIT_BLOCK;} #if TRACE_MSGS > 1 printf("\n"); #endif } splits[s_top] = splits[s_top-SPLIT_BLOCK+1]; s_top++; splits[s_top++] = N; splits[s_top++] = 1; #if TRACE_MSGS > 1 { int i; printf("local final splits: "); for (i=0; i<s_top; i+=SPLIT_BLOCK) { printf("[%d,%d],",splits[i],splits[i+1]);} printf("\n");} #endif ierr = MPIAllGatherIntV(splits,s_top,&Splits,&S_top,comm); CHKERRQ(ierr); Splits++; #ifdef TRACE_MSGS if (mytid==0) { int i; printf("global final splits: "); for (i=0; i<S_top; i+=SPLIT_BLOCK) { printf("[%d,%d:%d],",Splits[i],Splits[i+1],Splits[i+2]); } printf("\n");} #endif ierr = PetscFree(splits); CHKERRQ(ierr); ierr = PetscFree(joffsets); CHKERRQ(ierr); ierr = PetscFree(jsplits); CHKERRQ(ierr); for (isplit=0; isplit<nsave; isplit++) { ierr = ISDestroy(isses[isplit]); CHKERRQ(ierr); if (isplit>0 && (jsses[isplit]!=jsses[isplit-1])) { ierr = ISDestroy(jsses[isplit]); CHKERRQ(ierr); } } PetscFree(isses); PetscFree(jsses); ierr = MatDestroyMatrices(ntest,&subs); CHKERRQ(ierr); /*also frees subs*/ } { IS isret; ierr = ChainSplits(N,Splits,S_top,np,&isret); CHKERRQ(ierr); *splitpoints = isret; } Splits--; ierr = PetscFree(Splits); CHKERRQ(ierr); PetscFunctionReturn(0); }