Actual source code: ex250.c

  1: static char help[] = "Test Mat products \n\n";

  3: #include <petscmat.h>
  4: int main(int argc,char **args)
  5: {
  6:   Mat            A = NULL,B=NULL,C=NULL,D=NULL,E=NULL;
  7:   PetscInt       k;
  8:   const PetscInt M = 18,N = 18;
  9:   PetscMPIInt    rank;

 11:   /* A, B are 18 x 18 nonsymmetric matrices and have the same sparsity pattern but different values.
 12:      Big enough to have complex communication patterns but still small enough for debugging.
 13:   */
 14:   PetscInt Ai[] = {0, 0, 1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 5,  6, 6, 7, 7, 8, 8,  9,  9, 10, 10, 11, 11, 12, 12, 13, 13, 13, 14, 14, 15, 15, 16, 16, 17, 17};
 15:   PetscInt Aj[] = {0, 1, 2, 7, 3, 8, 4, 9, 5, 8, 2, 6, 11, 0, 7, 1, 6, 2, 4, 10, 16, 11, 15, 12, 17, 12, 13, 14, 15, 17, 11, 13,  3, 16,  9, 15, 11, 13};
 16:   PetscInt Bi[] = {0, 0, 1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 5,  6, 6, 7, 7, 8, 8,  9,  9, 10, 10, 11, 11, 12, 12, 13, 13, 13, 14, 14, 15, 15, 16, 16, 17, 17};
 17:   PetscInt Bj[] = {0, 1, 2, 7, 3, 8, 4, 9, 5, 8, 2, 6, 11, 0, 7, 1, 6, 2, 4, 10, 16, 11, 15, 12, 17, 12, 13, 14, 15, 17, 11, 13,  3, 16,  9, 15, 11, 13};

 19:   PetscInt Annz = sizeof(Ai)/sizeof(PetscInt);
 20:   PetscInt Bnnz = sizeof(Bi)/sizeof(PetscInt);

 22:   PetscInitialize(&argc,&args,(char*)0,help);
 23:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);

 25:   MatCreate(PETSC_COMM_WORLD,&A);
 26:   MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,M,N);
 27:   MatSetFromOptions(A);
 28:   MatSeqAIJSetPreallocation(A,2,NULL);
 29:   MatMPIAIJSetPreallocation(A,2,NULL,2,NULL);
 30:   MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);

 32:   if (rank == 0) {
 33:     for (k=0; k<Annz; k++) MatSetValue(A,Ai[k],Aj[k],Ai[k]+Aj[k]+1.0,INSERT_VALUES);
 34:   }

 36:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 37:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

 39:   MatCreate(PETSC_COMM_WORLD,&B);
 40:   MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,M,N);
 41:   MatSetFromOptions(B);
 42:   MatSeqAIJSetPreallocation(B,2,NULL);
 43:   MatMPIAIJSetPreallocation(B,2,NULL,2,NULL);
 44:   MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);

 46:   if (rank == 0) {
 47:     for (k=0; k<Bnnz; k++) MatSetValue(B,Bi[k],Bj[k],Bi[k]+Bj[k]+2.0,INSERT_VALUES);
 48:   }
 49:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
 50:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);

 52:   MatMatMult(A,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&C);
 53:   MatView(C,PETSC_VIEWER_STDOUT_WORLD);

 55:   /* B, A have the same nonzero pattern, so it is legitimate to do so */
 56:   MatMatMult(B,A,MAT_REUSE_MATRIX,PETSC_DEFAULT,&C);
 57:   MatView(C,PETSC_VIEWER_STDOUT_WORLD);

 59:   MatTransposeMatMult(A,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&D);
 60:   MatView(D, PETSC_VIEWER_STDOUT_WORLD);

 62:   MatPtAP(A,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&E);
 63:   MatView(E,PETSC_VIEWER_STDOUT_WORLD);

 65:   MatDestroy(&A);
 66:   MatDestroy(&B);
 67:   MatDestroy(&C);
 68:   MatDestroy(&D);
 69:   MatDestroy(&E);

 71:   PetscFinalize();
 72:   return 0;
 73: }

 75: /*TEST
 76:   testset:
 77:     filter: grep -ve type -ve "Mat Object"
 78:     output_file: output/ex250_1.out

 80:     test:
 81:       suffix: 1
 82:       nsize: {{1 3}}
 83:       args: -mat_type aij

 85:     test:
 86:       suffix: 2
 87:       nsize: {{3 4}}
 88:       args: -mat_type aij -matmatmult_via backend -matptap_via backend -mattransposematmult_via backend

 90:     test:
 91:       suffix: cuda
 92:       requires: cuda
 93:       nsize: {{1 3 4}}
 94:       args: -mat_type aijcusparse

 96:     test:
 97:       suffix: kok
 98:       requires: !sycl kokkos_kernels
 99:       nsize: {{1 3 4}}
100:       args: -mat_type aijkokkos

102: TEST*/