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;
8: PetscBool flg;
9: PetscMPIInt size;
10: PetscInt i, j;
11: PetscScalar v, sigma2;
12: PetscReal h2, sigma1 = 100.0;
13: PetscInt dim, Ii, J, n = 3, rstart, rend;
15: PetscFunctionBeginUser;
16: PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
17: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
18: PetscCall(PetscOptionsGetReal(NULL, NULL, "-sigma1", &sigma1, NULL));
19: PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
20: dim = n * n;
22: PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
23: PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, dim, dim));
24: PetscCall(MatSetType(A, MATAIJ));
25: PetscCall(MatSetFromOptions(A));
26: PetscCall(MatSetUp(A));
28: sigma2 = 10.0 * PETSC_i;
29: h2 = 1.0 / ((n + 1) * (n + 1));
31: PetscCall(MatGetOwnershipRange(A, &rstart, &rend));
32: for (Ii = rstart; Ii < rend; Ii++) {
33: v = -1.0;
34: i = Ii / n;
35: j = Ii - i * n;
36: if (i > 0) {
37: J = Ii - n;
38: PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES));
39: }
40: if (i < n - 1) {
41: J = Ii + n;
42: PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES));
43: }
44: if (j > 0) {
45: J = Ii - 1;
46: PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES));
47: }
48: if (j < n - 1) {
49: J = Ii + 1;
50: PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES));
51: }
52: v = 4.0 - sigma1 * h2;
53: PetscCall(MatSetValues(A, 1, &Ii, 1, &Ii, &v, ADD_VALUES));
54: }
55: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
56: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
58: /* Check whether A is symmetric */
59: PetscCall(PetscOptionsHasName(NULL, NULL, "-check_symmetric", &flg));
60: if (flg) {
61: PetscCall(MatIsSymmetric(A, 0.0, &flg));
62: PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_USER, "A is not symmetric");
63: }
64: PetscCall(MatSetOption(A, MAT_SYMMETRIC, PETSC_TRUE));
66: /* make A complex Hermitian */
67: Ii = 0;
68: J = dim - 1;
69: if (Ii >= rstart && Ii < rend) {
70: v = sigma2 * h2; /* RealPart(v) = 0.0 */
71: PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES));
72: v = -sigma2 * h2;
73: PetscCall(MatSetValues(A, 1, &J, 1, &Ii, &v, ADD_VALUES));
74: }
76: Ii = dim - 2;
77: J = dim - 1;
78: if (Ii >= rstart && Ii < rend) {
79: v = sigma2 * h2; /* RealPart(v) = 0.0 */
80: PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES));
81: v = -sigma2 * h2;
82: PetscCall(MatSetValues(A, 1, &J, 1, &Ii, &v, ADD_VALUES));
83: }
85: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
86: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
87: PetscCall(MatViewFromOptions(A, NULL, "-disp_mat"));
89: /* Check whether A is Hermitian, then set A->hermitian flag */
90: PetscCall(PetscOptionsHasName(NULL, NULL, "-check_Hermitian", &flg));
91: if (flg && size == 1) {
92: PetscCall(MatIsHermitian(A, 0.0, &flg));
93: PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_USER, "A is not Hermitian");
94: }
95: PetscCall(MatSetOption(A, MAT_HERMITIAN, PETSC_TRUE));
97: #if defined(PETSC_HAVE_SUPERLU_DIST)
98: /* Test Cholesky factorization */
99: PetscCall(PetscOptionsHasName(NULL, NULL, "-test_choleskyfactor", &flg));
100: if (flg) {
101: Mat F;
102: IS perm, iperm;
103: MatFactorInfo info;
104: PetscInt nneg, nzero, npos;
106: PetscCall(MatGetFactor(A, MATSOLVERSUPERLU_DIST, MAT_FACTOR_CHOLESKY, &F));
107: PetscCall(MatGetOrdering(A, MATORDERINGND, &perm, &iperm));
108: PetscCall(MatCholeskyFactorSymbolic(F, A, perm, &info));
109: PetscCall(MatCholeskyFactorNumeric(F, A, &info));
111: PetscCall(MatGetInertia(F, &nneg, &nzero, &npos));
112: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " MatInertia: nneg: %" PetscInt_FMT ", nzero: %" PetscInt_FMT ", npos: %" PetscInt_FMT "\n", nneg, nzero, npos));
113: PetscCall(MatDestroy(&F));
114: PetscCall(ISDestroy(&perm));
115: PetscCall(ISDestroy(&iperm));
116: }
117: #endif
119: /* Create a Hermitian matrix As in sbaij format */
120: PetscCall(MatConvert(A, MATSBAIJ, MAT_INITIAL_MATRIX, &As));
121: PetscCall(MatViewFromOptions(As, NULL, "-disp_mat"));
123: /* Test MatMult */
124: PetscCall(MatMultEqual(A, As, 10, &flg));
125: PetscCheck(flg, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "MatMult not equal");
126: PetscCall(MatMultAddEqual(A, As, 10, &flg));
127: PetscCheck(flg, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "MatMultAdd not equal");
129: /* Free spaces */
130: PetscCall(MatDestroy(&A));
131: PetscCall(MatDestroy(&As));
132: PetscCall(PetscFinalize());
133: return 0;
134: }
136: /*TEST
138: build:
139: requires: complex
141: test:
142: args: -n 1000
143: output_file: output/ex127.out
145: test:
146: suffix: 2
147: nsize: 3
148: args: -n 1000
149: output_file: output/ex127.out
151: test:
152: suffix: superlu_dist
153: nsize: 3
154: requires: superlu_dist
155: args: -test_choleskyfactor -mat_superlu_dist_rowperm NOROWPERM
156: output_file: output/ex127_superlu_dist.out
157: TEST*/