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