Actual source code: ex49.c


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

  4: #include <petscmat.h>

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

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

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

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

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

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

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

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

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

 94: /*TEST

 96:    test:

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

117: TEST*/