1 static char help[] = "Test MatSetPreallocationCOO and MatSetValuesCOO\n\n"; 2 3 #include <petscmat.h> 4 5 int main(int argc, char **args) 6 { 7 Mat A, At, AAt, T = NULL; 8 Vec x, y, z; 9 ISLocalToGlobalMapping rl2g, cl2g; 10 IS is; 11 PetscLayout rmap, cmap; 12 PetscInt *it, *jt; 13 PetscInt n1 = 11, n2 = 9; 14 PetscInt i1[] = {7, 6, 2, 0, 4, 1, 1, 0, 2, 2, 1, -1, -1}; 15 PetscInt j1[] = {1, 4, 3, 5, 3, 3, 4, 5, 0, 3, 1, -1, -1}; 16 PetscInt i2[] = {7, 6, 2, 0, 4, 1, 1, 2, 1, -1, -1}; 17 PetscInt j2[] = {1, 4, 3, 5, 3, 3, 4, 0, 1, -1, -1}; 18 PetscScalar v1[] = {-1., 1., 2., 3., 4., 5., 6., 7., 8., 9., 10., PETSC_MAX_REAL, PETSC_MAX_REAL}; 19 PetscScalar v2[] = {1., -1., -2., -3., -4., -5., -6., -7., -8., -9., -10., PETSC_MAX_REAL, PETSC_MAX_REAL}; 20 PetscInt N = 6, m = 8, M, rstart, cstart, i; 21 PetscMPIInt size; 22 PetscBool loc = PETSC_FALSE; 23 PetscBool locdiag = PETSC_TRUE; 24 PetscBool localapi = PETSC_FALSE; 25 PetscBool neg = PETSC_FALSE; 26 PetscBool ismatis, ismpiaij, ishypre; 27 28 PetscFunctionBeginUser; 29 PetscCall(PetscInitialize(&argc, &args, (char *)0, help)); 30 PetscCall(PetscOptionsGetBool(NULL, NULL, "-neg", &neg, NULL)); 31 PetscCall(PetscOptionsGetBool(NULL, NULL, "-loc", &loc, NULL)); 32 PetscCall(PetscOptionsGetBool(NULL, NULL, "-locdiag", &locdiag, NULL)); 33 PetscCall(PetscOptionsGetBool(NULL, NULL, "-localapi", &localapi, NULL)); 34 PetscCall(MatCreate(PETSC_COMM_WORLD, &A)); 35 if (loc) { 36 if (locdiag) { 37 PetscCall(MatSetSizes(A, m, N, PETSC_DECIDE, PETSC_DECIDE)); 38 } else { 39 PetscCall(MatSetSizes(A, m, m + N, PETSC_DECIDE, PETSC_DECIDE)); 40 } 41 } else { 42 PetscCall(MatSetSizes(A, m, PETSC_DECIDE, PETSC_DECIDE, N)); 43 } 44 PetscCall(MatSetFromOptions(A)); 45 PetscCall(MatGetLayouts(A, &rmap, &cmap)); 46 PetscCall(PetscLayoutSetUp(rmap)); 47 PetscCall(PetscLayoutSetUp(cmap)); 48 PetscCall(PetscLayoutGetRange(rmap, &rstart, NULL)); 49 PetscCall(PetscLayoutGetRange(cmap, &cstart, NULL)); 50 PetscCall(PetscLayoutGetSize(rmap, &M)); 51 PetscCall(PetscLayoutGetSize(cmap, &N)); 52 53 PetscCall(PetscObjectTypeCompare((PetscObject)A, MATIS, &ismatis)); 54 PetscCall(PetscObjectTypeCompare((PetscObject)A, MATHYPRE, &ishypre)); 55 56 /* create fake l2g maps to test the local API */ 57 PetscCall(ISCreateStride(PETSC_COMM_WORLD, M - rstart, rstart, 1, &is)); 58 PetscCall(ISLocalToGlobalMappingCreateIS(is, &rl2g)); 59 PetscCall(ISDestroy(&is)); 60 PetscCall(ISCreateStride(PETSC_COMM_WORLD, N, 0, 1, &is)); 61 PetscCall(ISLocalToGlobalMappingCreateIS(is, &cl2g)); 62 PetscCall(ISDestroy(&is)); 63 PetscCall(MatSetLocalToGlobalMapping(A, rl2g, cl2g)); 64 PetscCall(ISLocalToGlobalMappingDestroy(&rl2g)); 65 PetscCall(ISLocalToGlobalMappingDestroy(&cl2g)); 66 67 PetscCall(MatCreateVecs(A, &x, &y)); 68 PetscCall(MatCreateVecs(A, NULL, &z)); 69 PetscCall(VecSet(x, 1.)); 70 PetscCall(VecSet(z, 2.)); 71 if (!localapi) 72 for (i = 0; i < n1; i++) i1[i] += rstart; 73 if (!localapi) 74 for (i = 0; i < n2; i++) i2[i] += rstart; 75 if (loc) { 76 if (locdiag) { 77 for (i = 0; i < n1; i++) j1[i] += cstart; 78 for (i = 0; i < n2; i++) j2[i] += cstart; 79 } else { 80 for (i = 0; i < n1; i++) j1[i] += cstart + m; 81 for (i = 0; i < n2; i++) j2[i] += cstart + m; 82 } 83 } 84 if (neg) { 85 n1 += 2; 86 n2 += 2; 87 } 88 /* MatSetPreallocationCOOLocal maps the indices! */ 89 PetscCall(PetscMalloc2(PetscMax(n1, n2), &it, PetscMax(n1, n2), &jt)); 90 /* test with repeated entries */ 91 if (!localapi) { 92 PetscCall(MatSetPreallocationCOO(A, n1, i1, j1)); 93 } else { 94 PetscCall(PetscArraycpy(it, i1, n1)); 95 PetscCall(PetscArraycpy(jt, j1, n1)); 96 PetscCall(MatSetPreallocationCOOLocal(A, n1, it, jt)); 97 } 98 PetscCall(MatSetValuesCOO(A, v1, ADD_VALUES)); 99 PetscCall(MatMult(A, x, y)); 100 PetscCall(MatView(A, NULL)); 101 PetscCall(VecView(y, NULL)); 102 PetscCall(MatSetValuesCOO(A, v2, ADD_VALUES)); 103 PetscCall(MatMultAdd(A, x, y, y)); 104 PetscCall(MatView(A, NULL)); 105 PetscCall(VecView(y, NULL)); 106 T = A; 107 if (ishypre) PetscCall(MatConvert(A, MATAIJ, MAT_INITIAL_MATRIX, &T)); 108 PetscCall(MatTranspose(T, MAT_INITIAL_MATRIX, &At)); 109 if (!ismatis) { 110 PetscCall(MatMatMult(T, At, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &AAt)); 111 PetscCall(MatView(AAt, NULL)); 112 PetscCall(MatDestroy(&AAt)); 113 PetscCall(MatMatMult(At, T, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &AAt)); 114 PetscCall(MatView(AAt, NULL)); 115 PetscCall(MatDestroy(&AAt)); 116 } 117 PetscCall(MatDestroy(&At)); 118 if (ishypre) PetscCall(MatDestroy(&T)); 119 120 /* INSERT_VALUES will overwrite matrix entries but 121 still perform the sum of the repeated entries */ 122 PetscCall(MatSetValuesCOO(A, v2, INSERT_VALUES)); 123 PetscCall(MatView(A, NULL)); 124 125 /* test with unique entries */ 126 PetscCall(PetscArraycpy(it, i2, n2)); 127 PetscCall(PetscArraycpy(jt, j2, n2)); 128 if (!localapi) { 129 PetscCall(MatSetPreallocationCOO(A, n2, it, jt)); 130 } else { 131 PetscCall(MatSetPreallocationCOOLocal(A, n2, it, jt)); 132 } 133 PetscCall(MatSetValuesCOO(A, v1, ADD_VALUES)); 134 PetscCall(MatMult(A, x, y)); 135 PetscCall(MatView(A, NULL)); 136 PetscCall(VecView(y, NULL)); 137 PetscCall(MatSetValuesCOO(A, v2, ADD_VALUES)); 138 PetscCall(MatMultAdd(A, x, y, z)); 139 PetscCall(MatView(A, NULL)); 140 PetscCall(VecView(z, NULL)); 141 PetscCall(PetscArraycpy(it, i2, n2)); 142 PetscCall(PetscArraycpy(jt, j2, n2)); 143 if (!localapi) { 144 PetscCall(MatSetPreallocationCOO(A, n2, it, jt)); 145 } else { 146 PetscCall(MatSetPreallocationCOOLocal(A, n2, it, jt)); 147 } 148 PetscCall(MatSetValuesCOO(A, v1, INSERT_VALUES)); 149 PetscCall(MatMult(A, x, y)); 150 PetscCall(MatView(A, NULL)); 151 PetscCall(VecView(y, NULL)); 152 PetscCall(MatSetValuesCOO(A, v2, INSERT_VALUES)); 153 PetscCall(MatMultAdd(A, x, y, z)); 154 PetscCall(MatView(A, NULL)); 155 PetscCall(VecView(z, NULL)); 156 T = A; 157 if (ishypre) PetscCall(MatConvert(A, MATAIJ, MAT_INITIAL_MATRIX, &T)); 158 PetscCall(MatTranspose(T, MAT_INITIAL_MATRIX, &At)); 159 if (!ismatis) { 160 PetscCall(MatMatMult(T, At, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &AAt)); 161 PetscCall(MatView(AAt, NULL)); 162 PetscCall(MatDestroy(&AAt)); 163 PetscCall(MatMatMult(At, T, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &AAt)); 164 PetscCall(MatView(AAt, NULL)); 165 PetscCall(MatDestroy(&AAt)); 166 } 167 PetscCall(MatDestroy(&At)); 168 if (ishypre) PetscCall(MatDestroy(&T)); 169 170 /* test providing diagonal first, then offdiagonal */ 171 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size)); 172 PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIAIJ, &ismpiaij)); 173 if ((ismpiaij || ishypre) && size > 1) { 174 Mat lA, lB; 175 const PetscInt *garray, *iA, *jA, *iB, *jB; 176 const PetscScalar *vA, *vB; 177 PetscScalar *coo_v; 178 PetscInt *coo_i, *coo_j; 179 PetscInt i, j, nA, nB, nnz; 180 PetscBool flg; 181 182 T = A; 183 if (ishypre) PetscCall(MatConvert(A, MATAIJ, MAT_INITIAL_MATRIX, &T)); 184 PetscCall(MatMPIAIJGetSeqAIJ(T, &lA, &lB, &garray)); 185 PetscCall(MatSeqAIJGetArrayRead(lA, &vA)); 186 PetscCall(MatSeqAIJGetArrayRead(lB, &vB)); 187 PetscCall(MatGetRowIJ(lA, 0, PETSC_FALSE, PETSC_FALSE, &nA, &iA, &jA, &flg)); 188 PetscCall(MatGetRowIJ(lB, 0, PETSC_FALSE, PETSC_FALSE, &nB, &iB, &jB, &flg)); 189 nnz = iA[nA] + iB[nB]; 190 PetscCall(PetscMalloc3(nnz, &coo_i, nnz, &coo_j, nnz, &coo_v)); 191 nnz = 0; 192 for (i = 0; i < nA; i++) { 193 for (j = iA[i]; j < iA[i + 1]; j++, nnz++) { 194 coo_i[nnz] = i + rstart; 195 coo_j[nnz] = jA[j] + cstart; 196 coo_v[nnz] = vA[j]; 197 } 198 } 199 for (i = 0; i < nB; i++) { 200 for (j = iB[i]; j < iB[i + 1]; j++, nnz++) { 201 coo_i[nnz] = i + rstart; 202 coo_j[nnz] = garray[jB[j]]; 203 coo_v[nnz] = vB[j]; 204 } 205 } 206 PetscCall(MatRestoreRowIJ(lA, 0, PETSC_FALSE, PETSC_FALSE, &nA, &iA, &jA, &flg)); 207 PetscCall(MatRestoreRowIJ(lB, 0, PETSC_FALSE, PETSC_FALSE, &nB, &iB, &jB, &flg)); 208 PetscCall(MatSeqAIJRestoreArrayRead(lA, &vA)); 209 PetscCall(MatSeqAIJRestoreArrayRead(lB, &vB)); 210 if (ishypre) PetscCall(MatDestroy(&T)); 211 212 PetscCall(MatSetPreallocationCOO(A, nnz, coo_i, coo_j)); 213 PetscCall(MatSetValuesCOO(A, coo_v, ADD_VALUES)); 214 PetscCall(MatMult(A, x, y)); 215 PetscCall(MatView(A, NULL)); 216 PetscCall(VecView(y, NULL)); 217 PetscCall(MatSetValuesCOO(A, coo_v, INSERT_VALUES)); 218 PetscCall(MatMult(A, x, y)); 219 PetscCall(MatView(A, NULL)); 220 PetscCall(VecView(y, NULL)); 221 222 T = A; 223 if (ishypre) PetscCall(MatConvert(A, MATAIJ, MAT_INITIAL_MATRIX, &T)); 224 PetscCall(MatTranspose(T, MAT_INITIAL_MATRIX, &At)); 225 PetscCall(MatMatMult(T, At, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &AAt)); 226 PetscCall(MatView(AAt, NULL)); 227 PetscCall(MatDestroy(&AAt)); 228 PetscCall(MatMatMult(At, T, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &AAt)); 229 PetscCall(MatView(AAt, NULL)); 230 PetscCall(MatDestroy(&AAt)); 231 PetscCall(MatDestroy(&At)); 232 if (ishypre) PetscCall(MatDestroy(&T)); 233 234 PetscCall(PetscFree3(coo_i, coo_j, coo_v)); 235 } 236 PetscCall(PetscFree2(it, jt)); 237 PetscCall(VecDestroy(&z)); 238 PetscCall(VecDestroy(&x)); 239 PetscCall(VecDestroy(&y)); 240 PetscCall(MatDestroy(&A)); 241 PetscCall(PetscFinalize()); 242 return 0; 243 } 244 245 /*TEST 246 247 test: 248 suffix: 1 249 filter: grep -v type | grep -v "Mat Object" 250 diff_args: -j 251 args: -mat_type {{seqaij mpiaij}} -localapi {{0 1}} -neg {{0 1}} 252 253 test: 254 requires: hypre !defined(PETSC_HAVE_HYPRE_DEVICE) 255 suffix: 1_hypre 256 filter: grep -v type | grep -v "Mat Object" 257 diff_args: -j 258 args: -mat_type hypre -localapi {{0 1}} -neg {{0 1}} 259 output_file: output/ex123_1.out 260 261 test: 262 requires: cuda 263 suffix: 1_cuda 264 filter: grep -v type | grep -v "Mat Object" 265 diff_args: -j 266 args: -mat_type {{seqaijcusparse mpiaijcusparse}} -localapi {{0 1}} -neg {{0 1}} 267 output_file: output/ex123_1.out 268 269 test: 270 requires: kokkos_kernels 271 suffix: 1_kokkos 272 filter: grep -v type | grep -v "Mat Object" 273 diff_args: -j 274 args: -mat_type {{seqaijkokkos mpiaijkokkos}} -localapi {{0 1}} -neg {{0 1}} 275 output_file: output/ex123_1.out 276 277 test: 278 suffix: 2 279 nsize: 7 280 filter: grep -v type | grep -v "Mat Object" 281 diff_args: -j 282 args: -mat_type mpiaij -localapi {{0 1}} -neg {{0 1}} 283 284 test: 285 requires: hypre !defined(PETSC_HAVE_HYPRE_DEVICE) 286 suffix: 2_hypre 287 nsize: 7 288 filter: grep -v type | grep -v "Mat Object" 289 diff_args: -j 290 args: -mat_type hypre -localapi {{0 1}} -neg {{0 1}} 291 output_file: output/ex123_2.out 292 293 test: 294 requires: cuda 295 suffix: 2_cuda 296 nsize: 7 297 filter: grep -v type | grep -v "Mat Object" 298 diff_args: -j 299 args: -mat_type mpiaijcusparse -localapi {{0 1}} -neg {{0 1}} 300 output_file: output/ex123_2.out 301 302 test: 303 requires: kokkos_kernels 304 suffix: 2_kokkos 305 nsize: 7 306 filter: grep -v type | grep -v "Mat Object" 307 diff_args: -j 308 args: -mat_type mpiaijkokkos -localapi {{0 1}} -neg {{0 1}} 309 output_file: output/ex123_2.out 310 311 test: 312 suffix: 3 313 nsize: 3 314 filter: grep -v type | grep -v "Mat Object" 315 diff_args: -j 316 args: -mat_type mpiaij -loc -localapi {{0 1}} -neg {{0 1}} 317 318 test: 319 requires: hypre !defined(PETSC_HAVE_HYPRE_DEVICE) 320 suffix: 3_hypre 321 nsize: 3 322 filter: grep -v type | grep -v "Mat Object" 323 diff_args: -j 324 args: -mat_type hypre -loc -localapi {{0 1}} -neg {{0 1}} 325 output_file: output/ex123_3.out 326 327 test: 328 requires: cuda 329 suffix: 3_cuda 330 nsize: 3 331 filter: grep -v type | grep -v "Mat Object" 332 diff_args: -j 333 args: -mat_type mpiaijcusparse -loc -localapi {{0 1}} -neg {{0 1}} 334 output_file: output/ex123_3.out 335 336 test: 337 requires: kokkos_kernels 338 suffix: 3_kokkos 339 nsize: 3 340 filter: grep -v type | grep -v "Mat Object" 341 diff_args: -j 342 args: -mat_type aijkokkos -loc -localapi {{0 1}} -neg {{0 1}} 343 output_file: output/ex123_3.out 344 345 test: 346 suffix: 4 347 nsize: 4 348 filter: grep -v type | grep -v "Mat Object" 349 diff_args: -j 350 args: -mat_type mpiaij -loc -locdiag 0 -localapi {{0 1}} -neg {{0 1}} 351 352 test: 353 requires: hypre !defined(PETSC_HAVE_HYPRE_DEVICE) 354 suffix: 4_hypre 355 nsize: 4 356 filter: grep -v type | grep -v "Mat Object" 357 diff_args: -j 358 args: -mat_type hypre -loc -locdiag 0 -localapi {{0 1}} -neg {{0 1}} 359 output_file: output/ex123_4.out 360 361 test: 362 requires: cuda 363 suffix: 4_cuda 364 nsize: 4 365 filter: grep -v type | grep -v "Mat Object" 366 diff_args: -j 367 args: -mat_type mpiaijcusparse -loc -locdiag 0 -localapi {{0 1}} -neg {{0 1}} 368 output_file: output/ex123_4.out 369 370 test: 371 requires: kokkos_kernels 372 suffix: 4_kokkos 373 nsize: 4 374 filter: grep -v type | grep -v "Mat Object" 375 diff_args: -j 376 args: -mat_type aijkokkos -loc -locdiag 0 -localapi {{0 1}} -neg {{0 1}} 377 output_file: output/ex123_4.out 378 379 test: 380 suffix: matis 381 nsize: 3 382 filter: grep -v type | grep -v "Mat Object" 383 diff_args: -j 384 args: -mat_type is -localapi {{0 1}} -neg {{0 1}} 385 386 TEST*/ 387