1 2 /* 3 Support for the parallel AIJ matrix vector multiply 4 */ 5 #include <../src/mat/impls/aij/mpi/mpiaij.h> 6 #include <petsc-private/isimpl.h> /* needed because accesses data structure of ISLocalToGlobalMapping directly */ 7 8 #undef __FUNCT__ 9 #define __FUNCT__ "MatSetUpMultiply_MPIAIJ" 10 PetscErrorCode MatSetUpMultiply_MPIAIJ(Mat mat) 11 { 12 Mat_MPIAIJ *aij = (Mat_MPIAIJ*)mat->data; 13 Mat_SeqAIJ *B = (Mat_SeqAIJ*)(aij->B->data); 14 PetscErrorCode ierr; 15 PetscInt i,j,*aj = B->j,ec = 0,*garray; 16 IS from,to; 17 Vec gvec; 18 #if defined(PETSC_USE_CTABLE) 19 PetscTable gid1_lid1; 20 PetscTablePosition tpos; 21 PetscInt gid,lid; 22 #else 23 PetscInt N = mat->cmap->N,*indices; 24 #endif 25 26 PetscFunctionBegin; 27 #if defined(PETSC_USE_CTABLE) 28 /* use a table */ 29 ierr = PetscTableCreate(aij->B->rmap->n,mat->cmap->N+1,&gid1_lid1);CHKERRQ(ierr); 30 for (i=0; i<aij->B->rmap->n; i++) { 31 for (j=0; j<B->ilen[i]; j++) { 32 PetscInt data,gid1 = aj[B->i[i] + j] + 1; 33 ierr = PetscTableFind(gid1_lid1,gid1,&data);CHKERRQ(ierr); 34 if (!data) { 35 /* one based table */ 36 ierr = PetscTableAdd(gid1_lid1,gid1,++ec,INSERT_VALUES);CHKERRQ(ierr); 37 } 38 } 39 } 40 /* form array of columns we need */ 41 ierr = PetscMalloc1((ec+1),&garray);CHKERRQ(ierr); 42 ierr = PetscTableGetHeadPosition(gid1_lid1,&tpos);CHKERRQ(ierr); 43 while (tpos) { 44 ierr = PetscTableGetNext(gid1_lid1,&tpos,&gid,&lid);CHKERRQ(ierr); 45 gid--; 46 lid--; 47 garray[lid] = gid; 48 } 49 ierr = PetscSortInt(ec,garray);CHKERRQ(ierr); /* sort, and rebuild */ 50 ierr = PetscTableRemoveAll(gid1_lid1);CHKERRQ(ierr); 51 for (i=0; i<ec; i++) { 52 ierr = PetscTableAdd(gid1_lid1,garray[i]+1,i+1,INSERT_VALUES);CHKERRQ(ierr); 53 } 54 /* compact out the extra columns in B */ 55 for (i=0; i<aij->B->rmap->n; i++) { 56 for (j=0; j<B->ilen[i]; j++) { 57 PetscInt gid1 = aj[B->i[i] + j] + 1; 58 ierr = PetscTableFind(gid1_lid1,gid1,&lid);CHKERRQ(ierr); 59 lid--; 60 aj[B->i[i] + j] = lid; 61 } 62 } 63 aij->B->cmap->n = aij->B->cmap->N = ec; 64 aij->B->cmap->bs = 1; 65 66 ierr = PetscLayoutSetUp((aij->B->cmap));CHKERRQ(ierr); 67 ierr = PetscTableDestroy(&gid1_lid1);CHKERRQ(ierr); 68 #else 69 /* Make an array as long as the number of columns */ 70 /* mark those columns that are in aij->B */ 71 ierr = PetscCalloc1(N+1,&indices);CHKERRQ(ierr); 72 for (i=0; i<aij->B->rmap->n; i++) { 73 for (j=0; j<B->ilen[i]; j++) { 74 if (!indices[aj[B->i[i] + j]]) ec++; 75 indices[aj[B->i[i] + j]] = 1; 76 } 77 } 78 79 /* form array of columns we need */ 80 ierr = PetscMalloc1((ec+1),&garray);CHKERRQ(ierr); 81 ec = 0; 82 for (i=0; i<N; i++) { 83 if (indices[i]) garray[ec++] = i; 84 } 85 86 /* make indices now point into garray */ 87 for (i=0; i<ec; i++) { 88 indices[garray[i]] = i; 89 } 90 91 /* compact out the extra columns in B */ 92 for (i=0; i<aij->B->rmap->n; i++) { 93 for (j=0; j<B->ilen[i]; j++) { 94 aj[B->i[i] + j] = indices[aj[B->i[i] + j]]; 95 } 96 } 97 aij->B->cmap->n = aij->B->cmap->N = ec; 98 aij->B->cmap->bs = 1; 99 100 ierr = PetscLayoutSetUp((aij->B->cmap));CHKERRQ(ierr); 101 ierr = PetscFree(indices);CHKERRQ(ierr); 102 #endif 103 /* create local vector that is used to scatter into */ 104 ierr = VecCreateSeq(PETSC_COMM_SELF,ec,&aij->lvec);CHKERRQ(ierr); 105 106 /* create two temporary Index sets for build scatter gather */ 107 ierr = ISCreateGeneral(((PetscObject)mat)->comm,ec,garray,PETSC_COPY_VALUES,&from);CHKERRQ(ierr); 108 109 ierr = ISCreateStride(PETSC_COMM_SELF,ec,0,1,&to);CHKERRQ(ierr); 110 111 /* create temporary global vector to generate scatter context */ 112 /* This does not allocate the array's memory so is efficient */ 113 ierr = VecCreateMPIWithArray(PetscObjectComm((PetscObject)mat),1,mat->cmap->n,mat->cmap->N,NULL,&gvec);CHKERRQ(ierr); 114 115 /* generate the scatter context */ 116 ierr = VecScatterCreate(gvec,from,aij->lvec,to,&aij->Mvctx);CHKERRQ(ierr); 117 ierr = PetscLogObjectParent((PetscObject)mat,(PetscObject)aij->Mvctx);CHKERRQ(ierr); 118 ierr = PetscLogObjectParent((PetscObject)mat,(PetscObject)aij->lvec);CHKERRQ(ierr); 119 ierr = PetscLogObjectParent((PetscObject)mat,(PetscObject)from);CHKERRQ(ierr); 120 ierr = PetscLogObjectParent((PetscObject)mat,(PetscObject)to);CHKERRQ(ierr); 121 122 aij->garray = garray; 123 124 ierr = PetscLogObjectMemory((PetscObject)mat,(ec+1)*sizeof(PetscInt));CHKERRQ(ierr); 125 ierr = ISDestroy(&from);CHKERRQ(ierr); 126 ierr = ISDestroy(&to);CHKERRQ(ierr); 127 ierr = VecDestroy(&gvec);CHKERRQ(ierr); 128 PetscFunctionReturn(0); 129 } 130 131 132 #undef __FUNCT__ 133 #define __FUNCT__ "MatDisAssemble_MPIAIJ" 134 /* 135 Takes the local part of an already assembled MPIAIJ matrix 136 and disassembles it. This is to allow new nonzeros into the matrix 137 that require more communication in the matrix vector multiply. 138 Thus certain data-structures must be rebuilt. 139 140 Kind of slow! But that's what application programmers get when 141 they are sloppy. 142 */ 143 PetscErrorCode MatDisAssemble_MPIAIJ(Mat A) 144 { 145 Mat_MPIAIJ *aij = (Mat_MPIAIJ*)A->data; 146 Mat B = aij->B,Bnew; 147 Mat_SeqAIJ *Baij = (Mat_SeqAIJ*)B->data; 148 PetscErrorCode ierr; 149 PetscInt i,j,m = B->rmap->n,n = A->cmap->N,col,ct = 0,*garray = aij->garray,*nz,ec; 150 PetscScalar v; 151 152 PetscFunctionBegin; 153 /* free stuff related to matrix-vec multiply */ 154 ierr = VecGetSize(aij->lvec,&ec);CHKERRQ(ierr); /* needed for PetscLogObjectMemory below */ 155 ierr = VecDestroy(&aij->lvec);CHKERRQ(ierr); 156 ierr = VecScatterDestroy(&aij->Mvctx);CHKERRQ(ierr); 157 if (aij->colmap) { 158 #if defined(PETSC_USE_CTABLE) 159 ierr = PetscTableDestroy(&aij->colmap);CHKERRQ(ierr); 160 #else 161 ierr = PetscFree(aij->colmap);CHKERRQ(ierr); 162 ierr = PetscLogObjectMemory((PetscObject)A,-aij->B->cmap->n*sizeof(PetscInt));CHKERRQ(ierr); 163 #endif 164 } 165 166 /* make sure that B is assembled so we can access its values */ 167 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 168 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 169 170 /* invent new B and copy stuff over */ 171 ierr = PetscMalloc1((m+1),&nz);CHKERRQ(ierr); 172 for (i=0; i<m; i++) { 173 nz[i] = Baij->i[i+1] - Baij->i[i]; 174 } 175 ierr = MatCreate(PETSC_COMM_SELF,&Bnew);CHKERRQ(ierr); 176 ierr = MatSetSizes(Bnew,m,n,m,n);CHKERRQ(ierr); 177 ierr = MatSetBlockSizesFromMats(Bnew,A,A);CHKERRQ(ierr); 178 ierr = MatSetType(Bnew,((PetscObject)B)->type_name);CHKERRQ(ierr); 179 ierr = MatSeqAIJSetPreallocation(Bnew,0,nz);CHKERRQ(ierr); 180 181 ((Mat_SeqAIJ*)Bnew->data)->nonew = Baij->nonew; /* Inherit insertion error options. */ 182 Bnew->nonzerostate = B->nonzerostate; 183 184 ierr = PetscFree(nz);CHKERRQ(ierr); 185 for (i=0; i<m; i++) { 186 for (j=Baij->i[i]; j<Baij->i[i+1]; j++) { 187 col = garray[Baij->j[ct]]; 188 v = Baij->a[ct++]; 189 ierr = MatSetValues(Bnew,1,&i,1,&col,&v,B->insertmode);CHKERRQ(ierr); 190 } 191 } 192 ierr = PetscFree(aij->garray);CHKERRQ(ierr); 193 ierr = PetscLogObjectMemory((PetscObject)A,-ec*sizeof(PetscInt));CHKERRQ(ierr); 194 ierr = MatDestroy(&B);CHKERRQ(ierr); 195 ierr = PetscLogObjectParent((PetscObject)A,(PetscObject)Bnew);CHKERRQ(ierr); 196 197 aij->B = Bnew; 198 A->was_assembled = PETSC_FALSE; 199 PetscFunctionReturn(0); 200 } 201 202 /* ugly stuff added for Glenn someday we should fix this up */ 203 204 static PetscInt *auglyrmapd = 0,*auglyrmapo = 0; /* mapping from the local ordering to the "diagonal" and "off-diagonal" parts of the local matrix */ 205 static Vec auglydd = 0,auglyoo = 0; /* work vectors used to scale the two parts of the local matrix */ 206 207 208 #undef __FUNCT__ 209 #define __FUNCT__ "MatMPIAIJDiagonalScaleLocalSetUp" 210 PetscErrorCode MatMPIAIJDiagonalScaleLocalSetUp(Mat inA,Vec scale) 211 { 212 Mat_MPIAIJ *ina = (Mat_MPIAIJ*) inA->data; /*access private part of matrix */ 213 PetscErrorCode ierr; 214 PetscInt i,n,nt,cstart,cend,no,*garray = ina->garray,*lindices; 215 PetscInt *r_rmapd,*r_rmapo; 216 217 PetscFunctionBegin; 218 ierr = MatGetOwnershipRange(inA,&cstart,&cend);CHKERRQ(ierr); 219 ierr = MatGetSize(ina->A,NULL,&n);CHKERRQ(ierr); 220 ierr = PetscCalloc1(inA->rmap->mapping->n+1,&r_rmapd);CHKERRQ(ierr); 221 nt = 0; 222 for (i=0; i<inA->rmap->mapping->n; i++) { 223 if (inA->rmap->mapping->indices[i] >= cstart && inA->rmap->mapping->indices[i] < cend) { 224 nt++; 225 r_rmapd[i] = inA->rmap->mapping->indices[i] + 1; 226 } 227 } 228 if (nt != n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Hmm nt %D n %D",nt,n); 229 ierr = PetscMalloc1((n+1),&auglyrmapd);CHKERRQ(ierr); 230 for (i=0; i<inA->rmap->mapping->n; i++) { 231 if (r_rmapd[i]) { 232 auglyrmapd[(r_rmapd[i]-1)-cstart] = i; 233 } 234 } 235 ierr = PetscFree(r_rmapd);CHKERRQ(ierr); 236 ierr = VecCreateSeq(PETSC_COMM_SELF,n,&auglydd);CHKERRQ(ierr); 237 238 ierr = PetscCalloc1(inA->cmap->N+1,&lindices);CHKERRQ(ierr); 239 for (i=0; i<ina->B->cmap->n; i++) { 240 lindices[garray[i]] = i+1; 241 } 242 no = inA->rmap->mapping->n - nt; 243 ierr = PetscCalloc1(inA->rmap->mapping->n+1,&r_rmapo);CHKERRQ(ierr); 244 nt = 0; 245 for (i=0; i<inA->rmap->mapping->n; i++) { 246 if (lindices[inA->rmap->mapping->indices[i]]) { 247 nt++; 248 r_rmapo[i] = lindices[inA->rmap->mapping->indices[i]]; 249 } 250 } 251 if (nt > no) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Hmm nt %D no %D",nt,n); 252 ierr = PetscFree(lindices);CHKERRQ(ierr); 253 ierr = PetscMalloc1((nt+1),&auglyrmapo);CHKERRQ(ierr); 254 for (i=0; i<inA->rmap->mapping->n; i++) { 255 if (r_rmapo[i]) { 256 auglyrmapo[(r_rmapo[i]-1)] = i; 257 } 258 } 259 ierr = PetscFree(r_rmapo);CHKERRQ(ierr); 260 ierr = VecCreateSeq(PETSC_COMM_SELF,nt,&auglyoo);CHKERRQ(ierr); 261 PetscFunctionReturn(0); 262 } 263 264 #undef __FUNCT__ 265 #define __FUNCT__ "MatMPIAIJDiagonalScaleLocal" 266 PetscErrorCode MatMPIAIJDiagonalScaleLocal(Mat A,Vec scale) 267 { 268 /* This routine should really be abandoned as it duplicates MatDiagonalScaleLocal */ 269 PetscErrorCode ierr; 270 271 PetscFunctionBegin; 272 ierr = PetscTryMethod(A,"MatDiagonalScaleLocal_C",(Mat,Vec),(A,scale));CHKERRQ(ierr); 273 PetscFunctionReturn(0); 274 } 275 276 #undef __FUNCT__ 277 #define __FUNCT__ "MatDiagonalScaleLocal_MPIAIJ" 278 PetscErrorCode MatDiagonalScaleLocal_MPIAIJ(Mat A,Vec scale) 279 { 280 Mat_MPIAIJ *a = (Mat_MPIAIJ*) A->data; /*access private part of matrix */ 281 PetscErrorCode ierr; 282 PetscInt n,i; 283 PetscScalar *d,*o,*s; 284 285 PetscFunctionBegin; 286 if (!auglyrmapd) { 287 ierr = MatMPIAIJDiagonalScaleLocalSetUp(A,scale);CHKERRQ(ierr); 288 } 289 290 ierr = VecGetArray(scale,&s);CHKERRQ(ierr); 291 292 ierr = VecGetLocalSize(auglydd,&n);CHKERRQ(ierr); 293 ierr = VecGetArray(auglydd,&d);CHKERRQ(ierr); 294 for (i=0; i<n; i++) { 295 d[i] = s[auglyrmapd[i]]; /* copy "diagonal" (true local) portion of scale into dd vector */ 296 } 297 ierr = VecRestoreArray(auglydd,&d);CHKERRQ(ierr); 298 /* column scale "diagonal" portion of local matrix */ 299 ierr = MatDiagonalScale(a->A,NULL,auglydd);CHKERRQ(ierr); 300 301 ierr = VecGetLocalSize(auglyoo,&n);CHKERRQ(ierr); 302 ierr = VecGetArray(auglyoo,&o);CHKERRQ(ierr); 303 for (i=0; i<n; i++) { 304 o[i] = s[auglyrmapo[i]]; /* copy "off-diagonal" portion of scale into oo vector */ 305 } 306 ierr = VecRestoreArray(scale,&s);CHKERRQ(ierr); 307 ierr = VecRestoreArray(auglyoo,&o);CHKERRQ(ierr); 308 /* column scale "off-diagonal" portion of local matrix */ 309 ierr = MatDiagonalScale(a->B,NULL,auglyoo);CHKERRQ(ierr); 310 PetscFunctionReturn(0); 311 } 312 313 314