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, (char *)0, 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 23 PetscCall(PetscOptionsGetInt(NULL, NULL, "-M", &M, NULL)); 24 PetscCall(PetscOptionsGetInt(NULL, NULL, "-N", &N, NULL)); 25 PetscCall(MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, M, N, NULL, &A)); 26 PetscCall(MatSetBlockSize(A, bs)); 27 PetscCall(MatSetRandom(A, NULL)); 28 } else { 29 PetscViewer viewer; 30 31 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, filename, FILE_MODE_READ, &viewer)); 32 PetscCall(MatCreate(PETSC_COMM_WORLD, &A)); 33 PetscCall(MatSetBlockSize(A, bs)); 34 PetscCall(MatSetFromOptions(A)); 35 PetscCall(MatLoad(A, viewer)); 36 PetscCall(PetscViewerDestroy(&viewer)); 37 } 38 PetscCall(PetscObjectSetName((PetscObject)A, "Matrix")); 39 PetscCall(MatViewFromOptions(A, NULL, "-view_expl")); 40 41 PetscCall(MatComputeOperator(A, etype, &Ae)); 42 PetscCall(PetscObjectSetName((PetscObject)Ae, "Explicit matrix")); 43 PetscCall(MatViewFromOptions(Ae, NULL, "-view_expl")); 44 45 PetscCall(PetscOptionsGetBool(NULL, NULL, "-check", &check, NULL)); 46 if (check) { 47 Mat A2; 48 PetscReal err, tol = PETSC_SMALL; 49 50 PetscCall(PetscOptionsGetReal(NULL, NULL, "-tol", &tol, NULL)); 51 PetscCall(MatConvert(A, etype, MAT_INITIAL_MATRIX, &A2)); 52 PetscCall(MatAXPY(A2, -1.0, Ae, DIFFERENT_NONZERO_PATTERN)); 53 PetscCall(MatNorm(A2, NORM_FROBENIUS, &err)); 54 if (err > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error %g > %g (type %s)\n", (double)err, (double)tol, etype)); 55 PetscCall(MatDestroy(&A2)); 56 } 57 58 PetscCall(MatComputeOperatorTranspose(A, etype, &Aet)); 59 PetscCall(PetscObjectSetName((PetscObject)Aet, "Explicit matrix transpose")); 60 PetscCall(MatViewFromOptions(Aet, NULL, "-view_expl")); 61 62 PetscCall(PetscFree(etype)); 63 PetscCall(MatDestroy(&Ae)); 64 PetscCall(MatDestroy(&Aet)); 65 PetscCall(MatDestroy(&A)); 66 PetscCall(PetscFinalize()); 67 return 0; 68 } 69 70 /*TEST 71 72 test: 73 output_file: output/ex222_null.out 74 75 testset: 76 suffix: matexpl_rect 77 output_file: output/ex222_null.out 78 nsize: {{1 3}} 79 args: -expl_type {{dense aij baij}} 80 81 testset: 82 suffix: matexpl_square 83 output_file: output/ex222_null.out 84 nsize: {{1 3}} 85 args: -bs {{1 2 3}} -M 36 -N 36 -expl_type {{dense aij baij sbaij}} 86 87 TEST*/ 88