Actual source code: ex49.c

  1: static char help[] = "Tests MatTranspose(), MatNorm(), and MatAXPY().\n\n";

  3: #include <petscmat.h>

  5: int main(int argc, char **argv)
  6: {
  7:   Mat         mat, tmat = 0;
  8:   PetscInt    m = 4, n, i, j;
  9:   PetscMPIInt size, rank;
 10:   PetscInt    rstart, rend, rect = 0;
 11:   PetscBool   flg;
 12:   PetscScalar v;
 13:   PetscReal   normf, normi, norm1;
 14:   MatInfo     info;

 16:   PetscFunctionBeginUser;
 17:   PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
 18:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL));
 19:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
 20:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
 21:   n = m;
 22:   PetscCall(PetscOptionsHasName(NULL, NULL, "-rect1", &flg));
 23:   if (flg) {
 24:     n += 2;
 25:     rect = 1;
 26:   }
 27:   PetscCall(PetscOptionsHasName(NULL, NULL, "-rect2", &flg));
 28:   if (flg) {
 29:     n -= 2;
 30:     rect = 1;
 31:   }

 33:   /* Create and assemble matrix */
 34:   PetscCall(MatCreate(PETSC_COMM_WORLD, &mat));
 35:   PetscCall(MatSetSizes(mat, PETSC_DECIDE, PETSC_DECIDE, m, n));
 36:   PetscCall(MatSetFromOptions(mat));
 37:   PetscCall(MatSetUp(mat));
 38:   PetscCall(MatGetOwnershipRange(mat, &rstart, &rend));
 39:   for (i = rstart; i < rend; i++) {
 40:     for (j = 0; j < n; j++) {
 41:       v = 10 * i + j;
 42:       PetscCall(MatSetValues(mat, 1, &i, 1, &j, &v, INSERT_VALUES));
 43:     }
 44:   }
 45:   PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY));
 46:   PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY));

 48:   /* Print info about original matrix */
 49:   PetscCall(MatGetInfo(mat, MAT_GLOBAL_SUM, &info));
 50:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "original matrix nonzeros = %" PetscInt_FMT ", allocated nonzeros = %" PetscInt_FMT "\n", (PetscInt)info.nz_used, (PetscInt)info.nz_allocated));
 51:   PetscCall(MatNorm(mat, NORM_FROBENIUS, &normf));
 52:   PetscCall(MatNorm(mat, NORM_1, &norm1));
 53:   PetscCall(MatNorm(mat, NORM_INFINITY, &normi));
 54:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "original: Frobenious norm = %g, one norm = %g, infinity norm = %g\n", (double)normf, (double)norm1, (double)normi));
 55:   PetscCall(MatView(mat, PETSC_VIEWER_STDOUT_WORLD));

 57:   /* Form matrix transpose */
 58:   PetscCall(PetscOptionsHasName(NULL, NULL, "-in_place", &flg));
 59:   if (flg) {
 60:     PetscCall(MatTranspose(mat, MAT_INPLACE_MATRIX, &mat)); /* in-place transpose */
 61:     tmat = mat;
 62:     mat  = 0;
 63:   } else { /* out-of-place transpose */
 64:     PetscCall(MatTranspose(mat, MAT_INITIAL_MATRIX, &tmat));
 65:   }

 67:   /* Print info about transpose matrix */
 68:   PetscCall(MatGetInfo(tmat, MAT_GLOBAL_SUM, &info));
 69:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "transpose matrix nonzeros = %" PetscInt_FMT ", allocated nonzeros = %" PetscInt_FMT "\n", (PetscInt)info.nz_used, (PetscInt)info.nz_allocated));
 70:   PetscCall(MatNorm(tmat, NORM_FROBENIUS, &normf));
 71:   PetscCall(MatNorm(tmat, NORM_1, &norm1));
 72:   PetscCall(MatNorm(tmat, NORM_INFINITY, &normi));
 73:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "transpose: Frobenious norm = %g, one norm = %g, infinity norm = %g\n", (double)normf, (double)norm1, (double)normi));
 74:   PetscCall(MatView(tmat, PETSC_VIEWER_STDOUT_WORLD));

 76:   /* Test MatAXPY */
 77:   if (mat && !rect) {
 78:     PetscScalar alpha = 1.0;
 79:     PetscCall(PetscOptionsGetScalar(NULL, NULL, "-alpha", &alpha, NULL));
 80:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "matrix addition:  B = B + alpha * A\n"));
 81:     PetscCall(MatAXPY(tmat, alpha, mat, DIFFERENT_NONZERO_PATTERN));
 82:     PetscCall(MatView(tmat, PETSC_VIEWER_STDOUT_WORLD));
 83:   }

 85:   /* Free data structures */
 86:   PetscCall(MatDestroy(&tmat));
 87:   if (mat) PetscCall(MatDestroy(&mat));

 89:   PetscCall(PetscFinalize());
 90:   return 0;
 91: }

 93: /*TEST

 95:    test:

 97:    testset:
 98:      args: -rect1
 99:      test:
100:        suffix: r1
101:        output_file: output/ex49_r1.out
102:      test:
103:        suffix: r1_inplace
104:        args: -in_place
105:        output_file: output/ex49_r1.out
106:      test:
107:        suffix: r1_par
108:        nsize: 2
109:        output_file: output/ex49_r1_par.out
110:      test:
111:        suffix: r1_par_inplace
112:        args: -in_place
113:        nsize: 2
114:        output_file: output/ex49_r1_par.out

116: TEST*/