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 PetscErrorCode ierr; 116d0b6147SHong Zhang 126d0b6147SHong Zhang PetscFunctionBegin; 136d0b6147SHong Zhang ierr = MatDestroy(&matmatmatmult->BC);CHKERRQ(ierr); 146d0b6147SHong Zhang ierr = PetscFree(matmatmatmult);CHKERRQ(ierr); 156d0b6147SHong Zhang PetscFunctionReturn(0); 166d0b6147SHong Zhang } 1771063593SHong Zhang 184222ddf1SHong Zhang PetscErrorCode MatMatMatMultSymbolic_SeqAIJ_SeqAIJ_SeqAIJ(Mat A,Mat B,Mat C,PetscReal fill,Mat D) 1971063593SHong Zhang { 2071063593SHong Zhang PetscErrorCode ierr; 216d0b6147SHong Zhang Mat BC; 226d0b6147SHong Zhang Mat_MatMatMatMult *matmatmatmult; 236718818eSStefano Zampini char *alg; 2471063593SHong Zhang 2571063593SHong Zhang PetscFunctionBegin; 266718818eSStefano Zampini MatCheckProduct(D,5); 27*2c71b3e2SJacob Faibussowitsch PetscCheckFalse(D->product->data,PetscObjectComm((PetscObject)D),PETSC_ERR_PLIB,"Product data not empty"); 284222ddf1SHong Zhang ierr = MatCreate(PETSC_COMM_SELF,&BC);CHKERRQ(ierr); 294222ddf1SHong Zhang ierr = MatMatMultSymbolic_SeqAIJ_SeqAIJ(B,C,fill,BC);CHKERRQ(ierr); 306d0b6147SHong Zhang 316718818eSStefano Zampini ierr = PetscStrallocpy(D->product->alg,&alg);CHKERRQ(ierr); 324222ddf1SHong Zhang ierr = MatProductSetAlgorithm(D,"sorted");CHKERRQ(ierr); /* set alg for D = A*BC */ 334222ddf1SHong Zhang ierr = MatMatMultSymbolic_SeqAIJ_SeqAIJ(A,BC,fill,D);CHKERRQ(ierr); 347a3c3d58SStefano Zampini ierr = MatProductSetAlgorithm(D,alg);CHKERRQ(ierr); /* resume original algorithm */ 356718818eSStefano Zampini ierr = PetscFree(alg);CHKERRQ(ierr); 364222ddf1SHong Zhang 374222ddf1SHong Zhang /* create struct Mat_MatMatMatMult and attached it to D */ 38*2c71b3e2SJacob Faibussowitsch PetscCheckFalse(D->product->data,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not yet coded"); 39b00a9115SJed Brown ierr = PetscNew(&matmatmatmult);CHKERRQ(ierr); 406d0b6147SHong Zhang matmatmatmult->BC = BC; 416718818eSStefano Zampini D->product->data = matmatmatmult; 426718818eSStefano Zampini D->product->destroy = MatDestroy_SeqAIJ_MatMatMatMult; 436d0b6147SHong Zhang 444222ddf1SHong Zhang D->ops->matmatmultnumeric = MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqAIJ; 4571063593SHong Zhang PetscFunctionReturn(0); 4671063593SHong Zhang } 4771063593SHong Zhang 4871063593SHong Zhang PetscErrorCode MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqAIJ(Mat A,Mat B,Mat C,Mat D) 4971063593SHong Zhang { 5071063593SHong Zhang PetscErrorCode ierr; 516718818eSStefano Zampini Mat_MatMatMatMult *matmatmatmult; 526718818eSStefano Zampini Mat BC; 5371063593SHong Zhang 5471063593SHong Zhang PetscFunctionBegin; 556718818eSStefano Zampini MatCheckProduct(D,4); 56*2c71b3e2SJacob Faibussowitsch PetscCheckFalse(!D->product->data,PetscObjectComm((PetscObject)D),PETSC_ERR_PLIB,"Product data empty"); 576718818eSStefano Zampini matmatmatmult = (Mat_MatMatMatMult*)D->product->data; 586718818eSStefano Zampini BC = matmatmatmult->BC; 59*2c71b3e2SJacob Faibussowitsch PetscCheckFalse(!BC,PetscObjectComm((PetscObject)D),PETSC_ERR_PLIB,"Missing BC mat"); 60*2c71b3e2SJacob Faibussowitsch PetscCheckFalse(!BC->ops->matmultnumeric,PetscObjectComm((PetscObject)BC),PETSC_ERR_PLIB,"Missing numeric operation"); 616718818eSStefano Zampini ierr = (*BC->ops->matmultnumeric)(B,C,BC);CHKERRQ(ierr); 62*2c71b3e2SJacob Faibussowitsch PetscCheckFalse(!D->ops->matmultnumeric,PetscObjectComm((PetscObject)D),PETSC_ERR_PLIB,"Missing numeric operation"); 636718818eSStefano Zampini ierr = (*D->ops->matmultnumeric)(A,BC,D);CHKERRQ(ierr); 6471063593SHong Zhang PetscFunctionReturn(0); 6571063593SHong Zhang } 66