1 2 /* 3 Defines projective product routines where A is a SeqAIJ 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 <petscbt.h> 10 11 #undef __FUNCT__ 12 #define __FUNCT__ "MatPtAPSymbolic_SeqAIJ" 13 PetscErrorCode MatPtAPSymbolic_SeqAIJ(Mat A,Mat P,PetscReal fill,Mat *C) 14 { 15 PetscErrorCode ierr; 16 17 PetscFunctionBegin; 18 if (!P->ops->ptapsymbolic_seqaij) SETERRQ2(((PetscObject)A)->comm,PETSC_ERR_SUP,"Not implemented for A=%s and P=%s",((PetscObject)A)->type_name,((PetscObject)P)->type_name); 19 ierr = (*P->ops->ptapsymbolic_seqaij)(A,P,fill,C);CHKERRQ(ierr); 20 PetscFunctionReturn(0); 21 } 22 23 #undef __FUNCT__ 24 #define __FUNCT__ "MatPtAPNumeric_SeqAIJ" 25 PetscErrorCode MatPtAPNumeric_SeqAIJ(Mat A,Mat P,Mat C) 26 { 27 PetscErrorCode ierr; 28 29 PetscFunctionBegin; 30 if (!P->ops->ptapnumeric_seqaij) SETERRQ2(((PetscObject)A)->comm,PETSC_ERR_SUP,"Not implemented for A=%s and P=%s",((PetscObject)A)->type_name,((PetscObject)P)->type_name); 31 ierr = (*P->ops->ptapnumeric_seqaij)(A,P,C);CHKERRQ(ierr); 32 PetscFunctionReturn(0); 33 } 34 35 #undef __FUNCT__ 36 #define __FUNCT__ "MatDestroy_SeqAIJ_PtAP" 37 PetscErrorCode MatDestroy_SeqAIJ_PtAP(Mat A) 38 { 39 PetscErrorCode ierr; 40 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 41 Mat_PtAP *ptap = a->ptap; 42 43 PetscFunctionBegin; 44 /* free ptap, then A */ 45 ierr = PetscFree(ptap->apa);CHKERRQ(ierr); 46 ierr = PetscFree(ptap->api);CHKERRQ(ierr); 47 ierr = PetscFree(ptap->apj);CHKERRQ(ierr); 48 ierr = (ptap->destroy)(A);CHKERRQ(ierr); 49 ierr = PetscFree(ptap);CHKERRQ(ierr); 50 PetscFunctionReturn(0); 51 } 52 53 #undef __FUNCT__ 54 #define __FUNCT__ "MatPtAPSymbolic_SeqAIJ_SeqAIJ_SparseAxpy2" 55 PetscErrorCode MatPtAPSymbolic_SeqAIJ_SeqAIJ_SparseAxpy2(Mat A,Mat P,PetscReal fill,Mat *C) 56 { 57 PetscErrorCode ierr; 58 PetscFreeSpaceList free_space=PETSC_NULL,current_space=PETSC_NULL; 59 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*p = (Mat_SeqAIJ*)P->data,*c; 60 PetscInt *pti,*ptj,*ptJ,*ai=a->i,*aj=a->j,*ajj,*pi=p->i,*pj=p->j,*pjj; 61 PetscInt *ci,*cj,*ptadenserow,*ptasparserow,*ptaj,nspacedouble=0; 62 PetscInt an=A->cmap->N,am=A->rmap->N,pn=P->cmap->N; 63 PetscInt i,j,k,ptnzi,arow,anzj,ptanzi,prow,pnzj,cnzi,nlnk,*lnk; 64 MatScalar *ca; 65 PetscBT lnkbt; 66 67 PetscFunctionBegin; 68 /* Get ij structure of P^T */ 69 ierr = MatGetSymbolicTranspose_SeqAIJ(P,&pti,&ptj);CHKERRQ(ierr); 70 ptJ=ptj; 71 72 /* Allocate ci array, arrays for fill computation and */ 73 /* free space for accumulating nonzero column info */ 74 ierr = PetscMalloc((pn+1)*sizeof(PetscInt),&ci);CHKERRQ(ierr); 75 ci[0] = 0; 76 77 ierr = PetscMalloc((2*an+1)*sizeof(PetscInt),&ptadenserow);CHKERRQ(ierr); 78 ierr = PetscMemzero(ptadenserow,(2*an+1)*sizeof(PetscInt));CHKERRQ(ierr); 79 ptasparserow = ptadenserow + an; 80 81 /* create and initialize a linked list */ 82 nlnk = pn+1; 83 ierr = PetscLLCreate(pn,pn,nlnk,lnk,lnkbt);CHKERRQ(ierr); 84 85 /* Set initial free space to be fill*nnz(A). */ 86 /* This should be reasonable if sparsity of PtAP is similar to that of A. */ 87 ierr = PetscFreeSpaceGet((PetscInt)(fill*ai[am]),&free_space); 88 current_space = free_space; 89 90 /* Determine symbolic info for each row of C: */ 91 for (i=0;i<pn;i++) { 92 ptnzi = pti[i+1] - pti[i]; 93 ptanzi = 0; 94 /* Determine symbolic row of PtA: */ 95 for (j=0;j<ptnzi;j++) { 96 arow = *ptJ++; 97 anzj = ai[arow+1] - ai[arow]; 98 ajj = aj + ai[arow]; 99 for (k=0;k<anzj;k++) { 100 if (!ptadenserow[ajj[k]]) { 101 ptadenserow[ajj[k]] = -1; 102 ptasparserow[ptanzi++] = ajj[k]; 103 } 104 } 105 } 106 /* Using symbolic info for row of PtA, determine symbolic info for row of C: */ 107 ptaj = ptasparserow; 108 cnzi = 0; 109 for (j=0;j<ptanzi;j++) { 110 prow = *ptaj++; 111 pnzj = pi[prow+1] - pi[prow]; 112 pjj = pj + pi[prow]; 113 /* add non-zero cols of P into the sorted linked list lnk */ 114 ierr = PetscLLAdd(pnzj,pjj,pn,nlnk,lnk,lnkbt);CHKERRQ(ierr); 115 cnzi += nlnk; 116 } 117 118 /* If free space is not available, make more free space */ 119 /* Double the amount of total space in the list */ 120 if (current_space->local_remaining<cnzi) { 121 ierr = PetscFreeSpaceGet(cnzi+current_space->total_array_size,¤t_space);CHKERRQ(ierr); 122 nspacedouble++; 123 } 124 125 /* Copy data into free space, and zero out denserows */ 126 ierr = PetscLLClean(pn,pn,cnzi,lnk,current_space->array,lnkbt);CHKERRQ(ierr); 127 current_space->array += cnzi; 128 current_space->local_used += cnzi; 129 current_space->local_remaining -= cnzi; 130 131 for (j=0;j<ptanzi;j++) { 132 ptadenserow[ptasparserow[j]] = 0; 133 } 134 /* Aside: Perhaps we should save the pta info for the numerical factorization. */ 135 /* For now, we will recompute what is needed. */ 136 ci[i+1] = ci[i] + cnzi; 137 } 138 /* nnz is now stored in ci[ptm], column indices are in the list of free space */ 139 /* Allocate space for cj, initialize cj, and */ 140 /* destroy list of free space and other temporary array(s) */ 141 ierr = PetscMalloc((ci[pn]+1)*sizeof(PetscInt),&cj);CHKERRQ(ierr); 142 ierr = PetscFreeSpaceContiguous(&free_space,cj);CHKERRQ(ierr); 143 ierr = PetscFree(ptadenserow);CHKERRQ(ierr); 144 ierr = PetscLLDestroy(lnk,lnkbt);CHKERRQ(ierr); 145 146 /* Allocate space for ca */ 147 ierr = PetscMalloc((ci[pn]+1)*sizeof(MatScalar),&ca);CHKERRQ(ierr); 148 ierr = PetscMemzero(ca,(ci[pn]+1)*sizeof(MatScalar));CHKERRQ(ierr); 149 150 /* put together the new matrix */ 151 ierr = MatCreateSeqAIJWithArrays(((PetscObject)A)->comm,pn,pn,ci,cj,ca,C);CHKERRQ(ierr); 152 153 /* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */ 154 /* Since these are PETSc arrays, change flags to free them as necessary. */ 155 c = (Mat_SeqAIJ *)((*C)->data); 156 c->free_a = PETSC_TRUE; 157 c->free_ij = PETSC_TRUE; 158 c->nonew = 0; 159 A->ops->ptapnumeric = MatPtAPNumeric_SeqAIJ_SeqAIJ_SparseAxpy2; /* should use *C->ops until PtAP insterface is updated to double dispatch as MatMatMult() */ 160 161 /* Clean up. */ 162 ierr = MatRestoreSymbolicTranspose_SeqAIJ(P,&pti,&ptj);CHKERRQ(ierr); 163 #if defined(PETSC_USE_INFO) 164 if (ci[pn] != 0) { 165 PetscReal afill = ((PetscReal)ci[pn])/ai[am]; 166 if (afill < 1.0) afill = 1.0; 167 ierr = PetscInfo3((*C),"Reallocs %D; Fill ratio: given %G needed %G.\n",nspacedouble,fill,afill);CHKERRQ(ierr); 168 ierr = PetscInfo1((*C),"Use MatPtAP(A,P,MatReuse,%G,&C) for best performance.\n",afill);CHKERRQ(ierr); 169 } else { 170 ierr = PetscInfo((*C),"Empty matrix product\n");CHKERRQ(ierr); 171 } 172 #endif 173 PetscFunctionReturn(0); 174 } 175 176 #undef __FUNCT__ 177 #define __FUNCT__ "MatPtAPNumeric_SeqAIJ_SeqAIJ_SparseAxpy2" 178 PetscErrorCode MatPtAPNumeric_SeqAIJ_SeqAIJ_SparseAxpy2(Mat A,Mat P,Mat C) 179 { 180 PetscErrorCode ierr; 181 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 182 Mat_SeqAIJ *p = (Mat_SeqAIJ *) P->data; 183 Mat_SeqAIJ *c = (Mat_SeqAIJ *) C->data; 184 PetscInt *ai=a->i,*aj=a->j,*apj,*apjdense,*pi=p->i,*pj=p->j,*pJ=p->j,*pjj; 185 PetscInt *ci=c->i,*cj=c->j,*cjj; 186 PetscInt am=A->rmap->N,cn=C->cmap->N,cm=C->rmap->N; 187 PetscInt i,j,k,anzi,pnzi,apnzj,nextap,pnzj,prow,crow; 188 MatScalar *aa=a->a,*apa,*pa=p->a,*pA=p->a,*paj,*ca=c->a,*caj; 189 190 PetscFunctionBegin; 191 /* Allocate temporary array for storage of one row of A*P */ 192 ierr = PetscMalloc(cn*(sizeof(MatScalar)+sizeof(PetscInt))+c->rmax*sizeof(PetscInt),&apa);CHKERRQ(ierr); 193 194 apjdense = (PetscInt *)(apa + cn); 195 apj = apjdense + cn; 196 ierr = PetscMemzero(apa,cn*(sizeof(MatScalar)+sizeof(PetscInt)));CHKERRQ(ierr); 197 198 /* Clear old values in C */ 199 ierr = PetscMemzero(ca,ci[cm]*sizeof(MatScalar));CHKERRQ(ierr); 200 201 for (i=0; i<am; i++) { 202 /* Form sparse row of A*P */ 203 anzi = ai[i+1] - ai[i]; 204 apnzj = 0; 205 for (j=0; j<anzi; j++) { 206 prow = *aj++; 207 pnzj = pi[prow+1] - pi[prow]; 208 pjj = pj + pi[prow]; 209 paj = pa + pi[prow]; 210 for (k=0;k<pnzj;k++) { 211 if (!apjdense[pjj[k]]) { 212 apjdense[pjj[k]] = -1; 213 apj[apnzj++] = pjj[k]; 214 } 215 apa[pjj[k]] += (*aa)*paj[k]; 216 } 217 ierr = PetscLogFlops(2.0*pnzj);CHKERRQ(ierr); 218 aa++; 219 } 220 221 /* Sort the j index array for quick sparse axpy. */ 222 /* Note: a array does not need sorting as it is in dense storage locations. */ 223 ierr = PetscSortInt(apnzj,apj);CHKERRQ(ierr); 224 225 /* Compute P^T*A*P using outer product (P^T)[:,j]*(A*P)[j,:]. */ 226 pnzi = pi[i+1] - pi[i]; 227 for (j=0; j<pnzi; j++) { 228 nextap = 0; 229 crow = *pJ++; 230 cjj = cj + ci[crow]; 231 caj = ca + ci[crow]; 232 /* Perform sparse axpy operation. Note cjj includes apj. */ 233 for (k=0;nextap<apnzj;k++) { 234 #if defined(PETSC_USE_DEBUG) 235 if (k >= ci[crow+1] - ci[crow]) { 236 SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"k too large k %d, crow %d",k,crow); 237 } 238 #endif 239 if (cjj[k]==apj[nextap]) { 240 caj[k] += (*pA)*apa[apj[nextap++]]; 241 } 242 } 243 ierr = PetscLogFlops(2.0*apnzj);CHKERRQ(ierr); 244 pA++; 245 } 246 247 /* Zero the current row info for A*P */ 248 for (j=0; j<apnzj; j++) { 249 apa[apj[j]] = 0.; 250 apjdense[apj[j]] = 0; 251 } 252 } 253 254 /* Assemble the final matrix and clean up */ 255 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 256 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 257 ierr = PetscFree(apa);CHKERRQ(ierr); 258 PetscFunctionReturn(0); 259 } 260 261 #undef __FUNCT__ 262 #define __FUNCT__ "MatPtAPSymbolic_SeqAIJ_SeqAIJ" 263 PetscErrorCode MatPtAPSymbolic_SeqAIJ_SeqAIJ(Mat A,Mat P,PetscReal fill,Mat *C) 264 { 265 PetscErrorCode ierr; 266 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*p = (Mat_SeqAIJ*)P->data,*c; 267 PetscInt *pti,*ptj,*ptJ,*ai=a->i,*aj=a->j,*pi=p->i,*pj=p->j,*api,*apj; 268 PetscInt *ci,*cj,ndouble_ap,ndouble_ptap; 269 PetscInt an=A->cmap->N,am=A->rmap->N,pn=P->cmap->N; 270 MatScalar *ca; 271 Mat_PtAP *ptap; 272 PetscInt sparse_axpy=0; 273 274 PetscFunctionBegin; 275 /* flag 'sparse_axpy' determines which implementations to be used: 276 0: do dense axpy in MatPtAPNumeric() - fastest, but requires storage of struct A*P; (default) 277 1: do one sparse axpy - uses same memory as sparse_axpy=0 and might execute less flops 278 (apnz vs. cnz in the outerproduct), slower than case '0' when cnz is not too large than apnz; 279 2: do two sparse axpy in MatPtAPNumeric() - slowest, does not store structure of A*P. */ 280 ierr = PetscOptionsGetInt(PETSC_NULL,"-matptap_sparseaxpy",&sparse_axpy,PETSC_NULL);CHKERRQ(ierr); 281 if (sparse_axpy == 2){ 282 ierr = MatPtAPSymbolic_SeqAIJ_SeqAIJ_SparseAxpy2(A,P,fill,C);CHKERRQ(ierr); 283 PetscFunctionReturn(0); 284 } 285 286 /* Get ij structure of Pt = P^T */ 287 ierr = MatGetSymbolicTranspose_SeqAIJ(P,&pti,&ptj);CHKERRQ(ierr); 288 ptJ=ptj; 289 290 /* Get structure of AP = A*P */ 291 ierr = MatGetSymbolicMatMatMult_SeqAIJ_SeqAIJ(am,ai,aj,an,pn,pi,pj,fill,&api,&apj,&ndouble_ap);CHKERRQ(ierr); 292 293 /* Get structure of C = Pt*AP */ 294 ierr = MatGetSymbolicMatMatMult_SeqAIJ_SeqAIJ(pn,pti,ptj,am,pn,api,apj,fill,&ci,&cj,&ndouble_ptap);CHKERRQ(ierr); 295 296 /* Allocate space for ca */ 297 ierr = PetscMalloc((ci[pn]+1)*sizeof(MatScalar),&ca);CHKERRQ(ierr); 298 ierr = PetscMemzero(ca,(ci[pn]+1)*sizeof(MatScalar));CHKERRQ(ierr); 299 300 /* put together the new matrix */ 301 ierr = MatCreateSeqAIJWithArrays(((PetscObject)A)->comm,pn,pn,ci,cj,ca,C);CHKERRQ(ierr); 302 303 /* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */ 304 /* Since these are PETSc arrays, change flags to free them as necessary. */ 305 c = (Mat_SeqAIJ *)(*C)->data; 306 c->free_a = PETSC_TRUE; 307 c->free_ij = PETSC_TRUE; 308 c->nonew = 0; 309 310 /* create a supporting struct for reuse by MatPtAPNumeric() */ 311 ierr = PetscNew(Mat_PtAP,&ptap);CHKERRQ(ierr); 312 c->ptap = ptap; 313 ptap->destroy = (*C)->ops->destroy; 314 (*C)->ops->destroy = MatDestroy_SeqAIJ_PtAP; 315 316 /* Allocate temporary array for storage of one row of A*P */ 317 ierr = PetscMalloc((pn+1)*sizeof(PetscScalar),&ptap->apa);CHKERRQ(ierr); 318 ierr = PetscMemzero(ptap->apa,(pn+1)*sizeof(PetscScalar));CHKERRQ(ierr); 319 if (sparse_axpy == 1){ 320 A->ops->ptapnumeric = MatPtAPNumeric_SeqAIJ_SeqAIJ_SparseAxpy; 321 } else { 322 A->ops->ptapnumeric = MatPtAPNumeric_SeqAIJ_SeqAIJ; 323 } 324 ptap->api = api; 325 ptap->apj = apj; 326 327 /* Clean up. */ 328 ierr = MatRestoreSymbolicTranspose_SeqAIJ(P,&pti,&ptj);CHKERRQ(ierr); 329 #if defined(PETSC_USE_INFO) 330 if (ci[pn] != 0) { 331 PetscReal apfill,ptapfill; 332 apfill = ((PetscReal)api[am])/(ai[am]+pi[an]); 333 if (apfill < 1.0) apfill = 1.0; 334 ierr = PetscInfo3((*C),"A*P: Reallocs %D; Fill ratio: given %G needed %G.\n",ndouble_ap,fill,apfill);CHKERRQ(ierr); 335 ptapfill = ((PetscReal)ci[pn])/(pi[an]+api[am]); 336 if (ptapfill < 1.0) ptapfill = 1.0; 337 ierr = PetscInfo3((*C),"Pt*AP: Reallocs %D; Fill ratio: given %G needed %G.\n",ndouble_ptap,fill,ptapfill);CHKERRQ(ierr); 338 339 ierr = PetscInfo1((*C),"Use MatPtAP(A,P,MatReuse,%G,&C) for best performance.\n",PetscMax(apfill,ptapfill));CHKERRQ(ierr); 340 ierr = PetscInfo4((*C),"nonzeros: A %D, P %D, A*P %D, C=PtAP %D\n",ai[am],pi[an],api[am],ci[pn]);CHKERRQ(ierr); 341 } else { 342 ierr = PetscInfo((*C),"Empty matrix product\n");CHKERRQ(ierr); 343 } 344 #endif 345 PetscFunctionReturn(0); 346 } 347 348 /* #define PROFILE_MatPtAPNumeric */ 349 #undef __FUNCT__ 350 #define __FUNCT__ "MatPtAPNumeric_SeqAIJ_SeqAIJ" 351 PetscErrorCode MatPtAPNumeric_SeqAIJ_SeqAIJ(Mat A,Mat P,Mat C) 352 { 353 PetscErrorCode ierr; 354 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 355 Mat_SeqAIJ *p = (Mat_SeqAIJ *) P->data; 356 Mat_SeqAIJ *c = (Mat_SeqAIJ *) C->data; 357 PetscInt *ai=a->i,*aj=a->j,*pi=p->i,*pj=p->j,*ci=c->i,*cj=c->j; 358 PetscScalar *aa=a->a,*pa=p->a; 359 PetscInt *apj,*pcol,*cjj,cnz; 360 PetscInt am=A->rmap->N,cm=C->rmap->N; 361 PetscInt i,j,k,anz,apnz,pnz,prow,crow; 362 PetscScalar *apa,*pval,*ca=c->a,*caj,pvalj; 363 Mat_PtAP *ptap = c->ptap; 364 #if defined(PROFILE_MatPtAPNumeric) 365 PetscLogDouble t0,tf,time_Cseq0=0.0,time_Cseq1=0.0; 366 PetscInt flops0=0,flops1=0; 367 #endif 368 369 PetscFunctionBegin; 370 /* Get temporary array for storage of one row of A*P */ 371 apa = ptap->apa; 372 373 /* Clear old values in C */ 374 ierr = PetscMemzero(ca,ci[cm]*sizeof(MatScalar));CHKERRQ(ierr); 375 376 for (i=0;i<am;i++) { 377 /* Form sparse row of AP[i,:] = A[i,:]*P */ 378 #if defined(PROFILE_MatPtAPNumeric) 379 ierr = PetscGetTime(&t0);CHKERRQ(ierr); 380 #endif 381 anz = ai[i+1] - ai[i]; 382 apnz = 0; 383 for (j=0; j<anz; j++) { 384 prow = aj[j]; 385 pnz = pi[prow+1] - pi[prow]; 386 pcol = pj + pi[prow]; 387 pval = pa + pi[prow]; 388 for (k=0; k<pnz; k++) { 389 apa[pcol[k]] += aa[j]*pval[k]; 390 } 391 ierr = PetscLogFlops(2.0*pnz);CHKERRQ(ierr); 392 #if defined(PROFILE_MatPtAPNumeric) 393 flops0 += 2.0*pnz; 394 #endif 395 } 396 aj += anz; aa += anz; 397 #if defined(PROFILE_MatPtAPNumeric) 398 ierr = PetscGetTime(&tf);CHKERRQ(ierr); 399 time_Cseq0 += tf - t0; 400 #endif 401 402 /* Compute P^T*A*P using outer product P[i,:]^T*AP[i,:]. */ 403 #if defined(PROFILE_MatPtAPNumeric) 404 ierr = PetscGetTime(&t0);CHKERRQ(ierr); 405 #endif 406 apj = ptap->apj + ptap->api[i]; 407 apnz = ptap->api[i+1] - ptap->api[i]; 408 pnz = pi[i+1] - pi[i]; 409 pcol = pj + pi[i]; 410 pval = pa + pi[i]; 411 412 /* Perform dense axpy */ 413 for (j=0; j<pnz; j++) { 414 crow = pcol[j]; 415 cjj = cj + ci[crow]; 416 caj = ca + ci[crow]; 417 pvalj = pval[j]; 418 cnz = ci[crow+1] - ci[crow]; 419 for (k=0; k<cnz; k++){ 420 caj[k] += pvalj*apa[cjj[k]]; 421 } 422 ierr = PetscLogFlops(2.0*cnz);CHKERRQ(ierr); 423 #if defined(PROFILE_MatPtAPNumeric) 424 flops1 += 2.0*cnz; 425 #endif 426 } 427 #if defined(PROFILE_MatPtAPNumeric) 428 ierr = PetscGetTime(&tf);CHKERRQ(ierr); 429 time_Cseq1 += tf - t0; 430 #endif 431 432 /* Zero the current row info for A*P */ 433 for (j=0; j<apnz; j++) apa[apj[j]] = 0.0; 434 } 435 436 /* Assemble the final matrix and clean up */ 437 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 438 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 439 #if defined(PROFILE_MatPtAPNumeric) 440 printf("PtAPNumeric_SeqAIJ time %g + %g, flops %d %d\n",time_Cseq0,time_Cseq1,flops0,flops1); 441 #endif 442 PetscFunctionReturn(0); 443 } 444 445 #undef __FUNCT__ 446 #define __FUNCT__ "MatPtAPNumeric_SeqAIJ_SeqAIJ_SparseAxpy" 447 PetscErrorCode MatPtAPNumeric_SeqAIJ_SeqAIJ_SparseAxpy(Mat A,Mat P,Mat C) 448 { 449 PetscErrorCode ierr; 450 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 451 Mat_SeqAIJ *p = (Mat_SeqAIJ *) P->data; 452 Mat_SeqAIJ *c = (Mat_SeqAIJ *) C->data; 453 PetscInt *ai=a->i,*aj=a->j,*pi=p->i,*pj=p->j,*ci=c->i,*cj=c->j; 454 PetscScalar *aa=a->a,*pa=p->a; 455 PetscInt *apj,*pcol,*cjj; 456 PetscInt am=A->rmap->N,cm=C->rmap->N; 457 PetscInt i,j,k,anz,apnz,pnz,prow,crow,apcol,nextap; 458 PetscScalar *apa,*pval,*ca=c->a,*caj,pvalj; 459 Mat_PtAP *ptap = c->ptap; 460 #if defined(PROFILE_MatPtAPNumeric) 461 PetscLogDouble t0,tf,time_Cseq0=0.0,time_Cseq1=0.0; 462 PetscInt flops0=0,flops1=0; 463 #endif 464 465 PetscFunctionBegin; 466 /* Get temporary array for storage of one row of A*P */ 467 apa = ptap->apa; 468 469 /* Clear old values in C */ 470 ierr = PetscMemzero(ca,ci[cm]*sizeof(MatScalar));CHKERRQ(ierr); 471 472 for (i=0;i<am;i++) { 473 /* Form sparse row of AP[i,:] = A[i,:]*P */ 474 #if defined(PROFILE_MatPtAPNumeric) 475 ierr = PetscGetTime(&t0);CHKERRQ(ierr); 476 #endif 477 anz = ai[i+1] - ai[i]; 478 apnz = 0; 479 for (j=0; j<anz; j++) { 480 prow = aj[j]; 481 pnz = pi[prow+1] - pi[prow]; 482 pcol = pj + pi[prow]; 483 pval = pa + pi[prow]; 484 for (k=0; k<pnz; k++) { 485 apa[pcol[k]] += aa[j]*pval[k]; 486 } 487 ierr = PetscLogFlops(2.0*pnz);CHKERRQ(ierr); 488 #if defined(PROFILE_MatPtAPNumeric) 489 flops0 += 2.0*pnz; 490 #endif 491 } 492 aj += anz; aa += anz; 493 #if defined(PROFILE_MatPtAPNumeric) 494 ierr = PetscGetTime(&tf);CHKERRQ(ierr); 495 time_Cseq0 += tf - t0; 496 #endif 497 498 /* Compute P^T*A*P using outer product P[i,:]^T*AP[i,:]. */ 499 #if defined(PROFILE_MatPtAPNumeric) 500 ierr = PetscGetTime(&t0);CHKERRQ(ierr); 501 #endif 502 apj = ptap->apj + ptap->api[i]; 503 apnz = ptap->api[i+1] - ptap->api[i]; 504 pnz = pi[i+1] - pi[i]; 505 pcol = pj + pi[i]; 506 pval = pa + pi[i]; 507 508 /* Perform sparse axpy */ 509 for (j=0; j<pnz; j++) { 510 crow = pcol[j]; 511 cjj = cj + ci[crow]; 512 caj = ca + ci[crow]; 513 pvalj = pval[j]; 514 nextap = 0; 515 apcol = apj[nextap]; 516 for (k=0; nextap<apnz; k++) { 517 if (cjj[k] == apcol) { 518 caj[k] += pvalj*apa[apcol]; 519 apcol = apj[++nextap]; 520 } 521 } 522 ierr = PetscLogFlops(2.0*apnz);CHKERRQ(ierr); 523 #if defined(PROFILE_MatPtAPNumeric) 524 flops1 += 2.0*apnz; 525 #endif 526 } 527 #if defined(PROFILE_MatPtAPNumeric) 528 ierr = PetscGetTime(&tf);CHKERRQ(ierr); 529 time_Cseq1 += tf - t0; 530 #endif 531 532 /* Zero the current row info for A*P */ 533 for (j=0; j<apnz; j++) { 534 apcol = apj[j]; 535 apa[apcol] = 0.; 536 } 537 } 538 539 /* Assemble the final matrix and clean up */ 540 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 541 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 542 #if defined(PROFILE_MatPtAPNumeric) 543 printf("MatPtAPNumeric_SeqAIJ_SeqAIJ_SparseAxpy time %g + %g, flops %d %d\n",time_Cseq0,time_Cseq1,flops0,flops1); 544 #endif 545 PetscFunctionReturn(0); 546 } 547