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