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;
 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);

 24:   PetscRandomCreate(PETSC_COMM_WORLD,&r);
 25:   PetscRandomSetFromOptions(r);

 27:   /* create a aij matrix A */
 28:   MatCreate(PETSC_COMM_WORLD,&A);
 29:   MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,N,M);
 30:   if (size == 1){
 31:     MatSetType(A,MATSEQAIJ);
 32:   } else {
 33:     MatSetType(A,MATMPIAIJ);
 34:   }
 35:   nza  = (PetscInt)(.3*M); /* num of nozeros in each row of A */
 36:   MatSeqAIJSetPreallocation(A,nza,PETSC_NULL);
 37:   MatMPIAIJSetPreallocation(A,nza,PETSC_NULL,nza,PETSC_NULL);
 38:   MatGetOwnershipRange(A,&rstart,&rend);
 39:   PetscMalloc((nza+1)*sizeof(PetscScalar),&a);
 40:   for (i=rstart; i<rend; i++) {
 41:     for (j=0; j<nza; j++) {
 42:       PetscRandomGetValue(r,&a[j]);
 43:       col   = (PetscInt)(M*PetscRealPart(a[j]));
 44:       MatSetValues(A,1,&i,1,&col,&a[j],ADD_VALUES);
 45:     }
 46:   }
 47:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 48:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

 50:   /* create a dense matrix B */
 51:   MatGetLocalSize(A,&am,&an);
 52:   MatCreate(PETSC_COMM_WORLD,&B);
 53:   MatSetSizes(B,an,PETSC_DECIDE,PETSC_DECIDE,N);
 54:   if (size == 1){
 55:     MatSetType(B,MATSEQDENSE);
 56:   } else {
 57:     MatSetType(B,MATMPIDENSE);
 58:   }
 59:   MatSeqDenseSetPreallocation(B,PETSC_NULL);
 60:   MatMPIDenseSetPreallocation(B,PETSC_NULL);
 61:   MatGetLocalSize(B,&bm,&bn);
 62:   MatGetArray(B,&array);
 63:   k = 0;
 64:   for (j=0; j<N; j++){ /* local column-wise entries */
 65:     for (i=0; i<bm; i++){
 66:       PetscRandomGetValue(r,&array[k++]);
 67:     }
 68:   }
 69:   MatRestoreArray(B,&array);
 70:   PetscRandomDestroy(r);
 71:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
 72:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);

 74:   /* Used for testing with Matlab as checker 
 75:   PetscViewerSetFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_MATLAB);
 76:   MatView(A,PETSC_VIEWER_STDOUT_WORLD);  
 77:   MatView(B,PETSC_VIEWER_STDOUT_WORLD);  */

 79:   /* Test MatMatMult() */
 80:   MatMatMult(A,B,MAT_INITIAL_MATRIX,fill,&C);

 82:   /* Used for testing with Matlab as checker 
 83:   MatView(C,PETSC_VIEWER_STDOUT_WORLD);  
 84:   PetscPrintf(PETSC_COMM_WORLD,"max(max(abs(Mat_3 - Mat_1*Mat_2)))\n");*/

 86:   MatMatMultSymbolic(A,B,1.0,&D);
 87:   for (i=0; i<2; i++){
 88:     MatMatMultNumeric(A,B,D);
 89:   }
 90:   MatEqual(C,D,&equal);
 91:   if (!equal) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"C != D");

 93:   MatDestroy(D);
 94:   MatDestroy(C);
 95:   MatDestroy(B);
 96:   MatDestroy(A);
 97:   PetscFree(a);
 98:   PetscFinalize();
 99:   return(0);
100: }