Actual source code: ex123.c

  1: static char help[] = "Test MatMatMult() for Dense matrices.\n\n";

 3:  #include petscmat.h

  7: int main(int argc,char **argv)
  8: {
  9:   Mat            A,B,C,D;
 10:   PetscInt       m,i,k,M=10,N=5;
 11:   PetscScalar    *array;
 13:   PetscRandom    r;
 14:   PetscTruth     equal;
 15:   PetscReal      fill = 1.0;

 17:   PetscInitialize(&argc,&argv,(char *)0,help);
 18:   PetscOptionsGetInt(PETSC_NULL,"-M",&M,PETSC_NULL);
 19:   PetscOptionsGetInt(PETSC_NULL,"-N",&N,PETSC_NULL);
 20:   MatCreate(PETSC_COMM_WORLD,&A);
 21:   MatSetSizes(A,PETSC_DETERMINE,PETSC_DETERMINE,M,N);
 22:   MatSetType(A,MATDENSE);
 23:   MatSeqDenseSetPreallocation(A,PETSC_NULL);
 24:   PetscRandomCreate(PETSC_COMM_SELF,&r);
 25:   PetscRandomSetFromOptions(r);
 26:   MatGetLocalSize(A,&m,PETSC_NULL);
 27:   MatGetArray(A,&array);
 28:   k = 0;
 29:   for (i=0; i<m*N; i++){
 30:     PetscRandomGetValue(r,&array[k++]);
 31:   }
 32:   MatRestoreArray(A,&array);
 33:   PetscRandomDestroy(r);
 34:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 35:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

 37:   /* Test MatMatMult() */
 38:   MatTranspose(A,MAT_INITIAL_MATRIX,&B); /* B = A^T */
 39:   PetscViewerSetFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_MATLAB);
 40:   MatView(B,0);
 41:   MatView(A,0);
 42:   MatMatMult(B,A,MAT_INITIAL_MATRIX,fill,&C); /* C = B*A = A^T*A */
 43:   MatView(C,0);

 45:   MatMatMultSymbolic(B,A,fill,&D); /* D = B*A = A^T*A */
 46:   for (i=0; i<2; i++){
 47:     /* Repeat the numeric product to test reuse of the previous symbolic product */
 48:     MatMatMultNumeric(B,A,D);
 49:   }
 50:   MatEqual(C,D,&equal);
 51:   if (!equal) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"C != D");
 52:   MatDestroy(D);


 55:   MatDestroy(C);
 56:   MatDestroy(B);
 57:   MatDestroy(A);
 58: 
 59:   PetscFinalize();
 60:   return(0);
 61: }