1 2 static char help[] = "Tests MatTranspose(), MatNorm(), and MatAXPY().\n\n"; 3 4 #include <petscmat.h> 5 6 int main(int argc, char **argv) { 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; 15 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 } 32 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)); 47 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)); 56 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 } 66 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)); 75 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 } 84 85 /* Free data structures */ 86 PetscCall(MatDestroy(&tmat)); 87 if (mat) PetscCall(MatDestroy(&mat)); 88 89 PetscCall(PetscFinalize()); 90 return 0; 91 } 92 93 /*TEST 94 95 test: 96 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 115 116 TEST*/ 117