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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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 PetscCall((*BC->ops->matmultnumeric)(B, C, BC)); 58 PetscCall((*D->ops->matmultnumeric)(A, BC, D)); 59 PetscFunctionReturn(PETSC_SUCCESS); 60 } 61