1 static int MatIncreaseOverlap_MPIAIJ(Mat, int id_msx, IS *is, int ov); 2 { 3 Mat_MPIAIJ *a = (Mat_MPIAIJ *) A->data; 4 int shift, **table; **data , *xtable, *xdata;, **idx, *n, *w1; 5 int *w2, *w3, *w4, *rtable ; 6 /* assume overlap = 1 */ 7 8 size = a->size; 9 rank = a->rank; 10 m = a->M; 11 idx = PetscMalloc((is_max)*sizeof(int *)); 12 n = PetscMalloc((is_max)*sizeof(int )); 13 rtable = PetscMalloc((m+1)sizeof(int )); /* Hash table for maping row ->proc */ 14 15 for ( i=0 ; i<is_max ; ++i) { 16 ierr = ISGetIndices(is[i],&idx[i]); CHKERRQ(ierr); 17 ierr = ISGetLocalSize(is[i],&n[i]); CHKERRQ(ierr); 18 } 19 20 /* Create hash table for the mapping :row -> proc*/ 21 for( i=0 j=0; i< size; ++i) { 22 for (; j <rowners[i+1]; ++j) { 23 rtable[j] = i; 24 } 25 } 26 27 /* evaluate communication - mesg to who, length of mesg, and buffer space 28 required. Based on this, buffers are allocated, and data copied into them*/ 29 w1 = PetscMalloc((size)*2*sizeof(int )); /* foreach proc mesg size */ 30 w2 = w1 + size; /* if w2[i] marked, then a message to proc i*/ 31 w3 = w2 + size; /* no of IS that needs to be sent to proc i */ 32 w4 = w3 + size; /* temp work space used in determining w1, w2, w3 */ 33 PetscMemzero(w1,(size)*3*sizeof(int)); /* initialise work vector*/ 34 for ( i=0; i<is_max ; ++i) { 35 PetscMemzero(w4,(size)*sizeof(int)); /* initialise work vector*/ 36 for ( j =0 ; j < n[i] ; ++j) { 37 row = idx[i][j]; 38 proc = rtable[row]; 39 w4[proc]++; 40 } 41 for( j = 0; j < size; ++j){ 42 if( w4[j] ) { w1[j] += w4[j]; w3[j] += 1;} 43 } 44 } 45 46 mct = 0; /* no of outgoing messages */ 47 msz = 0; /* total mesg length (for all proc */ 48 w1[rank] = 0; /* no mesg sent to intself */ 49 for (j =0; j < size ; ++j) { 50 if (w1[j]) { w2[j] = 1; mct++;} /* there exists a message to proc i */ 51 } 52 /* Each message would have a header = 1 + no of IS + data */ 53 for (j = 0; j<size ; ++j) { 54 w1[j] += w2[j] + w3[j]; /* 0+0 or 1 + XX */ 55 msz += w1[j]; 56 } 57 58 /* MPI_Allreduce */ 59 60 /* Allocate Memory for outgoing messages */ 61 outdat = PetscMalloc( 2*size*sizeof(int*)); 62 outdat[0] = PetscMalloc(msz *sizeof (int)); 63 ptr = outdat +size; /* Pointers to the data in outgoing buffers */ 64 ctr = PetscMalloc( size*sizeof(int)); 65 66 for (i = 1; i < size ; ++i) { 67 if ( w1[i]) { outdat[i] = w1[ i-1];} 68 else { outdat[i] = PETSC_NULL; } 69 } 70 71 /* Form the outgoing messages */ 72 /*plug in the headers*/ 73 for ( i=0 ; i<size ; ++i) { 74 if (w3[i]) { 75 outdat[i][0] = 0; /* Updated later*/ 76 PetscMemzero(outdat[i]+1, w3[i]*sizeof(int)); 77 ptr[i] = outdat[i] + w3[i]; 78 } 79 else { ptr[i] = PETSC_NULL; } 80 } 81 82 /* ??? How to handle the local data computations? */ 83 for ( i=0 ; i<is_max ; ++i) { 84 PetscMemzero(ctr,size*sizeof(int)); 85 for( j=0; j<n[i]; ++j) { /* parse the indices of each IS */ 86 row = idx[i][j]; 87 proc = rtable[row]; 88 if (proc != rank) { 89 ++count[proc]; 90 *ptr[proc] = row; 91 ++ptr[proc]; 92 } 93 } 94 /* Update the headers*/ 95 for( i = 0; i<size; ++i) { 96 if (ctr[i]) { 97 outdat[i][0] ++; 98 outdat[i][ outdat[i][0]] = ctr[i]; 99 } 100 } 101 } 102 103 /* Check Validity */ 104 for ( i=0 ; i<size ; ++i) { 105 if (w3[i] != outdat[i][0]) {SETERRQ(1,"MatIncreaseOverlap_MPIAIJ: Blew it!\n"); } 106 } 107 108 109 110 111 table = PetscMalloc((ismax)*sizeof(int *)); /* Hash table for each IS */ 112 data = PetscMalloc((is_max)*sizeof(int *)); /* The overlap sol stored here */ 113 table[0] = PetscMalloc((m+1)*(is_max+1)*2*sizeof(int)); 114 data [0] = table[0] + (m+1)*(is_max+1); 115 116 for(i = i; i<is_max ; ++i) { 117 table[i] = table[0] + (m+1)*i; 118 data[i] = table[0] + (m+1)*i; 119 } 120 xdata = table[0] + (m+1)*i; 121 PetscMemzero(*table,(m+1)*(is_max+1)*sizeof(int)); 122 123 /* whew! already done? check again :) */ 124 return 0 125 } 126