1 2 static char help[] = "Tests reusing MPI parallel matrices and MatGetValues().\n\ 3 To test the parallel matrix assembly, this example intentionally lays out\n\ 4 the matrix across processors differently from the way it is assembled.\n\ 5 This example uses bilinear elements on the unit square. Input arguments are:\n\ 6 -m <size> : problem size\n\n"; 7 8 #include <petscmat.h> 9 10 PetscErrorCode FormElementStiffness(PetscReal H, PetscScalar *Ke) 11 { 12 PetscFunctionBegin; 13 Ke[0] = H / 6.0; 14 Ke[1] = -.125 * H; 15 Ke[2] = H / 12.0; 16 Ke[3] = -.125 * H; 17 Ke[4] = -.125 * H; 18 Ke[5] = H / 6.0; 19 Ke[6] = -.125 * H; 20 Ke[7] = H / 12.0; 21 Ke[8] = H / 12.0; 22 Ke[9] = -.125 * H; 23 Ke[10] = H / 6.0; 24 Ke[11] = -.125 * H; 25 Ke[12] = -.125 * H; 26 Ke[13] = H / 12.0; 27 Ke[14] = -.125 * H; 28 Ke[15] = H / 6.0; 29 PetscFunctionReturn(PETSC_SUCCESS); 30 } 31 32 int main(int argc, char **args) 33 { 34 Mat C; 35 Vec u, b; 36 PetscMPIInt size, rank; 37 PetscInt i, m = 5, N, start, end, M, idx[4]; 38 PetscInt j, nrsub, ncsub, *rsub, *csub, mystart, myend; 39 PetscBool flg; 40 PetscScalar one = 1.0, Ke[16], *vals; 41 PetscReal h, norm; 42 43 PetscFunctionBeginUser; 44 PetscCall(PetscInitialize(&argc, &args, (char *)0, help)); 45 PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL)); 46 47 N = (m + 1) * (m + 1); /* dimension of matrix */ 48 M = m * m; /* number of elements */ 49 h = 1.0 / m; /* mesh width */ 50 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank)); 51 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 52 53 /* Create stiffness matrix */ 54 PetscCall(MatCreate(PETSC_COMM_WORLD, &C)); 55 PetscCall(MatSetSizes(C, PETSC_DECIDE, PETSC_DECIDE, N, N)); 56 PetscCall(MatSetFromOptions(C)); 57 PetscCall(MatSetUp(C)); 58 59 start = rank * (M / size) + ((M % size) < rank ? (M % size) : rank); 60 end = start + M / size + ((M % size) > rank); 61 62 /* Form the element stiffness for the Laplacian */ 63 PetscCall(FormElementStiffness(h * h, Ke)); 64 for (i = start; i < end; i++) { 65 /* location of lower left corner of element */ 66 /* node numbers for the four corners of element */ 67 idx[0] = (m + 1) * (i / m) + (i % m); 68 idx[1] = idx[0] + 1; 69 idx[2] = idx[1] + m + 1; 70 idx[3] = idx[2] - 1; 71 PetscCall(MatSetValues(C, 4, idx, 4, idx, Ke, ADD_VALUES)); 72 } 73 PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY)); 74 PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY)); 75 76 /* Assemble the matrix again */ 77 PetscCall(MatZeroEntries(C)); 78 79 for (i = start; i < end; i++) { 80 /* location of lower left corner of element */ 81 /* node numbers for the four corners of element */ 82 idx[0] = (m + 1) * (i / m) + (i % m); 83 idx[1] = idx[0] + 1; 84 idx[2] = idx[1] + m + 1; 85 idx[3] = idx[2] - 1; 86 PetscCall(MatSetValues(C, 4, idx, 4, idx, Ke, ADD_VALUES)); 87 } 88 PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY)); 89 PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY)); 90 91 /* Create test vectors */ 92 PetscCall(VecCreate(PETSC_COMM_WORLD, &u)); 93 PetscCall(VecSetSizes(u, PETSC_DECIDE, N)); 94 PetscCall(VecSetFromOptions(u)); 95 PetscCall(VecDuplicate(u, &b)); 96 PetscCall(VecSet(u, one)); 97 98 /* Check error */ 99 PetscCall(MatMult(C, u, b)); 100 PetscCall(VecNorm(b, NORM_2, &norm)); 101 if (norm > PETSC_SQRT_MACHINE_EPSILON) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Norm of error b %g should be near 0\n", (double)norm)); 102 103 /* Now test MatGetValues() */ 104 PetscCall(PetscOptionsHasName(NULL, NULL, "-get_values", &flg)); 105 if (flg) { 106 PetscCall(MatGetOwnershipRange(C, &mystart, &myend)); 107 nrsub = myend - mystart; 108 ncsub = 4; 109 PetscCall(PetscMalloc1(nrsub * ncsub, &vals)); 110 PetscCall(PetscMalloc1(nrsub, &rsub)); 111 PetscCall(PetscMalloc1(ncsub, &csub)); 112 for (i = myend - 1; i >= mystart; i--) rsub[myend - i - 1] = i; 113 for (i = 0; i < ncsub; i++) csub[i] = 2 * (ncsub - i) + mystart; 114 PetscCall(MatGetValues(C, nrsub, rsub, ncsub, csub, vals)); 115 PetscCall(MatView(C, PETSC_VIEWER_STDOUT_WORLD)); 116 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "processor number %d: start=%" PetscInt_FMT ", end=%" PetscInt_FMT ", mystart=%" PetscInt_FMT ", myend=%" PetscInt_FMT "\n", rank, start, end, mystart, myend)); 117 for (i = 0; i < nrsub; i++) { 118 for (j = 0; j < ncsub; j++) { 119 if (PetscImaginaryPart(vals[i * ncsub + j]) != 0.0) { 120 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, " C[%" PetscInt_FMT ", %" PetscInt_FMT "] = %g + %g i\n", rsub[i], csub[j], (double)PetscRealPart(vals[i * ncsub + j]), (double)PetscImaginaryPart(vals[i * ncsub + j]))); 121 } else { 122 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, " C[%" PetscInt_FMT ", %" PetscInt_FMT "] = %g\n", rsub[i], csub[j], (double)PetscRealPart(vals[i * ncsub + j]))); 123 } 124 } 125 } 126 PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT)); 127 PetscCall(PetscFree(rsub)); 128 PetscCall(PetscFree(csub)); 129 PetscCall(PetscFree(vals)); 130 } 131 132 /* Free data structures */ 133 PetscCall(VecDestroy(&u)); 134 PetscCall(VecDestroy(&b)); 135 PetscCall(MatDestroy(&C)); 136 PetscCall(PetscFinalize()); 137 return 0; 138 } 139 140 /*TEST 141 142 test: 143 nsize: 4 144 145 TEST*/ 146