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, (char *)0, 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 for (i = 0; i < m; i++) data[j * m + i] = s2 * j + i + s1 + 1; 128 } 129 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 130 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 131 132 if (submat) { 133 Mat A2; 134 IS r, c; 135 PetscInt rst, ren, cst, cen; 136 137 PetscCall(MatGetOwnershipRange(A, &rst, &ren)); 138 PetscCall(MatGetOwnershipRangeColumn(A, &cst, &cen)); 139 PetscCall(ISCreateStride(PetscObjectComm((PetscObject)A), (ren - rst) / 2, rst, 1, &r)); 140 PetscCall(ISCreateStride(PetscObjectComm((PetscObject)A), (cen - cst) / 2, cst, 1, &c)); 141 PetscCall(MatCreateSubMatrix(A, r, c, MAT_INITIAL_MATRIX, &A2)); 142 PetscCall(ISDestroy(&r)); 143 PetscCall(ISDestroy(&c)); 144 PetscCall(MatDestroy(&A)); 145 A = A2; 146 } 147 148 PetscCall(MatGetSize(A, &M, &N)); 149 PetscCall(MatGetLocalSize(A, &m, &n)); 150 PetscCall(MatHasCongruentLayouts(A, &cong)); 151 152 PetscCall(MatConvert(A, MATAIJ, MAT_INPLACE_MATRIX, &A)); 153 PetscCall(MatSetOption(A, MAT_KEEP_NONZERO_PATTERN, keep)); 154 PetscCall(PetscObjectSetName((PetscObject)A, "initial")); 155 PetscCall(MatViewFromOptions(A, NULL, "-view_mat")); 156 157 PetscCall(PetscNew(&user)); 158 PetscCall(MatCreateShell(PETSC_COMM_WORLD, m, n, M, N, user, &S)); 159 PetscCall(MatShellSetOperation(S, MATOP_MULT, (void (*)(void))MatMult_User)); 160 PetscCall(MatShellSetOperation(S, MATOP_MULT_TRANSPOSE, (void (*)(void))MatMultTranspose_User)); 161 if (cong) PetscCall(MatShellSetOperation(S, MATOP_GET_DIAGONAL, (void (*)(void))MatGetDiagonal_User)); 162 PetscCall(MatShellSetOperation(S, MATOP_COPY, (void (*)(void))MatCopy_User)); 163 PetscCall(MatShellSetOperation(S, MATOP_DESTROY, (void (*)(void))MatDestroy_User)); 164 PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &user->B)); 165 166 /* Square and rows only scaling */ 167 ronl = cong ? ronl : PETSC_TRUE; 168 169 for (test = 0; test < ntest; test++) { 170 PetscReal err; 171 172 PetscCall(MatMultAddEqual(A, S, 10, &flg)); 173 if (!flg) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[test %" PetscInt_FMT "] Error mult add\n", test)); 174 PetscCall(MatMultTransposeAddEqual(A, S, 10, &flg)); 175 if (!flg) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[test %" PetscInt_FMT "] Error mult add (T)\n", test)); 176 if (testzerorows) { 177 Mat ST, B, C, BT, BTT; 178 IS zr; 179 Vec x = NULL, b1 = NULL, b2 = NULL; 180 PetscInt *idxs = NULL, nr = 0; 181 182 if (rank == (test % size)) { 183 nr = 1; 184 PetscCall(PetscMalloc1(nr, &idxs)); 185 if (test % 2) { 186 idxs[0] = (2 * M - 1 - test / 2) % M; 187 } else { 188 idxs[0] = (test / 2) % M; 189 } 190 idxs[0] = PetscMax(idxs[0], 0); 191 } 192 PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, nr, idxs, PETSC_OWN_POINTER, &zr)); 193 PetscCall(PetscObjectSetName((PetscObject)zr, "ZR")); 194 PetscCall(ISViewFromOptions(zr, NULL, "-view_is")); 195 PetscCall(MatCreateVecs(A, &x, &b1)); 196 if (randomize) { 197 PetscCall(VecSetRandom(x, NULL)); 198 PetscCall(VecSetRandom(b1, NULL)); 199 } else { 200 PetscCall(VecSet(x, 11.4)); 201 PetscCall(VecSet(b1, -14.2)); 202 } 203 PetscCall(VecDuplicate(b1, &b2)); 204 PetscCall(VecCopy(b1, b2)); 205 PetscCall(PetscObjectSetName((PetscObject)b1, "A_B1")); 206 PetscCall(PetscObjectSetName((PetscObject)b2, "A_B2")); 207 if (size > 1 && !cong) { /* MATMPIAIJ ZeroRows and ZeroRowsColumns are buggy in this case */ 208 PetscCall(VecDestroy(&b1)); 209 } 210 if (ronl) { 211 PetscCall(MatZeroRowsIS(A, zr, diag, x, b1)); 212 PetscCall(MatZeroRowsIS(S, zr, diag, x, b2)); 213 } else { 214 PetscCall(MatZeroRowsColumnsIS(A, zr, diag, x, b1)); 215 PetscCall(MatZeroRowsColumnsIS(S, zr, diag, x, b2)); 216 PetscCall(ISDestroy(&zr)); 217 /* Mix zerorows and zerorowscols */ 218 nr = 0; 219 idxs = NULL; 220 if (rank == 0) { 221 nr = 1; 222 PetscCall(PetscMalloc1(nr, &idxs)); 223 if (test % 2) { 224 idxs[0] = (3 * M - 2 - test / 2) % M; 225 } else { 226 idxs[0] = (test / 2 + 1) % M; 227 } 228 idxs[0] = PetscMax(idxs[0], 0); 229 } 230 PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, nr, idxs, PETSC_OWN_POINTER, &zr)); 231 PetscCall(PetscObjectSetName((PetscObject)zr, "ZR2")); 232 PetscCall(ISViewFromOptions(zr, NULL, "-view_is")); 233 PetscCall(MatZeroRowsIS(A, zr, diag * 2.0 + PETSC_SMALL, NULL, NULL)); 234 PetscCall(MatZeroRowsIS(S, zr, diag * 2.0 + PETSC_SMALL, NULL, NULL)); 235 } 236 PetscCall(ISDestroy(&zr)); 237 238 if (b1) { 239 Vec b; 240 241 PetscCall(VecViewFromOptions(b1, NULL, "-view_b")); 242 PetscCall(VecViewFromOptions(b2, NULL, "-view_b")); 243 PetscCall(VecDuplicate(b1, &b)); 244 PetscCall(VecCopy(b1, b)); 245 PetscCall(VecAXPY(b, -1.0, b2)); 246 PetscCall(VecNorm(b, NORM_INFINITY, &err)); 247 if (err >= tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[test %" PetscInt_FMT "] Error b %g\n", test, (double)err)); 248 PetscCall(VecDestroy(&b)); 249 } 250 PetscCall(VecDestroy(&b1)); 251 PetscCall(VecDestroy(&b2)); 252 PetscCall(VecDestroy(&x)); 253 PetscCall(MatConvert(S, MATDENSE, MAT_INITIAL_MATRIX, &B)); 254 255 PetscCall(MatCreateTranspose(S, &ST)); 256 PetscCall(MatComputeOperator(ST, MATDENSE, &BT)); 257 PetscCall(MatTranspose(BT, MAT_INITIAL_MATRIX, &BTT)); 258 PetscCall(PetscObjectSetName((PetscObject)B, "S")); 259 PetscCall(PetscObjectSetName((PetscObject)BTT, "STT")); 260 PetscCall(MatConvert(A, MATDENSE, MAT_INITIAL_MATRIX, &C)); 261 PetscCall(PetscObjectSetName((PetscObject)C, "A")); 262 263 PetscCall(MatViewFromOptions(C, NULL, "-view_mat")); 264 PetscCall(MatViewFromOptions(B, NULL, "-view_mat")); 265 PetscCall(MatViewFromOptions(BTT, NULL, "-view_mat")); 266 267 PetscCall(MatAXPY(C, -1.0, B, SAME_NONZERO_PATTERN)); 268 PetscCall(MatNorm(C, NORM_FROBENIUS, &err)); 269 if (err >= tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[test %" PetscInt_FMT "] Error mat mult after %s %g\n", test, ronl ? "MatZeroRows" : "MatZeroRowsColumns", (double)err)); 270 271 PetscCall(MatConvert(A, MATDENSE, MAT_REUSE_MATRIX, &C)); 272 PetscCall(MatAXPY(C, -1.0, BTT, SAME_NONZERO_PATTERN)); 273 PetscCall(MatNorm(C, NORM_FROBENIUS, &err)); 274 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)); 275 276 PetscCall(MatDestroy(&ST)); 277 PetscCall(MatDestroy(&BTT)); 278 PetscCall(MatDestroy(&BT)); 279 PetscCall(MatDestroy(&B)); 280 PetscCall(MatDestroy(&C)); 281 } 282 if (testdiagscale) { /* MatDiagonalScale() */ 283 Vec vr, vl; 284 285 PetscCall(MatCreateVecs(A, &vr, &vl)); 286 if (randomize) { 287 PetscCall(VecSetRandom(vr, NULL)); 288 PetscCall(VecSetRandom(vl, NULL)); 289 } else { 290 PetscCall(VecSet(vr, test % 2 ? 0.15 : 1.0 / 0.15)); 291 PetscCall(VecSet(vl, test % 2 ? -1.2 : 1.0 / -1.2)); 292 } 293 PetscCall(MatDiagonalScale(A, vl, vr)); 294 PetscCall(MatDiagonalScale(S, vl, vr)); 295 PetscCall(VecDestroy(&vr)); 296 PetscCall(VecDestroy(&vl)); 297 } 298 299 if (testscale) { /* MatScale() */ 300 PetscCall(MatScale(A, test % 2 ? 1.4 : 1.0 / 1.4)); 301 PetscCall(MatScale(S, test % 2 ? 1.4 : 1.0 / 1.4)); 302 } 303 304 if (testshift && cong) { /* MatShift() : MATSHELL shift is broken when row/cols layout are not congruent and left/right scaling have been applied */ 305 PetscCall(MatShift(A, test % 2 ? -77.5 : 77.5)); 306 PetscCall(MatShift(S, test % 2 ? -77.5 : 77.5)); 307 } 308 309 if (testgetdiag && cong) { /* MatGetDiagonal() */ 310 Vec dA, dS; 311 312 PetscCall(MatCreateVecs(A, &dA, NULL)); 313 PetscCall(MatCreateVecs(S, &dS, NULL)); 314 PetscCall(MatGetDiagonal(A, dA)); 315 PetscCall(MatGetDiagonal(S, dS)); 316 PetscCall(VecAXPY(dA, -1.0, dS)); 317 PetscCall(VecNorm(dA, NORM_INFINITY, &err)); 318 if (err >= tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[test %" PetscInt_FMT "] Error diag %g\n", test, (double)err)); 319 PetscCall(VecDestroy(&dA)); 320 PetscCall(VecDestroy(&dS)); 321 } 322 323 if (testdup && !test) { 324 Mat A2, S2; 325 326 PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &A2)); 327 PetscCall(MatDuplicate(S, MAT_COPY_VALUES, &S2)); 328 PetscCall(MatDestroy(&A)); 329 PetscCall(MatDestroy(&S)); 330 A = A2; 331 S = S2; 332 } 333 334 if (testsubmat) { 335 Mat sA, sS, dA, dS, At, St; 336 IS r, c; 337 PetscInt rst, ren, cst, cen; 338 339 PetscCall(MatGetOwnershipRange(A, &rst, &ren)); 340 PetscCall(MatGetOwnershipRangeColumn(A, &cst, &cen)); 341 PetscCall(ISCreateStride(PetscObjectComm((PetscObject)A), (ren - rst) / 2, rst, 1, &r)); 342 PetscCall(ISCreateStride(PetscObjectComm((PetscObject)A), (cen - cst) / 2, cst, 1, &c)); 343 PetscCall(MatCreateSubMatrix(A, r, c, MAT_INITIAL_MATRIX, &sA)); 344 PetscCall(MatCreateSubMatrix(S, r, c, MAT_INITIAL_MATRIX, &sS)); 345 PetscCall(MatMultAddEqual(sA, sS, 10, &flg)); 346 if (!flg) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[test %" PetscInt_FMT "] Error submatrix mult add\n", test)); 347 PetscCall(MatMultTransposeAddEqual(sA, sS, 10, &flg)); 348 if (!flg) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[test %" PetscInt_FMT "] Error submatrix mult add (T)\n", test)); 349 PetscCall(MatConvert(sA, MATDENSE, MAT_INITIAL_MATRIX, &dA)); 350 PetscCall(MatConvert(sS, MATDENSE, MAT_INITIAL_MATRIX, &dS)); 351 PetscCall(MatAXPY(dA, -1.0, dS, SAME_NONZERO_PATTERN)); 352 PetscCall(MatNorm(dA, NORM_FROBENIUS, &err)); 353 if (err >= tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[test %" PetscInt_FMT "] Error mat submatrix %g\n", test, (double)err)); 354 PetscCall(MatDestroy(&sA)); 355 PetscCall(MatDestroy(&sS)); 356 PetscCall(MatDestroy(&dA)); 357 PetscCall(MatDestroy(&dS)); 358 PetscCall(MatCreateTranspose(A, &At)); 359 PetscCall(MatCreateTranspose(S, &St)); 360 PetscCall(MatCreateSubMatrix(At, c, r, MAT_INITIAL_MATRIX, &sA)); 361 PetscCall(MatCreateSubMatrix(St, c, r, MAT_INITIAL_MATRIX, &sS)); 362 PetscCall(MatMultAddEqual(sA, sS, 10, &flg)); 363 if (!flg) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[test %" PetscInt_FMT "] Error submatrix (T) mult add\n", test)); 364 PetscCall(MatMultTransposeAddEqual(sA, sS, 10, &flg)); 365 if (!flg) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[test %" PetscInt_FMT "] Error submatrix (T) mult add (T)\n", test)); 366 PetscCall(MatConvert(sA, MATDENSE, MAT_INITIAL_MATRIX, &dA)); 367 PetscCall(MatConvert(sS, MATDENSE, MAT_INITIAL_MATRIX, &dS)); 368 PetscCall(MatAXPY(dA, -1.0, dS, SAME_NONZERO_PATTERN)); 369 PetscCall(MatNorm(dA, NORM_FROBENIUS, &err)); 370 if (err >= tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[test %" PetscInt_FMT "] Error mat submatrix (T) %g\n", test, (double)err)); 371 PetscCall(MatDestroy(&sA)); 372 PetscCall(MatDestroy(&sS)); 373 PetscCall(MatDestroy(&dA)); 374 PetscCall(MatDestroy(&dS)); 375 PetscCall(MatDestroy(&At)); 376 PetscCall(MatDestroy(&St)); 377 PetscCall(ISDestroy(&r)); 378 PetscCall(ISDestroy(&c)); 379 } 380 381 if (testaxpy) { 382 Mat tA, tS, dA, dS; 383 MatStructure str[3] = {SAME_NONZERO_PATTERN, SUBSET_NONZERO_PATTERN, DIFFERENT_NONZERO_PATTERN}; 384 385 PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &tA)); 386 if (testaxpyd && !(test % 2)) { 387 PetscCall(PetscObjectReference((PetscObject)tA)); 388 tS = tA; 389 } else { 390 PetscCall(PetscObjectReference((PetscObject)S)); 391 tS = S; 392 } 393 PetscCall(MatAXPY(A, 0.5, tA, str[test % 3])); 394 PetscCall(MatAXPY(S, 0.5, tS, str[test % 3])); 395 /* this will trigger an error the next MatMult or MatMultTranspose call for S */ 396 if (testaxpyerr) PetscCall(MatScale(tA, 0)); 397 PetscCall(MatDestroy(&tA)); 398 PetscCall(MatDestroy(&tS)); 399 PetscCall(MatMultAddEqual(A, S, 10, &flg)); 400 if (!flg) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[test %" PetscInt_FMT "] Error axpy mult add\n", test)); 401 PetscCall(MatMultTransposeAddEqual(A, S, 10, &flg)); 402 if (!flg) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[test %" PetscInt_FMT "] Error axpy mult add (T)\n", test)); 403 PetscCall(MatConvert(A, MATDENSE, MAT_INITIAL_MATRIX, &dA)); 404 PetscCall(MatConvert(S, MATDENSE, MAT_INITIAL_MATRIX, &dS)); 405 PetscCall(MatAXPY(dA, -1.0, dS, SAME_NONZERO_PATTERN)); 406 PetscCall(MatNorm(dA, NORM_FROBENIUS, &err)); 407 if (err >= tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[test %" PetscInt_FMT "] Error mat submatrix %g\n", test, (double)err)); 408 PetscCall(MatDestroy(&dA)); 409 PetscCall(MatDestroy(&dS)); 410 } 411 412 if (testreset && (ntest == 1 || test == ntest - 2)) { 413 /* reset MATSHELL */ 414 PetscCall(MatAssemblyBegin(S, MAT_FINAL_ASSEMBLY)); 415 PetscCall(MatAssemblyEnd(S, MAT_FINAL_ASSEMBLY)); 416 /* reset A */ 417 PetscCall(MatCopy(user->B, A, DIFFERENT_NONZERO_PATTERN)); 418 } 419 } 420 421 PetscCall(MatDestroy(&A)); 422 PetscCall(MatDestroy(&S)); 423 PetscCall(PetscFree(user)); 424 PetscCall(PetscFinalize()); 425 return 0; 426 } 427 428 /*TEST 429 430 testset: 431 suffix: rect 432 requires: !single 433 output_file: output/ex221_1.out 434 nsize: {{1 3}} 435 args: -loop 3 -keep {{0 1}} -M {{12 19}} -N {{19 12}} -submat {{0 1}} -test_axpy_different {{0 1}} 436 437 testset: 438 suffix: square 439 requires: !single 440 output_file: output/ex221_1.out 441 nsize: {{1 3}} 442 args: -M 21 -N 21 -loop 4 -rows_only {{0 1}} -keep {{0 1}} -submat {{0 1}} -test_axpy_different {{0 1}} 443 TEST*/ 444