1 static char help[] = "Test MatMult() for Hermitian matrix.\n\n"; 2 3 #include <petscmat.h> 4 5 int main(int argc, char **args) { 6 Mat A, As; 7 PetscBool flg; 8 PetscMPIInt size; 9 PetscInt i, j; 10 PetscScalar v, sigma2; 11 PetscReal h2, sigma1 = 100.0; 12 PetscInt dim, Ii, J, n = 3, rstart, rend; 13 14 PetscFunctionBeginUser; 15 PetscCall(PetscInitialize(&argc, &args, (char *)0, help)); 16 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 17 PetscCall(PetscOptionsGetReal(NULL, NULL, "-sigma1", &sigma1, NULL)); 18 PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL)); 19 dim = n * n; 20 21 PetscCall(MatCreate(PETSC_COMM_WORLD, &A)); 22 PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, dim, dim)); 23 PetscCall(MatSetType(A, MATAIJ)); 24 PetscCall(MatSetFromOptions(A)); 25 PetscCall(MatSetUp(A)); 26 27 sigma2 = 10.0 * PETSC_i; 28 h2 = 1.0 / ((n + 1) * (n + 1)); 29 30 PetscCall(MatGetOwnershipRange(A, &rstart, &rend)); 31 for (Ii = rstart; Ii < rend; Ii++) { 32 v = -1.0; 33 i = Ii / n; 34 j = Ii - i * n; 35 if (i > 0) { 36 J = Ii - n; 37 PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES)); 38 } 39 if (i < n - 1) { 40 J = Ii + n; 41 PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES)); 42 } 43 if (j > 0) { 44 J = Ii - 1; 45 PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES)); 46 } 47 if (j < n - 1) { 48 J = Ii + 1; 49 PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES)); 50 } 51 v = 4.0 - sigma1 * h2; 52 PetscCall(MatSetValues(A, 1, &Ii, 1, &Ii, &v, ADD_VALUES)); 53 } 54 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 55 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 56 57 /* Check whether A is symmetric */ 58 PetscCall(PetscOptionsHasName(NULL, NULL, "-check_symmetric", &flg)); 59 if (flg) { 60 PetscCall(MatIsSymmetric(A, 0.0, &flg)); 61 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_USER, "A is not symmetric"); 62 } 63 PetscCall(MatSetOption(A, MAT_SYMMETRIC, PETSC_TRUE)); 64 65 /* make A complex Hermitian */ 66 Ii = 0; 67 J = dim - 1; 68 if (Ii >= rstart && Ii < rend) { 69 v = sigma2 * h2; /* RealPart(v) = 0.0 */ 70 PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES)); 71 v = -sigma2 * h2; 72 PetscCall(MatSetValues(A, 1, &J, 1, &Ii, &v, ADD_VALUES)); 73 } 74 75 Ii = dim - 2; 76 J = dim - 1; 77 if (Ii >= rstart && Ii < rend) { 78 v = sigma2 * h2; /* RealPart(v) = 0.0 */ 79 PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES)); 80 v = -sigma2 * h2; 81 PetscCall(MatSetValues(A, 1, &J, 1, &Ii, &v, ADD_VALUES)); 82 } 83 84 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 85 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 86 PetscCall(MatViewFromOptions(A, NULL, "-disp_mat")); 87 88 /* Check whether A is Hermitian, then set A->hermitian flag */ 89 PetscCall(PetscOptionsHasName(NULL, NULL, "-check_Hermitian", &flg)); 90 if (flg && size == 1) { 91 PetscCall(MatIsHermitian(A, 0.0, &flg)); 92 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_USER, "A is not Hermitian"); 93 } 94 PetscCall(MatSetOption(A, MAT_HERMITIAN, PETSC_TRUE)); 95 96 #if defined(PETSC_HAVE_SUPERLU_DIST) 97 /* Test Cholesky factorization */ 98 PetscCall(PetscOptionsHasName(NULL, NULL, "-test_choleskyfactor", &flg)); 99 if (flg) { 100 Mat F; 101 IS perm, iperm; 102 MatFactorInfo info; 103 PetscInt nneg, nzero, npos; 104 105 PetscCall(MatGetFactor(A, MATSOLVERSUPERLU_DIST, MAT_FACTOR_CHOLESKY, &F)); 106 PetscCall(MatGetOrdering(A, MATORDERINGND, &perm, &iperm)); 107 PetscCall(MatCholeskyFactorSymbolic(F, A, perm, &info)); 108 PetscCall(MatCholeskyFactorNumeric(F, A, &info)); 109 110 PetscCall(MatGetInertia(F, &nneg, &nzero, &npos)); 111 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " MatInertia: nneg: %" PetscInt_FMT ", nzero: %" PetscInt_FMT ", npos: %" PetscInt_FMT "\n", nneg, nzero, npos)); 112 PetscCall(MatDestroy(&F)); 113 PetscCall(ISDestroy(&perm)); 114 PetscCall(ISDestroy(&iperm)); 115 } 116 #endif 117 118 /* Create a Hermitian matrix As in sbaij format */ 119 PetscCall(MatConvert(A, MATSBAIJ, MAT_INITIAL_MATRIX, &As)); 120 PetscCall(MatViewFromOptions(As, NULL, "-disp_mat")); 121 122 /* Test MatMult */ 123 PetscCall(MatMultEqual(A, As, 10, &flg)); 124 PetscCheck(flg, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "MatMult not equal"); 125 PetscCall(MatMultAddEqual(A, As, 10, &flg)); 126 PetscCheck(flg, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "MatMultAdd not equal"); 127 128 /* Free spaces */ 129 PetscCall(MatDestroy(&A)); 130 PetscCall(MatDestroy(&As)); 131 PetscCall(PetscFinalize()); 132 return 0; 133 } 134 135 /*TEST 136 137 build: 138 requires: complex 139 140 test: 141 args: -n 1000 142 output_file: output/ex127.out 143 144 test: 145 suffix: 2 146 nsize: 3 147 args: -n 1000 148 output_file: output/ex127.out 149 150 test: 151 suffix: superlu_dist 152 nsize: 3 153 requires: superlu_dist 154 args: -test_choleskyfactor -mat_superlu_dist_rowperm NOROWPERM 155 output_file: output/ex127_superlu_dist.out 156 TEST*/ 157