SALSA Analysis Modules
Defines | Functions | Variables
icmk.c File Reference

Functions for computing the ICMK structure of a matrix. More...

#include <stdlib.h>
#include "anamod.h"
#include "anamodsalsamodules.h"
#include "icmk.h"
#include "petscmat.h"
#include "petscis.h"
#include "petsc.h"
#include "anamodutils.h"
Include dependency graph for icmk.c:

Go to the source code of this file.

Defines

#define CHDIF(a)   rar[a+1]-rar[a]
#define SPLIT_BLOCK   3
#define VERY_FIRST_ROW   (isFirst && isplit==0)
#define VERY_LAST_ROW   (isLast && isplit==nsplits-1)
#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
#define IF_TOO_FAR_RIGHT_BREAK   if (j-ts>lastcol) break;
#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
#define OVERFLOW_DOWN   i+tsd-1>=last
#define OVERFLOW_UP   i-ts<first

Functions

static void iSpaces (int len)
int MatIsNull (Mat A, int *flag)
int MatSubMatIsNull (Mat A, int i1, int i2, int j1, int j2, PetscTruth *isnull)
static PetscErrorCode PrintArray (int *ar, int len, char *txt)
static PetscErrorCode CondenseChain (int *chain, int len, int tar, IS *israr)
static PetscErrorCode CanChain (int lev, int b, int N, int *splits, int nsplits, int *chain, int len, int q, IS *rar)
int CondenseSplits (int p, IS *splits)
static PetscErrorCode ChainSplits (int N, int *splits, int s_top, int p, IS *rar)
static PetscErrorCode GetUpBiDiagSplits (Mat A, int first, int last, int isFirst, int isLast, int **ar, int *n_ar, int **Ar, int *N_ar)
int MatSplitPoints (Mat A, int np, IS *splitpoints)
int MatRedistribute (Mat *A, IS splits)
int VecRedistribute (Vec *V, IS splits)
static PetscErrorCode compute_icm_splits (Mat A)
static PetscErrorCode NSplits (AnaModNumericalProblem prob, AnalysisItem *rv, int *dummy, PetscTruth *flg)
static PetscErrorCode Splits (AnaModNumericalProblem prob, AnalysisItem *rv, int *dummy, PetscTruth *flg)
PetscErrorCode RegisterICMKModules (void)

Variables

double mq
int nc
int ml

Detailed Description

Functions for computing the ICMK structure of a matrix.

Inverse Cuthill-McKee distribution

ICMK (`Inverse Cuthill-McKee') is a heuristic (Eijkhout[2001]) that, based on the sparsity structure of a matrix, tries to find a meaningful block structure. This is done without permuting the matrix. This block structure can be used in Block Jacobi preconditioners: distributing the matrix over parallel processors according to the block structure often improves convergence.

Usage

Activate this module with:

PetscErrorCode RegisterICMKModules();

Compute these elements with

ComputeQuantity("icmk","element",A,&res,&flg);

Available elements are:

References

@techreport{Eijk:auto-block,
author = {Victor Eijkhout},
title = {Automatic Determination of Matrix Blocks},
institution = {Department of Computer Science, University of Tennessee},
number = {ut-cs-01-458},
note = {Lapack Working Note 151},
year = {2001}
}

Definition in file icmk.c.


Define Documentation

#define CHDIF (   a)    rar[a+1]-rar[a]

Referenced by CondenseChain().

#define CORNERUP   16

Referenced by MatSplitPoints().

#define IF_TOO_FAR_RIGHT_BREAK   if (j-ts>lastcol) break;

Referenced by MatSplitPoints().

#define LEFTDOWN   1

Referenced by MatSplitPoints().

#define LEFTUP   4

Referenced by MatSplitPoints().

#define OVERFLOW_DOWN   i+tsd-1>=last

Referenced by MatSplitPoints().

#define OVERFLOW_UP   i-ts<first

Referenced by MatSplitPoints().

#define RIGHTDOWN   2

Referenced by MatSplitPoints().

#define RIGHTUP   8

Referenced by MatSplitPoints().

#define SET_J   j = Split_array[jsplit];

Referenced by MatSplitPoints().

#define SET_LAST_COL
Value:
ierr = MatGetRow(A,i,&ncols,&cols,PETSC_NULL); CHKERRQ(ierr); \
      lastcol = cols[ncols-1]; \
      ierr = MatRestoreRow(A,i,&ncols,&cols,PETSC_NULL); CHKERRQ(ierr);

Referenced by MatSplitPoints().

#define SET_STATUS (   isp,
  jsp,
 
)    STATUS(isp,jsp) += s

Referenced by MatSplitPoints().

#define SET_TS   ts = PetscMax(PetscMin(bs/10,i),1);

Referenced by MatSplitPoints().

#define SET_TSS   tsr = PetscMin(ts,N-j); tsd = PetscMin(ts,M-i); tsc = PetscMin(ts,M-i);

Referenced by MatSplitPoints().

#define SPLIT_BLOCK   3

Definition at line 145 of file icmk.c.

Referenced by CanChain(), and MatSplitPoints().

#define STATUS (   isp,
  jsp 
)    status[jsplits[isp]+(jsp-joffsets[isplit])]

Referenced by MatSplitPoints().

#define VERY_FIRST_ROW   (isFirst && isplit==0)

Definition at line 253 of file icmk.c.

Referenced by MatSplitPoints().

#define VERY_LAST_ROW   (isLast && isplit==nsplits-1)

Definition at line 254 of file icmk.c.

Referenced by MatSplitPoints().


Function Documentation

static PetscErrorCode CanChain ( int  lev,
int  b,
int  N,
int *  splits,
int  nsplits,
int *  chain,
int  len,
int  q,
IS *  rar 
) [static]

Definition at line 153 of file icmk.c.

References PrintArray(), and SPLIT_BLOCK.

Referenced by ChainSplits().

{
  int nnext,next,*conts,ierr;
  double dq=(1.*q)/(len-1); 
  PetscFunctionBegin;
  if (len>=nsplits) SETERRQ(1,"chain overflow");
  chain[len] = b;
#if TRACE_MSGS
  printf("[%d] ",lev); PrintArray(chain,len+1,"chain so far");
#endif
  if (b==N) {
    ierr = ISCreateGeneral(MPI_COMM_SELF,len+1,chain,rar); CHKERRQ(ierr);
    PetscFunctionReturn(0);}
  {
    int split_start,split_end,isplit,next_block,contd,level;
    for (split_start=0; splits[split_start]<b; split_start+=SPLIT_BLOCK) ;
    /*
    if (splits[split_start]>b) {
#if TRACE_MSGS
      printf(".. there is no block starting at %d\n",b);
#endif
      PetscFunctionReturn(0);}
    */
    for (isplit=split_start;
         splits[isplit]==b && isplit<nsplits; isplit+=SPLIT_BLOCK) ;
    split_end = isplit-SPLIT_BLOCK;
    if (isplit>=nsplits)
      next_block = N;
    else
      next_block = splits[isplit/* +1 ??? */];
    nnext = isplit-split_start+SPLIT_BLOCK;
    ierr = PetscMalloc(nnext*sizeof(int),&conts); CHKERRQ(ierr);
    contd = 0;
    for (level=3; level>=1; level--)
      for (isplit=split_end; isplit>=split_start; isplit-=SPLIT_BLOCK)
        if (splits[isplit+2]==level) {
          ierr = PetscMemcpy
            (conts+contd,splits+isplit,SPLIT_BLOCK*sizeof(int)); CHKERRQ(ierr);
          contd += SPLIT_BLOCK;}
    if (next_block==N) {
      nnext -= SPLIT_BLOCK;
    } else {
      conts[contd] = b; conts[contd+1] = next_block; conts[contd+2] = 1;}
#if TRACE_MSGS
    PrintArray(conts,nnext,"continuations");
#endif
  }
  for (next=0; next<nnext; next+=SPLIT_BLOCK) {
#ifdef TRACE_MSGS
    if (next==0)
      printf("go on to level %d\n",lev+1);
    else
      printf("next attempt at level %d\n",lev+1);
#endif
    ierr = CanChain
      (lev+1,conts[next+1],N,splits,nsplits,chain,len+1,q,rar); CHKERRQ(ierr);
    if (rar) break;
  }
  PetscFunctionReturn(0);
}

Here is the call graph for this function:

static PetscErrorCode ChainSplits ( int  N,
int *  splits,
int  s_top,
int  p,
IS *  rar 
) [static]

Definition at line 235 of file icmk.c.

References CanChain(), CondenseSplits(), ml, mq, and nc.

Referenced by MatSplitPoints().

{
  int *chain,ierr; IS isret;

  PetscFunctionBegin;
  ierr = PetscMalloc(s_top*sizeof(int),&chain); CHKERRQ(ierr);
  mq = 0.; ml = 0; nc = 0;
  ierr = CanChain(0,0,N,splits,s_top,chain,0,0,&isret); CHKERRQ(ierr);
  if (!isret) SETERRQ(1,"Did not deliver a chain");
  ierr = PetscFree(chain); CHKERRQ(ierr);
  if (p) {ierr = CondenseSplits(p,&isret); CHKERRQ(ierr);}
#if TRACE_MSGS
  printf("final chain\n"); ISView(isret,0);
#endif
  *rar = isret;
  PetscFunctionReturn(0);
}

Here is the call graph for this function:

static PetscErrorCode compute_icm_splits ( Mat  A) [static]

Definition at line 671 of file icmk.c.

References GetDataID(), HASTOEXIST, id, and MatSplitPoints().

Referenced by NSplits(), and Splits().

{
  MPI_Comm comm; IS splits; int *indices,ns,*ss,id,ntids; 
  PetscTruth has; PetscErrorCode ierr;
  PetscFunctionBegin;

  ierr = PetscObjectGetComm((PetscObject)A,&comm); CHKERRQ(ierr);
  ierr = MPI_Comm_size(comm,&ntids); CHKERRQ(ierr);
  ierr = MatSplitPoints(A,0,&splits); CHKERRQ(ierr);

  ierr = ISGetSize(splits,&ns); CHKERRQ(ierr);
  ierr = ISGetIndices(splits,(const PetscInt**)&indices); CHKERRQ(ierr);
  ierr = PetscMalloc((ns+1)*sizeof(int),&ss); CHKERRQ(ierr);
  ierr = PetscMemcpy(ss+1,indices,ns*sizeof(int)); CHKERRQ(ierr);
  ss[0] = ns;
  ierr = ISRestoreIndices(splits,(const PetscInt**)&indices); CHKERRQ(ierr);
  ierr = ISDestroy(splits); CHKERRQ(ierr);

  ierr = GetDataID("icmk","n-splits",&id,&has); CHKERRQ(ierr);
  HASTOEXIST(has);
  ierr = PetscObjectComposedDataSetInt
    ((PetscObject)A,id,ns); CHKERRQ(ierr);
  ierr = GetDataID("icmk","splits",&id,&has); CHKERRQ(ierr);
  HASTOEXIST(has);
  ierr = PetscObjectComposedDataSetIntstar
    ((PetscObject)A,id,ss); CHKERRQ(ierr);

  PetscFunctionReturn(0);
}

Here is the call graph for this function:

static PetscErrorCode CondenseChain ( int *  chain,
int  len,
int  tar,
IS *  israr 
) [static]

Definition at line 114 of file icmk.c.

References CHDIF, and PrintArray().

Referenced by CondenseSplits().

{
  int i,ierr; int *rar; IS isret;
#define CHDIF(a) rar[a+1]-rar[a]
  PetscFunctionBegin;
  if (len<tar) SETERRQ(1,"chain length less than target");
  /* nasty: we are not allowed to do anything to 'chain', so we need to
   * overallocate the return array; just by a few elements, though */
  ierr = PetscMalloc((len+1)*sizeof(int),&rar); CHKERRQ(ierr);
  for (i=0; i<=len; i++) {rar[i] = chain[i];}
  while (len>tar) {
    int lloc=0,lmin=CHDIF(0),l,shift_start;
    for (i=1; i<len; i++) {
      l = CHDIF(i); if (l<lmin) {lloc=i; lmin=l;}}
    if (lloc==0 || (lloc<len-1 && CHDIF(lloc+1)<CHDIF(lloc-1))) {
      shift_start = lloc+1; } else { shift_start = lloc;}
    for (i=shift_start; i<len; i++) rar[i] = rar[i+1];
    len--;
  }
#ifdef TRACE_MSGS
  ierr = PrintArray(rar,len+1,"Condensed chain"); CHKERRQ(ierr);
#endif
  ierr = ISCreateGeneral(MPI_COMM_SELF,len+1,rar,&isret); CHKERRQ(ierr);
  *israr = isret;
  PetscFunctionReturn(0);
}

Here is the call graph for this function:

int CondenseSplits ( int  p,
IS *  splits 
)

Definition at line 216 of file icmk.c.

References CondenseChain().

Referenced by ChainSplits().

{
  IS newsplits; int lchain,*chain,ierr;
  PetscFunctionBegin;
  ierr = ISGetSize(*splits,&lchain); CHKERRQ(ierr);
  ierr = ISGetIndices(*splits,(const PetscInt**)&chain); CHKERRQ(ierr);
  ierr = CondenseChain(chain,lchain-1,p,&newsplits); CHKERRQ(ierr);
  ierr = ISRestoreIndices(*splits,(const PetscInt**)&chain); CHKERRQ(ierr);
  ierr = ISDestroy(*splits); CHKERRQ(ierr);
  *splits = newsplits;
  PetscFunctionReturn(0);
}

Here is the call graph for this function:

static PetscErrorCode GetUpBiDiagSplits ( Mat  A,
int  first,
int  last,
int  isFirst,
int  isLast,
int **  ar,
int *  n_ar,
int **  Ar,
int *  N_ar 
) [static]

Definition at line 259 of file icmk.c.

References MPIAllGatherIntV(), and TRUTH.

Referenced by MatSplitPoints().

{
  MPI_Comm comm; PetscTruth candidate;
  int *split_array,nsplits,*Split_array,Nsplits, local_size,row,ierr;

  PetscFunctionBegin;
  ierr = PetscObjectGetComm((PetscObject)A,&comm); CHKERRQ(ierr);
  local_size = last-first;
  ierr = PetscMalloc(local_size*sizeof(int),&split_array); CHKERRQ(ierr);
  nsplits = 0;
  if (isFirst) candidate=PETSC_TRUE;
  for (row=first; row<last; row++) {
    int ncols,icol; const int *cols; 
    ierr = MatGetRow(A,row,&ncols,&cols,PETSC_NULL); CHKERRQ(ierr);
    for (icol=0; icol<ncols; icol++) {
      if (cols[icol]==row) {
        if (candidate && icol<ncols-1 && cols[icol+1]==row+1)
          split_array[nsplits++] = row;
        candidate = TRUTH(icol==ncols-1 || cols[icol+1]>row+1);
        break;
      }
      if (cols[icol]>row) break;
    }
    ierr = MatRestoreRow(A,row,&ncols,&cols,PETSC_NULL); CHKERRQ(ierr);
  }
  if (isLast) split_array[nsplits++] = last;

  ierr = MPIAllGatherIntV
    (split_array,nsplits,&Split_array,&Nsplits,comm); CHKERRQ(ierr);
#if TRACE_MSGS
  {
    int i;
    printf("local upbidiag splits: ");
    for (i=0; i<nsplits; i++) {printf("%d,",split_array[i]);}
    printf("\n");
    printf("global upbidiag splits: ");
    for (i=0; i<Nsplits; i++) {printf("%d,",Split_array[i+1]);}
    printf("\n");
  }
#endif

  *ar = split_array; *n_ar = nsplits; *Ar = Split_array; *N_ar = Nsplits;

  PetscFunctionReturn(0);
}

Here is the call graph for this function:

static void iSpaces ( int  len) [static]

Definition at line 50 of file icmk.c.

{int i;for (i=0;i<len;i++){printf(" ");}}
int MatIsNull ( Mat  A,
int *  flag 
)

Definition at line 57 of file icmk.c.

Referenced by MatSplitPoints().

{
  int size,row,ncol,ierr;
  PetscFunctionBegin;
  ierr = MatGetSize(A,&size,&ncol); CHKERRQ(ierr);
  *flag = 1;
  for (row=0; row<size; row++) {
    ierr = MatGetRow(A,row,&ncol,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr);
    if (ncol>0) *flag = 0;
    ierr = MatRestoreRow(A,row,&ncol,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr);
    if (*flag==0) break;
  }
  PetscFunctionReturn(0);
}
int MatRedistribute ( Mat *  A,
IS  splits 
)

Definition at line 582 of file icmk.c.

{
  MPI_Comm comm; Mat B;
  int ntids,mytid,local_size,ierr;

  PetscFunctionBegin;
  ierr = PetscObjectGetComm((PetscObject)*A,&comm); CHKERRQ(ierr);
  MPI_Comm_size(comm,&ntids); MPI_Comm_rank(comm,&mytid);
  {
    /* find the improved distribution, return the new local_size */
    int *indices,chck;
    ierr = ISGetLocalSize(splits,&chck); CHKERRQ(ierr);
    if (chck!=ntids+1) SETERRQ(1,"Set of split points wrong length");
    ierr = ISGetIndices(splits,(const PetscInt**)&indices); CHKERRQ(ierr);
    local_size = indices[mytid+1]-indices[mytid];
    ierr = ISRestoreIndices(splits,(const PetscInt**)&indices); CHKERRQ(ierr);
  }
  {
    /* redistribute the matrix, unless the distribution is already as wanted */
    int perfect,gperfect,first,last,row;
    ierr = MatGetOwnershipRange(*A,&first,&last); CHKERRQ(ierr);
    perfect = local_size==(last-first);
    MPI_Allreduce(&perfect,&gperfect,1,MPI_INT,MPI_PROD,comm);
    if (gperfect) PetscFunctionReturn(0);
    ierr = MatCreateMPIAIJ /* optimise the d_ and o_ parameters ! */
      (comm,local_size,local_size,PETSC_DECIDE,PETSC_DECIDE,
       5,0,5,0,&B); CHKERRQ(ierr);
    for (row=first; row<last; row++) {
      int ncols; const int *cols; const PetscScalar *vals;
      ierr = MatGetRow(*A,row,&ncols,&cols,&vals); CHKERRQ(ierr);
      ierr = MatSetValues
        (B,1,&row,ncols,cols,vals,INSERT_VALUES); CHKERRQ(ierr);
      ierr = MatRestoreRow(*A,row,&ncols,&cols,&vals); CHKERRQ(ierr);
    }
    ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
    ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
  }
  ierr = MatDestroy(*A); CHKERRQ(ierr);
  *A = B;
  PetscFunctionReturn(0);
}
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:

int MatSubMatIsNull ( Mat  A,
int  i1,
int  i2,
int  j1,
int  j2,
PetscTruth *  isnull 
)

Definition at line 78 of file icmk.c.

Referenced by MatSplitPoints().

{
  int first,last,row,ncol,j,ierr; const int *cols;
  PetscFunctionBegin;
  ierr = MatGetOwnershipRange(A,&first,&last); CHKERRQ(ierr);
  if (i1<first || i2>=last) SETERRQ(1,"Requesting (part) non-local submat");
  *isnull = PETSC_TRUE;
  for (row=i1; row<=i2; row++) {
    ierr = MatGetRow(A,row,&ncol,&cols,PETSC_NULL); CHKERRQ(ierr);
    if (ncol>=0)
      if (cols[0]<=j2 && cols[ncol-1]>=j1)
        for (j=0; j<ncol; j++)
          if (cols[j]>=j1 && cols[j]<=j2) *isnull = PETSC_FALSE;
    ierr = MatRestoreRow(A,row,&ncol,&cols,PETSC_NULL); CHKERRQ(ierr);
    if (!*isnull) break;
  }
  PetscFunctionReturn(0);
}
static PetscErrorCode NSplits ( AnaModNumericalProblem  prob,
AnalysisItem rv,
int *  dummy,
PetscTruth *  flg 
) [static]

Definition at line 704 of file icmk.c.

References compute_icm_splits(), GetDataID(), HASTOEXIST, AnalysisItem::i, and id.

Referenced by RegisterICMKModules().

{
  Mat A = (Mat)prob;
  int v; PetscTruth has; int id; PetscErrorCode ierr;
  PetscFunctionBegin;
  ierr = GetDataID("icmk","n-splits",&id,&has); CHKERRQ(ierr);
  HASTOEXIST(has);
  ierr = PetscObjectComposedDataGetInt
    ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
  if (*flg) rv->i = v;
  else {
    ierr = compute_icm_splits(A); CHKERRQ(ierr);
    ierr = PetscObjectComposedDataGetInt
      ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
    if (*flg) rv->i = v;
  }
  
  PetscFunctionReturn(0);
}

Here is the call graph for this function:

static PetscErrorCode PrintArray ( int *  ar,
int  len,
char *  txt 
) [static]

Definition at line 99 of file icmk.c.

Referenced by CanChain(), and CondenseChain().

{
  int i;
  PetscFunctionBegin;
  printf("%s: ",txt); for (i=0; i<len; i++) printf("%d,",ar[i]); printf("\n");
  PetscFunctionReturn(0);
}
PetscErrorCode RegisterICMKModules ( void  )

Definition at line 749 of file icmk.c.

References ANALYSISINTARRAY, ANALYSISINTEGER, NSplits(), RegisterModule(), and Splits().

{
  PetscErrorCode ierr;
  PetscFunctionBegin;

  ierr = RegisterModule
    ("icmk","n-splits",ANALYSISINTEGER,&NSplits); CHKERRQ(ierr);
  ierr = RegisterModule
    ("icmk","splits",ANALYSISINTARRAY,&Splits); CHKERRQ(ierr);

  PetscFunctionReturn(0);
}

Here is the call graph for this function:

static PetscErrorCode Splits ( AnaModNumericalProblem  prob,
AnalysisItem rv,
int *  dummy,
PetscTruth *  flg 
) [static]

Definition at line 727 of file icmk.c.

References compute_icm_splits(), GetDataID(), HASTOEXIST, id, and AnalysisItem::ii.

Referenced by MatSplitPoints(), and RegisterICMKModules().

{
  Mat A = (Mat)prob;
  int *v; PetscTruth has; int id; PetscErrorCode ierr;
  PetscFunctionBegin;
  ierr = GetDataID("icmk","splits",&id,&has); CHKERRQ(ierr);
  HASTOEXIST(has);
  ierr = PetscObjectComposedDataGetIntstar
    ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
  if (*flg) rv->ii = v;
  else {
    ierr = compute_icm_splits(A); CHKERRQ(ierr);
    ierr = PetscObjectComposedDataGetIntstar
      ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
    if (*flg) rv->ii = v;
  }
  
  PetscFunctionReturn(0);
}

Here is the call graph for this function:

int VecRedistribute ( Vec *  V,
IS  splits 
)

Definition at line 623 of file icmk.c.

{
  MPI_Comm comm; Vec W;
  int ntids,mytid,local_size,new_first,ierr;

  PetscFunctionBegin;
  ierr = PetscObjectGetComm((PetscObject)*V,&comm); CHKERRQ(ierr);
  MPI_Comm_size(comm,&ntids); MPI_Comm_rank(comm,&mytid);
  {
    /* find the new local_size */
    int *indices;
    ierr = ISGetIndices(splits,(const PetscInt**)&indices); CHKERRQ(ierr);
    local_size = indices[mytid+1]-indices[mytid]; new_first = indices[mytid];
    ierr = ISRestoreIndices(splits,(const PetscInt**)&indices); CHKERRQ(ierr);
  }
  {
    /* redistribute the matrix, unless the distribution is already as wanted */
    int perfect,gperfect,first,last;
    ierr = VecGetOwnershipRange(*V,&first,&last); CHKERRQ(ierr);
    perfect = local_size==(last-first);
    MPI_Allreduce(&perfect,&gperfect,1,MPI_INT,MPI_PROD,comm);
    if (gperfect) PetscFunctionReturn(0);
    ierr = VecCreateMPI
      (comm,local_size,PETSC_DECIDE,&W); CHKERRQ(ierr);
    {
      PetscScalar *elements; int *range,i;
      ierr = VecGetArray(*V,&elements); CHKERRQ(ierr);
      ierr = PetscMalloc(local_size*sizeof(int),&range); CHKERRQ(ierr);
      for (i=first; i<last; i++) range[i-first] = i;
      ierr = VecSetValues
        (W,last-first,range,elements,INSERT_VALUES); CHKERRQ(ierr);
      ierr = PetscFree(range); CHKERRQ(ierr);
      ierr = VecRestoreArray(*V,&elements); CHKERRQ(ierr);
    }
    ierr = VecAssemblyBegin(W); CHKERRQ(ierr);
    ierr = VecAssemblyEnd(W); CHKERRQ(ierr);
  }
  ierr = VecDestroy(*V); CHKERRQ(ierr);
  *V = W;

  PetscFunctionReturn(0);
}

Variable Documentation

int ml

Definition at line 144 of file icmk.c.

Referenced by ChainSplits().

double mq

Definition at line 144 of file icmk.c.

Referenced by ChainSplits().

int nc

Definition at line 144 of file icmk.c.

Referenced by ChainSplits().