1 /*$Id: matmatmult.c,v 1.15 2001/09/07 20:04:44 buschelm Exp $*/ 2 /* 3 Defines a matrix-matrix product for 2 SeqAIJ matrices 4 C = A * B 5 */ 6 7 #include "src/mat/impls/aij/seq/aij.h" 8 9 typedef struct _Space *FreeSpaceList; 10 typedef struct _Space { 11 FreeSpaceList more_space; 12 int *array; 13 int *array_head; 14 int total_array_size; 15 int local_used; 16 int local_remaining; 17 } FreeSpace; 18 19 #undef __FUNCT__ 20 #define __FUNCT__ "GetMoreSpace" 21 int GetMoreSpace(int size,FreeSpaceList *list) { 22 FreeSpaceList a; 23 int ierr; 24 25 PetscFunctionBegin; 26 ierr = PetscMalloc(sizeof(FreeSpace),&a);CHKERRQ(ierr); 27 ierr = PetscMalloc(size*sizeof(int),&(a->array_head));CHKERRQ(ierr); 28 a->array = a->array_head; 29 a->local_remaining = size; 30 a->local_used = 0; 31 a->total_array_size = 0; 32 a->more_space = NULL; 33 34 if (*list) { 35 (*list)->more_space = a; 36 a->total_array_size = (*list)->total_array_size; 37 } 38 39 a->total_array_size += size; 40 *list = a; 41 PetscFunctionReturn(0); 42 } 43 44 #undef __FUNCT__ 45 #define __FUNCT__ "MakeSpaceContiguous" 46 int MakeSpaceContiguous(int *space,FreeSpaceList *head) { 47 FreeSpaceList a; 48 int ierr; 49 50 PetscFunctionBegin; 51 while ((*head)!=NULL) { 52 a = (*head)->more_space; 53 ierr = PetscMemcpy(space,(*head)->array_head,((*head)->local_used)*sizeof(int));CHKERRQ(ierr); 54 space += (*head)->local_used; 55 ierr = PetscFree((*head)->array_head);CHKERRQ(ierr); 56 ierr = PetscFree(*head);CHKERRQ(ierr); 57 *head = a; 58 } 59 PetscFunctionReturn(0); 60 } 61 62 static int logkey_matmatmult_symbolic = 0; 63 static int logkey_matmatmult_numeric = 0; 64 65 /* 66 MatMatMult_SeqAIJ_SeqAIJ_Symbolic - Forms the symbolic product of two SeqAIJ matrices 67 C=A*B; 68 69 Note: C is assumed to be uninitialized. 70 If this is not the case, Destroy C before calling this routine. 71 */ 72 #undef __FUNCT__ 73 #define __FUNCT__ "MatMatMult_SeqAIJ_SeqAIJ_Symbolic" 74 int MatMatMult_SeqAIJ_SeqAIJ_Symbolic(Mat A,Mat B,Mat *C) 75 { 76 int ierr; 77 FreeSpaceList free_space=PETSC_NULL,current_space=PETSC_NULL; 78 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data,*b=(Mat_SeqAIJ*)B->data,*c; 79 int aishift=a->indexshift,bishift=b->indexshift; 80 int *ai=a->i,*aj=a->j,*bi=b->i,*bj=b->j,*bjj; 81 int *ci,*cj,*densefill,*sparsefill; 82 int an=A->N,am=A->M,bn=B->N,bm=B->M; 83 int i,j,k,anzi,brow,bnzj,cnzi; 84 MatScalar *ca; 85 86 PetscFunctionBegin; 87 /* some error checking which could be moved into interface layer */ 88 if (aishift || bishift) SETERRQ(PETSC_ERR_SUP,"Shifted matrix indices are not supported."); 89 if (an!=bm) SETERRQ2(PETSC_ERR_ARG_SIZ,"Matrix dimensions are incompatible, %d != %d",an,bm); 90 91 if (!logkey_matmatmult_symbolic) { 92 ierr = PetscLogEventRegister(&logkey_matmatmult_symbolic,"MatMatMult_Symbolic",MAT_COOKIE);CHKERRQ(ierr); 93 } 94 ierr = PetscLogEventBegin(logkey_matmatmult_symbolic,A,B,0,0);CHKERRQ(ierr); 95 96 /* Set up */ 97 /* Allocate ci array, arrays for fill computation and */ 98 /* free space for accumulating nonzero column info */ 99 ierr = PetscMalloc(((am+1)+1)*sizeof(int),&ci);CHKERRQ(ierr); 100 ci[0] = 0; 101 102 ierr = PetscMalloc((2*bn+1)*sizeof(int),&densefill);CHKERRQ(ierr); 103 ierr = PetscMemzero(densefill,(2*bn+1)*sizeof(int));CHKERRQ(ierr); 104 sparsefill = densefill + bn; 105 106 /* Initial FreeSpace size is nnz(B)=bi[bm] */ 107 ierr = GetMoreSpace(bi[bm],&free_space);CHKERRQ(ierr); 108 current_space = free_space; 109 110 /* Determine fill for each row: */ 111 for (i=0;i<am;i++) { 112 anzi = ai[i+1] - ai[i]; 113 cnzi = 0; 114 for (j=0;j<anzi;j++) { 115 brow = *aj++; 116 bnzj = bi[brow+1] - bi[brow]; 117 bjj = bj + bi[brow]; 118 for (k=0;k<bnzj;k++) { 119 /* If column is not marked, mark it in compressed and uncompressed locations. */ 120 /* For simplicity, leave uncompressed row unsorted until finished with row, */ 121 /* and increment nonzero count for this row. */ 122 if (!densefill[bjj[k]]) { 123 densefill[bjj[k]] = -1; 124 sparsefill[cnzi++] = bjj[k]; 125 } 126 } 127 } 128 129 /* sort sparsefill */ 130 ierr = PetscSortInt(cnzi,sparsefill);CHKERRQ(ierr); 131 132 /* If free space is not available, make more free space */ 133 /* Double the amount of total space in the list */ 134 if (current_space->local_remaining<cnzi) { 135 ierr = GetMoreSpace(current_space->total_array_size,¤t_space);CHKERRQ(ierr); 136 } 137 138 /* Copy data into free space, and zero out densefill */ 139 ierr = PetscMemcpy(current_space->array,sparsefill,cnzi*sizeof(int));CHKERRQ(ierr); 140 current_space->array += cnzi; 141 current_space->local_used += cnzi; 142 current_space->local_remaining -= cnzi; 143 for (j=0;j<cnzi;j++) { 144 densefill[sparsefill[j]] = 0; 145 } 146 ci[i+1] = ci[i] + cnzi; 147 } 148 149 /* nnz is now stored in ci[am], column indices are in the list of free space */ 150 /* Allocate space for cj, initialize cj, and */ 151 /* destroy list of free space and other temporary array(s) */ 152 ierr = PetscMalloc((ci[am]+1)*sizeof(int),&cj);CHKERRQ(ierr); 153 ierr = MakeSpaceContiguous(cj,&free_space);CHKERRQ(ierr); 154 ierr = PetscFree(densefill);CHKERRQ(ierr); 155 156 /* Allocate space for ca */ 157 ierr = PetscMalloc((ci[am]+1)*sizeof(MatScalar),&ca);CHKERRQ(ierr); 158 ierr = PetscMemzero(ca,(ci[am]+1)*sizeof(MatScalar));CHKERRQ(ierr); 159 160 /* put together the new matrix */ 161 ierr = MatCreateSeqAIJWithArrays(A->comm,am,bn,ci,cj,ca,C);CHKERRQ(ierr); 162 163 /* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */ 164 /* These are PETSc arrays, so change flags so arrays can be deleted by PETSc */ 165 c = (Mat_SeqAIJ *)((*C)->data); 166 c->freedata = PETSC_TRUE; 167 c->nonew = 0; 168 169 ierr = PetscLogEventEnd(logkey_matmatmult_symbolic,A,B,0,0);CHKERRQ(ierr); 170 PetscFunctionReturn(0); 171 } 172 173 /* 174 MatMatMult_SeqAIJ_SeqAIJ_Numeric - Forms the numeric product of two SeqAIJ matrices 175 C=A*B; 176 Note: C must have been created by calling MatMatMult_SeqAIJ_SeqAIJ_Symbolic. 177 */ 178 #undef __FUNCT__ 179 #define __FUNCT__ "MatMatMult_SeqAIJ_SeqAIJ_Numeric" 180 int MatMatMult_SeqAIJ_SeqAIJ_Numeric(Mat A,Mat B,Mat C) 181 { 182 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 183 Mat_SeqAIJ *b = (Mat_SeqAIJ *)B->data; 184 Mat_SeqAIJ *c = (Mat_SeqAIJ *)C->data; 185 int aishift=a->indexshift,bishift=b->indexshift,cishift=c->indexshift; 186 int *ai=a->i,*aj=a->j,*bi=b->i,*bj=b->j,*bjj,*ci=c->i,*cj=c->j; 187 int an=A->N,am=A->M,bn=B->N,bm=B->M,cn=C->N,cm=C->M; 188 int ierr,i,j,k,anzi,bnzi,cnzi,brow,flops; 189 MatScalar *aa=a->a,*ba=b->a,*baj,*ca=c->a,*temp; 190 191 PetscFunctionBegin; 192 193 /* This error checking should be unnecessary if the symbolic was performed */ 194 if (aishift || bishift || cishift) SETERRQ(PETSC_ERR_SUP,"Shifted matrix indices are not supported."); 195 if (am!=cm) SETERRQ2(PETSC_ERR_ARG_SIZ,"Matrix dimensions are incompatible, %d != %d",am,cm); 196 if (an!=bm) SETERRQ2(PETSC_ERR_ARG_SIZ,"Matrix dimensions are incompatible, %d != %d",an,bm); 197 if (bn!=cn) SETERRQ2(PETSC_ERR_ARG_SIZ,"Matrix dimensions are incompatible, %d != %d",bn,cn); 198 199 if (!logkey_matmatmult_numeric) { 200 ierr = PetscLogEventRegister(&logkey_matmatmult_numeric,"MatMatMult_Numeric",MAT_COOKIE);CHKERRQ(ierr); 201 } 202 ierr = PetscLogEventBegin(logkey_matmatmult_numeric,A,B,C,0);CHKERRQ(ierr); 203 flops = 0; 204 /* Allocate temp accumulation space to avoid searching for nonzero columns in C */ 205 ierr = PetscMalloc((cn+1)*sizeof(MatScalar),&temp);CHKERRQ(ierr); 206 ierr = PetscMemzero(temp,cn*sizeof(MatScalar));CHKERRQ(ierr); 207 /* Traverse A row-wise. */ 208 /* Build the ith row in C by summing over nonzero columns in A, */ 209 /* the rows of B corresponding to nonzeros of A. */ 210 for (i=0;i<am;i++) { 211 anzi = ai[i+1] - ai[i]; 212 for (j=0;j<anzi;j++) { 213 brow = *aj++; 214 bnzi = bi[brow+1] - bi[brow]; 215 bjj = bj + bi[brow]; 216 baj = ba + bi[brow]; 217 for (k=0;k<bnzi;k++) { 218 temp[bjj[k]] += (*aa)*baj[k]; 219 } 220 flops += 2*bnzi; 221 aa++; 222 } 223 /* Store row back into C, and re-zero temp */ 224 cnzi = ci[i+1] - ci[i]; 225 for (j=0;j<cnzi;j++) { 226 ca[j] = temp[cj[j]]; 227 temp[cj[j]] = 0.0; 228 } 229 ca += cnzi; 230 cj += cnzi; 231 } 232 /* Free temp */ 233 ierr = PetscFree(temp);CHKERRQ(ierr); 234 ierr = PetscLogFlops(flops);CHKERRQ(ierr); 235 ierr = PetscLogEventEnd(logkey_matmatmult_numeric,A,B,C,0);CHKERRQ(ierr); 236 PetscFunctionReturn(0); 237 } 238 239 #undef __FUNCT__ 240 #define __FUNCT__ "MatMatMult_SeqAIJ_SeqAIJ" 241 int MatMatMult_SeqAIJ_SeqAIJ(Mat A,Mat B,Mat *C) { 242 int ierr; 243 244 PetscFunctionBegin; 245 ierr = MatMatMult_SeqAIJ_SeqAIJ_Symbolic(A,B,C);CHKERRQ(ierr); 246 ierr = MatMatMult_SeqAIJ_SeqAIJ_Numeric(A,B,*C);CHKERRQ(ierr); 247 PetscFunctionReturn(0); 248 } 249 250 static int logkey_matapplyptap_symbolic = 0; 251 static int logkey_matapplyptap_numeric = 0; 252 253 #undef __FUNCT__ 254 #define __FUNCT__ "MatApplyPtAP_SeqAIJ_Symbolic" 255 int MatApplyPtAP_SeqAIJ_Symbolic(Mat A,Mat P,Mat *C) { 256 int ierr; 257 FreeSpaceList free_space=PETSC_NULL,current_space=PETSC_NULL; 258 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data,*p=(Mat_SeqAIJ*)P->data,*c; 259 int aishift=a->indexshift,pishift=p->indexshift; 260 int *pti,*ptj,*ptfill,*ai=a->i,*aj=a->j,*ajj,*pi=p->i,*pj=p->j,*pjj; 261 int *ci,*cj,*densefill,*sparsefill,*ptadensefill,*ptasparsefill,*ptaj; 262 int an=A->N,am=A->M,pn=P->N,pm=P->M; 263 int i,j,k,ptnzi,arow,anzj,ptanzi,prow,pnzj,cnzi; 264 MatScalar *ca; 265 266 PetscFunctionBegin; 267 268 /* some error checking which could be moved into interface layer */ 269 if (aishift || pishift) SETERRQ(PETSC_ERR_SUP,"Shifted matrix indices are not supported."); 270 if (pm!=an) SETERRQ2(PETSC_ERR_ARG_SIZ,"Matrix dimensions are incompatible, %d != %d",pm,an); 271 if (am!=an) SETERRQ2(PETSC_ERR_ARG_SIZ,"Matrix 'A' must be square, %d != %d",am, an); 272 273 if (!logkey_matapplyptap_symbolic) { 274 ierr = PetscLogEventRegister(&logkey_matapplyptap_symbolic,"MatApplyPtAP_Symbolic",MAT_COOKIE);CHKERRQ(ierr); 275 } 276 ierr = PetscLogEventBegin(logkey_matapplyptap_symbolic,A,P,0,0);CHKERRQ(ierr); 277 278 /* Create ij structure of P^T */ 279 /* Recall in P^T there are pn rows and pi[pm] nonzeros. */ 280 ierr = PetscMalloc((pn+1+pi[pm])*sizeof(int),&pti);CHKERRQ(ierr); 281 ierr = PetscMemzero(pti,(pn+1+pi[pm])*sizeof(int));CHKERRQ(ierr); 282 ptj = pti + pn+1; 283 284 /* Walk through pj and count ## of non-zeros in each row of P^T. */ 285 for (i=0;i<pi[pm]-1;i++) { 286 pti[pj[i]+1] += 1; 287 } 288 /* Form pti for csr format of P^T. */ 289 for (i=0;i<pm;i++) { 290 pti[i+1] += pti[i]; 291 } 292 293 /* Allocate temporary space for next insert location in each row of P^T. */ 294 ierr = PetscMalloc(pn*sizeof(int),&ptfill);CHKERRQ(ierr); 295 ierr = PetscMemcpy(ptfill,pti,pn*sizeof(int));CHKERRQ(ierr); 296 297 /* Walk through P row-wise and mark nonzero entries of P^T. */ 298 for (i=0;i<pm;i++) { 299 pnzj = pi[i+1] - pi[i]; 300 for (j=0;j<pnzj;j++) { 301 ptj[ptfill[j]] = i; 302 ptfill[j] += 1; 303 } 304 } 305 306 /* Clean-up temporary space. */ 307 ierr = PetscFree(ptfill);CHKERRQ(ierr); 308 309 /* Allocate ci array, arrays for fill computation and */ 310 /* free space for accumulating nonzero column info */ 311 ierr = PetscMalloc(((pn+1)*1)*sizeof(int),&ci);CHKERRQ(ierr); 312 ci[0] = 0; 313 314 ierr = PetscMalloc((2*pn+2*an+1)*sizeof(int),&ptadensefill);CHKERRQ(ierr); 315 ierr = PetscMemzero(ptadensefill,(2*pn+2*an+1)*sizeof(int));CHKERRQ(ierr); 316 ptasparsefill = ptadensefill + an; 317 densefill = ptasparsefill + an; 318 sparsefill = densefill + pn; 319 320 /* Set initial free space to be nnz(A) scaled by aspect ratio of P. */ 321 /* Reason: Take pn/pm = 1/2. */ 322 /* P^T*A*P will take A(NxN) and create C(N/2xN/2). */ 323 /* If C has same sparsity pattern as A, nnz(C)~1/2*nnz(A). */ 324 /* Is this reasonable???? */ 325 ierr = GetMoreSpace((ai[am]*pn)/pm,&free_space); 326 current_space = free_space; 327 328 /* Determine fill for each row of C: */ 329 for (i=0;i<pn;i++) { 330 ptnzi = pti[i+1] - pti[i]; 331 ptanzi = 0; 332 /* Determine fill for row of PtA: */ 333 for (j=0;j<ptnzi;j++) { 334 arow = *ptj++; 335 anzj = ai[arow+1] - ai[arow]; 336 ajj = aj + ai[arow]; 337 for (k=0;k<anzj;k++) { 338 if (!ptadensefill[ajj[k]]) { 339 ptadensefill[ajj[k]] = -1; 340 ptasparsefill[ptanzi++] = ajj[k]; 341 } 342 } 343 } 344 /* Using fill info for row of PtA, determine fill for row of C: */ 345 ptaj = ptasparsefill; 346 cnzi = 0; 347 for (j=0;j<ptanzi;j++) { 348 prow = *ptaj++; 349 pnzj = pi[prow+1] - pi[prow]; 350 pjj = pj + pi[prow]; 351 for (k=0;k<pnzj;k++) { 352 if (!densefill[pjj[k]]) { 353 densefill[pjj[k]] = -1; 354 sparsefill[cnzi++] = pjj[k]; 355 } 356 } 357 } 358 359 /* sort sparsefill */ 360 ierr = PetscSortInt(cnzi,sparsefill);CHKERRQ(ierr); 361 362 /* If free space is not available, make more free space */ 363 /* Double the amount of total space in the list */ 364 if (current_space->local_remaining<cnzi) { 365 ierr = GetMoreSpace(current_space->total_array_size,¤t_space);CHKERRQ(ierr); 366 } 367 368 /* Copy data into free space, and zero out densefills */ 369 ierr = PetscMemcpy(current_space->array,sparsefill,cnzi*sizeof(int));CHKERRQ(ierr); 370 current_space->array += cnzi; 371 current_space->local_used += cnzi; 372 current_space->local_remaining -= cnzi; 373 374 for (j=0;j<ptanzi;j++) { 375 ptadensefill[ptasparsefill[j]] = 0; 376 } 377 for (j=0;j<cnzi;j++) { 378 densefill[sparsefill[j]] = 0; 379 } 380 /* Aside: Perhaps we should save the pta info for the numerical factorization. */ 381 /* For now, we will recompute what is needed. */ 382 ci[i+1] = ci[i] + cnzi; 383 } 384 /* nnz is now stored in ci[ptm], column indices are in the list of free space */ 385 /* Allocate space for cj, initialize cj, and */ 386 /* destroy list of free space and other temporary array(s) */ 387 ierr = PetscMalloc((ci[pn]+1)*sizeof(int),&cj);CHKERRQ(ierr); 388 ierr = MakeSpaceContiguous(cj,&free_space);CHKERRQ(ierr); 389 ierr = PetscFree(ptadensefill);CHKERRQ(ierr); 390 391 /* Allocate space for ca */ 392 ierr = PetscMalloc((ci[pn]+1)*sizeof(MatScalar),&ca);CHKERRQ(ierr); 393 ierr = PetscMemzero(ca,(ci[pn]+1)*sizeof(MatScalar));CHKERRQ(ierr); 394 395 /* put together the new matrix */ 396 ierr = MatCreateSeqAIJWithArrays(A->comm,pn,pn,ci,cj,ca,C);CHKERRQ(ierr); 397 398 /* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */ 399 /* Since these are PETSc arrays, change flags to free them as necessary. */ 400 c = (Mat_SeqAIJ *)((*C)->data); 401 c->freedata = PETSC_TRUE; 402 c->nonew = 0; 403 404 /* Clean up. */ 405 /* Perhaps we should attach the (i,j) info for P^T to P for future use. */ 406 /* For now, we won't. */ 407 ierr = PetscFree(pti); 408 409 ierr = PetscLogEventEnd(logkey_matapplyptap_symbolic,A,P,0,0);CHKERRQ(ierr); 410 PetscFunctionReturn(0); 411 } 412 413 #undef __FUNCT__ 414 #define __FUNCT__ "MatApplyPtAP_SeqAIJ_Numeric" 415 int MatApplyPtAP_SeqAIJ_Numeric(Mat A,Mat P,Mat C) { 416 int ierr,flops; 417 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 418 Mat_SeqAIJ *p = (Mat_SeqAIJ *) P->data; 419 Mat_SeqAIJ *c = (Mat_SeqAIJ *) C->data; 420 int aishift=a->indexshift,pishift=p->indexshift,cishift=c->indexshift; 421 int *ai=a->i,*aj=a->j,*apj,*pi=p->i,*pj=p->j,*pJ=p->j,*pjj,*ci=c->i,*cj=c->j,*cjj; 422 int an=A->N,am=A->M,pn=P->N,pm=P->M,cn=C->N,cm=C->M; 423 int i,j,k,anzi,pnzi,apnzj,pnzj,cnzj,prow,crow; 424 MatScalar *aa=a->a,*apa,*pa=p->a,*pA=p->a,*paj,*ca=c->a,*caj; 425 426 PetscFunctionBegin; 427 428 /* This error checking should be unnecessary if the symbolic was performed */ 429 if (aishift || pishift || cishift) SETERRQ(PETSC_ERR_SUP,"Shifted matrix indices are not supported."); 430 if (pn!=cm) SETERRQ2(PETSC_ERR_ARG_SIZ,"Matrix dimensions are incompatible, %d != %d",pn,cm); 431 if (pm!=an) SETERRQ2(PETSC_ERR_ARG_SIZ,"Matrix dimensions are incompatible, %d != %d",pm,an); 432 if (am!=an) SETERRQ2(PETSC_ERR_ARG_SIZ,"Matrix 'A' must be square, %d != %d",am, an); 433 if (pn!=cn) SETERRQ2(PETSC_ERR_ARG_SIZ,"Matrix dimensions are incompatible, %d != %d",pn, cn); 434 435 if (!logkey_matapplyptap_numeric) { 436 ierr = PetscLogEventRegister(&logkey_matapplyptap_numeric,"MatApplyPtAP_Numeric",MAT_COOKIE);CHKERRQ(ierr); 437 } 438 ierr = PetscLogEventBegin(logkey_matapplyptap_numeric,A,P,C,0);CHKERRQ(ierr); 439 flops = 0; 440 441 ierr = PetscMalloc(cn*(sizeof(MatScalar)+sizeof(int)),&apa);CHKERRQ(ierr); 442 ierr = PetscMemzero(apa,cn*(sizeof(MatScalar)+sizeof(int)));CHKERRQ(ierr); 443 apj = (int *)(apa + cn); 444 ierr = PetscMemzero(ca,ci[cm]*sizeof(MatScalar));CHKERRQ(ierr); 445 446 for (i=0;i<am;i++) { 447 CHKMEMQ; 448 /* Form sparse row of A*P */ 449 anzi = ai[i+1] - ai[i]; 450 apnzj = 0; 451 for (j=0;j<anzi;j++) { 452 prow = *aj++; 453 pnzj = pi[prow+1] - pi[prow]; 454 pjj = pj + pi[prow]; 455 paj = pa + pi[prow]; 456 for (k=0;k<pnzj;k++) { 457 if (!apa[pjj[k]]) { 458 CHKMEMQ; 459 apj[apnzj++]=pjj[k]; 460 CHKMEMQ; 461 } 462 CHKMEMQ; 463 apa[pjj[k]] += (*aa)*paj[k]; 464 CHKMEMQ; 465 } 466 flops += 2*pnzj; 467 aa++; 468 } 469 470 /* Sort the j index array for quick sparse axpy. */ 471 ierr = PetscSortInt(apnzj,apj);CHKERRQ(ierr); 472 473 /* Compute P^T*A*P using outer product (P^T)[:,j]*(A*P)[j,:]. */ 474 pnzi = pi[i+1] - pi[i]; 475 for (j=0;j<pnzi;j++) { 476 int nextap=0; 477 crow = *pJ++; 478 cnzj = ci[crow+1] - ci[crow]; 479 cjj = cj + ci[crow]; 480 caj = ca + ci[crow]; 481 /* Perform the sparse axpy operation. Note cjj includes apj. */ 482 for (k=0;k<cnzj;k++) { 483 if (cjj[k]==apj[nextap]) { 484 caj[k] += (*pA)*apa[apj[nextap++]]; 485 } 486 } 487 flops += 2*apnzj; 488 pA++; 489 } 490 491 for (j=0;j<apnzj;j++) { 492 CHKMEMQ; 493 apa[apj[j]] = 0.; 494 CHKMEMQ; 495 } 496 } 497 ierr = PetscFree(apa);CHKERRQ(ierr); 498 ierr = PetscLogFlops(flops);CHKERRQ(ierr); 499 ierr = PetscLogEventEnd(logkey_matapplyptap_numeric,A,P,C,0);CHKERRQ(ierr); 500 PetscFunctionReturn(0); 501 } 502 503 #undef __FUNCT__ 504 #define __FUNCT__ "MatApplyPtAP_SeqAIJ" 505 int MatApplyPtAP_SeqAIJ(Mat A,Mat P,Mat *C) { 506 int ierr; 507 508 PetscFunctionBegin; 509 ierr = MatApplyPtAP_SeqAIJ_Symbolic(A,P,C);CHKERRQ(ierr); 510 ierr = MatApplyPtAP_SeqAIJ_Numeric(A,P,*C);CHKERRQ(ierr); 511 PetscFunctionReturn(0); 512 } 513