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 (  )     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().

00154 {
00155   int nnext,next,*conts,ierr;
00156   double dq=(1.*q)/(len-1); 
00157   PetscFunctionBegin;
00158   if (len>=nsplits) SETERRQ(1,"chain overflow");
00159   chain[len] = b;
00160 #if TRACE_MSGS
00161   printf("[%d] ",lev); PrintArray(chain,len+1,"chain so far");
00162 #endif
00163   if (b==N) {
00164     ierr = ISCreateGeneral(MPI_COMM_SELF,len+1,chain,rar); CHKERRQ(ierr);
00165     PetscFunctionReturn(0);}
00166   {
00167     int split_start,split_end,isplit,next_block,contd,level;
00168     for (split_start=0; splits[split_start]<b; split_start+=SPLIT_BLOCK) ;
00169     /*
00170     if (splits[split_start]>b) {
00171 #if TRACE_MSGS
00172       printf(".. there is no block starting at %d\n",b);
00173 #endif
00174       PetscFunctionReturn(0);}
00175     */
00176     for (isplit=split_start;
00177          splits[isplit]==b && isplit<nsplits; isplit+=SPLIT_BLOCK) ;
00178     split_end = isplit-SPLIT_BLOCK;
00179     if (isplit>=nsplits)
00180       next_block = N;
00181     else
00182       next_block = splits[isplit/* +1 ??? */];
00183     nnext = isplit-split_start+SPLIT_BLOCK;
00184     ierr = PetscMalloc(nnext*sizeof(int),&conts); CHKERRQ(ierr);
00185     contd = 0;
00186     for (level=3; level>=1; level--)
00187       for (isplit=split_end; isplit>=split_start; isplit-=SPLIT_BLOCK)
00188         if (splits[isplit+2]==level) {
00189           ierr = PetscMemcpy
00190             (conts+contd,splits+isplit,SPLIT_BLOCK*sizeof(int)); CHKERRQ(ierr);
00191           contd += SPLIT_BLOCK;}
00192     if (next_block==N) {
00193       nnext -= SPLIT_BLOCK;
00194     } else {
00195       conts[contd] = b; conts[contd+1] = next_block; conts[contd+2] = 1;}
00196 #if TRACE_MSGS
00197     PrintArray(conts,nnext,"continuations");
00198 #endif
00199   }
00200   for (next=0; next<nnext; next+=SPLIT_BLOCK) {
00201 #ifdef TRACE_MSGS
00202     if (next==0)
00203       printf("go on to level %d\n",lev+1);
00204     else
00205       printf("next attempt at level %d\n",lev+1);
00206 #endif
00207     ierr = CanChain
00208       (lev+1,conts[next+1],N,splits,nsplits,chain,len+1,q,rar); CHKERRQ(ierr);
00209     if (rar) break;
00210   }
00211   PetscFunctionReturn(0);
00212 }

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().

00236 {
00237   int *chain,ierr; IS isret;
00238 
00239   PetscFunctionBegin;
00240   ierr = PetscMalloc(s_top*sizeof(int),&chain); CHKERRQ(ierr);
00241   mq = 0.; ml = 0; nc = 0;
00242   ierr = CanChain(0,0,N,splits,s_top,chain,0,0,&isret); CHKERRQ(ierr);
00243   if (!isret) SETERRQ(1,"Did not deliver a chain");
00244   ierr = PetscFree(chain); CHKERRQ(ierr);
00245   if (p) {ierr = CondenseSplits(p,&isret); CHKERRQ(ierr);}
00246 #if TRACE_MSGS
00247   printf("final chain\n"); ISView(isret,0);
00248 #endif
00249   *rar = isret;
00250   PetscFunctionReturn(0);
00251 }

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().

00672 {
00673   MPI_Comm comm; IS splits; int *indices,ns,*ss,id,ntids; 
00674   PetscTruth has; PetscErrorCode ierr;
00675   PetscFunctionBegin;
00676 
00677   ierr = PetscObjectGetComm((PetscObject)A,&comm); CHKERRQ(ierr);
00678   ierr = MPI_Comm_size(comm,&ntids); CHKERRQ(ierr);
00679   ierr = MatSplitPoints(A,0,&splits); CHKERRQ(ierr);
00680 
00681   ierr = ISGetSize(splits,&ns); CHKERRQ(ierr);
00682   ierr = ISGetIndices(splits,(const PetscInt**)&indices); CHKERRQ(ierr);
00683   ierr = PetscMalloc((ns+1)*sizeof(int),&ss); CHKERRQ(ierr);
00684   ierr = PetscMemcpy(ss+1,indices,ns*sizeof(int)); CHKERRQ(ierr);
00685   ss[0] = ns;
00686   ierr = ISRestoreIndices(splits,(const PetscInt**)&indices); CHKERRQ(ierr);
00687   ierr = ISDestroy(splits); CHKERRQ(ierr);
00688 
00689   ierr = GetDataID("icmk","n-splits",&id,&has); CHKERRQ(ierr);
00690   HASTOEXIST(has);
00691   ierr = PetscObjectComposedDataSetInt
00692     ((PetscObject)A,id,ns); CHKERRQ(ierr);
00693   ierr = GetDataID("icmk","splits",&id,&has); CHKERRQ(ierr);
00694   HASTOEXIST(has);
00695   ierr = PetscObjectComposedDataSetIntstar
00696     ((PetscObject)A,id,ss); CHKERRQ(ierr);
00697 
00698   PetscFunctionReturn(0);
00699 }

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().

00115 {
00116   int i,ierr; int *rar; IS isret;
00117 #define CHDIF(a) rar[a+1]-rar[a]
00118   PetscFunctionBegin;
00119   if (len<tar) SETERRQ(1,"chain length less than target");
00120   /* nasty: we are not allowed to do anything to 'chain', so we need to
00121    * overallocate the return array; just by a few elements, though */
00122   ierr = PetscMalloc((len+1)*sizeof(int),&rar); CHKERRQ(ierr);
00123   for (i=0; i<=len; i++) {rar[i] = chain[i];}
00124   while (len>tar) {
00125     int lloc=0,lmin=CHDIF(0),l,shift_start;
00126     for (i=1; i<len; i++) {
00127       l = CHDIF(i); if (l<lmin) {lloc=i; lmin=l;}}
00128     if (lloc==0 || (lloc<len-1 && CHDIF(lloc+1)<CHDIF(lloc-1))) {
00129       shift_start = lloc+1; } else { shift_start = lloc;}
00130     for (i=shift_start; i<len; i++) rar[i] = rar[i+1];
00131     len--;
00132   }
00133 #ifdef TRACE_MSGS
00134   ierr = PrintArray(rar,len+1,"Condensed chain"); CHKERRQ(ierr);
00135 #endif
00136   ierr = ISCreateGeneral(MPI_COMM_SELF,len+1,rar,&isret); CHKERRQ(ierr);
00137   *israr = isret;
00138   PetscFunctionReturn(0);
00139 }

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().

00217 {
00218   IS newsplits; int lchain,*chain,ierr;
00219   PetscFunctionBegin;
00220   ierr = ISGetSize(*splits,&lchain); CHKERRQ(ierr);
00221   ierr = ISGetIndices(*splits,(const PetscInt**)&chain); CHKERRQ(ierr);
00222   ierr = CondenseChain(chain,lchain-1,p,&newsplits); CHKERRQ(ierr);
00223   ierr = ISRestoreIndices(*splits,(const PetscInt**)&chain); CHKERRQ(ierr);
00224   ierr = ISDestroy(*splits); CHKERRQ(ierr);
00225   *splits = newsplits;
00226   PetscFunctionReturn(0);
00227 }

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().

00261 {
00262   MPI_Comm comm; PetscTruth candidate;
00263   int *split_array,nsplits,*Split_array,Nsplits, local_size,row,ierr;
00264 
00265   PetscFunctionBegin;
00266   ierr = PetscObjectGetComm((PetscObject)A,&comm); CHKERRQ(ierr);
00267   local_size = last-first;
00268   ierr = PetscMalloc(local_size*sizeof(int),&split_array); CHKERRQ(ierr);
00269   nsplits = 0;
00270   if (isFirst) candidate=PETSC_TRUE;
00271   for (row=first; row<last; row++) {
00272     int ncols,icol; const int *cols; 
00273     ierr = MatGetRow(A,row,&ncols,&cols,PETSC_NULL); CHKERRQ(ierr);
00274     for (icol=0; icol<ncols; icol++) {
00275       if (cols[icol]==row) {
00276         if (candidate && icol<ncols-1 && cols[icol+1]==row+1)
00277           split_array[nsplits++] = row;
00278         candidate = TRUTH(icol==ncols-1 || cols[icol+1]>row+1);
00279         break;
00280       }
00281       if (cols[icol]>row) break;
00282     }
00283     ierr = MatRestoreRow(A,row,&ncols,&cols,PETSC_NULL); CHKERRQ(ierr);
00284   }
00285   if (isLast) split_array[nsplits++] = last;
00286 
00287   ierr = MPIAllGatherIntV
00288     (split_array,nsplits,&Split_array,&Nsplits,comm); CHKERRQ(ierr);
00289 #if TRACE_MSGS
00290   {
00291     int i;
00292     printf("local upbidiag splits: ");
00293     for (i=0; i<nsplits; i++) {printf("%d,",split_array[i]);}
00294     printf("\n");
00295     printf("global upbidiag splits: ");
00296     for (i=0; i<Nsplits; i++) {printf("%d,",Split_array[i+1]);}
00297     printf("\n");
00298   }
00299 #endif
00300 
00301   *ar = split_array; *n_ar = nsplits; *Ar = Split_array; *N_ar = Nsplits;
00302 
00303   PetscFunctionReturn(0);
00304 }

Here is the call graph for this function:

static void iSpaces ( int  len  )  [static]

Definition at line 50 of file icmk.c.

00050 {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().

00058 {
00059   int size,row,ncol,ierr;
00060   PetscFunctionBegin;
00061   ierr = MatGetSize(A,&size,&ncol); CHKERRQ(ierr);
00062   *flag = 1;
00063   for (row=0; row<size; row++) {
00064     ierr = MatGetRow(A,row,&ncol,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr);
00065     if (ncol>0) *flag = 0;
00066     ierr = MatRestoreRow(A,row,&ncol,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr);
00067     if (*flag==0) break;
00068   }
00069   PetscFunctionReturn(0);
00070 }

int MatRedistribute ( Mat *  A,
IS  splits 
)

Definition at line 582 of file icmk.c.

00583 {
00584   MPI_Comm comm; Mat B;
00585   int ntids,mytid,local_size,ierr;
00586 
00587   PetscFunctionBegin;
00588   ierr = PetscObjectGetComm((PetscObject)*A,&comm); CHKERRQ(ierr);
00589   MPI_Comm_size(comm,&ntids); MPI_Comm_rank(comm,&mytid);
00590   {
00591     /* find the improved distribution, return the new local_size */
00592     int *indices,chck;
00593     ierr = ISGetLocalSize(splits,&chck); CHKERRQ(ierr);
00594     if (chck!=ntids+1) SETERRQ(1,"Set of split points wrong length");
00595     ierr = ISGetIndices(splits,(const PetscInt**)&indices); CHKERRQ(ierr);
00596     local_size = indices[mytid+1]-indices[mytid];
00597     ierr = ISRestoreIndices(splits,(const PetscInt**)&indices); CHKERRQ(ierr);
00598   }
00599   {
00600     /* redistribute the matrix, unless the distribution is already as wanted */
00601     int perfect,gperfect,first,last,row;
00602     ierr = MatGetOwnershipRange(*A,&first,&last); CHKERRQ(ierr);
00603     perfect = local_size==(last-first);
00604     MPI_Allreduce(&perfect,&gperfect,1,MPI_INT,MPI_PROD,comm);
00605     if (gperfect) PetscFunctionReturn(0);
00606     ierr = MatCreateMPIAIJ /* optimise the d_ and o_ parameters ! */
00607       (comm,local_size,local_size,PETSC_DECIDE,PETSC_DECIDE,
00608        5,0,5,0,&B); CHKERRQ(ierr);
00609     for (row=first; row<last; row++) {
00610       int ncols; const int *cols; const PetscScalar *vals;
00611       ierr = MatGetRow(*A,row,&ncols,&cols,&vals); CHKERRQ(ierr);
00612       ierr = MatSetValues
00613         (B,1,&row,ncols,cols,vals,INSERT_VALUES); CHKERRQ(ierr);
00614       ierr = MatRestoreRow(*A,row,&ncols,&cols,&vals); CHKERRQ(ierr);
00615     }
00616     ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
00617     ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
00618   }
00619   ierr = MatDestroy(*A); CHKERRQ(ierr);
00620   *A = B;
00621   PetscFunctionReturn(0);
00622 }

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().

00309 {
00310   MPI_Comm comm;
00311   int M,N,ntids,mytid, first,last, isFirst,isLast,
00312     nsplits,Nsplits,*split_array,*Split_array,
00313     *Splits,S_top,
00314     ierr;
00315 
00316   PetscFunctionBegin;
00317 
00318   ierr = PetscObjectGetComm((PetscObject)A,&comm); CHKERRQ(ierr);
00319   MPI_Comm_size(comm,&ntids); MPI_Comm_rank(comm,&mytid);
00320   isFirst = (mytid==0); isLast = (mytid==ntids-1);
00321   ierr = MatGetSize(A,&M,&N); CHKERRQ(ierr);
00322   ierr = MatGetOwnershipRange(A,&first,&last); CHKERRQ(ierr);
00323 
00324   ierr = GetUpBiDiagSplits
00325     (A,first,last,isFirst,isLast,
00326      &split_array,&nsplits,&Split_array,&Nsplits); CHKERRQ(ierr);
00327   Split_array++; /* zero location is size */
00328 
00329   {
00330     IS *isses,*jsses; Mat *subs;
00331     int *splits,r_top,s_top, isplit,jsplit,*jsplits,*joffsets,
00332       ntest,nsave,*status;
00333 
00334 #define SET_J j = Split_array[jsplit];
00335 #define SET_TS  ts = PetscMax(PetscMin(bs/10,i),1);
00336 #define SET_TSS \
00337       tsr = PetscMin(ts,N-j); tsd = PetscMin(ts,M-i); tsc = PetscMin(ts,M-i);
00338 #define SET_LAST_COL \
00339       ierr = MatGetRow(A,i,&ncols,&cols,PETSC_NULL); CHKERRQ(ierr); \
00340       lastcol = cols[ncols-1]; \
00341       ierr = MatRestoreRow(A,i,&ncols,&cols,PETSC_NULL); CHKERRQ(ierr);
00342 #define IF_TOO_FAR_RIGHT_BREAK if (j-ts>lastcol) break;
00343     /* 
00344      * what is the number of splits for each of the local splits?
00345      * (we only need this for efficient alloation of the "status" array;
00346      * straightforward allocation of nsplits*Nsplits is far too much.)
00347      */
00348     ierr = PetscMalloc(nsplits*sizeof(int),&joffsets); CHKERRQ(ierr);
00349     ierr = PetscMemzero(joffsets,nsplits*sizeof(int)); CHKERRQ(ierr);
00350     ierr = PetscMalloc((nsplits+1)*sizeof(int),&jsplits); CHKERRQ(ierr);
00351     ierr = PetscMemzero(jsplits,(nsplits+1)*sizeof(int)); CHKERRQ(ierr);
00352     jsplits[0] = 0;
00353     for (isplit=0; isplit<nsplits; isplit++) {
00354       int i = split_array[isplit],ncols,lastcol,njsplits=0; const int *cols;
00355       if (!VERY_LAST_ROW) {
00356         SET_LAST_COL;
00357         for (jsplit=0; jsplit<Nsplits; jsplit++) {
00358           int j,bs,ts;
00359           SET_J; bs = j-i; SET_TS;
00360           if (bs<=0) {joffsets[isplit]++; continue;}
00361           IF_TOO_FAR_RIGHT_BREAK;
00362           njsplits++; /* count how many splits we actually use */
00363         }
00364       }
00365       jsplits[isplit+1] = jsplits[isplit]+njsplits;
00366     }
00367     ierr = PetscMalloc(jsplits[nsplits]*sizeof(int),&status); CHKERRQ(ierr);
00368     ierr = PetscMemzero(status,jsplits[nsplits]*sizeof(int)); CHKERRQ(ierr);
00369 #define LEFTDOWN 1
00370 #define RIGHTDOWN 2
00371 #define LEFTUP 4
00372 #define RIGHTUP 8
00373 #define CORNERUP 16
00374 #define STATUS(isp,jsp) status[jsplits[isp]+(jsp-joffsets[isplit])]
00375 #define SET_STATUS(isp,jsp,s) STATUS(isp,jsp) += s
00376     ntest = 3*jsplits[nsplits];
00377 #if TRACE_MSGS
00378     printf("total # tests: %d\n",ntest);
00379 #endif
00380     ierr = PetscMalloc(ntest*sizeof(IS),&isses); CHKERRQ(ierr);
00381     ierr = PetscMalloc(ntest*sizeof(IS),&jsses); CHKERRQ(ierr);
00382 #define OVERFLOW_DOWN i+tsd-1>=last
00383     /* why did that say >=last-1 ? */
00384 #define OVERFLOW_UP i-ts<first
00385     ntest = 0; nsave = 0;
00386     for (isplit=0; isplit<nsplits; isplit++) {
00387       int i=split_array[isplit],ncols,lastcol;
00388       const int *cols; PetscTruth flag;
00389       if (VERY_LAST_ROW) break;
00390 #if TRACE_MSGS > 1
00391       printf("from split %d=>%d: ",isplit,i);
00392 #endif
00393       SET_LAST_COL;
00394       for (jsplit=joffsets[isplit]; jsplit<Nsplits; jsplit++) {
00395         IS left,right,up,down,corner; int j,bs,ts,tsr,tsd,tsc;
00396         SET_J; bs = j-i; SET_TS; SET_TSS;
00397         IF_TOO_FAR_RIGHT_BREAK;
00398 #if TRACE_MSGS > 1
00399         printf("split %d=>%d, ",jsplit,j);
00400 #endif
00401         if ((OVERFLOW_DOWN) || (OVERFLOW_UP)) {
00402           ierr = ISCreateStride(MPI_COMM_SELF,tsr,j,1,&right); CHKERRQ(ierr);
00403           ierr = ISCreateStride(MPI_COMM_SELF,ts,j-ts,1,&left); CHKERRQ(ierr);}
00404         if (!VERY_LAST_ROW) {
00405           if (OVERFLOW_DOWN) {
00406             /*printf("saving down i=%d j=%d\n",i,j);*/
00407             ierr = ISCreateStride(MPI_COMM_SELF,tsd,i,1,&down); CHKERRQ(ierr);
00408             isses[nsave] = left; jsses[nsave] = down; nsave++;
00409             isses[nsave] = right; jsses[nsave] = down; nsave++;
00410           } else {
00411             ierr = MatSubMatIsNull
00412               (A,i,i+tsd-1,j-ts,j-1,&flag); CHKERRQ(ierr);
00413 #if TRACE_MSGS > 1
00414             printf("submat (%d,%d) x (%d,%d) : n=%d,  ",
00415                    i,i+tsd-1,j-ts,j-1,flag);
00416 #endif
00417             SET_STATUS(isplit,jsplit,flag*LEFTDOWN);
00418             ierr = MatSubMatIsNull
00419               (A,i,i+tsd-1,j,j+tsr-1,&flag); CHKERRQ(ierr);
00420 #if TRACE_MSGS > 1
00421             printf("submat (%d,%d) x (%d,%d) : n=%d,  ",
00422                    i,i+tsd-1,j,j+tsr-1,flag);
00423 #endif
00424             SET_STATUS(isplit,jsplit,flag*RIGHTDOWN);
00425           }
00426           ntest += 2;
00427         }
00428         if (!VERY_FIRST_ROW) {
00429           if (OVERFLOW_UP) {
00430             ierr = ISCreateStride(MPI_COMM_SELF,ts,i-ts,1,&up); CHKERRQ(ierr);
00431             ierr = ISCreateStride(MPI_COMM_SELF,tsc,i,1,&corner); CHKERRQ(ierr);
00432             /*printf("saving up i=%d j=%d\n",i,j);*/
00433             isses[nsave] = left; jsses[nsave] = up; nsave++;
00434             isses[nsave] = right; jsses[nsave] = up; nsave++;
00435             isses[nsave] = corner; jsses[nsave] = up; nsave++;
00436           } else {
00437             ierr = MatSubMatIsNull
00438               (A,i-ts,i-1,j-ts-1,j-1,&flag); CHKERRQ(ierr);
00439 #if TRACE_MSGS > 1
00440             printf("submat (%d,%d) x (%d,%d) : n=%d,  ",
00441                    i-ts,i-1,j-ts-1,j-1,flag);
00442 #endif
00443             SET_STATUS(isplit,jsplit,flag*LEFTUP);
00444             ierr = MatSubMatIsNull
00445               (A,i-ts,i-1,j,j+tsr-1,&flag); CHKERRQ(ierr);
00446             SET_STATUS(isplit,jsplit,flag*RIGHTUP);
00447             ierr = MatSubMatIsNull
00448               (A,i-ts,i-1,i,i+tsc-1,&flag); CHKERRQ(ierr);
00449             SET_STATUS(isplit,jsplit,flag*CORNERUP);
00450           }
00451           ntest += 3;
00452         }
00453       }
00454 #if TRACE_MSGS > 1
00455       printf("\n\n");
00456 #endif
00457     }
00458 
00459 #if TRACE_MSGS
00460     printf("saved %d straddling matrices out of %d tests\n",nsave,ntest);
00461 #endif
00462     ierr = MatGetSubMatrices
00463       (A,nsave,isses,jsses,MAT_INITIAL_MATRIX,&subs); CHKERRQ(ierr);
00464 
00465     ierr = PetscMalloc(SPLIT_BLOCK*ntest*sizeof(int),&splits); CHKERRQ(ierr);
00466     ntest = 0; r_top = 0; s_top = 0;
00467     for (isplit=0; isplit<nsplits; isplit++) {
00468       int i=split_array[isplit],ncols,lastcol,restart=1,split=1,c=0;
00469       const int *cols;
00470       if (VERY_LAST_ROW) break;
00471 #if TRACE_MSGS > 1
00472       printf("from split %d=>%d: ",isplit,i);
00473 #endif
00474       SET_LAST_COL;
00475       for (jsplit=joffsets[isplit]; jsplit<Nsplits; jsplit++) {
00476         int j,bs,ts,tsr,tsd,tsc,rank=0;
00477         SET_J; bs = j-i; SET_TS; SET_TSS;
00478 #if TRACE_MSGS > 1
00479         if (j-ts>lastcol) printf("%d:a ",j);
00480 #endif
00481         IF_TOO_FAR_RIGHT_BREAK;
00482         if (!VERY_LAST_ROW) {
00483           Mat leftdown,rightdown; int flag_leftdown,flag_rightdown;
00484           if (OVERFLOW_DOWN) {
00485             leftdown = subs[ntest++]; rightdown = subs[ntest++];
00486             ierr = MatIsNull(leftdown,&flag_leftdown); CHKERRQ(ierr);
00487             SET_STATUS(isplit,jsplit,flag_leftdown*LEFTDOWN);
00488             ierr = MatIsNull(rightdown,&flag_rightdown); CHKERRQ(ierr);
00489             SET_STATUS(isplit,jsplit,flag_rightdown*RIGHTDOWN);
00490           } else {
00491             flag_leftdown = STATUS(isplit,jsplit) & LEFTDOWN;
00492             flag_rightdown = STATUS(isplit,jsplit) & RIGHTDOWN;
00493           }
00494           restart = flag_leftdown && (!flag_rightdown || j==N);
00495 #if TRACE_MSGS > 1
00496           if (restart) printf("%d:dn:r, ",j); else printf("%d:dn:x, ",j);
00497 #endif
00498           if (restart) rank=1;
00499         }
00500         if (VERY_FIRST_ROW) {
00501           if (rank) rank=3;
00502         } else {
00503           Mat leftup,rightup,corner; int flag_leftup,flag_rightup,flag_corner;
00504           if (OVERFLOW_UP) {
00505             leftup = subs[ntest++]; rightup = subs[ntest++];
00506             ierr = MatIsNull(leftup,&flag_leftup); CHKERRQ(ierr);
00507             ierr = MatIsNull(rightup,&flag_rightup); CHKERRQ(ierr);
00508             corner = subs[ntest++];
00509             ierr = MatIsNull(corner,&flag_corner); CHKERRQ(ierr);
00510           } else {
00511             flag_leftup = STATUS(isplit,jsplit) & LEFTUP;
00512             flag_rightup = STATUS(isplit,jsplit) & RIGHTUP;
00513             flag_corner = STATUS(isplit,jsplit) & CORNERUP;
00514           }
00515           split = (flag_rightup || j==N) && !flag_leftup;
00516 #if TRACE_MSGS > 1
00517           if (split) printf("%d:up:s, ",j); else printf("%d:up:x, ",j);
00518 #endif
00519           if (rank && split) rank=2;
00520           if (rank==2 && flag_corner) rank=3;
00521         }
00522         if (rank>0) {
00523           c++;
00524           splits[s_top] = i; splits[s_top+1] = j; splits[s_top+2] = rank;
00525           s_top+=SPLIT_BLOCK;}
00526       }
00527       if (!c) { /* append next split as default case */
00528         splits[s_top] = i; splits[s_top+1] = split_array[isplit+1];
00529         splits[s_top+2] = 1; s_top+=SPLIT_BLOCK;}
00530 #if TRACE_MSGS > 1
00531       printf("\n");
00532 #endif
00533     }
00534     splits[s_top] = splits[s_top-SPLIT_BLOCK+1]; s_top++;
00535     splits[s_top++] = N; splits[s_top++] = 1;
00536 #if TRACE_MSGS > 1
00537     {
00538       int i;
00539       printf("local final splits: ");
00540       for (i=0; i<s_top; i+=SPLIT_BLOCK) {
00541         printf("[%d,%d],",splits[i],splits[i+1]);}
00542       printf("\n");}
00543 #endif
00544     ierr = MPIAllGatherIntV(splits,s_top,&Splits,&S_top,comm); CHKERRQ(ierr);
00545     Splits++;
00546 #ifdef TRACE_MSGS
00547     if (mytid==0) {
00548       int i;
00549       printf("global final splits: ");
00550       for (i=0; i<S_top; i+=SPLIT_BLOCK) {
00551         printf("[%d,%d:%d],",Splits[i],Splits[i+1],Splits[i+2]);
00552       }
00553       printf("\n");}
00554 #endif
00555     ierr = PetscFree(splits); CHKERRQ(ierr);
00556     ierr = PetscFree(joffsets); CHKERRQ(ierr);
00557     ierr = PetscFree(jsplits); CHKERRQ(ierr);
00558 
00559     for (isplit=0; isplit<nsave; isplit++) {
00560       ierr = ISDestroy(isses[isplit]); CHKERRQ(ierr);
00561       if (isplit>0 && (jsses[isplit]!=jsses[isplit-1])) {
00562         ierr = ISDestroy(jsses[isplit]); CHKERRQ(ierr);
00563       }
00564     }
00565     PetscFree(isses); PetscFree(jsses);
00566 
00567     ierr = MatDestroyMatrices(ntest,&subs); CHKERRQ(ierr); /*also frees subs*/
00568   }
00569   {
00570     IS isret;
00571     ierr = ChainSplits(N,Splits,S_top,np,&isret); CHKERRQ(ierr);
00572     *splitpoints = isret;
00573   }
00574   Splits--;
00575   ierr = PetscFree(Splits); CHKERRQ(ierr);
00576 
00577   PetscFunctionReturn(0);
00578 }

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().

00079 {
00080   int first,last,row,ncol,j,ierr; const int *cols;
00081   PetscFunctionBegin;
00082   ierr = MatGetOwnershipRange(A,&first,&last); CHKERRQ(ierr);
00083   if (i1<first || i2>=last) SETERRQ(1,"Requesting (part) non-local submat");
00084   *isnull = PETSC_TRUE;
00085   for (row=i1; row<=i2; row++) {
00086     ierr = MatGetRow(A,row,&ncol,&cols,PETSC_NULL); CHKERRQ(ierr);
00087     if (ncol>=0)
00088       if (cols[0]<=j2 && cols[ncol-1]>=j1)
00089         for (j=0; j<ncol; j++)
00090           if (cols[j]>=j1 && cols[j]<=j2) *isnull = PETSC_FALSE;
00091     ierr = MatRestoreRow(A,row,&ncol,&cols,PETSC_NULL); CHKERRQ(ierr);
00092     if (!*isnull) break;
00093   }
00094   PetscFunctionReturn(0);
00095 }

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().

00705 {
00706   Mat A = (Mat)prob;
00707   int v; PetscTruth has; int id; PetscErrorCode ierr;
00708   PetscFunctionBegin;
00709   ierr = GetDataID("icmk","n-splits",&id,&has); CHKERRQ(ierr);
00710   HASTOEXIST(has);
00711   ierr = PetscObjectComposedDataGetInt
00712     ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00713   if (*flg) rv->i = v;
00714   else {
00715     ierr = compute_icm_splits(A); CHKERRQ(ierr);
00716     ierr = PetscObjectComposedDataGetInt
00717       ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00718     if (*flg) rv->i = v;
00719   }
00720   
00721   PetscFunctionReturn(0);
00722 }

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().

00100 {
00101   int i;
00102   PetscFunctionBegin;
00103   printf("%s: ",txt); for (i=0; i<len; i++) printf("%d,",ar[i]); printf("\n");
00104   PetscFunctionReturn(0);
00105 }

PetscErrorCode RegisterICMKModules ( void   ) 

Definition at line 749 of file icmk.c.

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

00750 {
00751   PetscErrorCode ierr;
00752   PetscFunctionBegin;
00753 
00754   ierr = RegisterModule
00755     ("icmk","n-splits",ANALYSISINTEGER,&NSplits); CHKERRQ(ierr);
00756   ierr = RegisterModule
00757     ("icmk","splits",ANALYSISINTARRAY,&Splits); CHKERRQ(ierr);
00758 
00759   PetscFunctionReturn(0);
00760 }

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().

00728 {
00729   Mat A = (Mat)prob;
00730   int *v; PetscTruth has; int id; PetscErrorCode ierr;
00731   PetscFunctionBegin;
00732   ierr = GetDataID("icmk","splits",&id,&has); CHKERRQ(ierr);
00733   HASTOEXIST(has);
00734   ierr = PetscObjectComposedDataGetIntstar
00735     ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00736   if (*flg) rv->ii = v;
00737   else {
00738     ierr = compute_icm_splits(A); CHKERRQ(ierr);
00739     ierr = PetscObjectComposedDataGetIntstar
00740       ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00741     if (*flg) rv->ii = v;
00742   }
00743   
00744   PetscFunctionReturn(0);
00745 }

Here is the call graph for this function:

int VecRedistribute ( Vec *  V,
IS  splits 
)

Definition at line 623 of file icmk.c.

00624 {
00625   MPI_Comm comm; Vec W;
00626   int ntids,mytid,local_size,new_first,ierr;
00627 
00628   PetscFunctionBegin;
00629   ierr = PetscObjectGetComm((PetscObject)*V,&comm); CHKERRQ(ierr);
00630   MPI_Comm_size(comm,&ntids); MPI_Comm_rank(comm,&mytid);
00631   {
00632     /* find the new local_size */
00633     int *indices;
00634     ierr = ISGetIndices(splits,(const PetscInt**)&indices); CHKERRQ(ierr);
00635     local_size = indices[mytid+1]-indices[mytid]; new_first = indices[mytid];
00636     ierr = ISRestoreIndices(splits,(const PetscInt**)&indices); CHKERRQ(ierr);
00637   }
00638   {
00639     /* redistribute the matrix, unless the distribution is already as wanted */
00640     int perfect,gperfect,first,last;
00641     ierr = VecGetOwnershipRange(*V,&first,&last); CHKERRQ(ierr);
00642     perfect = local_size==(last-first);
00643     MPI_Allreduce(&perfect,&gperfect,1,MPI_INT,MPI_PROD,comm);
00644     if (gperfect) PetscFunctionReturn(0);
00645     ierr = VecCreateMPI
00646       (comm,local_size,PETSC_DECIDE,&W); CHKERRQ(ierr);
00647     {
00648       PetscScalar *elements; int *range,i;
00649       ierr = VecGetArray(*V,&elements); CHKERRQ(ierr);
00650       ierr = PetscMalloc(local_size*sizeof(int),&range); CHKERRQ(ierr);
00651       for (i=first; i<last; i++) range[i-first] = i;
00652       ierr = VecSetValues
00653         (W,last-first,range,elements,INSERT_VALUES); CHKERRQ(ierr);
00654       ierr = PetscFree(range); CHKERRQ(ierr);
00655       ierr = VecRestoreArray(*V,&elements); CHKERRQ(ierr);
00656     }
00657     ierr = VecAssemblyBegin(W); CHKERRQ(ierr);
00658     ierr = VecAssemblyEnd(W); CHKERRQ(ierr);
00659   }
00660   ierr = VecDestroy(*V); CHKERRQ(ierr);
00661   *V = W;
00662 
00663   PetscFunctionReturn(0);
00664 }


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().


Generated on Sun Oct 4 04:01:10 2009 for SALSA Analysis Modules by  doxygen 1.5.9