1 2 static char help[] = "Tests matrix factorization. Note that most users should\n\ 3 employ the KSP interface to the linear solvers instead of using the factorization\n\ 4 routines directly.\n\n"; 5 6 #include <petscmat.h> 7 8 int main(int argc, char **args) 9 { 10 Mat C, LU; 11 MatInfo info; 12 PetscInt i, j, m, n, Ii, J; 13 PetscScalar v, one = 1.0; 14 IS perm, iperm; 15 Vec x, u, b; 16 PetscReal norm, fill; 17 MatFactorInfo luinfo; 18 19 PetscFunctionBeginUser; 20 PetscCall(PetscInitialize(&argc, &args, (char *)0, help)); 21 22 PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Mat test ex7 options", "Mat"); 23 m = 3; 24 n = 3; 25 fill = 2.0; 26 PetscCall(PetscOptionsInt("-m", "Number of rows in grid", NULL, m, &m, NULL)); 27 PetscCall(PetscOptionsInt("-n", "Number of columns in grid", NULL, n, &n, NULL)); 28 PetscCall(PetscOptionsReal("-fill", "Expected fill ratio for factorization", NULL, fill, &fill, NULL)); 29 30 PetscOptionsEnd(); 31 32 /* Create the matrix for the five point stencil, YET AGAIN */ 33 PetscCall(MatCreate(PETSC_COMM_SELF, &C)); 34 PetscCall(MatSetSizes(C, PETSC_DECIDE, PETSC_DECIDE, m * n, m * n)); 35 PetscCall(MatSetFromOptions(C)); 36 PetscCall(MatSetUp(C)); 37 for (i = 0; i < m; i++) { 38 for (j = 0; j < n; j++) { 39 v = -1.0; 40 Ii = j + n * i; 41 if (i > 0) { 42 J = Ii - n; 43 PetscCall(MatSetValues(C, 1, &Ii, 1, &J, &v, INSERT_VALUES)); 44 } 45 if (i < m - 1) { 46 J = Ii + n; 47 PetscCall(MatSetValues(C, 1, &Ii, 1, &J, &v, INSERT_VALUES)); 48 } 49 if (j > 0) { 50 J = Ii - 1; 51 PetscCall(MatSetValues(C, 1, &Ii, 1, &J, &v, INSERT_VALUES)); 52 } 53 if (j < n - 1) { 54 J = Ii + 1; 55 PetscCall(MatSetValues(C, 1, &Ii, 1, &J, &v, INSERT_VALUES)); 56 } 57 v = 4.0; 58 PetscCall(MatSetValues(C, 1, &Ii, 1, &Ii, &v, INSERT_VALUES)); 59 } 60 } 61 PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY)); 62 PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY)); 63 PetscCall(MatGetOrdering(C, MATORDERINGRCM, &perm, &iperm)); 64 PetscCall(MatView(C, PETSC_VIEWER_STDOUT_WORLD)); 65 PetscCall(ISView(perm, PETSC_VIEWER_STDOUT_SELF)); 66 67 PetscCall(MatFactorInfoInitialize(&luinfo)); 68 69 luinfo.fill = fill; 70 luinfo.dtcol = 0.0; 71 luinfo.zeropivot = 1.e-14; 72 luinfo.pivotinblocks = 1.0; 73 74 PetscCall(MatGetFactor(C, MATSOLVERPETSC, MAT_FACTOR_LU, &LU)); 75 PetscCall(MatLUFactorSymbolic(LU, C, perm, iperm, &luinfo)); 76 PetscCall(MatLUFactorNumeric(LU, C, &luinfo)); 77 78 PetscCall(VecCreateSeq(PETSC_COMM_SELF, m * n, &u)); 79 PetscCall(VecSet(u, one)); 80 PetscCall(VecDuplicate(u, &x)); 81 PetscCall(VecDuplicate(u, &b)); 82 83 PetscCall(MatMult(C, u, b)); 84 PetscCall(MatSolve(LU, b, x)); 85 86 PetscCall(VecView(b, PETSC_VIEWER_STDOUT_SELF)); 87 PetscCall(VecView(x, PETSC_VIEWER_STDOUT_SELF)); 88 89 PetscCall(VecAXPY(x, -1.0, u)); 90 PetscCall(VecNorm(x, NORM_2, &norm)); 91 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Norm of error %g\n", (double)norm)); 92 93 PetscCall(MatGetInfo(C, MAT_LOCAL, &info)); 94 PetscCall(PetscPrintf(PETSC_COMM_SELF, "original matrix nonzeros = %" PetscInt_FMT "\n", (PetscInt)info.nz_used)); 95 PetscCall(MatGetInfo(LU, MAT_LOCAL, &info)); 96 PetscCall(PetscPrintf(PETSC_COMM_SELF, "factored matrix nonzeros = %" PetscInt_FMT "\n", (PetscInt)info.nz_used)); 97 98 PetscCall(VecDestroy(&u)); 99 PetscCall(VecDestroy(&b)); 100 PetscCall(VecDestroy(&x)); 101 PetscCall(ISDestroy(&perm)); 102 PetscCall(ISDestroy(&iperm)); 103 PetscCall(MatDestroy(&C)); 104 PetscCall(MatDestroy(&LU)); 105 106 PetscCall(PetscFinalize()); 107 return 0; 108 } 109 110 /*TEST 111 112 test: 113 suffix: 1 114 filter: grep -v " MPI process" 115 116 test: 117 suffix: 2 118 args: -m 1 -n 1 -fill 0.49 119 filter: grep -v " MPI process" 120 121 TEST*/ 122