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, PetscBool *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, PetscBool *flg)
static PetscErrorCode	Splits (AnaModNumericalProblem prob, AnalysisItem *rv, PetscBool *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:

• "nsplits" : number of split points found in the matrix
• "splits" : the sequence of split points, including 0 and N

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,
			s
	)		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(MPI_COMM_WORLD,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,PETSC_OWN_POINTER,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(MPI_COMM_WORLD,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; 
  PetscBool 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(MPI_COMM_WORLD,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,PETSC_OWN_POINTER,&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; PetscBool 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);
  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(MPI_COMM_WORLD,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);
    ierr = MPI_Allreduce(&perfect,&gperfect,1,MPI_INT,MPI_PROD,comm);
    if (gperfect) PetscFunctionReturn(0);
    ierr = MatCreateAIJ /* 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);
  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:

int MatSubMatIsNull	(	Mat	A,
		int	i1,
		int	i2,
		int	j1,
		int	j2,
		PetscBool *	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(MPI_COMM_WORLD,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,
		PetscBool *	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; PetscBool 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,
		PetscBool *	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; PetscBool 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);
  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);
    ierr = 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().