Actual source code: ex225.c

  1: static char help[] = "Test Hypre matrix APIs\n";

  3: #include <petscmathypre.h>

  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;

 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));

 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));

 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;

 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:   }

 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));

 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));

 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));

 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));

 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};

 85:     for (j = 0; j < 6; j++) vals[j] = ((PetscReal)j) / NP;

 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:   }

 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));

 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;

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:     }

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;

128:     PetscCall(MatGetValues(A, rend - rstart, rows, 6, cols, valuesA));
129:     PetscCall(MatGetValues(B, rend - rstart, rows, 6, cols, valuesB));

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:   }

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);

144:   PetscCall(MatDestroy(&A));
145:   PetscCall(MatDestroy(&B));
146:   PetscCall(MatDestroy(&C));

148:   PetscCall(PetscFinalize());
149:   return 0;
150: }

152: /*TEST

154:    build:
155:       requires: hypre

157:    test:
158:       suffix: 1
159:       requires: !defined(PETSC_HAVE_HYPRE_DEVICE)

161:    test:
162:       suffix: 2
163:       requires: !defined(PETSC_HAVE_HYPRE_DEVICE)
164:       output_file: output/ex225_1.out
165:       nsize: 2

167: TEST*/