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