Actual source code: ex275.c
1: static char help[] = "Tests MatNormApproximate() on the Brusselator matrix.\n\n"
2: "The command line options are:\n"
3: " -n <n>, where <n> = block dimension of the 2x2 block matrix.\n"
4: " -L <L>, where <L> = bifurcation parameter.\n"
5: " -alpha <alpha>, -beta <beta>, -delta1 <delta1>, -delta2 <delta2>,\n"
6: " where <alpha> <beta> <delta1> <delta2> = model parameters.\n\n";
8: #include <petscmat.h>
10: /*
11: The Brusselator matrix is
13: A = [ tau1*T+(beta-1)*I alpha^2*I
14: -beta*I tau2*T-alpha^2*I ],
16: where
18: T = tridiag{1,-2,1}
19: h = 1/(n+1)
20: tau1 = delta1/(h*L)^2
21: tau2 = delta2/(h*L)^2
22: */
24: int main(int argc, char **argv)
25: {
26: Mat A, T1, T2, D1, D2, mats[4], Ae;
27: PetscScalar alpha, beta, tau1, tau2, delta1, delta2, L, h;
28: PetscReal norm[3];
29: PetscInt N = 30, i, Istart, Iend, samples = PETSC_DECIDE;
30: PetscBool test_fwd = PETSC_FALSE;
32: PetscFunctionBeginUser;
33: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
34: PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &N, NULL));
36: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
37: Generate the matrix
38: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
40: alpha = 2.0;
41: beta = 5.45;
42: delta1 = 0.008;
43: delta2 = 0.004;
44: L = 0.51302;
46: PetscCall(PetscOptionsGetScalar(NULL, NULL, "-L", &L, NULL));
47: PetscCall(PetscOptionsGetScalar(NULL, NULL, "-alpha", &alpha, NULL));
48: PetscCall(PetscOptionsGetScalar(NULL, NULL, "-beta", &beta, NULL));
49: PetscCall(PetscOptionsGetScalar(NULL, NULL, "-delta1", &delta1, NULL));
50: PetscCall(PetscOptionsGetScalar(NULL, NULL, "-delta2", &delta2, NULL));
51: PetscCall(PetscOptionsGetInt(NULL, NULL, "-samples", &samples, NULL));
52: PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_fwd", &test_fwd, NULL));
54: h = 1.0 / (PetscReal)(N + 1);
55: tau1 = delta1 / ((h * L) * (h * L));
56: tau2 = delta2 / ((h * L) * (h * L));
58: /* Create matrices T1, T2 */
59: PetscCall(MatCreate(PETSC_COMM_WORLD, &T1));
60: PetscCall(MatSetSizes(T1, PETSC_DECIDE, PETSC_DECIDE, N, N));
61: PetscCall(MatSetFromOptions(T1));
62: PetscCall(MatGetOwnershipRange(T1, &Istart, &Iend));
63: for (i = Istart; i < Iend; i++) {
64: if (i > 0) PetscCall(MatSetValue(T1, i, i - 1, 1.0, INSERT_VALUES));
65: if (i < N - 1) PetscCall(MatSetValue(T1, i, i + 1, 1.0, INSERT_VALUES));
66: PetscCall(MatSetValue(T1, i, i, -2.0, INSERT_VALUES));
67: }
68: PetscCall(MatAssemblyBegin(T1, MAT_FINAL_ASSEMBLY));
69: PetscCall(MatAssemblyEnd(T1, MAT_FINAL_ASSEMBLY));
71: PetscCall(MatDuplicate(T1, MAT_COPY_VALUES, &T2));
72: PetscCall(MatScale(T1, tau1));
73: PetscCall(MatShift(T1, beta - 1.0));
74: PetscCall(MatScale(T2, tau2));
75: PetscCall(MatShift(T2, -alpha * alpha));
77: /* Create matrices D1, D2 */
78: PetscCall(MatCreateConstantDiagonal(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, N, N, alpha * alpha, &D1));
79: PetscCall(MatCreateConstantDiagonal(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, N, N, -beta, &D2));
81: /* Create the nest matrix */
82: mats[0] = T1;
83: mats[1] = D1;
84: mats[2] = D2;
85: mats[3] = T2;
86: PetscCall(MatCreateNest(PETSC_COMM_WORLD, 2, NULL, 2, NULL, mats, &A));
87: if (test_fwd) PetscCall(MatSetOperation(A, MATOP_MULT_TRANSPOSE, NULL));
89: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
90: Estimate the norm
91: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
93: PetscCall(MatComputeOperator(A, MATDENSE, &Ae));
94: PetscCall(MatViewFromOptions(Ae, NULL, "-view_explicit"));
96: PetscCall(MatNormApproximate(A, NORM_1, samples, norm));
97: PetscCall(MatNormApproximate(A, NORM_INFINITY, samples, norm + 1));
98: PetscCall(MatNormApproximate(A, NORM_2, samples, norm + 2));
99: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nBrusselator matrix, n=%" PetscInt_FMT " - estimated 1-norm = %g - estimated infinity-norm = %g - estimated 2-norm = %g\n\n", N, (double)norm[0], (double)norm[1], (double)norm[2]));
101: PetscCall(MatNorm(Ae, NORM_1, norm));
102: PetscCall(MatNorm(Ae, NORM_INFINITY, norm + 1));
103: PetscCall(MatNorm(Ae, NORM_2, norm + 2));
104: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nBrusselator matrix, n=%" PetscInt_FMT " - exact 1-norm = %g - exact infinity-norm = %g - exact 2-norm = %g\n\n", N, (double)norm[0], (double)norm[1], (double)norm[2]));
106: PetscCall(MatDestroy(&A));
107: PetscCall(MatDestroy(&Ae));
108: PetscCall(MatDestroy(&T1));
109: PetscCall(MatDestroy(&T2));
110: PetscCall(MatDestroy(&D1));
111: PetscCall(MatDestroy(&D2));
112: PetscCall(PetscFinalize());
113: return 0;
114: }
116: /*TEST
118: test:
119: suffix: 1
120: args: -test_fwd {{false true}shared output}
121: output_file: output/ex275_1.out
123: TEST*/