1 static char help[] = "Tests MatComputeOperator() and MatComputeOperatorTranspose()\n\n"; 2 3 #include <petscmat.h> 4 5 int main(int argc, char **argv) 6 { 7 Mat A, Ae, Aet; 8 char filename[PETSC_MAX_PATH_LEN]; 9 char expltype[128], *etype = NULL; 10 PetscInt bs = 1; 11 PetscBool flg, check = PETSC_TRUE; 12 13 PetscFunctionBeginUser; 14 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 15 16 PetscCall(PetscOptionsGetString(NULL, NULL, "-expl_type", expltype, sizeof(expltype), &flg)); 17 if (flg) PetscCall(PetscStrallocpy(expltype, &etype)); 18 PetscCall(PetscOptionsGetString(NULL, NULL, "-f", filename, sizeof(filename), &flg)); 19 PetscCall(PetscOptionsGetInt(NULL, NULL, "-bs", &bs, NULL)); 20 if (!flg) { 21 PetscInt M = 13, N = 6; 22 PetscBool sbaij; 23 24 PetscCall(PetscOptionsGetInt(NULL, NULL, "-M", &M, NULL)); 25 PetscCall(PetscOptionsGetInt(NULL, NULL, "-N", &N, NULL)); 26 PetscCall(MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, M, N, NULL, &A)); 27 PetscCall(MatSetBlockSize(A, bs)); 28 PetscCall(MatSetRandom(A, NULL)); 29 PetscCall(PetscStrcmp(etype, MATSBAIJ, &sbaij)); 30 if (sbaij) { 31 Mat At; 32 PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &At)); 33 PetscCall(MatAXPY(A, 1.0, At, DIFFERENT_NONZERO_PATTERN)); 34 PetscCall(MatDestroy(&At)); 35 } 36 } else { 37 PetscViewer viewer; 38 39 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, filename, FILE_MODE_READ, &viewer)); 40 PetscCall(MatCreate(PETSC_COMM_WORLD, &A)); 41 PetscCall(MatSetBlockSize(A, bs)); 42 PetscCall(MatSetFromOptions(A)); 43 PetscCall(MatLoad(A, viewer)); 44 PetscCall(PetscViewerDestroy(&viewer)); 45 } 46 PetscCall(PetscObjectSetName((PetscObject)A, "Matrix")); 47 PetscCall(MatViewFromOptions(A, NULL, "-view_expl")); 48 49 PetscCall(MatComputeOperator(A, etype, &Ae)); 50 PetscCall(PetscObjectSetName((PetscObject)Ae, "Explicit matrix")); 51 PetscCall(MatViewFromOptions(Ae, NULL, "-view_expl")); 52 53 PetscCall(PetscOptionsGetBool(NULL, NULL, "-check", &check, NULL)); 54 if (check) { 55 Mat A2; 56 PetscReal err, tol = PETSC_SMALL; 57 58 PetscCall(PetscOptionsGetReal(NULL, NULL, "-tol", &tol, NULL)); 59 PetscCall(MatConvert(A, etype, MAT_INITIAL_MATRIX, &A2)); 60 PetscCall(MatAXPY(A2, -1.0, Ae, DIFFERENT_NONZERO_PATTERN)); 61 PetscCall(MatNorm(A2, NORM_FROBENIUS, &err)); 62 if (err > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error %g > %g (type %s)\n", (double)err, (double)tol, etype)); 63 PetscCall(MatDestroy(&A2)); 64 } 65 66 PetscCall(MatComputeOperatorTranspose(A, etype, &Aet)); 67 PetscCall(PetscObjectSetName((PetscObject)Aet, "Explicit matrix transpose")); 68 PetscCall(MatViewFromOptions(Aet, NULL, "-view_expl")); 69 70 PetscCall(PetscFree(etype)); 71 PetscCall(MatDestroy(&Ae)); 72 PetscCall(MatDestroy(&Aet)); 73 PetscCall(MatDestroy(&A)); 74 PetscCall(PetscFinalize()); 75 return 0; 76 } 77 78 /*TEST 79 80 test: 81 output_file: output/empty.out 82 83 testset: 84 suffix: matexpl_rect 85 output_file: output/empty.out 86 nsize: {{1 3}} 87 args: -expl_type {{dense aij baij}} 88 89 testset: 90 suffix: matexpl_square 91 output_file: output/empty.out 92 nsize: {{1 3}} 93 args: -bs {{1 2 3}} -M 36 -N 36 -expl_type {{dense aij baij sbaij}} 94 95 TEST*/ 96