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) { 96 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Norm of error difference = %g\n",(double)norm)); 97 } 98 99 /* ------- Test MatMultTranspose(), MatMultTransposeAdd() ------- */ 100 101 for (i=rstart; i<rend; i++) { 102 v = one*((PetscReal)i); 103 PetscCall(VecSetValues(x,1,&i,&v,INSERT_VALUES)); 104 } 105 PetscCall(VecAssemblyBegin(x)); 106 PetscCall(VecAssemblyEnd(x)); 107 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"testing MatMultTranspose()\n")); 108 PetscCall(MatMultTranspose(C,x,y)); 109 PetscCall(VecView(y,PETSC_VIEWER_STDOUT_WORLD)); 110 111 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"testing MatMultTransposeAdd()\n")); 112 PetscCall(MatMultTransposeAdd(C,x,u,s)); 113 PetscCall(VecAXPY(y,one,u)); 114 PetscCall(VecAXPY(y,negone,s)); 115 PetscCall(VecNorm(y,NORM_2,&norm)); 116 if (norm > tol) { 117 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Norm of error difference = %g\n",(double)norm)); 118 } 119 120 /* -------------------- Test MatGetDiagonal() ------------------ */ 121 122 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"testing MatGetDiagonal(), MatDiagonalScale()\n")); 123 PetscCall(MatView(C,PETSC_VIEWER_STDOUT_WORLD)); 124 PetscCall(VecSet(x,one)); 125 PetscCall(MatGetDiagonal(C,x)); 126 PetscCall(VecView(x,PETSC_VIEWER_STDOUT_WORLD)); 127 for (i=vstart; i<vend; i++) { 128 v = one*((PetscReal)(i+1)); 129 PetscCall(VecSetValues(y,1,&i,&v,INSERT_VALUES)); 130 } 131 132 /* -------------------- Test () MatDiagonalScale ------------------ */ 133 PetscCall(PetscOptionsHasName(NULL,NULL,"-test_diagonalscale",&flg)); 134 if (flg) { 135 PetscCall(MatDiagonalScale(C,x,y)); 136 PetscCall(MatView(C,PETSC_VIEWER_STDOUT_WORLD)); 137 } 138 /* Free data structures */ 139 PetscCall(VecDestroy(&u)); PetscCall(VecDestroy(&s)); 140 PetscCall(VecDestroy(&w)); PetscCall(VecDestroy(&x)); 141 PetscCall(VecDestroy(&y)); PetscCall(VecDestroy(&z)); 142 PetscCall(MatDestroy(&C)); 143 144 PetscCall(PetscFinalize()); 145 return 0; 146 } 147 148 /*TEST 149 150 test: 151 suffix: 11_A 152 args: -mat_type seqaij -rectA 153 filter: grep -v type 154 155 test: 156 args: -mat_type seqdense -rectA 157 suffix: 12_A 158 159 test: 160 args: -mat_type seqaij -rectB 161 suffix: 11_B 162 filter: grep -v type 163 164 test: 165 args: -mat_type seqdense -rectB 166 suffix: 12_B 167 168 test: 169 suffix: 21 170 args: -mat_type mpiaij 171 filter: grep -v type 172 173 test: 174 suffix: 22 175 args: -mat_type mpidense 176 177 test: 178 suffix: 23 179 nsize: 3 180 args: -mat_type mpiaij 181 filter: grep -v type 182 183 test: 184 suffix: 24 185 nsize: 3 186 args: -mat_type mpidense 187 188 test: 189 suffix: 2_aijcusparse_1 190 args: -mat_type mpiaijcusparse -vec_type cuda 191 filter: grep -v type 192 output_file: output/ex5_21.out 193 requires: cuda 194 195 test: 196 nsize: 3 197 suffix: 2_aijcusparse_2 198 filter: grep -v type 199 args: -mat_type mpiaijcusparse -vec_type cuda 200 args: -sf_type {{basic neighbor}} 201 output_file: output/ex5_23.out 202 requires: cuda 203 204 test: 205 nsize: 3 206 suffix: 2_aijcusparse_3 207 filter: grep -v type 208 args: -mat_type mpiaijcusparse -vec_type cuda 209 args: -sf_type {{basic neighbor}} 210 output_file: output/ex5_23.out 211 requires: cuda defined(PETSC_HAVE_MPI_GPU_AWARE) 212 213 test: 214 suffix: 31 215 args: -mat_type mpiaij -test_diagonalscale 216 filter: grep -v type 217 218 test: 219 suffix: 32 220 args: -mat_type mpibaij -test_diagonalscale 221 filter: grep -v Mat_ 222 223 test: 224 suffix: 33 225 nsize: 3 226 args: -mat_type mpiaij -test_diagonalscale 227 filter: grep -v type 228 229 test: 230 suffix: 34 231 nsize: 3 232 args: -mat_type mpibaij -test_diagonalscale 233 filter: grep -v Mat_ 234 235 test: 236 suffix: 3_aijcusparse_1 237 args: -mat_type mpiaijcusparse -vec_type cuda -test_diagonalscale 238 filter: grep -v type 239 output_file: output/ex5_31.out 240 requires: cuda 241 242 test: 243 suffix: 3_aijcusparse_2 244 nsize: 3 245 args: -mat_type mpiaijcusparse -vec_type cuda -test_diagonalscale 246 filter: grep -v type 247 output_file: output/ex5_33.out 248 requires: cuda 249 250 test: 251 suffix: 3_kokkos 252 nsize: 3 253 args: -mat_type mpiaijkokkos -vec_type kokkos -test_diagonalscale 254 filter: grep -v type 255 output_file: output/ex5_33.out 256 requires: !sycl kokkos_kernels 257 258 test: 259 suffix: aijcusparse_1 260 args: -mat_type seqaijcusparse -vec_type cuda -rectA 261 filter: grep -v type 262 output_file: output/ex5_11_A.out 263 requires: cuda 264 265 test: 266 suffix: aijcusparse_2 267 args: -mat_type seqaijcusparse -vec_type cuda -rectB 268 filter: grep -v type 269 output_file: output/ex5_11_B.out 270 requires: cuda 271 272 test: 273 suffix: sell_1 274 args: -mat_type sell 275 output_file: output/ex5_41.out 276 277 test: 278 suffix: sell_2 279 nsize: 3 280 args: -mat_type sell 281 output_file: output/ex5_43.out 282 283 test: 284 suffix: sell_3 285 args: -mat_type sell -test_diagonalscale 286 output_file: output/ex5_51.out 287 288 test: 289 suffix: sell_4 290 nsize: 3 291 args: -mat_type sell -test_diagonalscale 292 output_file: output/ex5_53.out 293 294 TEST*/ 295