Actual source code: ex132.c

  1: static char help[] = "Test MatAXPY()\n\n";

  3: #include <petscmat.h>

  5: int main(int argc, char **args)
  6: {
  7:   Mat         C, C1, C2, CU;
  8:   PetscScalar v;
  9:   PetscInt    Ii, J, Istart, Iend;
 10:   PetscInt    i, j, m = 3, n;
 11:   PetscMPIInt size;
 12:   PetscBool   mat_nonsymmetric = PETSC_FALSE, flg;
 13:   MatInfo     info;

 15:   PetscFunctionBeginUser;
 16:   PetscCall(PetscInitialize(&argc, &args, NULL, help));
 17:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL));
 18:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
 19:   n = 2 * size;

 21:   /* Set flag if we are doing a nonsymmetric problem; the default is symmetric. */
 22:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-mat_nonsym", &mat_nonsymmetric, NULL));

 24:   PetscCall(MatCreate(PETSC_COMM_WORLD, &C));
 25:   PetscCall(MatSetSizes(C, PETSC_DECIDE, PETSC_DECIDE, m * n, m * n));
 26:   PetscCall(MatSetFromOptions(C));
 27:   PetscCall(MatSeqAIJSetPreallocation(C, 5, NULL));
 28:   PetscCall(MatMPIAIJSetPreallocation(C, 5, NULL, 5, NULL));

 30:   PetscCall(MatGetOwnershipRange(C, &Istart, &Iend));
 31:   for (Ii = Istart; Ii < Iend; Ii++) {
 32:     v = -1.0;
 33:     i = Ii / n;
 34:     j = Ii - i * n;
 35:     if (i > 0) {
 36:       J = Ii - n;
 37:       PetscCall(MatSetValues(C, 1, &Ii, 1, &J, &v, ADD_VALUES));
 38:     }
 39:     if (i < m - 1) {
 40:       J = Ii + n;
 41:       PetscCall(MatSetValues(C, 1, &Ii, 1, &J, &v, ADD_VALUES));
 42:     }
 43:     if (j > 0) {
 44:       J = Ii - 1;
 45:       PetscCall(MatSetValues(C, 1, &Ii, 1, &J, &v, ADD_VALUES));
 46:     }
 47:     if (j < n - 1) {
 48:       J = Ii + 1;
 49:       PetscCall(MatSetValues(C, 1, &Ii, 1, &J, &v, ADD_VALUES));
 50:     }
 51:     v = 4.0;
 52:     PetscCall(MatSetValues(C, 1, &Ii, 1, &Ii, &v, ADD_VALUES));
 53:   }

 55:   /* Make the matrix nonsymmetric if desired */
 56:   if (mat_nonsymmetric) {
 57:     for (Ii = Istart; Ii < Iend; Ii++) {
 58:       v = -1.5;
 59:       i = Ii / n;
 60:       if (i > 1) {
 61:         J = Ii - n - 1;
 62:         PetscCall(MatSetValues(C, 1, &Ii, 1, &J, &v, ADD_VALUES));
 63:       }
 64:     }
 65:   } else {
 66:     PetscCall(MatSetOption(C, MAT_SYMMETRIC, PETSC_TRUE));
 67:     PetscCall(MatSetOption(C, MAT_SYMMETRY_ETERNAL, PETSC_TRUE));
 68:   }
 69:   PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
 70:   PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
 71:   PetscCall(PetscObjectSetName((PetscObject)C, "C"));
 72:   PetscCall(MatViewFromOptions(C, NULL, "-view"));

 74:   /* C1 = 2.0*C1 + C, C1 is anti-diagonal and has different non-zeros than C */
 75:   PetscCall(MatCreate(PETSC_COMM_WORLD, &C1));
 76:   PetscCall(MatSetSizes(C1, PETSC_DECIDE, PETSC_DECIDE, m * n, m * n));
 77:   PetscCall(MatSetFromOptions(C1));
 78:   PetscCall(MatSeqAIJSetPreallocation(C1, 1, NULL));
 79:   PetscCall(MatMPIAIJSetPreallocation(C1, 1, NULL, 1, NULL));
 80:   for (Ii = Istart; Ii < Iend; Ii++) {
 81:     v = 1.0;
 82:     i = m * n - Ii - 1;
 83:     j = Ii;
 84:     PetscCall(MatSetValues(C1, 1, &i, 1, &j, &v, ADD_VALUES));
 85:   }
 86:   PetscCall(MatAssemblyBegin(C1, MAT_FINAL_ASSEMBLY));
 87:   PetscCall(MatAssemblyEnd(C1, MAT_FINAL_ASSEMBLY));
 88:   PetscCall(PetscObjectSetName((PetscObject)C1, "C1"));
 89:   PetscCall(MatViewFromOptions(C1, NULL, "-view"));
 90:   PetscCall(MatDuplicate(C1, MAT_COPY_VALUES, &CU));

 92:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, " MatAXPY(C1,2.0,C,DIFFERENT_NONZERO_PATTERN)...\n"));
 93:   PetscCall(MatAXPY(C1, 2.0, C, DIFFERENT_NONZERO_PATTERN));
 94:   PetscCall(MatAXPY(CU, 2.0, C, UNKNOWN_NONZERO_PATTERN));
 95:   PetscCall(MatGetInfo(C1, MAT_GLOBAL_SUM, &info));
 96:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, " C1: nz_allocated = %g; nz_used = %g; nz_unneeded = %g\n", info.nz_allocated, info.nz_used, info.nz_unneeded));
 97:   PetscCall(MatViewFromOptions(C1, NULL, "-view"));
 98:   PetscCall(MatMultEqual(CU, C1, 10, &flg));
 99:   if (!flg) {
100:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error UNKNOWN_NONZERO_PATTERN (supposedly DIFFERENT_NONZERO_PATTERN)\n"));
101:     PetscCall(MatViewFromOptions(CU, NULL, "-view"));
102:   }
103:   PetscCall(MatDestroy(&CU));

105:   /* Secondly, compute C1 = 2.0*C2 + C1, C2 has non-zero pattern of C */
106:   PetscCall(MatDuplicate(C, MAT_DO_NOT_COPY_VALUES, &C2));
107:   PetscCall(MatDuplicate(C1, MAT_COPY_VALUES, &CU));

109:   for (Ii = Istart; Ii < Iend; Ii++) {
110:     v = 1.0;
111:     PetscCall(MatSetValues(C2, 1, &Ii, 1, &Ii, &v, ADD_VALUES));
112:   }
113:   PetscCall(MatAssemblyBegin(C2, MAT_FINAL_ASSEMBLY));
114:   PetscCall(MatAssemblyEnd(C2, MAT_FINAL_ASSEMBLY));
115:   PetscCall(PetscObjectSetName((PetscObject)C2, "C2"));
116:   PetscCall(MatViewFromOptions(C2, NULL, "-view"));
117:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, " MatAXPY(C1,2.0,C2,SUBSET_NONZERO_PATTERN)...\n"));
118:   PetscCall(MatAXPY(C1, 2.0, C2, SUBSET_NONZERO_PATTERN));
119:   PetscCall(MatAXPY(CU, 2.0, C2, UNKNOWN_NONZERO_PATTERN));
120:   PetscCall(MatGetInfo(C1, MAT_GLOBAL_SUM, &info));
121:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, " C1: nz_allocated = %g; nz_used = %g; nz_unneeded = %g\n", info.nz_allocated, info.nz_used, info.nz_unneeded));
122:   PetscCall(MatViewFromOptions(C1, NULL, "-view"));
123:   PetscCall(MatMultEqual(CU, C1, 10, &flg));
124:   if (!flg) {
125:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error UNKNOWN_NONZERO_PATTERN (supposedly SUBSET_NONZERO_PATTERN)\n"));
126:     PetscCall(MatViewFromOptions(CU, NULL, "-view"));
127:   }
128:   PetscCall(MatDestroy(&CU));

130:   /* Test SAME_NONZERO_PATTERN computing C2 = C2 + 2.0 * C */
131:   PetscCall(MatDuplicate(C2, MAT_COPY_VALUES, &CU));
132:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, " MatAXPY(C2,2.0,C,SAME_NONZERO_PATTERN)...\n"));
133:   PetscCall(MatAXPY(C2, 2.0, C, SAME_NONZERO_PATTERN));
134:   PetscCall(MatAXPY(CU, 2.0, C, UNKNOWN_NONZERO_PATTERN));
135:   PetscCall(MatGetInfo(C2, MAT_GLOBAL_SUM, &info));
136:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, " C2: nz_allocated = %g; nz_used = %g; nz_unneeded = %g\n", info.nz_allocated, info.nz_used, info.nz_unneeded));
137:   PetscCall(MatViewFromOptions(C2, NULL, "-view"));
138:   PetscCall(MatMultEqual(CU, C2, 10, &flg));
139:   if (!flg) {
140:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error UNKNOWN_NONZERO_PATTERN (supposedly SUBSET_NONZERO_PATTERN)\n"));
141:     PetscCall(MatViewFromOptions(CU, NULL, "-view"));
142:   }
143:   PetscCall(MatDestroy(&CU));

145:   PetscCall(MatDestroy(&C1));
146:   PetscCall(MatDestroy(&C2));
147:   PetscCall(MatDestroy(&C));

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

153: /*TEST

155:    test:
156:      suffix: 1
157:      filter: grep -v " type:" | grep -v "Mat Object"
158:      args: -view
159:      diff_args: -j

161:    test:
162:      output_file: output/ex132_1.out
163:      requires: cuda
164:      suffix: 1_cuda
165:      filter: grep -v " type:" | grep -v "Mat Object"
166:      args: -view -mat_type aijcusparse
167:      diff_args: -j

169:    test:
170:      output_file: output/ex132_1.out
171:      requires: kokkos_kernels
172:      suffix: 1_kokkos
173:      filter: grep -v " type:" | grep -v "Mat Object"
174:      args: -view -mat_type aijkokkos
175:      diff_args: -j

177:    test:
178:      suffix: 2
179:      filter: grep -v " type:" | grep -v "Mat Object"
180:      args: -view -mat_nonsym
181:      diff_args: -j

183:    test:
184:      output_file: output/ex132_2.out
185:      requires: cuda
186:      suffix: 2_cuda
187:      filter: grep -v " type:" | grep -v "Mat Object"
188:      args: -view -mat_type aijcusparse -mat_nonsym
189:      diff_args: -j

191:    test:
192:      output_file: output/ex132_2.out
193:      requires: kokkos_kernels
194:      suffix: 2_kokkos
195:      filter: grep -v " type:" | grep -v "Mat Object"
196:      args: -view -mat_type aijkokkos -mat_nonsym
197:      diff_args: -j

199:    test:
200:      nsize: 2
201:      suffix: 1_par
202:      filter: grep -v " type:" | grep -v "Mat Object"
203:      args: -view
204:      diff_args: -j

206:    test:
207:      nsize: 2
208:      output_file: output/ex132_1_par.out
209:      requires: cuda
210:      suffix: 1_par_cuda
211:      filter: grep -v " type:" | grep -v "Mat Object"
212:      args: -view -mat_type aijcusparse
213:      diff_args: -j

215:    test:
216:      nsize: 2
217:      output_file: output/ex132_1_par.out
218:      requires: kokkos_kernels
219:      suffix: 1_par_kokkos
220:      filter: grep -v " type:" | grep -v "Mat Object"
221:      args: -view -mat_type aijkokkos
222:      diff_args: -j

224:    test:
225:      nsize: 2
226:      suffix: 2_par
227:      filter: grep -v " type:" | grep -v "Mat Object"
228:      args: -view -mat_nonsym
229:      diff_args: -j

231:    test:
232:      nsize: 2
233:      output_file: output/ex132_2_par.out
234:      requires: cuda
235:      suffix: 2_par_cuda
236:      filter: grep -v " type:" | grep -v "Mat Object"
237:      args: -view -mat_type aijcusparse -mat_nonsym
238:      diff_args: -j

240:    testset:
241:      nsize: 2
242:      output_file: output/ex132_2_par.out
243:      requires: kokkos_kernels
244:      filter: grep -v " type:" | grep -v "Mat Object"
245:      args: -view -mat_type aijkokkos -mat_nonsym
246:      diff_args: -j
247:      test:
248:        suffix: 2_par_kokkos_no_gpu_aware
249:        args: -use_gpu_aware_mpi 0
250:      test:
251:        requires: defined(HAVE_MPI_GPU_AWARE)
252:        suffix: 2_par_kokkos_gpu_aware
253:        args: -use_gpu_aware_mpi 1

255: TEST*/