Actual source code: ex127.c

  1: static char help[] = "Test MatMult() for Hermitian matrix.\n\n";

  3: #include <petscmat.h>

  5: int main(int argc,char **args)
  6: {
  7:   Mat            A,As;
  8:   PetscBool      flg;
  9:   PetscMPIInt    size;
 10:   PetscInt       i,j;
 11:   PetscScalar    v,sigma2;
 12:   PetscReal      h2,sigma1=100.0;
 13:   PetscInt       dim,Ii,J,n = 3,rstart,rend;

 15:   PetscInitialize(&argc,&args,(char*)0,help);
 16:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 17:   PetscOptionsGetReal(NULL,NULL,"-sigma1",&sigma1,NULL);
 18:   PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
 19:   dim  = n*n;

 21:   MatCreate(PETSC_COMM_WORLD,&A);
 22:   MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,dim,dim);
 23:   MatSetType(A,MATAIJ);
 24:   MatSetFromOptions(A);
 25:   MatSetUp(A);

 27:   sigma2 = 10.0*PETSC_i;
 28:   h2 = 1.0/((n+1)*(n+1));

 30:   MatGetOwnershipRange(A,&rstart,&rend);
 31:   for (Ii=rstart; Ii<rend; Ii++) {
 32:     v = -1.0; i = Ii/n; j = Ii - i*n;
 33:     if (i>0) {
 34:       J = Ii-n; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);
 35:     }
 36:     if (i<n-1) {
 37:       J = Ii+n; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);
 38:     }
 39:     if (j>0) {
 40:       J = Ii-1; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);
 41:     }
 42:     if (j<n-1) {
 43:       J = Ii+1; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);
 44:     }
 45:     v    = 4.0 - sigma1*h2;
 46:     MatSetValues(A,1,&Ii,1,&Ii,&v,ADD_VALUES);
 47:   }
 48:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 49:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

 51:   /* Check whether A is symmetric */
 52:   PetscOptionsHasName(NULL,NULL, "-check_symmetric", &flg);
 53:   if (flg) {
 54:     MatIsSymmetric(A,0.0,&flg);
 56:   }
 57:   MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);

 59:   /* make A complex Hermitian */
 60:   Ii = 0; J = dim-1;
 61:   if (Ii >= rstart && Ii < rend) {
 62:     v    = sigma2*h2; /* RealPart(v) = 0.0 */
 63:     MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);
 64:     v    = -sigma2*h2;
 65:     MatSetValues(A,1,&J,1,&Ii,&v,ADD_VALUES);
 66:   }

 68:   Ii = dim-2; J = dim-1;
 69:   if (Ii >= rstart && Ii < rend) {
 70:     v    = sigma2*h2; /* RealPart(v) = 0.0 */
 71:     MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);
 72:     v    = -sigma2*h2;
 73:     MatSetValues(A,1,&J,1,&Ii,&v,ADD_VALUES);
 74:   }

 76:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 77:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
 78:   MatViewFromOptions(A,NULL,"-disp_mat");

 80:   /* Check whether A is Hermitian, then set A->hermitian flag */
 81:   PetscOptionsHasName(NULL,NULL, "-check_Hermitian", &flg);
 82:   if (flg && size == 1) {
 83:     MatIsHermitian(A,0.0,&flg);
 85:   }
 86:   MatSetOption(A,MAT_HERMITIAN,PETSC_TRUE);

 88: #if defined(PETSC_HAVE_SUPERLU_DIST)
 89:   /* Test Cholesky factorization */
 90:   PetscOptionsHasName(NULL,NULL, "-test_choleskyfactor", &flg);
 91:   if (flg) {
 92:     Mat      F;
 93:     IS       perm,iperm;
 94:     MatFactorInfo info;
 95:     PetscInt nneg,nzero,npos;

 97:     MatGetFactor(A,MATSOLVERSUPERLU_DIST,MAT_FACTOR_CHOLESKY,&F);
 98:     MatGetOrdering(A,MATORDERINGND,&perm,&iperm);
 99:     MatCholeskyFactorSymbolic(F,A,perm,&info);
100:     MatCholeskyFactorNumeric(F,A,&info);

102:     MatGetInertia(F,&nneg,&nzero,&npos);
103:     PetscPrintf(PETSC_COMM_WORLD," MatInertia: nneg: %" PetscInt_FMT ", nzero: %" PetscInt_FMT ", npos: %" PetscInt_FMT "\n",nneg,nzero,npos);
104:     MatDestroy(&F);
105:     ISDestroy(&perm);
106:     ISDestroy(&iperm);
107:   }
108: #endif

110:   /* Create a Hermitian matrix As in sbaij format */
111:   MatConvert(A,MATSBAIJ,MAT_INITIAL_MATRIX,&As);
112:   MatViewFromOptions(As,NULL,"-disp_mat");

114:   /* Test MatMult */
115:   MatMultEqual(A,As,10,&flg);
117:   MatMultAddEqual(A,As,10,&flg);

120:   /* Free spaces */
121:   MatDestroy(&A);
122:   MatDestroy(&As);
123:   PetscFinalize();
124:   return 0;
125: }

127: /*TEST

129:    build:
130:       requires: complex

132:    test:
133:       args: -n 1000
134:       output_file: output/ex127.out

136:    test:
137:       suffix: 2
138:       nsize: 3
139:       args: -n 1000
140:       output_file: output/ex127.out

142:    test:
143:       suffix: superlu_dist
144:       nsize: 3
145:       requires: superlu_dist
146:       args: -test_choleskyfactor -mat_superlu_dist_rowperm NOROWPERM
147:       output_file: output/ex127_superlu_dist.out
148: TEST*/