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: }