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 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 232 test: 233 suffix: 3_aijcusparse_1 234 args: -mat_type mpiaijcusparse -vec_type cuda -test_diagonalscale 235 filter: grep -v type 236 output_file: output/ex5_31.out 237 requires: cuda 238 239 test: 240 suffix: 3_aijcusparse_2 241 nsize: 3 242 args: -mat_type mpiaijcusparse -vec_type cuda -test_diagonalscale 243 filter: grep -v type 244 output_file: output/ex5_33.out 245 requires: cuda 246 247 test: 248 suffix: 3_kokkos 249 nsize: 3 250 args: -mat_type mpiaijkokkos -vec_type kokkos -test_diagonalscale 251 filter: grep -v type 252 output_file: output/ex5_33.out 253 requires: kokkos_kernels 254 255 test: 256 suffix: aijcusparse_1 257 args: -mat_type seqaijcusparse -vec_type cuda -rectA 258 filter: grep -v type 259 output_file: output/ex5_11_A.out 260 requires: cuda 261 262 test: 263 suffix: aijcusparse_2 264 args: -mat_type seqaijcusparse -vec_type cuda -rectB 265 filter: grep -v type 266 output_file: output/ex5_11_B.out 267 requires: cuda 268 269 test: 270 suffix: sell_1 271 args: -mat_type sell 272 output_file: output/ex5_41.out 273 274 test: 275 suffix: sell_2 276 nsize: 3 277 args: -mat_type sell 278 output_file: output/ex5_43.out 279 280 test: 281 suffix: sell_3 282 args: -mat_type sell -test_diagonalscale 283 output_file: output/ex5_51.out 284 285 test: 286 suffix: sell_4 287 nsize: 3 288 args: -mat_type sell -test_diagonalscale 289 output_file: output/ex5_53.out 290 291 TEST*/ 292