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 PetscErrorCode ierr; 12 PetscInt i,j,m = 8,n,rstart,rend,vstart,vend; 13 PetscScalar one = 1.0,negone = -1.0,v,alpha=0.1; 14 PetscReal norm, tol = PETSC_SQRT_MACHINE_EPSILON; 15 PetscBool flg; 16 17 ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr; 18 CHKERRQ(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_COMMON)); 19 CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL)); 20 n = m; 21 CHKERRQ(PetscOptionsHasName(NULL,NULL,"-rectA",&flg)); 22 if (flg) n += 2; 23 CHKERRQ(PetscOptionsHasName(NULL,NULL,"-rectB",&flg)); 24 if (flg) n -= 2; 25 26 /* ---------- Assemble matrix and vectors ----------- */ 27 28 CHKERRQ(MatCreate(PETSC_COMM_WORLD,&C)); 29 CHKERRQ(MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,m,n)); 30 CHKERRQ(MatSetFromOptions(C)); 31 CHKERRQ(MatSetUp(C)); 32 CHKERRQ(MatGetOwnershipRange(C,&rstart,&rend)); 33 CHKERRQ(VecCreate(PETSC_COMM_WORLD,&x)); 34 CHKERRQ(VecSetSizes(x,PETSC_DECIDE,m)); 35 CHKERRQ(VecSetFromOptions(x)); 36 CHKERRQ(VecDuplicate(x,&z)); 37 CHKERRQ(VecDuplicate(x,&w)); 38 CHKERRQ(VecCreate(PETSC_COMM_WORLD,&y)); 39 CHKERRQ(VecSetSizes(y,PETSC_DECIDE,n)); 40 CHKERRQ(VecSetFromOptions(y)); 41 CHKERRQ(VecDuplicate(y,&u)); 42 CHKERRQ(VecDuplicate(y,&s)); 43 CHKERRQ(VecGetOwnershipRange(y,&vstart,&vend)); 44 45 /* Assembly */ 46 for (i=rstart; i<rend; i++) { 47 v = 100*(i+1); 48 CHKERRQ(VecSetValues(z,1,&i,&v,INSERT_VALUES)); 49 for (j=0; j<n; j++) { 50 v = 10*(i+1)+j+1; 51 CHKERRQ(MatSetValues(C,1,&i,1,&j,&v,INSERT_VALUES)); 52 } 53 } 54 55 /* Flush off proc Vec values and do more assembly */ 56 CHKERRQ(VecAssemblyBegin(z)); 57 for (i=vstart; i<vend; i++) { 58 v = one*((PetscReal)i); 59 CHKERRQ(VecSetValues(y,1,&i,&v,INSERT_VALUES)); 60 v = 100.0*i; 61 CHKERRQ(VecSetValues(u,1,&i,&v,INSERT_VALUES)); 62 } 63 64 /* Flush off proc Mat values and do more assembly */ 65 CHKERRQ(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 CHKERRQ(MatSetValues(C,1,&i,1,&j,&v,INSERT_VALUES)); 70 } 71 } 72 /* Try overlap Coomunication with the next stage XXXSetValues */ 73 CHKERRQ(VecAssemblyEnd(z)); 74 75 CHKERRQ(MatAssemblyEnd(C,MAT_FLUSH_ASSEMBLY)); 76 CHKMEMQ; 77 /* The Assembly for the second Stage */ 78 CHKERRQ(MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY)); 79 CHKERRQ(MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY)); 80 CHKERRQ(VecAssemblyBegin(y)); 81 CHKERRQ(VecAssemblyEnd(y)); 82 CHKERRQ(MatScale(C,alpha)); 83 CHKERRQ(VecAssemblyBegin(u)); 84 CHKERRQ(VecAssemblyEnd(u)); 85 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"testing MatMult()\n")); 86 CHKMEMQ; 87 CHKERRQ(MatMult(C,y,x)); 88 CHKMEMQ; 89 CHKERRQ(VecView(x,PETSC_VIEWER_STDOUT_WORLD)); 90 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"testing MatMultAdd()\n")); 91 CHKERRQ(MatMultAdd(C,y,z,w)); 92 CHKERRQ(VecAXPY(x,one,z)); 93 CHKERRQ(VecAXPY(x,negone,w)); 94 CHKERRQ(VecNorm(x,NORM_2,&norm)); 95 if (norm > tol) { 96 CHKERRQ(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 CHKERRQ(VecSetValues(x,1,&i,&v,INSERT_VALUES)); 104 } 105 CHKERRQ(VecAssemblyBegin(x)); 106 CHKERRQ(VecAssemblyEnd(x)); 107 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"testing MatMultTranspose()\n")); 108 CHKERRQ(MatMultTranspose(C,x,y)); 109 CHKERRQ(VecView(y,PETSC_VIEWER_STDOUT_WORLD)); 110 111 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"testing MatMultTransposeAdd()\n")); 112 CHKERRQ(MatMultTransposeAdd(C,x,u,s)); 113 CHKERRQ(VecAXPY(y,one,u)); 114 CHKERRQ(VecAXPY(y,negone,s)); 115 CHKERRQ(VecNorm(y,NORM_2,&norm)); 116 if (norm > tol) { 117 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"Norm of error difference = %g\n",(double)norm)); 118 } 119 120 /* -------------------- Test MatGetDiagonal() ------------------ */ 121 122 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"testing MatGetDiagonal(), MatDiagonalScale()\n")); 123 CHKERRQ(MatView(C,PETSC_VIEWER_STDOUT_WORLD)); 124 CHKERRQ(VecSet(x,one)); 125 CHKERRQ(MatGetDiagonal(C,x)); 126 CHKERRQ(VecView(x,PETSC_VIEWER_STDOUT_WORLD)); 127 for (i=vstart; i<vend; i++) { 128 v = one*((PetscReal)(i+1)); 129 CHKERRQ(VecSetValues(y,1,&i,&v,INSERT_VALUES)); 130 } 131 132 /* -------------------- Test () MatDiagonalScale ------------------ */ 133 CHKERRQ(PetscOptionsHasName(NULL,NULL,"-test_diagonalscale",&flg)); 134 if (flg) { 135 CHKERRQ(MatDiagonalScale(C,x,y)); 136 CHKERRQ(MatView(C,PETSC_VIEWER_STDOUT_WORLD)); 137 } 138 /* Free data structures */ 139 CHKERRQ(VecDestroy(&u)); CHKERRQ(VecDestroy(&s)); 140 CHKERRQ(VecDestroy(&w)); CHKERRQ(VecDestroy(&x)); 141 CHKERRQ(VecDestroy(&y)); CHKERRQ(VecDestroy(&z)); 142 CHKERRQ(MatDestroy(&C)); 143 144 ierr = PetscFinalize(); 145 return ierr; 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