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: PetscMPIInt size;
17: PetscInt rstart,rend,nza,col,am,an,bm,bn;
19: PetscInitialize(&argc,&argv,(char *)0,help);
20: MPI_Comm_size(PETSC_COMM_WORLD,&size);
21: PetscOptionsGetInt(PETSC_NULL,"-M",&M,PETSC_NULL);
22: PetscOptionsGetInt(PETSC_NULL,"-N",&N,PETSC_NULL);
23: PetscOptionsGetReal(PETSC_NULL,"-fill",&fill,PETSC_NULL);
25: PetscRandomCreate(PETSC_COMM_WORLD,&r);
26: PetscRandomSetFromOptions(r);
28: /* create a aij matrix A */
29: MatCreate(PETSC_COMM_WORLD,&A);
30: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,N,M);
31: if (size == 1){
32: MatSetType(A,MATSEQAIJ);
33: } else {
34: MatSetType(A,MATMPIAIJ);
35: }
36: nza = (PetscInt)(.3*M); /* num of nozeros in each row of A */
37: MatSeqAIJSetPreallocation(A,nza,PETSC_NULL);
38: MatMPIAIJSetPreallocation(A,nza,PETSC_NULL,nza,PETSC_NULL);
39: MatGetOwnershipRange(A,&rstart,&rend);
40: PetscMalloc((nza+1)*sizeof(PetscScalar),&a);
41: for (i=rstart; i<rend; i++) {
42: for (j=0; j<nza; j++) {
43: PetscRandomGetValue(r,&a[j]);
44: col = (PetscInt)(M*PetscRealPart(a[j]));
45: MatSetValues(A,1,&i,1,&col,&a[j],ADD_VALUES);
46: }
47: }
48: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
49: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
51: /* create a dense matrix B */
52: MatGetLocalSize(A,&am,&an);
53: MatCreate(PETSC_COMM_WORLD,&B);
54: MatSetSizes(B,an,PETSC_DECIDE,PETSC_DECIDE,N);
55: if (size == 1){
56: MatSetType(B,MATSEQDENSE);
57: } else {
58: MatSetType(B,MATMPIDENSE);
59: }
60: MatSeqDenseSetPreallocation(B,PETSC_NULL);
61: MatMPIDenseSetPreallocation(B,PETSC_NULL);
62: MatGetLocalSize(B,&bm,&bn);
63: MatGetArray(B,&array);
64: k = 0;
65: for (j=0; j<N; j++){ /* local column-wise entries */
66: for (i=0; i<bm; i++){
67: PetscRandomGetValue(r,&array[k++]);
68: }
69: }
70: MatRestoreArray(B,&array);
71: PetscRandomDestroy(r);
72: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
73: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
75: /* Test MatMatMult() */
76: MatMatMult(A,B,MAT_INITIAL_MATRIX,fill,&C);
78: MatMatMultSymbolic(A,B,fill,&D);
79: for (i=0; i<2; i++){
80: MatMatMultNumeric(A,B,D);
81: }
82: MatEqual(C,D,&equal);
83: if (!equal) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"C != D");
85: MatDestroy(D);
86: MatDestroy(C);
87: MatDestroy(B);
88: MatDestroy(A);
89: PetscFree(a);
90: PetscFinalize();
91: return(0);
92: }