Actual source code: ex101.c
petsc-3.13.6 2020-09-29
1: static char help[] = "Testing PtAP for SeqMAIJ matrix, P, with SeqAIJ matrix, A.\n\n";
3: #include <petscmat.h>
5: int main(int argc,char **argv)
6: {
7: Mat pA,P,aijP;
8: PetscScalar pa[]={1.,-1.,0.,0.,1.,-1.,0.,0.,1.};
9: PetscInt i,pij[]={0,1,2};
10: PetscInt aij[3][3]={{0,1,2},{3,4,5},{6,7,8}};
11: Mat A,mC,C;
12: PetscScalar one=1.;
14: PetscMPIInt size;
15: PetscBool flg;
17: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
18: MPI_Comm_size(PETSC_COMM_WORLD,&size);
19: if (size != 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"This is a uniprocessor example only!");
21: /* Create MAIJ matrix, P */
22: MatCreate(PETSC_COMM_SELF,&pA);
23: MatSetSizes(pA,3,3,3,3);
24: MatSetType(pA,MATSEQAIJ);
25: MatSetUp(pA);
26: MatSetOption(pA,MAT_IGNORE_ZERO_ENTRIES,PETSC_TRUE);
27: MatSetValues(pA,3,pij,3,pij,pa,ADD_VALUES);
28: MatAssemblyBegin(pA,MAT_FINAL_ASSEMBLY);
29: MatAssemblyEnd(pA,MAT_FINAL_ASSEMBLY);
30: MatCreateMAIJ(pA,3,&P);
31: MatDestroy(&pA);
33: /* Create AIJ equivalent matrix, aijP, for comparison testing */
34: MatConvert(P,MATSEQAIJ,MAT_INITIAL_MATRIX,&aijP);
36: /* Create AIJ matrix A */
37: MatCreate(PETSC_COMM_SELF,&A);
38: MatSetSizes(A,9,9,9,9);
39: MatSetType(A,MATSEQAIJ);
40: MatSetUp(A);
41: MatSetOption(A,MAT_IGNORE_ZERO_ENTRIES,PETSC_TRUE);
42: MatSetValues(A,3,aij[0],3,aij[0],pa,ADD_VALUES);
43: MatSetValues(A,3,aij[1],3,aij[1],pa,ADD_VALUES);
44: MatSetValues(A,3,aij[2],3,aij[2],pa,ADD_VALUES);
45: for (i=0; i<9; i++) {
46: MatSetValue(A,i,i,one,ADD_VALUES);
47: }
48: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
49: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
51: /* Perform PtAP_SeqAIJ_SeqAIJ for comparison testing */
52: MatPtAP(A,aijP,MAT_INITIAL_MATRIX,1.,&C);
54: /* Perform PtAP_SeqAIJ_SeqMAIJ */
55: /* Developer API */
56: MatProductCreate(A,P,NULL,&mC);
57: MatProductSetType(mC,MATPRODUCT_PtAP);
58: MatProductSetAlgorithm(mC,"default");
59: MatProductSetFill(mC,PETSC_DEFAULT);
60: MatProductSetFromOptions(mC);
61: MatProductSymbolic(mC);
62: MatProductNumeric(mC);
63: MatProductNumeric(mC);
65: /* Check mC = C */
66: MatEqual(C,mC,&flg);
67: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"MatProduct C != mC");
68: MatDestroy(&mC);
70: /* User API */
71: MatPtAP(A,P,MAT_INITIAL_MATRIX,1.,&mC);
72: MatPtAP(A,P,MAT_REUSE_MATRIX,1.,&mC);
74: /* Check mC = C */
75: MatEqual(C,mC,&flg);
76: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"MatPtAP C != mC");
77: MatDestroy(&mC);
79: /* Cleanup */
80: MatDestroy(&P);
81: MatDestroy(&aijP);
82: MatDestroy(&A);
83: MatDestroy(&C);
84: PetscFinalize();
85: return ierr;
86: }
88: /*TEST
90: test:
91: output_file: output/ex101.out
93: TEST*/