Actual source code: ex127.c
petsc-3.12.5 2020-03-29
1: static char help[] = "Test MatMult() for Hermitian matrix.\n\n";
2: /*
3: Example of usage
4: ./ex127 -check_Hermitian -display_mat -display_vec
5: mpiexec -n 2 ./ex127
6: */
8: #include <petscmat.h>
10: int main(int argc,char **args)
11: {
12: Mat A,As;
13: Vec x,y,ys;
14: PetscBool flg,disp_mat=PETSC_FALSE,disp_vec=PETSC_FALSE;
16: PetscMPIInt size,rank;
17: PetscInt m,i,j;
18: PetscScalar v,sigma2;
19: PetscRandom rctx = NULL;
20: PetscReal h2,sigma1=100.0,norm;
21: PetscInt dim,Ii,J,n = 3,rstart,rend;
23: PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
24: #if !defined(PETSC_USE_COMPLEX)
25: SETERRQ(PETSC_COMM_WORLD,1,"This example requires complex build");
26: #endif
28: MPI_Comm_size(PETSC_COMM_WORLD,&size);
29: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
30: PetscOptionsHasName(NULL,NULL, "-display_mat", &disp_mat);
31: PetscOptionsHasName(NULL,NULL, "-display_vec", &disp_vec);
33: PetscOptionsGetReal(NULL,NULL,"-sigma1",&sigma1,NULL);
34: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
35: dim = n*n;
37: MatCreate(PETSC_COMM_WORLD,&A);
38: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,dim,dim);
39: MatSetType(A,MATAIJ);
40: MatSetFromOptions(A);
41: MatSetUp(A);
43: PetscOptionsHasName(NULL,NULL,"-norandom",&flg);
44: if (!flg) {
45: PetscRandomCreate(PETSC_COMM_WORLD,&rctx);
46: PetscRandomSetFromOptions(rctx);
47: PetscRandomSetInterval(rctx,0.0,PETSC_i);
48: PetscRandomGetValue(rctx,&sigma2); /* RealPart(sigma2) == 0.0 */
49: } else {
50: sigma2 = 10.0*PETSC_i;
51: }
52: h2 = 1.0/((n+1)*(n+1));
54: MatGetOwnershipRange(A,&rstart,&rend);
55: for (Ii=rstart; Ii<rend; Ii++) {
56: v = -1.0; i = Ii/n; j = Ii - i*n;
57: if (i>0) {
58: J = Ii-n; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);
59: }
60: if (i<n-1) {
61: J = Ii+n; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);
62: }
63: if (j>0) {
64: J = Ii-1; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);
65: }
66: if (j<n-1) {
67: J = Ii+1; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);
68: }
69: v = 4.0 - sigma1*h2;
70: MatSetValues(A,1,&Ii,1,&Ii,&v,ADD_VALUES);
71: }
72: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
73: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
75: /* Check whether A is symmetric */
76: PetscOptionsHasName(NULL,NULL, "-check_symmetric", &flg);
77: if (flg) {
78: MatIsSymmetric(A,0.0,&flg);
79: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"A is not symmetric");
80: }
81: MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);
83: /* make A complex Hermitian */
84: Ii = 0; J = dim-1;
85: if (Ii >= rstart && Ii < rend) {
86: v = sigma2*h2; /* RealPart(v) = 0.0 */
87: MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);
88: v = -sigma2*h2;
89: MatSetValues(A,1,&J,1,&Ii,&v,ADD_VALUES);
90: }
92: Ii = dim-2; J = dim-1;
93: if (Ii >= rstart && Ii < rend) {
94: v = sigma2*h2; /* RealPart(v) = 0.0 */
95: MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);
96: v = -sigma2*h2;
97: MatSetValues(A,1,&J,1,&Ii,&v,ADD_VALUES);
98: }
100: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
101: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
103: /* Check whether A is Hermitian, then set A->hermitian flag */
104: PetscOptionsHasName(NULL,NULL, "-check_Hermitian", &flg);
105: if (flg && size == 1) {
106: MatIsHermitian(A,0.0,&flg);
107: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"A is not Hermitian");
108: }
109: MatSetOption(A,MAT_HERMITIAN,PETSC_TRUE);
111: #if defined(PETSC_HAVE_SUPERLU_DIST)
112: /* Test Cholesky factorization */
113: PetscOptionsHasName(NULL,NULL, "-test_choleskyfactor", &flg);
114: if (flg) {
115: Mat F;
116: IS perm,iperm;
117: MatFactorInfo info;
118: PetscInt nneg,nzero,npos;
120: MatGetFactor(A,MATSOLVERSUPERLU_DIST,MAT_FACTOR_CHOLESKY,&F);
121: MatGetOrdering(A,MATORDERINGND,&perm,&iperm);
122: MatCholeskyFactorSymbolic(F,A,perm,&info);
123: MatCholeskyFactorNumeric(F,A,&info);
125: MatGetInertia(F,&nneg,&nzero,&npos);
126: if (!rank) {
127: PetscPrintf(PETSC_COMM_SELF," MatInertia: nneg: %D, nzero: %D, npos: %D\n",nneg,nzero,npos);
128: }
129: MatDestroy(&F);
130: ISDestroy(&perm);
131: ISDestroy(&iperm);
132: }
133: #endif
135: /* Create a Hermitian matrix As in sbaij format */
136: MatConvert(A,MATSBAIJ,MAT_INITIAL_MATRIX,&As);
137: if (disp_mat) {
138: PetscPrintf(PETSC_COMM_WORLD," As:\n");
139: MatView(As,PETSC_VIEWER_STDOUT_WORLD);
140: }
142: MatGetLocalSize(A,&m,&n);
143: VecCreate(PETSC_COMM_WORLD,&x);
144: VecSetSizes(x,n,PETSC_DECIDE);
145: VecSetFromOptions(x);
146: if (rctx) {
147: VecSetRandom(x,rctx);
148: } else {
149: VecSet(x,1.0);
150: }
152: /* Create vectors */
153: VecCreate(PETSC_COMM_WORLD,&y);
154: VecSetSizes(y,m,PETSC_DECIDE);
155: VecSetFromOptions(y);
156: VecDuplicate(y,&ys);
158: /* Test MatMult */
159: MatMult(A,x,y);
160: MatMult(As,x,ys);
161: if (disp_vec) {
162: PetscPrintf(PETSC_COMM_WORLD,"y = A*x:\n");
163: VecView(y,PETSC_VIEWER_STDOUT_WORLD);
164: PetscPrintf(PETSC_COMM_WORLD,"ys = As*x:\n");
165: VecView(ys,PETSC_VIEWER_STDOUT_WORLD);
166: }
167: VecAXPY(y,-1.0,ys);
168: VecNorm(y,NORM_INFINITY,&norm);
169: if (norm > 100.0*PETSC_MACHINE_EPSILON || disp_vec) {
170: PetscPrintf(PETSC_COMM_WORLD,"|| A*x - As*x || = %g\n",(double)norm);
171: }
173: /* Free spaces */
174: PetscRandomDestroy(&rctx);
175: MatDestroy(&A);
176: MatDestroy(&As);
178: VecDestroy(&x);
179: VecDestroy(&y);
180: VecDestroy(&ys);
181: PetscFinalize();
182: return ierr;
183: }
186: /*TEST
188: build:
189: requires: complex
191: test:
192: args: -n 1000
193: output_file: output/ex127.out
195: test:
196: suffix: 2
197: nsize: 3
198: args: -n 1000
199: output_file: output/ex127.out
202: test:
203: suffix: superlu_dist
204: nsize: 3
205: requires: superlu_dist
206: args: -test_choleskyfactor -mat_superlu_dist_rowperm NOROWPERM
207: output_file: output/ex127_superlu_dist.out
208: TEST*/