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];i++) { 286 pti[pj[i]+1] += 1; 287 } 288 /* Form pti for csr format of P^T. */ 289 for (i=0;i<pn;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[*pj]] = i; 302 ptfill[*pj++] += 1; 303 } 304 } 305 pj = p->j; 306 307 /* Clean-up temporary space. */ 308 ierr = PetscFree(ptfill);CHKERRQ(ierr); 309 310 /* Allocate ci array, arrays for fill computation and */ 311 /* free space for accumulating nonzero column info */ 312 ierr = PetscMalloc(((pn+1)*1)*sizeof(int),&ci);CHKERRQ(ierr); 313 ci[0] = 0; 314 315 ierr = PetscMalloc((2*pn+2*an+1)*sizeof(int),&ptadensefill);CHKERRQ(ierr); 316 ierr = PetscMemzero(ptadensefill,(2*pn+2*an+1)*sizeof(int));CHKERRQ(ierr); 317 ptasparsefill = ptadensefill + an; 318 densefill = ptasparsefill + an; 319 sparsefill = densefill + pn; 320 321 /* Set initial free space to be nnz(A) scaled by aspect ratio of P. */ 322 /* Reason: Take pn/pm = 1/2. */ 323 /* P^T*A*P will take A(NxN) and create C(N/2xN/2). */ 324 /* If C has same sparsity pattern as A, nnz(C)~1/2*nnz(A). */ 325 /* Is this reasonable???? */ 326 ierr = GetMoreSpace((ai[am]*pn)/pm,&free_space); 327 current_space = free_space; 328 329 /* Determine fill for each row of C: */ 330 for (i=0;i<pn;i++) { 331 ptnzi = pti[i+1] - pti[i]; 332 ptanzi = 0; 333 /* Determine fill for row of PtA: */ 334 for (j=0;j<ptnzi;j++) { 335 arow = *ptj++; 336 anzj = ai[arow+1] - ai[arow]; 337 ajj = aj + ai[arow]; 338 for (k=0;k<anzj;k++) { 339 if (!ptadensefill[ajj[k]]) { 340 ptadensefill[ajj[k]] = -1; 341 ptasparsefill[ptanzi++] = ajj[k]; 342 } 343 } 344 } 345 /* Using fill info for row of PtA, determine fill for row of C: */ 346 ptaj = ptasparsefill; 347 cnzi = 0; 348 for (j=0;j<ptanzi;j++) { 349 prow = *ptaj++; 350 pnzj = pi[prow+1] - pi[prow]; 351 pjj = pj + pi[prow]; 352 for (k=0;k<pnzj;k++) { 353 if (!densefill[pjj[k]]) { 354 densefill[pjj[k]] = -1; 355 sparsefill[cnzi++] = pjj[k]; 356 } 357 } 358 } 359 360 /* sort sparsefill */ 361 ierr = PetscSortInt(cnzi,sparsefill);CHKERRQ(ierr); 362 363 /* If free space is not available, make more free space */ 364 /* Double the amount of total space in the list */ 365 if (current_space->local_remaining<cnzi) { 366 ierr = GetMoreSpace(current_space->total_array_size,¤t_space);CHKERRQ(ierr); 367 } 368 369 /* Copy data into free space, and zero out densefills */ 370 ierr = PetscMemcpy(current_space->array,sparsefill,cnzi*sizeof(int));CHKERRQ(ierr); 371 current_space->array += cnzi; 372 current_space->local_used += cnzi; 373 current_space->local_remaining -= cnzi; 374 375 for (j=0;j<ptanzi;j++) { 376 ptadensefill[ptasparsefill[j]] = 0; 377 } 378 for (j=0;j<cnzi;j++) { 379 densefill[sparsefill[j]] = 0; 380 } 381 /* Aside: Perhaps we should save the pta info for the numerical factorization. */ 382 /* For now, we will recompute what is needed. */ 383 ci[i+1] = ci[i] + cnzi; 384 } 385 /* nnz is now stored in ci[ptm], column indices are in the list of free space */ 386 /* Allocate space for cj, initialize cj, and */ 387 /* destroy list of free space and other temporary array(s) */ 388 ierr = PetscMalloc((ci[pn]+1)*sizeof(int),&cj);CHKERRQ(ierr); 389 ierr = MakeSpaceContiguous(cj,&free_space);CHKERRQ(ierr); 390 ierr = PetscFree(ptadensefill);CHKERRQ(ierr); 391 392 /* Allocate space for ca */ 393 ierr = PetscMalloc((ci[pn]+1)*sizeof(MatScalar),&ca);CHKERRQ(ierr); 394 ierr = PetscMemzero(ca,(ci[pn]+1)*sizeof(MatScalar));CHKERRQ(ierr); 395 396 /* put together the new matrix */ 397 ierr = MatCreateSeqAIJWithArrays(A->comm,pn,pn,ci,cj,ca,C);CHKERRQ(ierr); 398 399 /* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */ 400 /* Since these are PETSc arrays, change flags to free them as necessary. */ 401 c = (Mat_SeqAIJ *)((*C)->data); 402 c->freedata = PETSC_TRUE; 403 c->nonew = 0; 404 405 /* Clean up. */ 406 /* Perhaps we should attach the (i,j) info for P^T to P for future use. */ 407 /* For now, we won't. */ 408 ierr = PetscFree(pti); 409 410 ierr = PetscLogEventEnd(logkey_matapplyptap_symbolic,A,P,0,0);CHKERRQ(ierr); 411 PetscFunctionReturn(0); 412 } 413 414 #undef __FUNCT__ 415 #define __FUNCT__ "MatApplyPtAP_SeqAIJ_Numeric" 416 int MatApplyPtAP_SeqAIJ_Numeric(Mat A,Mat P,Mat C) { 417 int ierr,flops; 418 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 419 Mat_SeqAIJ *p = (Mat_SeqAIJ *) P->data; 420 Mat_SeqAIJ *c = (Mat_SeqAIJ *) C->data; 421 int aishift=a->indexshift,pishift=p->indexshift,cishift=c->indexshift; 422 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; 423 int an=A->N,am=A->M,pn=P->N,pm=P->M,cn=C->N,cm=C->M; 424 int i,j,k,anzi,pnzi,apnzj,pnzj,cnzj,prow,crow; 425 MatScalar *aa=a->a,*apa,*pa=p->a,*pA=p->a,*paj,*ca=c->a,*caj; 426 427 PetscFunctionBegin; 428 429 /* This error checking should be unnecessary if the symbolic was performed */ 430 if (aishift || pishift || cishift) SETERRQ(PETSC_ERR_SUP,"Shifted matrix indices are not supported."); 431 if (pn!=cm) SETERRQ2(PETSC_ERR_ARG_SIZ,"Matrix dimensions are incompatible, %d != %d",pn,cm); 432 if (pm!=an) SETERRQ2(PETSC_ERR_ARG_SIZ,"Matrix dimensions are incompatible, %d != %d",pm,an); 433 if (am!=an) SETERRQ2(PETSC_ERR_ARG_SIZ,"Matrix 'A' must be square, %d != %d",am, an); 434 if (pn!=cn) SETERRQ2(PETSC_ERR_ARG_SIZ,"Matrix dimensions are incompatible, %d != %d",pn, cn); 435 436 if (!logkey_matapplyptap_numeric) { 437 ierr = PetscLogEventRegister(&logkey_matapplyptap_numeric,"MatApplyPtAP_Numeric",MAT_COOKIE);CHKERRQ(ierr); 438 } 439 ierr = PetscLogEventBegin(logkey_matapplyptap_numeric,A,P,C,0);CHKERRQ(ierr); 440 flops = 0; 441 442 ierr = PetscMalloc(cn*(sizeof(MatScalar)+sizeof(int)),&apa);CHKERRQ(ierr); 443 ierr = PetscMemzero(apa,cn*(sizeof(MatScalar)+sizeof(int)));CHKERRQ(ierr); 444 apj = (int *)(apa + cn); 445 ierr = PetscMemzero(ca,ci[cm]*sizeof(MatScalar));CHKERRQ(ierr); 446 447 for (i=0;i<am;i++) { 448 CHKMEMQ; 449 /* Form sparse row of A*P */ 450 anzi = ai[i+1] - ai[i]; 451 apnzj = 0; 452 for (j=0;j<anzi;j++) { 453 prow = *aj++; 454 pnzj = pi[prow+1] - pi[prow]; 455 pjj = pj + pi[prow]; 456 paj = pa + pi[prow]; 457 for (k=0;k<pnzj;k++) { 458 if (!apa[pjj[k]]) { 459 CHKMEMQ; 460 apj[apnzj++]=pjj[k]; 461 CHKMEMQ; 462 } 463 CHKMEMQ; 464 apa[pjj[k]] += (*aa)*paj[k]; 465 CHKMEMQ; 466 } 467 flops += 2*pnzj; 468 aa++; 469 } 470 471 /* Sort the j index array for quick sparse axpy. */ 472 ierr = PetscSortInt(apnzj,apj);CHKERRQ(ierr); 473 474 /* Compute P^T*A*P using outer product (P^T)[:,j]*(A*P)[j,:]. */ 475 pnzi = pi[i+1] - pi[i]; 476 for (j=0;j<pnzi;j++) { 477 int nextap=0; 478 crow = *pJ++; 479 cnzj = ci[crow+1] - ci[crow]; 480 cjj = cj + ci[crow]; 481 caj = ca + ci[crow]; 482 /* Perform the sparse axpy operation. Note cjj includes apj. */ 483 for (k=0;k<cnzj;k++) { 484 if (cjj[k]==apj[nextap]) { 485 caj[k] += (*pA)*apa[apj[nextap++]]; 486 } 487 } 488 flops += 2*apnzj; 489 pA++; 490 } 491 492 for (j=0;j<apnzj;j++) { 493 CHKMEMQ; 494 apa[apj[j]] = 0.; 495 CHKMEMQ; 496 } 497 } 498 ierr = PetscFree(apa);CHKERRQ(ierr); 499 ierr = PetscLogFlops(flops);CHKERRQ(ierr); 500 ierr = PetscLogEventEnd(logkey_matapplyptap_numeric,A,P,C,0);CHKERRQ(ierr); 501 PetscFunctionReturn(0); 502 } 503 504 #undef __FUNCT__ 505 #define __FUNCT__ "MatApplyPtAP_SeqAIJ" 506 int MatApplyPtAP_SeqAIJ(Mat A,Mat P,Mat *C) { 507 int ierr; 508 509 PetscFunctionBegin; 510 ierr = MatApplyPtAP_SeqAIJ_Symbolic(A,P,C);CHKERRQ(ierr); 511 ierr = MatApplyPtAP_SeqAIJ_Numeric(A,P,*C);CHKERRQ(ierr); 512 PetscFunctionReturn(0); 513 } 514