Actual source code: ex225.c


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

  4: #include <petscmathypre.h>

  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;

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

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

 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;

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

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

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

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

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

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

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

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

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

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;

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

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;

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

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

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

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

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

153: /*TEST

155:    build:
156:       requires: hypre

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

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

168: TEST*/