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