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(void* data) 8 { 9 Mat_MatMatMatMult *matmatmatmult = (Mat_MatMatMatMult*)data; 10 PetscErrorCode ierr; 11 12 PetscFunctionBegin; 13 ierr = MatDestroy(&matmatmatmult->BC);CHKERRQ(ierr); 14 ierr = PetscFree(matmatmatmult);CHKERRQ(ierr); 15 PetscFunctionReturn(0); 16 } 17 18 PetscErrorCode MatMatMatMultSymbolic_SeqAIJ_SeqAIJ_SeqAIJ(Mat A,Mat B,Mat C,PetscReal fill,Mat D) 19 { 20 PetscErrorCode ierr; 21 Mat BC; 22 Mat_MatMatMatMult *matmatmatmult; 23 char *alg; 24 25 PetscFunctionBegin; 26 MatCheckProduct(D,5); 27 if (D->product->data) SETERRQ(PetscObjectComm((PetscObject)D),PETSC_ERR_PLIB,"Product data not empty"); 28 ierr = MatCreate(PETSC_COMM_SELF,&BC);CHKERRQ(ierr); 29 ierr = MatMatMultSymbolic_SeqAIJ_SeqAIJ(B,C,fill,BC);CHKERRQ(ierr); 30 31 ierr = PetscStrallocpy(D->product->alg,&alg);CHKERRQ(ierr); 32 ierr = MatProductSetAlgorithm(D,"sorted");CHKERRQ(ierr); /* set alg for D = A*BC */ 33 ierr = MatMatMultSymbolic_SeqAIJ_SeqAIJ(A,BC,fill,D);CHKERRQ(ierr); 34 ierr = MatProductSetAlgorithm(D,alg);CHKERRQ(ierr); /* resume original algorithm */ 35 ierr = PetscFree(alg);CHKERRQ(ierr); 36 37 /* create struct Mat_MatMatMatMult and attached it to D */ 38 if (D->product->data) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not yet coded"); 39 ierr = PetscNew(&matmatmatmult);CHKERRQ(ierr); 40 matmatmatmult->BC = BC; 41 D->product->data = matmatmatmult; 42 D->product->destroy = MatDestroy_SeqAIJ_MatMatMatMult; 43 44 D->ops->matmatmultnumeric = MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqAIJ; 45 PetscFunctionReturn(0); 46 } 47 48 PetscErrorCode MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqAIJ(Mat A,Mat B,Mat C,Mat D) 49 { 50 PetscErrorCode ierr; 51 Mat_MatMatMatMult *matmatmatmult; 52 Mat BC; 53 54 PetscFunctionBegin; 55 MatCheckProduct(D,4); 56 if (!D->product->data) SETERRQ(PetscObjectComm((PetscObject)D),PETSC_ERR_PLIB,"Product data empty"); 57 matmatmatmult = (Mat_MatMatMatMult*)D->product->data; 58 BC = matmatmatmult->BC; 59 if (!BC) SETERRQ(PetscObjectComm((PetscObject)D),PETSC_ERR_PLIB,"Missing BC mat"); 60 if (!BC->ops->matmultnumeric) SETERRQ(PetscObjectComm((PetscObject)BC),PETSC_ERR_PLIB,"Missing numeric operation"); 61 ierr = (*BC->ops->matmultnumeric)(B,C,BC);CHKERRQ(ierr); 62 if (!D->ops->matmultnumeric) SETERRQ(PetscObjectComm((PetscObject)D),PETSC_ERR_PLIB,"Missing numeric operation"); 63 ierr = (*D->ops->matmultnumeric)(A,BC,D);CHKERRQ(ierr); 64 PetscFunctionReturn(0); 65 } 66