1 2 /* 3 Defines matrix-matrix product routines for pairs of MPIAIJ matrices 4 C = A * B 5 */ 6 #include <../src/mat/impls/aij/seq/aij.h> /*I "petscmat.h" I*/ 7 #include <../src/mat/utils/freespace.h> 8 #include <../src/mat/impls/aij/mpi/mpiaij.h> 9 #include <petscbt.h> 10 #include <../src/mat/impls/dense/mpi/mpidense.h> 11 12 #undef __FUNCT__ 13 #define __FUNCT__ "MatMatMult_MPIAIJ_MPIAIJ" 14 PetscErrorCode MatMatMult_MPIAIJ_MPIAIJ(Mat A,Mat B,MatReuse scall,PetscReal fill, Mat *C) 15 { 16 PetscErrorCode ierr; 17 18 PetscFunctionBegin; 19 if (scall == MAT_INITIAL_MATRIX){ 20 ierr = PetscLogEventBegin(MAT_MatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 21 ierr = MatMatMultSymbolic_MPIAIJ_MPIAIJ(A,B,fill,C);CHKERRQ(ierr); 22 ierr = PetscLogEventEnd(MAT_MatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 23 } 24 25 ierr = PetscLogEventBegin(MAT_MatMultNumeric,A,B,0,0);CHKERRQ(ierr); 26 ierr = (*(*C)->ops->matmultnumeric)(A,B,*C);CHKERRQ(ierr); 27 ierr = PetscLogEventEnd(MAT_MatMultNumeric,A,B,0,0);CHKERRQ(ierr); 28 PetscFunctionReturn(0); 29 } 30 31 #undef __FUNCT__ 32 #define __FUNCT__ "MatDestroy_MPIAIJ_MatMatMult" 33 PetscErrorCode MatDestroy_MPIAIJ_MatMatMult(Mat A) 34 { 35 PetscErrorCode ierr; 36 Mat_MPIAIJ *a=(Mat_MPIAIJ*)A->data; 37 Mat_PtAPMPI *ptap=a->ptap; 38 39 PetscFunctionBegin; 40 ierr = PetscFree2(ptap->startsj_s,ptap->startsj_r);CHKERRQ(ierr); 41 ierr = PetscFree(ptap->bufa);CHKERRQ(ierr); 42 ierr = MatDestroy(&ptap->P_loc);CHKERRQ(ierr); 43 ierr = MatDestroy(&ptap->P_oth);CHKERRQ(ierr); 44 ierr = MatDestroy(&ptap->Pt);CHKERRQ(ierr); 45 ierr = PetscFree(ptap->api);CHKERRQ(ierr); 46 ierr = PetscFree(ptap->apj);CHKERRQ(ierr); 47 ierr = PetscFree(ptap->apa);CHKERRQ(ierr); 48 ierr = ptap->destroy(A);CHKERRQ(ierr); 49 ierr = PetscFree(ptap);CHKERRQ(ierr); 50 PetscFunctionReturn(0); 51 } 52 53 #undef __FUNCT__ 54 #define __FUNCT__ "MatDuplicate_MPIAIJ_MatMatMult" 55 PetscErrorCode MatDuplicate_MPIAIJ_MatMatMult(Mat A, MatDuplicateOption op, Mat *M) 56 { 57 PetscErrorCode ierr; 58 Mat_MPIAIJ *a=(Mat_MPIAIJ*)A->data; 59 Mat_PtAPMPI *ptap=a->ptap; 60 61 PetscFunctionBegin; 62 ierr = (*ptap->duplicate)(A,op,M);CHKERRQ(ierr); 63 (*M)->ops->destroy = ptap->destroy; /* = MatDestroy_MPIAIJ, *M doesn't duplicate A's special structure! */ 64 (*M)->ops->duplicate = ptap->duplicate; /* = MatDuplicate_MPIAIJ */ 65 PetscFunctionReturn(0); 66 } 67 68 #undef __FUNCT__ 69 #define __FUNCT__ "MatMatMultNumeric_MPIAIJ_MPIAIJ" 70 PetscErrorCode MatMatMultNumeric_MPIAIJ_MPIAIJ(Mat A,Mat P,Mat C) 71 { 72 PetscErrorCode ierr; 73 Mat_MPIAIJ *a=(Mat_MPIAIJ*)A->data,*c=(Mat_MPIAIJ*)C->data; 74 Mat_SeqAIJ *ad=(Mat_SeqAIJ*)(a->A)->data,*ao=(Mat_SeqAIJ*)(a->B)->data; 75 Mat_SeqAIJ *cd=(Mat_SeqAIJ*)(c->A)->data,*co=(Mat_SeqAIJ*)(c->B)->data; 76 PetscInt *adi=ad->i,*adj,*aoi=ao->i,*aoj; 77 PetscScalar *ada,*aoa,*cda=cd->a,*coa=co->a; 78 Mat_SeqAIJ *p_loc,*p_oth; 79 PetscInt *pi_loc,*pj_loc,*pi_oth,*pj_oth,*pj; 80 PetscScalar *pa_loc,*pa_oth,*pa,*apa,valtmp,*ca; 81 PetscInt cm=C->rmap->n,anz,pnz; 82 Mat_PtAPMPI *ptap=c->ptap; 83 PetscInt *api,*apj,*apJ,i,j,k,row; 84 PetscInt cstart=C->cmap->rstart; 85 PetscInt cdnz,conz,k0,k1; 86 87 PetscFunctionBegin; 88 /* 1) get P_oth = ptap->P_oth and P_loc = ptap->P_loc */ 89 /*-----------------------------------------------------*/ 90 /* update numerical values of P_oth and P_loc */ 91 ierr = MatGetBrowsOfAoCols_MPIAIJ(A,P,MAT_REUSE_MATRIX,&ptap->startsj_s,&ptap->startsj_r,&ptap->bufa,&ptap->P_oth);CHKERRQ(ierr); 92 ierr = MatMPIAIJGetLocalMat(P,MAT_REUSE_MATRIX,&ptap->P_loc);CHKERRQ(ierr); 93 94 /* 2) compute numeric C_loc = A_loc*P = Ad*P_loc + Ao*P_oth */ 95 /*----------------------------------------------------------*/ 96 /* get data from symbolic products */ 97 p_loc = (Mat_SeqAIJ*)(ptap->P_loc)->data; 98 p_oth = (Mat_SeqAIJ*)(ptap->P_oth)->data; 99 pi_loc=p_loc->i; pj_loc=p_loc->j; pa_loc=p_loc->a; 100 pi_oth=p_oth->i; pj_oth=p_oth->j; pa_oth=p_oth->a; 101 102 /* get apa for storing dense row A[i,:]*P */ 103 apa = ptap->apa; 104 105 api = ptap->api; 106 apj = ptap->apj; 107 for (i=0; i<cm; i++) { 108 /* diagonal portion of A */ 109 anz = adi[i+1] - adi[i]; 110 adj = ad->j + adi[i]; 111 ada = ad->a + adi[i]; 112 for (j=0; j<anz; j++) { 113 row = adj[j]; 114 pnz = pi_loc[row+1] - pi_loc[row]; 115 pj = pj_loc + pi_loc[row]; 116 pa = pa_loc + pi_loc[row]; 117 118 /* perform dense axpy */ 119 valtmp = ada[j]; 120 for (k=0; k<pnz; k++){ 121 apa[pj[k]] += valtmp*pa[k]; 122 } 123 ierr = PetscLogFlops(2.0*pnz);CHKERRQ(ierr); 124 } 125 126 /* off-diagonal portion of A */ 127 anz = aoi[i+1] - aoi[i]; 128 aoj = ao->j + aoi[i]; 129 aoa = ao->a + aoi[i]; 130 for (j=0; j<anz; j++) { 131 row = aoj[j]; 132 pnz = pi_oth[row+1] - pi_oth[row]; 133 pj = pj_oth + pi_oth[row]; 134 pa = pa_oth + pi_oth[row]; 135 136 /* perform dense axpy */ 137 valtmp = aoa[j]; 138 for (k=0; k<pnz; k++){ 139 apa[pj[k]] += valtmp*pa[k]; 140 } 141 ierr = PetscLogFlops(2.0*pnz);CHKERRQ(ierr); 142 } 143 144 /* set values in C */ 145 apJ = apj + api[i]; 146 cdnz = cd->i[i+1] - cd->i[i]; 147 conz = co->i[i+1] - co->i[i]; 148 149 /* 1st off-diagoanl part of C */ 150 ca = coa + co->i[i]; 151 k = 0; 152 for (k0=0; k0<conz; k0++){ 153 if (apJ[k] >= cstart) break; 154 ca[k0] = apa[apJ[k]]; 155 apa[apJ[k]] = 0.0; 156 k++; 157 } 158 159 /* diagonal part of C */ 160 ca = cda + cd->i[i]; 161 for (k1=0; k1<cdnz; k1++){ 162 ca[k1] = apa[apJ[k]]; 163 apa[apJ[k]] = 0.0; 164 k++; 165 } 166 167 /* 2nd off-diagoanl part of C */ 168 ca = coa + co->i[i]; 169 for (; k0<conz; k0++){ 170 ca[k0] = apa[apJ[k]]; 171 apa[apJ[k]] = 0.0; 172 k++; 173 } 174 } 175 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 176 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 177 PetscFunctionReturn(0); 178 } 179 180 #undef __FUNCT__ 181 #define __FUNCT__ "MatMatMultSymbolic_MPIAIJ_MPIAIJ" 182 PetscErrorCode MatMatMultSymbolic_MPIAIJ_MPIAIJ(Mat A,Mat P,PetscReal fill,Mat *C) 183 { 184 PetscErrorCode ierr; 185 MPI_Comm comm=((PetscObject)A)->comm; 186 Mat Cmpi; 187 Mat_PtAPMPI *ptap; 188 PetscFreeSpaceList free_space=PETSC_NULL,current_space=PETSC_NULL; 189 Mat_MPIAIJ *a=(Mat_MPIAIJ*)A->data,*c; 190 Mat_SeqAIJ *ad=(Mat_SeqAIJ*)(a->A)->data,*ao=(Mat_SeqAIJ*)(a->B)->data,*p_loc,*p_oth; 191 PetscInt *pi_loc,*pj_loc,*pi_oth,*pj_oth,*dnz,*onz; 192 PetscInt *adi=ad->i,*adj=ad->j,*aoi=ao->i,*aoj=ao->j,rstart=A->rmap->rstart; 193 PetscInt *lnk,i,pnz,row,*api,*apj,*Jptr,apnz,nspacedouble=0,j,nzi; 194 PetscInt am=A->rmap->n,pN=P->cmap->N,pn=P->cmap->n,pm=P->rmap->n; 195 PetscBT lnkbt; 196 PetscScalar *apa; 197 PetscReal afill; 198 PetscBool scalable=PETSC_TRUE; 199 PetscInt nlnk_max,armax,prmax; 200 201 PetscFunctionBegin; 202 if (A->cmap->rstart != P->rmap->rstart || A->cmap->rend != P->rmap->rend){ 203 SETERRQ4(comm,PETSC_ERR_ARG_SIZ,"Matrix local dimensions are incompatible, (%D, %D) != (%D,%D)",A->cmap->rstart,A->cmap->rend,P->rmap->rstart,P->rmap->rend); 204 } 205 206 ierr = PetscObjectOptionsBegin((PetscObject)A);CHKERRQ(ierr); 207 ierr = PetscOptionsBool("-matmatmult_scalable","Use a scalable but slower C=A*B","",scalable,&scalable,PETSC_NULL);CHKERRQ(ierr); 208 if (scalable){ 209 ierr = MatMatMultSymbolic_MPIAIJ_MPIAIJ_Scalable(A,P,fill,C);;CHKERRQ(ierr); 210 PetscFunctionReturn(0); 211 } 212 ierr = PetscOptionsEnd();CHKERRQ(ierr); 213 214 /* create struct Mat_PtAPMPI and attached it to C later */ 215 ierr = PetscNew(Mat_PtAPMPI,&ptap);CHKERRQ(ierr); 216 217 /* get P_oth by taking rows of P (= non-zero cols of local A) from other processors */ 218 ierr = MatGetBrowsOfAoCols_MPIAIJ(A,P,MAT_INITIAL_MATRIX,&ptap->startsj_s,&ptap->startsj_r,&ptap->bufa,&ptap->P_oth);CHKERRQ(ierr); 219 220 /* get P_loc by taking all local rows of P */ 221 ierr = MatMPIAIJGetLocalMat(P,MAT_INITIAL_MATRIX,&ptap->P_loc);CHKERRQ(ierr); 222 223 p_loc = (Mat_SeqAIJ*)(ptap->P_loc)->data; 224 p_oth = (Mat_SeqAIJ*)(ptap->P_oth)->data; 225 pi_loc = p_loc->i; pj_loc = p_loc->j; 226 pi_oth = p_oth->i; pj_oth = p_oth->j; 227 228 /* first, compute symbolic AP = A_loc*P = A_diag*P_loc + A_off*P_oth */ 229 /*-------------------------------------------------------------------*/ 230 ierr = PetscMalloc((am+2)*sizeof(PetscInt),&api);CHKERRQ(ierr); 231 ptap->api = api; 232 api[0] = 0; 233 234 /* create and initialize a linked list */ 235 armax = ad->rmax+ao->rmax; 236 prmax = PetscMax(p_loc->rmax,p_oth->rmax); 237 nlnk_max = armax*prmax; 238 if (!nlnk_max || nlnk_max > pN) nlnk_max = pN; 239 ierr = PetscLLCondensedCreate(nlnk_max,pN,&lnk,&lnkbt);CHKERRQ(ierr); 240 241 /* Initial FreeSpace size is fill*(nnz(A)+nnz(P)) */ 242 ierr = PetscFreeSpaceGet((PetscInt)(fill*(adi[am]+aoi[am]+pi_loc[pm])),&free_space);CHKERRQ(ierr); 243 current_space = free_space; 244 245 ierr = MatPreallocateInitialize(comm,am,pn,dnz,onz);CHKERRQ(ierr); 246 for (i=0; i<am; i++) { 247 /* diagonal portion of A */ 248 nzi = adi[i+1] - adi[i]; 249 for (j=0; j<nzi; j++){ 250 row = *adj++; 251 pnz = pi_loc[row+1] - pi_loc[row]; 252 Jptr = pj_loc + pi_loc[row]; 253 /* add non-zero cols of P into the sorted linked list lnk */ 254 ierr = PetscLLCondensedAddSorted(pnz,Jptr,lnk,lnkbt);CHKERRQ(ierr); 255 } 256 /* off-diagonal portion of A */ 257 nzi = aoi[i+1] - aoi[i]; 258 for (j=0; j<nzi; j++){ 259 row = *aoj++; 260 pnz = pi_oth[row+1] - pi_oth[row]; 261 Jptr = pj_oth + pi_oth[row]; 262 ierr = PetscLLCondensedAddSorted(pnz,Jptr,lnk,lnkbt);CHKERRQ(ierr); 263 } 264 265 apnz = lnk[0]; 266 api[i+1] = api[i] + apnz; 267 268 /* if free space is not available, double the total space in the list */ 269 if (current_space->local_remaining<apnz) { 270 ierr = PetscFreeSpaceGet(apnz+current_space->total_array_size,¤t_space);CHKERRQ(ierr); 271 nspacedouble++; 272 } 273 274 /* Copy data into free space, then initialize lnk */ 275 ierr = PetscLLCondensedClean(pN,apnz,current_space->array,lnk,lnkbt);CHKERRQ(ierr); 276 ierr = MatPreallocateSet(i+rstart,apnz,current_space->array,dnz,onz);CHKERRQ(ierr); 277 current_space->array += apnz; 278 current_space->local_used += apnz; 279 current_space->local_remaining -= apnz; 280 } 281 282 /* Allocate space for apj, initialize apj, and */ 283 /* destroy list of free space and other temporary array(s) */ 284 ierr = PetscMalloc((api[am]+1)*sizeof(PetscInt),&ptap->apj);CHKERRQ(ierr); 285 apj = ptap->apj; 286 ierr = PetscFreeSpaceContiguous(&free_space,ptap->apj);CHKERRQ(ierr); 287 ierr = PetscLLDestroy(lnk,lnkbt);CHKERRQ(ierr); 288 289 /* malloc apa to store dense row A[i,:]*P */ 290 ierr = PetscMalloc(pN*sizeof(PetscScalar),&apa);CHKERRQ(ierr); 291 ierr = PetscMemzero(apa,pN*sizeof(PetscScalar));CHKERRQ(ierr); 292 ptap->apa = apa; 293 294 /* create and assemble symbolic parallel matrix Cmpi */ 295 /*----------------------------------------------------*/ 296 ierr = MatCreate(comm,&Cmpi);CHKERRQ(ierr); 297 ierr = MatSetSizes(Cmpi,am,pn,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 298 ierr = MatSetBlockSizes(Cmpi,A->rmap->bs,P->cmap->bs);CHKERRQ(ierr); 299 300 ierr = MatSetType(Cmpi,MATMPIAIJ);CHKERRQ(ierr); 301 ierr = MatMPIAIJSetPreallocation(Cmpi,0,dnz,0,onz);CHKERRQ(ierr); 302 ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr); 303 for (i=0; i<am; i++){ 304 row = i + rstart; 305 apnz = api[i+1] - api[i]; 306 ierr = MatSetValues(Cmpi,1,&row,apnz,apj,apa,INSERT_VALUES);CHKERRQ(ierr); 307 apj += apnz; 308 } 309 ierr = MatAssemblyBegin(Cmpi,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 310 ierr = MatAssemblyEnd(Cmpi,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 311 312 ptap->destroy = Cmpi->ops->destroy; 313 ptap->duplicate = Cmpi->ops->duplicate; 314 Cmpi->ops->destroy = MatDestroy_MPIAIJ_MatMatMult; 315 Cmpi->ops->duplicate = MatDuplicate_MPIAIJ_MatMatMult; 316 317 /* attach the supporting struct to Cmpi for reuse */ 318 c = (Mat_MPIAIJ*)Cmpi->data; 319 c->ptap = ptap; 320 321 *C = Cmpi; 322 323 /* set MatInfo */ 324 afill = (PetscReal)api[am]/(adi[am]+aoi[am]+pi_loc[pm]+1) + 1.e-5; 325 if (afill < 1.0) afill = 1.0; 326 Cmpi->info.mallocs = nspacedouble; 327 Cmpi->info.fill_ratio_given = fill; 328 Cmpi->info.fill_ratio_needed = afill; 329 330 #if defined(PETSC_USE_INFO) 331 if (api[am]) { 332 ierr = PetscInfo3(Cmpi,"Reallocs %D; Fill ratio: given %G needed %G.\n",nspacedouble,fill,afill);CHKERRQ(ierr); 333 ierr = PetscInfo1(Cmpi,"Use MatMatMult(A,B,MatReuse,%G,&C) for best performance.;\n",afill);CHKERRQ(ierr); 334 } else { 335 ierr = PetscInfo(Cmpi,"Empty matrix product\n");CHKERRQ(ierr); 336 } 337 #endif 338 PetscFunctionReturn(0); 339 } 340 341 #undef __FUNCT__ 342 #define __FUNCT__ "MatMatMult_MPIAIJ_MPIDense" 343 PetscErrorCode MatMatMult_MPIAIJ_MPIDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 344 { 345 PetscErrorCode ierr; 346 347 PetscFunctionBegin; 348 if (scall == MAT_INITIAL_MATRIX){ 349 ierr = MatMatMultSymbolic_MPIAIJ_MPIDense(A,B,fill,C);CHKERRQ(ierr); 350 } 351 ierr = MatMatMultNumeric_MPIAIJ_MPIDense(A,B,*C);CHKERRQ(ierr); 352 PetscFunctionReturn(0); 353 } 354 355 typedef struct { 356 Mat workB; 357 PetscScalar *rvalues,*svalues; 358 MPI_Request *rwaits,*swaits; 359 } MPIAIJ_MPIDense; 360 361 #undef __FUNCT__ 362 #define __FUNCT__ "MatMPIAIJ_MPIDenseDestroy" 363 PetscErrorCode MatMPIAIJ_MPIDenseDestroy(void *ctx) 364 { 365 MPIAIJ_MPIDense *contents = (MPIAIJ_MPIDense*) ctx; 366 PetscErrorCode ierr; 367 368 PetscFunctionBegin; 369 ierr = MatDestroy(&contents->workB);CHKERRQ(ierr); 370 ierr = PetscFree4(contents->rvalues,contents->svalues,contents->rwaits,contents->swaits);CHKERRQ(ierr); 371 ierr = PetscFree(contents);CHKERRQ(ierr); 372 PetscFunctionReturn(0); 373 } 374 375 #undef __FUNCT__ 376 #define __FUNCT__ "MatMatMultSymbolic_MPIAIJ_MPIDense" 377 PetscErrorCode MatMatMultSymbolic_MPIAIJ_MPIDense(Mat A,Mat B,PetscReal fill,Mat *C) 378 { 379 PetscErrorCode ierr; 380 Mat_MPIAIJ *aij = (Mat_MPIAIJ*) A->data; 381 PetscInt nz = aij->B->cmap->n; 382 PetscContainer container; 383 MPIAIJ_MPIDense *contents; 384 VecScatter ctx = aij->Mvctx; 385 VecScatter_MPI_General *from = (VecScatter_MPI_General*) ctx->fromdata; 386 VecScatter_MPI_General *to = ( VecScatter_MPI_General*) ctx->todata; 387 PetscInt m=A->rmap->n,n=B->cmap->n; 388 389 PetscFunctionBegin; 390 ierr = MatCreate(((PetscObject)B)->comm,C);CHKERRQ(ierr); 391 ierr = MatSetSizes(*C,m,n,A->rmap->N,B->cmap->N);CHKERRQ(ierr); 392 ierr = MatSetBlockSizes(*C,A->rmap->bs,B->cmap->bs);CHKERRQ(ierr); 393 ierr = MatSetType(*C,MATMPIDENSE);CHKERRQ(ierr); 394 ierr = MatMPIDenseSetPreallocation(*C,PETSC_NULL);CHKERRQ(ierr); 395 ierr = MatAssemblyBegin(*C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 396 ierr = MatAssemblyEnd(*C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 397 (*C)->ops->matmultnumeric = MatMatMultNumeric_MPIAIJ_MPIDense; 398 399 ierr = PetscNew(MPIAIJ_MPIDense,&contents);CHKERRQ(ierr); 400 /* Create work matrix used to store off processor rows of B needed for local product */ 401 ierr = MatCreateSeqDense(PETSC_COMM_SELF,nz,B->cmap->N,PETSC_NULL,&contents->workB);CHKERRQ(ierr); 402 /* Create work arrays needed */ 403 ierr = PetscMalloc4(B->cmap->N*from->starts[from->n],PetscScalar,&contents->rvalues, 404 B->cmap->N*to->starts[to->n],PetscScalar,&contents->svalues, 405 from->n,MPI_Request,&contents->rwaits, 406 to->n,MPI_Request,&contents->swaits);CHKERRQ(ierr); 407 408 ierr = PetscContainerCreate(((PetscObject)A)->comm,&container);CHKERRQ(ierr); 409 ierr = PetscContainerSetPointer(container,contents);CHKERRQ(ierr); 410 ierr = PetscContainerSetUserDestroy(container,MatMPIAIJ_MPIDenseDestroy);CHKERRQ(ierr); 411 ierr = PetscObjectCompose((PetscObject)(*C),"workB",(PetscObject)container);CHKERRQ(ierr); 412 ierr = PetscContainerDestroy(&container);CHKERRQ(ierr); 413 PetscFunctionReturn(0); 414 } 415 416 #undef __FUNCT__ 417 #define __FUNCT__ "MatMPIDenseScatter" 418 /* 419 Performs an efficient scatter on the rows of B needed by this process; this is 420 a modification of the VecScatterBegin_() routines. 421 */ 422 PetscErrorCode MatMPIDenseScatter(Mat A,Mat B,Mat C,Mat *outworkB) 423 { 424 Mat_MPIAIJ *aij = (Mat_MPIAIJ*)A->data; 425 PetscErrorCode ierr; 426 PetscScalar *b,*w,*svalues,*rvalues; 427 VecScatter ctx = aij->Mvctx; 428 VecScatter_MPI_General *from = (VecScatter_MPI_General*) ctx->fromdata; 429 VecScatter_MPI_General *to = ( VecScatter_MPI_General*) ctx->todata; 430 PetscInt i,j,k; 431 PetscInt *sindices,*sstarts,*rindices,*rstarts; 432 PetscMPIInt *sprocs,*rprocs,nrecvs; 433 MPI_Request *swaits,*rwaits; 434 MPI_Comm comm = ((PetscObject)A)->comm; 435 PetscMPIInt tag = ((PetscObject)ctx)->tag,ncols = B->cmap->N, nrows = aij->B->cmap->n,imdex,nrowsB = B->rmap->n; 436 MPI_Status status; 437 MPIAIJ_MPIDense *contents; 438 PetscContainer container; 439 Mat workB; 440 441 PetscFunctionBegin; 442 ierr = PetscObjectQuery((PetscObject)C,"workB",(PetscObject*)&container);CHKERRQ(ierr); 443 if (!container) SETERRQ(comm,PETSC_ERR_PLIB,"Container does not exist"); 444 ierr = PetscContainerGetPointer(container,(void**)&contents);CHKERRQ(ierr); 445 446 workB = *outworkB = contents->workB; 447 if (nrows != workB->rmap->n) SETERRQ2(comm,PETSC_ERR_PLIB,"Number of rows of workB %D not equal to columns of aij->B %D",nrows,workB->cmap->n); 448 sindices = to->indices; 449 sstarts = to->starts; 450 sprocs = to->procs; 451 swaits = contents->swaits; 452 svalues = contents->svalues; 453 454 rindices = from->indices; 455 rstarts = from->starts; 456 rprocs = from->procs; 457 rwaits = contents->rwaits; 458 rvalues = contents->rvalues; 459 460 ierr = MatDenseGetArray(B,&b);CHKERRQ(ierr); 461 ierr = MatDenseGetArray(workB,&w);CHKERRQ(ierr); 462 463 for (i=0; i<from->n; i++) { 464 ierr = MPI_Irecv(rvalues+ncols*rstarts[i],ncols*(rstarts[i+1]-rstarts[i]),MPIU_SCALAR,rprocs[i],tag,comm,rwaits+i);CHKERRQ(ierr); 465 } 466 467 for (i=0; i<to->n; i++) { 468 /* pack a message at a time */ 469 CHKMEMQ; 470 for (j=0; j<sstarts[i+1]-sstarts[i]; j++){ 471 for (k=0; k<ncols; k++) { 472 svalues[ncols*(sstarts[i] + j) + k] = b[sindices[sstarts[i]+j] + nrowsB*k]; 473 } 474 } 475 CHKMEMQ; 476 ierr = MPI_Isend(svalues+ncols*sstarts[i],ncols*(sstarts[i+1]-sstarts[i]),MPIU_SCALAR,sprocs[i],tag,comm,swaits+i);CHKERRQ(ierr); 477 } 478 479 nrecvs = from->n; 480 while (nrecvs) { 481 ierr = MPI_Waitany(from->n,rwaits,&imdex,&status);CHKERRQ(ierr); 482 nrecvs--; 483 /* unpack a message at a time */ 484 CHKMEMQ; 485 for (j=0; j<rstarts[imdex+1]-rstarts[imdex]; j++){ 486 for (k=0; k<ncols; k++) { 487 w[rindices[rstarts[imdex]+j] + nrows*k] = rvalues[ncols*(rstarts[imdex] + j) + k]; 488 } 489 } 490 CHKMEMQ; 491 } 492 if (to->n) {ierr = MPI_Waitall(to->n,swaits,to->sstatus);CHKERRQ(ierr);} 493 494 ierr = MatDenseRestoreArray(B,&b);CHKERRQ(ierr); 495 ierr = MatDenseRestoreArray(workB,&w);CHKERRQ(ierr); 496 ierr = MatAssemblyBegin(workB,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 497 ierr = MatAssemblyEnd(workB,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 498 PetscFunctionReturn(0); 499 } 500 extern PetscErrorCode MatMatMultNumericAdd_SeqAIJ_SeqDense(Mat,Mat,Mat); 501 502 #undef __FUNCT__ 503 #define __FUNCT__ "MatMatMultNumeric_MPIAIJ_MPIDense" 504 PetscErrorCode MatMatMultNumeric_MPIAIJ_MPIDense(Mat A,Mat B,Mat C) 505 { 506 PetscErrorCode ierr; 507 Mat_MPIAIJ *aij = (Mat_MPIAIJ*)A->data; 508 Mat_MPIDense *bdense = (Mat_MPIDense*)B->data; 509 Mat_MPIDense *cdense = (Mat_MPIDense*)C->data; 510 Mat workB; 511 512 PetscFunctionBegin; 513 514 /* diagonal block of A times all local rows of B*/ 515 ierr = MatMatMultNumeric_SeqAIJ_SeqDense(aij->A,bdense->A,cdense->A);CHKERRQ(ierr); 516 517 /* get off processor parts of B needed to complete the product */ 518 ierr = MatMPIDenseScatter(A,B,C,&workB);CHKERRQ(ierr); 519 520 /* off-diagonal block of A times nonlocal rows of B */ 521 ierr = MatMatMultNumericAdd_SeqAIJ_SeqDense(aij->B,workB,cdense->A);CHKERRQ(ierr); 522 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 523 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 524 PetscFunctionReturn(0); 525 } 526 527 #undef __FUNCT__ 528 #define __FUNCT__ "MatMatMultNumeric_MPIAIJ_MPIAIJ_Scalable" 529 PetscErrorCode MatMatMultNumeric_MPIAIJ_MPIAIJ_Scalable(Mat A,Mat P,Mat C) 530 { 531 PetscErrorCode ierr; 532 Mat_MPIAIJ *a=(Mat_MPIAIJ*)A->data,*c=(Mat_MPIAIJ*)C->data; 533 Mat_SeqAIJ *ad=(Mat_SeqAIJ*)(a->A)->data,*ao=(Mat_SeqAIJ*)(a->B)->data; 534 Mat_SeqAIJ *cd=(Mat_SeqAIJ*)(c->A)->data,*co=(Mat_SeqAIJ*)(c->B)->data; 535 PetscInt *adi=ad->i,*adj,*aoi=ao->i,*aoj; 536 PetscScalar *ada,*aoa,*cda=cd->a,*coa=co->a; 537 Mat_SeqAIJ *p_loc,*p_oth; 538 PetscInt *pi_loc,*pj_loc,*pi_oth,*pj_oth,*pj; 539 PetscScalar *pa_loc,*pa_oth,*pa,valtmp,*ca; 540 PetscInt cm=C->rmap->n,anz,pnz; 541 Mat_PtAPMPI *ptap=c->ptap; 542 PetscScalar *apa_sparse=ptap->apa; 543 PetscInt *api,*apj,*apJ,i,j,k,row; 544 PetscInt cstart=C->cmap->rstart; 545 PetscInt cdnz,conz,k0,k1,nextp; 546 547 PetscFunctionBegin; 548 /* 1) get P_oth = ptap->P_oth and P_loc = ptap->P_loc */ 549 /*-----------------------------------------------------*/ 550 /* update numerical values of P_oth and P_loc */ 551 ierr = MatGetBrowsOfAoCols_MPIAIJ(A,P,MAT_REUSE_MATRIX,&ptap->startsj_s,&ptap->startsj_r,&ptap->bufa,&ptap->P_oth);CHKERRQ(ierr); 552 ierr = MatMPIAIJGetLocalMat(P,MAT_REUSE_MATRIX,&ptap->P_loc);CHKERRQ(ierr); 553 554 /* 2) compute numeric C_loc = A_loc*P = Ad*P_loc + Ao*P_oth */ 555 /*----------------------------------------------------------*/ 556 /* get data from symbolic products */ 557 p_loc = (Mat_SeqAIJ*)(ptap->P_loc)->data; 558 p_oth = (Mat_SeqAIJ*)(ptap->P_oth)->data; 559 pi_loc=p_loc->i; pj_loc=p_loc->j; pa_loc=p_loc->a; 560 pi_oth=p_oth->i; pj_oth=p_oth->j; pa_oth=p_oth->a; 561 562 api = ptap->api; 563 apj = ptap->apj; 564 for (i=0; i<cm; i++) { 565 apJ = apj + api[i]; 566 567 /* diagonal portion of A */ 568 anz = adi[i+1] - adi[i]; 569 adj = ad->j + adi[i]; 570 ada = ad->a + adi[i]; 571 for (j=0; j<anz; j++) { 572 row = adj[j]; 573 pnz = pi_loc[row+1] - pi_loc[row]; 574 pj = pj_loc + pi_loc[row]; 575 pa = pa_loc + pi_loc[row]; 576 /* perform sparse axpy */ 577 valtmp = ada[j]; 578 nextp = 0; 579 for (k=0; nextp<pnz; k++) { 580 if (apJ[k] == pj[nextp]) { /* column of AP == column of P */ 581 apa_sparse[k] += valtmp*pa[nextp++]; 582 } 583 } 584 ierr = PetscLogFlops(2.0*pnz);CHKERRQ(ierr); 585 } 586 587 /* off-diagonal portion of A */ 588 anz = aoi[i+1] - aoi[i]; 589 aoj = ao->j + aoi[i]; 590 aoa = ao->a + aoi[i]; 591 for (j=0; j<anz; j++) { 592 row = aoj[j]; 593 pnz = pi_oth[row+1] - pi_oth[row]; 594 pj = pj_oth + pi_oth[row]; 595 pa = pa_oth + pi_oth[row]; 596 /* perform sparse axpy */ 597 valtmp = aoa[j]; 598 nextp = 0; 599 for (k=0; nextp<pnz; k++) { 600 if (apJ[k] == pj[nextp]) { /* column of AP == column of P */ 601 apa_sparse[k] += valtmp*pa[nextp++]; 602 } 603 } 604 ierr = PetscLogFlops(2.0*pnz);CHKERRQ(ierr); 605 } 606 607 /* set values in C */ 608 cdnz = cd->i[i+1] - cd->i[i]; 609 conz = co->i[i+1] - co->i[i]; 610 611 /* 1st off-diagoanl part of C */ 612 ca = coa + co->i[i]; 613 k = 0; 614 for (k0=0; k0<conz; k0++){ 615 if (apJ[k] >= cstart) break; 616 ca[k0] = apa_sparse[k]; 617 apa_sparse[k] = 0.0; 618 k++; 619 } 620 621 /* diagonal part of C */ 622 ca = cda + cd->i[i]; 623 for (k1=0; k1<cdnz; k1++){ 624 ca[k1] = apa_sparse[k]; 625 apa_sparse[k] = 0.0; 626 k++; 627 } 628 629 /* 2nd off-diagoanl part of C */ 630 ca = coa + co->i[i]; 631 for (; k0<conz; k0++){ 632 ca[k0] = apa_sparse[k]; 633 apa_sparse[k] = 0.0; 634 k++; 635 } 636 } 637 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 638 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 639 PetscFunctionReturn(0); 640 } 641 642 /* same as MatMatMultSymbolic_MPIAIJ_MPIAIJ(), except using LLCondensed to avoid O(BN) memory requirement */ 643 #undef __FUNCT__ 644 #define __FUNCT__ "MatMatMultSymbolic_MPIAIJ_MPIAIJ_Scalable" 645 PetscErrorCode MatMatMultSymbolic_MPIAIJ_MPIAIJ_Scalable(Mat A,Mat P,PetscReal fill,Mat *C) 646 { 647 PetscErrorCode ierr; 648 MPI_Comm comm=((PetscObject)A)->comm; 649 Mat Cmpi; 650 Mat_PtAPMPI *ptap; 651 PetscFreeSpaceList free_space=PETSC_NULL,current_space=PETSC_NULL; 652 Mat_MPIAIJ *a=(Mat_MPIAIJ*)A->data,*c; 653 Mat_SeqAIJ *ad=(Mat_SeqAIJ*)(a->A)->data,*ao=(Mat_SeqAIJ*)(a->B)->data,*p_loc,*p_oth; 654 PetscInt *pi_loc,*pj_loc,*pi_oth,*pj_oth,*dnz,*onz; 655 PetscInt *adi=ad->i,*adj=ad->j,*aoi=ao->i,*aoj=ao->j,rstart=A->rmap->rstart; 656 PetscInt i,pnz,row,*api,*apj,*Jptr,apnz,nspacedouble=0,j,nzi,*lnk,apnz_max=0; 657 PetscInt am=A->rmap->n,pN=P->cmap->N,pn=P->cmap->n,pm=P->rmap->n; 658 PetscInt nlnk_max,armax,prmax; 659 PetscReal afill; 660 PetscScalar *apa; 661 662 PetscFunctionBegin; 663 /* create struct Mat_PtAPMPI and attached it to C later */ 664 ierr = PetscNew(Mat_PtAPMPI,&ptap);CHKERRQ(ierr); 665 666 /* get P_oth by taking rows of P (= non-zero cols of local A) from other processors */ 667 ierr = MatGetBrowsOfAoCols_MPIAIJ(A,P,MAT_INITIAL_MATRIX,&ptap->startsj_s,&ptap->startsj_r,&ptap->bufa,&ptap->P_oth);CHKERRQ(ierr); 668 669 /* get P_loc by taking all local rows of P */ 670 ierr = MatMPIAIJGetLocalMat(P,MAT_INITIAL_MATRIX,&ptap->P_loc);CHKERRQ(ierr); 671 672 p_loc = (Mat_SeqAIJ*)(ptap->P_loc)->data; 673 p_oth = (Mat_SeqAIJ*)(ptap->P_oth)->data; 674 pi_loc = p_loc->i; pj_loc = p_loc->j; 675 pi_oth = p_oth->i; pj_oth = p_oth->j; 676 677 /* first, compute symbolic AP = A_loc*P = A_diag*P_loc + A_off*P_oth */ 678 /*-------------------------------------------------------------------*/ 679 ierr = PetscMalloc((am+2)*sizeof(PetscInt),&api);CHKERRQ(ierr); 680 ptap->api = api; 681 api[0] = 0; 682 683 /* create and initialize a linked list */ 684 armax = ad->rmax+ao->rmax; 685 prmax = PetscMax(p_loc->rmax,p_oth->rmax); 686 nlnk_max = armax*prmax; 687 if (!nlnk_max || nlnk_max > pN) nlnk_max = pN; 688 ierr = PetscLLCondensedCreate_Scalable(nlnk_max,&lnk);CHKERRQ(ierr); 689 690 /* Initial FreeSpace size is fill*(nnz(A)+nnz(P)) */ 691 ierr = PetscFreeSpaceGet((PetscInt)(fill*(adi[am]+aoi[am]+pi_loc[pm])),&free_space);CHKERRQ(ierr); 692 current_space = free_space; 693 694 ierr = MatPreallocateInitialize(comm,am,pn,dnz,onz);CHKERRQ(ierr); 695 for (i=0; i<am; i++) { 696 /* diagonal portion of A */ 697 nzi = adi[i+1] - adi[i]; 698 for (j=0; j<nzi; j++){ 699 row = *adj++; 700 pnz = pi_loc[row+1] - pi_loc[row]; 701 Jptr = pj_loc + pi_loc[row]; 702 /* add non-zero cols of P into the sorted linked list lnk */ 703 ierr = PetscLLCondensedAddSorted_Scalable(pnz,Jptr,lnk);CHKERRQ(ierr); 704 } 705 /* off-diagonal portion of A */ 706 nzi = aoi[i+1] - aoi[i]; 707 for (j=0; j<nzi; j++){ 708 row = *aoj++; 709 pnz = pi_oth[row+1] - pi_oth[row]; 710 Jptr = pj_oth + pi_oth[row]; 711 ierr = PetscLLCondensedAddSorted_Scalable(pnz,Jptr,lnk);CHKERRQ(ierr); 712 } 713 714 apnz = *lnk; 715 api[i+1] = api[i] + apnz; 716 if (apnz > apnz_max) apnz_max = apnz; 717 718 /* if free space is not available, double the total space in the list */ 719 if (current_space->local_remaining<apnz) { 720 ierr = PetscFreeSpaceGet(apnz+current_space->total_array_size,¤t_space);CHKERRQ(ierr); 721 nspacedouble++; 722 } 723 724 /* Copy data into free space, then initialize lnk */ 725 ierr = PetscLLCondensedClean_Scalable(apnz,current_space->array,lnk);CHKERRQ(ierr); 726 ierr = MatPreallocateSet(i+rstart,apnz,current_space->array,dnz,onz);CHKERRQ(ierr); 727 current_space->array += apnz; 728 current_space->local_used += apnz; 729 current_space->local_remaining -= apnz; 730 } 731 732 /* Allocate space for apj, initialize apj, and */ 733 /* destroy list of free space and other temporary array(s) */ 734 ierr = PetscMalloc((api[am]+1)*sizeof(PetscInt),&ptap->apj);CHKERRQ(ierr); 735 apj = ptap->apj; 736 ierr = PetscFreeSpaceContiguous(&free_space,ptap->apj);CHKERRQ(ierr); 737 ierr = PetscLLCondensedDestroy_Scalable(lnk);CHKERRQ(ierr); 738 739 /* create and assemble symbolic parallel matrix Cmpi */ 740 /*----------------------------------------------------*/ 741 ierr = MatCreate(comm,&Cmpi);CHKERRQ(ierr); 742 ierr = MatSetSizes(Cmpi,am,pn,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 743 ierr = MatSetBlockSizes(Cmpi,A->rmap->bs,P->cmap->bs);CHKERRQ(ierr); 744 ierr = MatSetType(Cmpi,MATMPIAIJ);CHKERRQ(ierr); 745 ierr = MatMPIAIJSetPreallocation(Cmpi,0,dnz,0,onz);CHKERRQ(ierr); 746 ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr); 747 748 /* malloc apa for assembly Cmpi */ 749 ierr = PetscMalloc(apnz_max*sizeof(PetscScalar),&apa);CHKERRQ(ierr); 750 ierr = PetscMemzero(apa,apnz_max*sizeof(PetscScalar));CHKERRQ(ierr); 751 ptap->apa = apa; 752 for (i=0; i<am; i++){ 753 row = i + rstart; 754 apnz = api[i+1] - api[i]; 755 ierr = MatSetValues(Cmpi,1,&row,apnz,apj,apa,INSERT_VALUES);CHKERRQ(ierr); 756 apj += apnz; 757 } 758 ierr = MatAssemblyBegin(Cmpi,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 759 ierr = MatAssemblyEnd(Cmpi,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 760 761 ptap->destroy = Cmpi->ops->destroy; 762 ptap->duplicate = Cmpi->ops->duplicate; 763 Cmpi->ops->matmultnumeric = MatMatMultNumeric_MPIAIJ_MPIAIJ_Scalable; 764 Cmpi->ops->destroy = MatDestroy_MPIAIJ_MatMatMult; 765 Cmpi->ops->duplicate = MatDuplicate_MPIAIJ_MatMatMult; 766 767 /* attach the supporting struct to Cmpi for reuse */ 768 c = (Mat_MPIAIJ*)Cmpi->data; 769 c->ptap = ptap; 770 771 *C = Cmpi; 772 773 /* set MatInfo */ 774 afill = (PetscReal)api[am]/(adi[am]+aoi[am]+pi_loc[pm]+1) + 1.e-5; 775 if (afill < 1.0) afill = 1.0; 776 Cmpi->info.mallocs = nspacedouble; 777 Cmpi->info.fill_ratio_given = fill; 778 Cmpi->info.fill_ratio_needed = afill; 779 780 #if defined(PETSC_USE_INFO) 781 if (api[am]) { 782 ierr = PetscInfo3(Cmpi,"Reallocs %D; Fill ratio: given %G needed %G.\n",nspacedouble,fill,afill);CHKERRQ(ierr); 783 ierr = PetscInfo1(Cmpi,"Use MatMatMult(A,B,MatReuse,%G,&C) for best performance.;\n",afill);CHKERRQ(ierr); 784 } else { 785 ierr = PetscInfo(Cmpi,"Empty matrix product\n");CHKERRQ(ierr); 786 } 787 #endif 788 PetscFunctionReturn(0); 789 } 790 791 /*-------------------------------------------------------------------------*/ 792 #undef __FUNCT__ 793 #define __FUNCT__ "MatTransposeMatMult_MPIAIJ_MPIAIJ" 794 PetscErrorCode MatTransposeMatMult_MPIAIJ_MPIAIJ(Mat P,Mat A,MatReuse scall,PetscReal fill,Mat *C) 795 { 796 PetscErrorCode ierr; 797 PetscBool scalable=PETSC_TRUE,viamatmatmult=PETSC_FALSE; 798 799 PetscFunctionBegin; 800 ierr = PetscOptionsBool("-mattransposematmult_viamatmatmult","Use R=Pt and C=R*A","",viamatmatmult,&viamatmatmult,PETSC_NULL);CHKERRQ(ierr); 801 if (viamatmatmult){ 802 Mat Pt; 803 Mat_PtAPMPI *ptap; 804 Mat_MPIAIJ *c; 805 if (scall == MAT_INITIAL_MATRIX){ 806 ierr = MatTranspose(P,MAT_INITIAL_MATRIX,&Pt);CHKERRQ(ierr); 807 ierr = MatMatMult(Pt,A,MAT_INITIAL_MATRIX,fill,C);CHKERRQ(ierr); 808 809 c = (Mat_MPIAIJ*)(*C)->data; 810 ptap = c->ptap; 811 ptap->Pt = Pt; 812 } else if (scall == MAT_REUSE_MATRIX){ 813 c = (Mat_MPIAIJ*)(*C)->data; 814 ptap = c->ptap; 815 Pt = ptap->Pt; 816 ierr = MatTranspose(P,scall,&Pt);CHKERRQ(ierr); 817 ierr = MatMatMult(Pt,A,scall,fill,C);CHKERRQ(ierr); 818 } else { 819 SETERRQ(((PetscObject)A)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not supported"); 820 } 821 PetscFunctionReturn(0); 822 } 823 824 if (scall == MAT_INITIAL_MATRIX){ 825 ierr = PetscObjectOptionsBegin((PetscObject)A);CHKERRQ(ierr); 826 ierr = PetscOptionsBool("-mattransposematmult_scalable","Use a scalable but slower C=Pt*A","",scalable,&scalable,PETSC_NULL);CHKERRQ(ierr); 827 if (scalable){ 828 ierr = MatTransposeMatMultSymbolic_MPIAIJ_MPIAIJ_Scalable(P,A,fill,C);CHKERRQ(ierr); 829 } else { 830 ierr = MatTransposeMatMultSymbolic_MPIAIJ_MPIAIJ(P,A,fill,C);CHKERRQ(ierr); 831 } 832 ierr = PetscOptionsEnd();CHKERRQ(ierr); 833 } 834 ierr = (*(*C)->ops->mattransposemultnumeric)(P,A,*C);CHKERRQ(ierr); 835 PetscFunctionReturn(0); 836 } 837 838 #undef __FUNCT__ 839 #define __FUNCT__ "MatTransposeMatMultNumeric_MPIAIJ_MPIAIJ" 840 PetscErrorCode MatTransposeMatMultNumeric_MPIAIJ_MPIAIJ(Mat P,Mat A,Mat C) 841 { 842 PetscErrorCode ierr; 843 Mat_Merge_SeqsToMPI *merge; 844 Mat_MPIAIJ *p=(Mat_MPIAIJ*)P->data,*c=(Mat_MPIAIJ*)C->data; 845 Mat_SeqAIJ *pd=(Mat_SeqAIJ*)(p->A)->data,*po=(Mat_SeqAIJ*)(p->B)->data; 846 Mat_PtAPMPI *ptap; 847 PetscInt *adj,*aJ; 848 PetscInt i,j,k,anz,pnz,row,*cj; 849 MatScalar *ada,*aval,*ca,valtmp; 850 PetscInt am=A->rmap->n,cm=C->rmap->n,pon=(p->B)->cmap->n; 851 MPI_Comm comm=((PetscObject)C)->comm; 852 PetscMPIInt size,rank,taga,*len_s; 853 PetscInt *owners,proc,nrows,**buf_ri_k,**nextrow,**nextci; 854 PetscInt **buf_ri,**buf_rj; 855 PetscInt cnz=0,*bj_i,*bi,*bj,bnz,nextcj; /* bi,bj,ba: local array of C(mpi mat) */ 856 MPI_Request *s_waits,*r_waits; 857 MPI_Status *status; 858 MatScalar **abuf_r,*ba_i,*pA,*coa,*ba; 859 PetscInt *ai,*aj,*coi,*coj; 860 PetscInt *poJ,*pdJ; 861 Mat A_loc; 862 Mat_SeqAIJ *a_loc; 863 864 PetscFunctionBegin; 865 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 866 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 867 868 ptap = c->ptap; 869 merge = ptap->merge; 870 871 /* 2) compute numeric C_seq = P_loc^T*A_loc*P - dominating part */ 872 /*--------------------------------------------------------------*/ 873 /* get data from symbolic products */ 874 coi = merge->coi; coj = merge->coj; 875 ierr = PetscMalloc((coi[pon]+1)*sizeof(MatScalar),&coa);CHKERRQ(ierr); 876 ierr = PetscMemzero(coa,coi[pon]*sizeof(MatScalar));CHKERRQ(ierr); 877 878 bi = merge->bi; bj = merge->bj; 879 owners = merge->rowmap->range; 880 ierr = PetscMalloc((bi[cm]+1)*sizeof(MatScalar),&ba);CHKERRQ(ierr); 881 ierr = PetscMemzero(ba,bi[cm]*sizeof(MatScalar));CHKERRQ(ierr); 882 883 /* get A_loc by taking all local rows of A */ 884 A_loc = ptap->A_loc; 885 ierr = MatMPIAIJGetLocalMat(A,MAT_REUSE_MATRIX,&A_loc);CHKERRQ(ierr); 886 a_loc = (Mat_SeqAIJ*)(A_loc)->data; 887 ai = a_loc->i; 888 aj = a_loc->j; 889 890 ierr = PetscMalloc((A->cmap->N)*sizeof(PetscScalar),&aval);CHKERRQ(ierr); /* non-scalable!!! */ 891 ierr = PetscMemzero(aval,A->cmap->N*sizeof(PetscScalar));CHKERRQ(ierr); 892 893 for (i=0; i<am; i++) { 894 /* 2-a) put A[i,:] to dense array aval */ 895 anz = ai[i+1] - ai[i]; 896 adj = aj + ai[i]; 897 ada = a_loc->a + ai[i]; 898 for (j=0; j<anz; j++){ 899 aval[adj[j]] = ada[j]; 900 } 901 902 /* 2-b) Compute Cseq = P_loc[i,:]^T*A[i,:] using outer product */ 903 /*--------------------------------------------------------------*/ 904 /* put the value into Co=(p->B)^T*A (off-diagonal part, send to others) */ 905 pnz = po->i[i+1] - po->i[i]; 906 poJ = po->j + po->i[i]; 907 pA = po->a + po->i[i]; 908 for (j=0; j<pnz; j++){ 909 row = poJ[j]; 910 cnz = coi[row+1] - coi[row]; 911 cj = coj + coi[row]; 912 ca = coa + coi[row]; 913 /* perform dense axpy */ 914 valtmp = pA[j]; 915 for (k=0; k<cnz; k++) { 916 ca[k] += valtmp*aval[cj[k]]; 917 } 918 ierr = PetscLogFlops(2.0*cnz);CHKERRQ(ierr); 919 } 920 921 /* put the value into Cd (diagonal part) */ 922 pnz = pd->i[i+1] - pd->i[i]; 923 pdJ = pd->j + pd->i[i]; 924 pA = pd->a + pd->i[i]; 925 for (j=0; j<pnz; j++){ 926 row = pdJ[j]; 927 cnz = bi[row+1] - bi[row]; 928 cj = bj + bi[row]; 929 ca = ba + bi[row]; 930 /* perform dense axpy */ 931 valtmp = pA[j]; 932 for (k=0; k<cnz; k++) { 933 ca[k] += valtmp*aval[cj[k]]; 934 } 935 ierr = PetscLogFlops(2.0*cnz);CHKERRQ(ierr); 936 } 937 938 /* zero the current row of Pt*A */ 939 aJ = aj + ai[i]; 940 for (k=0; k<anz; k++) aval[aJ[k]] = 0.0; 941 } 942 943 /* 3) send and recv matrix values coa */ 944 /*------------------------------------*/ 945 buf_ri = merge->buf_ri; 946 buf_rj = merge->buf_rj; 947 len_s = merge->len_s; 948 ierr = PetscCommGetNewTag(comm,&taga);CHKERRQ(ierr); 949 ierr = PetscPostIrecvScalar(comm,taga,merge->nrecv,merge->id_r,merge->len_r,&abuf_r,&r_waits);CHKERRQ(ierr); 950 951 ierr = PetscMalloc2(merge->nsend+1,MPI_Request,&s_waits,size,MPI_Status,&status);CHKERRQ(ierr); 952 for (proc=0,k=0; proc<size; proc++){ 953 if (!len_s[proc]) continue; 954 i = merge->owners_co[proc]; 955 ierr = MPI_Isend(coa+coi[i],len_s[proc],MPIU_MATSCALAR,proc,taga,comm,s_waits+k);CHKERRQ(ierr); 956 k++; 957 } 958 if (merge->nrecv) {ierr = MPI_Waitall(merge->nrecv,r_waits,status);CHKERRQ(ierr);} 959 if (merge->nsend) {ierr = MPI_Waitall(merge->nsend,s_waits,status);CHKERRQ(ierr);} 960 961 ierr = PetscFree2(s_waits,status);CHKERRQ(ierr); 962 ierr = PetscFree(r_waits);CHKERRQ(ierr); 963 ierr = PetscFree(coa);CHKERRQ(ierr); 964 965 /* 4) insert local Cseq and received values into Cmpi */ 966 /*----------------------------------------------------*/ 967 ierr = PetscMalloc3(merge->nrecv,PetscInt**,&buf_ri_k,merge->nrecv,PetscInt*,&nextrow,merge->nrecv,PetscInt*,&nextci);CHKERRQ(ierr); 968 for (k=0; k<merge->nrecv; k++){ 969 buf_ri_k[k] = buf_ri[k]; /* beginning of k-th recved i-structure */ 970 nrows = *(buf_ri_k[k]); 971 nextrow[k] = buf_ri_k[k]+1; /* next row number of k-th recved i-structure */ 972 nextci[k] = buf_ri_k[k] + (nrows + 1);/* poins to the next i-structure of k-th recved i-structure */ 973 } 974 975 for (i=0; i<cm; i++) { 976 row = owners[rank] + i; /* global row index of C_seq */ 977 bj_i = bj + bi[i]; /* col indices of the i-th row of C */ 978 ba_i = ba + bi[i]; 979 bnz = bi[i+1] - bi[i]; 980 /* add received vals into ba */ 981 for (k=0; k<merge->nrecv; k++){ /* k-th received message */ 982 /* i-th row */ 983 if (i == *nextrow[k]) { 984 cnz = *(nextci[k]+1) - *nextci[k]; 985 cj = buf_rj[k] + *(nextci[k]); 986 ca = abuf_r[k] + *(nextci[k]); 987 nextcj = 0; 988 for (j=0; nextcj<cnz; j++){ 989 if (bj_i[j] == cj[nextcj]){ /* bcol == ccol */ 990 ba_i[j] += ca[nextcj++]; 991 } 992 } 993 nextrow[k]++; nextci[k]++; 994 ierr = PetscLogFlops(2.0*cnz);CHKERRQ(ierr); 995 } 996 } 997 ierr = MatSetValues(C,1,&row,bnz,bj_i,ba_i,INSERT_VALUES);CHKERRQ(ierr); 998 } 999 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1000 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1001 1002 ierr = PetscFree(ba);CHKERRQ(ierr); 1003 ierr = PetscFree(abuf_r[0]);CHKERRQ(ierr); 1004 ierr = PetscFree(abuf_r);CHKERRQ(ierr); 1005 ierr = PetscFree3(buf_ri_k,nextrow,nextci);CHKERRQ(ierr); 1006 ierr = PetscFree(aval);CHKERRQ(ierr); 1007 PetscFunctionReturn(0); 1008 } 1009 1010 /* This routine is modified from MatPtAPSymbolic_MPIAIJ_MPIAIJ() */ 1011 #undef __FUNCT__ 1012 #define __FUNCT__ "MatTransposeMatMultSymbolic_MPIAIJ_MPIAIJ" 1013 PetscErrorCode MatTransposeMatMultSymbolic_MPIAIJ_MPIAIJ(Mat P,Mat A,PetscReal fill,Mat *C) 1014 { 1015 PetscErrorCode ierr; 1016 Mat Cmpi,A_loc,POt,PDt; 1017 Mat_PtAPMPI *ptap; 1018 PetscFreeSpaceList free_space=PETSC_NULL,current_space=PETSC_NULL; 1019 Mat_MPIAIJ *p=(Mat_MPIAIJ*)P->data,*c; 1020 PetscInt *pdti,*pdtj,*poti,*potj,*ptJ; 1021 PetscInt nnz; 1022 PetscInt *lnk,*owners_co,*coi,*coj,i,k,pnz,row; 1023 PetscInt am=A->rmap->n,pn=P->cmap->n; 1024 PetscBT lnkbt; 1025 MPI_Comm comm=((PetscObject)A)->comm; 1026 PetscMPIInt size,rank,tagi,tagj,*len_si,*len_s,*len_ri; 1027 PetscInt **buf_rj,**buf_ri,**buf_ri_k; 1028 PetscInt len,proc,*dnz,*onz,*owners; 1029 PetscInt nzi,*bi,*bj; 1030 PetscInt nrows,*buf_s,*buf_si,*buf_si_i,**nextrow,**nextci; 1031 MPI_Request *swaits,*rwaits; 1032 MPI_Status *sstatus,rstatus; 1033 Mat_Merge_SeqsToMPI *merge; 1034 PetscInt *ai,*aj,*Jptr,anz,*prmap=p->garray,pon,nspacedouble=0,j; 1035 PetscReal afill=1.0,afill_tmp; 1036 PetscInt rstart = P->cmap->rstart,rmax,aN=A->cmap->N,Crmax; 1037 PetscScalar *vals; 1038 Mat_SeqAIJ *a_loc, *pdt,*pot; 1039 1040 PetscFunctionBegin; 1041 /* check if matrix local sizes are compatible */ 1042 if (A->rmap->rstart != P->rmap->rstart || A->rmap->rend != P->rmap->rend){ 1043 SETERRQ4(comm,PETSC_ERR_ARG_SIZ,"Matrix local dimensions are incompatible, A (%D, %D) != P (%D,%D)",A->rmap->rstart,A->rmap->rend,P->rmap->rstart,P->rmap->rend); 1044 } 1045 1046 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 1047 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 1048 1049 /* create struct Mat_PtAPMPI and attached it to C later */ 1050 ierr = PetscNew(Mat_PtAPMPI,&ptap);CHKERRQ(ierr); 1051 1052 /* get A_loc by taking all local rows of A */ 1053 ierr = MatMPIAIJGetLocalMat(A,MAT_INITIAL_MATRIX,&A_loc);CHKERRQ(ierr); 1054 ptap->A_loc = A_loc; 1055 a_loc = (Mat_SeqAIJ*)(A_loc)->data; 1056 ai = a_loc->i; 1057 aj = a_loc->j; 1058 1059 /* determine symbolic Co=(p->B)^T*A - send to others */ 1060 /*----------------------------------------------------*/ 1061 ierr = MatTransposeSymbolic_SeqAIJ(p->A,&PDt);CHKERRQ(ierr); 1062 pdt = (Mat_SeqAIJ*)PDt->data; 1063 pdti = pdt->i; pdtj = pdt->j; 1064 1065 ierr = MatTransposeSymbolic_SeqAIJ(p->B,&POt);CHKERRQ(ierr); 1066 pot = (Mat_SeqAIJ*)POt->data; 1067 poti = pot->i; potj = pot->j; 1068 1069 /* then, compute symbolic Co = (p->B)^T*A */ 1070 pon = (p->B)->cmap->n; /* total num of rows to be sent to other processors 1071 >= (num of nonzero rows of C_seq) - pn */ 1072 ierr = PetscMalloc((pon+1)*sizeof(PetscInt),&coi);CHKERRQ(ierr); 1073 coi[0] = 0; 1074 1075 /* set initial free space to be fill*(nnz(p->B) + nnz(A)) */ 1076 nnz = fill*(poti[pon] + ai[am]); 1077 ierr = PetscFreeSpaceGet(nnz,&free_space);CHKERRQ(ierr); 1078 current_space = free_space; 1079 1080 /* create and initialize a linked list */ 1081 i = PetscMax(pdt->rmax,pot->rmax); 1082 Crmax = i*a_loc->rmax*size; 1083 if (!Crmax || Crmax > aN) Crmax = aN; 1084 ierr = PetscLLCondensedCreate(Crmax,aN,&lnk,&lnkbt);CHKERRQ(ierr); 1085 1086 for (i=0; i<pon; i++) { 1087 pnz = poti[i+1] - poti[i]; 1088 ptJ = potj + poti[i]; 1089 for (j=0; j<pnz; j++){ 1090 row = ptJ[j]; /* row of A_loc == col of Pot */ 1091 anz = ai[row+1] - ai[row]; 1092 Jptr = aj + ai[row]; 1093 /* add non-zero cols of AP into the sorted linked list lnk */ 1094 ierr = PetscLLCondensedAddSorted(anz,Jptr,lnk,lnkbt);CHKERRQ(ierr); 1095 } 1096 nnz = lnk[0]; 1097 1098 /* If free space is not available, double the total space in the list */ 1099 if (current_space->local_remaining<nnz) { 1100 ierr = PetscFreeSpaceGet(nnz+current_space->total_array_size,¤t_space);CHKERRQ(ierr); 1101 nspacedouble++; 1102 } 1103 1104 /* Copy data into free space, and zero out denserows */ 1105 ierr = PetscLLCondensedClean(aN,nnz,current_space->array,lnk,lnkbt);CHKERRQ(ierr); 1106 current_space->array += nnz; 1107 current_space->local_used += nnz; 1108 current_space->local_remaining -= nnz; 1109 coi[i+1] = coi[i] + nnz; 1110 } 1111 1112 ierr = PetscMalloc((coi[pon]+1)*sizeof(PetscInt),&coj);CHKERRQ(ierr); 1113 ierr = PetscFreeSpaceContiguous(&free_space,coj);CHKERRQ(ierr); 1114 afill_tmp = (PetscReal)coi[pon]/(poti[pon] + ai[am]+1); 1115 if (afill_tmp > afill) afill = afill_tmp; 1116 1117 /* send j-array (coj) of Co to other processors */ 1118 /*----------------------------------------------*/ 1119 /* determine row ownership */ 1120 ierr = PetscNew(Mat_Merge_SeqsToMPI,&merge);CHKERRQ(ierr); 1121 ierr = PetscLayoutCreate(comm,&merge->rowmap);CHKERRQ(ierr); 1122 merge->rowmap->n = pn; 1123 merge->rowmap->bs = 1; 1124 ierr = PetscLayoutSetUp(merge->rowmap);CHKERRQ(ierr); 1125 owners = merge->rowmap->range; 1126 1127 /* determine the number of messages to send, their lengths */ 1128 ierr = PetscMalloc(size*sizeof(PetscMPIInt),&len_si);CHKERRQ(ierr); 1129 ierr = PetscMemzero(len_si,size*sizeof(PetscMPIInt));CHKERRQ(ierr); 1130 ierr = PetscMalloc(size*sizeof(PetscMPIInt),&merge->len_s);CHKERRQ(ierr); 1131 len_s = merge->len_s; 1132 merge->nsend = 0; 1133 1134 ierr = PetscMalloc((size+2)*sizeof(PetscInt),&owners_co);CHKERRQ(ierr); 1135 ierr = PetscMemzero(len_s,size*sizeof(PetscMPIInt));CHKERRQ(ierr); 1136 1137 proc = 0; 1138 for (i=0; i<pon; i++){ 1139 while (prmap[i] >= owners[proc+1]) proc++; 1140 len_si[proc]++; /* num of rows in Co to be sent to [proc] */ 1141 len_s[proc] += coi[i+1] - coi[i]; 1142 } 1143 1144 len = 0; /* max length of buf_si[] */ 1145 owners_co[0] = 0; 1146 for (proc=0; proc<size; proc++){ 1147 owners_co[proc+1] = owners_co[proc] + len_si[proc]; 1148 if (len_si[proc]){ 1149 merge->nsend++; 1150 len_si[proc] = 2*(len_si[proc] + 1); 1151 len += len_si[proc]; 1152 } 1153 } 1154 1155 /* determine the number and length of messages to receive for coi and coj */ 1156 ierr = PetscGatherNumberOfMessages(comm,PETSC_NULL,len_s,&merge->nrecv);CHKERRQ(ierr); 1157 ierr = PetscGatherMessageLengths2(comm,merge->nsend,merge->nrecv,len_s,len_si,&merge->id_r,&merge->len_r,&len_ri);CHKERRQ(ierr); 1158 1159 /* post the Irecv and Isend of coj */ 1160 ierr = PetscCommGetNewTag(comm,&tagj);CHKERRQ(ierr); 1161 ierr = PetscPostIrecvInt(comm,tagj,merge->nrecv,merge->id_r,merge->len_r,&buf_rj,&rwaits);CHKERRQ(ierr); 1162 ierr = PetscMalloc((merge->nsend+1)*sizeof(MPI_Request),&swaits);CHKERRQ(ierr); 1163 for (proc=0, k=0; proc<size; proc++){ 1164 if (!len_s[proc]) continue; 1165 i = owners_co[proc]; 1166 ierr = MPI_Isend(coj+coi[i],len_s[proc],MPIU_INT,proc,tagj,comm,swaits+k);CHKERRQ(ierr); 1167 k++; 1168 } 1169 1170 /* receives and sends of coj are complete */ 1171 ierr = PetscMalloc(size*sizeof(MPI_Status),&sstatus);CHKERRQ(ierr); 1172 for (i=0; i<merge->nrecv; i++){ 1173 PetscMPIInt icompleted; 1174 ierr = MPI_Waitany(merge->nrecv,rwaits,&icompleted,&rstatus);CHKERRQ(ierr); 1175 } 1176 ierr = PetscFree(rwaits);CHKERRQ(ierr); 1177 if (merge->nsend) {ierr = MPI_Waitall(merge->nsend,swaits,sstatus);CHKERRQ(ierr);} 1178 1179 /* send and recv coi */ 1180 /*-------------------*/ 1181 ierr = PetscCommGetNewTag(comm,&tagi);CHKERRQ(ierr); 1182 ierr = PetscPostIrecvInt(comm,tagi,merge->nrecv,merge->id_r,len_ri,&buf_ri,&rwaits);CHKERRQ(ierr); 1183 ierr = PetscMalloc((len+1)*sizeof(PetscInt),&buf_s);CHKERRQ(ierr); 1184 buf_si = buf_s; /* points to the beginning of k-th msg to be sent */ 1185 for (proc=0,k=0; proc<size; proc++){ 1186 if (!len_s[proc]) continue; 1187 /* form outgoing message for i-structure: 1188 buf_si[0]: nrows to be sent 1189 [1:nrows]: row index (global) 1190 [nrows+1:2*nrows+1]: i-structure index 1191 */ 1192 /*-------------------------------------------*/ 1193 nrows = len_si[proc]/2 - 1; 1194 buf_si_i = buf_si + nrows+1; 1195 buf_si[0] = nrows; 1196 buf_si_i[0] = 0; 1197 nrows = 0; 1198 for (i=owners_co[proc]; i<owners_co[proc+1]; i++){ 1199 nzi = coi[i+1] - coi[i]; 1200 buf_si_i[nrows+1] = buf_si_i[nrows] + nzi; /* i-structure */ 1201 buf_si[nrows+1] =prmap[i] -owners[proc]; /* local row index */ 1202 nrows++; 1203 } 1204 ierr = MPI_Isend(buf_si,len_si[proc],MPIU_INT,proc,tagi,comm,swaits+k);CHKERRQ(ierr); 1205 k++; 1206 buf_si += len_si[proc]; 1207 } 1208 i = merge->nrecv; 1209 while (i--) { 1210 PetscMPIInt icompleted; 1211 ierr = MPI_Waitany(merge->nrecv,rwaits,&icompleted,&rstatus);CHKERRQ(ierr); 1212 } 1213 ierr = PetscFree(rwaits);CHKERRQ(ierr); 1214 if (merge->nsend) {ierr = MPI_Waitall(merge->nsend,swaits,sstatus);CHKERRQ(ierr);} 1215 ierr = PetscFree(len_si);CHKERRQ(ierr); 1216 ierr = PetscFree(len_ri);CHKERRQ(ierr); 1217 ierr = PetscFree(swaits);CHKERRQ(ierr); 1218 ierr = PetscFree(sstatus);CHKERRQ(ierr); 1219 ierr = PetscFree(buf_s);CHKERRQ(ierr); 1220 1221 /* compute the local portion of C (mpi mat) */ 1222 /*------------------------------------------*/ 1223 /* allocate bi array and free space for accumulating nonzero column info */ 1224 ierr = PetscMalloc((pn+1)*sizeof(PetscInt),&bi);CHKERRQ(ierr); 1225 bi[0] = 0; 1226 1227 /* set initial free space to be fill*(nnz(P) + nnz(A)) */ 1228 nnz = fill*(pdti[pn] + poti[pon] + ai[am]); 1229 ierr = PetscFreeSpaceGet(nnz,&free_space);CHKERRQ(ierr); 1230 current_space = free_space; 1231 1232 ierr = PetscMalloc3(merge->nrecv,PetscInt**,&buf_ri_k,merge->nrecv,PetscInt*,&nextrow,merge->nrecv,PetscInt*,&nextci);CHKERRQ(ierr); 1233 for (k=0; k<merge->nrecv; k++){ 1234 buf_ri_k[k] = buf_ri[k]; /* beginning of k-th recved i-structure */ 1235 nrows = *buf_ri_k[k]; 1236 nextrow[k] = buf_ri_k[k] + 1; /* next row number of k-th recved i-structure */ 1237 nextci[k] = buf_ri_k[k] + (nrows + 1);/* poins to the next i-structure of k-th recved i-structure */ 1238 } 1239 1240 ierr = MatPreallocateInitialize(comm,pn,A->cmap->n,dnz,onz);CHKERRQ(ierr); 1241 rmax = 0; 1242 for (i=0; i<pn; i++) { 1243 /* add pdt[i,:]*AP into lnk */ 1244 pnz = pdti[i+1] - pdti[i]; 1245 ptJ = pdtj + pdti[i]; 1246 for (j=0; j<pnz; j++){ 1247 row = ptJ[j]; /* row of AP == col of Pt */ 1248 anz = ai[row+1] - ai[row]; 1249 Jptr = aj + ai[row]; 1250 /* add non-zero cols of AP into the sorted linked list lnk */ 1251 ierr = PetscLLCondensedAddSorted(anz,Jptr,lnk,lnkbt);CHKERRQ(ierr); 1252 } 1253 1254 /* add received col data into lnk */ 1255 for (k=0; k<merge->nrecv; k++){ /* k-th received message */ 1256 if (i == *nextrow[k]) { /* i-th row */ 1257 nzi = *(nextci[k]+1) - *nextci[k]; 1258 Jptr = buf_rj[k] + *nextci[k]; 1259 ierr = PetscLLCondensedAddSorted(nzi,Jptr,lnk,lnkbt);CHKERRQ(ierr); 1260 nextrow[k]++; nextci[k]++; 1261 } 1262 } 1263 nnz = lnk[0]; 1264 1265 /* if free space is not available, make more free space */ 1266 if (current_space->local_remaining<nnz) { 1267 ierr = PetscFreeSpaceGet(nnz+current_space->total_array_size,¤t_space);CHKERRQ(ierr); 1268 nspacedouble++; 1269 } 1270 /* copy data into free space, then initialize lnk */ 1271 ierr = PetscLLCondensedClean(aN,nnz,current_space->array,lnk,lnkbt);CHKERRQ(ierr); 1272 ierr = MatPreallocateSet(i+owners[rank],nnz,current_space->array,dnz,onz);CHKERRQ(ierr); 1273 current_space->array += nnz; 1274 current_space->local_used += nnz; 1275 current_space->local_remaining -= nnz; 1276 bi[i+1] = bi[i] + nnz; 1277 if (nnz > rmax) rmax = nnz; 1278 } 1279 ierr = PetscFree3(buf_ri_k,nextrow,nextci);CHKERRQ(ierr); 1280 1281 ierr = PetscMalloc((bi[pn]+1)*sizeof(PetscInt),&bj);CHKERRQ(ierr); 1282 ierr = PetscFreeSpaceContiguous(&free_space,bj);CHKERRQ(ierr); 1283 afill_tmp = (PetscReal)bi[pn]/(pdti[pn] + poti[pon] + ai[am]+1); 1284 if (afill_tmp > afill) afill = afill_tmp; 1285 ierr = PetscLLCondensedDestroy(lnk,lnkbt);CHKERRQ(ierr); 1286 ierr = MatDestroy(&POt);CHKERRQ(ierr); 1287 ierr = MatDestroy(&PDt);CHKERRQ(ierr); 1288 1289 /* create symbolic parallel matrix Cmpi - why cannot be assembled in Numeric part */ 1290 /*----------------------------------------------------------------------------------*/ 1291 ierr = PetscMalloc((rmax+1)*sizeof(PetscScalar),&vals);CHKERRQ(ierr); 1292 ierr = PetscMemzero(vals,rmax*sizeof(PetscScalar));CHKERRQ(ierr); 1293 1294 ierr = MatCreate(comm,&Cmpi);CHKERRQ(ierr); 1295 ierr = MatSetSizes(Cmpi,pn,A->cmap->n,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 1296 ierr = MatSetBlockSizes(Cmpi,P->cmap->bs,A->cmap->bs);CHKERRQ(ierr); 1297 ierr = MatSetType(Cmpi,MATMPIAIJ);CHKERRQ(ierr); 1298 ierr = MatMPIAIJSetPreallocation(Cmpi,0,dnz,0,onz);CHKERRQ(ierr); 1299 ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr); 1300 ierr = MatSetBlockSize(Cmpi,1);CHKERRQ(ierr); 1301 for (i=0; i<pn; i++){ 1302 row = i + rstart; 1303 nnz = bi[i+1] - bi[i]; 1304 Jptr = bj + bi[i]; 1305 ierr = MatSetValues(Cmpi,1,&row,nnz,Jptr,vals,INSERT_VALUES);CHKERRQ(ierr); 1306 } 1307 ierr = MatAssemblyBegin(Cmpi,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1308 ierr = MatAssemblyEnd(Cmpi,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1309 ierr = PetscFree(vals);CHKERRQ(ierr); 1310 1311 merge->bi = bi; 1312 merge->bj = bj; 1313 merge->coi = coi; 1314 merge->coj = coj; 1315 merge->buf_ri = buf_ri; 1316 merge->buf_rj = buf_rj; 1317 merge->owners_co = owners_co; 1318 merge->destroy = Cmpi->ops->destroy; 1319 merge->duplicate = Cmpi->ops->duplicate; 1320 1321 Cmpi->ops->mattransposemultnumeric = MatTransposeMatMultNumeric_MPIAIJ_MPIAIJ; 1322 Cmpi->ops->destroy = MatDestroy_MPIAIJ_PtAP; 1323 1324 /* attach the supporting struct to Cmpi for reuse */ 1325 c = (Mat_MPIAIJ*)Cmpi->data; 1326 c->ptap = ptap; 1327 ptap->api = PETSC_NULL; 1328 ptap->apj = PETSC_NULL; 1329 ptap->merge = merge; 1330 ptap->rmax = rmax; 1331 1332 *C = Cmpi; 1333 #if defined(PETSC_USE_INFO) 1334 if (bi[pn] != 0) { 1335 ierr = PetscInfo3(Cmpi,"Reallocs %D; Fill ratio: given %G needed %G.\n",nspacedouble,fill,afill);CHKERRQ(ierr); 1336 ierr = PetscInfo1(Cmpi,"Use MatTransposeMatMult(A,B,MatReuse,%G,&C) for best performance.\n",afill);CHKERRQ(ierr); 1337 } else { 1338 ierr = PetscInfo(Cmpi,"Empty matrix product\n");CHKERRQ(ierr); 1339 } 1340 #endif 1341 PetscFunctionReturn(0); 1342 } 1343 1344 #undef __FUNCT__ 1345 #define __FUNCT__ "MatTransposeMatMultNumeric_MPIAIJ_MPIAIJ_Scalable" 1346 PetscErrorCode MatTransposeMatMultNumeric_MPIAIJ_MPIAIJ_Scalable(Mat P,Mat A,Mat C) 1347 { 1348 PetscErrorCode ierr; 1349 Mat_Merge_SeqsToMPI *merge; 1350 Mat_MPIAIJ *p=(Mat_MPIAIJ*)P->data,*c=(Mat_MPIAIJ*)C->data; 1351 Mat_SeqAIJ *pd=(Mat_SeqAIJ*)(p->A)->data,*po=(Mat_SeqAIJ*)(p->B)->data; 1352 Mat_PtAPMPI *ptap; 1353 PetscInt *adj; 1354 PetscInt i,j,k,anz,pnz,row,*cj,nexta; 1355 MatScalar *ada,*ca,valtmp; 1356 PetscInt am=A->rmap->n,cm=C->rmap->n,pon=(p->B)->cmap->n; 1357 MPI_Comm comm=((PetscObject)C)->comm; 1358 PetscMPIInt size,rank,taga,*len_s; 1359 PetscInt *owners,proc,nrows,**buf_ri_k,**nextrow,**nextci; 1360 PetscInt **buf_ri,**buf_rj; 1361 PetscInt cnz=0,*bj_i,*bi,*bj,bnz,nextcj; /* bi,bj,ba: local array of C(mpi mat) */ 1362 MPI_Request *s_waits,*r_waits; 1363 MPI_Status *status; 1364 MatScalar **abuf_r,*ba_i,*pA,*coa,*ba; 1365 PetscInt *ai,*aj,*coi,*coj; 1366 PetscInt *poJ,*pdJ; 1367 Mat A_loc; 1368 Mat_SeqAIJ *a_loc; 1369 1370 PetscFunctionBegin; 1371 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 1372 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 1373 1374 ptap = c->ptap; 1375 merge = ptap->merge; 1376 1377 /* 2) compute numeric C_seq = P_loc^T*A_loc */ 1378 /*------------------------------------------*/ 1379 /* get data from symbolic products */ 1380 coi = merge->coi; coj = merge->coj; 1381 ierr = PetscMalloc((coi[pon]+1)*sizeof(MatScalar),&coa);CHKERRQ(ierr); 1382 ierr = PetscMemzero(coa,coi[pon]*sizeof(MatScalar));CHKERRQ(ierr); 1383 bi = merge->bi; bj = merge->bj; 1384 owners = merge->rowmap->range; 1385 ierr = PetscMalloc((bi[cm]+1)*sizeof(MatScalar),&ba);CHKERRQ(ierr); 1386 ierr = PetscMemzero(ba,bi[cm]*sizeof(MatScalar));CHKERRQ(ierr); 1387 1388 /* get A_loc by taking all local rows of A */ 1389 A_loc = ptap->A_loc; 1390 ierr = MatMPIAIJGetLocalMat(A,MAT_REUSE_MATRIX,&A_loc);CHKERRQ(ierr); 1391 a_loc = (Mat_SeqAIJ*)(A_loc)->data; 1392 ai = a_loc->i; 1393 aj = a_loc->j; 1394 1395 for (i=0; i<am; i++) { 1396 anz = ai[i+1] - ai[i]; 1397 adj = aj + ai[i]; 1398 ada = a_loc->a + ai[i]; 1399 1400 /* 2-b) Compute Cseq = P_loc[i,:]^T*A[i,:] using outer product */ 1401 /*-------------------------------------------------------------*/ 1402 /* put the value into Co=(p->B)^T*A (off-diagonal part, send to others) */ 1403 pnz = po->i[i+1] - po->i[i]; 1404 poJ = po->j + po->i[i]; 1405 pA = po->a + po->i[i]; 1406 for (j=0; j<pnz; j++){ 1407 row = poJ[j]; 1408 cj = coj + coi[row]; 1409 ca = coa + coi[row]; 1410 /* perform sparse axpy */ 1411 nexta = 0; 1412 valtmp = pA[j]; 1413 for (k=0; nexta<anz; k++) { 1414 if (cj[k] == adj[nexta]){ 1415 ca[k] += valtmp*ada[nexta]; 1416 nexta++; 1417 } 1418 } 1419 ierr = PetscLogFlops(2.0*anz);CHKERRQ(ierr); 1420 } 1421 1422 /* put the value into Cd (diagonal part) */ 1423 pnz = pd->i[i+1] - pd->i[i]; 1424 pdJ = pd->j + pd->i[i]; 1425 pA = pd->a + pd->i[i]; 1426 for (j=0; j<pnz; j++){ 1427 row = pdJ[j]; 1428 cj = bj + bi[row]; 1429 ca = ba + bi[row]; 1430 /* perform sparse axpy */ 1431 nexta = 0; 1432 valtmp = pA[j]; 1433 for (k=0; nexta<anz; k++) { 1434 if (cj[k] == adj[nexta]){ 1435 ca[k] += valtmp*ada[nexta]; 1436 nexta++; 1437 } 1438 } 1439 ierr = PetscLogFlops(2.0*anz);CHKERRQ(ierr); 1440 } 1441 } 1442 1443 /* 3) send and recv matrix values coa */ 1444 /*------------------------------------*/ 1445 buf_ri = merge->buf_ri; 1446 buf_rj = merge->buf_rj; 1447 len_s = merge->len_s; 1448 ierr = PetscCommGetNewTag(comm,&taga);CHKERRQ(ierr); 1449 ierr = PetscPostIrecvScalar(comm,taga,merge->nrecv,merge->id_r,merge->len_r,&abuf_r,&r_waits);CHKERRQ(ierr); 1450 1451 ierr = PetscMalloc2(merge->nsend+1,MPI_Request,&s_waits,size,MPI_Status,&status);CHKERRQ(ierr); 1452 for (proc=0,k=0; proc<size; proc++){ 1453 if (!len_s[proc]) continue; 1454 i = merge->owners_co[proc]; 1455 ierr = MPI_Isend(coa+coi[i],len_s[proc],MPIU_MATSCALAR,proc,taga,comm,s_waits+k);CHKERRQ(ierr); 1456 k++; 1457 } 1458 if (merge->nrecv) {ierr = MPI_Waitall(merge->nrecv,r_waits,status);CHKERRQ(ierr);} 1459 if (merge->nsend) {ierr = MPI_Waitall(merge->nsend,s_waits,status);CHKERRQ(ierr);} 1460 1461 ierr = PetscFree2(s_waits,status);CHKERRQ(ierr); 1462 ierr = PetscFree(r_waits);CHKERRQ(ierr); 1463 ierr = PetscFree(coa);CHKERRQ(ierr); 1464 1465 /* 4) insert local Cseq and received values into Cmpi */ 1466 /*----------------------------------------------------*/ 1467 ierr = PetscMalloc3(merge->nrecv,PetscInt**,&buf_ri_k,merge->nrecv,PetscInt*,&nextrow,merge->nrecv,PetscInt*,&nextci);CHKERRQ(ierr); 1468 for (k=0; k<merge->nrecv; k++){ 1469 buf_ri_k[k] = buf_ri[k]; /* beginning of k-th recved i-structure */ 1470 nrows = *(buf_ri_k[k]); 1471 nextrow[k] = buf_ri_k[k]+1; /* next row number of k-th recved i-structure */ 1472 nextci[k] = buf_ri_k[k] + (nrows + 1);/* poins to the next i-structure of k-th recved i-structure */ 1473 } 1474 1475 for (i=0; i<cm; i++) { 1476 row = owners[rank] + i; /* global row index of C_seq */ 1477 bj_i = bj + bi[i]; /* col indices of the i-th row of C */ 1478 ba_i = ba + bi[i]; 1479 bnz = bi[i+1] - bi[i]; 1480 /* add received vals into ba */ 1481 for (k=0; k<merge->nrecv; k++){ /* k-th received message */ 1482 /* i-th row */ 1483 if (i == *nextrow[k]) { 1484 cnz = *(nextci[k]+1) - *nextci[k]; 1485 cj = buf_rj[k] + *(nextci[k]); 1486 ca = abuf_r[k] + *(nextci[k]); 1487 nextcj = 0; 1488 for (j=0; nextcj<cnz; j++){ 1489 if (bj_i[j] == cj[nextcj]){ /* bcol == ccol */ 1490 ba_i[j] += ca[nextcj++]; 1491 } 1492 } 1493 nextrow[k]++; nextci[k]++; 1494 ierr = PetscLogFlops(2.0*cnz);CHKERRQ(ierr); 1495 } 1496 } 1497 ierr = MatSetValues(C,1,&row,bnz,bj_i,ba_i,INSERT_VALUES);CHKERRQ(ierr); 1498 } 1499 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1500 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1501 1502 ierr = PetscFree(ba);CHKERRQ(ierr); 1503 ierr = PetscFree(abuf_r[0]);CHKERRQ(ierr); 1504 ierr = PetscFree(abuf_r);CHKERRQ(ierr); 1505 ierr = PetscFree3(buf_ri_k,nextrow,nextci);CHKERRQ(ierr); 1506 PetscFunctionReturn(0); 1507 } 1508 1509 /* This routine is modified from MatPtAPSymbolic_MPIAIJ_MPIAIJ() */ 1510 #undef __FUNCT__ 1511 #define __FUNCT__ "MatTransposeMatMultSymbolic_MPIAIJ_MPIAIJ_Scalable" 1512 PetscErrorCode MatTransposeMatMultSymbolic_MPIAIJ_MPIAIJ_Scalable(Mat P,Mat A,PetscReal fill,Mat *C) 1513 { 1514 PetscErrorCode ierr; 1515 Mat Cmpi,A_loc,POt,PDt; 1516 Mat_PtAPMPI *ptap; 1517 PetscFreeSpaceList free_space=PETSC_NULL,current_space=PETSC_NULL; 1518 Mat_MPIAIJ *p=(Mat_MPIAIJ*)P->data,*c; 1519 PetscInt *pdti,*pdtj,*poti,*potj,*ptJ; 1520 PetscInt nnz; 1521 PetscInt *lnk,*owners_co,*coi,*coj,i,k,pnz,row; 1522 PetscInt am=A->rmap->n,pn=P->cmap->n; 1523 MPI_Comm comm=((PetscObject)A)->comm; 1524 PetscMPIInt size,rank,tagi,tagj,*len_si,*len_s,*len_ri; 1525 PetscInt **buf_rj,**buf_ri,**buf_ri_k; 1526 PetscInt len,proc,*dnz,*onz,*owners; 1527 PetscInt nzi,*bi,*bj; 1528 PetscInt nrows,*buf_s,*buf_si,*buf_si_i,**nextrow,**nextci; 1529 MPI_Request *swaits,*rwaits; 1530 MPI_Status *sstatus,rstatus; 1531 Mat_Merge_SeqsToMPI *merge; 1532 PetscInt *ai,*aj,*Jptr,anz,*prmap=p->garray,pon,nspacedouble=0,j; 1533 PetscReal afill=1.0,afill_tmp; 1534 PetscInt rstart = P->cmap->rstart,rmax,aN=A->cmap->N,Crmax; 1535 PetscScalar *vals; 1536 Mat_SeqAIJ *a_loc, *pdt,*pot; 1537 1538 PetscFunctionBegin; 1539 /* check if matrix local sizes are compatible */ 1540 if (A->rmap->rstart != P->rmap->rstart || A->rmap->rend != P->rmap->rend){ 1541 SETERRQ4(comm,PETSC_ERR_ARG_SIZ,"Matrix local dimensions are incompatible, A (%D, %D) != P (%D,%D)",A->rmap->rstart,A->rmap->rend,P->rmap->rstart,P->rmap->rend); 1542 } 1543 1544 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 1545 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 1546 1547 /* create struct Mat_PtAPMPI and attached it to C later */ 1548 ierr = PetscNew(Mat_PtAPMPI,&ptap);CHKERRQ(ierr); 1549 1550 /* get A_loc by taking all local rows of A */ 1551 ierr = MatMPIAIJGetLocalMat(A,MAT_INITIAL_MATRIX,&A_loc);CHKERRQ(ierr); 1552 ptap->A_loc = A_loc; 1553 a_loc = (Mat_SeqAIJ*)(A_loc)->data; 1554 ai = a_loc->i; 1555 aj = a_loc->j; 1556 1557 /* determine symbolic Co=(p->B)^T*A - send to others */ 1558 /*----------------------------------------------------*/ 1559 ierr = MatTransposeSymbolic_SeqAIJ(p->A,&PDt);CHKERRQ(ierr); 1560 pdt = (Mat_SeqAIJ*)PDt->data; 1561 pdti = pdt->i; pdtj = pdt->j; 1562 1563 ierr = MatTransposeSymbolic_SeqAIJ(p->B,&POt);CHKERRQ(ierr); 1564 pot = (Mat_SeqAIJ*)POt->data; 1565 poti = pot->i; potj = pot->j; 1566 1567 /* then, compute symbolic Co = (p->B)^T*A */ 1568 pon = (p->B)->cmap->n; /* total num of rows to be sent to other processors 1569 >= (num of nonzero rows of C_seq) - pn */ 1570 ierr = PetscMalloc((pon+1)*sizeof(PetscInt),&coi);CHKERRQ(ierr); 1571 coi[0] = 0; 1572 1573 /* set initial free space to be fill*(nnz(p->B) + nnz(A)) */ 1574 nnz = fill*(poti[pon] + ai[am]); 1575 ierr = PetscFreeSpaceGet(nnz,&free_space);CHKERRQ(ierr); 1576 current_space = free_space; 1577 1578 /* create and initialize a linked list */ 1579 i = PetscMax(pdt->rmax,pot->rmax); 1580 Crmax = i*a_loc->rmax*size; /* non-scalable! */ 1581 if (!Crmax || Crmax > aN) Crmax = aN; 1582 ierr = PetscLLCondensedCreate_Scalable(Crmax,&lnk);CHKERRQ(ierr); 1583 1584 for (i=0; i<pon; i++) { 1585 pnz = poti[i+1] - poti[i]; 1586 ptJ = potj + poti[i]; 1587 for (j=0; j<pnz; j++){ 1588 row = ptJ[j]; /* row of A_loc == col of Pot */ 1589 anz = ai[row+1] - ai[row]; 1590 Jptr = aj + ai[row]; 1591 /* add non-zero cols of AP into the sorted linked list lnk */ 1592 ierr = PetscLLCondensedAddSorted_Scalable(anz,Jptr,lnk);CHKERRQ(ierr); 1593 } 1594 nnz = lnk[0]; 1595 1596 /* If free space is not available, double the total space in the list */ 1597 if (current_space->local_remaining<nnz) { 1598 ierr = PetscFreeSpaceGet(nnz+current_space->total_array_size,¤t_space);CHKERRQ(ierr); 1599 nspacedouble++; 1600 } 1601 1602 /* Copy data into free space, and zero out denserows */ 1603 ierr = PetscLLCondensedClean_Scalable(nnz,current_space->array,lnk);CHKERRQ(ierr); 1604 current_space->array += nnz; 1605 current_space->local_used += nnz; 1606 current_space->local_remaining -= nnz; 1607 coi[i+1] = coi[i] + nnz; 1608 } 1609 1610 ierr = PetscMalloc((coi[pon]+1)*sizeof(PetscInt),&coj);CHKERRQ(ierr); 1611 ierr = PetscFreeSpaceContiguous(&free_space,coj);CHKERRQ(ierr); 1612 afill_tmp = (PetscReal)coi[pon]/(poti[pon] + ai[am]+1); 1613 if (afill_tmp > afill) afill = afill_tmp; 1614 1615 /* send j-array (coj) of Co to other processors */ 1616 /*----------------------------------------------*/ 1617 /* determine row ownership */ 1618 ierr = PetscNew(Mat_Merge_SeqsToMPI,&merge);CHKERRQ(ierr); 1619 ierr = PetscLayoutCreate(comm,&merge->rowmap);CHKERRQ(ierr); 1620 merge->rowmap->n = pn; 1621 merge->rowmap->bs = 1; 1622 ierr = PetscLayoutSetUp(merge->rowmap);CHKERRQ(ierr); 1623 owners = merge->rowmap->range; 1624 1625 /* determine the number of messages to send, their lengths */ 1626 ierr = PetscMalloc(size*sizeof(PetscMPIInt),&len_si);CHKERRQ(ierr); 1627 ierr = PetscMemzero(len_si,size*sizeof(PetscMPIInt));CHKERRQ(ierr); 1628 ierr = PetscMalloc(size*sizeof(PetscMPIInt),&merge->len_s);CHKERRQ(ierr); 1629 len_s = merge->len_s; 1630 merge->nsend = 0; 1631 1632 ierr = PetscMalloc((size+2)*sizeof(PetscInt),&owners_co);CHKERRQ(ierr); 1633 ierr = PetscMemzero(len_s,size*sizeof(PetscMPIInt));CHKERRQ(ierr); 1634 1635 proc = 0; 1636 for (i=0; i<pon; i++){ 1637 while (prmap[i] >= owners[proc+1]) proc++; 1638 len_si[proc]++; /* num of rows in Co to be sent to [proc] */ 1639 len_s[proc] += coi[i+1] - coi[i]; 1640 } 1641 1642 len = 0; /* max length of buf_si[] */ 1643 owners_co[0] = 0; 1644 for (proc=0; proc<size; proc++){ 1645 owners_co[proc+1] = owners_co[proc] + len_si[proc]; 1646 if (len_si[proc]){ 1647 merge->nsend++; 1648 len_si[proc] = 2*(len_si[proc] + 1); 1649 len += len_si[proc]; 1650 } 1651 } 1652 1653 /* determine the number and length of messages to receive for coi and coj */ 1654 ierr = PetscGatherNumberOfMessages(comm,PETSC_NULL,len_s,&merge->nrecv);CHKERRQ(ierr); 1655 ierr = PetscGatherMessageLengths2(comm,merge->nsend,merge->nrecv,len_s,len_si,&merge->id_r,&merge->len_r,&len_ri);CHKERRQ(ierr); 1656 1657 /* post the Irecv and Isend of coj */ 1658 ierr = PetscCommGetNewTag(comm,&tagj);CHKERRQ(ierr); 1659 ierr = PetscPostIrecvInt(comm,tagj,merge->nrecv,merge->id_r,merge->len_r,&buf_rj,&rwaits);CHKERRQ(ierr); 1660 ierr = PetscMalloc((merge->nsend+1)*sizeof(MPI_Request),&swaits);CHKERRQ(ierr); 1661 for (proc=0, k=0; proc<size; proc++){ 1662 if (!len_s[proc]) continue; 1663 i = owners_co[proc]; 1664 ierr = MPI_Isend(coj+coi[i],len_s[proc],MPIU_INT,proc,tagj,comm,swaits+k);CHKERRQ(ierr); 1665 k++; 1666 } 1667 1668 /* receives and sends of coj are complete */ 1669 ierr = PetscMalloc(size*sizeof(MPI_Status),&sstatus);CHKERRQ(ierr); 1670 for (i=0; i<merge->nrecv; i++){ 1671 PetscMPIInt icompleted; 1672 ierr = MPI_Waitany(merge->nrecv,rwaits,&icompleted,&rstatus);CHKERRQ(ierr); 1673 } 1674 ierr = PetscFree(rwaits);CHKERRQ(ierr); 1675 if (merge->nsend) {ierr = MPI_Waitall(merge->nsend,swaits,sstatus);CHKERRQ(ierr);} 1676 1677 /* send and recv coi */ 1678 /*-------------------*/ 1679 ierr = PetscCommGetNewTag(comm,&tagi);CHKERRQ(ierr); 1680 ierr = PetscPostIrecvInt(comm,tagi,merge->nrecv,merge->id_r,len_ri,&buf_ri,&rwaits);CHKERRQ(ierr); 1681 ierr = PetscMalloc((len+1)*sizeof(PetscInt),&buf_s);CHKERRQ(ierr); 1682 buf_si = buf_s; /* points to the beginning of k-th msg to be sent */ 1683 for (proc=0,k=0; proc<size; proc++){ 1684 if (!len_s[proc]) continue; 1685 /* form outgoing message for i-structure: 1686 buf_si[0]: nrows to be sent 1687 [1:nrows]: row index (global) 1688 [nrows+1:2*nrows+1]: i-structure index 1689 */ 1690 /*-------------------------------------------*/ 1691 nrows = len_si[proc]/2 - 1; 1692 buf_si_i = buf_si + nrows+1; 1693 buf_si[0] = nrows; 1694 buf_si_i[0] = 0; 1695 nrows = 0; 1696 for (i=owners_co[proc]; i<owners_co[proc+1]; i++){ 1697 nzi = coi[i+1] - coi[i]; 1698 buf_si_i[nrows+1] = buf_si_i[nrows] + nzi; /* i-structure */ 1699 buf_si[nrows+1] =prmap[i] -owners[proc]; /* local row index */ 1700 nrows++; 1701 } 1702 ierr = MPI_Isend(buf_si,len_si[proc],MPIU_INT,proc,tagi,comm,swaits+k);CHKERRQ(ierr); 1703 k++; 1704 buf_si += len_si[proc]; 1705 } 1706 i = merge->nrecv; 1707 while (i--) { 1708 PetscMPIInt icompleted; 1709 ierr = MPI_Waitany(merge->nrecv,rwaits,&icompleted,&rstatus);CHKERRQ(ierr); 1710 } 1711 ierr = PetscFree(rwaits);CHKERRQ(ierr); 1712 if (merge->nsend) {ierr = MPI_Waitall(merge->nsend,swaits,sstatus);CHKERRQ(ierr);} 1713 ierr = PetscFree(len_si);CHKERRQ(ierr); 1714 ierr = PetscFree(len_ri);CHKERRQ(ierr); 1715 ierr = PetscFree(swaits);CHKERRQ(ierr); 1716 ierr = PetscFree(sstatus);CHKERRQ(ierr); 1717 ierr = PetscFree(buf_s);CHKERRQ(ierr); 1718 1719 /* compute the local portion of C (mpi mat) */ 1720 /*------------------------------------------*/ 1721 /* allocate bi array and free space for accumulating nonzero column info */ 1722 ierr = PetscMalloc((pn+1)*sizeof(PetscInt),&bi);CHKERRQ(ierr); 1723 bi[0] = 0; 1724 1725 /* set initial free space to be fill*(nnz(P) + nnz(AP)) */ 1726 nnz = fill*(pdti[pn] + poti[pon] + ai[am]); 1727 ierr = PetscFreeSpaceGet(nnz,&free_space);CHKERRQ(ierr); 1728 current_space = free_space; 1729 1730 ierr = PetscMalloc3(merge->nrecv,PetscInt**,&buf_ri_k,merge->nrecv,PetscInt*,&nextrow,merge->nrecv,PetscInt*,&nextci);CHKERRQ(ierr); 1731 for (k=0; k<merge->nrecv; k++){ 1732 buf_ri_k[k] = buf_ri[k]; /* beginning of k-th recved i-structure */ 1733 nrows = *buf_ri_k[k]; 1734 nextrow[k] = buf_ri_k[k] + 1; /* next row number of k-th recved i-structure */ 1735 nextci[k] = buf_ri_k[k] + (nrows + 1);/* poins to the next i-structure of k-th recved i-structure */ 1736 } 1737 1738 ierr = MatPreallocateInitialize(comm,pn,A->cmap->n,dnz,onz);CHKERRQ(ierr); 1739 rmax = 0; 1740 for (i=0; i<pn; i++) { 1741 /* add pdt[i,:]*AP into lnk */ 1742 pnz = pdti[i+1] - pdti[i]; 1743 ptJ = pdtj + pdti[i]; 1744 for (j=0; j<pnz; j++){ 1745 row = ptJ[j]; /* row of AP == col of Pt */ 1746 anz = ai[row+1] - ai[row]; 1747 Jptr = aj + ai[row]; 1748 /* add non-zero cols of AP into the sorted linked list lnk */ 1749 ierr = PetscLLCondensedAddSorted_Scalable(anz,Jptr,lnk);CHKERRQ(ierr); 1750 } 1751 1752 /* add received col data into lnk */ 1753 for (k=0; k<merge->nrecv; k++){ /* k-th received message */ 1754 if (i == *nextrow[k]) { /* i-th row */ 1755 nzi = *(nextci[k]+1) - *nextci[k]; 1756 Jptr = buf_rj[k] + *nextci[k]; 1757 ierr = PetscLLCondensedAddSorted_Scalable(nzi,Jptr,lnk);CHKERRQ(ierr); 1758 nextrow[k]++; nextci[k]++; 1759 } 1760 } 1761 nnz = lnk[0]; 1762 1763 /* if free space is not available, make more free space */ 1764 if (current_space->local_remaining<nnz) { 1765 ierr = PetscFreeSpaceGet(nnz+current_space->total_array_size,¤t_space);CHKERRQ(ierr); 1766 nspacedouble++; 1767 } 1768 /* copy data into free space, then initialize lnk */ 1769 ierr = PetscLLCondensedClean_Scalable(nnz,current_space->array,lnk);CHKERRQ(ierr); 1770 ierr = MatPreallocateSet(i+owners[rank],nnz,current_space->array,dnz,onz);CHKERRQ(ierr); 1771 current_space->array += nnz; 1772 current_space->local_used += nnz; 1773 current_space->local_remaining -= nnz; 1774 bi[i+1] = bi[i] + nnz; 1775 if (nnz > rmax) rmax = nnz; 1776 } 1777 ierr = PetscFree3(buf_ri_k,nextrow,nextci);CHKERRQ(ierr); 1778 1779 ierr = PetscMalloc((bi[pn]+1)*sizeof(PetscInt),&bj);CHKERRQ(ierr); 1780 ierr = PetscFreeSpaceContiguous(&free_space,bj);CHKERRQ(ierr); 1781 afill_tmp = (PetscReal)bi[pn]/(pdti[pn] + poti[pon] + ai[am]+1); 1782 if (afill_tmp > afill) afill = afill_tmp; 1783 ierr = PetscLLCondensedDestroy_Scalable(lnk);CHKERRQ(ierr); 1784 ierr = MatDestroy(&POt);CHKERRQ(ierr); 1785 ierr = MatDestroy(&PDt);CHKERRQ(ierr); 1786 1787 /* create symbolic parallel matrix Cmpi - why cannot be assembled in Numeric part */ 1788 /*----------------------------------------------------------------------------------*/ 1789 ierr = PetscMalloc((rmax+1)*sizeof(PetscScalar),&vals);CHKERRQ(ierr); 1790 ierr = PetscMemzero(vals,rmax*sizeof(PetscScalar));CHKERRQ(ierr); 1791 1792 ierr = MatCreate(comm,&Cmpi);CHKERRQ(ierr); 1793 ierr = MatSetSizes(Cmpi,pn,A->cmap->n,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 1794 ierr = MatSetBlockSizes(Cmpi,P->cmap->bs,A->cmap->bs);CHKERRQ(ierr); 1795 ierr = MatSetType(Cmpi,MATMPIAIJ);CHKERRQ(ierr); 1796 ierr = MatMPIAIJSetPreallocation(Cmpi,0,dnz,0,onz);CHKERRQ(ierr); 1797 ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr); 1798 ierr = MatSetBlockSize(Cmpi,1);CHKERRQ(ierr); 1799 for (i=0; i<pn; i++){ 1800 row = i + rstart; 1801 nnz = bi[i+1] - bi[i]; 1802 Jptr = bj + bi[i]; 1803 ierr = MatSetValues(Cmpi,1,&row,nnz,Jptr,vals,INSERT_VALUES);CHKERRQ(ierr); 1804 } 1805 ierr = MatAssemblyBegin(Cmpi,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1806 ierr = MatAssemblyEnd(Cmpi,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1807 ierr = PetscFree(vals);CHKERRQ(ierr); 1808 1809 merge->bi = bi; 1810 merge->bj = bj; 1811 merge->coi = coi; 1812 merge->coj = coj; 1813 merge->buf_ri = buf_ri; 1814 merge->buf_rj = buf_rj; 1815 merge->owners_co = owners_co; 1816 merge->destroy = Cmpi->ops->destroy; 1817 merge->duplicate = Cmpi->ops->duplicate; 1818 1819 Cmpi->ops->mattransposemultnumeric = MatTransposeMatMultNumeric_MPIAIJ_MPIAIJ_Scalable; 1820 Cmpi->ops->destroy = MatDestroy_MPIAIJ_PtAP; 1821 1822 /* attach the supporting struct to Cmpi for reuse */ 1823 c = (Mat_MPIAIJ*)Cmpi->data; 1824 c->ptap = ptap; 1825 ptap->api = PETSC_NULL; 1826 ptap->apj = PETSC_NULL; 1827 ptap->merge = merge; 1828 ptap->rmax = rmax; 1829 ptap->apa = PETSC_NULL; 1830 1831 *C = Cmpi; 1832 #if defined(PETSC_USE_INFO) 1833 if (bi[pn] != 0) { 1834 ierr = PetscInfo3(Cmpi,"Reallocs %D; Fill ratio: given %G needed %G.\n",nspacedouble,fill,afill);CHKERRQ(ierr); 1835 ierr = PetscInfo1(Cmpi,"Use MatTransposeMatMult(A,B,MatReuse,%G,&C) for best performance.\n",afill);CHKERRQ(ierr); 1836 } else { 1837 ierr = PetscInfo(Cmpi,"Empty matrix product\n");CHKERRQ(ierr); 1838 } 1839 #endif 1840 PetscFunctionReturn(0); 1841 } 1842