1 2 static char help[] = "Tests sequential and parallel MatMatMult() and MatPtAP(), MatTransposeMatMult(), sequential MatMatTransposeMult(), MatRARt()\n\ 3 Input arguments are:\n\ 4 -f0 <input_file> -f1 <input_file> -f2 <input_file> -f3 <input_file> : file to load\n\n"; 5 /* Example of usage: 6 ./ex94 -f0 <A_binary> -f1 <B_binary> -matmatmult_mat_view ascii::ascii_info -matmatmulttr_mat_view 7 mpiexec -n 3 ./ex94 -f0 medium -f1 medium -f2 arco1 -f3 arco1 -matmatmult_mat_view 8 */ 9 10 #include <petscmat.h> 11 12 /* 13 B = A - B 14 norm = norm(B) 15 */ 16 PetscErrorCode MatNormDifference(Mat A, Mat B, PetscReal *norm) 17 { 18 PetscFunctionBegin; 19 PetscCall(MatAXPY(B, -1.0, A, DIFFERENT_NONZERO_PATTERN)); 20 PetscCall(MatNorm(B, NORM_FROBENIUS, norm)); 21 PetscFunctionReturn(PETSC_SUCCESS); 22 } 23 24 int main(int argc, char **args) 25 { 26 Mat A, A_save, B, AT, ATT, BT, BTT, P, R, C, C1; 27 Vec x, v1, v2, v3, v4; 28 PetscViewer viewer; 29 PetscMPIInt size, rank; 30 PetscInt i, m, n, j, *idxn, M, N, nzp, rstart, rend; 31 PetscReal norm, norm_abs, norm_tmp, fill = 4.0; 32 PetscRandom rdm; 33 char file[4][128]; 34 PetscBool flg, preload = PETSC_TRUE; 35 PetscScalar *a, rval, alpha, none = -1.0; 36 PetscBool Test_MatMatMult = PETSC_TRUE, Test_MatMatTr = PETSC_TRUE, Test_MatPtAP = PETSC_TRUE, Test_MatRARt = PETSC_TRUE, Test_MatMatMatMult = PETSC_TRUE; 37 PetscBool Test_MatAXPY = PETSC_FALSE, view = PETSC_FALSE; 38 PetscInt pm, pn, pM, pN; 39 MatInfo info; 40 PetscBool seqaij; 41 MatType mattype; 42 Mat Cdensetest, Pdense, Cdense, Adense; 43 44 PetscFunctionBeginUser; 45 PetscCall(PetscInitialize(&argc, &args, (char *)0, help)); 46 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 47 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank)); 48 49 PetscCall(PetscOptionsGetReal(NULL, NULL, "-fill", &fill, NULL)); 50 PetscCall(PetscOptionsGetBool(NULL, NULL, "-matops_view", &view, NULL)); 51 if (view) PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_INFO)); 52 53 /* Load the matrices A_save and B */ 54 PetscCall(PetscOptionsGetString(NULL, NULL, "-f0", file[0], sizeof(file[0]), &flg)); 55 PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_USER, "Must indicate a file name for small matrix A with the -f0 option."); 56 PetscCall(PetscOptionsGetString(NULL, NULL, "-f1", file[1], sizeof(file[1]), &flg)); 57 PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_USER, "Must indicate a file name for small matrix B with the -f1 option."); 58 PetscCall(PetscOptionsGetString(NULL, NULL, "-f2", file[2], sizeof(file[2]), &flg)); 59 if (!flg) { 60 preload = PETSC_FALSE; 61 } else { 62 PetscCall(PetscOptionsGetString(NULL, NULL, "-f3", file[3], sizeof(file[3]), &flg)); 63 PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_USER, "Must indicate a file name for test matrix B with the -f3 option."); 64 } 65 66 PetscPreLoadBegin(preload, "Load system"); 67 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, file[2 * PetscPreLoadIt], FILE_MODE_READ, &viewer)); 68 PetscCall(MatCreate(PETSC_COMM_WORLD, &A_save)); 69 PetscCall(MatSetFromOptions(A_save)); 70 PetscCall(MatLoad(A_save, viewer)); 71 PetscCall(PetscViewerDestroy(&viewer)); 72 73 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, file[2 * PetscPreLoadIt + 1], FILE_MODE_READ, &viewer)); 74 PetscCall(MatCreate(PETSC_COMM_WORLD, &B)); 75 PetscCall(MatSetFromOptions(B)); 76 PetscCall(MatLoad(B, viewer)); 77 PetscCall(PetscViewerDestroy(&viewer)); 78 79 PetscCall(MatGetType(B, &mattype)); 80 81 PetscCall(MatGetSize(B, &M, &N)); 82 nzp = PetscMax((PetscInt)(0.1 * M), 5); 83 PetscCall(PetscMalloc((nzp + 1) * (sizeof(PetscInt) + sizeof(PetscScalar)), &idxn)); 84 a = (PetscScalar *)(idxn + nzp); 85 86 /* Create vectors v1 and v2 that are compatible with A_save */ 87 PetscCall(VecCreate(PETSC_COMM_WORLD, &v1)); 88 PetscCall(MatGetLocalSize(A_save, &m, NULL)); 89 PetscCall(VecSetSizes(v1, m, PETSC_DECIDE)); 90 PetscCall(VecSetFromOptions(v1)); 91 PetscCall(VecDuplicate(v1, &v2)); 92 93 PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rdm)); 94 PetscCall(PetscRandomSetFromOptions(rdm)); 95 PetscCall(PetscOptionsGetReal(NULL, NULL, "-fill", &fill, NULL)); 96 97 /* Test MatAXPY() */ 98 /*-------------------*/ 99 PetscCall(PetscOptionsHasName(NULL, NULL, "-test_MatAXPY", &Test_MatAXPY)); 100 if (Test_MatAXPY) { 101 Mat Btmp; 102 PetscCall(MatDuplicate(A_save, MAT_COPY_VALUES, &A)); 103 PetscCall(MatDuplicate(B, MAT_COPY_VALUES, &Btmp)); 104 PetscCall(MatAXPY(A, -1.0, B, DIFFERENT_NONZERO_PATTERN)); /* A = -B + A_save */ 105 106 PetscCall(MatScale(A, -1.0)); /* A = -A = B - A_save */ 107 PetscCall(MatAXPY(Btmp, -1.0, A, DIFFERENT_NONZERO_PATTERN)); /* Btmp = -A + B = A_save */ 108 PetscCall(MatMultEqual(A_save, Btmp, 10, &flg)); 109 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatAXPY() is incorrect"); 110 PetscCall(MatDestroy(&A)); 111 PetscCall(MatDestroy(&Btmp)); 112 113 Test_MatMatMult = PETSC_FALSE; 114 Test_MatMatTr = PETSC_FALSE; 115 Test_MatPtAP = PETSC_FALSE; 116 Test_MatRARt = PETSC_FALSE; 117 Test_MatMatMatMult = PETSC_FALSE; 118 } 119 120 /* 1) Test MatMatMult() */ 121 /* ---------------------*/ 122 if (Test_MatMatMult) { 123 PetscCall(MatDuplicate(A_save, MAT_COPY_VALUES, &A)); 124 PetscCall(MatCreateTranspose(A, &AT)); 125 PetscCall(MatCreateTranspose(AT, &ATT)); 126 PetscCall(MatCreateTranspose(B, &BT)); 127 PetscCall(MatCreateTranspose(BT, &BTT)); 128 129 PetscCall(MatMatMult(AT, B, MAT_INITIAL_MATRIX, fill, &C)); 130 PetscCall(MatMatMultEqual(AT, B, C, 10, &flg)); 131 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Error: MatMatMult() for C=AT*B"); 132 PetscCall(MatDestroy(&C)); 133 134 PetscCall(MatMatMult(ATT, B, MAT_INITIAL_MATRIX, fill, &C)); 135 PetscCall(MatMatMultEqual(ATT, B, C, 10, &flg)); 136 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Error: MatMatMult() for C=ATT*B"); 137 PetscCall(MatDestroy(&C)); 138 139 PetscCall(MatMatMult(A, B, MAT_INITIAL_MATRIX, fill, &C)); 140 PetscCall(MatMatMultEqual(A, B, C, 10, &flg)); 141 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Error: MatMatMult() for reuse C=A*B"); 142 /* ATT has different matrix type as A (although they have same internal data structure), 143 we cannot call MatProductReplaceMats(ATT,NULL,NULL,C) and MatMatMult(ATT,B,MAT_REUSE_MATRIX,fill,&C) */ 144 PetscCall(MatDestroy(&C)); 145 146 PetscCall(MatMatMult(A, BTT, MAT_INITIAL_MATRIX, fill, &C)); 147 PetscCall(MatMatMultEqual(A, BTT, C, 10, &flg)); 148 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Error: MatMatMult() for C=A*BTT"); 149 PetscCall(MatDestroy(&C)); 150 151 PetscCall(MatMatMult(ATT, BTT, MAT_INITIAL_MATRIX, fill, &C)); 152 PetscCall(MatMatMultEqual(A, B, C, 10, &flg)); 153 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Error: MatMatMult()"); 154 PetscCall(MatDestroy(&C)); 155 156 PetscCall(MatDestroy(&BTT)); 157 PetscCall(MatDestroy(&BT)); 158 PetscCall(MatDestroy(&ATT)); 159 PetscCall(MatDestroy(&AT)); 160 161 PetscCall(MatMatMult(A, B, MAT_INITIAL_MATRIX, fill, &C)); 162 PetscCall(MatSetOptionsPrefix(C, "matmatmult_")); /* enable option '-matmatmult_' for matrix C */ 163 PetscCall(MatGetInfo(C, MAT_GLOBAL_SUM, &info)); 164 165 /* Test MAT_REUSE_MATRIX - reuse symbolic C */ 166 alpha = 1.0; 167 for (i = 0; i < 2; i++) { 168 alpha -= 0.1; 169 PetscCall(MatScale(A, alpha)); 170 PetscCall(MatMatMult(A, B, MAT_REUSE_MATRIX, fill, &C)); 171 } 172 PetscCall(MatMatMultEqual(A, B, C, 10, &flg)); 173 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Error: MatMatMult()"); 174 PetscCall(MatDestroy(&A)); 175 176 /* Test MatDuplicate() of C=A*B */ 177 PetscCall(MatDuplicate(C, MAT_COPY_VALUES, &C1)); 178 PetscCall(MatDestroy(&C1)); 179 PetscCall(MatDestroy(&C)); 180 } /* if (Test_MatMatMult) */ 181 182 /* 2) Test MatTransposeMatMult() and MatMatTransposeMult() */ 183 /* ------------------------------------------------------- */ 184 if (Test_MatMatTr) { 185 /* Create P */ 186 PetscInt PN, rstart, rend; 187 PN = M / 2; 188 nzp = 5; /* num of nonzeros in each row of P */ 189 PetscCall(MatCreate(PETSC_COMM_WORLD, &P)); 190 PetscCall(MatSetSizes(P, PETSC_DECIDE, PETSC_DECIDE, M, PN)); 191 PetscCall(MatSetType(P, mattype)); 192 PetscCall(MatSeqAIJSetPreallocation(P, nzp, NULL)); 193 PetscCall(MatMPIAIJSetPreallocation(P, nzp, NULL, nzp, NULL)); 194 PetscCall(MatGetOwnershipRange(P, &rstart, &rend)); 195 for (i = 0; i < nzp; i++) PetscCall(PetscRandomGetValue(rdm, &a[i])); 196 for (i = rstart; i < rend; i++) { 197 for (j = 0; j < nzp; j++) { 198 PetscCall(PetscRandomGetValue(rdm, &rval)); 199 idxn[j] = (PetscInt)(PetscRealPart(rval) * PN); 200 } 201 PetscCall(MatSetValues(P, 1, &i, nzp, idxn, a, ADD_VALUES)); 202 } 203 PetscCall(MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY)); 204 PetscCall(MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY)); 205 206 /* Create R = P^T */ 207 PetscCall(MatTranspose(P, MAT_INITIAL_MATRIX, &R)); 208 209 { /* Test R = P^T, C1 = R*B */ 210 PetscCall(MatMatMult(R, B, MAT_INITIAL_MATRIX, fill, &C1)); 211 PetscCall(MatTranspose(P, MAT_REUSE_MATRIX, &R)); 212 PetscCall(MatMatMult(R, B, MAT_REUSE_MATRIX, fill, &C1)); 213 PetscCall(MatDestroy(&C1)); 214 } 215 216 /* C = P^T*B */ 217 PetscCall(MatTransposeMatMult(P, B, MAT_INITIAL_MATRIX, fill, &C)); 218 PetscCall(MatGetInfo(C, MAT_GLOBAL_SUM, &info)); 219 220 /* Test MAT_REUSE_MATRIX - reuse symbolic C */ 221 PetscCall(MatTransposeMatMult(P, B, MAT_REUSE_MATRIX, fill, &C)); 222 if (view) { 223 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "C = P^T * B:\n")); 224 PetscCall(MatView(C, PETSC_VIEWER_STDOUT_WORLD)); 225 } 226 PetscCall(MatProductClear(C)); 227 if (view) { 228 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nC = P^T * B after MatProductClear():\n")); 229 PetscCall(MatView(C, PETSC_VIEWER_STDOUT_WORLD)); 230 } 231 232 /* Compare P^T*B and R*B */ 233 PetscCall(MatMatMult(R, B, MAT_INITIAL_MATRIX, fill, &C1)); 234 PetscCall(MatNormDifference(C, C1, &norm)); 235 PetscCheck(norm <= PETSC_SMALL, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Error in MatTransposeMatMult(): %g", (double)norm); 236 PetscCall(MatDestroy(&C1)); 237 238 /* Test MatDuplicate() of C=P^T*B */ 239 PetscCall(MatDuplicate(C, MAT_COPY_VALUES, &C1)); 240 PetscCall(MatDestroy(&C1)); 241 PetscCall(MatDestroy(&C)); 242 243 /* C = B*R^T */ 244 PetscCall(PetscObjectTypeCompare((PetscObject)B, MATSEQAIJ, &seqaij)); 245 if (size == 1 && seqaij) { 246 PetscCall(MatMatTransposeMult(B, R, MAT_INITIAL_MATRIX, fill, &C)); 247 PetscCall(MatSetOptionsPrefix(C, "matmatmulttr_")); /* enable '-matmatmulttr_' for matrix C */ 248 PetscCall(MatGetInfo(C, MAT_GLOBAL_SUM, &info)); 249 250 /* Test MAT_REUSE_MATRIX - reuse symbolic C */ 251 PetscCall(MatMatTransposeMult(B, R, MAT_REUSE_MATRIX, fill, &C)); 252 253 /* Check */ 254 PetscCall(MatMatMult(B, P, MAT_INITIAL_MATRIX, fill, &C1)); 255 PetscCall(MatNormDifference(C, C1, &norm)); 256 PetscCheck(norm <= PETSC_SMALL, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Error in MatMatTransposeMult() %g", (double)norm); 257 PetscCall(MatDestroy(&C1)); 258 PetscCall(MatDestroy(&C)); 259 } 260 PetscCall(MatDestroy(&P)); 261 PetscCall(MatDestroy(&R)); 262 } 263 264 /* 3) Test MatPtAP() */ 265 /*-------------------*/ 266 if (Test_MatPtAP) { 267 PetscInt PN; 268 Mat Cdup; 269 270 PetscCall(MatDuplicate(A_save, MAT_COPY_VALUES, &A)); 271 PetscCall(MatGetSize(A, &M, &N)); 272 PetscCall(MatGetLocalSize(A, &m, &n)); 273 274 PN = M / 2; 275 nzp = (PetscInt)(0.1 * PN + 1); /* num of nonzeros in each row of P */ 276 PetscCall(MatCreate(PETSC_COMM_WORLD, &P)); 277 PetscCall(MatSetSizes(P, PETSC_DECIDE, PETSC_DECIDE, N, PN)); 278 PetscCall(MatSetType(P, mattype)); 279 PetscCall(MatSeqAIJSetPreallocation(P, nzp, NULL)); 280 PetscCall(MatMPIAIJSetPreallocation(P, nzp, NULL, nzp, NULL)); 281 for (i = 0; i < nzp; i++) PetscCall(PetscRandomGetValue(rdm, &a[i])); 282 PetscCall(MatGetOwnershipRange(P, &rstart, &rend)); 283 for (i = rstart; i < rend; i++) { 284 for (j = 0; j < nzp; j++) { 285 PetscCall(PetscRandomGetValue(rdm, &rval)); 286 idxn[j] = (PetscInt)(PetscRealPart(rval) * PN); 287 } 288 PetscCall(MatSetValues(P, 1, &i, nzp, idxn, a, ADD_VALUES)); 289 } 290 PetscCall(MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY)); 291 PetscCall(MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY)); 292 293 /* PetscCall(MatView(P,PETSC_VIEWER_STDOUT_WORLD)); */ 294 PetscCall(MatGetSize(P, &pM, &pN)); 295 PetscCall(MatGetLocalSize(P, &pm, &pn)); 296 PetscCall(MatPtAP(A, P, MAT_INITIAL_MATRIX, fill, &C)); 297 298 /* Test MAT_REUSE_MATRIX - reuse symbolic C */ 299 alpha = 1.0; 300 for (i = 0; i < 2; i++) { 301 alpha -= 0.1; 302 PetscCall(MatScale(A, alpha)); 303 PetscCall(MatPtAP(A, P, MAT_REUSE_MATRIX, fill, &C)); 304 } 305 306 /* Test PtAP ops with P Dense and A either AIJ or SeqDense (it assumes MatPtAP_XAIJ_XAIJ is fine) */ 307 PetscCall(PetscObjectTypeCompare((PetscObject)A, MATSEQAIJ, &seqaij)); 308 if (seqaij) { 309 PetscCall(MatConvert(C, MATSEQDENSE, MAT_INITIAL_MATRIX, &Cdensetest)); 310 PetscCall(MatConvert(P, MATSEQDENSE, MAT_INITIAL_MATRIX, &Pdense)); 311 } else { 312 PetscCall(MatConvert(C, MATMPIDENSE, MAT_INITIAL_MATRIX, &Cdensetest)); 313 PetscCall(MatConvert(P, MATMPIDENSE, MAT_INITIAL_MATRIX, &Pdense)); 314 } 315 316 /* test with A(AIJ), Pdense -- call MatPtAP_Basic() when np>1 */ 317 PetscCall(MatPtAP(A, Pdense, MAT_INITIAL_MATRIX, fill, &Cdense)); 318 PetscCall(MatPtAP(A, Pdense, MAT_REUSE_MATRIX, fill, &Cdense)); 319 PetscCall(MatPtAPMultEqual(A, Pdense, Cdense, 10, &flg)); 320 PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Error in MatPtAP with A AIJ and P Dense"); 321 PetscCall(MatDestroy(&Cdense)); 322 323 /* test with A SeqDense */ 324 if (seqaij) { 325 PetscCall(MatConvert(A, MATSEQDENSE, MAT_INITIAL_MATRIX, &Adense)); 326 PetscCall(MatPtAP(Adense, Pdense, MAT_INITIAL_MATRIX, fill, &Cdense)); 327 PetscCall(MatPtAP(Adense, Pdense, MAT_REUSE_MATRIX, fill, &Cdense)); 328 PetscCall(MatPtAPMultEqual(Adense, Pdense, Cdense, 10, &flg)); 329 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Error in MatPtAP with A SeqDense and P SeqDense"); 330 PetscCall(MatDestroy(&Cdense)); 331 PetscCall(MatDestroy(&Adense)); 332 } 333 PetscCall(MatDestroy(&Cdensetest)); 334 PetscCall(MatDestroy(&Pdense)); 335 336 /* Test MatDuplicate() of C=PtAP and MatView(Cdup,...) */ 337 PetscCall(MatDuplicate(C, MAT_COPY_VALUES, &Cdup)); 338 if (view) { 339 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nC = P^T * A * P:\n")); 340 PetscCall(MatView(C, PETSC_VIEWER_STDOUT_WORLD)); 341 342 PetscCall(MatProductClear(C)); 343 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nC = P^T * A * P after MatProductClear():\n")); 344 PetscCall(MatView(C, PETSC_VIEWER_STDOUT_WORLD)); 345 346 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nCdup:\n")); 347 PetscCall(MatView(Cdup, PETSC_VIEWER_STDOUT_WORLD)); 348 } 349 PetscCall(MatDestroy(&Cdup)); 350 351 if (size > 1 || !seqaij) Test_MatRARt = PETSC_FALSE; 352 /* 4) Test MatRARt() */ 353 /* ----------------- */ 354 if (Test_MatRARt) { 355 Mat R, RARt, Rdense, RARtdense; 356 PetscCall(MatTranspose(P, MAT_INITIAL_MATRIX, &R)); 357 358 /* Test MatRARt_Basic(), MatMatMatMult_Basic() */ 359 PetscCall(MatConvert(R, MATDENSE, MAT_INITIAL_MATRIX, &Rdense)); 360 PetscCall(MatRARt(A, Rdense, MAT_INITIAL_MATRIX, 2.0, &RARtdense)); 361 PetscCall(MatRARt(A, Rdense, MAT_REUSE_MATRIX, 2.0, &RARtdense)); 362 363 PetscCall(MatConvert(RARtdense, MATAIJ, MAT_INITIAL_MATRIX, &RARt)); 364 PetscCall(MatNormDifference(C, RARt, &norm)); 365 PetscCheck(norm <= PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_PLIB, "|PtAP - RARtdense| = %g", (double)norm); 366 PetscCall(MatDestroy(&Rdense)); 367 PetscCall(MatDestroy(&RARtdense)); 368 PetscCall(MatDestroy(&RARt)); 369 370 /* Test MatRARt() for aij matrices */ 371 PetscCall(MatRARt(A, R, MAT_INITIAL_MATRIX, 2.0, &RARt)); 372 PetscCall(MatRARt(A, R, MAT_REUSE_MATRIX, 2.0, &RARt)); 373 PetscCall(MatNormDifference(C, RARt, &norm)); 374 PetscCheck(norm <= PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_PLIB, "|PtAP - RARt| = %g", (double)norm); 375 PetscCall(MatDestroy(&R)); 376 PetscCall(MatDestroy(&RARt)); 377 } 378 379 if (Test_MatMatMatMult && size == 1) { 380 Mat R, RAP; 381 PetscCall(MatTranspose(P, MAT_INITIAL_MATRIX, &R)); 382 PetscCall(MatMatMatMult(R, A, P, MAT_INITIAL_MATRIX, 2.0, &RAP)); 383 PetscCall(MatMatMatMult(R, A, P, MAT_REUSE_MATRIX, 2.0, &RAP)); 384 PetscCall(MatNormDifference(C, RAP, &norm)); 385 PetscCheck(norm <= PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_PLIB, "PtAP != RAP %g", (double)norm); 386 PetscCall(MatDestroy(&R)); 387 PetscCall(MatDestroy(&RAP)); 388 } 389 390 /* Create vector x that is compatible with P */ 391 PetscCall(VecCreate(PETSC_COMM_WORLD, &x)); 392 PetscCall(MatGetLocalSize(P, &m, &n)); 393 PetscCall(VecSetSizes(x, n, PETSC_DECIDE)); 394 PetscCall(VecSetFromOptions(x)); 395 396 PetscCall(VecCreate(PETSC_COMM_WORLD, &v3)); 397 PetscCall(VecSetSizes(v3, n, PETSC_DECIDE)); 398 PetscCall(VecSetFromOptions(v3)); 399 PetscCall(VecDuplicate(v3, &v4)); 400 401 norm = 0.0; 402 for (i = 0; i < 10; i++) { 403 PetscCall(VecSetRandom(x, rdm)); 404 PetscCall(MatMult(P, x, v1)); 405 PetscCall(MatMult(A, v1, v2)); /* v2 = A*P*x */ 406 407 PetscCall(MatMultTranspose(P, v2, v3)); /* v3 = Pt*A*P*x */ 408 PetscCall(MatMult(C, x, v4)); /* v3 = C*x */ 409 PetscCall(VecNorm(v4, NORM_2, &norm_abs)); 410 PetscCall(VecAXPY(v4, none, v3)); 411 PetscCall(VecNorm(v4, NORM_2, &norm_tmp)); 412 413 norm_tmp /= norm_abs; 414 if (norm_tmp > norm) norm = norm_tmp; 415 } 416 PetscCheck(norm < PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Error: MatPtAP(), |v1 - v2|: %g", (double)norm); 417 418 PetscCall(MatDestroy(&A)); 419 PetscCall(MatDestroy(&P)); 420 PetscCall(MatDestroy(&C)); 421 PetscCall(VecDestroy(&v3)); 422 PetscCall(VecDestroy(&v4)); 423 PetscCall(VecDestroy(&x)); 424 } 425 426 /* Destroy objects */ 427 PetscCall(VecDestroy(&v1)); 428 PetscCall(VecDestroy(&v2)); 429 PetscCall(PetscRandomDestroy(&rdm)); 430 PetscCall(PetscFree(idxn)); 431 432 PetscCall(MatDestroy(&A_save)); 433 PetscCall(MatDestroy(&B)); 434 435 PetscPreLoadEnd(); 436 PetscCall(PetscFinalize()); 437 return 0; 438 } 439 440 /*TEST 441 442 test: 443 suffix: 2_mattransposematmult_matmatmult 444 nsize: 3 445 requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES) 446 args: -f0 ${DATAFILESPATH}/matrices/medium -f1 ${DATAFILESPATH}/matrices/medium -mattransposematmult_via at*b> ex94_2.tmp 2>&1 447 448 test: 449 suffix: 2_mattransposematmult_scalable 450 nsize: 3 451 requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES) 452 args: -f0 ${DATAFILESPATH}/matrices/medium -f1 ${DATAFILESPATH}/matrices/medium -mattransposematmult_via scalable> ex94_2.tmp 2>&1 453 output_file: output/ex94_1.out 454 455 test: 456 suffix: axpy_mpiaij 457 requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES) 458 nsize: 8 459 args: -f0 ${DATAFILESPATH}/matrices/poisson_2d5p -f1 ${DATAFILESPATH}/matrices/poisson_2d13p -test_MatAXPY 460 output_file: output/ex94_1.out 461 462 test: 463 suffix: axpy_mpibaij 464 requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES) 465 nsize: 8 466 args: -f0 ${DATAFILESPATH}/matrices/poisson_2d5p -f1 ${DATAFILESPATH}/matrices/poisson_2d13p -test_MatAXPY -mat_type baij 467 output_file: output/ex94_1.out 468 469 test: 470 suffix: axpy_mpisbaij 471 requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES) 472 nsize: 8 473 args: -f0 ${DATAFILESPATH}/matrices/poisson_2d5p -f1 ${DATAFILESPATH}/matrices/poisson_2d13p -test_MatAXPY -mat_type sbaij 474 output_file: output/ex94_1.out 475 476 test: 477 suffix: matmatmult 478 requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES) 479 args: -f0 ${DATAFILESPATH}/matrices/arco1 -f1 ${DATAFILESPATH}/matrices/arco1 -viewer_binary_skip_info 480 output_file: output/ex94_1.out 481 482 test: 483 suffix: matmatmult_2 484 requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES) 485 args: -f0 ${DATAFILESPATH}/matrices/arco1 -f1 ${DATAFILESPATH}/matrices/arco1 -mat_type mpiaij -viewer_binary_skip_info 486 output_file: output/ex94_1.out 487 488 test: 489 suffix: matmatmult_scalable 490 nsize: 4 491 requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES) 492 args: -f0 ${DATAFILESPATH}/matrices/arco1 -f1 ${DATAFILESPATH}/matrices/arco1 -matmatmult_via scalable 493 output_file: output/ex94_1.out 494 495 test: 496 suffix: ptap 497 nsize: 3 498 requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES) 499 args: -f0 ${DATAFILESPATH}/matrices/medium -f1 ${DATAFILESPATH}/matrices/medium -matptap_via scalable 500 output_file: output/ex94_1.out 501 502 test: 503 suffix: rap 504 nsize: 3 505 requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES) 506 args: -f0 ${DATAFILESPATH}/matrices/medium -f1 ${DATAFILESPATH}/matrices/medium 507 output_file: output/ex94_1.out 508 509 test: 510 suffix: scalable0 511 requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES) 512 args: -f0 ${DATAFILESPATH}/matrices/arco1 -f1 ${DATAFILESPATH}/matrices/arco1 -viewer_binary_skip_info 513 output_file: output/ex94_1.out 514 515 test: 516 suffix: scalable1 517 requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES) 518 args: -f0 ${DATAFILESPATH}/matrices/arco1 -f1 ${DATAFILESPATH}/matrices/arco1 -viewer_binary_skip_info -matptap_via scalable 519 output_file: output/ex94_1.out 520 521 test: 522 suffix: view 523 nsize: 2 524 requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES) 525 args: -f0 ${DATAFILESPATH}/matrices/tiny -f1 ${DATAFILESPATH}/matrices/tiny -viewer_binary_skip_info -matops_view 526 output_file: output/ex94_2.out 527 528 TEST*/ 529