Actual source code: ex127.c
petsc-3.3-p6 2013-02-11
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>
12: PetscInt main(PetscInt argc,char **args)
13: {
14: Mat A,As;
15: Vec x,y,ys;
16: PetscBool flg,disp_mat=PETSC_FALSE,disp_vec=PETSC_FALSE;
18: PetscMPIInt size,rank;
19: PetscInt m,i,j;
20: PetscScalar v,sigma2;
21: PetscRandom rctx;
22: PetscReal h2,sigma1=100.0,norm;
23: PetscInt dim,Ii,J,n = 3,use_random,rstart,rend;
24:
25: PetscInitialize(&argc,&args,(char *)0,help);
26: #if !defined(PETSC_USE_COMPLEX)
27: SETERRQ(PETSC_COMM_WORLD,1,"This example requires complex numbers");
28: #endif
29: MPI_Comm_size(PETSC_COMM_WORLD,&size);
30: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
31: PetscOptionsHasName(PETSC_NULL, "-display_mat", &disp_mat);
32: PetscOptionsHasName(PETSC_NULL, "-display_vec", &disp_vec);
34: PetscOptionsGetReal(PETSC_NULL,"-sigma1",&sigma1,PETSC_NULL);
35: PetscOptionsGetInt(PETSC_NULL,"-n",&n,PETSC_NULL);
36: dim = n*n;
38: MatCreate(PETSC_COMM_WORLD,&A);
39: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,dim,dim);
40: MatSetType(A,MATAIJ);
41: MatSetFromOptions(A);
43: PetscOptionsHasName(PETSC_NULL,"-norandom",&flg);
44: if (flg) use_random = 0;
45: else use_random = 1;
46: if (use_random) {
47: PetscRandomCreate(PETSC_COMM_WORLD,&rctx);
48: PetscRandomSetFromOptions(rctx);
49: PetscRandomSetInterval(rctx,0.0,PETSC_i);
50: PetscRandomGetValue(rctx,&sigma2); /* RealPart(sigma2) == 0.0 */
51: } else {
52: sigma2 = 10.0*PETSC_i;
53: }
54: h2 = 1.0/((n+1)*(n+1));
56: MatGetOwnershipRange(A,&rstart,&rend);
57: for (Ii=rstart; Ii<rend; Ii++) {
58: v = -1.0; i = Ii/n; j = Ii - i*n;
59: if (i>0) {
60: J = Ii-n; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);}
61: if (i<n-1) {
62: J = Ii+n; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);}
63: if (j>0) {
64: J = Ii-1; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);}
65: if (j<n-1) {
66: J = Ii+1; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);}
67: v = 4.0 - sigma1*h2;
68: MatSetValues(A,1,&Ii,1,&Ii,&v,ADD_VALUES);
69: }
70: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
71: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
73: /* Check whether A is symmetric */
74: PetscOptionsHasName(PETSC_NULL, "-check_symmetric", &flg);
75: if (flg) {
76: Mat Trans;
77: MatTranspose(A,MAT_INITIAL_MATRIX, &Trans);
78: MatEqual(A, Trans, &flg);
79: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"A is not symmetric");
80: MatDestroy(&Trans);
81: }
82: MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);
84: /* make A complex Hermitian */
85: Ii = 0; J = dim-1;
86: if (Ii >= rstart && Ii < rend){
87: v = sigma2*h2; /* RealPart(v) = 0.0 */
88: MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);
89: v = -sigma2*h2;
90: MatSetValues(A,1,&J,1,&Ii,&v,ADD_VALUES);
91: }
93: Ii = dim-2; J = dim-1;
94: if (Ii >= rstart && Ii < rend){
95: v = sigma2*h2; /* RealPart(v) = 0.0 */
96: MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);
97: v = -sigma2*h2;
98: MatSetValues(A,1,&J,1,&Ii,&v,ADD_VALUES);
99: }
101: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
102: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
104: /* Check whether A is Hermitian */
105: PetscOptionsHasName(PETSC_NULL, "-check_Hermitian", &flg);
106: if (flg) {
107: Mat Hermit;
108: if (disp_mat){
109: if (!rank) printf(" A:\n");
110: MatView(A,PETSC_VIEWER_STDOUT_WORLD);
111: }
112: MatHermitianTranspose(A,MAT_INITIAL_MATRIX, &Hermit);
113: if (disp_mat){
114: if (!rank) printf(" A_Hermitian:\n");
115: MatView(Hermit,PETSC_VIEWER_STDOUT_WORLD);
116: }
117: MatEqual(A, Hermit, &flg);
118: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"A is not Hermitian");
119: MatDestroy(&Hermit);
120: }
121: MatSetOption(A,MAT_HERMITIAN,PETSC_TRUE);
122:
123: /* Create a Hermitian matrix As in sbaij format */
124: MatConvert(A,MATSBAIJ,MAT_INITIAL_MATRIX,&As);
125: if (disp_mat){
126: if (!rank) {PetscPrintf(PETSC_COMM_SELF," As:\n");}
127: MatView(As,PETSC_VIEWER_STDOUT_WORLD);
128: }
130: MatGetLocalSize(A,&m,&n);
131: VecCreate(PETSC_COMM_WORLD,&x);
132: VecSetSizes(x,n,PETSC_DECIDE);
133: VecSetFromOptions(x);
134: if (use_random){
135: VecSetRandom(x,rctx);
136: } else {
137: VecSet(x,1.0);
138: }
140: /* Create vectors */
141: VecCreate(PETSC_COMM_WORLD,&y);
142: VecSetSizes(y,m,PETSC_DECIDE);
143: VecSetFromOptions(y);
144: VecDuplicate(y,&ys);
146: /* Test MatMult */
147: MatMult(A,x,y);
148: MatMult(As,x,ys);
149: if (disp_vec){
150: printf("y = A*x:\n");
151: VecView(y,PETSC_VIEWER_STDOUT_WORLD);
152: PetscPrintf(PETSC_COMM_WORLD,"ys = As*x:\n");
153: VecView(ys,PETSC_VIEWER_STDOUT_WORLD);
154: }
155: VecAXPY(y,-1.0,ys);
156: VecNorm(y,NORM_INFINITY,&norm);
157: if (norm > 1.e-12 || disp_vec){
158: printf("|| A*x - As*x || = %G\n",norm);
159: }
161: /* Free spaces */
162: if (use_random) {PetscRandomDestroy(&rctx);}
163: MatDestroy(&A);
164: MatDestroy(&As);
165:
166: VecDestroy(&x);
167: VecDestroy(&y);
168: VecDestroy(&ys);
169: PetscFinalize();
170: return 0;
171: }