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