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 #if defined(PETSC_HAVE_HYPRE) 23 PETSC_INTERN PetscErrorCode MatMatMatMultSymbolic_AIJ_AIJ_AIJ_wHYPRE(Mat,Mat,Mat,PetscReal,Mat*); 24 #endif 25 26 #undef __FUNCT__ 27 #define __FUNCT__ "MatMatMatMult_SeqAIJ_SeqAIJ_SeqAIJ" 28 PETSC_INTERN PetscErrorCode MatMatMatMult_SeqAIJ_SeqAIJ_SeqAIJ(Mat A,Mat B,Mat C,MatReuse scall,PetscReal fill,Mat *D) 29 { 30 PetscErrorCode ierr; 31 32 PetscFunctionBegin; 33 if (scall == MAT_INITIAL_MATRIX) { 34 #if defined(PETSC_HAVE_HYPRE) 35 PetscBool hypre = PETSC_FALSE; 36 37 ierr = PetscObjectOptionsBegin((PetscObject)B);CHKERRQ(ierr); 38 ierr = PetscOptionsBool("-matmatmatmult_viahypre","Use HYPRE to perform D=A*B*C","",hypre,&hypre,NULL);CHKERRQ(ierr); 39 ierr = PetscOptionsEnd();CHKERRQ(ierr); 40 #endif 41 ierr = PetscLogEventBegin(MAT_MatMatMultSymbolic,A,B,C,0);CHKERRQ(ierr); 42 #if !defined(PETSC_HAVE_HYPRE) 43 ierr = MatMatMatMultSymbolic_SeqAIJ_SeqAIJ_SeqAIJ(A,B,C,fill,D);CHKERRQ(ierr); 44 #else 45 if (hypre) { 46 ierr = MatMatMatMultSymbolic_AIJ_AIJ_AIJ_wHYPRE(A,B,C,fill,D);CHKERRQ(ierr); 47 } else { 48 ierr = MatMatMatMultSymbolic_SeqAIJ_SeqAIJ_SeqAIJ(A,B,C,fill,D);CHKERRQ(ierr); 49 } 50 #endif 51 ierr = PetscLogEventEnd(MAT_MatMatMultSymbolic,A,B,C,0);CHKERRQ(ierr); 52 } 53 ierr = PetscLogEventBegin(MAT_MatMatMultNumeric,A,B,C,0);CHKERRQ(ierr); 54 ierr = ((*D)->ops->matmatmultnumeric)(A,B,C,*D);CHKERRQ(ierr); 55 ierr = PetscLogEventEnd(MAT_MatMatMultNumeric,A,B,C,0);CHKERRQ(ierr); 56 PetscFunctionReturn(0); 57 } 58 59 #undef __FUNCT__ 60 #define __FUNCT__ "MatMatMatMultSymbolic_SeqAIJ_SeqAIJ_SeqAIJ" 61 PetscErrorCode MatMatMatMultSymbolic_SeqAIJ_SeqAIJ_SeqAIJ(Mat A,Mat B,Mat C,PetscReal fill,Mat *D) 62 { 63 PetscErrorCode ierr; 64 Mat BC; 65 Mat_MatMatMatMult *matmatmatmult; 66 Mat_SeqAIJ *d; 67 PetscBool scalable=PETSC_TRUE; 68 69 PetscFunctionBegin; 70 ierr = PetscObjectOptionsBegin((PetscObject)B);CHKERRQ(ierr); 71 ierr = PetscOptionsBool("-matmatmatmult_scalable","Use a scalable but slower D=A*B*C","",scalable,&scalable,NULL);CHKERRQ(ierr); 72 ierr = PetscOptionsEnd();CHKERRQ(ierr); 73 if (scalable) { 74 ierr = MatMatMultSymbolic_SeqAIJ_SeqAIJ_Scalable(B,C,fill,&BC);CHKERRQ(ierr); 75 ierr = MatMatMultSymbolic_SeqAIJ_SeqAIJ_Scalable(A,BC,fill,D);CHKERRQ(ierr); 76 } else { 77 ierr = MatMatMultSymbolic_SeqAIJ_SeqAIJ(B,C,fill,&BC);CHKERRQ(ierr); 78 ierr = MatMatMultSymbolic_SeqAIJ_SeqAIJ(A,BC,fill,D);CHKERRQ(ierr); 79 } 80 81 /* create struct Mat_MatMatMatMult and attached it to *D */ 82 ierr = PetscNew(&matmatmatmult);CHKERRQ(ierr); 83 84 matmatmatmult->BC = BC; 85 matmatmatmult->destroy = (*D)->ops->destroy; 86 d = (Mat_SeqAIJ*)(*D)->data; 87 d->matmatmatmult = matmatmatmult; 88 89 (*D)->ops->matmatmultnumeric = MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqAIJ; 90 (*D)->ops->destroy = MatDestroy_SeqAIJ_MatMatMatMult; 91 PetscFunctionReturn(0); 92 } 93 94 #undef __FUNCT__ 95 #define __FUNCT__ "MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqAIJ" 96 PetscErrorCode MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqAIJ(Mat A,Mat B,Mat C,Mat D) 97 { 98 PetscErrorCode ierr; 99 Mat_SeqAIJ *d =(Mat_SeqAIJ*)D->data; 100 Mat_MatMatMatMult *matmatmatmult=d->matmatmatmult; 101 Mat BC = matmatmatmult->BC; 102 103 PetscFunctionBegin; 104 ierr = (BC->ops->matmultnumeric)(B,C,BC);CHKERRQ(ierr); 105 ierr = (D->ops->matmultnumeric)(A,BC,D);CHKERRQ(ierr); 106 PetscFunctionReturn(0); 107 } 108