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