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*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatDestroy(&matmatmatmult->BC)); 13*5f80ce2aSJacob 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); 252c71b3e2SJacob Faibussowitsch PetscCheckFalse(D->product->data,PetscObjectComm((PetscObject)D),PETSC_ERR_PLIB,"Product data not empty"); 26*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatCreate(PETSC_COMM_SELF,&BC)); 27*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatMatMultSymbolic_SeqAIJ_SeqAIJ(B,C,fill,BC)); 286d0b6147SHong Zhang 29*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscStrallocpy(D->product->alg,&alg)); 30*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatProductSetAlgorithm(D,"sorted")); /* set alg for D = A*BC */ 31*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatMatMultSymbolic_SeqAIJ_SeqAIJ(A,BC,fill,D)); 32*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatProductSetAlgorithm(D,alg)); /* resume original algorithm */ 33*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscFree(alg)); 344222ddf1SHong Zhang 354222ddf1SHong Zhang /* create struct Mat_MatMatMatMult and attached it to D */ 362c71b3e2SJacob Faibussowitsch PetscCheckFalse(D->product->data,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not yet coded"); 37*5f80ce2aSJacob 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); 532c71b3e2SJacob Faibussowitsch PetscCheckFalse(!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; 562c71b3e2SJacob Faibussowitsch PetscCheckFalse(!BC,PetscObjectComm((PetscObject)D),PETSC_ERR_PLIB,"Missing BC mat"); 572c71b3e2SJacob Faibussowitsch PetscCheckFalse(!BC->ops->matmultnumeric,PetscObjectComm((PetscObject)BC),PETSC_ERR_PLIB,"Missing numeric operation"); 58*5f80ce2aSJacob Faibussowitsch CHKERRQ((*BC->ops->matmultnumeric)(B,C,BC)); 592c71b3e2SJacob Faibussowitsch PetscCheckFalse(!D->ops->matmultnumeric,PetscObjectComm((PetscObject)D),PETSC_ERR_PLIB,"Missing numeric operation"); 60*5f80ce2aSJacob Faibussowitsch CHKERRQ((*D->ops->matmultnumeric)(A,BC,D)); 6171063593SHong Zhang PetscFunctionReturn(0); 6271063593SHong Zhang } 63