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 { 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; 14 15 PetscFunctionBeginUser; 16 PetscCall(PetscInitialize(&argc, &args, (char *)0, 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; 21 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)); 27 28 sigma2 = 10.0 * PETSC_i; 29 h2 = 1.0 / ((n + 1) * (n + 1)); 30 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)); 57 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)); 65 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 } 75 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 } 84 85 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 86 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 87 PetscCall(MatViewFromOptions(A, NULL, "-disp_mat")); 88 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)); 96 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; 105 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)); 110 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 118 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")); 122 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"); 128 129 /* Free spaces */ 130 PetscCall(MatDestroy(&A)); 131 PetscCall(MatDestroy(&As)); 132 PetscCall(PetscFinalize()); 133 return 0; 134 } 135 136 /*TEST 137 138 build: 139 requires: complex 140 141 test: 142 args: -n 1000 143 output_file: output/ex127.out 144 145 test: 146 suffix: 2 147 nsize: 3 148 args: -n 1000 149 output_file: output/ex127.out 150 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*/ 158