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;
16: PetscInitialize(&argc, &args, (char *)0, help);
17: MPI_Comm_size(PETSC_COMM_WORLD, &size);
18: PetscOptionsGetReal(NULL, NULL, "-sigma1", &sigma1, NULL);
19: PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL);
20: dim = n * n;
22: MatCreate(PETSC_COMM_WORLD, &A);
23: MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, dim, dim);
24: MatSetType(A, MATAIJ);
25: MatSetFromOptions(A);
26: MatSetUp(A);
28: sigma2 = 10.0 * PETSC_i;
29: h2 = 1.0 / ((n + 1) * (n + 1));
31: 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: MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES);
39: }
40: if (i < n - 1) {
41: J = Ii + n;
42: MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES);
43: }
44: if (j > 0) {
45: J = Ii - 1;
46: MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES);
47: }
48: if (j < n - 1) {
49: J = Ii + 1;
50: MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES);
51: }
52: v = 4.0 - sigma1 * h2;
53: MatSetValues(A, 1, &Ii, 1, &Ii, &v, ADD_VALUES);
54: }
55: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
56: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
58: /* Check whether A is symmetric */
59: PetscOptionsHasName(NULL, NULL, "-check_symmetric", &flg);
60: if (flg) {
61: MatIsSymmetric(A, 0.0, &flg);
63: }
64: 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: MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES);
72: v = -sigma2 * h2;
73: 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: MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES);
81: v = -sigma2 * h2;
82: MatSetValues(A, 1, &J, 1, &Ii, &v, ADD_VALUES);
83: }
85: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
86: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
87: MatViewFromOptions(A, NULL, "-disp_mat");
89: /* Check whether A is Hermitian, then set A->hermitian flag */
90: PetscOptionsHasName(NULL, NULL, "-check_Hermitian", &flg);
91: if (flg && size == 1) {
92: MatIsHermitian(A, 0.0, &flg);
94: }
95: MatSetOption(A, MAT_HERMITIAN, PETSC_TRUE);
97: #if defined(PETSC_HAVE_SUPERLU_DIST)
98: /* Test Cholesky factorization */
99: 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: MatGetFactor(A, MATSOLVERSUPERLU_DIST, MAT_FACTOR_CHOLESKY, &F);
107: MatGetOrdering(A, MATORDERINGND, &perm, &iperm);
108: MatCholeskyFactorSymbolic(F, A, perm, &info);
109: MatCholeskyFactorNumeric(F, A, &info);
111: MatGetInertia(F, &nneg, &nzero, &npos);
112: PetscPrintf(PETSC_COMM_WORLD, " MatInertia: nneg: %" PetscInt_FMT ", nzero: %" PetscInt_FMT ", npos: %" PetscInt_FMT "\n", nneg, nzero, npos);
113: MatDestroy(&F);
114: ISDestroy(&perm);
115: ISDestroy(&iperm);
116: }
117: #endif
119: /* Create a Hermitian matrix As in sbaij format */
120: MatConvert(A, MATSBAIJ, MAT_INITIAL_MATRIX, &As);
121: MatViewFromOptions(As, NULL, "-disp_mat");
123: /* Test MatMult */
124: MatMultEqual(A, As, 10, &flg);
126: MatMultAddEqual(A, As, 10, &flg);
129: /* Free spaces */
130: MatDestroy(&A);
131: MatDestroy(&As);
132: 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*/