1 2 static char help[] = "Tests MatMult(), MatMultAdd(), MatMultTranspose().\n\ 3 Also MatMultTransposeAdd(), MatScale(), MatGetDiagonal(), MatDiagonalScale(), MatZeroEntries() and MatDuplicate().\n\n"; 4 5 #include <petscmat.h> 6 7 int main(int argc, char **args) 8 { 9 Mat C; 10 Vec s, u, w, x, y, z; 11 PetscInt i, j, m = 8, n, rstart, rend, vstart, vend; 12 PetscScalar one = 1.0, negone = -1.0, v, alpha = 0.1; 13 PetscReal norm, tol = PETSC_SQRT_MACHINE_EPSILON; 14 PetscBool flg; 15 16 PetscFunctionBeginUser; 17 PetscCall(PetscInitialize(&argc, &args, (char *)0, help)); 18 PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_COMMON)); 19 PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL)); 20 n = m; 21 PetscCall(PetscOptionsHasName(NULL, NULL, "-rectA", &flg)); 22 if (flg) n += 2; 23 PetscCall(PetscOptionsHasName(NULL, NULL, "-rectB", &flg)); 24 if (flg) n -= 2; 25 26 /* ---------- Assemble matrix and vectors ----------- */ 27 28 PetscCall(MatCreate(PETSC_COMM_WORLD, &C)); 29 PetscCall(MatSetSizes(C, PETSC_DECIDE, PETSC_DECIDE, m, n)); 30 PetscCall(MatSetFromOptions(C)); 31 PetscCall(MatSetUp(C)); 32 PetscCall(MatGetOwnershipRange(C, &rstart, &rend)); 33 PetscCall(VecCreate(PETSC_COMM_WORLD, &x)); 34 PetscCall(VecSetSizes(x, PETSC_DECIDE, m)); 35 PetscCall(VecSetFromOptions(x)); 36 PetscCall(VecDuplicate(x, &z)); 37 PetscCall(VecDuplicate(x, &w)); 38 PetscCall(VecCreate(PETSC_COMM_WORLD, &y)); 39 PetscCall(VecSetSizes(y, PETSC_DECIDE, n)); 40 PetscCall(VecSetFromOptions(y)); 41 PetscCall(VecDuplicate(y, &u)); 42 PetscCall(VecDuplicate(y, &s)); 43 PetscCall(VecGetOwnershipRange(y, &vstart, &vend)); 44 45 /* Assembly */ 46 for (i = rstart; i < rend; i++) { 47 v = 100 * (i + 1); 48 PetscCall(VecSetValues(z, 1, &i, &v, INSERT_VALUES)); 49 for (j = 0; j < n; j++) { 50 v = 10 * (i + 1) + j + 1; 51 PetscCall(MatSetValues(C, 1, &i, 1, &j, &v, INSERT_VALUES)); 52 } 53 } 54 55 /* Flush off proc Vec values and do more assembly */ 56 PetscCall(VecAssemblyBegin(z)); 57 for (i = vstart; i < vend; i++) { 58 v = one * ((PetscReal)i); 59 PetscCall(VecSetValues(y, 1, &i, &v, INSERT_VALUES)); 60 v = 100.0 * i; 61 PetscCall(VecSetValues(u, 1, &i, &v, INSERT_VALUES)); 62 } 63 64 /* Flush off proc Mat values and do more assembly */ 65 PetscCall(MatAssemblyBegin(C, MAT_FLUSH_ASSEMBLY)); 66 for (i = rstart; i < rend; i++) { 67 for (j = 0; j < n; j++) { 68 v = 10 * (i + 1) + j + 1; 69 PetscCall(MatSetValues(C, 1, &i, 1, &j, &v, INSERT_VALUES)); 70 } 71 } 72 /* Try overlap Coomunication with the next stage XXXSetValues */ 73 PetscCall(VecAssemblyEnd(z)); 74 75 PetscCall(MatAssemblyEnd(C, MAT_FLUSH_ASSEMBLY)); 76 CHKMEMQ; 77 /* The Assembly for the second Stage */ 78 PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY)); 79 PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY)); 80 PetscCall(VecAssemblyBegin(y)); 81 PetscCall(VecAssemblyEnd(y)); 82 PetscCall(MatScale(C, alpha)); 83 PetscCall(VecAssemblyBegin(u)); 84 PetscCall(VecAssemblyEnd(u)); 85 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "testing MatMult()\n")); 86 CHKMEMQ; 87 PetscCall(MatMult(C, y, x)); 88 CHKMEMQ; 89 PetscCall(VecView(x, PETSC_VIEWER_STDOUT_WORLD)); 90 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "testing MatMultAdd()\n")); 91 PetscCall(MatMultAdd(C, y, z, w)); 92 PetscCall(VecAXPY(x, one, z)); 93 PetscCall(VecAXPY(x, negone, w)); 94 PetscCall(VecNorm(x, NORM_2, &norm)); 95 if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Norm of error difference = %g\n", (double)norm)); 96 97 /* ------- Test MatMultTranspose(), MatMultTransposeAdd() ------- */ 98 99 for (i = rstart; i < rend; i++) { 100 v = one * ((PetscReal)i); 101 PetscCall(VecSetValues(x, 1, &i, &v, INSERT_VALUES)); 102 } 103 PetscCall(VecAssemblyBegin(x)); 104 PetscCall(VecAssemblyEnd(x)); 105 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "testing MatMultTranspose()\n")); 106 PetscCall(MatMultTranspose(C, x, y)); 107 PetscCall(VecView(y, PETSC_VIEWER_STDOUT_WORLD)); 108 109 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "testing MatMultTransposeAdd()\n")); 110 PetscCall(MatMultTransposeAdd(C, x, u, s)); 111 PetscCall(VecAXPY(y, one, u)); 112 PetscCall(VecAXPY(y, negone, s)); 113 PetscCall(VecNorm(y, NORM_2, &norm)); 114 if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Norm of error difference = %g\n", (double)norm)); 115 116 /* -------------------- Test MatGetDiagonal() ------------------ */ 117 118 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "testing MatGetDiagonal(), MatDiagonalScale()\n")); 119 PetscCall(MatView(C, PETSC_VIEWER_STDOUT_WORLD)); 120 PetscCall(VecSet(x, one)); 121 PetscCall(MatGetDiagonal(C, x)); 122 PetscCall(VecView(x, PETSC_VIEWER_STDOUT_WORLD)); 123 for (i = vstart; i < vend; i++) { 124 v = one * ((PetscReal)(i + 1)); 125 PetscCall(VecSetValues(y, 1, &i, &v, INSERT_VALUES)); 126 } 127 128 /* -------------------- Test () MatDiagonalScale ------------------ */ 129 PetscCall(PetscOptionsHasName(NULL, NULL, "-test_diagonalscale", &flg)); 130 if (flg) { 131 PetscCall(MatDiagonalScale(C, x, y)); 132 PetscCall(MatView(C, PETSC_VIEWER_STDOUT_WORLD)); 133 } 134 /* -------------------- Test () MatZeroEntries() and MatDuplicate() ------------------ */ 135 PetscCall(PetscOptionsHasName(NULL, NULL, "-test_zeroentries", &flg)); 136 if (flg) { 137 Mat D; 138 PetscCall(MatDuplicate(C, MAT_COPY_VALUES, &D)); 139 PetscCall(MatZeroEntries(D)); 140 PetscCall(MatView(D, PETSC_VIEWER_STDOUT_WORLD)); 141 PetscCall(MatDestroy(&D)); 142 } 143 /* Free data structures */ 144 PetscCall(VecDestroy(&u)); 145 PetscCall(VecDestroy(&s)); 146 PetscCall(VecDestroy(&w)); 147 PetscCall(VecDestroy(&x)); 148 PetscCall(VecDestroy(&y)); 149 PetscCall(VecDestroy(&z)); 150 PetscCall(MatDestroy(&C)); 151 152 PetscCall(PetscFinalize()); 153 return 0; 154 } 155 156 /*TEST 157 158 test: 159 suffix: 11_A 160 args: -mat_type seqaij -rectA 161 filter: grep -v type 162 163 test: 164 args: -mat_type seqdense -rectA 165 suffix: 12_A 166 167 test: 168 args: -mat_type seqaij -rectB 169 suffix: 11_B 170 filter: grep -v type 171 172 test: 173 args: -mat_type seqdense -rectB 174 suffix: 12_B 175 176 test: 177 suffix: 21 178 args: -mat_type mpiaij 179 filter: grep -v type 180 181 test: 182 suffix: 22 183 args: -mat_type mpidense 184 185 test: 186 suffix: 23 187 nsize: 3 188 args: -mat_type mpiaij 189 filter: grep -v type 190 191 test: 192 suffix: 24 193 nsize: 3 194 args: -mat_type mpidense 195 196 test: 197 suffix: 2_aijcusparse_1 198 args: -mat_type mpiaijcusparse -vec_type cuda 199 filter: grep -v type 200 output_file: output/ex5_21.out 201 requires: cuda 202 203 test: 204 nsize: 3 205 suffix: 2_aijcusparse_2 206 filter: grep -v type 207 args: -mat_type mpiaijcusparse -vec_type cuda 208 args: -sf_type {{basic neighbor}} 209 output_file: output/ex5_23.out 210 requires: cuda 211 212 test: 213 nsize: 3 214 suffix: 2_aijcusparse_3 215 filter: grep -v type 216 args: -mat_type mpiaijcusparse -vec_type cuda 217 args: -sf_type {{basic neighbor}} 218 output_file: output/ex5_23.out 219 requires: cuda defined(PETSC_HAVE_MPI_GPU_AWARE) 220 221 test: 222 suffix: 31 223 args: -mat_type mpiaij -test_diagonalscale 224 filter: grep -v type 225 226 test: 227 suffix: 32 228 args: -mat_type mpibaij -test_diagonalscale 229 230 test: 231 suffix: 33 232 nsize: 3 233 args: -mat_type mpiaij -test_diagonalscale 234 filter: grep -v type 235 236 test: 237 suffix: 34 238 nsize: 3 239 args: -mat_type mpibaij -test_diagonalscale 240 241 test: 242 suffix: 3_aijcusparse_1 243 args: -mat_type mpiaijcusparse -vec_type cuda -test_diagonalscale 244 filter: grep -v type 245 output_file: output/ex5_31.out 246 requires: cuda 247 248 test: 249 suffix: 3_aijcusparse_2 250 nsize: 3 251 args: -mat_type mpiaijcusparse -vec_type cuda -test_diagonalscale 252 filter: grep -v type 253 output_file: output/ex5_33.out 254 requires: cuda 255 256 test: 257 suffix: 3_kokkos 258 nsize: 3 259 args: -mat_type mpiaijkokkos -vec_type kokkos -test_diagonalscale 260 filter: grep -v type 261 output_file: output/ex5_33.out 262 requires: kokkos_kernels 263 264 test: 265 suffix: aijcusparse_1 266 args: -mat_type seqaijcusparse -vec_type cuda -rectA 267 filter: grep -v type 268 output_file: output/ex5_11_A.out 269 requires: cuda 270 271 test: 272 suffix: aijcusparse_2 273 args: -mat_type seqaijcusparse -vec_type cuda -rectB 274 filter: grep -v type 275 output_file: output/ex5_11_B.out 276 requires: cuda 277 278 test: 279 suffix: sell_1 280 args: -mat_type sell -mat_sell_slice_height 8 281 output_file: output/ex5_41.out 282 283 test: 284 suffix: sell_2 285 nsize: 3 286 args: -mat_type sell -mat_sell_slice_height 8 287 output_file: output/ex5_43.out 288 289 test: 290 suffix: sell_3 291 args: -mat_type sell -test_diagonalscale -mat_sell_slice_height 8 292 output_file: output/ex5_51.out 293 294 test: 295 suffix: sell_4 296 nsize: 3 297 args: -mat_type sell -test_diagonalscale -mat_sell_slice_height 8 298 output_file: output/ex5_53.out 299 300 test: 301 suffix: sell_5 302 nsize: 3 303 args: -mat_type sellcuda -vec_type cuda -test_diagonalscale -test_zeroentries 304 output_file: output/ex5_55.out 305 requires: cuda !complex 306 307 test: 308 suffix: sell_6 309 nsize: 3 310 args: -mat_type sellcuda -vec_type cuda -mat_sell_spmv_cuda_kernel {{1 2 3 4 5 6}} 311 output_file: output/ex5_56.out 312 requires: cuda !complex 313 314 test: 315 suffix: sell_7 316 args: -m 32 -mat_type sellcuda -vec_type cuda -mat_sell_spmv_cuda_kernel {{0 7 9}} -mat_sell_spmv_cuda_blocky {{2 4 8 16 32}} 317 output_file: output/ex5_57.out 318 requires: cuda !complex !single 319 TEST*/ 320