1 static char help[] = "Test Hypre matrix APIs\n"; 2 3 #include <petscmathypre.h> 4 5 int main(int argc, char **args) 6 { 7 Mat A, B, C; 8 PetscReal err; 9 PetscInt i, j, M = 20; 10 PetscMPIInt NP; 11 MPI_Comm comm; 12 PetscInt *rows; 13 14 PetscFunctionBeginUser; 15 PetscCall(PetscInitialize(&argc, &args, (char *)0, help)); 16 comm = PETSC_COMM_WORLD; 17 PetscCallMPI(MPI_Comm_size(comm, &NP)); 18 PetscCall(PetscOptionsGetInt(NULL, NULL, "-M", &M, NULL)); 19 PetscCheck(M >= 6, PETSC_COMM_WORLD, PETSC_ERR_SUP, "Matrix has to have more than 6 columns"); 20 /* Hypre matrix */ 21 PetscCall(MatCreate(comm, &B)); 22 PetscCall(MatSetSizes(B, PETSC_DECIDE, PETSC_DECIDE, M, M)); 23 PetscCall(MatSetType(B, MATHYPRE)); 24 PetscCall(MatHYPRESetPreallocation(B, 9, NULL, 9, NULL)); 25 26 /* PETSc AIJ matrix */ 27 PetscCall(MatCreate(comm, &A)); 28 PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, M, M)); 29 PetscCall(MatSetType(A, MATAIJ)); 30 PetscCall(MatSeqAIJSetPreallocation(A, 9, NULL)); 31 PetscCall(MatMPIAIJSetPreallocation(A, 9, NULL, 9, NULL)); 32 33 /*Set Values */ 34 for (i = 0; i < M; i++) { 35 PetscInt cols[] = {0, 1, 2, 3, 4, 5}; 36 PetscScalar vals[6] = {0}; 37 PetscScalar value[] = {100}; 38 for (j = 0; j < 6; j++) vals[j] = ((PetscReal)j) / NP; 39 40 PetscCall(MatSetValues(B, 1, &i, 6, cols, vals, ADD_VALUES)); 41 PetscCall(MatSetValues(B, 1, &i, 1, &i, value, ADD_VALUES)); 42 PetscCall(MatSetValues(A, 1, &i, 6, cols, vals, ADD_VALUES)); 43 PetscCall(MatSetValues(A, 1, &i, 1, &i, value, ADD_VALUES)); 44 } 45 46 /* MAT_FLUSH_ASSEMBLY currently not supported */ 47 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 48 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 49 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 50 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 51 52 /* Compare A and B */ 53 PetscCall(MatConvert(B, MATAIJ, MAT_INITIAL_MATRIX, &C)); 54 PetscCall(MatAXPY(C, -1., A, SAME_NONZERO_PATTERN)); 55 PetscCall(MatNorm(C, NORM_INFINITY, &err)); 56 PetscCheck(err <= PETSC_SMALL, PetscObjectComm((PetscObject)B), PETSC_ERR_PLIB, "Error MatSetValues %g", err); 57 PetscCall(MatDestroy(&C)); 58 59 /* MatZeroRows */ 60 PetscCall(PetscMalloc1(M, &rows)); 61 for (i = 0; i < M; i++) rows[i] = i; 62 PetscCall(MatZeroRows(B, M, rows, 10.0, NULL, NULL)); 63 PetscCall(MatSetOption(A, MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE)); 64 PetscCall(MatZeroRows(A, M, rows, 10.0, NULL, NULL)); 65 PetscCall(MatConvert(B, MATAIJ, MAT_INITIAL_MATRIX, &C)); 66 PetscCall(MatAXPY(C, -1., A, SAME_NONZERO_PATTERN)); 67 PetscCall(MatNorm(C, NORM_INFINITY, &err)); 68 PetscCheck(err <= PETSC_SMALL, PetscObjectComm((PetscObject)B), PETSC_ERR_PLIB, "Error MatZeroRows %g", err); 69 PetscCall(MatDestroy(&C)); 70 PetscCall(PetscFree(rows)); 71 72 /* Test MatZeroEntries */ 73 PetscCall(MatZeroEntries(B)); 74 PetscCall(MatConvert(B, MATAIJ, MAT_INITIAL_MATRIX, &C)); 75 PetscCall(MatNorm(C, NORM_INFINITY, &err)); 76 PetscCheck(err <= PETSC_SMALL, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Error MatZeroEntries %g", err); 77 PetscCall(MatDestroy(&C)); 78 79 /* Insert Values */ 80 for (i = 0; i < M; i++) { 81 PetscInt cols[] = {0, 1, 2, 3, 4, 5}; 82 PetscScalar vals[6] = {0}; 83 PetscScalar value[] = {100}; 84 85 for (j = 0; j < 6; j++) vals[j] = ((PetscReal)j) / NP; 86 87 PetscCall(MatSetValues(B, 1, &i, 6, cols, vals, INSERT_VALUES)); 88 PetscCall(MatSetValues(B, 1, &i, 1, &i, value, INSERT_VALUES)); 89 PetscCall(MatSetValues(A, 1, &i, 6, cols, vals, INSERT_VALUES)); 90 PetscCall(MatSetValues(A, 1, &i, 1, &i, value, INSERT_VALUES)); 91 } 92 93 /* MAT_FLUSH_ASSEMBLY currently not supported */ 94 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 95 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 96 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 97 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 98 99 /* Rows are not sorted with HYPRE so we need an intermediate sort 100 They use a temporary buffer, so we can sort inplace the const memory */ 101 { 102 const PetscInt *idxA, *idxB; 103 const PetscScalar *vA, *vB; 104 PetscInt rstart, rend, nzA, nzB; 105 PetscInt cols[] = {0, 1, 2, 3, 4, -5}; 106 PetscInt *rows; 107 PetscScalar *valuesA, *valuesB; 108 PetscBool flg; 109 110 PetscCall(MatGetOwnershipRange(A, &rstart, &rend)); 111 for (i = rstart; i < rend; i++) { 112 PetscCall(MatGetRow(A, i, &nzA, &idxA, &vA)); 113 PetscCall(MatGetRow(B, i, &nzB, &idxB, &vB)); 114 PetscCheck(nzA == nzB, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Error MatGetRow %" PetscInt_FMT, nzA - nzB); 115 PetscCall(PetscSortIntWithScalarArray(nzB, (PetscInt *)idxB, (PetscScalar *)vB)); 116 PetscCall(PetscArraycmp(idxA, idxB, nzA, &flg)); 117 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Error MatGetRow %" PetscInt_FMT " (indices)", i); 118 PetscCall(PetscArraycmp(vA, vB, nzA, &flg)); 119 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Error MatGetRow %" PetscInt_FMT " (values)", i); 120 PetscCall(MatRestoreRow(A, i, &nzA, &idxA, &vA)); 121 PetscCall(MatRestoreRow(B, i, &nzB, &idxB, &vB)); 122 } 123 124 PetscCall(MatGetOwnershipRange(A, &rstart, &rend)); 125 PetscCall(PetscCalloc3((rend - rstart) * 6, &valuesA, (rend - rstart) * 6, &valuesB, rend - rstart, &rows)); 126 for (i = rstart; i < rend; i++) rows[i - rstart] = i; 127 128 PetscCall(MatGetValues(A, rend - rstart, rows, 6, cols, valuesA)); 129 PetscCall(MatGetValues(B, rend - rstart, rows, 6, cols, valuesB)); 130 131 for (i = 0; i < (rend - rstart); i++) { 132 PetscCall(PetscArraycmp(valuesA + 6 * i, valuesB + 6 * i, 6, &flg)); 133 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Error MatGetValues %" PetscInt_FMT, i + rstart); 134 } 135 PetscCall(PetscFree3(valuesA, valuesB, rows)); 136 } 137 138 /* Compare A and B */ 139 PetscCall(MatConvert(B, MATAIJ, MAT_INITIAL_MATRIX, &C)); 140 PetscCall(MatAXPY(C, -1., A, SAME_NONZERO_PATTERN)); 141 PetscCall(MatNorm(C, NORM_INFINITY, &err)); 142 PetscCheck(err <= PETSC_SMALL, PetscObjectComm((PetscObject)B), PETSC_ERR_PLIB, "Error MatSetValues with INSERT_VALUES %g", err); 143 144 PetscCall(MatDestroy(&A)); 145 PetscCall(MatDestroy(&B)); 146 PetscCall(MatDestroy(&C)); 147 148 PetscCall(PetscFinalize()); 149 return 0; 150 } 151 152 /*TEST 153 154 build: 155 requires: hypre 156 157 test: 158 suffix: 1 159 requires: !defined(PETSC_HAVE_HYPRE_DEVICE) 160 161 test: 162 suffix: 2 163 requires: !defined(PETSC_HAVE_HYPRE_DEVICE) 164 output_file: output/ex225_1.out 165 nsize: 2 166 167 TEST*/ 168