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;

 16:   PetscInitialize(&argc, &args, (char *)0, help);
 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;
 34:     i = Ii / n;
 35:     j = Ii - i * n;
 36:     if (i > 0) {
 37:       J = Ii - n;
 38:       MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES);
 39:     }
 40:     if (i < n - 1) {
 41:       J = Ii + n;
 42:       MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES);
 43:     }
 44:     if (j > 0) {
 45:       J = Ii - 1;
 46:       MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES);
 47:     }
 48:     if (j < n - 1) {
 49:       J = Ii + 1;
 50:       MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES);
 51:     }
 52:     v = 4.0 - sigma1 * h2;
 53:     MatSetValues(A, 1, &Ii, 1, &Ii, &v, ADD_VALUES);
 54:   }
 55:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
 56:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);

 58:   /* Check whether A is symmetric */
 59:   PetscOptionsHasName(NULL, NULL, "-check_symmetric", &flg);
 60:   if (flg) {
 61:     MatIsSymmetric(A, 0.0, &flg);
 63:   }
 64:   MatSetOption(A, MAT_SYMMETRIC, PETSC_TRUE);

 66:   /* make A complex Hermitian */
 67:   Ii = 0;
 68:   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:   Ii = dim - 2;
 77:   J  = dim - 1;
 78:   if (Ii >= rstart && Ii < rend) {
 79:     v = sigma2 * h2; /* RealPart(v) = 0.0 */
 80:     MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES);
 81:     v = -sigma2 * h2;
 82:     MatSetValues(A, 1, &J, 1, &Ii, &v, ADD_VALUES);
 83:   }

 85:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
 86:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
 87:   MatViewFromOptions(A, NULL, "-disp_mat");

 89:   /* Check whether A is Hermitian, then set A->hermitian flag */
 90:   PetscOptionsHasName(NULL, NULL, "-check_Hermitian", &flg);
 91:   if (flg && size == 1) {
 92:     MatIsHermitian(A, 0.0, &flg);
 94:   }
 95:   MatSetOption(A, MAT_HERMITIAN, PETSC_TRUE);

 97: #if defined(PETSC_HAVE_SUPERLU_DIST)
 98:   /* Test Cholesky factorization */
 99:   PetscOptionsHasName(NULL, NULL, "-test_choleskyfactor", &flg);
100:   if (flg) {
101:     Mat           F;
102:     IS            perm, iperm;
103:     MatFactorInfo info;
104:     PetscInt      nneg, nzero, npos;

106:     MatGetFactor(A, MATSOLVERSUPERLU_DIST, MAT_FACTOR_CHOLESKY, &F);
107:     MatGetOrdering(A, MATORDERINGND, &perm, &iperm);
108:     MatCholeskyFactorSymbolic(F, A, perm, &info);
109:     MatCholeskyFactorNumeric(F, A, &info);

111:     MatGetInertia(F, &nneg, &nzero, &npos);
112:     PetscPrintf(PETSC_COMM_WORLD, " MatInertia: nneg: %" PetscInt_FMT ", nzero: %" PetscInt_FMT ", npos: %" PetscInt_FMT "\n", nneg, nzero, npos);
113:     MatDestroy(&F);
114:     ISDestroy(&perm);
115:     ISDestroy(&iperm);
116:   }
117: #endif

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

123:   /* Test MatMult */
124:   MatMultEqual(A, As, 10, &flg);
126:   MatMultAddEqual(A, As, 10, &flg);

129:   /* Free spaces */
130:   MatDestroy(&A);
131:   MatDestroy(&As);
132:   PetscFinalize();
133:   return 0;
134: }

136: /*TEST

138:    build:
139:       requires: complex

141:    test:
142:       args: -n 1000
143:       output_file: output/ex127.out

145:    test:
146:       suffix: 2
147:       nsize: 3
148:       args: -n 1000
149:       output_file: output/ex127.out

151:    test:
152:       suffix: superlu_dist
153:       nsize: 3
154:       requires: superlu_dist
155:       args: -test_choleskyfactor -mat_superlu_dist_rowperm NOROWPERM
156:       output_file: output/ex127_superlu_dist.out
157: TEST*/