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; 125f80ce2aSJacob Faibussowitsch CHKERRQ(MatDestroy(&matmatmatmult->BC)); 135f80ce2aSJacob Faibussowitsch CHKERRQ(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); 25*28b400f6SJacob Faibussowitsch PetscCheck(!D->product->data,PetscObjectComm((PetscObject)D),PETSC_ERR_PLIB,"Product data not empty"); 265f80ce2aSJacob Faibussowitsch CHKERRQ(MatCreate(PETSC_COMM_SELF,&BC)); 275f80ce2aSJacob Faibussowitsch CHKERRQ(MatMatMultSymbolic_SeqAIJ_SeqAIJ(B,C,fill,BC)); 286d0b6147SHong Zhang 295f80ce2aSJacob Faibussowitsch CHKERRQ(PetscStrallocpy(D->product->alg,&alg)); 305f80ce2aSJacob Faibussowitsch CHKERRQ(MatProductSetAlgorithm(D,"sorted")); /* set alg for D = A*BC */ 315f80ce2aSJacob Faibussowitsch CHKERRQ(MatMatMultSymbolic_SeqAIJ_SeqAIJ(A,BC,fill,D)); 325f80ce2aSJacob Faibussowitsch CHKERRQ(MatProductSetAlgorithm(D,alg)); /* resume original algorithm */ 335f80ce2aSJacob Faibussowitsch CHKERRQ(PetscFree(alg)); 344222ddf1SHong Zhang 354222ddf1SHong Zhang /* create struct Mat_MatMatMatMult and attached it to D */ 36*28b400f6SJacob Faibussowitsch PetscCheck(!D->product->data,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not yet coded"); 375f80ce2aSJacob Faibussowitsch CHKERRQ(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); 53*28b400f6SJacob 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; 56*28b400f6SJacob Faibussowitsch PetscCheck(BC,PetscObjectComm((PetscObject)D),PETSC_ERR_PLIB,"Missing BC mat"); 57*28b400f6SJacob Faibussowitsch PetscCheck(BC->ops->matmultnumeric,PetscObjectComm((PetscObject)BC),PETSC_ERR_PLIB,"Missing numeric operation"); 585f80ce2aSJacob Faibussowitsch CHKERRQ((*BC->ops->matmultnumeric)(B,C,BC)); 59*28b400f6SJacob Faibussowitsch PetscCheck(D->ops->matmultnumeric,PetscObjectComm((PetscObject)D),PETSC_ERR_PLIB,"Missing numeric operation"); 605f80ce2aSJacob Faibussowitsch CHKERRQ((*D->ops->matmultnumeric)(A,BC,D)); 6171063593SHong Zhang PetscFunctionReturn(0); 6271063593SHong Zhang } 63