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 76d0b6147SHong Zhang PetscErrorCode MatDestroy_SeqAIJ_MatMatMatMult(Mat A) 86d0b6147SHong Zhang { 96d0b6147SHong Zhang Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 106d0b6147SHong Zhang Mat_MatMatMatMult *matmatmatmult=a->matmatmatmult; 116d0b6147SHong Zhang PetscErrorCode ierr; 126d0b6147SHong Zhang 136d0b6147SHong Zhang PetscFunctionBegin; 146d0b6147SHong Zhang ierr = MatDestroy(&matmatmatmult->BC);CHKERRQ(ierr); 156d0b6147SHong Zhang ierr = matmatmatmult->destroy(A);CHKERRQ(ierr); 166d0b6147SHong Zhang ierr = PetscFree(matmatmatmult);CHKERRQ(ierr); 176d0b6147SHong Zhang PetscFunctionReturn(0); 186d0b6147SHong Zhang } 1971063593SHong Zhang 204222ddf1SHong Zhang PetscErrorCode MatMatMatMultSymbolic_SeqAIJ_SeqAIJ_SeqAIJ(Mat A,Mat B,Mat C,PetscReal fill,Mat D) 2171063593SHong Zhang { 2271063593SHong Zhang PetscErrorCode ierr; 236d0b6147SHong Zhang Mat BC; 246d0b6147SHong Zhang Mat_MatMatMatMult *matmatmatmult; 256d0b6147SHong Zhang Mat_SeqAIJ *d; 264222ddf1SHong Zhang Mat_Product *product = D->product; 274222ddf1SHong Zhang MatProductAlgorithm alg=product->alg; 2871063593SHong Zhang 2971063593SHong Zhang PetscFunctionBegin; 304222ddf1SHong Zhang if (!product) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NULL,"Data struc Mat_Product is not created, call MatProductCreate() first"); 314222ddf1SHong Zhang ierr = MatCreate(PETSC_COMM_SELF,&BC);CHKERRQ(ierr); 324222ddf1SHong Zhang ierr = MatMatMultSymbolic_SeqAIJ_SeqAIJ(B,C,fill,BC);CHKERRQ(ierr); 336d0b6147SHong Zhang 344222ddf1SHong Zhang ierr = MatProductSetAlgorithm(D,"sorted");CHKERRQ(ierr); /* set alg for D = A*BC */ 354222ddf1SHong Zhang ierr = MatMatMultSymbolic_SeqAIJ_SeqAIJ(A,BC,fill,D);CHKERRQ(ierr); 36*7a3c3d58SStefano Zampini ierr = MatProductSetAlgorithm(D,alg);CHKERRQ(ierr); /* resume original algorithm */ 374222ddf1SHong Zhang 384222ddf1SHong Zhang /* create struct Mat_MatMatMatMult and attached it to D */ 39b00a9115SJed Brown ierr = PetscNew(&matmatmatmult);CHKERRQ(ierr); 402205254eSKarl Rupp 416d0b6147SHong Zhang matmatmatmult->BC = BC; 424222ddf1SHong Zhang matmatmatmult->destroy = D->ops->destroy; 434222ddf1SHong Zhang d = (Mat_SeqAIJ*)D->data; 446d0b6147SHong Zhang d->matmatmatmult = matmatmatmult; 456d0b6147SHong Zhang 464222ddf1SHong Zhang D->ops->matmatmultnumeric = MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqAIJ; 474222ddf1SHong Zhang D->ops->destroy = MatDestroy_SeqAIJ_MatMatMatMult; 4871063593SHong Zhang PetscFunctionReturn(0); 4971063593SHong Zhang } 5071063593SHong Zhang 5171063593SHong Zhang PetscErrorCode MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqAIJ(Mat A,Mat B,Mat C,Mat D) 5271063593SHong Zhang { 5371063593SHong Zhang PetscErrorCode ierr; 546d0b6147SHong Zhang Mat_SeqAIJ *d =(Mat_SeqAIJ*)D->data; 556d0b6147SHong Zhang Mat_MatMatMatMult *matmatmatmult=d->matmatmatmult; 566d0b6147SHong Zhang Mat BC = matmatmatmult->BC; 5771063593SHong Zhang 5871063593SHong Zhang PetscFunctionBegin; 596d0b6147SHong Zhang ierr = (BC->ops->matmultnumeric)(B,C,BC);CHKERRQ(ierr); 606d0b6147SHong Zhang ierr = (D->ops->matmultnumeric)(A,BC,D);CHKERRQ(ierr); 6171063593SHong Zhang PetscFunctionReturn(0); 6271063593SHong Zhang } 63