Actual source code: ex101.c
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.;
13: PetscMPIInt size;
14: PetscBool flg;
16: PetscInitialize(&argc,&argv,(char*)0,help);
17: MPI_Comm_size(PETSC_COMM_WORLD,&size);
20: /* Create MAIJ matrix, P */
21: MatCreate(PETSC_COMM_SELF,&pA);
22: MatSetSizes(pA,3,3,3,3);
23: MatSetType(pA,MATSEQAIJ);
24: MatSetUp(pA);
25: MatSetOption(pA,MAT_IGNORE_ZERO_ENTRIES,PETSC_TRUE);
26: MatSetValues(pA,3,pij,3,pij,pa,ADD_VALUES);
27: MatAssemblyBegin(pA,MAT_FINAL_ASSEMBLY);
28: MatAssemblyEnd(pA,MAT_FINAL_ASSEMBLY);
29: MatCreateMAIJ(pA,3,&P);
30: MatDestroy(&pA);
32: /* Create AIJ equivalent matrix, aijP, for comparison testing */
33: MatConvert(P,MATSEQAIJ,MAT_INITIAL_MATRIX,&aijP);
35: /* Create AIJ matrix A */
36: MatCreate(PETSC_COMM_SELF,&A);
37: MatSetSizes(A,9,9,9,9);
38: MatSetType(A,MATSEQAIJ);
39: MatSetUp(A);
40: MatSetOption(A,MAT_IGNORE_ZERO_ENTRIES,PETSC_TRUE);
41: MatSetValues(A,3,aij[0],3,aij[0],pa,ADD_VALUES);
42: MatSetValues(A,3,aij[1],3,aij[1],pa,ADD_VALUES);
43: MatSetValues(A,3,aij[2],3,aij[2],pa,ADD_VALUES);
44: for (i=0; i<9; i++) {
45: MatSetValue(A,i,i,one,ADD_VALUES);
46: }
47: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
48: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
50: /* Perform PtAP_SeqAIJ_SeqAIJ for comparison testing */
51: MatPtAP(A,aijP,MAT_INITIAL_MATRIX,1.,&C);
53: /* Perform PtAP_SeqAIJ_SeqMAIJ */
54: /* Developer API */
55: MatProductCreate(A,P,NULL,&mC);
56: MatProductSetType(mC,MATPRODUCT_PtAP);
57: MatProductSetAlgorithm(mC,"default");
58: MatProductSetFill(mC,PETSC_DEFAULT);
59: MatProductSetFromOptions(mC);
60: MatProductSymbolic(mC);
61: MatProductNumeric(mC);
62: MatProductNumeric(mC);
64: /* Check mC = C */
65: MatEqual(C,mC,&flg);
67: MatDestroy(&mC);
69: /* User API */
70: MatPtAP(A,P,MAT_INITIAL_MATRIX,1.,&mC);
71: MatPtAP(A,P,MAT_REUSE_MATRIX,1.,&mC);
73: /* Check mC = C */
74: MatEqual(C,mC,&flg);
76: MatDestroy(&mC);
78: /* Cleanup */
79: MatDestroy(&P);
80: MatDestroy(&aijP);
81: MatDestroy(&A);
82: MatDestroy(&C);
83: PetscFinalize();
84: return 0;
85: }
87: /*TEST
89: test:
90: output_file: output/ex101.out
92: TEST*/