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: }