icmk.h File Reference

#include "petscmat.h"
#include "petscis.h"

Include dependency graph for icmk.h:

This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Functions

int MatSplitPoints (Mat A, int np, IS *splitpoints)


Function Documentation

int MatSplitPoints ( Mat  A,
int  np,
IS *  splitpoints 
)

Definition at line 308 of file icmk.c.

References ChainSplits(), CORNERUP, GetUpBiDiagSplits(), IF_TOO_FAR_RIGHT_BREAK, LEFTDOWN, LEFTUP, MatIsNull(), MatSubMatIsNull(), MPIAllGatherIntV(), OVERFLOW_DOWN, OVERFLOW_UP, RIGHTDOWN, RIGHTUP, SET_J, SET_LAST_COL, SET_STATUS, SET_TS, SET_TSS, SPLIT_BLOCK, Splits(), STATUS, VERY_FIRST_ROW, and VERY_LAST_ROW.

Referenced by compute_icm_splits().

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:


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