Actual source code: ex36.c
petsc-3.3-p6 2013-02-11
1: static char help[] = "Test MatGetInertia() for Hermitian matrix. \n\n";
2: /*
3: Example of usage
4: ./ex36 -check_Hermitian -display_mat -display_vec
5: mpiexec -n 2 ./ex36
7: This example is modified from src/mat/examples/tests/ex127.c
8: */
10: #include <petscksp.h>
14: PetscInt main(PetscInt argc,char **args)
15: {
16: Mat A,As;
17: PetscBool flg,disp_mat=PETSC_FALSE;
19: PetscMPIInt size,rank;
20: PetscInt i,j;
21: PetscScalar v,sigma2;
22: PetscRandom rctx;
23: PetscReal h2,sigma1=100.0;
24: PetscInt dim,Ii,J,n = 3,use_random,rstart,rend;
25: KSP ksp;
26: PC pc;
27: Mat F;
28: PetscInt nneg, nzero, npos;
29:
30: PetscInitialize(&argc,&args,(char *)0,help);
31: #if !defined(PETSC_USE_COMPLEX)
32: SETERRQ(PETSC_COMM_WORLD,1,"This example requires complex numbers");
33: #endif
34: MPI_Comm_size(PETSC_COMM_WORLD,&size);
35: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
36: PetscOptionsHasName(PETSC_NULL, "-display_mat", &disp_mat);
38: PetscOptionsGetReal(PETSC_NULL,"-sigma1",&sigma1,PETSC_NULL);
39: PetscOptionsGetInt(PETSC_NULL,"-n",&n,PETSC_NULL);
40: dim = n*n;
42: MatCreate(PETSC_COMM_WORLD,&A);
43: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,dim,dim);
44: MatSetType(A,MATAIJ);
45: MatSetFromOptions(A);
47: PetscOptionsHasName(PETSC_NULL,"-norandom",&flg);
48: if (flg) use_random = 0;
49: else use_random = 1;
50: if (use_random) {
51: PetscRandomCreate(PETSC_COMM_WORLD,&rctx);
52: PetscRandomSetFromOptions(rctx);
53: PetscRandomSetInterval(rctx,0.0,PETSC_i);
54: PetscRandomGetValue(rctx,&sigma2); /* RealPart(sigma2) == 0.0 */
55: } else {
56: sigma2 = 10.0*PETSC_i;
57: }
58: h2 = 1.0/((n+1)*(n+1));
60: MatGetOwnershipRange(A,&rstart,&rend);
61: for (Ii=rstart; Ii<rend; Ii++) {
62: v = -1.0; i = Ii/n; j = Ii - i*n;
63: if (i>0) {
64: J = Ii-n; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);}
65: if (i<n-1) {
66: J = Ii+n; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);}
67: if (j>0) {
68: J = Ii-1; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);}
69: if (j<n-1) {
70: J = Ii+1; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);}
71: v = 4.0 - sigma1*h2;
72: MatSetValues(A,1,&Ii,1,&Ii,&v,ADD_VALUES);
73: }
74: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
75: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
77: /* Check whether A is symmetric */
78: PetscOptionsHasName(PETSC_NULL, "-check_symmetric", &flg);
79: if (flg) {
80: Mat Trans;
81: MatTranspose(A,MAT_INITIAL_MATRIX, &Trans);
82: MatEqual(A, Trans, &flg);
83: if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"A is not symmetric");
84: MatDestroy(&Trans);
85: }
86: MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);
88: /* make A complex Hermitian */
89: Ii = 0; J = dim-1;
90: if (Ii >= rstart && Ii < rend){
91: v = sigma2*h2; /* RealPart(v) = 0.0 */
92: MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);
93: v = -sigma2*h2;
94: MatSetValues(A,1,&J,1,&Ii,&v,ADD_VALUES);
95: }
97: Ii = dim-2; J = dim-1;
98: if (Ii >= rstart && Ii < rend){
99: v = sigma2*h2; /* RealPart(v) = 0.0 */
100: MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);
101: v = -sigma2*h2;
102: MatSetValues(A,1,&J,1,&Ii,&v,ADD_VALUES);
103: }
105: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
106: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
108: /* Check whether A is Hermitian */
109: PetscOptionsHasName(PETSC_NULL, "-check_Hermitian", &flg);
110: if (flg) {
111: Mat Hermit;
112: if (disp_mat){
113: if (!rank) printf(" A:\n");
114: MatView(A,PETSC_VIEWER_STDOUT_WORLD);
115: }
116: MatHermitianTranspose(A,MAT_INITIAL_MATRIX, &Hermit);
117: if (disp_mat){
118: if (!rank) printf(" A_Hermitian:\n");
119: MatView(Hermit,PETSC_VIEWER_STDOUT_WORLD);
120: }
121: MatEqual(A, Hermit, &flg);
122: if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"A is not Hermitian");
123: MatDestroy(&Hermit);
124: }
125: MatSetOption(A,MAT_HERMITIAN,PETSC_TRUE);
126:
127: /* Create a Hermitian matrix As in sbaij format */
128: MatConvert(A,MATSBAIJ,MAT_INITIAL_MATRIX,&As);
129: if (disp_mat){
130: if (!rank) {PetscPrintf(PETSC_COMM_SELF," As:\n");}
131: MatView(As,PETSC_VIEWER_STDOUT_WORLD);
132: }
134: /* Test MatGetInertia() */
135: KSPCreate(PETSC_COMM_WORLD,&ksp);
136: KSPSetType(ksp,KSPPREONLY);
137: KSPSetOperators(ksp,As,As,DIFFERENT_NONZERO_PATTERN);
139: KSPGetPC(ksp,&pc);
140: PCSetType(pc,PCCHOLESKY);
141: PCSetFromOptions(pc);
143: PCSetUp(pc);
144: PCFactorGetMatrix(pc,&F);
145: MatGetInertia(F,&nneg,&nzero,&npos);
146: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
147: if (!rank){
148: PetscPrintf(PETSC_COMM_SELF," MatInertia: nneg: %D, nzero: %D, npos: %D\n",nneg,nzero,npos);
149: }
151: /* Free spaces */
152: KSPDestroy(&ksp);
153: if (use_random) {PetscRandomDestroy(&rctx);}
154: MatDestroy(&A);
155: MatDestroy(&As);
156: PetscFinalize();
157: return 0;
158: }