Actual source code: ex222.c

  1: static char help[] = "Tests MatComputeOperator() and MatComputeOperatorTranspose()\n\n";

  3: #include <petscmat.h>

  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;

 13:   PetscFunctionBeginUser;
 14:   PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));

 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;

 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;

 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"));

 49:   PetscCall(MatComputeOperator(A, etype, &Ae));
 50:   PetscCall(PetscObjectSetName((PetscObject)Ae, "Explicit matrix"));
 51:   PetscCall(MatViewFromOptions(Ae, NULL, "-view_expl"));

 53:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-check", &check, NULL));
 54:   if (check) {
 55:     Mat       A2;
 56:     PetscReal err, tol = PETSC_SMALL;

 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:   }

 66:   PetscCall(MatComputeOperatorTranspose(A, etype, &Aet));
 67:   PetscCall(PetscObjectSetName((PetscObject)Aet, "Explicit matrix transpose"));
 68:   PetscCall(MatViewFromOptions(Aet, NULL, "-view_expl"));

 70:   PetscCall(PetscFree(etype));
 71:   PetscCall(MatDestroy(&Ae));
 72:   PetscCall(MatDestroy(&Aet));
 73:   PetscCall(MatDestroy(&A));
 74:   PetscCall(PetscFinalize());
 75:   return 0;
 76: }

 78: /*TEST

 80:    test:
 81:      output_file: output/ex222_null.out

 83:    testset:
 84:      suffix: matexpl_rect
 85:      output_file: output/ex222_null.out
 86:      nsize: {{1 3}}
 87:      args: -expl_type {{dense aij baij}}

 89:    testset:
 90:      suffix: matexpl_square
 91:      output_file: output/ex222_null.out
 92:      nsize: {{1 3}}
 93:      args: -bs {{1 2 3}} -M 36 -N 36 -expl_type {{dense aij baij sbaij}}

 95: TEST*/