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:   PetscFunctionBeginUser;
 16:   PetscCall(PetscInitialize(&argc, &args, NULL, help));
 17:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
 18:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-sigma1", &sigma1, NULL));
 19:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
 20:   dim = n * n;

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

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

 31:   PetscCall(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:       PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES));
 39:     }
 40:     if (i < n - 1) {
 41:       J = Ii + n;
 42:       PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES));
 43:     }
 44:     if (j > 0) {
 45:       J = Ii - 1;
 46:       PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES));
 47:     }
 48:     if (j < n - 1) {
 49:       J = Ii + 1;
 50:       PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES));
 51:     }
 52:     v = 4.0 - sigma1 * h2;
 53:     PetscCall(MatSetValues(A, 1, &Ii, 1, &Ii, &v, ADD_VALUES));
 54:   }
 55:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
 56:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));

 58:   /* Check whether A is symmetric */
 59:   PetscCall(PetscOptionsHasName(NULL, NULL, "-check_symmetric", &flg));
 60:   if (flg) {
 61:     PetscCall(MatIsSymmetric(A, 0.0, &flg));
 62:     PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_USER, "A is not symmetric");
 63:   }
 64:   PetscCall(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:     PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES));
 72:     v = -sigma2 * h2;
 73:     PetscCall(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:     PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES));
 81:     v = -sigma2 * h2;
 82:     PetscCall(MatSetValues(A, 1, &J, 1, &Ii, &v, ADD_VALUES));
 83:   }

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

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

 97: #if defined(PETSC_HAVE_SUPERLU_DIST)
 98:   /* Test Cholesky factorization */
 99:   PetscCall(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:     PetscCall(MatGetFactor(A, MATSOLVERSUPERLU_DIST, MAT_FACTOR_CHOLESKY, &F));
107:     PetscCall(MatGetOrdering(A, MATORDERINGND, &perm, &iperm));
108:     PetscCall(MatCholeskyFactorSymbolic(F, A, perm, &info));
109:     PetscCall(MatCholeskyFactorNumeric(F, A, &info));

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

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

123:   /* Test MatMult */
124:   PetscCall(MatMultEqual(A, As, 10, &flg));
125:   PetscCheck(flg, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "MatMult not equal");
126:   PetscCall(MatMultAddEqual(A, As, 10, &flg));
127:   PetscCheck(flg, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "MatMultAdd not equal");

129:   /* Free spaces */
130:   PetscCall(MatDestroy(&A));
131:   PetscCall(MatDestroy(&As));
132:   PetscCall(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*/