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 76718818eSStefano Zampini PetscErrorCode MatDestroy_SeqAIJ_MatMatMatMult(void* data) 86d0b6147SHong Zhang { 96718818eSStefano Zampini Mat_MatMatMatMult *matmatmatmult = (Mat_MatMatMatMult*)data; 106d0b6147SHong Zhang 116d0b6147SHong Zhang PetscFunctionBegin; 12*9566063dSJacob Faibussowitsch PetscCall(MatDestroy(&matmatmatmult->BC)); 13*9566063dSJacob Faibussowitsch PetscCall(PetscFree(matmatmatmult)); 146d0b6147SHong Zhang PetscFunctionReturn(0); 156d0b6147SHong Zhang } 1671063593SHong Zhang 174222ddf1SHong Zhang PetscErrorCode MatMatMatMultSymbolic_SeqAIJ_SeqAIJ_SeqAIJ(Mat A,Mat B,Mat C,PetscReal fill,Mat D) 1871063593SHong Zhang { 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"); 26*9566063dSJacob Faibussowitsch PetscCall(MatCreate(PETSC_COMM_SELF,&BC)); 27*9566063dSJacob Faibussowitsch PetscCall(MatMatMultSymbolic_SeqAIJ_SeqAIJ(B,C,fill,BC)); 286d0b6147SHong Zhang 29*9566063dSJacob Faibussowitsch PetscCall(PetscStrallocpy(D->product->alg,&alg)); 30*9566063dSJacob Faibussowitsch PetscCall(MatProductSetAlgorithm(D,"sorted")); /* set alg for D = A*BC */ 31*9566063dSJacob Faibussowitsch PetscCall(MatMatMultSymbolic_SeqAIJ_SeqAIJ(A,BC,fill,D)); 32*9566063dSJacob Faibussowitsch PetscCall(MatProductSetAlgorithm(D,alg)); /* resume original algorithm */ 33*9566063dSJacob 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"); 37*9566063dSJacob 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 4671063593SHong Zhang PetscErrorCode MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqAIJ(Mat A,Mat B,Mat C,Mat D) 4771063593SHong Zhang { 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"); 5728b400f6SJacob Faibussowitsch PetscCheck(BC->ops->matmultnumeric,PetscObjectComm((PetscObject)BC),PETSC_ERR_PLIB,"Missing numeric operation"); 58*9566063dSJacob Faibussowitsch PetscCall((*BC->ops->matmultnumeric)(B,C,BC)); 5928b400f6SJacob Faibussowitsch PetscCheck(D->ops->matmultnumeric,PetscObjectComm((PetscObject)D),PETSC_ERR_PLIB,"Missing numeric operation"); 60*9566063dSJacob Faibussowitsch PetscCall((*D->ops->matmultnumeric)(A,BC,D)); 6171063593SHong Zhang PetscFunctionReturn(0); 6271063593SHong Zhang } 63