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:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 25:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
 26:   PetscOptionsHasName(NULL,NULL, "-display_mat", &disp_mat);
 27:   PetscOptionsHasName(NULL,NULL, "-display_vec", &disp_vec);

 29:   PetscOptionsGetReal(NULL,NULL,"-sigma1",&sigma1,NULL);
 30:   PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
 31:   dim  = n*n;

 33:   MatCreate(PETSC_COMM_WORLD,&A);
 34:   MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,dim,dim);
 35:   MatSetType(A,MATAIJ);
 36:   MatSetFromOptions(A);
 37:   MatSetUp(A);

 39:   PetscOptionsHasName(NULL,NULL,"-norandom",&flg);
 40:   if (!flg) {
 41:     PetscRandomCreate(PETSC_COMM_WORLD,&rctx);
 42:     PetscRandomSetFromOptions(rctx);
 43:     PetscRandomSetInterval(rctx,0.0,PETSC_i);
 44:     PetscRandomGetValue(rctx,&sigma2); /* RealPart(sigma2) == 0.0 */
 45:   } else {
 46:     sigma2 = 10.0*PETSC_i;
 47:   }
 48:   h2 = 1.0/((n+1)*(n+1));

 50:   MatGetOwnershipRange(A,&rstart,&rend);
 51:   for (Ii=rstart; Ii<rend; Ii++) {
 52:     v = -1.0; i = Ii/n; j = Ii - i*n;
 53:     if (i>0) {
 54:       J = Ii-n; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);
 55:     }
 56:     if (i<n-1) {
 57:       J = Ii+n; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);
 58:     }
 59:     if (j>0) {
 60:       J = Ii-1; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);
 61:     }
 62:     if (j<n-1) {
 63:       J = Ii+1; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);
 64:     }
 65:     v    = 4.0 - sigma1*h2;
 66:     MatSetValues(A,1,&Ii,1,&Ii,&v,ADD_VALUES);
 67:   }
 68:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 69:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

 71:   /* Check whether A is symmetric */
 72:   PetscOptionsHasName(NULL,NULL, "-check_symmetric", &flg);
 73:   if (flg) {
 74:     Mat Trans;
 75:     MatTranspose(A,MAT_INITIAL_MATRIX, &Trans);
 76:     MatEqual(A, Trans, &flg);
 77:     if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"A is not symmetric");
 78:     MatDestroy(&Trans);
 79:   }
 80:   MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);

 82:   /* make A complex Hermitian */
 83:   Ii = 0; J = dim-1;
 84:   if (Ii >= rstart && Ii < rend) {
 85:     v    = sigma2*h2; /* RealPart(v) = 0.0 */
 86:     MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);
 87:     v    = -sigma2*h2;
 88:     MatSetValues(A,1,&J,1,&Ii,&v,ADD_VALUES);
 89:   }

 91:   Ii = dim-2; J = dim-1;
 92:   if (Ii >= rstart && Ii < rend) {
 93:     v    = sigma2*h2; /* RealPart(v) = 0.0 */
 94:     MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);
 95:     v    = -sigma2*h2;
 96:     MatSetValues(A,1,&J,1,&Ii,&v,ADD_VALUES);
 97:   }

 99:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
100:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

102:   /* Check whether A is Hermitian */
103:   PetscOptionsHasName(NULL,NULL, "-check_Hermitian", &flg);
104:   if (flg) {
105:     Mat Hermit;
106:     if (disp_mat) {
107:       PetscPrintf(PETSC_COMM_WORLD," A:\n");
108:       MatView(A,PETSC_VIEWER_STDOUT_WORLD);
109:     }
110:     MatHermitianTranspose(A,MAT_INITIAL_MATRIX, &Hermit);
111:     if (disp_mat) {
112:       PetscPrintf(PETSC_COMM_WORLD," A_Hermitian:\n");
113:       MatView(Hermit,PETSC_VIEWER_STDOUT_WORLD);
114:     }
115:     MatEqual(A, Hermit, &flg);
116:     if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"A is not Hermitian");
117:     MatDestroy(&Hermit);
118:   }
119:   MatSetOption(A,MAT_HERMITIAN,PETSC_TRUE);

121:   /* Create a Hermitian matrix As in sbaij format */
122:   MatConvert(A,MATSBAIJ,MAT_INITIAL_MATRIX,&As);
123:   if (disp_mat) {
124:     PetscPrintf(PETSC_COMM_WORLD," As:\n");
125:     MatView(As,PETSC_VIEWER_STDOUT_WORLD);
126:   }

128:   MatGetLocalSize(A,&m,&n);
129:   VecCreate(PETSC_COMM_WORLD,&x);
130:   VecSetSizes(x,n,PETSC_DECIDE);
131:   VecSetFromOptions(x);
132:   if (rctx) {
133:     VecSetRandom(x,rctx);
134:   } else {
135:     VecSet(x,1.0);
136:   }

138:   /* Create vectors */
139:   VecCreate(PETSC_COMM_WORLD,&y);
140:   VecSetSizes(y,m,PETSC_DECIDE);
141:   VecSetFromOptions(y);
142:   VecDuplicate(y,&ys);

144:   /* Test MatMult */
145:   MatMult(A,x,y);
146:   MatMult(As,x,ys);
147:   if (disp_vec) {
148:     PetscPrintf(PETSC_COMM_WORLD,"y = A*x:\n");
149:     VecView(y,PETSC_VIEWER_STDOUT_WORLD);
150:     PetscPrintf(PETSC_COMM_WORLD,"ys = As*x:\n");
151:     VecView(ys,PETSC_VIEWER_STDOUT_WORLD);
152:   }
153:   VecAXPY(y,-1.0,ys);
154:   VecNorm(y,NORM_INFINITY,&norm);
155:   if (norm > 100.0*PETSC_MACHINE_EPSILON || disp_vec) {
156:     PetscPrintf(PETSC_COMM_WORLD,"|| A*x - As*x || = %g\n",(double)norm);
157:   }

159:   /* Free spaces */
160:   PetscRandomDestroy(&rctx);
161:   MatDestroy(&A);
162:   MatDestroy(&As);

164:   VecDestroy(&x);
165:   VecDestroy(&y);
166:   VecDestroy(&ys);
167:   PetscFinalize();
168:   return ierr;
169: }