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