1 /* 2 Defines matrix-matrix-matrix product routines for SeqAIJ matrices 3 D = A * B * C 4 */ 5 #include <../src/mat/impls/aij/seq/aij.h> /*I "petscmat.h" I*/ 6 7 #undef __FUNCT__ 8 #define __FUNCT__ "MatDestroy_SeqAIJ_MatMatMatMult" 9 PetscErrorCode MatDestroy_SeqAIJ_MatMatMatMult(Mat A) 10 { 11 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 12 Mat_MatMatMatMult *matmatmatmult=a->matmatmatmult; 13 PetscErrorCode ierr; 14 15 PetscFunctionBegin; 16 ierr = MatDestroy(&matmatmatmult->BC);CHKERRQ(ierr); 17 ierr = matmatmatmult->destroy(A);CHKERRQ(ierr); 18 ierr = PetscFree(matmatmatmult);CHKERRQ(ierr); 19 PetscFunctionReturn(0); 20 } 21 22 #undef __FUNCT__ 23 #define __FUNCT__ "MatMatMatMult_SeqAIJ_SeqAIJ_SeqAIJ" 24 PetscErrorCode MatMatMatMult_SeqAIJ_SeqAIJ_SeqAIJ(Mat A,Mat B,Mat C,MatReuse scall,PetscReal fill,Mat *D) 25 { 26 PetscErrorCode ierr; 27 28 PetscFunctionBegin; 29 if (scall == MAT_INITIAL_MATRIX){ 30 ierr = PetscLogEventBegin(MAT_MatMatMultSymbolic,A,B,C,0);CHKERRQ(ierr); 31 ierr = MatMatMatMultSymbolic_SeqAIJ_SeqAIJ_SeqAIJ(A,B,C,fill,D);CHKERRQ(ierr); 32 ierr = PetscLogEventEnd(MAT_MatMatMultSymbolic,A,B,C,0);CHKERRQ(ierr); 33 } 34 ierr = PetscLogEventBegin(MAT_MatMatMultNumeric,A,B,C,0);CHKERRQ(ierr); 35 ierr = MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqAIJ(A,B,C,*D);CHKERRQ(ierr); 36 ierr = PetscLogEventEnd(MAT_MatMatMultNumeric,A,B,C,0);CHKERRQ(ierr); 37 PetscFunctionReturn(0); 38 } 39 40 #undef __FUNCT__ 41 #define __FUNCT__ "MatMatMatMultSymbolic_SeqAIJ_SeqAIJ_SeqAIJ" 42 PetscErrorCode MatMatMatMultSymbolic_SeqAIJ_SeqAIJ_SeqAIJ(Mat A,Mat B,Mat C,PetscReal fill,Mat *D) 43 { 44 PetscErrorCode ierr; 45 Mat BC; 46 Mat_MatMatMatMult *matmatmatmult; 47 Mat_SeqAIJ *d; 48 PetscBool scalable=PETSC_TRUE; 49 PetscLogDouble t0,t1,t2; 50 51 PetscFunctionBegin; 52 ierr = PetscObjectOptionsBegin((PetscObject)B);CHKERRQ(ierr); 53 ierr = PetscOptionsBool("-matmatmatmult_scalable","Use a scalable but slower D=A*B*C","",scalable,&scalable,PETSC_NULL);CHKERRQ(ierr); 54 ierr = PetscOptionsEnd();CHKERRQ(ierr); 55 if (scalable){ 56 ierr = PetscGetTime(&t0);CHKERRQ(ierr); 57 ierr = MatMatMultSymbolic_SeqAIJ_SeqAIJ_Scalable(B,C,fill,&BC);CHKERRQ(ierr); 58 ierr = PetscGetTime(&t1);CHKERRQ(ierr); 59 ierr = MatMatMultSymbolic_SeqAIJ_SeqAIJ_Scalable(A,BC,fill,D);CHKERRQ(ierr); 60 ierr = PetscGetTime(&t2);CHKERRQ(ierr); 61 printf(" Mat %d %d, 3MultSymbolic_SeqAIJ_Scalable time: %g + %g = %g\n",A->rmap->N,A->cmap->N,t1-t0,t2-t1,t2-t0); 62 } else { 63 ierr = PetscGetTime(&t0);CHKERRQ(ierr); 64 ierr = MatMatMultSymbolic_SeqAIJ_SeqAIJ(B,C,fill,&BC);CHKERRQ(ierr); 65 ierr = PetscGetTime(&t1);CHKERRQ(ierr); 66 ierr = MatMatMultSymbolic_SeqAIJ_SeqAIJ(A,BC,fill,D);CHKERRQ(ierr); 67 ierr = PetscGetTime(&t2);CHKERRQ(ierr); 68 printf(" Mat %d %d, 3MultSymbolic_SeqAIJ time: %g + %g = %g\n",A->rmap->N,A->cmap->N,t1-t0,t2-t1,t2-t0); 69 } 70 71 /* create struct Mat_MatMatMatMult and attached it to *D */ 72 ierr = PetscNew(Mat_MatMatMatMult,&matmatmatmult);CHKERRQ(ierr); 73 matmatmatmult->BC = BC; 74 matmatmatmult->destroy = (*D)->ops->destroy; 75 d = (Mat_SeqAIJ*)(*D)->data; 76 d->matmatmatmult = matmatmatmult; 77 78 (*D)->ops->matmatmultnumeric = MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqAIJ; 79 (*D)->ops->destroy = MatDestroy_SeqAIJ_MatMatMatMult; 80 PetscFunctionReturn(0); 81 } 82 83 #undef __FUNCT__ 84 #define __FUNCT__ "MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqAIJ" 85 PetscErrorCode MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqAIJ(Mat A,Mat B,Mat C,Mat D) 86 { 87 PetscErrorCode ierr; 88 Mat_SeqAIJ *d=(Mat_SeqAIJ*)D->data; 89 Mat_MatMatMatMult *matmatmatmult=d->matmatmatmult; 90 Mat BC= matmatmatmult->BC; 91 PetscLogDouble t0,t1,t2; 92 93 PetscFunctionBegin; 94 ierr = PetscGetTime(&t0);CHKERRQ(ierr); 95 ierr = (BC->ops->matmultnumeric)(B,C,BC);CHKERRQ(ierr); 96 ierr = PetscGetTime(&t1);CHKERRQ(ierr); 97 ierr = (D->ops->matmultnumeric)(A,BC,D);CHKERRQ(ierr); 98 ierr = PetscGetTime(&t2);CHKERRQ(ierr); 99 printf(" 3MultNumeric_SeqAIJ time: %g + %g = %g\n",t1-t0,t2-t1,t2-t0); 100 PetscFunctionReturn(0); 101 } 102