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