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