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 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) { 95 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Norm of error difference = %g\n",(double)norm)); 96 } 97 98 /* ------- Test MatMultTranspose(), MatMultTransposeAdd() ------- */ 99 100 for (i=rstart; i<rend; i++) { 101 v = one*((PetscReal)i); 102 PetscCall(VecSetValues(x,1,&i,&v,INSERT_VALUES)); 103 } 104 PetscCall(VecAssemblyBegin(x)); 105 PetscCall(VecAssemblyEnd(x)); 106 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"testing MatMultTranspose()\n")); 107 PetscCall(MatMultTranspose(C,x,y)); 108 PetscCall(VecView(y,PETSC_VIEWER_STDOUT_WORLD)); 109 110 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"testing MatMultTransposeAdd()\n")); 111 PetscCall(MatMultTransposeAdd(C,x,u,s)); 112 PetscCall(VecAXPY(y,one,u)); 113 PetscCall(VecAXPY(y,negone,s)); 114 PetscCall(VecNorm(y,NORM_2,&norm)); 115 if (norm > tol) { 116 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Norm of error difference = %g\n",(double)norm)); 117 } 118 119 /* -------------------- Test MatGetDiagonal() ------------------ */ 120 121 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"testing MatGetDiagonal(), MatDiagonalScale()\n")); 122 PetscCall(MatView(C,PETSC_VIEWER_STDOUT_WORLD)); 123 PetscCall(VecSet(x,one)); 124 PetscCall(MatGetDiagonal(C,x)); 125 PetscCall(VecView(x,PETSC_VIEWER_STDOUT_WORLD)); 126 for (i=vstart; i<vend; i++) { 127 v = one*((PetscReal)(i+1)); 128 PetscCall(VecSetValues(y,1,&i,&v,INSERT_VALUES)); 129 } 130 131 /* -------------------- Test () MatDiagonalScale ------------------ */ 132 PetscCall(PetscOptionsHasName(NULL,NULL,"-test_diagonalscale",&flg)); 133 if (flg) { 134 PetscCall(MatDiagonalScale(C,x,y)); 135 PetscCall(MatView(C,PETSC_VIEWER_STDOUT_WORLD)); 136 } 137 /* Free data structures */ 138 PetscCall(VecDestroy(&u)); PetscCall(VecDestroy(&s)); 139 PetscCall(VecDestroy(&w)); PetscCall(VecDestroy(&x)); 140 PetscCall(VecDestroy(&y)); 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: !sycl 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