Actual source code: ex264.c
1: static char help[] = "Test MatConvert() with a MATNEST with scaled and shifted MATTRANSPOSEVIRTUAL blocks.\n\n";
3: #include <petscmat.h>
5: /*
6: This example builds the matrix
8: H = [ R C
9: alpha C^H + beta I gamma R^T + delta I ],
11: where R is Hermitian and C is complex symmetric. In particular, R and C have the
12: following Toeplitz structure:
14: R = pentadiag{a,b,c,conj(b),conj(a)}
15: C = tridiag{b,d,b}
17: where a,b,d are complex scalars, and c is real.
18: */
20: int main(int argc, char **argv)
21: {
22: Mat block[4], H, R, C, M;
23: PetscScalar a, b, c, d;
24: PetscInt n = 13, Istart, Iend, i;
25: PetscBool flg;
27: PetscFunctionBeginUser;
28: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
30: PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
32: a = PetscCMPLX(-0.1, 0.2);
33: b = PetscCMPLX(1.0, 0.5);
34: c = 4.5;
35: d = PetscCMPLX(2.0, 0.2);
37: PetscCall(MatCreate(PETSC_COMM_WORLD, &R));
38: PetscCall(MatSetSizes(R, PETSC_DECIDE, PETSC_DECIDE, n, n));
39: PetscCall(MatSetFromOptions(R));
41: PetscCall(MatCreate(PETSC_COMM_WORLD, &C));
42: PetscCall(MatSetSizes(C, PETSC_DECIDE, PETSC_DECIDE, n, n));
43: PetscCall(MatSetFromOptions(C));
45: PetscCall(MatGetOwnershipRange(R, &Istart, &Iend));
46: for (i = Istart; i < Iend; i++) {
47: if (i > 1) PetscCall(MatSetValue(R, i, i - 2, a, INSERT_VALUES));
48: if (i > 0) PetscCall(MatSetValue(R, i, i - 1, b, INSERT_VALUES));
49: PetscCall(MatSetValue(R, i, i, c, INSERT_VALUES));
50: if (i < n - 1) PetscCall(MatSetValue(R, i, i + 1, PetscConj(b), INSERT_VALUES));
51: if (i < n - 2) PetscCall(MatSetValue(R, i, i + 2, PetscConj(a), INSERT_VALUES));
52: }
54: PetscCall(MatGetOwnershipRange(C, &Istart, &Iend));
55: for (i = Istart; i < Iend; i++) {
56: if (i > 0) PetscCall(MatSetValue(C, i, i - 1, b, INSERT_VALUES));
57: PetscCall(MatSetValue(C, i, i, d, INSERT_VALUES));
58: if (i < n - 1) PetscCall(MatSetValue(C, i, i + 1, b, INSERT_VALUES));
59: }
61: PetscCall(MatAssemblyBegin(R, MAT_FINAL_ASSEMBLY));
62: PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
63: PetscCall(MatAssemblyEnd(R, MAT_FINAL_ASSEMBLY));
64: PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
66: block[0] = R;
67: block[1] = C;
69: PetscCall(MatCreateHermitianTranspose(C, &block[2]));
70: PetscCall(MatScale(block[2], PetscConj(b)));
71: PetscCall(MatShift(block[2], d));
72: PetscCall(MatCreateTranspose(R, &block[3]));
73: PetscCall(MatScale(block[3], PetscConj(d)));
74: PetscCall(MatShift(block[3], b));
75: PetscCall(MatCreateNest(PetscObjectComm((PetscObject)R), 2, NULL, 2, NULL, block, &H));
76: PetscCall(MatDestroy(&block[2]));
77: PetscCall(MatDestroy(&block[3]));
79: PetscCall(MatConvert(H, MATAIJ, MAT_INITIAL_MATRIX, &M));
80: PetscCall(MatMultEqual(H, M, 20, &flg));
81: PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "MatNest != MatAIJ");
83: PetscCall(MatDestroy(&R));
84: PetscCall(MatDestroy(&C));
85: PetscCall(MatDestroy(&H));
86: PetscCall(MatDestroy(&M));
87: PetscCall(PetscFinalize());
88: return 0;
89: }
91: /*TEST
93: build:
94: requires: complex
96: test:
97: output_file: output/ex109.out
98: nsize: {{1 4}}
100: TEST*/