1 static char help[] = "Tests various routines for MATSHELL\n\n"; 2 3 #include <petscmat.h> 4 5 typedef struct _n_User *User; 6 struct _n_User { 7 Mat B; 8 }; 9 10 static PetscErrorCode MatGetDiagonal_User(Mat A, Vec X) 11 { 12 User user; 13 14 PetscFunctionBegin; 15 PetscCall(MatShellGetContext(A, &user)); 16 PetscCall(MatGetDiagonal(user->B, X)); 17 PetscFunctionReturn(PETSC_SUCCESS); 18 } 19 20 static PetscErrorCode MatMult_User(Mat A, Vec X, Vec Y) 21 { 22 User user; 23 24 PetscFunctionBegin; 25 PetscCall(MatShellGetContext(A, &user)); 26 PetscCall(MatMult(user->B, X, Y)); 27 PetscFunctionReturn(PETSC_SUCCESS); 28 } 29 30 static PetscErrorCode MatMultTranspose_User(Mat A, Vec X, Vec Y) 31 { 32 User user; 33 34 PetscFunctionBegin; 35 PetscCall(MatShellGetContext(A, &user)); 36 PetscCall(MatMultTranspose(user->B, X, Y)); 37 PetscFunctionReturn(PETSC_SUCCESS); 38 } 39 40 static PetscErrorCode MatCopy_User(Mat A, Mat X, MatStructure str) 41 { 42 User user, userX; 43 44 PetscFunctionBegin; 45 PetscCall(MatShellGetContext(A, &user)); 46 PetscCall(MatShellGetContext(X, &userX)); 47 PetscCheck(user == userX, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "This should not happen"); 48 PetscCall(PetscObjectReference((PetscObject)user->B)); 49 PetscFunctionReturn(PETSC_SUCCESS); 50 } 51 52 static PetscErrorCode MatDestroy_User(Mat A) 53 { 54 User user; 55 56 PetscFunctionBegin; 57 PetscCall(MatShellGetContext(A, &user)); 58 PetscCall(PetscObjectDereference((PetscObject)user->B)); 59 PetscFunctionReturn(PETSC_SUCCESS); 60 } 61 62 int main(int argc, char **args) 63 { 64 User user; 65 Mat A, S; 66 PetscScalar *data, diag = 1.3; 67 PetscReal tol = PETSC_SMALL; 68 PetscInt i, j, m = PETSC_DECIDE, n = PETSC_DECIDE, M = 17, N = 15, s1, s2; 69 PetscInt test, ntest = 2; 70 PetscMPIInt rank, size; 71 PetscBool nc = PETSC_FALSE, cong, flg; 72 PetscBool ronl = PETSC_TRUE; 73 PetscBool randomize = PETSC_FALSE, submat = PETSC_FALSE; 74 PetscBool keep = PETSC_FALSE; 75 PetscBool testzerorows = PETSC_TRUE, testdiagscale = PETSC_TRUE, testgetdiag = PETSC_TRUE, testsubmat = PETSC_TRUE; 76 PetscBool testshift = PETSC_TRUE, testscale = PETSC_TRUE, testdup = PETSC_TRUE, testreset = PETSC_TRUE; 77 PetscBool testaxpy = PETSC_TRUE, testaxpyd = PETSC_TRUE, testaxpyerr = PETSC_FALSE; 78 79 PetscFunctionBeginUser; 80 PetscCall(PetscInitialize(&argc, &args, NULL, help)); 81 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank)); 82 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 83 PetscCall(PetscOptionsGetInt(NULL, NULL, "-M", &M, NULL)); 84 PetscCall(PetscOptionsGetInt(NULL, NULL, "-N", &N, NULL)); 85 PetscCall(PetscOptionsGetInt(NULL, NULL, "-ml", &m, NULL)); 86 PetscCall(PetscOptionsGetInt(NULL, NULL, "-nl", &n, NULL)); 87 PetscCall(PetscOptionsGetBool(NULL, NULL, "-square_nc", &nc, NULL)); 88 PetscCall(PetscOptionsGetBool(NULL, NULL, "-rows_only", &ronl, NULL)); 89 PetscCall(PetscOptionsGetBool(NULL, NULL, "-randomize", &randomize, NULL)); 90 PetscCall(PetscOptionsGetBool(NULL, NULL, "-submat", &submat, NULL)); 91 PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_zerorows", &testzerorows, NULL)); 92 PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_diagscale", &testdiagscale, NULL)); 93 PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_getdiag", &testgetdiag, NULL)); 94 PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_shift", &testshift, NULL)); 95 PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_scale", &testscale, NULL)); 96 PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_dup", &testdup, NULL)); 97 PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_reset", &testreset, NULL)); 98 PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_submat", &testsubmat, NULL)); 99 PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_axpy", &testaxpy, NULL)); 100 PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_axpy_different", &testaxpyd, NULL)); 101 PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_axpy_error", &testaxpyerr, NULL)); 102 PetscCall(PetscOptionsGetInt(NULL, NULL, "-loop", &ntest, NULL)); 103 PetscCall(PetscOptionsGetReal(NULL, NULL, "-tol", &tol, NULL)); 104 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-diag", &diag, NULL)); 105 PetscCall(PetscOptionsGetBool(NULL, NULL, "-keep", &keep, NULL)); 106 /* This tests square matrices with different row/col layout */ 107 if (nc && size > 1) { 108 M = PetscMax(PetscMax(N, M), 1); 109 N = M; 110 m = n = 0; 111 if (rank == 0) { 112 m = M - 1; 113 n = 1; 114 } else if (rank == 1) { 115 m = 1; 116 n = N - 1; 117 } 118 } 119 PetscCall(MatCreateDense(PETSC_COMM_WORLD, m, n, M, N, NULL, &A)); 120 PetscCall(MatGetLocalSize(A, &m, &n)); 121 PetscCall(MatGetSize(A, &M, &N)); 122 PetscCall(MatGetOwnershipRange(A, &s1, NULL)); 123 s2 = 1; 124 while (s2 < M) s2 *= 10; 125 PetscCall(MatDenseGetArray(A, &data)); 126 for (j = 0; j < N; j++) { 127 #if defined(PETSC_USE_COMPLEX) 128 for (i = 0; i < m; i++) data[j * m + i] = s2 * j + i + s1 + 1 + PETSC_i * (s1 - 1); 129 #else 130 for (i = 0; i < m; i++) data[j * m + i] = s2 * j + i + s1 + 1; 131 #endif 132 } 133 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 134 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 135 136 if (submat) { 137 Mat A2; 138 IS r, c; 139 PetscInt rst, ren, cst, cen; 140 141 PetscCall(MatGetOwnershipRange(A, &rst, &ren)); 142 PetscCall(MatGetOwnershipRangeColumn(A, &cst, &cen)); 143 PetscCall(ISCreateStride(PetscObjectComm((PetscObject)A), (ren - rst) / 2, rst, 1, &r)); 144 PetscCall(ISCreateStride(PetscObjectComm((PetscObject)A), (cen - cst) / 2, cst, 1, &c)); 145 PetscCall(MatCreateSubMatrix(A, r, c, MAT_INITIAL_MATRIX, &A2)); 146 PetscCall(ISDestroy(&r)); 147 PetscCall(ISDestroy(&c)); 148 PetscCall(MatDestroy(&A)); 149 A = A2; 150 } 151 152 PetscCall(MatGetSize(A, &M, &N)); 153 PetscCall(MatGetLocalSize(A, &m, &n)); 154 PetscCall(MatHasCongruentLayouts(A, &cong)); 155 156 PetscCall(MatConvert(A, MATAIJ, MAT_INPLACE_MATRIX, &A)); 157 PetscCall(MatSetOption(A, MAT_KEEP_NONZERO_PATTERN, keep)); 158 PetscCall(PetscObjectSetName((PetscObject)A, "initial")); 159 PetscCall(MatViewFromOptions(A, NULL, "-view_mat")); 160 161 PetscCall(PetscNew(&user)); 162 PetscCall(MatCreateShell(PETSC_COMM_WORLD, m, n, M, N, user, &S)); 163 PetscCall(MatShellSetOperation(S, MATOP_MULT, (PetscErrorCodeFn *)MatMult_User)); 164 PetscCall(MatShellSetOperation(S, MATOP_MULT_TRANSPOSE, (PetscErrorCodeFn *)MatMultTranspose_User)); 165 if (cong) PetscCall(MatShellSetOperation(S, MATOP_GET_DIAGONAL, (PetscErrorCodeFn *)MatGetDiagonal_User)); 166 PetscCall(MatShellSetOperation(S, MATOP_COPY, (PetscErrorCodeFn *)MatCopy_User)); 167 PetscCall(MatShellSetOperation(S, MATOP_DESTROY, (PetscErrorCodeFn *)MatDestroy_User)); 168 PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &user->B)); 169 170 /* Square and rows only scaling */ 171 ronl = cong ? ronl : PETSC_TRUE; 172 173 for (test = 0; test < ntest; test++) { 174 PetscReal err; 175 176 PetscCall(MatMultAddEqual(A, S, 10, &flg)); 177 if (!flg) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[test %" PetscInt_FMT "] Error mult add\n", test)); 178 PetscCall(MatMultTransposeAddEqual(A, S, 10, &flg)); 179 if (!flg) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[test %" PetscInt_FMT "] Error mult transpose add\n", test)); 180 PetscCall(MatMultHermitianTransposeAddEqual(A, S, 10, &flg)); 181 if (!flg) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[test %" PetscInt_FMT "] Error mult hermitian transpose add\n", test)); 182 if (testzerorows) { 183 Mat ST, B, C, BT, BTT; 184 IS zr; 185 Vec x = NULL, b1 = NULL, b2 = NULL; 186 PetscInt *idxs = NULL, nr = 0; 187 188 if (rank == (test % size)) { 189 nr = 1; 190 PetscCall(PetscMalloc1(nr, &idxs)); 191 if (test % 2) { 192 idxs[0] = (2 * M - 1 - test / 2) % M; 193 } else { 194 idxs[0] = (test / 2) % M; 195 } 196 idxs[0] = PetscMax(idxs[0], 0); 197 } 198 PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, nr, idxs, PETSC_OWN_POINTER, &zr)); 199 PetscCall(PetscObjectSetName((PetscObject)zr, "ZR")); 200 PetscCall(ISViewFromOptions(zr, NULL, "-view_is")); 201 PetscCall(MatCreateVecs(A, &x, &b1)); 202 if (randomize) { 203 PetscCall(VecSetRandom(x, NULL)); 204 PetscCall(VecSetRandom(b1, NULL)); 205 } else { 206 PetscCall(VecSet(x, 11.4)); 207 PetscCall(VecSet(b1, -14.2)); 208 } 209 PetscCall(VecDuplicate(b1, &b2)); 210 PetscCall(VecCopy(b1, b2)); 211 PetscCall(PetscObjectSetName((PetscObject)b1, "A_B1")); 212 PetscCall(PetscObjectSetName((PetscObject)b2, "A_B2")); 213 if (size > 1 && !cong) { /* MATMPIAIJ ZeroRows and ZeroRowsColumns are buggy in this case */ 214 PetscCall(VecDestroy(&b1)); 215 } 216 if (ronl) { 217 PetscCall(MatZeroRowsIS(A, zr, diag, x, b1)); 218 PetscCall(MatZeroRowsIS(S, zr, diag, x, b2)); 219 } else { 220 PetscCall(MatZeroRowsColumnsIS(A, zr, diag, x, b1)); 221 PetscCall(MatZeroRowsColumnsIS(S, zr, diag, x, b2)); 222 PetscCall(ISDestroy(&zr)); 223 /* Mix zerorows and zerorowscols */ 224 nr = 0; 225 idxs = NULL; 226 if (rank == 0) { 227 nr = 1; 228 PetscCall(PetscMalloc1(nr, &idxs)); 229 if (test % 2) { 230 idxs[0] = (3 * M - 2 - test / 2) % M; 231 } else { 232 idxs[0] = (test / 2 + 1) % M; 233 } 234 idxs[0] = PetscMax(idxs[0], 0); 235 } 236 PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, nr, idxs, PETSC_OWN_POINTER, &zr)); 237 PetscCall(PetscObjectSetName((PetscObject)zr, "ZR2")); 238 PetscCall(ISViewFromOptions(zr, NULL, "-view_is")); 239 PetscCall(MatZeroRowsIS(A, zr, diag * 2.0 + PETSC_SMALL, NULL, NULL)); 240 PetscCall(MatZeroRowsIS(S, zr, diag * 2.0 + PETSC_SMALL, NULL, NULL)); 241 } 242 PetscCall(ISDestroy(&zr)); 243 244 if (b1) { 245 Vec b; 246 247 PetscCall(VecViewFromOptions(b1, NULL, "-view_b")); 248 PetscCall(VecViewFromOptions(b2, NULL, "-view_b")); 249 PetscCall(VecDuplicate(b1, &b)); 250 PetscCall(VecCopy(b1, b)); 251 PetscCall(VecAXPY(b, -1.0, b2)); 252 PetscCall(VecNorm(b, NORM_INFINITY, &err)); 253 if (err >= tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[test %" PetscInt_FMT "] Error b %g\n", test, (double)err)); 254 PetscCall(VecDestroy(&b)); 255 } 256 PetscCall(VecDestroy(&b1)); 257 PetscCall(VecDestroy(&b2)); 258 PetscCall(VecDestroy(&x)); 259 PetscCall(MatConvert(S, MATDENSE, MAT_INITIAL_MATRIX, &B)); 260 261 PetscCall(MatCreateTranspose(S, &ST)); 262 PetscCall(MatComputeOperator(ST, MATDENSE, &BT)); 263 PetscCall(MatTranspose(BT, MAT_INITIAL_MATRIX, &BTT)); 264 PetscCall(PetscObjectSetName((PetscObject)B, "S")); 265 PetscCall(PetscObjectSetName((PetscObject)BTT, "STT")); 266 PetscCall(MatConvert(A, MATDENSE, MAT_INITIAL_MATRIX, &C)); 267 PetscCall(PetscObjectSetName((PetscObject)C, "A")); 268 269 PetscCall(MatViewFromOptions(C, NULL, "-view_mat")); 270 PetscCall(MatViewFromOptions(B, NULL, "-view_mat")); 271 PetscCall(MatViewFromOptions(BTT, NULL, "-view_mat")); 272 273 PetscCall(MatAXPY(C, -1.0, B, SAME_NONZERO_PATTERN)); 274 PetscCall(MatNorm(C, NORM_FROBENIUS, &err)); 275 if (err >= tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[test %" PetscInt_FMT "] Error mat mult after %s %g\n", test, ronl ? "MatZeroRows" : "MatZeroRowsColumns", (double)err)); 276 277 PetscCall(MatConvert(A, MATDENSE, MAT_REUSE_MATRIX, &C)); 278 PetscCall(MatAXPY(C, -1.0, BTT, SAME_NONZERO_PATTERN)); 279 PetscCall(MatNorm(C, NORM_FROBENIUS, &err)); 280 if (err >= tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[test %" PetscInt_FMT "] Error mat mult transpose after %s %g\n", test, ronl ? "MatZeroRows" : "MatZeroRowsColumns", (double)err)); 281 282 PetscCall(MatDestroy(&ST)); 283 PetscCall(MatDestroy(&BTT)); 284 PetscCall(MatDestroy(&BT)); 285 PetscCall(MatDestroy(&B)); 286 PetscCall(MatDestroy(&C)); 287 } 288 if (testdiagscale) { /* MatDiagonalScale() */ 289 Vec vr, vl; 290 291 PetscCall(MatCreateVecs(A, &vr, &vl)); 292 if (randomize) { 293 PetscCall(VecSetRandom(vr, NULL)); 294 PetscCall(VecSetRandom(vl, NULL)); 295 } else { 296 PetscCall(VecSet(vr, test % 2 ? 0.15 : 1.0 / 0.15)); 297 PetscCall(VecSet(vl, test % 2 ? -1.2 : 1.0 / -1.2)); 298 } 299 PetscCall(MatDiagonalScale(A, vl, vr)); 300 PetscCall(MatDiagonalScale(S, vl, vr)); 301 PetscCall(VecDestroy(&vr)); 302 PetscCall(VecDestroy(&vl)); 303 } 304 305 if (testscale) { /* MatScale() */ 306 PetscCall(MatScale(A, test % 2 ? 1.4 : 1.0 / 1.4)); 307 PetscCall(MatScale(S, test % 2 ? 1.4 : 1.0 / 1.4)); 308 } 309 310 if (testshift && cong) { /* MatShift() : MATSHELL shift is broken when row/cols layout are not congruent and left/right scaling have been applied */ 311 PetscCall(MatShift(A, test % 2 ? -77.5 : 77.5)); 312 PetscCall(MatShift(S, test % 2 ? -77.5 : 77.5)); 313 } 314 315 if (testgetdiag && cong) { /* MatGetDiagonal() */ 316 Vec dA, dS; 317 318 PetscCall(MatCreateVecs(A, &dA, NULL)); 319 PetscCall(MatCreateVecs(S, &dS, NULL)); 320 PetscCall(MatGetDiagonal(A, dA)); 321 PetscCall(MatGetDiagonal(S, dS)); 322 PetscCall(VecAXPY(dA, -1.0, dS)); 323 PetscCall(VecNorm(dA, NORM_INFINITY, &err)); 324 if (err >= tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[test %" PetscInt_FMT "] Error diag %g\n", test, (double)err)); 325 PetscCall(VecDestroy(&dA)); 326 PetscCall(VecDestroy(&dS)); 327 } 328 329 if (testdup && !test) { 330 Mat A2, S2; 331 332 PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &A2)); 333 PetscCall(MatDuplicate(S, MAT_COPY_VALUES, &S2)); 334 PetscCall(MatDestroy(&A)); 335 PetscCall(MatDestroy(&S)); 336 A = A2; 337 S = S2; 338 } 339 340 if (testsubmat) { 341 Mat sA, sS, dA, dS, At, St; 342 IS r, c; 343 PetscInt rst, ren, cst, cen; 344 345 PetscCall(MatGetOwnershipRange(A, &rst, &ren)); 346 PetscCall(MatGetOwnershipRangeColumn(A, &cst, &cen)); 347 PetscCall(ISCreateStride(PetscObjectComm((PetscObject)A), (ren - rst) / 2, rst, 1, &r)); 348 PetscCall(ISCreateStride(PetscObjectComm((PetscObject)A), (cen - cst) / 2, cst, 1, &c)); 349 PetscCall(MatCreateSubMatrix(A, r, c, MAT_INITIAL_MATRIX, &sA)); 350 PetscCall(MatCreateSubMatrix(S, r, c, MAT_INITIAL_MATRIX, &sS)); 351 PetscCall(MatMultAddEqual(sA, sS, 10, &flg)); 352 if (!flg) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[test %" PetscInt_FMT "] Error submatrix mult add\n", test)); 353 PetscCall(MatMultTransposeAddEqual(sA, sS, 10, &flg)); 354 if (!flg) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[test %" PetscInt_FMT "] Error submatrix mult add (T)\n", test)); 355 PetscCall(MatConvert(sA, MATDENSE, MAT_INITIAL_MATRIX, &dA)); 356 PetscCall(MatConvert(sS, MATDENSE, MAT_INITIAL_MATRIX, &dS)); 357 PetscCall(MatAXPY(dA, -1.0, dS, SAME_NONZERO_PATTERN)); 358 PetscCall(MatNorm(dA, NORM_FROBENIUS, &err)); 359 if (err >= tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[test %" PetscInt_FMT "] Error mat submatrix %g\n", test, (double)err)); 360 PetscCall(MatDestroy(&sA)); 361 PetscCall(MatDestroy(&sS)); 362 PetscCall(MatDestroy(&dA)); 363 PetscCall(MatDestroy(&dS)); 364 PetscCall(MatCreateTranspose(A, &At)); 365 PetscCall(MatCreateTranspose(S, &St)); 366 PetscCall(MatCreateSubMatrix(At, c, r, MAT_INITIAL_MATRIX, &sA)); 367 PetscCall(MatCreateSubMatrix(St, c, r, MAT_INITIAL_MATRIX, &sS)); 368 PetscCall(MatMultAddEqual(sA, sS, 10, &flg)); 369 if (!flg) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[test %" PetscInt_FMT "] Error submatrix (T) mult add\n", test)); 370 PetscCall(MatMultTransposeAddEqual(sA, sS, 10, &flg)); 371 if (!flg) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[test %" PetscInt_FMT "] Error submatrix (T) mult add (T)\n", test)); 372 PetscCall(MatConvert(sA, MATDENSE, MAT_INITIAL_MATRIX, &dA)); 373 PetscCall(MatConvert(sS, MATDENSE, MAT_INITIAL_MATRIX, &dS)); 374 PetscCall(MatAXPY(dA, -1.0, dS, SAME_NONZERO_PATTERN)); 375 PetscCall(MatNorm(dA, NORM_FROBENIUS, &err)); 376 if (err >= tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[test %" PetscInt_FMT "] Error mat submatrix (T) %g\n", test, (double)err)); 377 PetscCall(MatDestroy(&sA)); 378 PetscCall(MatDestroy(&sS)); 379 PetscCall(MatDestroy(&dA)); 380 PetscCall(MatDestroy(&dS)); 381 PetscCall(MatDestroy(&At)); 382 PetscCall(MatDestroy(&St)); 383 PetscCall(ISDestroy(&r)); 384 PetscCall(ISDestroy(&c)); 385 } 386 387 if (testaxpy) { 388 Mat tA, tS, dA, dS; 389 MatStructure str[3] = {SAME_NONZERO_PATTERN, SUBSET_NONZERO_PATTERN, DIFFERENT_NONZERO_PATTERN}; 390 391 PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &tA)); 392 if (testaxpyd && !(test % 2)) { 393 PetscCall(PetscObjectReference((PetscObject)tA)); 394 tS = tA; 395 } else { 396 PetscCall(PetscObjectReference((PetscObject)S)); 397 tS = S; 398 } 399 PetscCall(MatAXPY(A, 0.5, tA, str[test % 3])); 400 PetscCall(MatAXPY(S, 0.5, tS, str[test % 3])); 401 /* this will trigger an error the next MatMult or MatMultTranspose call for S */ 402 if (testaxpyerr) PetscCall(MatScale(tA, 0)); 403 PetscCall(MatDestroy(&tA)); 404 PetscCall(MatDestroy(&tS)); 405 PetscCall(MatMultAddEqual(A, S, 10, &flg)); 406 if (!flg) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[test %" PetscInt_FMT "] Error axpy mult add\n", test)); 407 PetscCall(MatMultTransposeAddEqual(A, S, 10, &flg)); 408 if (!flg) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[test %" PetscInt_FMT "] Error axpy mult add (T)\n", test)); 409 PetscCall(MatConvert(A, MATDENSE, MAT_INITIAL_MATRIX, &dA)); 410 PetscCall(MatConvert(S, MATDENSE, MAT_INITIAL_MATRIX, &dS)); 411 PetscCall(MatAXPY(dA, -1.0, dS, SAME_NONZERO_PATTERN)); 412 PetscCall(MatNorm(dA, NORM_FROBENIUS, &err)); 413 if (err >= tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[test %" PetscInt_FMT "] Error mat submatrix %g\n", test, (double)err)); 414 PetscCall(MatDestroy(&dA)); 415 PetscCall(MatDestroy(&dS)); 416 } 417 418 if (testreset && (ntest == 1 || test == ntest - 2)) { 419 /* reset MATSHELL */ 420 PetscCall(MatAssemblyBegin(S, MAT_FINAL_ASSEMBLY)); 421 PetscCall(MatAssemblyEnd(S, MAT_FINAL_ASSEMBLY)); 422 /* reset A */ 423 PetscCall(MatCopy(user->B, A, DIFFERENT_NONZERO_PATTERN)); 424 } 425 } 426 427 PetscCall(MatDestroy(&A)); 428 PetscCall(MatDestroy(&S)); 429 PetscCall(PetscFree(user)); 430 PetscCall(PetscFinalize()); 431 return 0; 432 } 433 434 /*TEST 435 436 testset: 437 suffix: rect 438 requires: !single 439 output_file: output/empty.out 440 nsize: {{1 3}} 441 args: -loop 3 -keep {{0 1}} -M {{12 19}} -N {{19 12}} -submat {{0 1}} -test_axpy_different {{0 1}} 442 443 testset: 444 suffix: square 445 requires: !single 446 output_file: output/empty.out 447 nsize: {{1 3}} 448 args: -M 21 -N 21 -loop 4 -rows_only {{0 1}} -keep {{0 1}} -submat {{0 1}} -test_axpy_different {{0 1}} 449 TEST*/ 450