Actual source code: ex109.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,an,PETSC_DECIDE,PETSC_DECIDE,N);
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);
65: /* Test MatMatMult() */
66: MatMatMult(A,B,MAT_INITIAL_MATRIX,fill,&C);
68: MatMatMultSymbolic(A,B,fill,&D);
69: for (i=0; i<2; i++){
70: MatMatMultNumeric(A,B,D);
71: }
72: MatEqual(C,D,&equal);
73: if (!equal) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"C != D");
75: MatDestroy(D);
76: MatDestroy(C);
77: MatDestroy(B);
78: MatDestroy(A);
79: PetscFree(a);
80: PetscFinalize();
81: return(0);
82: }