1 static char help[] = "Test MatMult() for Hermitian matrix.\n\n";
2
3 #include <petscmat.h>
4
main(int argc,char ** args)5 int main(int argc, char **args)
6 {
7 Mat A, As, B;
8 Vec l, r;
9 PetscBool flg;
10 PetscMPIInt size;
11 PetscInt i, j;
12 PetscScalar v, sigma2;
13 PetscReal h2, sigma1 = 100.0;
14 PetscInt dim, Ii, J, n = 3, rstart, rend;
15
16 PetscFunctionBeginUser;
17 PetscCall(PetscInitialize(&argc, &args, NULL, help));
18 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
19 PetscCall(PetscOptionsGetReal(NULL, NULL, "-sigma1", &sigma1, NULL));
20 PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
21 dim = n * n;
22
23 PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
24 PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, dim, dim));
25 PetscCall(MatSetType(A, MATAIJ));
26 PetscCall(MatSetFromOptions(A));
27 PetscCall(MatSetUp(A));
28
29 sigma2 = 10.0 * PETSC_i;
30 h2 = 1.0 / ((n + 1) * (n + 1));
31
32 PetscCall(MatGetOwnershipRange(A, &rstart, &rend));
33 for (Ii = rstart; Ii < rend; Ii++) {
34 v = -1.0;
35 i = Ii / n;
36 j = Ii - i * n;
37 if (i > 0) {
38 J = Ii - n;
39 PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES));
40 }
41 if (i < n - 1) {
42 J = Ii + n;
43 PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES));
44 }
45 if (j > 0) {
46 J = Ii - 1;
47 PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES));
48 }
49 if (j < n - 1) {
50 J = Ii + 1;
51 PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES));
52 }
53 v = 4.0 - sigma1 * h2;
54 PetscCall(MatSetValues(A, 1, &Ii, 1, &Ii, &v, ADD_VALUES));
55 }
56 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
57 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
58
59 /* Check whether A is symmetric */
60 PetscCall(PetscOptionsHasName(NULL, NULL, "-check_symmetric", &flg));
61 if (flg) {
62 PetscCall(MatIsSymmetric(A, 0.0, &flg));
63 PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_USER, "A is not symmetric");
64 }
65 PetscCall(MatSetOption(A, MAT_SYMMETRIC, PETSC_TRUE));
66
67 /* make A complex Hermitian */
68 Ii = 0;
69 J = dim - 1;
70 if (Ii >= rstart && Ii < rend) {
71 v = sigma2 * h2; /* RealPart(v) = 0.0 */
72 PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES));
73 v = -sigma2 * h2;
74 PetscCall(MatSetValues(A, 1, &J, 1, &Ii, &v, ADD_VALUES));
75 }
76
77 Ii = dim - 2;
78 J = dim - 1;
79 if (Ii >= rstart && Ii < rend) {
80 v = sigma2 * h2; /* RealPart(v) = 0.0 */
81 PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES));
82 v = -sigma2 * h2;
83 PetscCall(MatSetValues(A, 1, &J, 1, &Ii, &v, ADD_VALUES));
84 }
85
86 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
87 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
88 PetscCall(MatViewFromOptions(A, NULL, "-disp_mat"));
89
90 if (flg) {
91 PetscCall(MatIsSymmetric(A, 0.0, &flg));
92 PetscCheck(!flg, PETSC_COMM_WORLD, PETSC_ERR_USER, "A is symmetric");
93 }
94 PetscCall(MatSetOption(A, MAT_SYMMETRIC, PETSC_FALSE));
95 /* Check whether A is Hermitian, then set A->hermitian flag */
96 PetscCall(PetscOptionsHasName(NULL, NULL, "-check_Hermitian", &flg));
97 if (flg && size == 1) {
98 PetscCall(MatIsHermitian(A, 0.0, &flg));
99 PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_USER, "A is not Hermitian");
100 }
101 PetscCall(MatSetOption(A, MAT_HERMITIAN, PETSC_TRUE));
102
103 #if defined(PETSC_HAVE_SUPERLU_DIST)
104 /* Test Cholesky factorization */
105 PetscCall(PetscOptionsHasName(NULL, NULL, "-test_choleskyfactor", &flg));
106 if (flg) {
107 Mat F;
108 IS perm, iperm;
109 MatFactorInfo info;
110 PetscInt nneg, nzero, npos;
111
112 PetscCall(MatGetFactor(A, MATSOLVERSUPERLU_DIST, MAT_FACTOR_CHOLESKY, &F));
113 PetscCall(MatGetOrdering(A, MATORDERINGND, &perm, &iperm));
114 PetscCall(MatCholeskyFactorSymbolic(F, A, perm, &info));
115 PetscCall(MatCholeskyFactorNumeric(F, A, &info));
116
117 PetscCall(MatGetInertia(F, &nneg, &nzero, &npos));
118 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " MatInertia: nneg: %" PetscInt_FMT ", nzero: %" PetscInt_FMT ", npos: %" PetscInt_FMT "\n", nneg, nzero, npos));
119 PetscCall(MatDestroy(&F));
120 PetscCall(ISDestroy(&perm));
121 PetscCall(ISDestroy(&iperm));
122 }
123 #endif
124
125 /* Create a Hermitian matrix As in sbaij format */
126 PetscCall(MatConvert(A, MATSBAIJ, MAT_INITIAL_MATRIX, &As));
127 PetscCall(MatViewFromOptions(As, NULL, "-disp_mat"));
128
129 /* Test MatMult */
130 PetscCall(MatMultEqual(A, As, 10, &flg));
131 PetscCheck(flg, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "MatMult not equal");
132 PetscCall(MatMultAddEqual(A, As, 10, &flg));
133 PetscCheck(flg, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "MatMultAdd not equal");
134
135 PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &B));
136 PetscCall(MatRealPart(B));
137 PetscCall(MatSetOption(B, MAT_SYMMETRIC, PETSC_TRUE));
138
139 PetscCall(MatCreateVecs(A, &r, &l));
140 PetscCall(VecSetRandom(r, NULL));
141 PetscCall(VecCopy(r, l));
142 PetscCall(MatDiagonalScale(A, r, l));
143 PetscCall(PetscOptionsHasName(NULL, NULL, "-check_Hermitian", &flg));
144 if (flg && size == 1) {
145 PetscCall(MatIsHermitian(A, 0.0, &flg));
146 PetscCheck(!flg, PETSC_COMM_WORLD, PETSC_ERR_USER, "A is Hermitian");
147 }
148 PetscCall(MatSetOption(A, MAT_HERMITIAN, PETSC_FALSE));
149
150 /* Free spaces */
151 PetscCall(MatDestroy(&A));
152 PetscCall(MatDestroy(&As));
153 PetscCall(MatDestroy(&B));
154 PetscCall(VecDestroy(&r));
155 PetscCall(VecDestroy(&l));
156 PetscCall(PetscFinalize());
157 return 0;
158 }
159
160 /*TEST
161
162 build:
163 requires: complex
164
165 test:
166 args: -n 1000 -check_symmetric -check_Hermitian
167 output_file: output/empty.out
168
169 test:
170 suffix: 2
171 nsize: 3
172 args: -n 1000
173 output_file: output/empty.out
174
175 test:
176 suffix: superlu_dist
177 nsize: 3
178 requires: superlu_dist
179 args: -test_choleskyfactor -mat_superlu_dist_rowperm NOROWPERM
180 output_file: output/ex127_superlu_dist.out
181 TEST*/
182