Actual source code: ex127.c
1: static char help[] = "Test MatMult() for Hermitian matrix.\n\n";
3: #include <petscmat.h>
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;
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;
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));
29: sigma2 = 10.0 * PETSC_i;
30: h2 = 1.0 / ((n + 1) * (n + 1));
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));
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));
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: }
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: }
86: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
87: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
88: PetscCall(MatViewFromOptions(A, NULL, "-disp_mat"));
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));
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;
112: PetscCall(MatGetFactor(A, MATSOLVERSUPERLU_DIST, MAT_FACTOR_CHOLESKY, &F));
113: PetscCall(MatGetOrdering(A, MATORDERINGND, &perm, &iperm));
114: PetscCall(MatFactorInfoInitialize(&info));
115: PetscCall(MatCholeskyFactorSymbolic(F, A, perm, &info));
116: PetscCall(MatCholeskyFactorNumeric(F, A, &info));
118: PetscCall(MatGetInertia(F, &nneg, &nzero, &npos));
119: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " MatInertia: nneg: %" PetscInt_FMT ", nzero: %" PetscInt_FMT ", npos: %" PetscInt_FMT "\n", nneg, nzero, npos));
120: PetscCall(MatDestroy(&F));
121: PetscCall(ISDestroy(&perm));
122: PetscCall(ISDestroy(&iperm));
123: }
124: #endif
126: /* Create a Hermitian matrix As in sbaij format */
127: PetscCall(MatConvert(A, MATSBAIJ, MAT_INITIAL_MATRIX, &As));
128: PetscCall(MatViewFromOptions(As, NULL, "-disp_mat"));
130: /* Test MatMult */
131: PetscCall(MatMultEqual(A, As, 10, &flg));
132: PetscCheck(flg, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "MatMult not equal");
133: PetscCall(MatMultAddEqual(A, As, 10, &flg));
134: PetscCheck(flg, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "MatMultAdd not equal");
136: PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &B));
137: PetscCall(MatRealPart(B));
138: PetscCall(MatSetOption(B, MAT_SYMMETRIC, PETSC_TRUE));
140: PetscCall(MatCreateVecs(A, &r, &l));
141: PetscCall(VecSetRandom(r, NULL));
142: PetscCall(VecCopy(r, l));
143: PetscCall(MatDiagonalScale(A, r, l));
144: PetscCall(PetscOptionsHasName(NULL, NULL, "-check_Hermitian", &flg));
145: if (flg && size == 1) {
146: PetscCall(MatIsHermitian(A, 0.0, &flg));
147: PetscCheck(!flg, PETSC_COMM_WORLD, PETSC_ERR_USER, "A is Hermitian");
148: }
149: PetscCall(MatSetOption(A, MAT_HERMITIAN, PETSC_FALSE));
151: /* Free spaces */
152: PetscCall(MatDestroy(&A));
153: PetscCall(MatDestroy(&As));
154: PetscCall(MatDestroy(&B));
155: PetscCall(VecDestroy(&r));
156: PetscCall(VecDestroy(&l));
157: PetscCall(PetscFinalize());
158: return 0;
159: }
161: /*TEST
163: build:
164: requires: complex
166: test:
167: args: -n 1000 -check_symmetric -check_Hermitian
168: output_file: output/empty.out
170: test:
171: suffix: 2
172: nsize: 3
173: args: -n 1000
174: output_file: output/empty.out
176: test:
177: suffix: superlu_dist
178: nsize: 3
179: requires: superlu_dist
180: args: -test_choleskyfactor -mat_superlu_dist_rowperm NOROWPERM
181: output_file: output/ex127_superlu_dist.out
182: TEST*/