1 static char help[] = "Mini-app to benchmark matrix--matrix multiplication\n\n"; 2 3 /* 4 See the paper below for more information 5 6 "KSPHPDDM and PCHPDDM: Extending PETSc with Robust Overlapping Schwarz Preconditioners and Advanced Krylov Methods", 7 P. Jolivet, J. E. Roman, and S. Zampini (2020). 8 */ 9 10 #include <petsc.h> 11 12 #if defined(PETSC_HAVE_MKL_SPARSE_OPTIMIZE) 13 #include <mkl.h> 14 #define PetscCallMKLSparse(func, args) \ 15 do { \ 16 sparse_status_t __ierr; \ 17 PetscStackPushExternal(#func); \ 18 __ierr = func args; \ 19 PetscStackPop; \ 20 PetscCheck(__ierr == SPARSE_STATUS_SUCCESS, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in %s(): error code %d", #func, (int)__ierr); \ 21 } while (0) 22 #else 23 #define PetscCallMKLSparse(func, args) \ 24 do { \ 25 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "No MKL support"); \ 26 } while (0) 27 #endif 28 29 int main(int argc, char **argv) 30 { 31 Mat A, C, D, E; 32 PetscInt nbs = 10, ntype = 10, nN = 8, m, M, trial = 5; 33 PetscViewer viewer; 34 PetscInt bs[10], N[8]; 35 char *type[10]; 36 PetscMPIInt size; 37 PetscBool flg, cuda, maij = PETSC_FALSE, check = PETSC_FALSE, trans = PETSC_FALSE, convert = PETSC_FALSE, mkl; 38 char file[PETSC_MAX_PATH_LEN]; 39 40 PetscFunctionBeginUser; 41 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 42 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 43 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only"); 44 PetscCall(PetscOptionsGetString(NULL, NULL, "-f", file, PETSC_MAX_PATH_LEN, &flg)); 45 PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_USER_INPUT, "Must indicate binary file with the -f option"); 46 PetscCall(PetscOptionsGetInt(NULL, NULL, "-trial", &trial, NULL)); 47 PetscCall(PetscOptionsGetIntArray(NULL, NULL, "-bs", bs, &nbs, &flg)); 48 if (!flg) { 49 nbs = 1; 50 bs[0] = 1; 51 } 52 PetscCall(PetscOptionsGetIntArray(NULL, NULL, "-N", N, &nN, &flg)); 53 if (!flg) { 54 nN = 8; 55 N[0] = 1; 56 N[1] = 2; 57 N[2] = 4; 58 N[3] = 8; 59 N[4] = 16; 60 N[5] = 32; 61 N[6] = 64; 62 N[7] = 128; 63 } 64 PetscCall(PetscOptionsGetStringArray(NULL, NULL, "-type", type, &ntype, &flg)); 65 if (!flg) { 66 ntype = 1; 67 PetscCall(PetscStrallocpy(MATSEQAIJ, &type[0])); 68 } 69 PetscCall(PetscOptionsGetBool(NULL, NULL, "-check", &check, NULL)); 70 PetscCall(PetscOptionsGetBool(NULL, NULL, "-trans", &trans, NULL)); 71 PetscCall(PetscOptionsGetBool(NULL, NULL, "-convert_aij", &convert, NULL)); 72 for (PetscInt j = 0; j < nbs; ++j) { 73 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, file, FILE_MODE_READ, &viewer)); 74 PetscCall(MatCreate(PETSC_COMM_WORLD, &A)); 75 PetscCall(MatSetFromOptions(A)); 76 PetscCall(MatLoad(A, viewer)); 77 PetscCall(PetscViewerDestroy(&viewer)); 78 PetscCall(PetscObjectTypeCompareAny((PetscObject)A, &flg, MATSEQAIJ, MATMPIAIJ, "")); 79 PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_USER_INPUT, "Must indicate a MatAIJ input matrix"); 80 PetscCall(MatGetSize(A, &m, &M)); 81 if (m == M) { 82 Mat oA; 83 PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &oA)); 84 PetscCall(MatAXPY(A, 1.0, oA, DIFFERENT_NONZERO_PATTERN)); 85 PetscCall(MatDestroy(&oA)); 86 } 87 PetscCall(MatGetLocalSize(A, &m, NULL)); 88 PetscCall(MatGetSize(A, &M, NULL)); 89 if (bs[j] > 1) { 90 Mat T, Tt, B; 91 const PetscScalar *ptr; 92 PetscScalar *val, *Aa; 93 const PetscInt *Ai, *Aj; 94 PetscInt An, i, k; 95 PetscBool done; 96 97 PetscCall(MatCreateDense(PETSC_COMM_SELF, bs[j], bs[j], bs[j], bs[j], NULL, &T)); 98 PetscCall(MatSetRandom(T, NULL)); 99 PetscCall(MatTranspose(T, MAT_INITIAL_MATRIX, &Tt)); 100 PetscCall(MatAXPY(T, 1.0, Tt, SAME_NONZERO_PATTERN)); 101 PetscCall(MatDestroy(&Tt)); 102 PetscCall(MatDenseGetArrayRead(T, &ptr)); 103 PetscCall(MatGetRowIJ(A, 0, PETSC_FALSE, PETSC_FALSE, &An, &Ai, &Aj, &done)); 104 PetscCheck(done && An == m, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Inconsistent sizes"); 105 PetscCall(MatSeqAIJGetArray(A, &Aa)); 106 PetscCall(MatCreate(PETSC_COMM_WORLD, &B)); 107 PetscCall(MatSetType(B, MATSEQBAIJ)); 108 PetscCall(MatSetSizes(B, bs[j] * An, bs[j] * An, PETSC_DECIDE, PETSC_DECIDE)); 109 PetscCall(PetscMalloc1(Ai[An] * bs[j] * bs[j], &val)); 110 for (i = 0; i < Ai[An]; ++i) 111 for (k = 0; k < bs[j] * bs[j]; ++k) val[i * bs[j] * bs[j] + k] = Aa[i] * ptr[k]; 112 PetscCall(MatSetOption(B, MAT_ROW_ORIENTED, PETSC_FALSE)); 113 PetscCall(MatSeqBAIJSetPreallocationCSR(B, bs[j], Ai, Aj, val)); 114 PetscCall(PetscFree(val)); 115 PetscCall(MatSeqAIJRestoreArray(A, &Aa)); 116 PetscCall(MatRestoreRowIJ(A, 0, PETSC_FALSE, PETSC_FALSE, &An, &Ai, &Aj, &done)); 117 PetscCall(MatDenseRestoreArrayRead(T, &ptr)); 118 PetscCall(MatDestroy(&T)); 119 PetscCall(MatDestroy(&A)); 120 A = B; 121 } 122 /* reconvert back to SeqAIJ before converting to the desired type later */ 123 if (!convert) E = A; 124 PetscCall(MatConvert(A, MATSEQAIJ, convert ? MAT_INITIAL_MATRIX : MAT_INPLACE_MATRIX, &E)); 125 PetscCall(MatSetOption(E, MAT_SYMMETRIC, PETSC_TRUE)); 126 for (PetscInt i = 0; i < ntype; ++i) { 127 char *tmp = NULL; 128 PetscInt *ia_ptr, *ja_ptr, k; 129 PetscScalar *a_ptr; 130 #if defined(PETSC_HAVE_MKL_SPARSE_OPTIMIZE) 131 struct matrix_descr descr; 132 sparse_matrix_t spr; 133 descr.type = SPARSE_MATRIX_TYPE_GENERAL; 134 descr.diag = SPARSE_DIAG_NON_UNIT; 135 #endif 136 if (convert) PetscCall(MatDestroy(&A)); 137 PetscCall(PetscStrstr(type[i], "mkl", &tmp)); 138 if (tmp) { 139 size_t mlen, tlen; 140 char base[256]; 141 142 mkl = PETSC_TRUE; 143 PetscCall(PetscStrlen(tmp, &mlen)); 144 PetscCall(PetscStrlen(type[i], &tlen)); 145 PetscCall(PetscStrncpy(base, type[i], tlen - mlen + 1)); 146 PetscCall(MatConvert(E, base, convert ? MAT_INITIAL_MATRIX : MAT_INPLACE_MATRIX, &A)); 147 } else { 148 mkl = PETSC_FALSE; 149 PetscCall(PetscStrstr(type[i], "maij", &tmp)); 150 if (!tmp) { 151 PetscCall(MatConvert(E, type[i], convert ? MAT_INITIAL_MATRIX : MAT_INPLACE_MATRIX, &A)); 152 } else { 153 PetscCall(MatConvert(E, MATAIJ, convert ? MAT_INITIAL_MATRIX : MAT_INPLACE_MATRIX, &A)); 154 maij = PETSC_TRUE; 155 } 156 } 157 PetscCall(PetscObjectTypeCompareAny((PetscObject)A, &cuda, MATSEQAIJCUSPARSE, MATMPIAIJCUSPARSE, "")); 158 if (mkl) { 159 const PetscInt *Ai, *Aj; 160 PetscInt An; 161 PetscBool done; 162 163 PetscCall(PetscObjectTypeCompareAny((PetscObject)A, &flg, MATSEQAIJ, MATSEQBAIJ, MATSEQSBAIJ, "")); 164 PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_USER_INPUT, "Not implemented"); 165 PetscCall(PetscObjectTypeCompare((PetscObject)A, MATSEQAIJ, &flg)); 166 PetscCall(MatGetRowIJ(A, 0, PETSC_FALSE, flg ? PETSC_FALSE : PETSC_TRUE, &An, &Ai, &Aj, &done)); 167 PetscCheck(done, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Inconsistent sizes"); 168 PetscCall(PetscMalloc1(An + 1, &ia_ptr)); 169 PetscCall(PetscMalloc1(Ai[An], &ja_ptr)); 170 if (flg) { /* SeqAIJ */ 171 for (k = 0; k < An + 1; ++k) ia_ptr[k] = Ai[k]; 172 for (k = 0; k < Ai[An]; ++k) ja_ptr[k] = Aj[k]; 173 PetscCall(MatSeqAIJGetArray(A, &a_ptr)); 174 PetscCallMKLSparse(mkl_sparse_d_create_csr, (&spr, SPARSE_INDEX_BASE_ZERO, An, An, ia_ptr, ia_ptr + 1, ja_ptr, a_ptr)); 175 } else { 176 PetscCall(PetscObjectTypeCompare((PetscObject)A, MATSEQBAIJ, &flg)); 177 if (flg) { 178 for (k = 0; k < An + 1; ++k) ia_ptr[k] = Ai[k] + 1; /* Fortran indexing to maximize cases covered by _mm routines */ 179 for (k = 0; k < Ai[An]; ++k) ja_ptr[k] = Aj[k] + 1; /* Fortran indexing to maximize cases covered by _mm routines */ 180 PetscCall(MatSeqBAIJGetArray(A, &a_ptr)); 181 PetscCallMKLSparse(mkl_sparse_d_create_bsr, (&spr, SPARSE_INDEX_BASE_ONE, SPARSE_LAYOUT_COLUMN_MAJOR, An, An, bs[j], ia_ptr, ia_ptr + 1, ja_ptr, a_ptr)); 182 } else { 183 PetscCall(PetscObjectTypeCompare((PetscObject)A, MATSEQSBAIJ, &flg)); 184 if (flg) { 185 for (k = 0; k < An + 1; ++k) ia_ptr[k] = Ai[k] + 1; /* Fortran indexing to maximize cases covered by _mm routines */ 186 for (k = 0; k < Ai[An]; ++k) ja_ptr[k] = Aj[k] + 1; /* Fortran indexing to maximize cases covered by _mm routines */ 187 PetscCall(MatSeqSBAIJGetArray(A, &a_ptr)); 188 PetscCallMKLSparse(mkl_sparse_d_create_bsr, (&spr, SPARSE_INDEX_BASE_ONE, SPARSE_LAYOUT_COLUMN_MAJOR, An, An, bs[j], ia_ptr, ia_ptr + 1, ja_ptr, a_ptr)); 189 #if defined(PETSC_HAVE_MKL_SPARSE_OPTIMIZE) 190 descr.type = SPARSE_MATRIX_TYPE_SYMMETRIC; 191 descr.mode = SPARSE_FILL_MODE_UPPER; 192 descr.diag = SPARSE_DIAG_NON_UNIT; 193 #endif 194 } 195 } 196 } 197 PetscCall(PetscObjectTypeCompare((PetscObject)A, MATSEQAIJ, &flg)); 198 PetscCall(MatRestoreRowIJ(A, 0, PETSC_FALSE, flg ? PETSC_FALSE : PETSC_TRUE, &An, &Ai, &Aj, &done)); 199 } 200 201 PetscCall(MatViewFromOptions(A, NULL, "-A_view")); 202 203 for (k = 0; k < nN; ++k) { 204 MatType Atype, Ctype; 205 PetscInt AM, AN, CM, CN, t; 206 PetscLogStage stage, tstage; 207 char stage_s[256]; 208 209 PetscCall(MatCreateDense(PETSC_COMM_WORLD, bs[j] * m, PETSC_DECIDE, bs[j] * M, N[k], NULL, &C)); 210 PetscCall(MatCreateDense(PETSC_COMM_WORLD, bs[j] * m, PETSC_DECIDE, bs[j] * M, N[k], NULL, &D)); 211 PetscCall(MatSetRandom(C, NULL)); 212 if (cuda) { /* convert to GPU if needed */ 213 PetscCall(MatConvert(C, MATDENSECUDA, MAT_INPLACE_MATRIX, &C)); 214 PetscCall(MatConvert(D, MATDENSECUDA, MAT_INPLACE_MATRIX, &D)); 215 } 216 if (mkl) { 217 if (N[k] > 1) PetscCallMKLSparse(mkl_sparse_set_mm_hint, (spr, SPARSE_OPERATION_NON_TRANSPOSE, descr, SPARSE_LAYOUT_COLUMN_MAJOR, N[k], 1 + trial)); 218 else PetscCallMKLSparse(mkl_sparse_set_mv_hint, (spr, SPARSE_OPERATION_NON_TRANSPOSE, descr, 1 + trial)); 219 PetscCallMKLSparse(mkl_sparse_set_memory_hint, (spr, SPARSE_MEMORY_AGGRESSIVE)); 220 PetscCallMKLSparse(mkl_sparse_optimize, (spr)); 221 } 222 PetscCall(MatGetType(A, &Atype)); 223 PetscCall(MatGetType(C, &Ctype)); 224 PetscCall(MatGetSize(A, &AM, &AN)); 225 PetscCall(MatGetSize(C, &CM, &CN)); 226 227 if (!maij || N[k] > 1) { 228 PetscCall(PetscSNPrintf(stage_s, sizeof(stage_s), "type_%s-bs_%" PetscInt_FMT "-N_%02d", type[i], bs[j], (int)N[k])); 229 PetscCall(PetscLogStageRegister(stage_s, &stage)); 230 } 231 if (trans && N[k] > 1) { 232 PetscCall(PetscSNPrintf(stage_s, sizeof(stage_s), "trans_type_%s-bs_%" PetscInt_FMT "-N_%02d", type[i], bs[j], (int)N[k])); 233 PetscCall(PetscLogStageRegister(stage_s, &tstage)); 234 } 235 /* A*B */ 236 if (N[k] > 1) { 237 if (!maij) { 238 PetscCall(MatProductCreateWithMat(A, C, NULL, D)); 239 PetscCall(MatProductSetType(D, MATPRODUCT_AB)); 240 PetscCall(MatProductSetFromOptions(D)); 241 PetscCall(MatProductSymbolic(D)); 242 } 243 244 if (!mkl) { 245 if (!maij) { 246 PetscCall(MatProductNumeric(D)); 247 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Benchmarking MatProduct %s: with A %s %" PetscInt_FMT "x%" PetscInt_FMT " and B %s %" PetscInt_FMT "x%" PetscInt_FMT "\n", MatProductTypes[MATPRODUCT_AB], Atype, AM, AN, Ctype, CM, CN)); 248 PetscCall(PetscLogStagePush(stage)); 249 for (t = 0; t < trial; ++t) PetscCall(MatProductNumeric(D)); 250 PetscCall(PetscLogStagePop()); 251 } else { 252 Mat E, Ct, Dt; 253 Vec cC, cD; 254 const PetscScalar *c_ptr; 255 PetscScalar *d_ptr; 256 PetscCall(MatCreateMAIJ(A, N[k], &E)); 257 PetscCall(MatDenseGetLocalMatrix(C, &Ct)); 258 PetscCall(MatDenseGetLocalMatrix(D, &Dt)); 259 PetscCall(MatTranspose(Ct, MAT_INPLACE_MATRIX, &Ct)); 260 PetscCall(MatTranspose(Dt, MAT_INPLACE_MATRIX, &Dt)); 261 PetscCall(MatDenseGetArrayRead(Ct, &c_ptr)); 262 PetscCall(MatDenseGetArrayWrite(Dt, &d_ptr)); 263 PetscCall(VecCreateMPIWithArray(PETSC_COMM_WORLD, 1, AM * N[k], PETSC_DECIDE, c_ptr, &cC)); 264 PetscCall(VecCreateMPIWithArray(PETSC_COMM_WORLD, 1, AM * N[k], PETSC_DECIDE, d_ptr, &cD)); 265 PetscCall(MatMult(E, cC, cD)); 266 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Benchmarking MatMult: with A %s %" PetscInt_FMT "x%" PetscInt_FMT " and B %s %" PetscInt_FMT "x%" PetscInt_FMT "\n", MATMAIJ, AM, AN, VECMPI, AM * N[k], 1)); 267 PetscCall(PetscLogStagePush(stage)); 268 for (t = 0; t < trial; ++t) PetscCall(MatMult(E, cC, cD)); 269 PetscCall(PetscLogStagePop()); 270 PetscCall(VecDestroy(&cD)); 271 PetscCall(VecDestroy(&cC)); 272 PetscCall(MatDestroy(&E)); 273 PetscCall(MatDenseRestoreArrayWrite(Dt, &d_ptr)); 274 PetscCall(MatDenseRestoreArrayRead(Ct, &c_ptr)); 275 PetscCall(MatTranspose(Ct, MAT_INPLACE_MATRIX, &Ct)); 276 PetscCall(MatTranspose(Dt, MAT_INPLACE_MATRIX, &Dt)); 277 } 278 } else { 279 const PetscScalar *c_ptr; 280 PetscScalar *d_ptr; 281 282 PetscCall(MatDenseGetArrayRead(C, &c_ptr)); 283 PetscCall(MatDenseGetArrayWrite(D, &d_ptr)); 284 PetscCallMKLSparse(mkl_sparse_d_mm, (SPARSE_OPERATION_NON_TRANSPOSE, 1.0, spr, descr, SPARSE_LAYOUT_COLUMN_MAJOR, c_ptr, CN, CM, 0.0, d_ptr, CM)); 285 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Benchmarking mkl_sparse_d_mm (COLUMN_MAJOR): with A %s %" PetscInt_FMT "x%" PetscInt_FMT " and B %s %" PetscInt_FMT "x%" PetscInt_FMT "\n", Atype, AM, AN, Ctype, CM, CN)); 286 PetscCall(PetscLogStagePush(stage)); 287 for (t = 0; t < trial; ++t) PetscCallMKLSparse(mkl_sparse_d_mm, (SPARSE_OPERATION_NON_TRANSPOSE, 1.0, spr, descr, SPARSE_LAYOUT_COLUMN_MAJOR, c_ptr, CN, CM, 0.0, d_ptr, CM)); 288 PetscCall(PetscLogStagePop()); 289 PetscCall(MatDenseRestoreArrayWrite(D, &d_ptr)); 290 PetscCall(MatDenseRestoreArrayRead(C, &c_ptr)); 291 } 292 } else if (maij) { 293 PetscCall(MatDestroy(&C)); 294 PetscCall(MatDestroy(&D)); 295 continue; 296 } else if (!mkl) { 297 Vec cC, cD; 298 299 PetscCall(MatDenseGetColumnVecRead(C, 0, &cC)); 300 PetscCall(MatDenseGetColumnVecWrite(D, 0, &cD)); 301 PetscCall(MatMult(A, cC, cD)); 302 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Benchmarking MatMult: with A %s %" PetscInt_FMT "x%" PetscInt_FMT "\n", Atype, AM, AN)); 303 PetscCall(PetscLogStagePush(stage)); 304 for (t = 0; t < trial; ++t) PetscCall(MatMult(A, cC, cD)); 305 PetscCall(PetscLogStagePop()); 306 PetscCall(MatDenseRestoreColumnVecRead(C, 0, &cC)); 307 PetscCall(MatDenseRestoreColumnVecWrite(D, 0, &cD)); 308 } else { 309 const PetscScalar *c_ptr; 310 PetscScalar *d_ptr; 311 312 PetscCall(MatDenseGetArrayRead(C, &c_ptr)); 313 PetscCall(MatDenseGetArrayWrite(D, &d_ptr)); 314 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Benchmarking mkl_sparse_d_mv: with A %s %" PetscInt_FMT "x%" PetscInt_FMT "\n", Atype, AM, AN)); 315 PetscCallMKLSparse(mkl_sparse_d_mv, (SPARSE_OPERATION_NON_TRANSPOSE, 1.0, spr, descr, c_ptr, 0.0, d_ptr)); 316 PetscCall(PetscLogStagePush(stage)); 317 for (t = 0; t < trial; ++t) PetscCallMKLSparse(mkl_sparse_d_mv, (SPARSE_OPERATION_NON_TRANSPOSE, 1.0, spr, descr, c_ptr, 0.0, d_ptr)); 318 PetscCall(PetscLogStagePop()); 319 PetscCall(MatDenseRestoreArrayWrite(D, &d_ptr)); 320 PetscCall(MatDenseRestoreArrayRead(C, &c_ptr)); 321 } 322 323 if (check) { 324 PetscCall(MatMatMultEqual(A, C, D, 10, &flg)); 325 if (!flg) { 326 MatType Dtype; 327 328 PetscCall(MatGetType(D, &Dtype)); 329 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error with A %s%s, C %s, D %s, Nk %" PetscInt_FMT "\n", Atype, mkl ? "mkl" : "", Ctype, Dtype, N[k])); 330 } 331 } 332 333 /* MKL implementation seems buggy for ABt */ 334 /* A*Bt */ 335 if (!mkl && trans && N[k] > 1) { 336 PetscCall(PetscObjectTypeCompareAny((PetscObject)A, &flg, MATSEQAIJ, MATMPIAIJ, "")); 337 if (flg) { 338 PetscCall(MatTranspose(C, MAT_INPLACE_MATRIX, &C)); 339 PetscCall(MatGetType(C, &Ctype)); 340 if (!mkl) { 341 PetscCall(MatProductCreateWithMat(A, C, NULL, D)); 342 PetscCall(MatProductSetType(D, MATPRODUCT_ABt)); 343 PetscCall(MatProductSetFromOptions(D)); 344 PetscCall(MatProductSymbolic(D)); 345 PetscCall(MatProductNumeric(D)); 346 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Benchmarking MatProduct %s: with A %s %" PetscInt_FMT "x%" PetscInt_FMT " and Bt %s %" PetscInt_FMT "x%" PetscInt_FMT "\n", MatProductTypes[MATPRODUCT_ABt], Atype, AM, AN, Ctype, CM, CN)); 347 PetscCall(PetscLogStagePush(tstage)); 348 for (t = 0; t < trial; ++t) PetscCall(MatProductNumeric(D)); 349 PetscCall(PetscLogStagePop()); 350 } else { 351 const PetscScalar *c_ptr; 352 PetscScalar *d_ptr; 353 354 PetscCallMKLSparse(mkl_sparse_set_mm_hint, (spr, SPARSE_OPERATION_NON_TRANSPOSE, descr, SPARSE_LAYOUT_ROW_MAJOR, N[k], 1 + trial)); 355 PetscCallMKLSparse(mkl_sparse_optimize, (spr)); 356 PetscCall(MatDenseGetArrayRead(C, &c_ptr)); 357 PetscCall(MatDenseGetArrayWrite(D, &d_ptr)); 358 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Benchmarking mkl_sparse_d_mm (ROW_MAJOR): with A %s %" PetscInt_FMT "x%" PetscInt_FMT " and B %s %" PetscInt_FMT "x%" PetscInt_FMT "\n", Atype, AM, AN, Ctype, CM, CN)); 359 PetscCallMKLSparse(mkl_sparse_d_mm, (SPARSE_OPERATION_NON_TRANSPOSE, 1.0, spr, descr, SPARSE_LAYOUT_ROW_MAJOR, c_ptr, CN, CM, 0.0, d_ptr, CM)); 360 PetscCall(PetscLogStagePush(stage)); 361 for (t = 0; t < trial; ++t) PetscCallMKLSparse(mkl_sparse_d_mm, (SPARSE_OPERATION_NON_TRANSPOSE, 1.0, spr, descr, SPARSE_LAYOUT_ROW_MAJOR, c_ptr, CN, CM, 0.0, d_ptr, CM)); 362 PetscCall(PetscLogStagePop()); 363 PetscCall(MatDenseRestoreArrayWrite(D, &d_ptr)); 364 PetscCall(MatDenseRestoreArrayRead(C, &c_ptr)); 365 } 366 } 367 } 368 369 if (!mkl && trans && N[k] > 1 && flg && check) { 370 PetscCall(MatMatTransposeMultEqual(A, C, D, 10, &flg)); 371 if (!flg) { 372 MatType Dtype; 373 PetscCall(MatGetType(D, &Dtype)); 374 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error with A %s%s, C %s, D %s, Nk %" PetscInt_FMT "\n", Atype, mkl ? "mkl" : "", Ctype, Dtype, N[k])); 375 } 376 } 377 PetscCall(MatDestroy(&C)); 378 PetscCall(MatDestroy(&D)); 379 } 380 if (mkl) { 381 PetscCallMKLSparse(mkl_sparse_destroy, (spr)); 382 PetscCall(PetscFree(ia_ptr)); 383 PetscCall(PetscFree(ja_ptr)); 384 } 385 if (cuda && i != ntype - 1) { 386 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "AIJCUSPARSE must be last, otherwise MatConvert() to another MatType is too slow\n")); 387 break; 388 } 389 } 390 if (E != A) PetscCall(MatDestroy(&E)); 391 PetscCall(MatDestroy(&A)); 392 } 393 for (m = 0; m < ntype; ++m) PetscCall(PetscFree(type[m])); 394 PetscCall(PetscFinalize()); 395 return 0; 396 } 397 398 /*TEST 399 build: 400 requires: double !complex !defined(PETSC_USE_64BIT_INDICES) 401 402 testset: 403 nsize: 1 404 filter: sed "/Benchmarking/d" 405 args: -f ${wPETSC_DIR}/share/petsc/datafiles/matrices/spd-real-int32-float64 -bs 1,2,3 -N 1,2,18 -check -trans -convert_aij {{false true}shared output} 406 test: 407 suffix: basic 408 args: -type aij,sbaij,baij 409 output_file: output/ex237.out 410 test: 411 suffix: maij 412 args: -type aij,maij 413 output_file: output/ex237.out 414 test: 415 suffix: cuda 416 requires: cuda 417 args: -type aij,aijcusparse 418 output_file: output/ex237.out 419 test: 420 suffix: mkl 421 requires: mkl_sparse_optimize 422 args: -type aij,aijmkl,baijmkl,sbaijmkl 423 output_file: output/ex237.out 424 425 TEST*/ 426