SALSA Analysis Modules
|
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"
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 |
Functions for computing the ICMK structure of a matrix.
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.
Activate this module with:
PetscErrorCode RegisterICMKModules();
Compute these elements with
ComputeQuantity("icmk","element",A,&res,&flg);
Available elements are:
@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 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 |
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().
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); }
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); }
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); }
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); }
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); }
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); }
static void iSpaces | ( | int | len | ) | [static] |
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); }
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); }
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); }
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); }
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); }
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().