171063593SHong Zhang /* 26d0b6147SHong Zhang Defines matrix-matrix-matrix product routines for SeqAIJ matrices 371063593SHong Zhang D = A * B * C 471063593SHong Zhang */ 571063593SHong Zhang #include <../src/mat/impls/aij/seq/aij.h> /*I "petscmat.h" I*/ 66d0b6147SHong Zhang 7*cc1eb50dSBarry Smith static PetscErrorCode MatProductCtxDestroy_SeqAIJ_MatMatMatMult(void **data) 8d71ae5a4SJacob Faibussowitsch { 9*cc1eb50dSBarry Smith MatProductCtx_MatMatMatMult *matmatmatmult = *(MatProductCtx_MatMatMatMult **)data; 106d0b6147SHong Zhang 116d0b6147SHong Zhang PetscFunctionBegin; 129566063dSJacob Faibussowitsch PetscCall(MatDestroy(&matmatmatmult->BC)); 139566063dSJacob Faibussowitsch PetscCall(PetscFree(matmatmatmult)); 143ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 156d0b6147SHong Zhang } 1671063593SHong Zhang 17d71ae5a4SJacob Faibussowitsch PetscErrorCode MatMatMatMultSymbolic_SeqAIJ_SeqAIJ_SeqAIJ(Mat A, Mat B, Mat C, PetscReal fill, Mat D) 18d71ae5a4SJacob Faibussowitsch { 196d0b6147SHong Zhang Mat BC; 20*cc1eb50dSBarry Smith MatProductCtx_MatMatMatMult *matmatmatmult; 216718818eSStefano Zampini char *alg; 2271063593SHong Zhang 2371063593SHong Zhang PetscFunctionBegin; 246718818eSStefano Zampini MatCheckProduct(D, 5); 2528b400f6SJacob Faibussowitsch PetscCheck(!D->product->data, PetscObjectComm((PetscObject)D), PETSC_ERR_PLIB, "Product data not empty"); 269566063dSJacob Faibussowitsch PetscCall(MatCreate(PETSC_COMM_SELF, &BC)); 279566063dSJacob Faibussowitsch PetscCall(MatMatMultSymbolic_SeqAIJ_SeqAIJ(B, C, fill, BC)); 286d0b6147SHong Zhang 299566063dSJacob Faibussowitsch PetscCall(PetscStrallocpy(D->product->alg, &alg)); 309566063dSJacob Faibussowitsch PetscCall(MatProductSetAlgorithm(D, "sorted")); /* set alg for D = A*BC */ 319566063dSJacob Faibussowitsch PetscCall(MatMatMultSymbolic_SeqAIJ_SeqAIJ(A, BC, fill, D)); 329566063dSJacob Faibussowitsch PetscCall(MatProductSetAlgorithm(D, alg)); /* resume original algorithm */ 339566063dSJacob Faibussowitsch PetscCall(PetscFree(alg)); 344222ddf1SHong Zhang 35*cc1eb50dSBarry Smith /* create struct MatProductCtx_MatMatMatMult and attached it to D */ 3628b400f6SJacob Faibussowitsch PetscCheck(!D->product->data, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Not yet coded"); 379566063dSJacob Faibussowitsch PetscCall(PetscNew(&matmatmatmult)); 386d0b6147SHong Zhang matmatmatmult->BC = BC; 396718818eSStefano Zampini D->product->data = matmatmatmult; 40*cc1eb50dSBarry Smith D->product->destroy = MatProductCtxDestroy_SeqAIJ_MatMatMatMult; 416d0b6147SHong Zhang 424222ddf1SHong Zhang D->ops->matmatmultnumeric = MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqAIJ; 433ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 4471063593SHong Zhang } 4571063593SHong Zhang 46d71ae5a4SJacob Faibussowitsch PetscErrorCode MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqAIJ(Mat A, Mat B, Mat C, Mat D) 47d71ae5a4SJacob Faibussowitsch { 48*cc1eb50dSBarry Smith MatProductCtx_MatMatMatMult *matmatmatmult; 496718818eSStefano Zampini Mat BC; 5071063593SHong Zhang 5171063593SHong Zhang PetscFunctionBegin; 526718818eSStefano Zampini MatCheckProduct(D, 4); 5328b400f6SJacob Faibussowitsch PetscCheck(D->product->data, PetscObjectComm((PetscObject)D), PETSC_ERR_PLIB, "Product data empty"); 54*cc1eb50dSBarry Smith matmatmatmult = (MatProductCtx_MatMatMatMult *)D->product->data; 556718818eSStefano Zampini BC = matmatmatmult->BC; 5628b400f6SJacob Faibussowitsch PetscCheck(BC, PetscObjectComm((PetscObject)D), PETSC_ERR_PLIB, "Missing BC mat"); 579566063dSJacob Faibussowitsch PetscCall((*BC->ops->matmultnumeric)(B, C, BC)); 589566063dSJacob Faibussowitsch PetscCall((*D->ops->matmultnumeric)(A, BC, D)); 593ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 6071063593SHong Zhang } 61