1 static char help[] = "Test MatSetPreallocationCOO and MatSetValuesCOO\n\n"; 2 3 #include <petscmat.h> 4 #define MyMatView(a,b) (PetscPrintf(PetscObjectComm((PetscObject)(a)),"LINE %d\n",__LINE__),MatView(a,b)); 5 #define MyVecView(a,b) (PetscPrintf(PetscObjectComm((PetscObject)(a)),"LINE %d\n",__LINE__),VecView(a,b)); 6 int main(int argc,char **args) 7 { 8 Mat A,At,AAt; 9 Vec x,y,z; 10 ISLocalToGlobalMapping rl2g,cl2g; 11 IS is; 12 PetscLayout rmap,cmap; 13 PetscInt *it,*jt; 14 PetscInt n1 = 11, n2 = 9; 15 PetscInt i1[] = { 7, 6, 2, 0, 4, 1, 1, 0, 2, 2, 1 , -1, -1}; 16 PetscInt j1[] = { 1, 4, 3, 5, 3, 3, 4, 5, 0, 3, 1 , -1, -1}; 17 PetscInt i2[] = { 7, 6, 2, 0, 4, 1, 1, 2, 1, -1, -1}; 18 PetscInt j2[] = { 1, 4, 3, 5, 3, 3, 4, 0, 1, -1, -1}; 19 PetscScalar v1[] = { -1., 1., 2., 3., 4., 5., 6., 7., 8., 9., 10., PETSC_MAX_REAL, PETSC_MAX_REAL}; 20 PetscScalar v2[] = { 1.,-1.,-2.,-3.,-4.,-5.,-6.,-7.,-8.,-9.,-10., PETSC_MAX_REAL, PETSC_MAX_REAL}; 21 PetscInt N = 6, m = 8, M, rstart, cstart, i; 22 PetscMPIInt size; 23 PetscBool loc = PETSC_FALSE; 24 PetscBool locdiag = PETSC_TRUE; 25 PetscBool localapi = PETSC_FALSE; 26 PetscBool neg = PETSC_FALSE; 27 PetscBool ismatis, ismpiaij; 28 PetscErrorCode ierr; 29 30 ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr; 31 ierr = PetscOptionsGetBool(NULL,NULL,"-neg",&neg,NULL);CHKERRQ(ierr); 32 ierr = PetscOptionsGetBool(NULL,NULL,"-loc",&loc,NULL);CHKERRQ(ierr); 33 ierr = PetscOptionsGetBool(NULL,NULL,"-locdiag",&locdiag,NULL);CHKERRQ(ierr); 34 ierr = PetscOptionsGetBool(NULL,NULL,"-localapi",&localapi,NULL);CHKERRQ(ierr); 35 ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr); 36 if (loc) { 37 if (locdiag) { 38 ierr = MatSetSizes(A,m,N,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr); 39 } else { 40 ierr = MatSetSizes(A,m,m+N,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr); 41 } 42 } else { 43 ierr = MatSetSizes(A,m,PETSC_DECIDE,PETSC_DECIDE,N);CHKERRQ(ierr); 44 } 45 ierr = MatSetFromOptions(A);CHKERRQ(ierr); 46 ierr = MatGetLayouts(A,&rmap,&cmap);CHKERRQ(ierr); 47 ierr = PetscLayoutSetUp(rmap);CHKERRQ(ierr); 48 ierr = PetscLayoutSetUp(cmap);CHKERRQ(ierr); 49 ierr = PetscLayoutGetRange(rmap,&rstart,NULL);CHKERRQ(ierr); 50 ierr = PetscLayoutGetRange(cmap,&cstart,NULL);CHKERRQ(ierr); 51 ierr = PetscLayoutGetSize(rmap,&M);CHKERRQ(ierr); 52 ierr = PetscLayoutGetSize(cmap,&N);CHKERRQ(ierr); 53 54 ierr = PetscObjectTypeCompare((PetscObject)A,MATIS,&ismatis);CHKERRQ(ierr); 55 56 /* create fake l2g maps to test the local API */ 57 ierr = ISCreateStride(PETSC_COMM_WORLD,M-rstart,rstart,1,&is);CHKERRQ(ierr); 58 ierr = ISLocalToGlobalMappingCreateIS(is,&rl2g);CHKERRQ(ierr); 59 ierr = ISDestroy(&is);CHKERRQ(ierr); 60 ierr = ISCreateStride(PETSC_COMM_WORLD,N,0,1,&is);CHKERRQ(ierr); 61 ierr = ISLocalToGlobalMappingCreateIS(is,&cl2g);CHKERRQ(ierr); 62 ierr = ISDestroy(&is);CHKERRQ(ierr); 63 ierr = MatSetLocalToGlobalMapping(A,rl2g,cl2g);CHKERRQ(ierr); 64 ierr = ISLocalToGlobalMappingDestroy(&rl2g);CHKERRQ(ierr); 65 ierr = ISLocalToGlobalMappingDestroy(&cl2g);CHKERRQ(ierr); 66 67 ierr = MatCreateVecs(A,&x,&y);CHKERRQ(ierr); 68 ierr = MatCreateVecs(A,NULL,&z);CHKERRQ(ierr); 69 ierr = VecSet(x,1.);CHKERRQ(ierr); 70 ierr = VecSet(z,2.);CHKERRQ(ierr); 71 if (!localapi) for (i = 0; i < n1; i++) i1[i] += rstart; 72 if (!localapi) for (i = 0; i < n2; i++) i2[i] += rstart; 73 if (loc) { 74 if (locdiag) { 75 for (i = 0; i < n1; i++) j1[i] += cstart; 76 for (i = 0; i < n2; i++) j2[i] += cstart; 77 } else { 78 for (i = 0; i < n1; i++) j1[i] += cstart + m; 79 for (i = 0; i < n2; i++) j2[i] += cstart + m; 80 } 81 } 82 if (neg) { n1 += 2; n2 += 2; } 83 /* MatSetPreallocationCOOLocal maps the indices! */ 84 ierr = PetscMalloc2(PetscMax(n1,n2),&it,PetscMax(n1,n2),&jt);CHKERRQ(ierr); 85 /* test with repeated entries */ 86 if (!localapi) { 87 ierr = MatSetPreallocationCOO(A,n1,i1,j1);CHKERRQ(ierr); 88 } else { 89 ierr = PetscArraycpy(it,i1,n1);CHKERRQ(ierr); 90 ierr = PetscArraycpy(jt,j1,n1);CHKERRQ(ierr); 91 ierr = MatSetPreallocationCOOLocal(A,n1,it,jt);CHKERRQ(ierr); 92 } 93 ierr = MatSetValuesCOO(A,v1,ADD_VALUES);CHKERRQ(ierr); 94 ierr = MatMult(A,x,y);CHKERRQ(ierr); 95 ierr = MyMatView(A,NULL);CHKERRQ(ierr); 96 ierr = MyVecView(y,NULL);CHKERRQ(ierr); 97 ierr = MatSetValuesCOO(A,v2,ADD_VALUES);CHKERRQ(ierr); 98 ierr = MatMultAdd(A,x,y,y);CHKERRQ(ierr); 99 ierr = MyMatView(A,NULL);CHKERRQ(ierr); 100 ierr = MyVecView(y,NULL);CHKERRQ(ierr); 101 ierr = MatTranspose(A,MAT_INITIAL_MATRIX,&At);CHKERRQ(ierr); 102 if (!ismatis) { 103 ierr = MatMatMult(A,At,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&AAt);CHKERRQ(ierr); 104 ierr = MyMatView(AAt,NULL);CHKERRQ(ierr); 105 ierr = MatDestroy(&AAt);CHKERRQ(ierr); 106 ierr = MatMatMult(At,A,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&AAt);CHKERRQ(ierr); 107 ierr = MyMatView(AAt,NULL);CHKERRQ(ierr); 108 ierr = MatDestroy(&AAt);CHKERRQ(ierr); 109 } 110 ierr = MatDestroy(&At);CHKERRQ(ierr); 111 112 /* INSERT_VALUES will overwrite matrix entries but 113 still perform the sum of the repeated entries */ 114 ierr = MatSetValuesCOO(A,v2,INSERT_VALUES);CHKERRQ(ierr); 115 ierr = MyMatView(A,NULL);CHKERRQ(ierr); 116 117 /* test with unique entries */ 118 if (!localapi) { 119 ierr = MatSetPreallocationCOO(A,n2,i2,j2);CHKERRQ(ierr); 120 } else { 121 ierr = PetscArraycpy(it,i2,n2);CHKERRQ(ierr); 122 ierr = PetscArraycpy(jt,j2,n2);CHKERRQ(ierr); 123 ierr = MatSetPreallocationCOOLocal(A,n2,it,jt);CHKERRQ(ierr); 124 } 125 ierr = MatSetValuesCOO(A,v1,ADD_VALUES);CHKERRQ(ierr); 126 ierr = MatMult(A,x,y);CHKERRQ(ierr); 127 ierr = MyMatView(A,NULL);CHKERRQ(ierr); 128 ierr = MyVecView(y,NULL);CHKERRQ(ierr); 129 ierr = MatSetValuesCOO(A,v2,ADD_VALUES);CHKERRQ(ierr); 130 ierr = MatMultAdd(A,x,y,z);CHKERRQ(ierr); 131 ierr = MyMatView(A,NULL);CHKERRQ(ierr); 132 ierr = MyVecView(z,NULL);CHKERRQ(ierr); 133 if (!localapi) { 134 ierr = MatSetPreallocationCOO(A,n2,i2,j2);CHKERRQ(ierr); 135 } else { 136 ierr = PetscArraycpy(it,i2,n2);CHKERRQ(ierr); 137 ierr = PetscArraycpy(jt,j2,n2);CHKERRQ(ierr); 138 ierr = MatSetPreallocationCOOLocal(A,n2,it,jt);CHKERRQ(ierr); 139 } 140 ierr = MatSetValuesCOO(A,v1,INSERT_VALUES);CHKERRQ(ierr); 141 ierr = MatMult(A,x,y);CHKERRQ(ierr); 142 ierr = MyMatView(A,NULL);CHKERRQ(ierr); 143 ierr = MyVecView(y,NULL);CHKERRQ(ierr); 144 ierr = MatSetValuesCOO(A,v2,INSERT_VALUES);CHKERRQ(ierr); 145 ierr = MatMultAdd(A,x,y,z);CHKERRQ(ierr); 146 ierr = MyMatView(A,NULL);CHKERRQ(ierr); 147 ierr = MyVecView(z,NULL);CHKERRQ(ierr); 148 ierr = MatTranspose(A,MAT_INITIAL_MATRIX,&At);CHKERRQ(ierr); 149 if (!ismatis) { 150 ierr = MatMatMult(A,At,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&AAt);CHKERRQ(ierr); 151 ierr = MyMatView(AAt,NULL);CHKERRQ(ierr); 152 ierr = MatDestroy(&AAt);CHKERRQ(ierr); 153 ierr = MatMatMult(At,A,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&AAt);CHKERRQ(ierr); 154 ierr = MyMatView(AAt,NULL);CHKERRQ(ierr); 155 ierr = MatDestroy(&AAt);CHKERRQ(ierr); 156 } 157 ierr = MatDestroy(&At);CHKERRQ(ierr); 158 159 /* test providing diagonal first, then offdiagonal */ 160 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)A),&size);CHKERRMPI(ierr); 161 ierr = PetscObjectBaseTypeCompare((PetscObject)A,MATMPIAIJ,&ismpiaij);CHKERRQ(ierr); 162 if (ismpiaij && size > 1) { 163 Mat lA,lB; 164 const PetscInt *garray,*iA,*jA,*iB,*jB; 165 const PetscScalar *vA,*vB; 166 PetscScalar *coo_v; 167 PetscInt *coo_i,*coo_j; 168 PetscInt i,j,nA,nB,nnz; 169 PetscBool flg; 170 171 ierr = MatMPIAIJGetSeqAIJ(A,&lA,&lB,&garray);CHKERRQ(ierr); 172 ierr = MatSeqAIJGetArrayRead(lA,&vA);CHKERRQ(ierr); 173 ierr = MatSeqAIJGetArrayRead(lB,&vB);CHKERRQ(ierr); 174 ierr = MatGetRowIJ(lA,0,PETSC_FALSE,PETSC_FALSE,&nA,&iA,&jA,&flg);CHKERRQ(ierr); 175 ierr = MatGetRowIJ(lB,0,PETSC_FALSE,PETSC_FALSE,&nB,&iB,&jB,&flg);CHKERRQ(ierr); 176 nnz = iA[nA] + iB[nB]; 177 ierr = PetscMalloc3(nnz,&coo_i,nnz,&coo_j,nnz,&coo_v);CHKERRQ(ierr); 178 nnz = 0; 179 for (i=0;i<nA;i++) { 180 for (j=iA[i];j<iA[i+1];j++,nnz++) { 181 coo_i[nnz] = i+rstart; 182 coo_j[nnz] = jA[j]+cstart; 183 coo_v[nnz] = vA[j]; 184 } 185 } 186 for (i=0;i<nB;i++) { 187 for (j=iB[i];j<iB[i+1];j++,nnz++) { 188 coo_i[nnz] = i+rstart; 189 coo_j[nnz] = garray[jB[j]]; 190 coo_v[nnz] = vB[j]; 191 } 192 } 193 ierr = MatRestoreRowIJ(lA,0,PETSC_FALSE,PETSC_FALSE,&nA,&iA,&jA,&flg);CHKERRQ(ierr); 194 ierr = MatRestoreRowIJ(lB,0,PETSC_FALSE,PETSC_FALSE,&nB,&iB,&jB,&flg);CHKERRQ(ierr); 195 ierr = MatSeqAIJRestoreArrayRead(lA,&vA);CHKERRQ(ierr); 196 ierr = MatSeqAIJRestoreArrayRead(lB,&vB);CHKERRQ(ierr); 197 198 ierr = MatSetPreallocationCOO(A,nnz,coo_i,coo_j);CHKERRQ(ierr); 199 ierr = MatSetValuesCOO(A,coo_v,ADD_VALUES);CHKERRQ(ierr); 200 ierr = MatMult(A,x,y);CHKERRQ(ierr); 201 ierr = MyMatView(A,NULL);CHKERRQ(ierr); 202 ierr = MyVecView(y,NULL);CHKERRQ(ierr); 203 ierr = MatSetValuesCOO(A,coo_v,INSERT_VALUES);CHKERRQ(ierr); 204 ierr = MatMult(A,x,y);CHKERRQ(ierr); 205 ierr = MyMatView(A,NULL);CHKERRQ(ierr); 206 ierr = MyVecView(y,NULL);CHKERRQ(ierr); 207 ierr = MatTranspose(A,MAT_INITIAL_MATRIX,&At);CHKERRQ(ierr); 208 ierr = MatMatMult(A,At,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&AAt);CHKERRQ(ierr); 209 ierr = MyMatView(AAt,NULL);CHKERRQ(ierr); 210 ierr = MatDestroy(&AAt);CHKERRQ(ierr); 211 ierr = MatMatMult(At,A,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&AAt);CHKERRQ(ierr); 212 ierr = MyMatView(AAt,NULL);CHKERRQ(ierr); 213 ierr = MatDestroy(&AAt);CHKERRQ(ierr); 214 ierr = MatDestroy(&At);CHKERRQ(ierr); 215 216 ierr = PetscFree3(coo_i,coo_j,coo_v);CHKERRQ(ierr); 217 } 218 ierr = PetscFree2(it,jt);CHKERRQ(ierr); 219 ierr = VecDestroy(&z);CHKERRQ(ierr); 220 ierr = VecDestroy(&x);CHKERRQ(ierr); 221 ierr = VecDestroy(&y);CHKERRQ(ierr); 222 ierr = MatDestroy(&A);CHKERRQ(ierr); 223 ierr = PetscFinalize(); 224 return ierr; 225 } 226 227 /*TEST 228 229 test: 230 suffix: 1 231 filter: grep -v type 232 diff_args: -j 233 args: -mat_type {{seqaij mpiaij}} -localapi {{0 1}} -neg {{0 1}} 234 235 test: 236 requires: cuda 237 suffix: 1_cuda 238 filter: grep -v type 239 diff_args: -j 240 args: -mat_type {{seqaijcusparse mpiaijcusparse}} -localapi {{0 1}} -neg {{0 1}} 241 output_file: output/ex123_1.out 242 243 test: 244 requires: kokkos_kernels !sycl 245 suffix: 1_kokkos 246 filter: grep -v type 247 diff_args: -j 248 args: -mat_type {{seqaijkokkos mpiaijkokkos}} -localapi {{0 1}} -neg {{0 1}} 249 output_file: output/ex123_1.out 250 251 test: 252 suffix: 2 253 nsize: 7 254 filter: grep -v type 255 diff_args: -j 256 args: -mat_type mpiaij -localapi {{0 1}} -neg {{0 1}} 257 258 test: 259 requires: cuda 260 suffix: 2_cuda 261 nsize: 7 262 filter: grep -v type 263 diff_args: -j 264 args: -mat_type mpiaijcusparse -localapi {{0 1}} -neg {{0 1}} 265 output_file: output/ex123_2.out 266 267 test: 268 requires: kokkos_kernels !sycl 269 suffix: 2_kokkos 270 nsize: 7 271 filter: grep -v type 272 diff_args: -j 273 args: -mat_type mpiaijkokkos -localapi {{0 1}} -neg {{0 1}} 274 output_file: output/ex123_2.out 275 276 test: 277 suffix: 3 278 nsize: 3 279 filter: grep -v type 280 diff_args: -j 281 args: -mat_type mpiaij -loc -localapi {{0 1}} -neg {{0 1}} 282 283 test: 284 requires: cuda 285 suffix: 3_cuda 286 nsize: 3 287 filter: grep -v type 288 diff_args: -j 289 args: -mat_type mpiaijcusparse -loc -localapi {{0 1}} -neg {{0 1}} 290 output_file: output/ex123_3.out 291 292 test: 293 requires: !sycl kokkos_kernels 294 suffix: 3_kokkos 295 nsize: 3 296 filter: grep -v type 297 diff_args: -j 298 args: -mat_type aijkokkos -loc -localapi {{0 1}} -neg {{0 1}} 299 output_file: output/ex123_3.out 300 301 test: 302 suffix: 4 303 nsize: 4 304 filter: grep -v type 305 diff_args: -j 306 args: -mat_type mpiaij -loc -locdiag 0 -localapi {{0 1}} -neg {{0 1}} 307 308 test: 309 requires: cuda 310 suffix: 4_cuda 311 nsize: 4 312 filter: grep -v type 313 diff_args: -j 314 args: -mat_type mpiaijcusparse -loc -locdiag 0 -localapi {{0 1}} -neg {{0 1}} 315 output_file: output/ex123_4.out 316 317 test: 318 requires: !sycl kokkos_kernels 319 suffix: 4_kokkos 320 nsize: 4 321 filter: grep -v type 322 diff_args: -j 323 args: -mat_type aijkokkos -loc -locdiag 0 -localapi {{0 1}} -neg {{0 1}} 324 output_file: output/ex123_4.out 325 326 test: 327 suffix: matis 328 nsize: 3 329 filter: grep -v type 330 diff_args: -j 331 args: -mat_type is -localapi {{0 1}} -neg {{0 1}} 332 333 TEST*/ 334