1 2 /* 3 Defines projective product routines where A is a MPIAIJ matrix 4 C = P^T * A * P 5 */ 6 7 #include <../src/mat/impls/aij/seq/aij.h> /*I "petscmat.h" I*/ 8 #include <../src/mat/utils/freespace.h> 9 #include <../src/mat/impls/aij/mpi/mpiaij.h> 10 #include <petscbt.h> 11 #include <petsctime.h> 12 #include <petsc/private/hashmapiv.h> 13 #include <petsc/private/hashseti.h> 14 #include <petscsf.h> 15 16 /* #define PTAP_PROFILE */ 17 18 PetscErrorCode MatView_MPIAIJ_PtAP(Mat A,PetscViewer viewer) 19 { 20 PetscErrorCode ierr; 21 Mat_MPIAIJ *a=(Mat_MPIAIJ*)A->data; 22 Mat_APMPI *ptap=a->ap; 23 PetscBool iascii; 24 PetscViewerFormat format; 25 26 PetscFunctionBegin; 27 if (!ptap) { 28 /* hack: MatDuplicate() sets oldmat->ops->view to newmat which is a base mat class with null ptpa! */ 29 A->ops->view = MatView_MPIAIJ; 30 ierr = (A->ops->view)(A,viewer);CHKERRQ(ierr); 31 PetscFunctionReturn(0); 32 } 33 34 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 35 if (iascii) { 36 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 37 if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) { 38 if (ptap->algType == 0) { 39 ierr = PetscViewerASCIIPrintf(viewer,"using scalable MatPtAP() implementation\n");CHKERRQ(ierr); 40 } else if (ptap->algType == 1) { 41 ierr = PetscViewerASCIIPrintf(viewer,"using nonscalable MatPtAP() implementation\n");CHKERRQ(ierr); 42 } else if (ptap->algType == 2) { 43 ierr = PetscViewerASCIIPrintf(viewer,"using allatonce MatPtAP() implementation\n");CHKERRQ(ierr); 44 } else if (ptap->algType == 3) { 45 ierr = PetscViewerASCIIPrintf(viewer,"using merged allatonce MatPtAP() implementation\n");CHKERRQ(ierr); 46 } 47 } 48 } 49 ierr = (ptap->view)(A,viewer);CHKERRQ(ierr); 50 PetscFunctionReturn(0); 51 } 52 53 PetscErrorCode MatFreeIntermediateDataStructures_MPIAIJ_AP(Mat A) 54 { 55 PetscErrorCode ierr; 56 Mat_MPIAIJ *a=(Mat_MPIAIJ*)A->data; 57 Mat_APMPI *ptap=a->ap; 58 Mat_Merge_SeqsToMPI *merge; 59 60 PetscFunctionBegin; 61 if (!ptap) PetscFunctionReturn(0); 62 63 ierr = PetscFree2(ptap->startsj_s,ptap->startsj_r);CHKERRQ(ierr); 64 ierr = PetscFree(ptap->bufa);CHKERRQ(ierr); 65 ierr = MatDestroy(&ptap->P_loc);CHKERRQ(ierr); 66 ierr = MatDestroy(&ptap->P_oth);CHKERRQ(ierr); 67 ierr = MatDestroy(&ptap->A_loc);CHKERRQ(ierr); /* used by MatTransposeMatMult() */ 68 ierr = MatDestroy(&ptap->Rd);CHKERRQ(ierr); 69 ierr = MatDestroy(&ptap->Ro);CHKERRQ(ierr); 70 if (ptap->AP_loc) { /* used by alg_rap */ 71 Mat_SeqAIJ *ap = (Mat_SeqAIJ*)(ptap->AP_loc)->data; 72 ierr = PetscFree(ap->i);CHKERRQ(ierr); 73 ierr = PetscFree2(ap->j,ap->a);CHKERRQ(ierr); 74 ierr = MatDestroy(&ptap->AP_loc);CHKERRQ(ierr); 75 } else { /* used by alg_ptap */ 76 ierr = PetscFree(ptap->api);CHKERRQ(ierr); 77 ierr = PetscFree(ptap->apj);CHKERRQ(ierr); 78 } 79 ierr = MatDestroy(&ptap->C_loc);CHKERRQ(ierr); 80 ierr = MatDestroy(&ptap->C_oth);CHKERRQ(ierr); 81 if (ptap->apa) {ierr = PetscFree(ptap->apa);CHKERRQ(ierr);} 82 83 ierr = MatDestroy(&ptap->Pt);CHKERRQ(ierr); 84 85 merge=ptap->merge; 86 if (merge) { /* used by alg_ptap */ 87 ierr = PetscFree(merge->id_r);CHKERRQ(ierr); 88 ierr = PetscFree(merge->len_s);CHKERRQ(ierr); 89 ierr = PetscFree(merge->len_r);CHKERRQ(ierr); 90 ierr = PetscFree(merge->bi);CHKERRQ(ierr); 91 ierr = PetscFree(merge->bj);CHKERRQ(ierr); 92 ierr = PetscFree(merge->buf_ri[0]);CHKERRQ(ierr); 93 ierr = PetscFree(merge->buf_ri);CHKERRQ(ierr); 94 ierr = PetscFree(merge->buf_rj[0]);CHKERRQ(ierr); 95 ierr = PetscFree(merge->buf_rj);CHKERRQ(ierr); 96 ierr = PetscFree(merge->coi);CHKERRQ(ierr); 97 ierr = PetscFree(merge->coj);CHKERRQ(ierr); 98 ierr = PetscFree(merge->owners_co);CHKERRQ(ierr); 99 ierr = PetscLayoutDestroy(&merge->rowmap);CHKERRQ(ierr); 100 ierr = PetscFree(ptap->merge);CHKERRQ(ierr); 101 } 102 ierr = ISLocalToGlobalMappingDestroy(&ptap->ltog);CHKERRQ(ierr); 103 104 ierr = PetscSFDestroy(&ptap->sf);CHKERRQ(ierr); 105 ierr = PetscFree(ptap->c_othi);CHKERRQ(ierr); 106 ierr = PetscFree(ptap->c_rmti);CHKERRQ(ierr); 107 PetscFunctionReturn(0); 108 } 109 110 PetscErrorCode MatDestroy_MPIAIJ_PtAP(Mat A) 111 { 112 PetscErrorCode ierr; 113 Mat_MPIAIJ *a=(Mat_MPIAIJ*)A->data; 114 Mat_APMPI *ptap=a->ap; 115 116 PetscFunctionBegin; 117 ierr = (*A->ops->freeintermediatedatastructures)(A);CHKERRQ(ierr); 118 ierr = (*ptap->destroy)(A);CHKERRQ(ierr); /* MatDestroy_MPIAIJ(A) */ 119 ierr = PetscFree(ptap);CHKERRQ(ierr); 120 PetscFunctionReturn(0); 121 } 122 123 PETSC_INTERN PetscErrorCode MatPtAP_MPIAIJ_MPIAIJ(Mat A,Mat P,MatReuse scall,PetscReal fill,Mat *C) 124 { 125 PetscErrorCode ierr; 126 PetscBool flg; 127 MPI_Comm comm; 128 #if !defined(PETSC_HAVE_HYPRE) 129 const char *algTypes[4] = {"scalable","nonscalable","allatonce","allatonce_merged"}; 130 PetscInt nalg=4; 131 #else 132 const char *algTypes[5] = {"scalable","nonscalable","allatonce","allatonce_merged","hypre"}; 133 PetscInt nalg=5; 134 #endif 135 PetscInt pN=P->cmap->N,alg=1; /* set default algorithm */ 136 137 PetscFunctionBegin; 138 /* check if matrix local sizes are compatible */ 139 ierr = PetscObjectGetComm((PetscObject)A,&comm);CHKERRQ(ierr); 140 if (A->rmap->rstart != P->rmap->rstart || A->rmap->rend != P->rmap->rend) SETERRQ4(comm,PETSC_ERR_ARG_SIZ,"Matrix local dimensions are incompatible, Arow (%D, %D) != Prow (%D,%D)",A->rmap->rstart,A->rmap->rend,P->rmap->rstart,P->rmap->rend); 141 if (A->cmap->rstart != P->rmap->rstart || A->cmap->rend != P->rmap->rend) SETERRQ4(comm,PETSC_ERR_ARG_SIZ,"Matrix local dimensions are incompatible, Acol (%D, %D) != Prow (%D,%D)",A->cmap->rstart,A->cmap->rend,P->rmap->rstart,P->rmap->rend); 142 143 if (scall == MAT_INITIAL_MATRIX) { 144 /* pick an algorithm */ 145 ierr = PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"MatPtAP","Mat");CHKERRQ(ierr); 146 ierr = PetscOptionsEList("-matptap_via","Algorithmic approach","MatPtAP",algTypes,nalg,algTypes[alg],&alg,&flg);CHKERRQ(ierr); 147 ierr = PetscOptionsEnd();CHKERRQ(ierr); 148 149 if (!flg && pN > 100000) { /* may switch to scalable algorithm as default */ 150 MatInfo Ainfo,Pinfo; 151 PetscInt nz_local; 152 PetscBool alg_scalable_loc=PETSC_FALSE,alg_scalable; 153 154 ierr = MatGetInfo(A,MAT_LOCAL,&Ainfo);CHKERRQ(ierr); 155 ierr = MatGetInfo(P,MAT_LOCAL,&Pinfo);CHKERRQ(ierr); 156 nz_local = (PetscInt)(Ainfo.nz_allocated + Pinfo.nz_allocated); 157 158 if (pN > fill*nz_local) alg_scalable_loc = PETSC_TRUE; 159 ierr = MPIU_Allreduce(&alg_scalable_loc,&alg_scalable,1,MPIU_BOOL,MPI_LOR,comm);CHKERRQ(ierr); 160 161 if (alg_scalable) { 162 alg = 0; /* scalable algorithm would 50% slower than nonscalable algorithm */ 163 } 164 } 165 166 switch (alg) { 167 case 1: 168 /* do R=P^T locally, then C=R*A*P -- nonscalable */ 169 ierr = PetscLogEventBegin(MAT_PtAPSymbolic,A,P,0,0);CHKERRQ(ierr); 170 ierr = MatPtAPSymbolic_MPIAIJ_MPIAIJ(A,P,fill,C);CHKERRQ(ierr); 171 ierr = PetscLogEventEnd(MAT_PtAPSymbolic,A,P,0,0);CHKERRQ(ierr); 172 break; 173 case 2: 174 /* compute C=P^T*A*P allatonce */ 175 ierr = PetscLogEventBegin(MAT_PtAPSymbolic,A,P,0,0);CHKERRQ(ierr); 176 ierr = MatPtAPSymbolic_MPIAIJ_MPIAIJ_allatonce(A,P,fill,C);CHKERRQ(ierr); 177 ierr = PetscLogEventEnd(MAT_PtAPSymbolic,A,P,0,0);CHKERRQ(ierr); 178 break; 179 case 3: 180 /* compute C=P^T*A*P allatonce */ 181 ierr = PetscLogEventBegin(MAT_PtAPSymbolic,A,P,0,0);CHKERRQ(ierr); 182 ierr = MatPtAPSymbolic_MPIAIJ_MPIAIJ_allatonce_merged(A,P,fill,C);CHKERRQ(ierr); 183 ierr = PetscLogEventEnd(MAT_PtAPSymbolic,A,P,0,0);CHKERRQ(ierr); 184 break; 185 #if defined(PETSC_HAVE_HYPRE) 186 case 4: 187 /* Use boomerAMGBuildCoarseOperator */ 188 ierr = PetscLogEventBegin(MAT_PtAPSymbolic,A,P,0,0);CHKERRQ(ierr); 189 ierr = MatPtAPSymbolic_AIJ_AIJ_wHYPRE(A,P,fill,C);CHKERRQ(ierr); 190 ierr = PetscLogEventEnd(MAT_PtAPSymbolic,A,P,0,0);CHKERRQ(ierr); 191 break; 192 #endif 193 default: 194 /* do R=P^T locally, then C=R*A*P */ 195 ierr = PetscLogEventBegin(MAT_PtAPSymbolic,A,P,0,0);CHKERRQ(ierr); 196 ierr = MatPtAPSymbolic_MPIAIJ_MPIAIJ_scalable(A,P,fill,C);CHKERRQ(ierr); 197 ierr = PetscLogEventEnd(MAT_PtAPSymbolic,A,P,0,0);CHKERRQ(ierr); 198 break; 199 } 200 201 if (alg == 0 || alg == 1 || alg == 2 || alg == 3) { 202 Mat_MPIAIJ *c = (Mat_MPIAIJ*)(*C)->data; 203 Mat_APMPI *ap = c->ap; 204 ierr = PetscOptionsBegin(PetscObjectComm((PetscObject)(*C)),((PetscObject)(*C))->prefix,"MatFreeIntermediateDataStructures","Mat");CHKERRQ(ierr); 205 ap->freestruct = PETSC_FALSE; 206 ierr = PetscOptionsBool("-mat_freeintermediatedatastructures","Free intermediate data structures", "MatFreeIntermediateDataStructures",ap->freestruct,&ap->freestruct, NULL);CHKERRQ(ierr); 207 ierr = PetscOptionsEnd();CHKERRQ(ierr); 208 } 209 } 210 211 ierr = PetscLogEventBegin(MAT_PtAPNumeric,A,P,0,0);CHKERRQ(ierr); 212 ierr = (*(*C)->ops->ptapnumeric)(A,P,*C);CHKERRQ(ierr); 213 ierr = PetscLogEventEnd(MAT_PtAPNumeric,A,P,0,0);CHKERRQ(ierr); 214 PetscFunctionReturn(0); 215 } 216 217 PetscErrorCode MatPtAPNumeric_MPIAIJ_MPIAIJ_scalable(Mat A,Mat P,Mat C) 218 { 219 PetscErrorCode ierr; 220 Mat_MPIAIJ *a=(Mat_MPIAIJ*)A->data,*p=(Mat_MPIAIJ*)P->data,*c=(Mat_MPIAIJ*)C->data; 221 Mat_SeqAIJ *ad=(Mat_SeqAIJ*)(a->A)->data,*ao=(Mat_SeqAIJ*)(a->B)->data; 222 Mat_SeqAIJ *ap,*p_loc,*p_oth,*c_seq; 223 Mat_APMPI *ptap = c->ap; 224 Mat AP_loc,C_loc,C_oth; 225 PetscInt i,rstart,rend,cm,ncols,row,*api,*apj,am = A->rmap->n,apnz,nout; 226 PetscScalar *apa; 227 const PetscInt *cols; 228 const PetscScalar *vals; 229 230 PetscFunctionBegin; 231 if (!ptap) { 232 MPI_Comm comm; 233 ierr = PetscObjectGetComm((PetscObject)C,&comm);CHKERRQ(ierr); 234 SETERRQ(comm,PETSC_ERR_ARG_WRONGSTATE,"PtAP cannot be reused. Do not call MatFreeIntermediateDataStructures() or use '-mat_freeintermediatedatastructures'"); 235 } 236 237 ierr = MatZeroEntries(C);CHKERRQ(ierr); 238 239 /* 1) get R = Pd^T,Ro = Po^T */ 240 if (ptap->reuse == MAT_REUSE_MATRIX) { 241 ierr = MatTranspose(p->A,MAT_REUSE_MATRIX,&ptap->Rd);CHKERRQ(ierr); 242 ierr = MatTranspose(p->B,MAT_REUSE_MATRIX,&ptap->Ro);CHKERRQ(ierr); 243 } 244 245 /* 2) get AP_loc */ 246 AP_loc = ptap->AP_loc; 247 ap = (Mat_SeqAIJ*)AP_loc->data; 248 249 /* 2-1) get P_oth = ptap->P_oth and P_loc = ptap->P_loc */ 250 /*-----------------------------------------------------*/ 251 if (ptap->reuse == MAT_REUSE_MATRIX) { 252 /* P_oth and P_loc are obtained in MatPtASymbolic() when reuse == MAT_INITIAL_MATRIX */ 253 ierr = MatGetBrowsOfAoCols_MPIAIJ(A,P,MAT_REUSE_MATRIX,&ptap->startsj_s,&ptap->startsj_r,&ptap->bufa,&ptap->P_oth);CHKERRQ(ierr); 254 ierr = MatMPIAIJGetLocalMat(P,MAT_REUSE_MATRIX,&ptap->P_loc);CHKERRQ(ierr); 255 } 256 257 /* 2-2) compute numeric A_loc*P - dominating part */ 258 /* ---------------------------------------------- */ 259 /* get data from symbolic products */ 260 p_loc = (Mat_SeqAIJ*)(ptap->P_loc)->data; 261 if (ptap->P_oth) p_oth = (Mat_SeqAIJ*)(ptap->P_oth)->data; 262 else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"ptap->P_oth is NULL. Cannot proceed!"); 263 264 api = ap->i; 265 apj = ap->j; 266 ierr = ISLocalToGlobalMappingApply(ptap->ltog,api[AP_loc->rmap->n],apj,apj);CHKERRQ(ierr); 267 for (i=0; i<am; i++) { 268 /* AP[i,:] = A[i,:]*P = Ad*P_loc Ao*P_oth */ 269 apnz = api[i+1] - api[i]; 270 apa = ap->a + api[i]; 271 ierr = PetscMemzero(apa,sizeof(PetscScalar)*apnz);CHKERRQ(ierr); 272 AProw_scalable(i,ad,ao,p_loc,p_oth,api,apj,apa); 273 ierr = PetscLogFlops(2.0*apnz);CHKERRQ(ierr); 274 } 275 ierr = ISGlobalToLocalMappingApply(ptap->ltog,IS_GTOLM_DROP,api[AP_loc->rmap->n],apj,&nout,apj);CHKERRQ(ierr); 276 if (api[AP_loc->rmap->n] != nout) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Incorrect mapping %D != %D\n",api[AP_loc->rmap->n],nout); 277 278 /* 3) C_loc = Rd*AP_loc, C_oth = Ro*AP_loc */ 279 /* Always use scalable version since we are in the MPI scalable version */ 280 ierr = MatMatMultNumeric_SeqAIJ_SeqAIJ_Scalable(ptap->Rd,AP_loc,ptap->C_loc);CHKERRQ(ierr); 281 ierr = MatMatMultNumeric_SeqAIJ_SeqAIJ_Scalable(ptap->Ro,AP_loc,ptap->C_oth);CHKERRQ(ierr); 282 283 C_loc = ptap->C_loc; 284 C_oth = ptap->C_oth; 285 286 /* add C_loc and Co to to C */ 287 ierr = MatGetOwnershipRange(C,&rstart,&rend);CHKERRQ(ierr); 288 289 /* C_loc -> C */ 290 cm = C_loc->rmap->N; 291 c_seq = (Mat_SeqAIJ*)C_loc->data; 292 cols = c_seq->j; 293 vals = c_seq->a; 294 ierr = ISLocalToGlobalMappingApply(ptap->ltog,c_seq->i[C_loc->rmap->n],c_seq->j,c_seq->j);CHKERRQ(ierr); 295 296 /* The (fast) MatSetValues_MPIAIJ_CopyFromCSRFormat function can only be used when C->was_assembled is PETSC_FALSE and */ 297 /* when there are no off-processor parts. */ 298 /* If was_assembled is true, then the statement aj[rowstart_diag+dnz_row] = mat_j[col] - cstart; in MatSetValues_MPIAIJ_CopyFromCSRFormat */ 299 /* is no longer true. Then the more complex function MatSetValues_MPIAIJ() has to be used, where the column index is looked up from */ 300 /* a table, and other, more complex stuff has to be done. */ 301 if (C->assembled) { 302 C->was_assembled = PETSC_TRUE; 303 C->assembled = PETSC_FALSE; 304 } 305 if (C->was_assembled) { 306 for (i=0; i<cm; i++) { 307 ncols = c_seq->i[i+1] - c_seq->i[i]; 308 row = rstart + i; 309 ierr = MatSetValues_MPIAIJ(C,1,&row,ncols,cols,vals,ADD_VALUES);CHKERRQ(ierr); 310 cols += ncols; vals += ncols; 311 } 312 } else { 313 ierr = MatSetValues_MPIAIJ_CopyFromCSRFormat(C,c_seq->j,c_seq->i,c_seq->a);CHKERRQ(ierr); 314 } 315 ierr = ISGlobalToLocalMappingApply(ptap->ltog,IS_GTOLM_DROP,c_seq->i[C_loc->rmap->n],c_seq->j,&nout,c_seq->j);CHKERRQ(ierr); 316 if (c_seq->i[C_loc->rmap->n] != nout) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Incorrect mapping %D != %D\n",c_seq->i[C_loc->rmap->n],nout); 317 318 /* Co -> C, off-processor part */ 319 cm = C_oth->rmap->N; 320 c_seq = (Mat_SeqAIJ*)C_oth->data; 321 cols = c_seq->j; 322 vals = c_seq->a; 323 ierr = ISLocalToGlobalMappingApply(ptap->ltog,c_seq->i[C_oth->rmap->n],c_seq->j,c_seq->j);CHKERRQ(ierr); 324 for (i=0; i<cm; i++) { 325 ncols = c_seq->i[i+1] - c_seq->i[i]; 326 row = p->garray[i]; 327 ierr = MatSetValues(C,1,&row,ncols,cols,vals,ADD_VALUES);CHKERRQ(ierr); 328 cols += ncols; vals += ncols; 329 } 330 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 331 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 332 333 ptap->reuse = MAT_REUSE_MATRIX; 334 335 ierr = ISGlobalToLocalMappingApply(ptap->ltog,IS_GTOLM_DROP,c_seq->i[C_oth->rmap->n],c_seq->j,&nout,c_seq->j);CHKERRQ(ierr); 336 if (c_seq->i[C_oth->rmap->n] != nout) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Incorrect mapping %D != %D\n",c_seq->i[C_loc->rmap->n],nout); 337 338 /* supporting struct ptap consumes almost same amount of memory as C=PtAP, release it if C will not be updated by A and P */ 339 if (ptap->freestruct) { 340 ierr = MatFreeIntermediateDataStructures(C);CHKERRQ(ierr); 341 } 342 PetscFunctionReturn(0); 343 } 344 345 PetscErrorCode MatPtAPSymbolic_MPIAIJ_MPIAIJ_scalable(Mat A,Mat P,PetscReal fill,Mat *C) 346 { 347 PetscErrorCode ierr; 348 Mat_APMPI *ptap; 349 Mat_MPIAIJ *a=(Mat_MPIAIJ*)A->data,*p=(Mat_MPIAIJ*)P->data,*c; 350 MPI_Comm comm; 351 PetscMPIInt size,rank; 352 Mat Cmpi,P_loc,P_oth; 353 PetscFreeSpaceList free_space=NULL,current_space=NULL; 354 PetscInt am=A->rmap->n,pm=P->rmap->n,pN=P->cmap->N,pn=P->cmap->n; 355 PetscInt *lnk,i,k,pnz,row,nsend; 356 PetscMPIInt tagi,tagj,*len_si,*len_s,*len_ri,icompleted=0,nrecv; 357 PetscInt **buf_rj,**buf_ri,**buf_ri_k; 358 PetscInt len,proc,*dnz,*onz,*owners,nzi,nspacedouble; 359 PetscInt nrows,*buf_s,*buf_si,*buf_si_i,**nextrow,**nextci; 360 MPI_Request *swaits,*rwaits; 361 MPI_Status *sstatus,rstatus; 362 PetscLayout rowmap; 363 PetscInt *owners_co,*coi,*coj; /* i and j array of (p->B)^T*A*P - used in the communication */ 364 PetscMPIInt *len_r,*id_r; /* array of length of comm->size, store send/recv matrix values */ 365 PetscInt *api,*apj,*Jptr,apnz,*prmap=p->garray,con,j,Crmax,*aj,*ai,*pi,nout; 366 Mat_SeqAIJ *p_loc,*p_oth=NULL,*ad=(Mat_SeqAIJ*)(a->A)->data,*ao=NULL,*c_loc,*c_oth; 367 PetscScalar *apv; 368 PetscTable ta; 369 MatType mtype; 370 const char *prefix; 371 #if defined(PETSC_USE_INFO) 372 PetscReal apfill; 373 #endif 374 375 PetscFunctionBegin; 376 ierr = PetscObjectGetComm((PetscObject)A,&comm);CHKERRQ(ierr); 377 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 378 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 379 380 if (size > 1) ao = (Mat_SeqAIJ*)(a->B)->data; 381 382 /* create symbolic parallel matrix Cmpi */ 383 ierr = MatCreate(comm,&Cmpi);CHKERRQ(ierr); 384 ierr = MatGetType(A,&mtype);CHKERRQ(ierr); 385 ierr = MatSetType(Cmpi,mtype);CHKERRQ(ierr); 386 387 /* create struct Mat_APMPI and attached it to C later */ 388 ierr = PetscNew(&ptap);CHKERRQ(ierr); 389 ptap->reuse = MAT_INITIAL_MATRIX; 390 ptap->algType = 0; 391 392 /* get P_oth by taking rows of P (= non-zero cols of local A) from other processors */ 393 ierr = MatGetBrowsOfAoCols_MPIAIJ(A,P,MAT_INITIAL_MATRIX,&ptap->startsj_s,&ptap->startsj_r,&ptap->bufa,&P_oth);CHKERRQ(ierr); 394 /* get P_loc by taking all local rows of P */ 395 ierr = MatMPIAIJGetLocalMat(P,MAT_INITIAL_MATRIX,&P_loc);CHKERRQ(ierr); 396 397 ptap->P_loc = P_loc; 398 ptap->P_oth = P_oth; 399 400 /* (0) compute Rd = Pd^T, Ro = Po^T */ 401 /* --------------------------------- */ 402 ierr = MatTranspose(p->A,MAT_INITIAL_MATRIX,&ptap->Rd);CHKERRQ(ierr); 403 ierr = MatTranspose(p->B,MAT_INITIAL_MATRIX,&ptap->Ro);CHKERRQ(ierr); 404 405 /* (1) compute symbolic AP = A_loc*P = Ad*P_loc + Ao*P_oth (api,apj) */ 406 /* ----------------------------------------------------------------- */ 407 p_loc = (Mat_SeqAIJ*)P_loc->data; 408 if (P_oth) p_oth = (Mat_SeqAIJ*)P_oth->data; 409 410 /* create and initialize a linked list */ 411 ierr = PetscTableCreate(pn,pN,&ta);CHKERRQ(ierr); /* for compute AP_loc and Cmpi */ 412 MatRowMergeMax_SeqAIJ(p_loc,P_loc->rmap->N,ta); 413 MatRowMergeMax_SeqAIJ(p_oth,P_oth->rmap->N,ta); 414 ierr = PetscTableGetCount(ta,&Crmax);CHKERRQ(ierr); /* Crmax = nnz(sum of Prows) */ 415 416 ierr = PetscLLCondensedCreate_Scalable(Crmax,&lnk);CHKERRQ(ierr); 417 418 /* Initial FreeSpace size is fill*(nnz(A) + nnz(P)) */ 419 if (ao) { 420 ierr = PetscFreeSpaceGet(PetscRealIntMultTruncate(fill,PetscIntSumTruncate(ad->i[am],PetscIntSumTruncate(ao->i[am],p_loc->i[pm]))),&free_space);CHKERRQ(ierr); 421 } else { 422 ierr = PetscFreeSpaceGet(PetscRealIntMultTruncate(fill,PetscIntSumTruncate(ad->i[am],p_loc->i[pm])),&free_space);CHKERRQ(ierr); 423 } 424 current_space = free_space; 425 nspacedouble = 0; 426 427 ierr = PetscMalloc1(am+1,&api);CHKERRQ(ierr); 428 api[0] = 0; 429 for (i=0; i<am; i++) { 430 /* diagonal portion: Ad[i,:]*P */ 431 ai = ad->i; pi = p_loc->i; 432 nzi = ai[i+1] - ai[i]; 433 aj = ad->j + ai[i]; 434 for (j=0; j<nzi; j++) { 435 row = aj[j]; 436 pnz = pi[row+1] - pi[row]; 437 Jptr = p_loc->j + pi[row]; 438 /* add non-zero cols of P into the sorted linked list lnk */ 439 ierr = PetscLLCondensedAddSorted_Scalable(pnz,Jptr,lnk);CHKERRQ(ierr); 440 } 441 /* off-diagonal portion: Ao[i,:]*P */ 442 if (ao) { 443 ai = ao->i; pi = p_oth->i; 444 nzi = ai[i+1] - ai[i]; 445 aj = ao->j + ai[i]; 446 for (j=0; j<nzi; j++) { 447 row = aj[j]; 448 pnz = pi[row+1] - pi[row]; 449 Jptr = p_oth->j + pi[row]; 450 ierr = PetscLLCondensedAddSorted_Scalable(pnz,Jptr,lnk);CHKERRQ(ierr); 451 } 452 } 453 apnz = lnk[0]; 454 api[i+1] = api[i] + apnz; 455 456 /* if free space is not available, double the total space in the list */ 457 if (current_space->local_remaining<apnz) { 458 ierr = PetscFreeSpaceGet(PetscIntSumTruncate(apnz,current_space->total_array_size),¤t_space);CHKERRQ(ierr); 459 nspacedouble++; 460 } 461 462 /* Copy data into free space, then initialize lnk */ 463 ierr = PetscLLCondensedClean_Scalable(apnz,current_space->array,lnk);CHKERRQ(ierr); 464 465 current_space->array += apnz; 466 current_space->local_used += apnz; 467 current_space->local_remaining -= apnz; 468 } 469 /* Allocate space for apj and apv, initialize apj, and */ 470 /* destroy list of free space and other temporary array(s) */ 471 ierr = PetscCalloc2(api[am],&apj,api[am],&apv);CHKERRQ(ierr); 472 ierr = PetscFreeSpaceContiguous(&free_space,apj);CHKERRQ(ierr); 473 ierr = PetscLLCondensedDestroy_Scalable(lnk);CHKERRQ(ierr); 474 475 /* Create AP_loc for reuse */ 476 ierr = MatCreateSeqAIJWithArrays(PETSC_COMM_SELF,am,pN,api,apj,apv,&ptap->AP_loc);CHKERRQ(ierr); 477 ierr = MatSeqAIJCompactOutExtraColumns_SeqAIJ(ptap->AP_loc, &ptap->ltog);CHKERRQ(ierr); 478 479 #if defined(PETSC_USE_INFO) 480 if (ao) { 481 apfill = (PetscReal)api[am]/(ad->i[am]+ao->i[am]+p_loc->i[pm]+1); 482 } else { 483 apfill = (PetscReal)api[am]/(ad->i[am]+p_loc->i[pm]+1); 484 } 485 ptap->AP_loc->info.mallocs = nspacedouble; 486 ptap->AP_loc->info.fill_ratio_given = fill; 487 ptap->AP_loc->info.fill_ratio_needed = apfill; 488 489 if (api[am]) { 490 ierr = PetscInfo3(ptap->AP_loc,"Scalable algorithm, AP_loc reallocs %D; Fill ratio: given %g needed %g.\n",nspacedouble,(double)fill,(double)apfill);CHKERRQ(ierr); 491 ierr = PetscInfo1(ptap->AP_loc,"Use MatPtAP(A,B,MatReuse,%g,&C) for best AP_loc performance.;\n",(double)apfill);CHKERRQ(ierr); 492 } else { 493 ierr = PetscInfo(ptap->AP_loc,"Scalable algorithm, AP_loc is empty \n");CHKERRQ(ierr); 494 } 495 #endif 496 497 /* (2-1) compute symbolic Co = Ro*AP_loc */ 498 /* ------------------------------------ */ 499 ierr = MatGetOptionsPrefix(A,&prefix);CHKERRQ(ierr); 500 ierr = MatSetOptionsPrefix(ptap->Ro,prefix);CHKERRQ(ierr); 501 ierr = MatAppendOptionsPrefix(ptap->Ro,"inner_offdiag_");CHKERRQ(ierr); 502 ierr = MatMatMultSymbolic_SeqAIJ_SeqAIJ(ptap->Ro,ptap->AP_loc,fill,&ptap->C_oth);CHKERRQ(ierr); 503 504 /* (3) send coj of C_oth to other processors */ 505 /* ------------------------------------------ */ 506 /* determine row ownership */ 507 ierr = PetscLayoutCreate(comm,&rowmap);CHKERRQ(ierr); 508 rowmap->n = pn; 509 rowmap->bs = 1; 510 ierr = PetscLayoutSetUp(rowmap);CHKERRQ(ierr); 511 owners = rowmap->range; 512 513 /* determine the number of messages to send, their lengths */ 514 ierr = PetscMalloc4(size,&len_s,size,&len_si,size,&sstatus,size+2,&owners_co);CHKERRQ(ierr); 515 ierr = PetscMemzero(len_s,size*sizeof(PetscMPIInt));CHKERRQ(ierr); 516 ierr = PetscMemzero(len_si,size*sizeof(PetscMPIInt));CHKERRQ(ierr); 517 518 c_oth = (Mat_SeqAIJ*)ptap->C_oth->data; 519 coi = c_oth->i; coj = c_oth->j; 520 con = ptap->C_oth->rmap->n; 521 proc = 0; 522 ierr = ISLocalToGlobalMappingApply(ptap->ltog,coi[con],coj,coj);CHKERRQ(ierr); 523 for (i=0; i<con; i++) { 524 while (prmap[i] >= owners[proc+1]) proc++; 525 len_si[proc]++; /* num of rows in Co(=Pt*AP) to be sent to [proc] */ 526 len_s[proc] += coi[i+1] - coi[i]; /* num of nonzeros in Co to be sent to [proc] */ 527 } 528 529 len = 0; /* max length of buf_si[], see (4) */ 530 owners_co[0] = 0; 531 nsend = 0; 532 for (proc=0; proc<size; proc++) { 533 owners_co[proc+1] = owners_co[proc] + len_si[proc]; 534 if (len_s[proc]) { 535 nsend++; 536 len_si[proc] = 2*(len_si[proc] + 1); /* length of buf_si to be sent to [proc] */ 537 len += len_si[proc]; 538 } 539 } 540 541 /* determine the number and length of messages to receive for coi and coj */ 542 ierr = PetscGatherNumberOfMessages(comm,NULL,len_s,&nrecv);CHKERRQ(ierr); 543 ierr = PetscGatherMessageLengths2(comm,nsend,nrecv,len_s,len_si,&id_r,&len_r,&len_ri);CHKERRQ(ierr); 544 545 /* post the Irecv and Isend of coj */ 546 ierr = PetscCommGetNewTag(comm,&tagj);CHKERRQ(ierr); 547 ierr = PetscPostIrecvInt(comm,tagj,nrecv,id_r,len_r,&buf_rj,&rwaits);CHKERRQ(ierr); 548 ierr = PetscMalloc1(nsend+1,&swaits);CHKERRQ(ierr); 549 for (proc=0, k=0; proc<size; proc++) { 550 if (!len_s[proc]) continue; 551 i = owners_co[proc]; 552 ierr = MPI_Isend(coj+coi[i],len_s[proc],MPIU_INT,proc,tagj,comm,swaits+k);CHKERRQ(ierr); 553 k++; 554 } 555 556 /* (2-2) compute symbolic C_loc = Rd*AP_loc */ 557 /* ---------------------------------------- */ 558 ierr = MatSetOptionsPrefix(ptap->Rd,prefix);CHKERRQ(ierr); 559 ierr = MatAppendOptionsPrefix(ptap->Rd,"inner_diag_");CHKERRQ(ierr); 560 ierr = MatMatMultSymbolic_SeqAIJ_SeqAIJ(ptap->Rd,ptap->AP_loc,fill,&ptap->C_loc);CHKERRQ(ierr); 561 c_loc = (Mat_SeqAIJ*)ptap->C_loc->data; 562 ierr = ISLocalToGlobalMappingApply(ptap->ltog,c_loc->i[ptap->C_loc->rmap->n],c_loc->j,c_loc->j);CHKERRQ(ierr); 563 564 /* receives coj are complete */ 565 for (i=0; i<nrecv; i++) { 566 ierr = MPI_Waitany(nrecv,rwaits,&icompleted,&rstatus);CHKERRQ(ierr); 567 } 568 ierr = PetscFree(rwaits);CHKERRQ(ierr); 569 if (nsend) {ierr = MPI_Waitall(nsend,swaits,sstatus);CHKERRQ(ierr);} 570 571 /* add received column indices into ta to update Crmax */ 572 for (k=0; k<nrecv; k++) {/* k-th received message */ 573 Jptr = buf_rj[k]; 574 for (j=0; j<len_r[k]; j++) { 575 ierr = PetscTableAdd(ta,*(Jptr+j)+1,1,INSERT_VALUES);CHKERRQ(ierr); 576 } 577 } 578 ierr = PetscTableGetCount(ta,&Crmax);CHKERRQ(ierr); 579 ierr = PetscTableDestroy(&ta);CHKERRQ(ierr); 580 581 /* (4) send and recv coi */ 582 /*-----------------------*/ 583 ierr = PetscCommGetNewTag(comm,&tagi);CHKERRQ(ierr); 584 ierr = PetscPostIrecvInt(comm,tagi,nrecv,id_r,len_ri,&buf_ri,&rwaits);CHKERRQ(ierr); 585 ierr = PetscMalloc1(len+1,&buf_s);CHKERRQ(ierr); 586 buf_si = buf_s; /* points to the beginning of k-th msg to be sent */ 587 for (proc=0,k=0; proc<size; proc++) { 588 if (!len_s[proc]) continue; 589 /* form outgoing message for i-structure: 590 buf_si[0]: nrows to be sent 591 [1:nrows]: row index (global) 592 [nrows+1:2*nrows+1]: i-structure index 593 */ 594 /*-------------------------------------------*/ 595 nrows = len_si[proc]/2 - 1; /* num of rows in Co to be sent to [proc] */ 596 buf_si_i = buf_si + nrows+1; 597 buf_si[0] = nrows; 598 buf_si_i[0] = 0; 599 nrows = 0; 600 for (i=owners_co[proc]; i<owners_co[proc+1]; i++) { 601 nzi = coi[i+1] - coi[i]; 602 buf_si_i[nrows+1] = buf_si_i[nrows] + nzi; /* i-structure */ 603 buf_si[nrows+1] = prmap[i] -owners[proc]; /* local row index */ 604 nrows++; 605 } 606 ierr = MPI_Isend(buf_si,len_si[proc],MPIU_INT,proc,tagi,comm,swaits+k);CHKERRQ(ierr); 607 k++; 608 buf_si += len_si[proc]; 609 } 610 for (i=0; i<nrecv; i++) { 611 ierr = MPI_Waitany(nrecv,rwaits,&icompleted,&rstatus);CHKERRQ(ierr); 612 } 613 ierr = PetscFree(rwaits);CHKERRQ(ierr); 614 if (nsend) {ierr = MPI_Waitall(nsend,swaits,sstatus);CHKERRQ(ierr);} 615 616 ierr = PetscFree4(len_s,len_si,sstatus,owners_co);CHKERRQ(ierr); 617 ierr = PetscFree(len_ri);CHKERRQ(ierr); 618 ierr = PetscFree(swaits);CHKERRQ(ierr); 619 ierr = PetscFree(buf_s);CHKERRQ(ierr); 620 621 /* (5) compute the local portion of Cmpi */ 622 /* ------------------------------------------ */ 623 /* set initial free space to be Crmax, sufficient for holding nozeros in each row of Cmpi */ 624 ierr = PetscFreeSpaceGet(Crmax,&free_space);CHKERRQ(ierr); 625 current_space = free_space; 626 627 ierr = PetscMalloc3(nrecv,&buf_ri_k,nrecv,&nextrow,nrecv,&nextci);CHKERRQ(ierr); 628 for (k=0; k<nrecv; k++) { 629 buf_ri_k[k] = buf_ri[k]; /* beginning of k-th recved i-structure */ 630 nrows = *buf_ri_k[k]; 631 nextrow[k] = buf_ri_k[k] + 1; /* next row number of k-th recved i-structure */ 632 nextci[k] = buf_ri_k[k] + (nrows + 1); /* poins to the next i-structure of k-th recved i-structure */ 633 } 634 635 ierr = MatPreallocateInitialize(comm,pn,pn,dnz,onz);CHKERRQ(ierr); 636 ierr = PetscLLCondensedCreate_Scalable(Crmax,&lnk);CHKERRQ(ierr); 637 for (i=0; i<pn; i++) { 638 /* add C_loc into Cmpi */ 639 nzi = c_loc->i[i+1] - c_loc->i[i]; 640 Jptr = c_loc->j + c_loc->i[i]; 641 ierr = PetscLLCondensedAddSorted_Scalable(nzi,Jptr,lnk);CHKERRQ(ierr); 642 643 /* add received col data into lnk */ 644 for (k=0; k<nrecv; k++) { /* k-th received message */ 645 if (i == *nextrow[k]) { /* i-th row */ 646 nzi = *(nextci[k]+1) - *nextci[k]; 647 Jptr = buf_rj[k] + *nextci[k]; 648 ierr = PetscLLCondensedAddSorted_Scalable(nzi,Jptr,lnk);CHKERRQ(ierr); 649 nextrow[k]++; nextci[k]++; 650 } 651 } 652 nzi = lnk[0]; 653 654 /* copy data into free space, then initialize lnk */ 655 ierr = PetscLLCondensedClean_Scalable(nzi,current_space->array,lnk);CHKERRQ(ierr); 656 ierr = MatPreallocateSet(i+owners[rank],nzi,current_space->array,dnz,onz);CHKERRQ(ierr); 657 } 658 ierr = PetscFree3(buf_ri_k,nextrow,nextci);CHKERRQ(ierr); 659 ierr = PetscLLCondensedDestroy_Scalable(lnk);CHKERRQ(ierr); 660 ierr = PetscFreeSpaceDestroy(free_space);CHKERRQ(ierr); 661 662 /* local sizes and preallocation */ 663 ierr = MatSetSizes(Cmpi,pn,pn,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 664 ierr = MatSetBlockSizes(Cmpi,PetscAbs(P->cmap->bs),PetscAbs(P->cmap->bs));CHKERRQ(ierr); 665 ierr = MatMPIAIJSetPreallocation(Cmpi,0,dnz,0,onz);CHKERRQ(ierr); 666 ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr); 667 668 /* members in merge */ 669 ierr = PetscFree(id_r);CHKERRQ(ierr); 670 ierr = PetscFree(len_r);CHKERRQ(ierr); 671 ierr = PetscFree(buf_ri[0]);CHKERRQ(ierr); 672 ierr = PetscFree(buf_ri);CHKERRQ(ierr); 673 ierr = PetscFree(buf_rj[0]);CHKERRQ(ierr); 674 ierr = PetscFree(buf_rj);CHKERRQ(ierr); 675 ierr = PetscLayoutDestroy(&rowmap);CHKERRQ(ierr); 676 677 /* attach the supporting struct to Cmpi for reuse */ 678 c = (Mat_MPIAIJ*)Cmpi->data; 679 c->ap = ptap; 680 ptap->duplicate = Cmpi->ops->duplicate; 681 ptap->destroy = Cmpi->ops->destroy; 682 ptap->view = Cmpi->ops->view; 683 684 /* Cmpi is not ready for use - assembly will be done by MatPtAPNumeric() */ 685 Cmpi->assembled = PETSC_FALSE; 686 Cmpi->ops->ptapnumeric = MatPtAPNumeric_MPIAIJ_MPIAIJ_scalable; 687 Cmpi->ops->destroy = MatDestroy_MPIAIJ_PtAP; 688 Cmpi->ops->view = MatView_MPIAIJ_PtAP; 689 Cmpi->ops->freeintermediatedatastructures = MatFreeIntermediateDataStructures_MPIAIJ_AP; 690 *C = Cmpi; 691 692 nout = 0; 693 ierr = ISGlobalToLocalMappingApply(ptap->ltog,IS_GTOLM_DROP,c_oth->i[ptap->C_oth->rmap->n],c_oth->j,&nout,c_oth->j);CHKERRQ(ierr); 694 if (c_oth->i[ptap->C_oth->rmap->n] != nout) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Incorrect mapping %D != %D\n",c_oth->i[ptap->C_oth->rmap->n],nout); 695 ierr = ISGlobalToLocalMappingApply(ptap->ltog,IS_GTOLM_DROP,c_loc->i[ptap->C_loc->rmap->n],c_loc->j,&nout,c_loc->j);CHKERRQ(ierr); 696 if (c_loc->i[ptap->C_loc->rmap->n] != nout) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Incorrect mapping %D != %D\n",c_loc->i[ptap->C_loc->rmap->n],nout); 697 698 PetscFunctionReturn(0); 699 } 700 701 PETSC_STATIC_INLINE PetscErrorCode MatPtAPSymbolicComputeOneRowOfAP_private(Mat A,Mat P,Mat P_oth,PetscInt i,PetscHSetI dht,PetscHSetI oht) 702 { 703 Mat_MPIAIJ *a=(Mat_MPIAIJ*)A->data,*p=(Mat_MPIAIJ*)P->data; 704 Mat_SeqAIJ *ad=(Mat_SeqAIJ*)(a->A)->data,*ao=(Mat_SeqAIJ*)(a->B)->data,*p_oth=(Mat_SeqAIJ*)P_oth->data,*pd=(Mat_SeqAIJ*)p->A->data,*po=(Mat_SeqAIJ*)p->B->data; 705 PetscInt *ai,nzi,j,*aj,row,col,*pi,*pj,pnz,nzpi,*p_othcols,k; 706 PetscInt pcstart,pcend,column; 707 PetscErrorCode ierr; 708 709 PetscFunctionBegin; 710 711 pcstart = P->cmap->rstart; 712 pcend = P->cmap->rend; 713 /* diagonal portion: Ad[i,:]*P */ 714 ai = ad->i; 715 nzi = ai[i+1] - ai[i]; 716 aj = ad->j + ai[i]; 717 for (j=0; j<nzi; j++) { 718 row = aj[j]; 719 nzpi = pd->i[row+1] - pd->i[row]; 720 pj = pd->j + pd->i[row]; 721 for (k=0; k<nzpi; k++) { 722 ierr = PetscHSetIAdd(dht,pj[k]+pcstart);CHKERRQ(ierr); 723 } 724 } 725 for (j=0; j<nzi; j++) { 726 row = aj[j]; 727 nzpi = po->i[row+1] - po->i[row]; 728 pj = po->j + po->i[row]; 729 for (k=0; k<nzpi; k++) { 730 ierr = PetscHSetIAdd(oht,p->garray[pj[k]]);CHKERRQ(ierr); 731 } 732 } 733 734 /* off diagonal part: Ao[i, :]*P_oth */ 735 if (ao) { 736 ai = ao->i; 737 pi = p_oth->i; 738 nzi = ai[i+1] - ai[i]; 739 aj = ao->j + ai[i]; 740 for (j=0; j<nzi; j++) { 741 row = aj[j]; 742 pnz = pi[row+1] - pi[row]; 743 p_othcols = p_oth->j + pi[row]; 744 for (col=0; col<pnz; col++) { 745 column = p_othcols[col]; 746 if (column>=pcstart && column<pcend) { 747 ierr = PetscHSetIAdd(dht,column);CHKERRQ(ierr); 748 } else { 749 ierr = PetscHSetIAdd(oht,column);CHKERRQ(ierr); 750 } 751 } 752 } 753 } /* end if (ao) */ 754 PetscFunctionReturn(0); 755 } 756 757 PETSC_STATIC_INLINE PetscErrorCode MatPtAPNumericComputeOneRowOfAP_private(Mat A,Mat P,Mat P_oth,PetscInt i,PetscHMapIV hmap) 758 { 759 Mat_MPIAIJ *a=(Mat_MPIAIJ*)A->data,*p=(Mat_MPIAIJ*)P->data; 760 Mat_SeqAIJ *ad=(Mat_SeqAIJ*)(a->A)->data,*ao=(Mat_SeqAIJ*)(a->B)->data,*p_oth=(Mat_SeqAIJ*)P_oth->data,*pd=(Mat_SeqAIJ*)p->A->data,*po=(Mat_SeqAIJ*)p->B->data; 761 PetscInt *ai,nzi,j,*aj,row,col,*pi,pnz,*p_othcols,pcstart,*pj,k,nzpi; 762 PetscScalar ra,*aa,*pa; 763 PetscErrorCode ierr; 764 765 PetscFunctionBegin; 766 pcstart = P->cmap->rstart; 767 /* diagonal portion: Ad[i,:]*P */ 768 ai = ad->i; 769 nzi = ai[i+1] - ai[i]; 770 aj = ad->j + ai[i]; 771 aa = ad->a + ai[i]; 772 773 for (j=0; j<nzi; j++) { 774 ra = aa[j]; 775 row = aj[j]; 776 nzpi = pd->i[row+1] - pd->i[row]; 777 pj = pd->j + pd->i[row]; 778 pa = pd->a + pd->i[row]; 779 for (k=0; k<nzpi; k++) { 780 ierr = PetscHMapIVAdd(hmap,pj[k]+pcstart,ra*pa[k]);CHKERRQ(ierr); 781 } 782 } 783 for (j=0; j<nzi; j++) { 784 ra = aa[j]; 785 row = aj[j]; 786 nzpi = po->i[row+1] - po->i[row]; 787 pj = po->j + po->i[row]; 788 pa = po->a + po->i[row]; 789 for (k=0; k<nzpi; k++) { 790 ierr = PetscHMapIVAdd(hmap,p->garray[pj[k]],ra*pa[k]);CHKERRQ(ierr); 791 } 792 } 793 794 795 /* off diagonal part: Ao[i, :]*P_oth */ 796 if (ao) { 797 ai = ao->i; 798 pi = p_oth->i; 799 nzi = ai[i+1] - ai[i]; 800 aj = ao->j + ai[i]; 801 aa = ao->a + ai[i]; 802 for (j=0; j<nzi; j++) { 803 row = aj[j]; 804 ra = aa[j]; 805 pnz = pi[row+1] - pi[row]; 806 p_othcols = p_oth->j + pi[row]; 807 pa = p_oth->a + pi[row]; 808 for (col=0; col<pnz; col++) { 809 ierr = PetscHMapIVAdd(hmap,p_othcols[col],ra*pa[col]);CHKERRQ(ierr); 810 } 811 } 812 } /* end if (ao) */ 813 PetscFunctionReturn(0); 814 } 815 816 PetscErrorCode MatPtAPNumeric_MPIAIJ_MPIAIJ_allatonce(Mat A,Mat P,Mat C) 817 { 818 PetscErrorCode ierr; 819 Mat_MPIAIJ *p=(Mat_MPIAIJ*)P->data,*c=(Mat_MPIAIJ*)C->data; 820 Mat_SeqAIJ *cd,*co,*po,*pd; 821 Mat_APMPI *ptap = c->ap; 822 PetscHMapIV hmap; 823 PetscInt i,j,jj,kk,nzi,*c_rmtj,voff,*c_othj,pn,pon,pcstart,pcend,ccstart,ccend,row,am,*poj,*pdj,*apindices,cmaxr,*c_rmtc,*c_rmtjj,*dcc,*occ,loc; 824 PetscScalar *c_rmta,*c_otha,*poa,*pda,*apvalues,*apvaluestmp,*c_rmtaa; 825 MPI_Comm comm; 826 827 PetscFunctionBegin; 828 ierr = PetscObjectGetComm((PetscObject)A,&comm);CHKERRQ(ierr); 829 if (!ptap) SETERRQ(comm,PETSC_ERR_ARG_WRONGSTATE,"PtAP cannot be reused. Do not call MatFreeIntermediateDataStructures() or use '-mat_freeintermediatedatastructures'"); 830 831 ierr = MatZeroEntries(C);CHKERRQ(ierr); 832 833 /* Get P_oth = ptap->P_oth and P_loc = ptap->P_loc */ 834 /*-----------------------------------------------------*/ 835 if (ptap->reuse == MAT_REUSE_MATRIX) { 836 /* P_oth and P_loc are obtained in MatPtASymbolic() when reuse == MAT_INITIAL_MATRIX */ 837 ierr = MatGetBrowsOfAoCols_MPIAIJ(A,P,MAT_REUSE_MATRIX,&ptap->startsj_s,&ptap->startsj_r,&ptap->bufa,&ptap->P_oth);CHKERRQ(ierr); 838 } 839 840 po = (Mat_SeqAIJ*) p->B->data; 841 pd = (Mat_SeqAIJ*) p->A->data; 842 843 ierr = MatGetLocalSize(p->B,NULL,&pon);CHKERRQ(ierr); 844 845 ierr = PetscCalloc2(ptap->c_rmti[pon],&c_rmtj,ptap->c_rmti[pon],&c_rmta);CHKERRQ(ierr); 846 ierr = MatGetLocalSize(A,&am,NULL);CHKERRQ(ierr); 847 cmaxr = 0; 848 for (i=0; i<pon; i++) { 849 cmaxr = PetscMax(cmaxr,ptap->c_rmti[i+1]-ptap->c_rmti[i]); 850 } 851 ierr = PetscCalloc4(cmaxr,&apindices,cmaxr,&apvalues,cmaxr,&apvaluestmp,pon,&c_rmtc);CHKERRQ(ierr); 852 ierr = PetscHMapIVCreate(&hmap);CHKERRQ(ierr); 853 ierr = PetscHMapIVResize(hmap,cmaxr);CHKERRQ(ierr); 854 for (i=0; i<am && pon; i++) { 855 ierr = PetscHMapIVClear(hmap);CHKERRQ(ierr); 856 nzi = po->i[i+1] - po->i[i]; 857 if (!nzi) continue; 858 ierr = MatPtAPNumericComputeOneRowOfAP_private(A,P,ptap->P_oth,i,hmap);CHKERRQ(ierr); 859 voff = 0; 860 ierr = PetscHMapIVGetPairs(hmap,&voff,apindices,apvalues);CHKERRQ(ierr); 861 if (!voff) continue; 862 /*ierr = PetscMemzero(c_rmtc,sizeof(PetscInt)*pon);CHKERRQ(ierr);*/ 863 864 /* Form C(ii, :) */ 865 poj = po->j + po->i[i]; 866 poa = po->a + po->i[i]; 867 for (j=0; j<nzi; j++) { 868 c_rmtjj = c_rmtj + ptap->c_rmti[poj[j]]; 869 c_rmtaa = c_rmta + ptap->c_rmti[poj[j]]; 870 for (jj=0; jj<voff; jj++) { 871 apvaluestmp[jj] = apvalues[jj]*poa[j]; 872 /*If the row is empty */ 873 if (!c_rmtc[poj[j]]) { 874 c_rmtjj[jj] = apindices[jj]; 875 c_rmtaa[jj] = apvaluestmp[jj]; 876 c_rmtc[poj[j]]++; 877 } else { 878 ierr = PetscFindInt(apindices[jj],c_rmtc[poj[j]],c_rmtjj,&loc);CHKERRQ(ierr); 879 if (loc>=0){ /* hit */ 880 c_rmtaa[loc] += apvaluestmp[jj]; 881 } else { /* new element */ 882 loc = -(loc+1); 883 /* Move data backward */ 884 for (kk=c_rmtc[poj[j]]; kk>loc; kk--) { 885 c_rmtjj[kk] = c_rmtjj[kk-1]; 886 c_rmtaa[kk] = c_rmtaa[kk-1]; 887 }/* End kk */ 888 c_rmtjj[loc] = apindices[jj]; 889 c_rmtaa[loc] = apvaluestmp[jj]; 890 c_rmtc[poj[j]]++; 891 } 892 } 893 } /* End jj */ 894 } /* End j */ 895 } /* End i */ 896 897 ierr = PetscFree4(apindices,apvalues,apvaluestmp,c_rmtc);CHKERRQ(ierr); 898 899 ierr = MatGetLocalSize(P,NULL,&pn);CHKERRQ(ierr); 900 ierr = PetscCalloc2(ptap->c_othi[pn],&c_othj,ptap->c_othi[pn],&c_otha);CHKERRQ(ierr); 901 902 ierr = PetscSFReduceBegin(ptap->sf,MPIU_INT,c_rmtj,c_othj,MPIU_REPLACE);CHKERRQ(ierr); 903 ierr = PetscSFReduceBegin(ptap->sf,MPIU_SCALAR,c_rmta,c_otha,MPIU_REPLACE);CHKERRQ(ierr); 904 ierr = MatGetOwnershipRangeColumn(P,&pcstart,&pcend);CHKERRQ(ierr); 905 cd = (Mat_SeqAIJ*)(c->A)->data; 906 co = (Mat_SeqAIJ*)(c->B)->data; 907 908 cmaxr = 0; 909 for (i=0; i<pn; i++) { 910 cmaxr = PetscMax(cmaxr,(cd->i[i+1]-cd->i[i])+(co->i[i+1]-co->i[i])); 911 } 912 ierr = PetscCalloc5(cmaxr,&apindices,cmaxr,&apvalues,cmaxr,&apvaluestmp,pn,&dcc,pn,&occ);CHKERRQ(ierr); 913 ierr = PetscHMapIVCreate(&hmap);CHKERRQ(ierr); 914 ierr = PetscHMapIVResize(hmap,cmaxr);CHKERRQ(ierr); 915 for (i=0; i<am && pn; i++) { 916 ierr = PetscHMapIVClear(hmap);CHKERRQ(ierr); 917 nzi = pd->i[i+1] - pd->i[i]; 918 if (!nzi) continue; 919 ierr = MatPtAPNumericComputeOneRowOfAP_private(A,P,ptap->P_oth,i,hmap);CHKERRQ(ierr); 920 voff = 0; 921 ierr = PetscHMapIVGetPairs(hmap,&voff,apindices,apvalues);CHKERRQ(ierr); 922 if (!voff) continue; 923 /* Form C(ii, :) */ 924 pdj = pd->j + pd->i[i]; 925 pda = pd->a + pd->i[i]; 926 for (j=0; j<nzi; j++) { 927 row = pcstart + pdj[j]; 928 for (jj=0; jj<voff; jj++) { 929 apvaluestmp[jj] = apvalues[jj]*pda[j]; 930 } 931 ierr = MatSetValues(C,1,&row,voff,apindices,apvaluestmp,ADD_VALUES);CHKERRQ(ierr); 932 } 933 } 934 935 ierr = MatGetOwnershipRangeColumn(C,&ccstart,&ccend);CHKERRQ(ierr); 936 ierr = PetscFree5(apindices,apvalues,apvaluestmp,dcc,occ);CHKERRQ(ierr); 937 ierr = PetscHMapIVDestroy(&hmap);CHKERRQ(ierr); 938 ierr = PetscSFReduceEnd(ptap->sf,MPIU_INT,c_rmtj,c_othj,MPI_SUM);CHKERRQ(ierr); 939 ierr = PetscSFReduceEnd(ptap->sf,MPIU_SCALAR,c_rmta,c_otha,MPIU_REPLACE);CHKERRQ(ierr); 940 ierr = PetscFree2(c_rmtj,c_rmta);CHKERRQ(ierr); 941 942 /* Add contributions from remote */ 943 for (i = 0; i < pn; i++) { 944 row = i + pcstart; 945 ierr = MatSetValues(C,1,&row,ptap->c_othi[i+1]-ptap->c_othi[i],c_othj+ptap->c_othi[i],c_otha+ptap->c_othi[i],ADD_VALUES);CHKERRQ(ierr); 946 } 947 ierr = PetscFree2(c_othj,c_otha);CHKERRQ(ierr); 948 949 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 950 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 951 952 ptap->reuse = MAT_REUSE_MATRIX; 953 954 /* supporting struct ptap consumes almost same amount of memory as C=PtAP, release it if C will not be updated by A and P */ 955 if (ptap->freestruct) { 956 ierr = MatFreeIntermediateDataStructures(C);CHKERRQ(ierr); 957 } 958 PetscFunctionReturn(0); 959 } 960 961 PetscErrorCode MatPtAPNumeric_MPIAIJ_MPIAIJ_allatonce_merged(Mat A,Mat P,Mat C) 962 { 963 PetscErrorCode ierr; 964 Mat_MPIAIJ *p=(Mat_MPIAIJ*)P->data,*c=(Mat_MPIAIJ*)C->data; 965 Mat_SeqAIJ *cd,*co,*po,*pd; 966 Mat_APMPI *ptap = c->ap; 967 PetscHMapIV hmap; 968 PetscInt i,j,jj,kk,nzi,dnzi,*c_rmtj,voff,*c_othj,pn,pon,pcstart,pcend,row,am,*poj,*pdj,*apindices,cmaxr,*c_rmtc,*c_rmtjj,loc; 969 PetscScalar *c_rmta,*c_otha,*poa,*pda,*apvalues,*apvaluestmp,*c_rmtaa; 970 MPI_Comm comm; 971 972 PetscFunctionBegin; 973 ierr = PetscObjectGetComm((PetscObject)A,&comm);CHKERRQ(ierr); 974 if (!ptap) SETERRQ(comm,PETSC_ERR_ARG_WRONGSTATE,"PtAP cannot be reused. Do not call MatFreeIntermediateDataStructures() or use '-mat_freeintermediatedatastructures'"); 975 976 ierr = MatZeroEntries(C);CHKERRQ(ierr); 977 978 /* Get P_oth = ptap->P_oth and P_loc = ptap->P_loc */ 979 /*-----------------------------------------------------*/ 980 if (ptap->reuse == MAT_REUSE_MATRIX) { 981 /* P_oth and P_loc are obtained in MatPtASymbolic() when reuse == MAT_INITIAL_MATRIX */ 982 ierr = MatGetBrowsOfAoCols_MPIAIJ(A,P,MAT_REUSE_MATRIX,&ptap->startsj_s,&ptap->startsj_r,&ptap->bufa,&ptap->P_oth);CHKERRQ(ierr); 983 } 984 985 po = (Mat_SeqAIJ*) p->B->data; 986 pd = (Mat_SeqAIJ*) p->A->data; 987 988 ierr = MatGetLocalSize(p->B,NULL,&pon);CHKERRQ(ierr); 989 ierr = MatGetLocalSize(P,NULL,&pn);CHKERRQ(ierr); 990 991 ierr = PetscCalloc2(ptap->c_rmti[pon],&c_rmtj,ptap->c_rmti[pon],&c_rmta);CHKERRQ(ierr); 992 ierr = MatGetLocalSize(A,&am,NULL);CHKERRQ(ierr); 993 ierr = MatGetOwnershipRangeColumn(P,&pcstart,&pcend);CHKERRQ(ierr); 994 cmaxr = 0; 995 for (i=0; i<pon; i++) { 996 cmaxr = PetscMax(cmaxr,ptap->c_rmti[i+1]-ptap->c_rmti[i]); 997 } 998 cd = (Mat_SeqAIJ*)(c->A)->data; 999 co = (Mat_SeqAIJ*)(c->B)->data; 1000 for (i=0; i<pn; i++) { 1001 cmaxr = PetscMax(cmaxr,(cd->i[i+1]-cd->i[i])+(co->i[i+1]-co->i[i])); 1002 } 1003 ierr = PetscCalloc4(cmaxr,&apindices,cmaxr,&apvalues,cmaxr,&apvaluestmp,pon,&c_rmtc);CHKERRQ(ierr); 1004 ierr = PetscHMapIVCreate(&hmap);CHKERRQ(ierr); 1005 ierr = PetscHMapIVResize(hmap,cmaxr);CHKERRQ(ierr); 1006 for (i=0; i<am && (pon || pn); i++) { 1007 ierr = PetscHMapIVClear(hmap);CHKERRQ(ierr); 1008 nzi = po->i[i+1] - po->i[i]; 1009 dnzi = pd->i[i+1] - pd->i[i]; 1010 if (!nzi && !dnzi) continue; 1011 ierr = MatPtAPNumericComputeOneRowOfAP_private(A,P,ptap->P_oth,i,hmap);CHKERRQ(ierr); 1012 voff = 0; 1013 ierr = PetscHMapIVGetPairs(hmap,&voff,apindices,apvalues);CHKERRQ(ierr); 1014 if (!voff) continue; 1015 1016 /* Form remote C(ii, :) */ 1017 poj = po->j + po->i[i]; 1018 poa = po->a + po->i[i]; 1019 for (j=0; j<nzi; j++) { 1020 c_rmtjj = c_rmtj + ptap->c_rmti[poj[j]]; 1021 c_rmtaa = c_rmta + ptap->c_rmti[poj[j]]; 1022 for (jj=0; jj<voff; jj++) { 1023 apvaluestmp[jj] = apvalues[jj]*poa[j]; 1024 /*If the row is empty */ 1025 if (!c_rmtc[poj[j]]) { 1026 c_rmtjj[jj] = apindices[jj]; 1027 c_rmtaa[jj] = apvaluestmp[jj]; 1028 c_rmtc[poj[j]]++; 1029 } else { 1030 ierr = PetscFindInt(apindices[jj],c_rmtc[poj[j]],c_rmtjj,&loc);CHKERRQ(ierr); 1031 if (loc>=0){ /* hit */ 1032 c_rmtaa[loc] += apvaluestmp[jj]; 1033 } else { /* new element */ 1034 loc = -(loc+1); 1035 /* Move data backward */ 1036 for (kk=c_rmtc[poj[j]]; kk>loc; kk--) { 1037 c_rmtjj[kk] = c_rmtjj[kk-1]; 1038 c_rmtaa[kk] = c_rmtaa[kk-1]; 1039 }/* End kk */ 1040 c_rmtjj[loc] = apindices[jj]; 1041 c_rmtaa[loc] = apvaluestmp[jj]; 1042 c_rmtc[poj[j]]++; 1043 } 1044 } 1045 } /* End jj */ 1046 } /* End j */ 1047 1048 /* Form local C(ii, :) */ 1049 pdj = pd->j + pd->i[i]; 1050 pda = pd->a + pd->i[i]; 1051 for (j=0; j<dnzi; j++) { 1052 row = pcstart + pdj[j]; 1053 for (jj=0; jj<voff; jj++) { 1054 apvaluestmp[jj] = apvalues[jj]*pda[j]; 1055 }/* End kk */ 1056 ierr = MatSetValues(C,1,&row,voff,apindices,apvaluestmp,ADD_VALUES);CHKERRQ(ierr); 1057 }/* End j */ 1058 } /* End i */ 1059 1060 ierr = PetscFree4(apindices,apvalues,apvaluestmp,c_rmtc);CHKERRQ(ierr); 1061 ierr = PetscHMapIVDestroy(&hmap);CHKERRQ(ierr); 1062 ierr = PetscCalloc2(ptap->c_othi[pn],&c_othj,ptap->c_othi[pn],&c_otha);CHKERRQ(ierr); 1063 1064 ierr = PetscSFReduceBegin(ptap->sf,MPIU_INT,c_rmtj,c_othj,MPIU_REPLACE);CHKERRQ(ierr); 1065 ierr = PetscSFReduceBegin(ptap->sf,MPIU_SCALAR,c_rmta,c_otha,MPIU_REPLACE);CHKERRQ(ierr); 1066 ierr = PetscSFReduceEnd(ptap->sf,MPIU_INT,c_rmtj,c_othj,MPI_SUM);CHKERRQ(ierr); 1067 ierr = PetscSFReduceEnd(ptap->sf,MPIU_SCALAR,c_rmta,c_otha,MPIU_REPLACE);CHKERRQ(ierr); 1068 ierr = PetscFree2(c_rmtj,c_rmta);CHKERRQ(ierr); 1069 1070 /* Add contributions from remote */ 1071 for (i = 0; i < pn; i++) { 1072 row = i + pcstart; 1073 ierr = MatSetValues(C,1,&row,ptap->c_othi[i+1]-ptap->c_othi[i],c_othj+ptap->c_othi[i],c_otha+ptap->c_othi[i],ADD_VALUES);CHKERRQ(ierr); 1074 } 1075 ierr = PetscFree2(c_othj,c_otha);CHKERRQ(ierr); 1076 1077 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1078 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1079 1080 ptap->reuse = MAT_REUSE_MATRIX; 1081 1082 /* supporting struct ptap consumes almost same amount of memory as C=PtAP, release it if C will not be updated by A and P */ 1083 if (ptap->freestruct) { 1084 ierr = MatFreeIntermediateDataStructures(C);CHKERRQ(ierr); 1085 } 1086 PetscFunctionReturn(0); 1087 } 1088 1089 PetscErrorCode MatPtAPSymbolic_MPIAIJ_MPIAIJ_allatonce(Mat A,Mat P,PetscReal fill,Mat *C) 1090 { 1091 Mat_APMPI *ptap; 1092 Mat_MPIAIJ *p=(Mat_MPIAIJ*)P->data,*c; 1093 MPI_Comm comm; 1094 Mat Cmpi; 1095 Mat_SeqAIJ *pd,*po; 1096 MatType mtype; 1097 PetscSF sf; 1098 PetscSFNode *iremote; 1099 PetscInt rootspacesize,*rootspace,*rootspaceoffsets,nleaves; 1100 const PetscInt *rootdegrees; 1101 PetscHSetI ht,oht,*hta,*hto; 1102 PetscInt pn,pon,*c_rmtc,i,j,nzi,htsize,htosize,*c_rmtj,off,*c_othj,rcvncols,sendncols,*c_rmtoffsets; 1103 PetscInt owner,lidx,*rdj,col,pcstart,pcend,*dnz,*onz,am,arstart,arend,*poj,*pdj; 1104 PetscInt nalg=2,alg=0; 1105 PetscBool flg; 1106 const char *algTypes[2] = {"overlapping","merged"}; 1107 PetscErrorCode ierr; 1108 1109 PetscFunctionBegin; 1110 ierr = PetscObjectGetComm((PetscObject)A,&comm);CHKERRQ(ierr); 1111 1112 /* Create symbolic parallel matrix Cmpi */ 1113 ierr = MatGetLocalSize(P,NULL,&pn);CHKERRQ(ierr); 1114 ierr = MatCreate(comm,&Cmpi);CHKERRQ(ierr); 1115 ierr = MatGetType(A,&mtype);CHKERRQ(ierr); 1116 ierr = MatSetType(Cmpi,mtype);CHKERRQ(ierr); 1117 ierr = MatSetSizes(Cmpi,pn,pn,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 1118 1119 ierr = PetscNew(&ptap);CHKERRQ(ierr); 1120 ptap->reuse = MAT_INITIAL_MATRIX; 1121 ptap->algType = 2; 1122 1123 /* Get P_oth by taking rows of P (= non-zero cols of local A) from other processors */ 1124 ierr = MatGetBrowsOfAoCols_MPIAIJ(A,P,MAT_INITIAL_MATRIX,&ptap->startsj_s,&ptap->startsj_r,&ptap->bufa,&ptap->P_oth);CHKERRQ(ierr); 1125 1126 po = (Mat_SeqAIJ*)p->B->data; 1127 pd = (Mat_SeqAIJ*)p->A->data; 1128 1129 /* This equals to the number of offdiag columns in P */ 1130 ierr = MatGetLocalSize(p->B,NULL,&pon);CHKERRQ(ierr); 1131 /* offsets */ 1132 ierr = PetscMalloc1(pon+1,&ptap->c_rmti);CHKERRQ(ierr); 1133 /* The number of columns we will send to remote ranks */ 1134 ierr = PetscMalloc1(pon,&c_rmtc);CHKERRQ(ierr); 1135 ierr = PetscMalloc1(pon,&hta);CHKERRQ(ierr); 1136 for (i=0; i<pon; i++) { 1137 ierr = PetscHSetICreate(&hta[i]);CHKERRQ(ierr); 1138 } 1139 ierr = MatGetLocalSize(A,&am,NULL);CHKERRQ(ierr); 1140 ierr = MatGetOwnershipRange(A,&arstart,&arend);CHKERRQ(ierr); 1141 /* Create hash table to merge all columns for C(i, :) */ 1142 ierr = PetscHSetICreate(&ht);CHKERRQ(ierr); 1143 1144 ptap->c_rmti[0] = 0; 1145 /* 2) Pass 1: calculate the size for C_rmt (a matrix need to be sent to other processors) */ 1146 for (i=0; i<am && pon; i++) { 1147 /* Form one row of AP */ 1148 ierr = PetscHSetIClear(ht);CHKERRQ(ierr); 1149 /* If the off diag is empty, we should not do any calculation */ 1150 nzi = po->i[i+1] - po->i[i]; 1151 if (!nzi) continue; 1152 1153 ierr = MatPtAPSymbolicComputeOneRowOfAP_private(A,P,ptap->P_oth,i,ht,ht);CHKERRQ(ierr); 1154 ierr = PetscHSetIGetSize(ht,&htsize);CHKERRQ(ierr); 1155 /* If AP is empty, just continue */ 1156 if (!htsize) continue; 1157 /* Form C(ii, :) */ 1158 poj = po->j + po->i[i]; 1159 for (j=0; j<nzi; j++) { 1160 ierr = PetscHSetIAppend(hta[poj[j]],ht);CHKERRQ(ierr); 1161 } 1162 } 1163 1164 for (i=0; i<pon; i++) { 1165 ierr = PetscHSetIGetSize(hta[i],&htsize);CHKERRQ(ierr); 1166 ptap->c_rmti[i+1] = ptap->c_rmti[i] + htsize; 1167 c_rmtc[i] = htsize; 1168 } 1169 1170 ierr = PetscMalloc1(ptap->c_rmti[pon],&c_rmtj);CHKERRQ(ierr); 1171 1172 for (i=0; i<pon; i++) { 1173 off = 0; 1174 ierr = PetscHSetIGetElems(hta[i],&off,c_rmtj+ptap->c_rmti[i]);CHKERRQ(ierr); 1175 ierr = PetscHSetIDestroy(&hta[i]);CHKERRQ(ierr); 1176 } 1177 ierr = PetscFree(hta);CHKERRQ(ierr); 1178 1179 ierr = PetscCalloc1(pon,&iremote);CHKERRQ(ierr); 1180 for (i=0; i<pon; i++) { 1181 ierr = PetscLayoutFindOwnerIndex(P->cmap,p->garray[i],&owner,&lidx);CHKERRQ(ierr); 1182 iremote[i].index = lidx; 1183 iremote[i].rank = owner; 1184 } 1185 1186 ierr = PetscSFCreate(comm,&sf);CHKERRQ(ierr); 1187 ierr = PetscSFSetGraph(sf,pn,pon,NULL,PETSC_OWN_POINTER,iremote,PETSC_OWN_POINTER);CHKERRQ(ierr); 1188 /* Reorder ranks properly so that the data handled by gather and scatter have the same order */ 1189 ierr = PetscSFSetRankOrder(sf,PETSC_TRUE);CHKERRQ(ierr); 1190 ierr = PetscSFSetFromOptions(sf);CHKERRQ(ierr); 1191 ierr = PetscSFSetUp(sf);CHKERRQ(ierr); 1192 /* How many neighbors have contributions to my rows? */ 1193 ierr = PetscSFComputeDegreeBegin(sf,&rootdegrees);CHKERRQ(ierr); 1194 ierr = PetscSFComputeDegreeEnd(sf,&rootdegrees);CHKERRQ(ierr); 1195 rootspacesize = 0; 1196 for (i = 0; i < pn; i++) { 1197 rootspacesize += rootdegrees[i]; 1198 } 1199 ierr = PetscCalloc1(rootspacesize,&rootspace);CHKERRQ(ierr); 1200 ierr = PetscCalloc1(rootspacesize+1,&rootspaceoffsets);CHKERRQ(ierr); 1201 /* Get information from leaves 1202 * Number of columns other people contribute to my rows 1203 * */ 1204 ierr = PetscSFGatherBegin(sf,MPIU_INT,c_rmtc,rootspace);CHKERRQ(ierr); 1205 ierr = PetscSFGatherEnd(sf,MPIU_INT,c_rmtc,rootspace);CHKERRQ(ierr); 1206 ierr = PetscFree(c_rmtc);CHKERRQ(ierr); 1207 ierr = PetscCalloc1(pn+1,&ptap->c_othi);CHKERRQ(ierr); 1208 /* The number of columns is received for each row */ 1209 ptap->c_othi[0] = 0; 1210 rootspacesize = 0; 1211 rootspaceoffsets[0] = 0; 1212 for (i = 0; i < pn; i++) { 1213 rcvncols = 0; 1214 for (j = 0; j<rootdegrees[i]; j++) { 1215 rcvncols += rootspace[rootspacesize]; 1216 rootspaceoffsets[rootspacesize+1] = rootspaceoffsets[rootspacesize] + rootspace[rootspacesize]; 1217 rootspacesize++; 1218 } 1219 ptap->c_othi[i+1] = ptap->c_othi[i] + rcvncols; 1220 } 1221 ierr = PetscFree(rootspace);CHKERRQ(ierr); 1222 1223 ierr = PetscCalloc1(pon,&c_rmtoffsets);CHKERRQ(ierr); 1224 ierr = PetscSFScatterBegin(sf,MPIU_INT,rootspaceoffsets,c_rmtoffsets);CHKERRQ(ierr); 1225 ierr = PetscSFScatterEnd(sf,MPIU_INT,rootspaceoffsets,c_rmtoffsets);CHKERRQ(ierr); 1226 ierr = PetscSFDestroy(&sf);CHKERRQ(ierr); 1227 ierr = PetscFree(rootspaceoffsets);CHKERRQ(ierr); 1228 1229 ierr = PetscCalloc1(ptap->c_rmti[pon],&iremote);CHKERRQ(ierr); 1230 nleaves = 0; 1231 for (i = 0; i<pon; i++) { 1232 ierr = PetscLayoutFindOwnerIndex(P->cmap,p->garray[i],&owner,&lidx);CHKERRQ(ierr); 1233 sendncols = ptap->c_rmti[i+1] - ptap->c_rmti[i]; 1234 for (j=0; j<sendncols; j++) { 1235 iremote[nleaves].rank = owner; 1236 iremote[nleaves++].index = c_rmtoffsets[i] + j; 1237 } 1238 } 1239 ierr = PetscFree(c_rmtoffsets);CHKERRQ(ierr); 1240 ierr = PetscCalloc1(ptap->c_othi[pn],&c_othj);CHKERRQ(ierr); 1241 1242 ierr = PetscSFCreate(comm,&ptap->sf);CHKERRQ(ierr); 1243 ierr = PetscSFSetGraph(ptap->sf,ptap->c_othi[pn],nleaves,NULL,PETSC_OWN_POINTER,iremote,PETSC_OWN_POINTER);CHKERRQ(ierr); 1244 ierr = PetscSFSetFromOptions(ptap->sf);CHKERRQ(ierr); 1245 /* One to one map */ 1246 ierr = PetscSFReduceBegin(ptap->sf,MPIU_INT,c_rmtj,c_othj,MPIU_REPLACE);CHKERRQ(ierr); 1247 1248 ierr = PetscCalloc2(pn,&dnz,pn,&onz);CHKERRQ(ierr); 1249 ierr = PetscHSetICreate(&oht);CHKERRQ(ierr); 1250 ierr = MatGetOwnershipRangeColumn(P,&pcstart,&pcend);CHKERRQ(ierr); 1251 ierr = PetscMalloc2(pn,&hta,pn,&hto);CHKERRQ(ierr); 1252 for (i=0; i<pn; i++) { 1253 ierr = PetscHSetICreate(&hta[i]);CHKERRQ(ierr); 1254 ierr = PetscHSetICreate(&hto[i]);CHKERRQ(ierr); 1255 } 1256 /* Work on local part */ 1257 /* 4) Pass 1: Estimate memory for C_loc */ 1258 for (i=0; i<am && pn; i++) { 1259 ierr = PetscHSetIClear(ht);CHKERRQ(ierr); 1260 ierr = PetscHSetIClear(oht);CHKERRQ(ierr); 1261 nzi = pd->i[i+1] - pd->i[i]; 1262 if (!nzi) continue; 1263 1264 ierr = MatPtAPSymbolicComputeOneRowOfAP_private(A,P,ptap->P_oth,i,ht,oht);CHKERRQ(ierr); 1265 ierr = PetscHSetIGetSize(ht,&htsize);CHKERRQ(ierr); 1266 ierr = PetscHSetIGetSize(oht,&htosize);CHKERRQ(ierr); 1267 if (!(htsize+htosize)) continue; 1268 /* Form C(ii, :) */ 1269 pdj = pd->j + pd->i[i]; 1270 for (j=0; j<nzi; j++) { 1271 ierr = PetscHSetIAppend(hta[pdj[j]],ht);CHKERRQ(ierr); 1272 ierr = PetscHSetIAppend(hto[pdj[j]],oht);CHKERRQ(ierr); 1273 } 1274 } 1275 1276 ierr = PetscHSetIDestroy(&ht);CHKERRQ(ierr); 1277 ierr = PetscHSetIDestroy(&oht);CHKERRQ(ierr); 1278 1279 /* Get remote data */ 1280 ierr = PetscSFReduceEnd(ptap->sf,MPIU_INT,c_rmtj,c_othj,MPIU_REPLACE);CHKERRQ(ierr); 1281 ierr = PetscFree(c_rmtj);CHKERRQ(ierr); 1282 1283 for (i = 0; i < pn; i++) { 1284 nzi = ptap->c_othi[i+1] - ptap->c_othi[i]; 1285 rdj = c_othj + ptap->c_othi[i]; 1286 for (j = 0; j < nzi; j++) { 1287 col = rdj[j]; 1288 /* diag part */ 1289 if (col>=pcstart && col<pcend) { 1290 ierr = PetscHSetIAdd(hta[i],col);CHKERRQ(ierr); 1291 } else { /* off diag */ 1292 ierr = PetscHSetIAdd(hto[i],col);CHKERRQ(ierr); 1293 } 1294 } 1295 ierr = PetscHSetIGetSize(hta[i],&htsize);CHKERRQ(ierr); 1296 dnz[i] = htsize; 1297 ierr = PetscHSetIDestroy(&hta[i]);CHKERRQ(ierr); 1298 ierr = PetscHSetIGetSize(hto[i],&htsize);CHKERRQ(ierr); 1299 onz[i] = htsize; 1300 ierr = PetscHSetIDestroy(&hto[i]);CHKERRQ(ierr); 1301 } 1302 1303 ierr = PetscFree2(hta,hto);CHKERRQ(ierr); 1304 ierr = PetscFree(c_othj);CHKERRQ(ierr); 1305 1306 /* local sizes and preallocation */ 1307 ierr = MatSetSizes(Cmpi,pn,pn,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 1308 ierr = MatSetBlockSizes(Cmpi,PetscAbs(P->cmap->bs),PetscAbs(P->cmap->bs));CHKERRQ(ierr); 1309 ierr = MatMPIAIJSetPreallocation(Cmpi,0,dnz,0,onz);CHKERRQ(ierr); 1310 ierr = PetscFree2(dnz,onz);CHKERRQ(ierr); 1311 1312 /* attach the supporting struct to Cmpi for reuse */ 1313 c = (Mat_MPIAIJ*)Cmpi->data; 1314 c->ap = ptap; 1315 ptap->duplicate = Cmpi->ops->duplicate; 1316 ptap->destroy = Cmpi->ops->destroy; 1317 ptap->view = Cmpi->ops->view; 1318 1319 /* Cmpi is not ready for use - assembly will be done by MatPtAPNumeric() */ 1320 Cmpi->assembled = PETSC_FALSE; 1321 /* pick an algorithm */ 1322 ierr = PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"MatPtAP","Mat");CHKERRQ(ierr); 1323 alg = 0; 1324 ierr = PetscOptionsEList("-matptap_allatonce_via","PtAP allatonce numeric approach","MatPtAP",algTypes,nalg,algTypes[alg],&alg,&flg);CHKERRQ(ierr); 1325 ierr = PetscOptionsEnd();CHKERRQ(ierr); 1326 switch (alg) { 1327 case 0: 1328 Cmpi->ops->ptapnumeric = MatPtAPNumeric_MPIAIJ_MPIAIJ_allatonce; 1329 break; 1330 case 1: 1331 Cmpi->ops->ptapnumeric = MatPtAPNumeric_MPIAIJ_MPIAIJ_allatonce_merged; 1332 break; 1333 default: 1334 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG," Unsupported allatonce numerical algorithm \n"); 1335 } 1336 Cmpi->ops->destroy = MatDestroy_MPIAIJ_PtAP; 1337 Cmpi->ops->view = MatView_MPIAIJ_PtAP; 1338 Cmpi->ops->freeintermediatedatastructures = MatFreeIntermediateDataStructures_MPIAIJ_AP; 1339 *C = Cmpi; 1340 PetscFunctionReturn(0); 1341 } 1342 1343 PetscErrorCode MatPtAPSymbolic_MPIAIJ_MPIAIJ_allatonce_merged(Mat A,Mat P,PetscReal fill,Mat *C) 1344 { 1345 Mat_APMPI *ptap; 1346 Mat_MPIAIJ *p=(Mat_MPIAIJ*)P->data,*c; 1347 MPI_Comm comm; 1348 Mat Cmpi; 1349 Mat_SeqAIJ *pd,*po; 1350 MatType mtype; 1351 PetscSF sf; 1352 PetscSFNode *iremote; 1353 PetscInt rootspacesize,*rootspace,*rootspaceoffsets,nleaves; 1354 const PetscInt *rootdegrees; 1355 PetscHSetI ht,oht,*hta,*hto,*htd; 1356 PetscInt pn,pon,*c_rmtc,i,j,nzi,dnzi,htsize,htosize,*c_rmtj,off,*c_othj,rcvncols,sendncols,*c_rmtoffsets; 1357 PetscInt owner,lidx,*rdj,col,pcstart,pcend,*dnz,*onz,am,arstart,arend,*poj,*pdj; 1358 PetscInt nalg=2,alg=0; 1359 PetscBool flg; 1360 const char *algTypes[2] = {"merged","overlapping"}; 1361 PetscErrorCode ierr; 1362 1363 PetscFunctionBegin; 1364 ierr = PetscObjectGetComm((PetscObject)A,&comm);CHKERRQ(ierr); 1365 1366 /* Create symbolic parallel matrix Cmpi */ 1367 ierr = MatGetLocalSize(P,NULL,&pn);CHKERRQ(ierr); 1368 ierr = MatCreate(comm,&Cmpi);CHKERRQ(ierr); 1369 ierr = MatGetType(A,&mtype);CHKERRQ(ierr); 1370 ierr = MatSetType(Cmpi,mtype);CHKERRQ(ierr); 1371 ierr = MatSetSizes(Cmpi,pn,pn,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 1372 1373 ierr = PetscNew(&ptap);CHKERRQ(ierr); 1374 ptap->reuse = MAT_INITIAL_MATRIX; 1375 ptap->algType = 3; 1376 1377 /* 0) Get P_oth by taking rows of P (= non-zero cols of local A) from other processors */ 1378 ierr = MatGetBrowsOfAoCols_MPIAIJ(A,P,MAT_INITIAL_MATRIX,&ptap->startsj_s,&ptap->startsj_r,&ptap->bufa,&ptap->P_oth);CHKERRQ(ierr); 1379 1380 po = (Mat_SeqAIJ*)p->B->data; 1381 pd = (Mat_SeqAIJ*)p->A->data; 1382 1383 /* This equals to the number of offdiag columns in P */ 1384 ierr = MatGetLocalSize(p->B,NULL,&pon);CHKERRQ(ierr); 1385 /* offsets */ 1386 ierr = PetscMalloc1(pon+1,&ptap->c_rmti);CHKERRQ(ierr); 1387 /* The number of columns we will send to remote ranks */ 1388 ierr = PetscMalloc1(pon,&c_rmtc);CHKERRQ(ierr); 1389 ierr = PetscMalloc1(pon,&hta);CHKERRQ(ierr); 1390 for (i=0; i<pon; i++) { 1391 ierr = PetscHSetICreate(&hta[i]);CHKERRQ(ierr); 1392 } 1393 ierr = MatGetLocalSize(A,&am,NULL);CHKERRQ(ierr); 1394 ierr = MatGetOwnershipRange(A,&arstart,&arend);CHKERRQ(ierr); 1395 /* Create hash table to merge all columns for C(i, :) */ 1396 ierr = PetscHSetICreate(&ht);CHKERRQ(ierr); 1397 ierr = PetscHSetICreate(&oht);CHKERRQ(ierr); 1398 ierr = PetscMalloc2(pn,&htd,pn,&hto);CHKERRQ(ierr); 1399 for (i=0; i<pn; i++) { 1400 ierr = PetscHSetICreate(&htd[i]);CHKERRQ(ierr); 1401 ierr = PetscHSetICreate(&hto[i]);CHKERRQ(ierr); 1402 } 1403 ptap->c_rmti[0] = 0; 1404 /* 2) Pass 1: calculate the size for C_rmt (a matrix need to be sent to other processors) */ 1405 for (i=0; i<am && (pon || pn); i++) { 1406 /* Form one row of AP */ 1407 ierr = PetscHSetIClear(ht);CHKERRQ(ierr); 1408 ierr = PetscHSetIClear(oht);CHKERRQ(ierr); 1409 /* If the off diag is empty, we should not do any calculation */ 1410 nzi = po->i[i+1] - po->i[i]; 1411 dnzi = pd->i[i+1] - pd->i[i]; 1412 if (!nzi && !dnzi) continue; 1413 1414 ierr = MatPtAPSymbolicComputeOneRowOfAP_private(A,P,ptap->P_oth,i,ht,oht);CHKERRQ(ierr); 1415 ierr = PetscHSetIGetSize(ht,&htsize);CHKERRQ(ierr); 1416 ierr = PetscHSetIGetSize(oht,&htosize);CHKERRQ(ierr); 1417 /* If AP is empty, just continue */ 1418 if (!(htsize+htosize)) continue; 1419 1420 /* Form remote C(ii, :) */ 1421 poj = po->j + po->i[i]; 1422 for (j=0; j<nzi; j++) { 1423 ierr = PetscHSetIAppend(hta[poj[j]],ht);CHKERRQ(ierr); 1424 ierr = PetscHSetIAppend(hta[poj[j]],oht);CHKERRQ(ierr); 1425 } 1426 1427 /* Form local C(ii, :) */ 1428 pdj = pd->j + pd->i[i]; 1429 for (j=0; j<dnzi; j++) { 1430 ierr = PetscHSetIAppend(htd[pdj[j]],ht);CHKERRQ(ierr); 1431 ierr = PetscHSetIAppend(hto[pdj[j]],oht);CHKERRQ(ierr); 1432 } 1433 } 1434 1435 ierr = PetscHSetIDestroy(&ht);CHKERRQ(ierr); 1436 ierr = PetscHSetIDestroy(&oht);CHKERRQ(ierr); 1437 1438 for (i=0; i<pon; i++) { 1439 ierr = PetscHSetIGetSize(hta[i],&htsize);CHKERRQ(ierr); 1440 ptap->c_rmti[i+1] = ptap->c_rmti[i] + htsize; 1441 c_rmtc[i] = htsize; 1442 } 1443 1444 ierr = PetscMalloc1(ptap->c_rmti[pon],&c_rmtj);CHKERRQ(ierr); 1445 1446 for (i=0; i<pon; i++) { 1447 off = 0; 1448 ierr = PetscHSetIGetElems(hta[i],&off,c_rmtj+ptap->c_rmti[i]);CHKERRQ(ierr); 1449 ierr = PetscHSetIDestroy(&hta[i]);CHKERRQ(ierr); 1450 } 1451 ierr = PetscFree(hta);CHKERRQ(ierr); 1452 1453 ierr = PetscCalloc1(pon,&iremote);CHKERRQ(ierr); 1454 for (i=0; i<pon; i++) { 1455 ierr = PetscLayoutFindOwnerIndex(P->cmap,p->garray[i],&owner,&lidx);CHKERRQ(ierr); 1456 iremote[i].index = lidx; 1457 iremote[i].rank = owner; 1458 } 1459 1460 ierr = PetscSFCreate(comm,&sf);CHKERRQ(ierr); 1461 ierr = PetscSFSetGraph(sf,pn,pon,NULL,PETSC_OWN_POINTER,iremote,PETSC_OWN_POINTER);CHKERRQ(ierr); 1462 /* Reorder ranks properly so that the data handled by gather and scatter have the same order */ 1463 ierr = PetscSFSetRankOrder(sf,PETSC_TRUE);CHKERRQ(ierr); 1464 ierr = PetscSFSetFromOptions(sf);CHKERRQ(ierr); 1465 ierr = PetscSFSetUp(sf);CHKERRQ(ierr); 1466 /* How many neighbors have contributions to my rows? */ 1467 ierr = PetscSFComputeDegreeBegin(sf,&rootdegrees);CHKERRQ(ierr); 1468 ierr = PetscSFComputeDegreeEnd(sf,&rootdegrees);CHKERRQ(ierr); 1469 rootspacesize = 0; 1470 for (i = 0; i < pn; i++) { 1471 rootspacesize += rootdegrees[i]; 1472 } 1473 ierr = PetscCalloc1(rootspacesize,&rootspace);CHKERRQ(ierr); 1474 ierr = PetscCalloc1(rootspacesize+1,&rootspaceoffsets);CHKERRQ(ierr); 1475 /* Get information from leaves 1476 * Number of columns other people contribute to my rows 1477 * */ 1478 ierr = PetscSFGatherBegin(sf,MPIU_INT,c_rmtc,rootspace);CHKERRQ(ierr); 1479 ierr = PetscSFGatherEnd(sf,MPIU_INT,c_rmtc,rootspace);CHKERRQ(ierr); 1480 ierr = PetscFree(c_rmtc);CHKERRQ(ierr); 1481 ierr = PetscCalloc1(pn+1,&ptap->c_othi);CHKERRQ(ierr); 1482 /* The number of columns is received for each row */ 1483 ptap->c_othi[0] = 0; 1484 rootspacesize = 0; 1485 rootspaceoffsets[0] = 0; 1486 for (i = 0; i < pn; i++) { 1487 rcvncols = 0; 1488 for (j = 0; j<rootdegrees[i]; j++) { 1489 rcvncols += rootspace[rootspacesize]; 1490 rootspaceoffsets[rootspacesize+1] = rootspaceoffsets[rootspacesize] + rootspace[rootspacesize]; 1491 rootspacesize++; 1492 } 1493 ptap->c_othi[i+1] = ptap->c_othi[i] + rcvncols; 1494 } 1495 ierr = PetscFree(rootspace);CHKERRQ(ierr); 1496 1497 ierr = PetscCalloc1(pon,&c_rmtoffsets);CHKERRQ(ierr); 1498 ierr = PetscSFScatterBegin(sf,MPIU_INT,rootspaceoffsets,c_rmtoffsets);CHKERRQ(ierr); 1499 ierr = PetscSFScatterEnd(sf,MPIU_INT,rootspaceoffsets,c_rmtoffsets);CHKERRQ(ierr); 1500 ierr = PetscSFDestroy(&sf);CHKERRQ(ierr); 1501 ierr = PetscFree(rootspaceoffsets);CHKERRQ(ierr); 1502 1503 ierr = PetscCalloc1(ptap->c_rmti[pon],&iremote);CHKERRQ(ierr); 1504 nleaves = 0; 1505 for (i = 0; i<pon; i++) { 1506 ierr = PetscLayoutFindOwnerIndex(P->cmap,p->garray[i],&owner,&lidx);CHKERRQ(ierr); 1507 sendncols = ptap->c_rmti[i+1] - ptap->c_rmti[i]; 1508 for (j=0; j<sendncols; j++) { 1509 iremote[nleaves].rank = owner; 1510 iremote[nleaves++].index = c_rmtoffsets[i] + j; 1511 } 1512 } 1513 ierr = PetscFree(c_rmtoffsets);CHKERRQ(ierr); 1514 ierr = PetscCalloc1(ptap->c_othi[pn],&c_othj);CHKERRQ(ierr); 1515 1516 ierr = PetscSFCreate(comm,&ptap->sf);CHKERRQ(ierr); 1517 ierr = PetscSFSetGraph(ptap->sf,ptap->c_othi[pn],nleaves,NULL,PETSC_OWN_POINTER,iremote,PETSC_OWN_POINTER);CHKERRQ(ierr); 1518 ierr = PetscSFSetFromOptions(ptap->sf);CHKERRQ(ierr); 1519 /* One to one map */ 1520 ierr = PetscSFReduceBegin(ptap->sf,MPIU_INT,c_rmtj,c_othj,MPIU_REPLACE);CHKERRQ(ierr); 1521 /* Get remote data */ 1522 ierr = PetscSFReduceEnd(ptap->sf,MPIU_INT,c_rmtj,c_othj,MPIU_REPLACE);CHKERRQ(ierr); 1523 ierr = PetscFree(c_rmtj);CHKERRQ(ierr); 1524 ierr = PetscCalloc2(pn,&dnz,pn,&onz);CHKERRQ(ierr); 1525 ierr = MatGetOwnershipRangeColumn(P,&pcstart,&pcend);CHKERRQ(ierr); 1526 1527 for (i = 0; i < pn; i++) { 1528 nzi = ptap->c_othi[i+1] - ptap->c_othi[i]; 1529 rdj = c_othj + ptap->c_othi[i]; 1530 for (j = 0; j < nzi; j++) { 1531 col = rdj[j]; 1532 /* diag part */ 1533 if (col>=pcstart && col<pcend) { 1534 ierr = PetscHSetIAdd(htd[i],col);CHKERRQ(ierr); 1535 } else { /* off diag */ 1536 ierr = PetscHSetIAdd(hto[i],col);CHKERRQ(ierr); 1537 } 1538 } 1539 ierr = PetscHSetIGetSize(htd[i],&htsize);CHKERRQ(ierr); 1540 dnz[i] = htsize; 1541 ierr = PetscHSetIDestroy(&htd[i]);CHKERRQ(ierr); 1542 ierr = PetscHSetIGetSize(hto[i],&htsize);CHKERRQ(ierr); 1543 onz[i] = htsize; 1544 ierr = PetscHSetIDestroy(&hto[i]);CHKERRQ(ierr); 1545 } 1546 1547 ierr = PetscFree2(htd,hto);CHKERRQ(ierr); 1548 ierr = PetscFree(c_othj);CHKERRQ(ierr); 1549 1550 /* local sizes and preallocation */ 1551 ierr = MatSetSizes(Cmpi,pn,pn,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 1552 ierr = MatSetBlockSizes(Cmpi,PetscAbs(P->cmap->bs),PetscAbs(P->cmap->bs));CHKERRQ(ierr); 1553 ierr = MatMPIAIJSetPreallocation(Cmpi,0,dnz,0,onz);CHKERRQ(ierr); 1554 ierr = PetscFree2(dnz,onz);CHKERRQ(ierr); 1555 1556 /* attach the supporting struct to Cmpi for reuse */ 1557 c = (Mat_MPIAIJ*)Cmpi->data; 1558 c->ap = ptap; 1559 ptap->duplicate = Cmpi->ops->duplicate; 1560 ptap->destroy = Cmpi->ops->destroy; 1561 ptap->view = Cmpi->ops->view; 1562 1563 /* Cmpi is not ready for use - assembly will be done by MatPtAPNumeric() */ 1564 Cmpi->assembled = PETSC_FALSE; 1565 /* pick an algorithm */ 1566 ierr = PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"MatPtAP","Mat");CHKERRQ(ierr); 1567 alg = 0; 1568 ierr = PetscOptionsEList("-matptap_allatonce_via","PtAP allatonce numeric approach","MatPtAP",algTypes,nalg,algTypes[alg],&alg,&flg);CHKERRQ(ierr); 1569 ierr = PetscOptionsEnd();CHKERRQ(ierr); 1570 switch (alg) { 1571 case 0: 1572 Cmpi->ops->ptapnumeric = MatPtAPNumeric_MPIAIJ_MPIAIJ_allatonce_merged; 1573 break; 1574 case 1: 1575 Cmpi->ops->ptapnumeric = MatPtAPNumeric_MPIAIJ_MPIAIJ_allatonce; 1576 break; 1577 default: 1578 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG," Unsupported allatonce numerical algorithm \n"); 1579 } 1580 Cmpi->ops->destroy = MatDestroy_MPIAIJ_PtAP; 1581 Cmpi->ops->view = MatView_MPIAIJ_PtAP; 1582 Cmpi->ops->freeintermediatedatastructures = MatFreeIntermediateDataStructures_MPIAIJ_AP; 1583 *C = Cmpi; 1584 PetscFunctionReturn(0); 1585 } 1586 1587 PetscErrorCode MatPtAPSymbolic_MPIAIJ_MPIAIJ(Mat A,Mat P,PetscReal fill,Mat *C) 1588 { 1589 PetscErrorCode ierr; 1590 Mat_APMPI *ptap; 1591 Mat_MPIAIJ *a=(Mat_MPIAIJ*)A->data,*p=(Mat_MPIAIJ*)P->data,*c; 1592 MPI_Comm comm; 1593 PetscMPIInt size,rank; 1594 Mat Cmpi; 1595 PetscFreeSpaceList free_space=NULL,current_space=NULL; 1596 PetscInt am=A->rmap->n,pm=P->rmap->n,pN=P->cmap->N,pn=P->cmap->n; 1597 PetscInt *lnk,i,k,pnz,row,nsend; 1598 PetscBT lnkbt; 1599 PetscMPIInt tagi,tagj,*len_si,*len_s,*len_ri,icompleted=0,nrecv; 1600 PetscInt **buf_rj,**buf_ri,**buf_ri_k; 1601 PetscInt len,proc,*dnz,*onz,*owners,nzi,nspacedouble; 1602 PetscInt nrows,*buf_s,*buf_si,*buf_si_i,**nextrow,**nextci; 1603 MPI_Request *swaits,*rwaits; 1604 MPI_Status *sstatus,rstatus; 1605 PetscLayout rowmap; 1606 PetscInt *owners_co,*coi,*coj; /* i and j array of (p->B)^T*A*P - used in the communication */ 1607 PetscMPIInt *len_r,*id_r; /* array of length of comm->size, store send/recv matrix values */ 1608 PetscInt *api,*apj,*Jptr,apnz,*prmap=p->garray,con,j,ap_rmax=0,Crmax,*aj,*ai,*pi; 1609 Mat_SeqAIJ *p_loc,*p_oth=NULL,*ad=(Mat_SeqAIJ*)(a->A)->data,*ao=NULL,*c_loc,*c_oth; 1610 PetscScalar *apv; 1611 PetscTable ta; 1612 MatType mtype; 1613 const char *prefix; 1614 #if defined(PETSC_USE_INFO) 1615 PetscReal apfill; 1616 #endif 1617 1618 PetscFunctionBegin; 1619 ierr = PetscObjectGetComm((PetscObject)A,&comm);CHKERRQ(ierr); 1620 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 1621 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 1622 1623 if (size > 1) ao = (Mat_SeqAIJ*)(a->B)->data; 1624 1625 /* create symbolic parallel matrix Cmpi */ 1626 ierr = MatCreate(comm,&Cmpi);CHKERRQ(ierr); 1627 ierr = MatGetType(A,&mtype);CHKERRQ(ierr); 1628 ierr = MatSetType(Cmpi,mtype);CHKERRQ(ierr); 1629 1630 /* Do dense axpy in MatPtAPNumeric_MPIAIJ_MPIAIJ() */ 1631 Cmpi->ops->ptapnumeric = MatPtAPNumeric_MPIAIJ_MPIAIJ; 1632 1633 /* create struct Mat_APMPI and attached it to C later */ 1634 ierr = PetscNew(&ptap);CHKERRQ(ierr); 1635 ptap->reuse = MAT_INITIAL_MATRIX; 1636 ptap->algType = 1; 1637 1638 /* get P_oth by taking rows of P (= non-zero cols of local A) from other processors */ 1639 ierr = MatGetBrowsOfAoCols_MPIAIJ(A,P,MAT_INITIAL_MATRIX,&ptap->startsj_s,&ptap->startsj_r,&ptap->bufa,&ptap->P_oth);CHKERRQ(ierr); 1640 /* get P_loc by taking all local rows of P */ 1641 ierr = MatMPIAIJGetLocalMat(P,MAT_INITIAL_MATRIX,&ptap->P_loc);CHKERRQ(ierr); 1642 1643 /* (0) compute Rd = Pd^T, Ro = Po^T */ 1644 /* --------------------------------- */ 1645 ierr = MatTranspose(p->A,MAT_INITIAL_MATRIX,&ptap->Rd);CHKERRQ(ierr); 1646 ierr = MatTranspose(p->B,MAT_INITIAL_MATRIX,&ptap->Ro);CHKERRQ(ierr); 1647 1648 /* (1) compute symbolic AP = A_loc*P = Ad*P_loc + Ao*P_oth (api,apj) */ 1649 /* ----------------------------------------------------------------- */ 1650 p_loc = (Mat_SeqAIJ*)(ptap->P_loc)->data; 1651 if (ptap->P_oth) p_oth = (Mat_SeqAIJ*)(ptap->P_oth)->data; 1652 1653 /* create and initialize a linked list */ 1654 ierr = PetscTableCreate(pn,pN,&ta);CHKERRQ(ierr); /* for compute AP_loc and Cmpi */ 1655 MatRowMergeMax_SeqAIJ(p_loc,ptap->P_loc->rmap->N,ta); 1656 MatRowMergeMax_SeqAIJ(p_oth,ptap->P_oth->rmap->N,ta); 1657 ierr = PetscTableGetCount(ta,&Crmax);CHKERRQ(ierr); /* Crmax = nnz(sum of Prows) */ 1658 /* printf("[%d] est %d, Crmax %d; pN %d\n",rank,5*(p_loc->rmax+p_oth->rmax + (PetscInt)(1.e-2*pN)),Crmax,pN); */ 1659 1660 ierr = PetscLLCondensedCreate(Crmax,pN,&lnk,&lnkbt);CHKERRQ(ierr); 1661 1662 /* Initial FreeSpace size is fill*(nnz(A) + nnz(P)) */ 1663 if (ao) { 1664 ierr = PetscFreeSpaceGet(PetscRealIntMultTruncate(fill,PetscIntSumTruncate(ad->i[am],PetscIntSumTruncate(ao->i[am],p_loc->i[pm]))),&free_space);CHKERRQ(ierr); 1665 } else { 1666 ierr = PetscFreeSpaceGet(PetscRealIntMultTruncate(fill,PetscIntSumTruncate(ad->i[am],p_loc->i[pm])),&free_space);CHKERRQ(ierr); 1667 } 1668 current_space = free_space; 1669 nspacedouble = 0; 1670 1671 ierr = PetscMalloc1(am+1,&api);CHKERRQ(ierr); 1672 api[0] = 0; 1673 for (i=0; i<am; i++) { 1674 /* diagonal portion: Ad[i,:]*P */ 1675 ai = ad->i; pi = p_loc->i; 1676 nzi = ai[i+1] - ai[i]; 1677 aj = ad->j + ai[i]; 1678 for (j=0; j<nzi; j++) { 1679 row = aj[j]; 1680 pnz = pi[row+1] - pi[row]; 1681 Jptr = p_loc->j + pi[row]; 1682 /* add non-zero cols of P into the sorted linked list lnk */ 1683 ierr = PetscLLCondensedAddSorted(pnz,Jptr,lnk,lnkbt);CHKERRQ(ierr); 1684 } 1685 /* off-diagonal portion: Ao[i,:]*P */ 1686 if (ao) { 1687 ai = ao->i; pi = p_oth->i; 1688 nzi = ai[i+1] - ai[i]; 1689 aj = ao->j + ai[i]; 1690 for (j=0; j<nzi; j++) { 1691 row = aj[j]; 1692 pnz = pi[row+1] - pi[row]; 1693 Jptr = p_oth->j + pi[row]; 1694 ierr = PetscLLCondensedAddSorted(pnz,Jptr,lnk,lnkbt);CHKERRQ(ierr); 1695 } 1696 } 1697 apnz = lnk[0]; 1698 api[i+1] = api[i] + apnz; 1699 if (ap_rmax < apnz) ap_rmax = apnz; 1700 1701 /* if free space is not available, double the total space in the list */ 1702 if (current_space->local_remaining<apnz) { 1703 ierr = PetscFreeSpaceGet(PetscIntSumTruncate(apnz,current_space->total_array_size),¤t_space);CHKERRQ(ierr); 1704 nspacedouble++; 1705 } 1706 1707 /* Copy data into free space, then initialize lnk */ 1708 ierr = PetscLLCondensedClean(pN,apnz,current_space->array,lnk,lnkbt);CHKERRQ(ierr); 1709 1710 current_space->array += apnz; 1711 current_space->local_used += apnz; 1712 current_space->local_remaining -= apnz; 1713 } 1714 /* Allocate space for apj and apv, initialize apj, and */ 1715 /* destroy list of free space and other temporary array(s) */ 1716 ierr = PetscMalloc2(api[am],&apj,api[am],&apv);CHKERRQ(ierr); 1717 ierr = PetscFreeSpaceContiguous(&free_space,apj);CHKERRQ(ierr); 1718 ierr = PetscLLDestroy(lnk,lnkbt);CHKERRQ(ierr); 1719 1720 /* Create AP_loc for reuse */ 1721 ierr = MatCreateSeqAIJWithArrays(PETSC_COMM_SELF,am,pN,api,apj,apv,&ptap->AP_loc);CHKERRQ(ierr); 1722 1723 #if defined(PETSC_USE_INFO) 1724 if (ao) { 1725 apfill = (PetscReal)api[am]/(ad->i[am]+ao->i[am]+p_loc->i[pm]+1); 1726 } else { 1727 apfill = (PetscReal)api[am]/(ad->i[am]+p_loc->i[pm]+1); 1728 } 1729 ptap->AP_loc->info.mallocs = nspacedouble; 1730 ptap->AP_loc->info.fill_ratio_given = fill; 1731 ptap->AP_loc->info.fill_ratio_needed = apfill; 1732 1733 if (api[am]) { 1734 ierr = PetscInfo3(ptap->AP_loc,"Nonscalable algorithm, AP_loc reallocs %D; Fill ratio: given %g needed %g.\n",nspacedouble,(double)fill,(double)apfill);CHKERRQ(ierr); 1735 ierr = PetscInfo1(ptap->AP_loc,"Use MatPtAP(A,B,MatReuse,%g,&C) for best AP_loc performance.;\n",(double)apfill);CHKERRQ(ierr); 1736 } else { 1737 ierr = PetscInfo(ptap->AP_loc,"Nonscalable algorithm, AP_loc is empty \n");CHKERRQ(ierr); 1738 } 1739 #endif 1740 1741 /* (2-1) compute symbolic Co = Ro*AP_loc */ 1742 /* ------------------------------------ */ 1743 ierr = MatGetOptionsPrefix(A,&prefix);CHKERRQ(ierr); 1744 ierr = MatSetOptionsPrefix(ptap->Ro,prefix);CHKERRQ(ierr); 1745 ierr = MatAppendOptionsPrefix(ptap->Ro,"inner_offdiag_");CHKERRQ(ierr); 1746 ierr = MatMatMultSymbolic_SeqAIJ_SeqAIJ(ptap->Ro,ptap->AP_loc,fill,&ptap->C_oth);CHKERRQ(ierr); 1747 1748 /* (3) send coj of C_oth to other processors */ 1749 /* ------------------------------------------ */ 1750 /* determine row ownership */ 1751 ierr = PetscLayoutCreate(comm,&rowmap);CHKERRQ(ierr); 1752 rowmap->n = pn; 1753 rowmap->bs = 1; 1754 ierr = PetscLayoutSetUp(rowmap);CHKERRQ(ierr); 1755 owners = rowmap->range; 1756 1757 /* determine the number of messages to send, their lengths */ 1758 ierr = PetscMalloc4(size,&len_s,size,&len_si,size,&sstatus,size+2,&owners_co);CHKERRQ(ierr); 1759 ierr = PetscMemzero(len_s,size*sizeof(PetscMPIInt));CHKERRQ(ierr); 1760 ierr = PetscMemzero(len_si,size*sizeof(PetscMPIInt));CHKERRQ(ierr); 1761 1762 c_oth = (Mat_SeqAIJ*)ptap->C_oth->data; 1763 coi = c_oth->i; coj = c_oth->j; 1764 con = ptap->C_oth->rmap->n; 1765 proc = 0; 1766 for (i=0; i<con; i++) { 1767 while (prmap[i] >= owners[proc+1]) proc++; 1768 len_si[proc]++; /* num of rows in Co(=Pt*AP) to be sent to [proc] */ 1769 len_s[proc] += coi[i+1] - coi[i]; /* num of nonzeros in Co to be sent to [proc] */ 1770 } 1771 1772 len = 0; /* max length of buf_si[], see (4) */ 1773 owners_co[0] = 0; 1774 nsend = 0; 1775 for (proc=0; proc<size; proc++) { 1776 owners_co[proc+1] = owners_co[proc] + len_si[proc]; 1777 if (len_s[proc]) { 1778 nsend++; 1779 len_si[proc] = 2*(len_si[proc] + 1); /* length of buf_si to be sent to [proc] */ 1780 len += len_si[proc]; 1781 } 1782 } 1783 1784 /* determine the number and length of messages to receive for coi and coj */ 1785 ierr = PetscGatherNumberOfMessages(comm,NULL,len_s,&nrecv);CHKERRQ(ierr); 1786 ierr = PetscGatherMessageLengths2(comm,nsend,nrecv,len_s,len_si,&id_r,&len_r,&len_ri);CHKERRQ(ierr); 1787 1788 /* post the Irecv and Isend of coj */ 1789 ierr = PetscCommGetNewTag(comm,&tagj);CHKERRQ(ierr); 1790 ierr = PetscPostIrecvInt(comm,tagj,nrecv,id_r,len_r,&buf_rj,&rwaits);CHKERRQ(ierr); 1791 ierr = PetscMalloc1(nsend+1,&swaits);CHKERRQ(ierr); 1792 for (proc=0, k=0; proc<size; proc++) { 1793 if (!len_s[proc]) continue; 1794 i = owners_co[proc]; 1795 ierr = MPI_Isend(coj+coi[i],len_s[proc],MPIU_INT,proc,tagj,comm,swaits+k);CHKERRQ(ierr); 1796 k++; 1797 } 1798 1799 /* (2-2) compute symbolic C_loc = Rd*AP_loc */ 1800 /* ---------------------------------------- */ 1801 ierr = MatSetOptionsPrefix(ptap->Rd,prefix);CHKERRQ(ierr); 1802 ierr = MatAppendOptionsPrefix(ptap->Rd,"inner_diag_");CHKERRQ(ierr); 1803 ierr = MatMatMultSymbolic_SeqAIJ_SeqAIJ(ptap->Rd,ptap->AP_loc,fill,&ptap->C_loc);CHKERRQ(ierr); 1804 c_loc = (Mat_SeqAIJ*)ptap->C_loc->data; 1805 1806 /* receives coj are complete */ 1807 for (i=0; i<nrecv; i++) { 1808 ierr = MPI_Waitany(nrecv,rwaits,&icompleted,&rstatus);CHKERRQ(ierr); 1809 } 1810 ierr = PetscFree(rwaits);CHKERRQ(ierr); 1811 if (nsend) {ierr = MPI_Waitall(nsend,swaits,sstatus);CHKERRQ(ierr);} 1812 1813 /* add received column indices into ta to update Crmax */ 1814 for (k=0; k<nrecv; k++) {/* k-th received message */ 1815 Jptr = buf_rj[k]; 1816 for (j=0; j<len_r[k]; j++) { 1817 ierr = PetscTableAdd(ta,*(Jptr+j)+1,1,INSERT_VALUES);CHKERRQ(ierr); 1818 } 1819 } 1820 ierr = PetscTableGetCount(ta,&Crmax);CHKERRQ(ierr); 1821 ierr = PetscTableDestroy(&ta);CHKERRQ(ierr); 1822 1823 /* (4) send and recv coi */ 1824 /*-----------------------*/ 1825 ierr = PetscCommGetNewTag(comm,&tagi);CHKERRQ(ierr); 1826 ierr = PetscPostIrecvInt(comm,tagi,nrecv,id_r,len_ri,&buf_ri,&rwaits);CHKERRQ(ierr); 1827 ierr = PetscMalloc1(len+1,&buf_s);CHKERRQ(ierr); 1828 buf_si = buf_s; /* points to the beginning of k-th msg to be sent */ 1829 for (proc=0,k=0; proc<size; proc++) { 1830 if (!len_s[proc]) continue; 1831 /* form outgoing message for i-structure: 1832 buf_si[0]: nrows to be sent 1833 [1:nrows]: row index (global) 1834 [nrows+1:2*nrows+1]: i-structure index 1835 */ 1836 /*-------------------------------------------*/ 1837 nrows = len_si[proc]/2 - 1; /* num of rows in Co to be sent to [proc] */ 1838 buf_si_i = buf_si + nrows+1; 1839 buf_si[0] = nrows; 1840 buf_si_i[0] = 0; 1841 nrows = 0; 1842 for (i=owners_co[proc]; i<owners_co[proc+1]; i++) { 1843 nzi = coi[i+1] - coi[i]; 1844 buf_si_i[nrows+1] = buf_si_i[nrows] + nzi; /* i-structure */ 1845 buf_si[nrows+1] = prmap[i] -owners[proc]; /* local row index */ 1846 nrows++; 1847 } 1848 ierr = MPI_Isend(buf_si,len_si[proc],MPIU_INT,proc,tagi,comm,swaits+k);CHKERRQ(ierr); 1849 k++; 1850 buf_si += len_si[proc]; 1851 } 1852 for (i=0; i<nrecv; i++) { 1853 ierr = MPI_Waitany(nrecv,rwaits,&icompleted,&rstatus);CHKERRQ(ierr); 1854 } 1855 ierr = PetscFree(rwaits);CHKERRQ(ierr); 1856 if (nsend) {ierr = MPI_Waitall(nsend,swaits,sstatus);CHKERRQ(ierr);} 1857 1858 ierr = PetscFree4(len_s,len_si,sstatus,owners_co);CHKERRQ(ierr); 1859 ierr = PetscFree(len_ri);CHKERRQ(ierr); 1860 ierr = PetscFree(swaits);CHKERRQ(ierr); 1861 ierr = PetscFree(buf_s);CHKERRQ(ierr); 1862 1863 /* (5) compute the local portion of Cmpi */ 1864 /* ------------------------------------------ */ 1865 /* set initial free space to be Crmax, sufficient for holding nozeros in each row of Cmpi */ 1866 ierr = PetscFreeSpaceGet(Crmax,&free_space);CHKERRQ(ierr); 1867 current_space = free_space; 1868 1869 ierr = PetscMalloc3(nrecv,&buf_ri_k,nrecv,&nextrow,nrecv,&nextci);CHKERRQ(ierr); 1870 for (k=0; k<nrecv; k++) { 1871 buf_ri_k[k] = buf_ri[k]; /* beginning of k-th recved i-structure */ 1872 nrows = *buf_ri_k[k]; 1873 nextrow[k] = buf_ri_k[k] + 1; /* next row number of k-th recved i-structure */ 1874 nextci[k] = buf_ri_k[k] + (nrows + 1); /* poins to the next i-structure of k-th recved i-structure */ 1875 } 1876 1877 ierr = MatPreallocateInitialize(comm,pn,pn,dnz,onz);CHKERRQ(ierr); 1878 ierr = PetscLLCondensedCreate(Crmax,pN,&lnk,&lnkbt);CHKERRQ(ierr); 1879 for (i=0; i<pn; i++) { 1880 /* add C_loc into Cmpi */ 1881 nzi = c_loc->i[i+1] - c_loc->i[i]; 1882 Jptr = c_loc->j + c_loc->i[i]; 1883 ierr = PetscLLCondensedAddSorted(nzi,Jptr,lnk,lnkbt);CHKERRQ(ierr); 1884 1885 /* add received col data into lnk */ 1886 for (k=0; k<nrecv; k++) { /* k-th received message */ 1887 if (i == *nextrow[k]) { /* i-th row */ 1888 nzi = *(nextci[k]+1) - *nextci[k]; 1889 Jptr = buf_rj[k] + *nextci[k]; 1890 ierr = PetscLLCondensedAddSorted(nzi,Jptr,lnk,lnkbt);CHKERRQ(ierr); 1891 nextrow[k]++; nextci[k]++; 1892 } 1893 } 1894 nzi = lnk[0]; 1895 1896 /* copy data into free space, then initialize lnk */ 1897 ierr = PetscLLCondensedClean(pN,nzi,current_space->array,lnk,lnkbt);CHKERRQ(ierr); 1898 ierr = MatPreallocateSet(i+owners[rank],nzi,current_space->array,dnz,onz);CHKERRQ(ierr); 1899 } 1900 ierr = PetscFree3(buf_ri_k,nextrow,nextci);CHKERRQ(ierr); 1901 ierr = PetscLLDestroy(lnk,lnkbt);CHKERRQ(ierr); 1902 ierr = PetscFreeSpaceDestroy(free_space);CHKERRQ(ierr); 1903 1904 /* local sizes and preallocation */ 1905 ierr = MatSetSizes(Cmpi,pn,pn,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 1906 ierr = MatSetBlockSizes(Cmpi,PetscAbs(P->cmap->bs),PetscAbs(P->cmap->bs));CHKERRQ(ierr); 1907 ierr = MatMPIAIJSetPreallocation(Cmpi,0,dnz,0,onz);CHKERRQ(ierr); 1908 ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr); 1909 1910 /* members in merge */ 1911 ierr = PetscFree(id_r);CHKERRQ(ierr); 1912 ierr = PetscFree(len_r);CHKERRQ(ierr); 1913 ierr = PetscFree(buf_ri[0]);CHKERRQ(ierr); 1914 ierr = PetscFree(buf_ri);CHKERRQ(ierr); 1915 ierr = PetscFree(buf_rj[0]);CHKERRQ(ierr); 1916 ierr = PetscFree(buf_rj);CHKERRQ(ierr); 1917 ierr = PetscLayoutDestroy(&rowmap);CHKERRQ(ierr); 1918 1919 /* attach the supporting struct to Cmpi for reuse */ 1920 c = (Mat_MPIAIJ*)Cmpi->data; 1921 c->ap = ptap; 1922 ptap->duplicate = Cmpi->ops->duplicate; 1923 ptap->destroy = Cmpi->ops->destroy; 1924 ptap->view = Cmpi->ops->view; 1925 ierr = PetscCalloc1(pN,&ptap->apa);CHKERRQ(ierr); 1926 1927 /* Cmpi is not ready for use - assembly will be done by MatPtAPNumeric() */ 1928 Cmpi->assembled = PETSC_FALSE; 1929 Cmpi->ops->destroy = MatDestroy_MPIAIJ_PtAP; 1930 Cmpi->ops->view = MatView_MPIAIJ_PtAP; 1931 Cmpi->ops->freeintermediatedatastructures = MatFreeIntermediateDataStructures_MPIAIJ_AP; 1932 *C = Cmpi; 1933 PetscFunctionReturn(0); 1934 } 1935 1936 PetscErrorCode MatPtAPNumeric_MPIAIJ_MPIAIJ(Mat A,Mat P,Mat C) 1937 { 1938 PetscErrorCode ierr; 1939 Mat_MPIAIJ *a=(Mat_MPIAIJ*)A->data,*p=(Mat_MPIAIJ*)P->data,*c=(Mat_MPIAIJ*)C->data; 1940 Mat_SeqAIJ *ad=(Mat_SeqAIJ*)(a->A)->data,*ao=(Mat_SeqAIJ*)(a->B)->data; 1941 Mat_SeqAIJ *ap,*p_loc,*p_oth=NULL,*c_seq; 1942 Mat_APMPI *ptap = c->ap; 1943 Mat AP_loc,C_loc,C_oth; 1944 PetscInt i,rstart,rend,cm,ncols,row; 1945 PetscInt *api,*apj,am = A->rmap->n,j,col,apnz; 1946 PetscScalar *apa; 1947 const PetscInt *cols; 1948 const PetscScalar *vals; 1949 1950 PetscFunctionBegin; 1951 if (!ptap) { 1952 MPI_Comm comm; 1953 ierr = PetscObjectGetComm((PetscObject)C,&comm);CHKERRQ(ierr); 1954 SETERRQ(comm,PETSC_ERR_ARG_WRONGSTATE,"PtAP cannot be reused. Do not call MatFreeIntermediateDataStructures() or use '-mat_freeintermediatedatastructures'"); 1955 } 1956 1957 ierr = MatZeroEntries(C);CHKERRQ(ierr); 1958 /* 1) get R = Pd^T,Ro = Po^T */ 1959 if (ptap->reuse == MAT_REUSE_MATRIX) { 1960 ierr = MatTranspose(p->A,MAT_REUSE_MATRIX,&ptap->Rd);CHKERRQ(ierr); 1961 ierr = MatTranspose(p->B,MAT_REUSE_MATRIX,&ptap->Ro);CHKERRQ(ierr); 1962 } 1963 1964 /* 2) get AP_loc */ 1965 AP_loc = ptap->AP_loc; 1966 ap = (Mat_SeqAIJ*)AP_loc->data; 1967 1968 /* 2-1) get P_oth = ptap->P_oth and P_loc = ptap->P_loc */ 1969 /*-----------------------------------------------------*/ 1970 if (ptap->reuse == MAT_REUSE_MATRIX) { 1971 /* P_oth and P_loc are obtained in MatPtASymbolic() when reuse == MAT_INITIAL_MATRIX */ 1972 ierr = MatGetBrowsOfAoCols_MPIAIJ(A,P,MAT_REUSE_MATRIX,&ptap->startsj_s,&ptap->startsj_r,&ptap->bufa,&ptap->P_oth);CHKERRQ(ierr); 1973 ierr = MatMPIAIJGetLocalMat(P,MAT_REUSE_MATRIX,&ptap->P_loc);CHKERRQ(ierr); 1974 } 1975 1976 /* 2-2) compute numeric A_loc*P - dominating part */ 1977 /* ---------------------------------------------- */ 1978 /* get data from symbolic products */ 1979 p_loc = (Mat_SeqAIJ*)(ptap->P_loc)->data; 1980 if (ptap->P_oth) { 1981 p_oth = (Mat_SeqAIJ*)(ptap->P_oth)->data; 1982 } 1983 apa = ptap->apa; 1984 api = ap->i; 1985 apj = ap->j; 1986 for (i=0; i<am; i++) { 1987 /* AP[i,:] = A[i,:]*P = Ad*P_loc Ao*P_oth */ 1988 AProw_nonscalable(i,ad,ao,p_loc,p_oth,apa); 1989 apnz = api[i+1] - api[i]; 1990 for (j=0; j<apnz; j++) { 1991 col = apj[j+api[i]]; 1992 ap->a[j+ap->i[i]] = apa[col]; 1993 apa[col] = 0.0; 1994 } 1995 ierr = PetscLogFlops(2.0*apnz);CHKERRQ(ierr); 1996 } 1997 1998 /* 3) C_loc = Rd*AP_loc, C_oth = Ro*AP_loc */ 1999 ierr = ((ptap->C_loc)->ops->matmultnumeric)(ptap->Rd,AP_loc,ptap->C_loc);CHKERRQ(ierr); 2000 ierr = ((ptap->C_oth)->ops->matmultnumeric)(ptap->Ro,AP_loc,ptap->C_oth);CHKERRQ(ierr); 2001 C_loc = ptap->C_loc; 2002 C_oth = ptap->C_oth; 2003 2004 /* add C_loc and Co to to C */ 2005 ierr = MatGetOwnershipRange(C,&rstart,&rend);CHKERRQ(ierr); 2006 2007 /* C_loc -> C */ 2008 cm = C_loc->rmap->N; 2009 c_seq = (Mat_SeqAIJ*)C_loc->data; 2010 cols = c_seq->j; 2011 vals = c_seq->a; 2012 2013 2014 /* The (fast) MatSetValues_MPIAIJ_CopyFromCSRFormat function can only be used when C->was_assembled is PETSC_FALSE and */ 2015 /* when there are no off-processor parts. */ 2016 /* If was_assembled is true, then the statement aj[rowstart_diag+dnz_row] = mat_j[col] - cstart; in MatSetValues_MPIAIJ_CopyFromCSRFormat */ 2017 /* is no longer true. Then the more complex function MatSetValues_MPIAIJ() has to be used, where the column index is looked up from */ 2018 /* a table, and other, more complex stuff has to be done. */ 2019 if (C->assembled) { 2020 C->was_assembled = PETSC_TRUE; 2021 C->assembled = PETSC_FALSE; 2022 } 2023 if (C->was_assembled) { 2024 for (i=0; i<cm; i++) { 2025 ncols = c_seq->i[i+1] - c_seq->i[i]; 2026 row = rstart + i; 2027 ierr = MatSetValues_MPIAIJ(C,1,&row,ncols,cols,vals,ADD_VALUES);CHKERRQ(ierr); 2028 cols += ncols; vals += ncols; 2029 } 2030 } else { 2031 ierr = MatSetValues_MPIAIJ_CopyFromCSRFormat(C,c_seq->j,c_seq->i,c_seq->a);CHKERRQ(ierr); 2032 } 2033 2034 /* Co -> C, off-processor part */ 2035 cm = C_oth->rmap->N; 2036 c_seq = (Mat_SeqAIJ*)C_oth->data; 2037 cols = c_seq->j; 2038 vals = c_seq->a; 2039 for (i=0; i<cm; i++) { 2040 ncols = c_seq->i[i+1] - c_seq->i[i]; 2041 row = p->garray[i]; 2042 ierr = MatSetValues(C,1,&row,ncols,cols,vals,ADD_VALUES);CHKERRQ(ierr); 2043 cols += ncols; vals += ncols; 2044 } 2045 2046 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2047 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2048 2049 ptap->reuse = MAT_REUSE_MATRIX; 2050 2051 /* supporting struct ptap consumes almost same amount of memory as C=PtAP, release it if C will not be updated by A and P */ 2052 if (ptap->freestruct) { 2053 ierr = MatFreeIntermediateDataStructures(C);CHKERRQ(ierr); 2054 } 2055 PetscFunctionReturn(0); 2056 } 2057