#include "petscmat.h"#include "petscis.h"
Include dependency graph for icmk.h:
This graph shows which files directly or indirectly include this file:
Go to the source code of this file.
Functions | |
| int | MatSplitPoints (Mat A, int np, IS *splitpoints) |
Function Documentation
| 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);
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; PetscBool 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);
}
Here is the call graph for this function:
