1 static char help[] = "Test MatConvert() with a MATNEST with scaled and shifted MATTRANSPOSEVIRTUAL blocks.\n\n";
2
3 #include <petscmat.h>
4
5 /*
6 This example builds the matrix
7
8 H = [ R C
9 alpha C^H + beta I gamma R^T + delta I ],
10
11 where R is Hermitian and C is complex symmetric. In particular, R and C have the
12 following Toeplitz structure:
13
14 R = pentadiag{a,b,c,conj(b),conj(a)}
15 C = tridiag{b,d,b}
16
17 where a,b,d are complex scalars, and c is real.
18 */
19
main(int argc,char ** argv)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;
26
27 PetscFunctionBeginUser;
28 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
29
30 PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
31
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);
36
37 PetscCall(MatCreate(PETSC_COMM_WORLD, &R));
38 PetscCall(MatSetSizes(R, PETSC_DECIDE, PETSC_DECIDE, n, n));
39 PetscCall(MatSetFromOptions(R));
40
41 PetscCall(MatCreate(PETSC_COMM_WORLD, &C));
42 PetscCall(MatSetSizes(C, PETSC_DECIDE, PETSC_DECIDE, n, n));
43 PetscCall(MatSetFromOptions(C));
44
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 }
53
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 }
60
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));
65
66 block[0] = R;
67 block[1] = C;
68
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]));
78
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");
82
83 PetscCall(MatDestroy(&R));
84 PetscCall(MatDestroy(&C));
85 PetscCall(MatDestroy(&H));
86 PetscCall(MatDestroy(&M));
87 PetscCall(PetscFinalize());
88 return 0;
89 }
90
91 /*TEST
92
93 build:
94 requires: complex
95
96 test:
97 output_file: output/empty.out
98 nsize: {{1 4}}
99
100 TEST*/
101