SALSA Analysis Modules
Functions
icmk.h File Reference
#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);
  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);
}

Here is the call graph for this function: