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 PetscErrorCode MatDestroy_SeqAIJ_MatMatMatMult(Mat A) 8 { 9 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 10 Mat_MatMatMatMult *matmatmatmult=a->matmatmatmult; 11 PetscErrorCode ierr; 12 13 PetscFunctionBegin; 14 ierr = MatDestroy(&matmatmatmult->BC);CHKERRQ(ierr); 15 ierr = matmatmatmult->destroy(A);CHKERRQ(ierr); 16 ierr = PetscFree(matmatmatmult);CHKERRQ(ierr); 17 PetscFunctionReturn(0); 18 } 19 20 PETSC_INTERN PetscErrorCode MatMatMatMult_SeqAIJ_SeqAIJ_SeqAIJ(Mat A,Mat B,Mat C,MatReuse scall,PetscReal fill,Mat *D) 21 { 22 PetscErrorCode ierr; 23 24 PetscFunctionBegin; 25 if (scall == MAT_INITIAL_MATRIX) { 26 ierr = PetscLogEventBegin(MAT_MatMatMultSymbolic,A,B,C,0);CHKERRQ(ierr); 27 ierr = MatMatMatMultSymbolic_SeqAIJ_SeqAIJ_SeqAIJ(A,B,C,fill,D);CHKERRQ(ierr); 28 ierr = PetscLogEventEnd(MAT_MatMatMultSymbolic,A,B,C,0);CHKERRQ(ierr); 29 } 30 ierr = PetscLogEventBegin(MAT_MatMatMultNumeric,A,B,C,0);CHKERRQ(ierr); 31 ierr = ((*D)->ops->matmatmultnumeric)(A,B,C,*D);CHKERRQ(ierr); 32 ierr = PetscLogEventEnd(MAT_MatMatMultNumeric,A,B,C,0);CHKERRQ(ierr); 33 PetscFunctionReturn(0); 34 } 35 36 PetscErrorCode MatMatMatMultSymbolic_SeqAIJ_SeqAIJ_SeqAIJ(Mat A,Mat B,Mat C,PetscReal fill,Mat *D) 37 { 38 PetscErrorCode ierr; 39 Mat BC; 40 Mat_MatMatMatMult *matmatmatmult; 41 Mat_SeqAIJ *d; 42 PetscBool scalable=PETSC_FALSE; 43 44 PetscFunctionBegin; 45 ierr = PetscObjectOptionsBegin((PetscObject)B);CHKERRQ(ierr); 46 ierr = PetscOptionsBool("-matmatmatmult_scalable","Use a scalable but slower D=A*B*C","",scalable,&scalable,NULL);CHKERRQ(ierr); 47 ierr = PetscOptionsEnd();CHKERRQ(ierr); 48 if (scalable) { 49 ierr = MatMatMultSymbolic_SeqAIJ_SeqAIJ_Scalable(B,C,fill,&BC);CHKERRQ(ierr); 50 ierr = MatMatMultSymbolic_SeqAIJ_SeqAIJ_Scalable(A,BC,fill,D);CHKERRQ(ierr); 51 } else { 52 ierr = MatMatMultSymbolic_SeqAIJ_SeqAIJ(B,C,fill,&BC);CHKERRQ(ierr); 53 ierr = MatMatMultSymbolic_SeqAIJ_SeqAIJ(A,BC,fill,D);CHKERRQ(ierr); 54 } 55 56 /* create struct Mat_MatMatMatMult and attached it to *D */ 57 ierr = PetscNew(&matmatmatmult);CHKERRQ(ierr); 58 59 matmatmatmult->BC = BC; 60 matmatmatmult->destroy = (*D)->ops->destroy; 61 d = (Mat_SeqAIJ*)(*D)->data; 62 d->matmatmatmult = matmatmatmult; 63 64 (*D)->ops->matmatmultnumeric = MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqAIJ; 65 (*D)->ops->destroy = MatDestroy_SeqAIJ_MatMatMatMult; 66 PetscFunctionReturn(0); 67 } 68 69 PetscErrorCode MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqAIJ(Mat A,Mat B,Mat C,Mat D) 70 { 71 PetscErrorCode ierr; 72 Mat_SeqAIJ *d =(Mat_SeqAIJ*)D->data; 73 Mat_MatMatMatMult *matmatmatmult=d->matmatmatmult; 74 Mat BC = matmatmatmult->BC; 75 76 PetscFunctionBegin; 77 ierr = (BC->ops->matmultnumeric)(B,C,BC);CHKERRQ(ierr); 78 ierr = (D->ops->matmultnumeric)(A,BC,D);CHKERRQ(ierr); 79 PetscFunctionReturn(0); 80 } 81