xref: /petsc/src/mat/tests/ex101.c (revision 28b400f66ebc7ae0049166a2294dfcd3df27e64b)
1c4762a1bSJed Brown static char help[] = "Testing PtAP for SeqMAIJ matrix, P, with SeqAIJ matrix, A.\n\n";
2c4762a1bSJed Brown 
3c4762a1bSJed Brown #include <petscmat.h>
4c4762a1bSJed Brown 
5c4762a1bSJed Brown int main(int argc,char **argv)
6c4762a1bSJed Brown {
7c4762a1bSJed Brown   Mat            pA,P,aijP;
8c4762a1bSJed Brown   PetscScalar    pa[]={1.,-1.,0.,0.,1.,-1.,0.,0.,1.};
9c4762a1bSJed Brown   PetscInt       i,pij[]={0,1,2};
10c4762a1bSJed Brown   PetscInt       aij[3][3]={{0,1,2},{3,4,5},{6,7,8}};
11c4762a1bSJed Brown   Mat            A,mC,C;
12c4762a1bSJed Brown   PetscScalar    one=1.;
13c4762a1bSJed Brown   PetscErrorCode ierr;
14c4762a1bSJed Brown   PetscMPIInt    size;
15c20d7725SJed Brown   PetscBool      flg;
16c4762a1bSJed Brown 
17c4762a1bSJed Brown   ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
185f80ce2aSJacob Faibussowitsch   CHKERRMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size));
192c71b3e2SJacob Faibussowitsch   PetscCheckFalse(size != 1,PETSC_COMM_WORLD,PETSC_ERR_SUP,"This is a uniprocessor example only!");
20c4762a1bSJed Brown 
21c4762a1bSJed Brown   /* Create MAIJ matrix, P */
225f80ce2aSJacob Faibussowitsch   CHKERRQ(MatCreate(PETSC_COMM_SELF,&pA));
235f80ce2aSJacob Faibussowitsch   CHKERRQ(MatSetSizes(pA,3,3,3,3));
245f80ce2aSJacob Faibussowitsch   CHKERRQ(MatSetType(pA,MATSEQAIJ));
255f80ce2aSJacob Faibussowitsch   CHKERRQ(MatSetUp(pA));
265f80ce2aSJacob Faibussowitsch   CHKERRQ(MatSetOption(pA,MAT_IGNORE_ZERO_ENTRIES,PETSC_TRUE));
275f80ce2aSJacob Faibussowitsch   CHKERRQ(MatSetValues(pA,3,pij,3,pij,pa,ADD_VALUES));
285f80ce2aSJacob Faibussowitsch   CHKERRQ(MatAssemblyBegin(pA,MAT_FINAL_ASSEMBLY));
295f80ce2aSJacob Faibussowitsch   CHKERRQ(MatAssemblyEnd(pA,MAT_FINAL_ASSEMBLY));
305f80ce2aSJacob Faibussowitsch   CHKERRQ(MatCreateMAIJ(pA,3,&P));
315f80ce2aSJacob Faibussowitsch   CHKERRQ(MatDestroy(&pA));
32c4762a1bSJed Brown 
33c4762a1bSJed Brown   /* Create AIJ equivalent matrix, aijP, for comparison testing */
345f80ce2aSJacob Faibussowitsch   CHKERRQ(MatConvert(P,MATSEQAIJ,MAT_INITIAL_MATRIX,&aijP));
35c4762a1bSJed Brown 
36c20d7725SJed Brown   /* Create AIJ matrix A */
375f80ce2aSJacob Faibussowitsch   CHKERRQ(MatCreate(PETSC_COMM_SELF,&A));
385f80ce2aSJacob Faibussowitsch   CHKERRQ(MatSetSizes(A,9,9,9,9));
395f80ce2aSJacob Faibussowitsch   CHKERRQ(MatSetType(A,MATSEQAIJ));
405f80ce2aSJacob Faibussowitsch   CHKERRQ(MatSetUp(A));
415f80ce2aSJacob Faibussowitsch   CHKERRQ(MatSetOption(A,MAT_IGNORE_ZERO_ENTRIES,PETSC_TRUE));
425f80ce2aSJacob Faibussowitsch   CHKERRQ(MatSetValues(A,3,aij[0],3,aij[0],pa,ADD_VALUES));
435f80ce2aSJacob Faibussowitsch   CHKERRQ(MatSetValues(A,3,aij[1],3,aij[1],pa,ADD_VALUES));
445f80ce2aSJacob Faibussowitsch   CHKERRQ(MatSetValues(A,3,aij[2],3,aij[2],pa,ADD_VALUES));
45c4762a1bSJed Brown   for (i=0; i<9; i++) {
465f80ce2aSJacob Faibussowitsch     CHKERRQ(MatSetValue(A,i,i,one,ADD_VALUES));
47c4762a1bSJed Brown   }
485f80ce2aSJacob Faibussowitsch   CHKERRQ(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY));
495f80ce2aSJacob Faibussowitsch   CHKERRQ(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY));
50c4762a1bSJed Brown 
51c4762a1bSJed Brown   /* Perform PtAP_SeqAIJ_SeqAIJ for comparison testing */
525f80ce2aSJacob Faibussowitsch   CHKERRQ(MatPtAP(A,aijP,MAT_INITIAL_MATRIX,1.,&C));
53c20d7725SJed Brown 
54c20d7725SJed Brown   /* Perform PtAP_SeqAIJ_SeqMAIJ */
55c20d7725SJed Brown   /* Developer API */
565f80ce2aSJacob Faibussowitsch   CHKERRQ(MatProductCreate(A,P,NULL,&mC));
575f80ce2aSJacob Faibussowitsch   CHKERRQ(MatProductSetType(mC,MATPRODUCT_PtAP));
585f80ce2aSJacob Faibussowitsch   CHKERRQ(MatProductSetAlgorithm(mC,"default"));
595f80ce2aSJacob Faibussowitsch   CHKERRQ(MatProductSetFill(mC,PETSC_DEFAULT));
605f80ce2aSJacob Faibussowitsch   CHKERRQ(MatProductSetFromOptions(mC));
615f80ce2aSJacob Faibussowitsch   CHKERRQ(MatProductSymbolic(mC));
625f80ce2aSJacob Faibussowitsch   CHKERRQ(MatProductNumeric(mC));
635f80ce2aSJacob Faibussowitsch   CHKERRQ(MatProductNumeric(mC));
64c4762a1bSJed Brown 
65c4762a1bSJed Brown   /* Check mC = C */
665f80ce2aSJacob Faibussowitsch   CHKERRQ(MatEqual(C,mC,&flg));
67*28b400f6SJacob Faibussowitsch   PetscCheck(flg,PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"MatProduct C != mC");
685f80ce2aSJacob Faibussowitsch   CHKERRQ(MatDestroy(&mC));
69c20d7725SJed Brown 
70c20d7725SJed Brown   /* User API */
715f80ce2aSJacob Faibussowitsch   CHKERRQ(MatPtAP(A,P,MAT_INITIAL_MATRIX,1.,&mC));
725f80ce2aSJacob Faibussowitsch   CHKERRQ(MatPtAP(A,P,MAT_REUSE_MATRIX,1.,&mC));
73c20d7725SJed Brown 
74c20d7725SJed Brown   /* Check mC = C */
755f80ce2aSJacob Faibussowitsch   CHKERRQ(MatEqual(C,mC,&flg));
76*28b400f6SJacob Faibussowitsch   PetscCheck(flg,PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"MatPtAP C != mC");
775f80ce2aSJacob Faibussowitsch   CHKERRQ(MatDestroy(&mC));
78c4762a1bSJed Brown 
79c4762a1bSJed Brown   /* Cleanup */
805f80ce2aSJacob Faibussowitsch   CHKERRQ(MatDestroy(&P));
815f80ce2aSJacob Faibussowitsch   CHKERRQ(MatDestroy(&aijP));
825f80ce2aSJacob Faibussowitsch   CHKERRQ(MatDestroy(&A));
835f80ce2aSJacob Faibussowitsch   CHKERRQ(MatDestroy(&C));
84c4762a1bSJed Brown   ierr = PetscFinalize();
85c4762a1bSJed Brown   return ierr;
86c4762a1bSJed Brown }
87c4762a1bSJed Brown 
88c4762a1bSJed Brown /*TEST
89c4762a1bSJed Brown 
90c4762a1bSJed Brown    test:
91c4762a1bSJed Brown       output_file: output/ex101.out
92c4762a1bSJed Brown 
93c4762a1bSJed Brown TEST*/
94