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

 3:  #include petscmat.h

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

 17:   PetscInitialize(&argc,&argv,(char *)0,help);
 18:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 19:   if (size != 1) SETERRQ(PETSC_ERR_SUP,"This is a uniprocessor example only!");
 20:   PetscOptionsGetInt(PETSC_NULL,"-M",&M,PETSC_NULL);
 21:   PetscOptionsGetInt(PETSC_NULL,"-N",&N,PETSC_NULL);
 22:   MatCreate(PETSC_COMM_SELF,&A);
 23:   MatSetSizes(A,M,N,M,N);
 24:   MatSetType(A,MATSEQDENSE);
 25:   MatSeqDenseSetPreallocation(A,PETSC_NULL);
 26:   PetscRandomCreate(PETSC_COMM_SELF,&r);
 27:   PetscRandomSetFromOptions(r);
 28:   MatGetArray(A,&array);
 29:   k = 0;
 30:   for (i=0; i<M; i++){
 31:     for (j=0; j<N; j++){
 32:       PetscRandomGetValue(r,&array[k++]);
 33:     }
 34:   }
 35:   MatRestoreArray(A,&array);
 36:   PetscRandomDestroy(r);
 37:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 38:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

 40:   /* Test MatMatMult() */
 41:   MatTranspose(A,&B); /* B = A^T */
 42:   MatMatMult(B,A,MAT_INITIAL_MATRIX,fill,&C); /* C = B*A = A^T*A */

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

 53:   /* Test MatMatMultTranspose() */
 54:   MatMatMultTranspose(A,A,MAT_INITIAL_MATRIX,fill,&D); /* D = A^T*A */
 55:   if (!equal) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"C != D");
 56:   MatDestroy(D);

 58:   MatDestroy(C);
 59:   MatDestroy(B);
 60:   MatDestroy(A);
 61: 
 62:   PetscFinalize();
 63:   return(0);
 64: }