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*d71ae5a4SJacob Faibussowitsch PetscErrorCode MatDestroy_SeqAIJ_MatMatMatMult(void *data) 8*d71ae5a4SJacob Faibussowitsch { 96718818eSStefano Zampini Mat_MatMatMatMult *matmatmatmult = (Mat_MatMatMatMult *)data; 106d0b6147SHong Zhang 116d0b6147SHong Zhang PetscFunctionBegin; 129566063dSJacob Faibussowitsch PetscCall(MatDestroy(&matmatmatmult->BC)); 139566063dSJacob Faibussowitsch PetscCall(PetscFree(matmatmatmult)); 146d0b6147SHong Zhang PetscFunctionReturn(0); 156d0b6147SHong Zhang } 1671063593SHong Zhang 17*d71ae5a4SJacob Faibussowitsch PetscErrorCode MatMatMatMultSymbolic_SeqAIJ_SeqAIJ_SeqAIJ(Mat A, Mat B, Mat C, PetscReal fill, Mat D) 18*d71ae5a4SJacob Faibussowitsch { 196d0b6147SHong Zhang Mat BC; 206d0b6147SHong Zhang Mat_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 354222ddf1SHong Zhang /* create struct Mat_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; 406718818eSStefano Zampini D->product->destroy = MatDestroy_SeqAIJ_MatMatMatMult; 416d0b6147SHong Zhang 424222ddf1SHong Zhang D->ops->matmatmultnumeric = MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqAIJ; 4371063593SHong Zhang PetscFunctionReturn(0); 4471063593SHong Zhang } 4571063593SHong Zhang 46*d71ae5a4SJacob Faibussowitsch PetscErrorCode MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqAIJ(Mat A, Mat B, Mat C, Mat D) 47*d71ae5a4SJacob Faibussowitsch { 486718818eSStefano Zampini Mat_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"); 546718818eSStefano Zampini matmatmatmult = (Mat_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)); 5971063593SHong Zhang PetscFunctionReturn(0); 6071063593SHong Zhang } 61