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 { 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; 16 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 } 33 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)); 48 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)); 57 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 } 67 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)); 76 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 } 85 86 /* Free data structures */ 87 PetscCall(MatDestroy(&tmat)); 88 if (mat) PetscCall(MatDestroy(&mat)); 89 90 PetscCall(PetscFinalize()); 91 return 0; 92 } 93 94 /*TEST 95 96 test: 97 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 116 117 TEST*/ 118