Actual source code: ex123.c
1: static char help[] = "Test MatMatMult() for Dense matrices.\n\n";
3: #include petscmat.h
7: int main(int argc,char **argv)
8: {
9: Mat A,B,C,D;
10: PetscInt m,i,k,M=10,N=5;
11: PetscScalar *array;
13: PetscRandom r;
14: PetscTruth equal;
15: PetscReal fill = 1.0;
17: PetscInitialize(&argc,&argv,(char *)0,help);
18: PetscOptionsGetInt(PETSC_NULL,"-M",&M,PETSC_NULL);
19: PetscOptionsGetInt(PETSC_NULL,"-N",&N,PETSC_NULL);
20: MatCreate(PETSC_COMM_WORLD,&A);
21: MatSetSizes(A,PETSC_DETERMINE,PETSC_DETERMINE,M,N);
22: MatSetType(A,MATDENSE);
23: MatSeqDenseSetPreallocation(A,PETSC_NULL);
24: PetscRandomCreate(PETSC_COMM_SELF,&r);
25: PetscRandomSetFromOptions(r);
26: MatGetLocalSize(A,&m,PETSC_NULL);
27: MatGetArray(A,&array);
28: k = 0;
29: for (i=0; i<m*N; i++){
30: PetscRandomGetValue(r,&array[k++]);
31: }
32: MatRestoreArray(A,&array);
33: PetscRandomDestroy(r);
34: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
35: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
37: /* Test MatMatMult() */
38: MatTranspose(A,MAT_INITIAL_MATRIX,&B); /* B = A^T */
39: PetscViewerSetFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_MATLAB);
40: MatView(B,0);
41: MatView(A,0);
42: MatMatMult(B,A,MAT_INITIAL_MATRIX,fill,&C); /* C = B*A = A^T*A */
43: MatView(C,0);
45: MatMatMultSymbolic(B,A,fill,&D); /* D = B*A = A^T*A */
46: for (i=0; i<2; i++){
47: /* Repeat the numeric product to test reuse of the previous symbolic product */
48: MatMatMultNumeric(B,A,D);
49: }
50: MatEqual(C,D,&equal);
51: if (!equal) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"C != D");
52: MatDestroy(D);
55: MatDestroy(C);
56: MatDestroy(B);
57: MatDestroy(A);
58:
59: PetscFinalize();
60: return(0);
61: }