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