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 { 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 /* Free data structures */ 135 PetscCall(VecDestroy(&u)); 136 PetscCall(VecDestroy(&s)); 137 PetscCall(VecDestroy(&w)); 138 PetscCall(VecDestroy(&x)); 139 PetscCall(VecDestroy(&y)); 140 PetscCall(VecDestroy(&z)); 141 PetscCall(MatDestroy(&C)); 142 143 PetscCall(PetscFinalize()); 144 return 0; 145 } 146 147 /*TEST 148 149 test: 150 suffix: 11_A 151 args: -mat_type seqaij -rectA 152 filter: grep -v type 153 154 test: 155 args: -mat_type seqdense -rectA 156 suffix: 12_A 157 158 test: 159 args: -mat_type seqaij -rectB 160 suffix: 11_B 161 filter: grep -v type 162 163 test: 164 args: -mat_type seqdense -rectB 165 suffix: 12_B 166 167 test: 168 suffix: 21 169 args: -mat_type mpiaij 170 filter: grep -v type 171 172 test: 173 suffix: 22 174 args: -mat_type mpidense 175 176 test: 177 suffix: 23 178 nsize: 3 179 args: -mat_type mpiaij 180 filter: grep -v type 181 182 test: 183 suffix: 24 184 nsize: 3 185 args: -mat_type mpidense 186 187 test: 188 suffix: 2_aijcusparse_1 189 args: -mat_type mpiaijcusparse -vec_type cuda 190 filter: grep -v type 191 output_file: output/ex5_21.out 192 requires: cuda 193 194 test: 195 nsize: 3 196 suffix: 2_aijcusparse_2 197 filter: grep -v type 198 args: -mat_type mpiaijcusparse -vec_type cuda 199 args: -sf_type {{basic neighbor}} 200 output_file: output/ex5_23.out 201 requires: cuda 202 203 test: 204 nsize: 3 205 suffix: 2_aijcusparse_3 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 defined(PETSC_HAVE_MPI_GPU_AWARE) 211 212 test: 213 suffix: 31 214 args: -mat_type mpiaij -test_diagonalscale 215 filter: grep -v type 216 217 test: 218 suffix: 32 219 args: -mat_type mpibaij -test_diagonalscale 220 filter: grep -v Mat_ 221 222 test: 223 suffix: 33 224 nsize: 3 225 args: -mat_type mpiaij -test_diagonalscale 226 filter: grep -v type 227 228 test: 229 suffix: 34 230 nsize: 3 231 args: -mat_type mpibaij -test_diagonalscale 232 filter: grep -v Mat_ 233 234 test: 235 suffix: 3_aijcusparse_1 236 args: -mat_type mpiaijcusparse -vec_type cuda -test_diagonalscale 237 filter: grep -v type 238 output_file: output/ex5_31.out 239 requires: cuda 240 241 test: 242 suffix: 3_aijcusparse_2 243 nsize: 3 244 args: -mat_type mpiaijcusparse -vec_type cuda -test_diagonalscale 245 filter: grep -v type 246 output_file: output/ex5_33.out 247 requires: cuda 248 249 test: 250 suffix: 3_kokkos 251 nsize: 3 252 args: -mat_type mpiaijkokkos -vec_type kokkos -test_diagonalscale 253 filter: grep -v type 254 output_file: output/ex5_33.out 255 requires: kokkos_kernels 256 257 test: 258 suffix: aijcusparse_1 259 args: -mat_type seqaijcusparse -vec_type cuda -rectA 260 filter: grep -v type 261 output_file: output/ex5_11_A.out 262 requires: cuda 263 264 test: 265 suffix: aijcusparse_2 266 args: -mat_type seqaijcusparse -vec_type cuda -rectB 267 filter: grep -v type 268 output_file: output/ex5_11_B.out 269 requires: cuda 270 271 test: 272 suffix: sell_1 273 args: -mat_type sell 274 output_file: output/ex5_41.out 275 276 test: 277 suffix: sell_2 278 nsize: 3 279 args: -mat_type sell 280 output_file: output/ex5_43.out 281 282 test: 283 suffix: sell_3 284 args: -mat_type sell -test_diagonalscale 285 output_file: output/ex5_51.out 286 287 test: 288 suffix: sell_4 289 nsize: 3 290 args: -mat_type sell -test_diagonalscale 291 output_file: output/ex5_53.out 292 293 TEST*/ 294