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;
 10:   PetscMPIInt    size;
 11:   PetscInt       i,j;
 12:   PetscScalar    v,sigma2;
 13:   PetscReal      h2,sigma1=100.0;
 14:   PetscInt       dim,Ii,J,n = 3,rstart,rend;

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

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

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

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

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

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

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

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

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

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

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

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

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

115:   /* Test MatMult */
116:   MatMultEqual(A,As,10,&flg);
117:   if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"MatMult not equal");
118:   MatMultAddEqual(A,As,10,&flg);
119:   if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"MatMultAdd not equal");

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


129: /*TEST

131:    build:
132:       requires: complex

134:    test:
135:       args: -n 1000
136:       output_file: output/ex127.out

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


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