Actual source code: ex122.c
1: static char help[] = "Test MatMatMult() for AIJ and Dense matrices.\n\n";
3: #include petscmat.h
7: int main(int argc,char **argv)
8: {
9: Mat A,B,C,D;
10: PetscInt i,j,k,M=10,N=5;
11: PetscScalar *array,*a;
13: PetscRandom r;
14: PetscTruth equal=PETSC_FALSE;
15: PetscReal fill = 1.0;
16: PetscInt rstart,rend,nza,col,am,an,bm,bn;
18: PetscInitialize(&argc,&argv,(char *)0,help);
19: PetscOptionsGetInt(PETSC_NULL,"-M",&M,PETSC_NULL);
20: PetscOptionsGetInt(PETSC_NULL,"-N",&N,PETSC_NULL);
21: PetscOptionsGetReal(PETSC_NULL,"-fill",&fill,PETSC_NULL);
23: PetscRandomCreate(PETSC_COMM_WORLD,&r);
24: PetscRandomSetFromOptions(r);
26: /* create a aij matrix A */
27: MatCreate(PETSC_COMM_WORLD,&A);
28: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,N,M);
29: MatSetType(A,MATAIJ);
30: nza = (PetscInt)(.3*M); /* num of nozeros in each row of A */
31: MatSeqAIJSetPreallocation(A,nza,PETSC_NULL);
32: MatMPIAIJSetPreallocation(A,nza,PETSC_NULL,nza,PETSC_NULL);
33: MatGetOwnershipRange(A,&rstart,&rend);
34: PetscMalloc((nza+1)*sizeof(PetscScalar),&a);
35: for (i=rstart; i<rend; i++) {
36: for (j=0; j<nza; j++) {
37: PetscRandomGetValue(r,&a[j]);
38: col = (PetscInt)(M*PetscRealPart(a[j]));
39: MatSetValues(A,1,&i,1,&col,&a[j],ADD_VALUES);
40: }
41: }
42: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
43: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
45: /* create a dense matrix B */
46: MatGetLocalSize(A,&am,&an);
47: MatCreate(PETSC_COMM_WORLD,&B);
48: MatSetSizes(B,PETSC_DECIDE,am,N,PETSC_DECIDE);
49: MatSetType(B,MATDENSE);
50: MatSeqDenseSetPreallocation(B,PETSC_NULL);
51: MatMPIDenseSetPreallocation(B,PETSC_NULL);
52: MatGetLocalSize(B,&bm,&bn);
53: MatGetArray(B,&array);
54: k = 0;
55: for (j=0; j<N; j++){ /* local column-wise entries */
56: for (i=0; i<bm; i++){
57: PetscRandomGetValue(r,&array[k++]);
58: }
59: }
60: MatRestoreArray(B,&array);
61: PetscRandomDestroy(r);
62: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
63: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
66: /* Test MatMatMult() */
67: MatMatMult(B,A,MAT_INITIAL_MATRIX,fill,&C);
69: /*
70: PetscViewerSetFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_MATLAB);
71: MatView(C,0);
72: MatView(B,0);
73: MatView(A,0);
74: */
75:
77: MatMatMultSymbolic(B,A,fill,&D);
78: for (i=0; i<2; i++){
79: MatMatMultNumeric(B,A,D);
80: }
81: MatEqual(C,D,&equal);
82: if (!equal) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"C != D");
84: MatDestroy(D);
85: MatDestroy(C);
86: MatDestroy(B);
87: MatDestroy(A);
88: PetscFree(a);
89: PetscFinalize();
90: return(0);
91: }