1 #ifndef lint 2 static char vcid[] = "$Id: mpiov.c,v 1.2 1996/01/23 00:44:57 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 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 for (j =0; j < size ; ++j) { 56 if (w1[j]) { w2[j] = 1; mct++;} /* there exists a message to proc i */ 57 } 58 /* Each message would have a header = 1 + no of IS + data */ 59 for (j = 0; j<size ; ++j) { 60 w1[j] += w2[j] + w3[j]; /* 0+0 or 1 + XX */ 61 msz += w1[j]; 62 } 63 64 /* MPI_Allreduce */ 65 66 /* Allocate Memory for outgoing messages */ 67 outdat = (int **)PetscMalloc( 2*size*sizeof(int*)); 68 outdat[0] = (int *)PetscMalloc(msz *sizeof (int)); 69 ptr = outdat +size; /* Pointers to the data in outgoing buffers */ 70 ctr = (int *)PetscMalloc( size*sizeof(int)); 71 72 for (i = 1; i < size ; ++i) { 73 if ( w1[i]) { outdat[i] = outdat[i-1] + w1[ i-1];} 74 else { outdat[i] = PETSC_NULL; } 75 } 76 77 /* Form the outgoing messages */ 78 /*plug in the headers*/ 79 for ( i=0 ; i<size ; ++i) { 80 if (w3[i]) { 81 outdat[i][0] = 0; /* Updated later*/ 82 PetscMemzero(outdat[i]+1, w3[i]*sizeof(int)); 83 ptr[i] = outdat[i] + w3[i]; 84 } 85 else { ptr[i] = PETSC_NULL; } 86 } 87 88 /* ??? How to handle the local data computations? */ 89 for ( i=0 ; i<is_max ; ++i) { 90 PetscMemzero(ctr,size*sizeof(int)); 91 for( j=0; j<n[i]; ++j) { /* parse the indices of each IS */ 92 row = idx[i][j]; 93 proc = rtable[row]; 94 if (proc != rank) { 95 ++ctr[proc]; 96 *ptr[proc] = row; 97 ++ptr[proc]; 98 } 99 } 100 /* Update the headers*/ 101 for( i = 0; i<size; ++i) { 102 if (ctr[i]) { 103 outdat[i][0] ++; 104 outdat[i][ outdat[i][0]] = ctr[i]; 105 } 106 } 107 } 108 109 /* Check Validity */ 110 for ( i=0 ; i<size ; ++i) { 111 if (w3[i] != outdat[i][0]) {SETERRQ(1,"MatIncreaseOverlap_MPIAIJ: Blew it!\n"); } 112 } 113 114 115 116 117 /* table = PetscMalloc((ismax)*sizeof(int *)); 118 data = PetscMalloc((is_max)*sizeof(int *)); 119 table[0] = PetscMalloc((m+1)*(is_max+1)*2*sizeof(int)); 120 data [0] = table[0] + (m+1)*(is_max+1); 121 122 for(i = i; i<is_max ; ++i) { 123 table[i] = table[0] + (m+1)*i; 124 data[i] = table[0] + (m+1)*i; 125 } 126 xdata = table[0] + (m+1)*i; 127 PetscMemzero(*table,(m+1)*(is_max+1)*sizeof(int));*/ 128 129 /* whew! already done? check again :) */ 130 PetscFree(idx); 131 PetscFree(n); 132 PetscFree(rtable); 133 PetscFree(w1); 134 PetscFree(outdat[0]); 135 PetscFree( outdat); 136 PetscFree(ctr); 137 return 0; 138 } 139 140