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 PetscErrorCode MatMatMatMultSymbolic_SeqAIJ_SeqAIJ_SeqAIJ(Mat A,Mat B,Mat C,PetscReal fill,Mat D) 21 { 22 PetscErrorCode ierr; 23 Mat BC; 24 Mat_MatMatMatMult *matmatmatmult; 25 Mat_SeqAIJ *d; 26 Mat_Product *product = D->product; 27 MatProductAlgorithm alg=product->alg; 28 29 PetscFunctionBegin; 30 if (!product) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NULL,"Data struc Mat_Product is not created, call MatProductCreate() first"); 31 ierr = MatCreate(PETSC_COMM_SELF,&BC);CHKERRQ(ierr); 32 ierr = MatMatMultSymbolic_SeqAIJ_SeqAIJ(B,C,fill,BC);CHKERRQ(ierr); 33 34 ierr = MatProductSetAlgorithm(D,"sorted");CHKERRQ(ierr); /* set alg for D = A*BC */ 35 ierr = MatMatMultSymbolic_SeqAIJ_SeqAIJ(A,BC,fill,D);CHKERRQ(ierr); 36 D->product->alg = alg; /* resume original algorithm for D */ 37 38 /* create struct Mat_MatMatMatMult and attached it to D */ 39 ierr = PetscNew(&matmatmatmult);CHKERRQ(ierr); 40 41 matmatmatmult->BC = BC; 42 matmatmatmult->destroy = D->ops->destroy; 43 d = (Mat_SeqAIJ*)D->data; 44 d->matmatmatmult = matmatmatmult; 45 46 D->ops->matmatmultnumeric = MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqAIJ; 47 D->ops->destroy = MatDestroy_SeqAIJ_MatMatMatMult; 48 PetscFunctionReturn(0); 49 } 50 51 PetscErrorCode MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqAIJ(Mat A,Mat B,Mat C,Mat D) 52 { 53 PetscErrorCode ierr; 54 Mat_SeqAIJ *d =(Mat_SeqAIJ*)D->data; 55 Mat_MatMatMatMult *matmatmatmult=d->matmatmatmult; 56 Mat BC = matmatmatmult->BC; 57 58 PetscFunctionBegin; 59 ierr = (BC->ops->matmultnumeric)(B,C,BC);CHKERRQ(ierr); 60 ierr = (D->ops->matmultnumeric)(A,BC,D);CHKERRQ(ierr); 61 PetscFunctionReturn(0); 62 } 63 64