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*6718818eSStefano Zampini PetscErrorCode MatDestroy_SeqAIJ_MatMatMatMult(void* data) 86d0b6147SHong Zhang { 9*6718818eSStefano 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; 23*6718818eSStefano Zampini char *alg; 2471063593SHong Zhang 2571063593SHong Zhang PetscFunctionBegin; 26*6718818eSStefano Zampini MatCheckProduct(D,5); 27*6718818eSStefano Zampini if (D->product->data) SETERRQ(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 31*6718818eSStefano 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 */ 35*6718818eSStefano Zampini ierr = PetscFree(alg);CHKERRQ(ierr); 364222ddf1SHong Zhang 374222ddf1SHong Zhang /* create struct Mat_MatMatMatMult and attached it to D */ 38*6718818eSStefano Zampini if (D->product->data) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not yet coded"); 39b00a9115SJed Brown ierr = PetscNew(&matmatmatmult);CHKERRQ(ierr); 406d0b6147SHong Zhang matmatmatmult->BC = BC; 41*6718818eSStefano Zampini D->product->data = matmatmatmult; 42*6718818eSStefano 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; 51*6718818eSStefano Zampini Mat_MatMatMatMult *matmatmatmult; 52*6718818eSStefano Zampini Mat BC; 5371063593SHong Zhang 5471063593SHong Zhang PetscFunctionBegin; 55*6718818eSStefano Zampini MatCheckProduct(D,4); 56*6718818eSStefano Zampini if (!D->product->data) SETERRQ(PetscObjectComm((PetscObject)D),PETSC_ERR_PLIB,"Product data empty"); 57*6718818eSStefano Zampini matmatmatmult = (Mat_MatMatMatMult*)D->product->data; 58*6718818eSStefano Zampini BC = matmatmatmult->BC; 59*6718818eSStefano Zampini if (!BC) SETERRQ(PetscObjectComm((PetscObject)D),PETSC_ERR_PLIB,"Missing BC mat"); 60*6718818eSStefano Zampini if (!BC->ops->matmultnumeric) SETERRQ(PetscObjectComm((PetscObject)BC),PETSC_ERR_PLIB,"Missing numeric operation"); 61*6718818eSStefano Zampini ierr = (*BC->ops->matmultnumeric)(B,C,BC);CHKERRQ(ierr); 62*6718818eSStefano Zampini if (!D->ops->matmultnumeric) SETERRQ(PetscObjectComm((PetscObject)D),PETSC_ERR_PLIB,"Missing numeric operation"); 63*6718818eSStefano Zampini ierr = (*D->ops->matmultnumeric)(A,BC,D);CHKERRQ(ierr); 6471063593SHong Zhang PetscFunctionReturn(0); 6571063593SHong Zhang } 66