Actual source code: ex101.c
petsc-3.3-p6 2013-02-11
1: static char help[] = "Testing PtAP for SeqMAIJ matrix, P, with SeqAIJ matrix, A.\n\n";
3: #include <petscmat.h>
7: int main(int argc,char **argv) {
8: Mat pA,P,aijP;
9: PetscScalar pa[]={1.,-1.,0.,0.,1.,-1.,0.,0.,1.};
10: PetscInt pij[]={0,1,2};
11: PetscInt aij[3][3]={{0,1,2},{3,4,5},{6,7,8}};
12: Mat A,mC,C;
13: PetscScalar one=1.;
16: PetscInitialize(&argc,&argv,(char *)0,help);
18: /* Create MAIJ matrix, P */
19: MatCreate(PETSC_COMM_SELF,&pA);
20: MatSetSizes(pA,3,3,3,3);
21: MatSetType(pA,MATSEQAIJ);
22: MatSetUp(pA);
23: MatSetOption(pA,MAT_IGNORE_ZERO_ENTRIES,PETSC_TRUE);
24: MatSetValues(pA,3,pij,3,pij,pa,ADD_VALUES);
25: MatAssemblyBegin(pA,MAT_FINAL_ASSEMBLY);
26: MatAssemblyEnd(pA,MAT_FINAL_ASSEMBLY);
27: MatCreateMAIJ(pA,3,&P);
28: MatDestroy(&pA);
30: /* Create AIJ equivalent matrix, aijP, for comparison testing */
31: MatConvert(P,MATSEQAIJ,MAT_INITIAL_MATRIX,&aijP);
33: /* Create AIJ matrix, A */
34: MatCreate(PETSC_COMM_SELF,&A);
35: MatSetSizes(A,9,9,9,9);
36: MatSetType(A,MATSEQAIJ);
37: MatSetUp(A);
38: MatSetOption(A,MAT_IGNORE_ZERO_ENTRIES,PETSC_TRUE);
39: MatSetValues(A,3,aij[0],3,aij[0],pa,ADD_VALUES);
40: MatSetValues(A,3,aij[1],3,aij[1],pa,ADD_VALUES);
41: MatSetValues(A,3,aij[2],3,aij[2],pa,ADD_VALUES);
42: {int i;
43: for (i=0;i<9;i++) {
44: MatSetValue(A,i,i,one,ADD_VALUES);
45: }
46: }
47: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
48: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
49:
50: /* Perform SeqAIJ_SeqMAIJ PtAP */
51: MatPtAP(A,P,MAT_INITIAL_MATRIX,1.,&mC);
52: MatView(mC,PETSC_VIEWER_STDOUT_SELF);
54: /* Perform SeqAIJ_SeqAIJ PtAP for comparison testing */
55: MatPtAP(A,aijP,MAT_INITIAL_MATRIX,1.,&C);
56: MatView(C,PETSC_VIEWER_STDOUT_SELF);
58: /* Perform diff of two matrices */
59: MatAXPY(C,-1.0,mC,DIFFERENT_NONZERO_PATTERN);
60: /* Note: We should be able to use SAME_NONZERO_PATTERN on the line above, */
61: /* but don't because this flag doesn't assist testing. */
62: MatView(C,PETSC_VIEWER_STDOUT_SELF);
64: /* Cleanup */
65: MatDestroy(&P);
66: MatDestroy(&aijP);
67: MatDestroy(&A);
68: MatDestroy(&C);
69: MatDestroy(&mC);
70: PetscFinalize();
71: return(0);
72: }