1: static char help[] = "Tests block sparse-dense matrix products.\n\n";

  3: #include <petscmat.h>

  5: static PetscErrorCode CheckEqual(Mat A, Mat B)
  6: {
  7:   Mat       D;
  8:   PetscReal norm;

 10:   PetscFunctionBegin;
 11:   PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &D));
 12:   PetscCall(MatAXPY(D, -1.0, B, SAME_NONZERO_PATTERN));
 13:   PetscCall(MatNorm(D, NORM_FROBENIUS, &norm));
 14:   PetscCheck(norm <= 100.0 * PETSC_MACHINE_EPSILON, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Matrix product error %g", (double)norm);
 15:   PetscCall(MatDestroy(&D));
 16:   PetscFunctionReturn(PETSC_SUCCESS);
 17: }

 19: int main(int argc, char **argv)
 20: {
 21:   Mat         A, Alibxsmm = NULL, AlibxsmmDuplicate = NULL, Aref, B, Cbaij, Clibxsmm = NULL, Cref;
 22:   MatType     type, actual;
 23:   PetscInt    rstart, rend, M, ncols = 5, bs = 2;
 24:   PetscMPIInt size;
 25:   PetscBool   match, testlibxsmm = PETSC_FALSE;
 26:   PetscScalar value;

 28:   PetscFunctionBeginUser;
 29:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
 30:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_libxsmm", &testlibxsmm, NULL));
 31:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
 32:   M = 4 * size;
 33:   PetscCall(MatCreateBAIJ(PETSC_COMM_WORLD, bs, PETSC_DECIDE, PETSC_DECIDE, M, M, 3, NULL, 3, NULL, &A));
 34:   PetscCall(MatGetOwnershipRange(A, &rstart, &rend));
 35:   for (PetscInt i = rstart; i < rend; i++) {
 36:     PetscInt brow = i / bs;

 38:     for (PetscInt bcol = PetscMax(brow - 1, 0); bcol <= PetscMin(brow + 1, M / bs - 1); bcol++) {
 39:       for (PetscInt j = bcol * bs; j < (bcol + 1) * bs; j++) {
 40:         value = 1.0 + i + 0.25 * j;
 41:         PetscCall(MatSetValue(A, i, j, value, INSERT_VALUES));
 42:       }
 43:     }
 44:   }
 45:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
 46:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));

 48:   PetscCall(MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, M, ncols, NULL, &B));
 49:   PetscCall(MatGetOwnershipRange(B, &rstart, &rend));
 50:   for (PetscInt i = rstart; i < rend; i++) {
 51:     for (PetscInt j = 0; j < ncols; j++) {
 52:       value = 0.5 + 0.125 * i - 0.25 * j;
 53:       PetscCall(MatSetValue(B, i, j, value, INSERT_VALUES));
 54:     }
 55:   }
 56:   PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
 57:   PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));

 59:   type = size == 1 ? MATSEQBAIJ : MATMPIBAIJ;
 60:   PetscCall(PetscObjectTypeCompare((PetscObject)A, type, &match));
 61:   PetscCall(MatGetType(A, &actual));
 62:   PetscCheck(match, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Unexpected BAIJ matrix type %s", actual);
 63:   PetscCall(MatConvert(A, MATAIJ, MAT_INITIAL_MATRIX, &Aref));
 64:   if (testlibxsmm) {
 65:     PetscCall(MatConvert(A, MATBAIJLIBXSMM, MAT_INITIAL_MATRIX, &Alibxsmm));
 66:     type = size == 1 ? MATSEQBAIJLIBXSMM : MATMPIBAIJLIBXSMM;
 67:     PetscCall(PetscObjectTypeCompare((PetscObject)Alibxsmm, type, &match));
 68:     PetscCall(MatGetType(Alibxsmm, &actual));
 69:     PetscCheck(match, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Unexpected converted matrix type %s", actual);
 70:     PetscCall(MatDuplicate(Alibxsmm, MAT_COPY_VALUES, &AlibxsmmDuplicate));
 71:     PetscCall(PetscObjectTypeCompare((PetscObject)AlibxsmmDuplicate, type, &match));
 72:     PetscCall(MatGetType(AlibxsmmDuplicate, &actual));
 73:     PetscCheck(match, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Unexpected duplicated matrix type %s", actual);
 74:   }

 76:   PetscCall(MatMatMult(Aref, B, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &Cref));
 77:   PetscCall(MatMatMult(A, B, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &Cbaij));
 78:   PetscCall(CheckEqual(Cbaij, Cref));
 79:   if (testlibxsmm) {
 80:     PetscCall(MatMatMult(AlibxsmmDuplicate, B, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &Clibxsmm));
 81:     PetscCall(CheckEqual(Clibxsmm, Cref));
 82:   }
 83:   PetscCall(MatScale(B, -0.5));
 84:   PetscCall(MatMatMult(Aref, B, MAT_REUSE_MATRIX, PETSC_DETERMINE, &Cref));
 85:   PetscCall(MatMatMult(A, B, MAT_REUSE_MATRIX, PETSC_DETERMINE, &Cbaij));
 86:   PetscCall(CheckEqual(Cbaij, Cref));
 87:   if (testlibxsmm) {
 88:     PetscCall(MatMatMult(AlibxsmmDuplicate, B, MAT_REUSE_MATRIX, PETSC_DETERMINE, &Clibxsmm));
 89:     PetscCall(CheckEqual(Clibxsmm, Cref));
 90:   }

 92:   if (testlibxsmm) {
 93:     type = size == 1 ? MATSEQBAIJ : MATMPIBAIJ;
 94:     PetscCall(MatConvert(Alibxsmm, type, MAT_INPLACE_MATRIX, &Alibxsmm));
 95:     PetscCall(PetscObjectTypeCompare((PetscObject)Alibxsmm, type, &match));
 96:     PetscCall(MatGetType(Alibxsmm, &actual));
 97:     PetscCheck(match, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Unexpected reverted matrix type %s", actual);
 98:   }

100:   PetscCall(MatDestroy(&Clibxsmm));
101:   PetscCall(MatDestroy(&Cbaij));
102:   PetscCall(MatDestroy(&Cref));
103:   PetscCall(MatDestroy(&B));
104:   PetscCall(MatDestroy(&Aref));
105:   PetscCall(MatDestroy(&AlibxsmmDuplicate));
106:   PetscCall(MatDestroy(&Alibxsmm));
107:   PetscCall(MatDestroy(&A));
108:   PetscCall(PetscFinalize());
109:   return 0;
110: }

112: /*TEST

114:   testset:
115:     output_file: output/empty.out

117:     test:
118:       suffix: baij_seq

120:     test:
121:       suffix: baij_mpi
122:       nsize: 2

124:     test:
125:       suffix: baij_mpi_batch
126:       nsize: 2
127:       args: -matproduct_batch_size 3

129:   testset:
130:     requires: libxsmm !complex
131:     args: -test_libxsmm
132:     output_file: output/empty.out

134:     test:
135:       suffix: libxsmm_seq

137:     test:
138:       suffix: libxsmm_mpi
139:       nsize: 2

141:     test:
142:       suffix: libxsmm_mpi_batch
143:       nsize: 2
144:       args: -matproduct_batch_size 3

146: TEST*/