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