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