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(MatCholeskyFactorSymbolic(F, A, perm, &info));
115: PetscCall(MatCholeskyFactorNumeric(F, A, &info));
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
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"));
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");
135: PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &B));
136: PetscCall(MatRealPart(B));
137: PetscCall(MatSetOption(B, MAT_SYMMETRIC, PETSC_TRUE));
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));
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: }
160: /*TEST
162: build:
163: requires: complex
165: test:
166: args: -n 1000 -check_symmetric -check_Hermitian
167: output_file: output/empty.out
169: test:
170: suffix: 2
171: nsize: 3
172: args: -n 1000
173: output_file: output/empty.out
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*/