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