1 /* 2 Defines projective product routines where A is a SeqAIJ matrix 3 C = P^T * A * P 4 */ 5 6 #include "src/mat/impls/aij/seq/aij.h" 7 #include "src/mat/utils/freespace.h" 8 9 int MatSeqAIJPtAP(Mat,Mat,Mat*); 10 int MatSeqAIJPtAPSymbolic(Mat,Mat,Mat*); 11 int MatSeqAIJPtAPNumeric(Mat,Mat,Mat); 12 13 static int MATSeqAIJ_PtAP = 0; 14 static int MATSeqAIJ_PtAPSymbolic = 0; 15 static int MATSeqAIJ_PtAPNumeric = 0; 16 17 /* 18 MatSeqAIJPtAP - Creates the SeqAIJ matrix product, C, 19 of SeqAIJ matrix A and matrix P: 20 C = P^T * A * P; 21 22 Note: C is assumed to be uncreated. 23 If this is not the case, Destroy C before calling this routine. 24 */ 25 #undef __FUNCT__ 26 #define __FUNCT__ "MatSeqAIJPtAP" 27 int MatSeqAIJPtAP(Mat A,Mat P,Mat *C) { 28 int ierr; 29 char funct[80]; 30 31 PetscFunctionBegin; 32 33 ierr = PetscLogEventBegin(MATSeqAIJ_PtAP,A,P,0,0);CHKERRQ(ierr); 34 35 ierr = MatSeqAIJPtAPSymbolic(A,P,C);CHKERRQ(ierr); 36 37 /* Avoid additional error checking included in */ 38 /* ierr = MatSeqAIJApplyPtAPNumeric(A,P,*C);CHKERRQ(ierr); */ 39 40 /* Query A for ApplyPtAPNumeric implementation based on types of P */ 41 ierr = PetscStrcpy(funct,"MatApplyPtAPNumeric_seqaij_");CHKERRQ(ierr); 42 ierr = PetscStrcat(funct,P->type_name);CHKERRQ(ierr); 43 ierr = PetscTryMethod(A,funct,(Mat,Mat,Mat),(A,P,*C));CHKERRQ(ierr); 44 45 ierr = PetscLogEventEnd(MATSeqAIJ_PtAP,A,P,0,0);CHKERRQ(ierr); 46 47 PetscFunctionReturn(0); 48 } 49 50 /* 51 MatSeqAIJPtAPSymbolic - Creates the (i,j) structure of the SeqAIJ matrix product, C, 52 of SeqAIJ matrix A and matrix P, according to: 53 C = P^T * A * P; 54 55 Note: C is assumed to be uncreated. 56 If this is not the case, Destroy C before calling this routine. 57 */ 58 #undef __FUNCT__ 59 #define __FUNCT__ "MatSeqAIJPtAPSymbolic" 60 int MatSeqAIJPtAPSymbolic(Mat A,Mat P,Mat *C) { 61 int ierr; 62 char funct[80]; 63 64 PetscFunctionBegin; 65 66 PetscValidPointer(C); 67 68 PetscValidHeaderSpecific(A,MAT_COOKIE); 69 PetscValidType(A); 70 MatPreallocated(A); 71 if (!A->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 72 if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 73 74 PetscValidHeaderSpecific(P,MAT_COOKIE); 75 PetscValidType(P); 76 MatPreallocated(P); 77 if (!P->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 78 if (P->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 79 80 if (P->M!=A->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Matrix dimensions are incompatible, %d != %d",P->M,A->N); 81 if (A->M!=A->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Matrix 'A' must be square, %d != %d",A->M,A->N); 82 83 /* Query A for ApplyPtAP implementation based on types of P */ 84 ierr = PetscStrcpy(funct,"MatApplyPtAPSymbolic_seqaij_");CHKERRQ(ierr); 85 ierr = PetscStrcat(funct,P->type_name);CHKERRQ(ierr); 86 ierr = PetscTryMethod(A,funct,(Mat,Mat,Mat*),(A,P,C));CHKERRQ(ierr); 87 88 PetscFunctionReturn(0); 89 } 90 91 EXTERN_C_BEGIN 92 #undef __FUNCT__ 93 #define __FUNCT__ "MatApplyPtAPSymbolic_SeqAIJ_SeqAIJ" 94 int MatApplyPtAPSymbolic_SeqAIJ_SeqAIJ(Mat A,Mat P,Mat *C) { 95 int ierr; 96 FreeSpaceList free_space=PETSC_NULL,current_space=PETSC_NULL; 97 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data,*p=(Mat_SeqAIJ*)P->data,*c; 98 int *pti,*ptj,*ptJ,*ai=a->i,*aj=a->j,*ajj,*pi=p->i,*pj=p->j,*pjj; 99 int *ci,*cj,*denserow,*sparserow,*ptadenserow,*ptasparserow,*ptaj; 100 int an=A->N,am=A->M,pn=P->N,pm=P->M; 101 int i,j,k,ptnzi,arow,anzj,ptanzi,prow,pnzj,cnzi; 102 MatScalar *ca; 103 104 PetscFunctionBegin; 105 106 /* Start timer */ 107 ierr = PetscLogEventBegin(MATSeqAIJ_PtAPSymbolic,A,P,0,0);CHKERRQ(ierr); 108 109 /* Get ij structure of P^T */ 110 ierr = MatGetSymbolicTranspose_SeqAIJ(P,&pti,&ptj);CHKERRQ(ierr); 111 ptJ=ptj; 112 113 /* Allocate ci array, arrays for fill computation and */ 114 /* free space for accumulating nonzero column info */ 115 ierr = PetscMalloc((pn+1)*sizeof(int),&ci);CHKERRQ(ierr); 116 ci[0] = 0; 117 118 ierr = PetscMalloc((2*pn+2*an+1)*sizeof(int),&ptadenserow);CHKERRQ(ierr); 119 ierr = PetscMemzero(ptadenserow,(2*pn+2*an+1)*sizeof(int));CHKERRQ(ierr); 120 ptasparserow = ptadenserow + an; 121 denserow = ptasparserow + an; 122 sparserow = denserow + pn; 123 124 /* Set initial free space to be nnz(A) scaled by aspect ratio of P. */ 125 /* This should be reasonable if sparsity of PtAP is similar to that of A. */ 126 ierr = GetMoreSpace((ai[am]/pm)*pn,&free_space); 127 current_space = free_space; 128 129 /* Determine symbolic info for each row of C: */ 130 for (i=0;i<pn;i++) { 131 ptnzi = pti[i+1] - pti[i]; 132 ptanzi = 0; 133 /* Determine symbolic row of PtA: */ 134 for (j=0;j<ptnzi;j++) { 135 arow = *ptJ++; 136 anzj = ai[arow+1] - ai[arow]; 137 ajj = aj + ai[arow]; 138 for (k=0;k<anzj;k++) { 139 if (!ptadenserow[ajj[k]]) { 140 ptadenserow[ajj[k]] = -1; 141 ptasparserow[ptanzi++] = ajj[k]; 142 } 143 } 144 } 145 /* Using symbolic info for row of PtA, determine symbolic info for row of C: */ 146 ptaj = ptasparserow; 147 cnzi = 0; 148 for (j=0;j<ptanzi;j++) { 149 prow = *ptaj++; 150 pnzj = pi[prow+1] - pi[prow]; 151 pjj = pj + pi[prow]; 152 for (k=0;k<pnzj;k++) { 153 if (!denserow[pjj[k]]) { 154 denserow[pjj[k]] = -1; 155 sparserow[cnzi++] = pjj[k]; 156 } 157 } 158 } 159 160 /* sort sparserow */ 161 ierr = PetscSortInt(cnzi,sparserow);CHKERRQ(ierr); 162 163 /* If free space is not available, make more free space */ 164 /* Double the amount of total space in the list */ 165 if (current_space->local_remaining<cnzi) { 166 ierr = GetMoreSpace(current_space->total_array_size,¤t_space);CHKERRQ(ierr); 167 } 168 169 /* Copy data into free space, and zero out denserows */ 170 ierr = PetscMemcpy(current_space->array,sparserow,cnzi*sizeof(int));CHKERRQ(ierr); 171 current_space->array += cnzi; 172 current_space->local_used += cnzi; 173 current_space->local_remaining -= cnzi; 174 175 for (j=0;j<ptanzi;j++) { 176 ptadenserow[ptasparserow[j]] = 0; 177 } 178 for (j=0;j<cnzi;j++) { 179 denserow[sparserow[j]] = 0; 180 } 181 /* Aside: Perhaps we should save the pta info for the numerical factorization. */ 182 /* For now, we will recompute what is needed. */ 183 ci[i+1] = ci[i] + cnzi; 184 } 185 /* nnz is now stored in ci[ptm], column indices are in the list of free space */ 186 /* Allocate space for cj, initialize cj, and */ 187 /* destroy list of free space and other temporary array(s) */ 188 ierr = PetscMalloc((ci[pn]+1)*sizeof(int),&cj);CHKERRQ(ierr); 189 ierr = MakeSpaceContiguous(&free_space,cj);CHKERRQ(ierr); 190 ierr = PetscFree(ptadenserow);CHKERRQ(ierr); 191 192 /* Allocate space for ca */ 193 ierr = PetscMalloc((ci[pn]+1)*sizeof(MatScalar),&ca);CHKERRQ(ierr); 194 ierr = PetscMemzero(ca,(ci[pn]+1)*sizeof(MatScalar));CHKERRQ(ierr); 195 196 /* put together the new matrix */ 197 ierr = MatCreateSeqAIJWithArrays(A->comm,pn,pn,ci,cj,ca,C);CHKERRQ(ierr); 198 199 /* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */ 200 /* Since these are PETSc arrays, change flags to free them as necessary. */ 201 c = (Mat_SeqAIJ *)((*C)->data); 202 c->freedata = PETSC_TRUE; 203 c->nonew = 0; 204 205 /* Clean up. */ 206 ierr = MatRestoreSymbolicTranspose_SeqAIJ(P,&pti,&ptj);CHKERRQ(ierr); 207 208 ierr = PetscLogEventEnd(MATSeqAIJ_PtAPSymbolic,A,P,0,0);CHKERRQ(ierr); 209 PetscFunctionReturn(0); 210 } 211 EXTERN_C_END 212 213 #include "src/mat/impls/maij/maij.h" 214 EXTERN_C_BEGIN 215 #undef __FUNCT__ 216 #define __FUNCT__ "MatApplyPtAPSymbolic_SeqAIJ_SeqMAIJ" 217 int MatApplyPtAPSymbolic_SeqAIJ_SeqMAIJ(Mat A,Mat PP,Mat *C) { 218 int ierr; 219 FreeSpaceList free_space=PETSC_NULL,current_space=PETSC_NULL; 220 Mat_SeqMAIJ *pp=(Mat_SeqMAIJ*)PP->data; 221 Mat P=pp->AIJ; 222 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data,*p=(Mat_SeqAIJ*)P->data,*c; 223 int *pti,*ptj,*ptJ,*ai=a->i,*aj=a->j,*ajj,*pi=p->i,*pj=p->j,*pjj; 224 int *ci,*cj,*denserow,*sparserow,*ptadenserow,*ptasparserow,*ptaj; 225 int an=A->N,am=A->M,pn=P->N,pm=P->M,ppdof=pp->dof; 226 int i,j,k,dof,pdof,ptnzi,arow,anzj,ptanzi,prow,pnzj,cnzi; 227 MatScalar *ca; 228 229 PetscFunctionBegin; 230 /* Start timer */ 231 ierr = PetscLogEventBegin(MATSeqAIJ_PtAPSymbolic,A,PP,0,0);CHKERRQ(ierr); 232 233 /* Get ij structure of P^T */ 234 ierr = MatGetSymbolicTranspose_SeqAIJ(P,&pti,&ptj);CHKERRQ(ierr); 235 236 /* Allocate ci array, arrays for fill computation and */ 237 /* free space for accumulating nonzero column info */ 238 ierr = PetscMalloc((pn+1)*sizeof(int),&ci);CHKERRQ(ierr); 239 ci[0] = 0; 240 241 ierr = PetscMalloc((2*pn+2*an+1)*sizeof(int),&ptadenserow);CHKERRQ(ierr); 242 ierr = PetscMemzero(ptadenserow,(2*pn+2*an+1)*sizeof(int));CHKERRQ(ierr); 243 ptasparserow = ptadenserow + an; 244 denserow = ptasparserow + an; 245 sparserow = denserow + pn; 246 247 /* Set initial free space to be nnz(A) scaled by aspect ratio of P. */ 248 /* This should be reasonable if sparsity of PtAP is similar to that of A. */ 249 ierr = GetMoreSpace((ai[am]/pm)*pn,&free_space); 250 current_space = free_space; 251 252 /* Determine symbolic info for each row of C: */ 253 for (i=0;i<pn/ppdof;i++) { 254 ptnzi = pti[i+1] - pti[i]; 255 ptanzi = 0; 256 ptJ = ptj + pti[i]; 257 for (dof=0;dof<ppdof;dof++) { 258 /* Determine symbolic row of PtA: */ 259 for (j=0;j<ptnzi;j++) { 260 arow = ptJ[j] + dof; 261 anzj = ai[arow+1] - ai[arow]; 262 ajj = aj + ai[arow]; 263 for (k=0;k<anzj;k++) { 264 if (!ptadenserow[ajj[k]]) { 265 ptadenserow[ajj[k]] = -1; 266 ptasparserow[ptanzi++] = ajj[k]; 267 } 268 } 269 } 270 /* Using symbolic info for row of PtA, determine symbolic info for row of C: */ 271 ptaj = ptasparserow; 272 cnzi = 0; 273 for (j=0;j<ptanzi;j++) { 274 pdof = *ptaj%dof; 275 prow = (*ptaj++)/dof; 276 pnzj = pi[prow+1] - pi[prow]; 277 pjj = pj + pi[prow]; 278 for (k=0;k<pnzj;k++) { 279 if (!denserow[pjj[k]+pdof]) { 280 denserow[pjj[k]+pdof] = -1; 281 sparserow[cnzi++] = pjj[k]+pdof; 282 } 283 } 284 } 285 286 /* sort sparserow */ 287 ierr = PetscSortInt(cnzi,sparserow);CHKERRQ(ierr); 288 289 /* If free space is not available, make more free space */ 290 /* Double the amount of total space in the list */ 291 if (current_space->local_remaining<cnzi) { 292 ierr = GetMoreSpace(current_space->total_array_size,¤t_space);CHKERRQ(ierr); 293 } 294 295 /* Copy data into free space, and zero out denserows */ 296 ierr = PetscMemcpy(current_space->array,sparserow,cnzi*sizeof(int));CHKERRQ(ierr); 297 current_space->array += cnzi; 298 current_space->local_used += cnzi; 299 current_space->local_remaining -= cnzi; 300 301 for (j=0;j<ptanzi;j++) { 302 ptadenserow[ptasparserow[j]] = 0; 303 } 304 for (j=0;j<cnzi;j++) { 305 denserow[sparserow[j]] = 0; 306 } 307 /* Aside: Perhaps we should save the pta info for the numerical factorization. */ 308 /* For now, we will recompute what is needed. */ 309 ci[i+1+dof] = ci[i+dof] + cnzi; 310 } 311 } 312 /* nnz is now stored in ci[ptm], column indices are in the list of free space */ 313 /* Allocate space for cj, initialize cj, and */ 314 /* destroy list of free space and other temporary array(s) */ 315 ierr = PetscMalloc((ci[pn]+1)*sizeof(int),&cj);CHKERRQ(ierr); 316 ierr = MakeSpaceContiguous(&free_space,cj);CHKERRQ(ierr); 317 ierr = PetscFree(ptadenserow);CHKERRQ(ierr); 318 319 /* Allocate space for ca */ 320 ierr = PetscMalloc((ci[pn]+1)*sizeof(MatScalar),&ca);CHKERRQ(ierr); 321 ierr = PetscMemzero(ca,(ci[pn]+1)*sizeof(MatScalar));CHKERRQ(ierr); 322 323 /* put together the new matrix */ 324 ierr = MatCreateSeqAIJWithArrays(A->comm,pn,pn,ci,cj,ca,C);CHKERRQ(ierr); 325 326 /* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */ 327 /* Since these are PETSc arrays, change flags to free them as necessary. */ 328 c = (Mat_SeqAIJ *)((*C)->data); 329 c->freedata = PETSC_TRUE; 330 c->nonew = 0; 331 332 /* Clean up. */ 333 ierr = MatRestoreSymbolicTranspose_SeqAIJ(P,&pti,&ptj);CHKERRQ(ierr); 334 335 ierr = PetscLogEventEnd(MATSeqAIJ_PtAPSymbolic,A,PP,0,0);CHKERRQ(ierr); 336 PetscFunctionReturn(0); 337 } 338 EXTERN_C_END 339 340 /* 341 MatSeqAIJPtAPNumeric - Computes the SeqAIJ matrix product, C, 342 of SeqAIJ matrix A and matrix P, according to: 343 C = P^T * A * P 344 Note: C must have been created by calling MatSeqAIJApplyPtAPSymbolic. 345 */ 346 #undef __FUNCT__ 347 #define __FUNCT__ "MatSeqAIJPtAPNumeric" 348 int MatSeqAIJPtAPNumeric(Mat A,Mat P,Mat C) { 349 int ierr; 350 char funct[80]; 351 352 PetscFunctionBegin; 353 354 PetscValidHeaderSpecific(A,MAT_COOKIE); 355 PetscValidType(A); 356 MatPreallocated(A); 357 if (!A->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 358 if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 359 360 PetscValidHeaderSpecific(P,MAT_COOKIE); 361 PetscValidType(P); 362 MatPreallocated(P); 363 if (!P->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 364 if (P->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 365 366 PetscValidHeaderSpecific(C,MAT_COOKIE); 367 PetscValidType(C); 368 MatPreallocated(C); 369 if (!C->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 370 if (C->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 371 372 if (P->N!=C->M) SETERRQ2(PETSC_ERR_ARG_SIZ,"Matrix dimensions are incompatible, %d != %d",P->N,C->M); 373 if (P->M!=A->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Matrix dimensions are incompatible, %d != %d",P->M,A->N); 374 if (A->M!=A->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Matrix 'A' must be square, %d != %d",A->M,A->N); 375 if (P->N!=C->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Matrix dimensions are incompatible, %d != %d",P->N,C->N); 376 377 /* Query A for ApplyPtAP implementation based on types of P */ 378 ierr = PetscStrcpy(funct,"MatApplyPtAPNumeric_seqaij_");CHKERRQ(ierr); 379 ierr = PetscStrcat(funct,P->type_name);CHKERRQ(ierr); 380 ierr = PetscTryMethod(A,funct,(Mat,Mat,Mat),(A,P,C));CHKERRQ(ierr); 381 382 PetscFunctionReturn(0); 383 } 384 385 EXTERN_C_BEGIN 386 #undef __FUNCT__ 387 #define __FUNCT__ "MatApplyPtAPNumeric_SeqAIJ_SeqAIJ" 388 int MatApplyPtAPNumeric_SeqAIJ_SeqAIJ(Mat A,Mat P,Mat C) { 389 int ierr,flops=0; 390 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 391 Mat_SeqAIJ *p = (Mat_SeqAIJ *) P->data; 392 Mat_SeqAIJ *c = (Mat_SeqAIJ *) C->data; 393 int *ai=a->i,*aj=a->j,*apj,*apjdense,*pi=p->i,*pj=p->j,*pJ=p->j,*pjj; 394 int *ci=c->i,*cj=c->j,*cjj; 395 int am=A->M,cn=C->N,cm=C->M; 396 int i,j,k,anzi,pnzi,apnzj,nextap,pnzj,prow,crow; 397 MatScalar *aa=a->a,*apa,*pa=p->a,*pA=p->a,*paj,*ca=c->a,*caj; 398 399 PetscFunctionBegin; 400 ierr = PetscLogEventBegin(MATSeqAIJ_PtAPNumeric,A,P,C,0);CHKERRQ(ierr); 401 402 /* Allocate temporary array for storage of one row of A*P */ 403 ierr = PetscMalloc(cn*(sizeof(MatScalar)+2*sizeof(int)),&apa);CHKERRQ(ierr); 404 ierr = PetscMemzero(apa,cn*(sizeof(MatScalar)+2*sizeof(int)));CHKERRQ(ierr); 405 406 apj = (int *)(apa + cn); 407 apjdense = apj + cn; 408 409 /* Clear old values in C */ 410 ierr = PetscMemzero(ca,ci[cm]*sizeof(MatScalar));CHKERRQ(ierr); 411 412 for (i=0;i<am;i++) { 413 /* Form sparse row of A*P */ 414 anzi = ai[i+1] - ai[i]; 415 apnzj = 0; 416 for (j=0;j<anzi;j++) { 417 prow = *aj++; 418 pnzj = pi[prow+1] - pi[prow]; 419 pjj = pj + pi[prow]; 420 paj = pa + pi[prow]; 421 for (k=0;k<pnzj;k++) { 422 if (!apjdense[pjj[k]]) { 423 apjdense[pjj[k]] = -1; 424 apj[apnzj++] = pjj[k]; 425 } 426 apa[pjj[k]] += (*aa)*paj[k]; 427 } 428 flops += 2*pnzj; 429 aa++; 430 } 431 432 /* Sort the j index array for quick sparse axpy. */ 433 ierr = PetscSortInt(apnzj,apj);CHKERRQ(ierr); 434 435 /* Compute P^T*A*P using outer product (P^T)[:,j]*(A*P)[j,:]. */ 436 pnzi = pi[i+1] - pi[i]; 437 for (j=0;j<pnzi;j++) { 438 nextap = 0; 439 crow = *pJ++; 440 cjj = cj + ci[crow]; 441 caj = ca + ci[crow]; 442 /* Perform sparse axpy operation. Note cjj includes apj. */ 443 for (k=0;nextap<apnzj;k++) { 444 if (cjj[k]==apj[nextap]) { 445 caj[k] += (*pA)*apa[apj[nextap++]]; 446 } 447 } 448 flops += 2*apnzj; 449 pA++; 450 } 451 452 /* Zero the current row info for A*P */ 453 for (j=0;j<apnzj;j++) { 454 apa[apj[j]] = 0.; 455 apjdense[apj[j]] = 0; 456 } 457 } 458 459 /* Assemble the final matrix and clean up */ 460 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 461 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 462 ierr = PetscFree(apa);CHKERRQ(ierr); 463 ierr = PetscLogFlops(flops);CHKERRQ(ierr); 464 ierr = PetscLogEventEnd(MATSeqAIJ_PtAPNumeric,A,P,C,0);CHKERRQ(ierr); 465 466 PetscFunctionReturn(0); 467 } 468 EXTERN_C_END 469 470 #undef __FUNCT__ 471 #define __FUNCT__ "RegisterApplyPtAPRoutines_Private" 472 int RegisterApplyPtAPRoutines_Private(Mat A) { 473 int ierr; 474 475 PetscFunctionBegin; 476 477 if (!MATSeqAIJ_PtAP) { 478 ierr = PetscLogEventRegister(&MATSeqAIJ_PtAP,"MatSeqAIJApplyPtAP",MAT_COOKIE);CHKERRQ(ierr); 479 } 480 481 if (!MATSeqAIJ_PtAPSymbolic) { 482 ierr = PetscLogEventRegister(&MATSeqAIJ_PtAPSymbolic,"MatSeqAIJApplyPtAPSymbolic",MAT_COOKIE);CHKERRQ(ierr); 483 } 484 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatApplyPtAPSymbolic_seqaij_seqaij", 485 "MatApplyPtAPSymbolic_SeqAIJ_SeqAIJ", 486 MatApplyPtAPSymbolic_SeqAIJ_SeqAIJ);CHKERRQ(ierr); 487 488 if (!MATSeqAIJ_PtAPNumeric) { 489 ierr = PetscLogEventRegister(&MATSeqAIJ_PtAPNumeric,"MatSeqAIJApplyPtAPNumeric",MAT_COOKIE);CHKERRQ(ierr); 490 } 491 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatApplyPtAPNumeric_seqaij_seqaij", 492 "MatApplyPtAPNumeric_SeqAIJ_SeqAIJ", 493 MatApplyPtAPNumeric_SeqAIJ_SeqAIJ);CHKERRQ(ierr); 494 PetscFunctionReturn(0); 495 } 496