1 2 static char help[] = "Tests the use of interface functions for MATIS matrices and conversion routines.\n"; 3 4 #include <petscmat.h> 5 6 PetscErrorCode TestMatZeroRows(Mat,Mat,PetscBool,IS,PetscScalar); 7 PetscErrorCode CheckMat(Mat,Mat,PetscBool,const char*); 8 9 int main(int argc,char **args) 10 { 11 Mat A,B,A2,B2,T; 12 Mat Aee,Aeo,Aoe,Aoo; 13 Mat *mats,*Asub,*Bsub; 14 Vec x,y; 15 MatInfo info; 16 ISLocalToGlobalMapping cmap,rmap; 17 IS is,is2,reven,rodd,ceven,codd; 18 IS *rows,*cols; 19 IS irow[2], icol[2]; 20 PetscLayout rlayout, clayout; 21 const PetscInt *rrange, *crange; 22 MatType lmtype; 23 PetscScalar diag = 2.; 24 PetscInt n,m,i,lm,ln; 25 PetscInt rst,ren,cst,cen,nr,nc; 26 PetscMPIInt rank,size,lrank,rrank; 27 PetscBool testT,squaretest,isaij; 28 PetscBool permute = PETSC_FALSE, negmap = PETSC_FALSE, repmap = PETSC_FALSE; 29 PetscBool diffmap = PETSC_TRUE, symmetric = PETSC_FALSE, issymmetric; 30 31 PetscCall(PetscInitialize(&argc,&args,(char*)0,help)); 32 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank)); 33 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 34 m = n = 2*size; 35 PetscCall(PetscOptionsGetBool(NULL,NULL,"-symmetric",&symmetric,NULL)); 36 PetscCall(PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL)); 37 PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL)); 38 PetscCall(PetscOptionsGetBool(NULL,NULL,"-negmap",&negmap,NULL)); 39 PetscCall(PetscOptionsGetBool(NULL,NULL,"-repmap",&repmap,NULL)); 40 PetscCall(PetscOptionsGetBool(NULL,NULL,"-permmap",&permute,NULL)); 41 PetscCall(PetscOptionsGetBool(NULL,NULL,"-diffmap",&diffmap,NULL)); 42 PetscCheck(size == 1 || m >= 4,PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"Number of rows should be larger or equal 4 for parallel runs"); 43 PetscCheck(size != 1 || m >= 2,PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"Number of rows should be larger or equal 2 for uniprocessor runs"); 44 PetscCheck(n >= 2,PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"Number of cols should be larger or equal 2"); 45 46 /* create a MATIS matrix */ 47 PetscCall(MatCreate(PETSC_COMM_WORLD,&A)); 48 PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m,n)); 49 PetscCall(MatSetType(A,MATIS)); 50 PetscCall(MatSetFromOptions(A)); 51 if (!negmap && !repmap) { 52 /* This is not the proper setting for MATIS for finite elements, it is just used to test the routines 53 Here we use a one-to-one correspondence between local row/column spaces and global row/column spaces 54 Equivalent to passing NULL for the mapping */ 55 PetscCall(ISCreateStride(PETSC_COMM_WORLD,n,0,1,&is)); 56 } else if (negmap && !repmap) { /* non repeated but with negative indices */ 57 PetscCall(ISCreateStride(PETSC_COMM_WORLD,n+2,-2,1,&is)); 58 } else if (!negmap && repmap) { /* non negative but repeated indices */ 59 IS isl[2]; 60 61 PetscCall(ISCreateStride(PETSC_COMM_WORLD,n,0,1,&isl[0])); 62 PetscCall(ISCreateStride(PETSC_COMM_WORLD,n,n-1,-1,&isl[1])); 63 PetscCall(ISConcatenate(PETSC_COMM_WORLD,2,isl,&is)); 64 PetscCall(ISDestroy(&isl[0])); 65 PetscCall(ISDestroy(&isl[1])); 66 } else { /* negative and repeated indices */ 67 IS isl[2]; 68 69 PetscCall(ISCreateStride(PETSC_COMM_WORLD,n+1,-1,1,&isl[0])); 70 PetscCall(ISCreateStride(PETSC_COMM_WORLD,n+1,n-1,-1,&isl[1])); 71 PetscCall(ISConcatenate(PETSC_COMM_WORLD,2,isl,&is)); 72 PetscCall(ISDestroy(&isl[0])); 73 PetscCall(ISDestroy(&isl[1])); 74 } 75 PetscCall(ISLocalToGlobalMappingCreateIS(is,&cmap)); 76 PetscCall(ISDestroy(&is)); 77 78 if (m != n || diffmap) { 79 PetscCall(ISCreateStride(PETSC_COMM_WORLD,m,permute ? m-1 : 0,permute ? -1 : 1,&is)); 80 PetscCall(ISLocalToGlobalMappingCreateIS(is,&rmap)); 81 PetscCall(ISDestroy(&is)); 82 } else { 83 PetscCall(PetscObjectReference((PetscObject)cmap)); 84 rmap = cmap; 85 } 86 87 PetscCall(MatSetLocalToGlobalMapping(A,rmap,cmap)); 88 PetscCall(MatISStoreL2L(A,PETSC_FALSE)); 89 PetscCall(MatISSetPreallocation(A,3,NULL,3,NULL)); 90 PetscCall(MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,(PetscBool)!(repmap || negmap))); /* I do not want to precompute the pattern */ 91 PetscCall(ISLocalToGlobalMappingGetSize(rmap,&lm)); 92 PetscCall(ISLocalToGlobalMappingGetSize(cmap,&ln)); 93 for (i=0; i<lm; i++) { 94 PetscScalar v[3]; 95 PetscInt cols[3]; 96 97 cols[0] = (i-1+n)%n; 98 cols[1] = i%n; 99 cols[2] = (i+1)%n; 100 v[0] = -1.*(symmetric ? PetscMin(i+1,cols[0]+1) : i+1); 101 v[1] = 2.*(symmetric ? PetscMin(i+1,cols[1]+1) : i+1); 102 v[2] = -1.*(symmetric ? PetscMin(i+1,cols[2]+1) : i+1); 103 PetscCall(ISGlobalToLocalMappingApply(cmap,IS_GTOLM_MASK,3,cols,NULL,cols)); 104 PetscCall(MatSetValuesLocal(A,1,&i,3,cols,v,ADD_VALUES)); 105 } 106 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 107 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 108 109 /* activate tests for square matrices with same maps only */ 110 PetscCall(MatHasCongruentLayouts(A,&squaretest)); 111 if (squaretest && rmap != cmap) { 112 PetscInt nr, nc; 113 114 PetscCall(ISLocalToGlobalMappingGetSize(rmap,&nr)); 115 PetscCall(ISLocalToGlobalMappingGetSize(cmap,&nc)); 116 if (nr != nc) squaretest = PETSC_FALSE; 117 else { 118 const PetscInt *idxs1,*idxs2; 119 120 PetscCall(ISLocalToGlobalMappingGetIndices(rmap,&idxs1)); 121 PetscCall(ISLocalToGlobalMappingGetIndices(cmap,&idxs2)); 122 PetscCall(PetscArraycmp(idxs1,idxs2,nr,&squaretest)); 123 PetscCall(ISLocalToGlobalMappingRestoreIndices(rmap,&idxs1)); 124 PetscCall(ISLocalToGlobalMappingRestoreIndices(cmap,&idxs2)); 125 } 126 PetscCall(MPIU_Allreduce(MPI_IN_PLACE,&squaretest,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A))); 127 } 128 129 /* test MatISGetLocalMat */ 130 PetscCall(MatISGetLocalMat(A,&B)); 131 PetscCall(MatGetType(B,&lmtype)); 132 133 /* test MatGetInfo */ 134 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Test MatGetInfo\n")); 135 PetscCall(MatGetInfo(A,MAT_LOCAL,&info)); 136 PetscCall(PetscViewerASCIIPushSynchronized(PETSC_VIEWER_STDOUT_WORLD)); 137 PetscCall(PetscViewerASCIISynchronizedPrintf(PETSC_VIEWER_STDOUT_WORLD,"Process %2d: %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n",PetscGlobalRank,(PetscInt)info.nz_used, 138 (PetscInt)info.nz_allocated,(PetscInt)info.nz_unneeded,(PetscInt)info.assemblies,(PetscInt)info.mallocs)); 139 PetscCall(PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD)); 140 PetscCall(MatGetInfo(A,MAT_GLOBAL_MAX,&info)); 141 PetscCall(PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_WORLD,"GlobalMax : %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n",(PetscInt)info.nz_used, 142 (PetscInt)info.nz_allocated,(PetscInt)info.nz_unneeded,(PetscInt)info.assemblies,(PetscInt)info.mallocs)); 143 PetscCall(MatGetInfo(A,MAT_GLOBAL_SUM,&info)); 144 PetscCall(PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_WORLD,"GlobalSum : %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n",(PetscInt)info.nz_used, 145 (PetscInt)info.nz_allocated,(PetscInt)info.nz_unneeded,(PetscInt)info.assemblies,(PetscInt)info.mallocs)); 146 147 /* test MatIsSymmetric */ 148 PetscCall(MatIsSymmetric(A,0.0,&issymmetric)); 149 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Test MatIsSymmetric: %d\n",issymmetric)); 150 151 /* Create a MPIAIJ matrix, same as A */ 152 PetscCall(MatCreate(PETSC_COMM_WORLD,&B)); 153 PetscCall(MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,m,n)); 154 PetscCall(MatSetType(B,MATAIJ)); 155 PetscCall(MatSetFromOptions(B)); 156 PetscCall(MatSetUp(B)); 157 PetscCall(MatSetLocalToGlobalMapping(B,rmap,cmap)); 158 PetscCall(MatMPIAIJSetPreallocation(B,3,NULL,3,NULL)); 159 PetscCall(MatMPIBAIJSetPreallocation(B,1,3,NULL,3,NULL)); 160 #if defined(PETSC_HAVE_HYPRE) 161 PetscCall(MatHYPRESetPreallocation(B,3,NULL,3,NULL)); 162 #endif 163 PetscCall(MatISSetPreallocation(B,3,NULL,3,NULL)); 164 PetscCall(MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,(PetscBool)!(repmap || negmap))); /* I do not want to precompute the pattern */ 165 for (i=0; i<lm; i++) { 166 PetscScalar v[3]; 167 PetscInt cols[3]; 168 169 cols[0] = (i-1+n)%n; 170 cols[1] = i%n; 171 cols[2] = (i+1)%n; 172 v[0] = -1.*(symmetric ? PetscMin(i+1,cols[0]+1) : i+1); 173 v[1] = 2.*(symmetric ? PetscMin(i+1,cols[1]+1) : i+1); 174 v[2] = -1.*(symmetric ? PetscMin(i+1,cols[2]+1) : i+1); 175 PetscCall(ISGlobalToLocalMappingApply(cmap,IS_GTOLM_MASK,3,cols,NULL,cols)); 176 PetscCall(MatSetValuesLocal(B,1,&i,3,cols,v,ADD_VALUES)); 177 } 178 PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 179 PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 180 181 /* test MatView */ 182 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Test MatView\n")); 183 PetscCall(MatView(A,NULL)); 184 PetscCall(MatView(B,NULL)); 185 186 /* test CheckMat */ 187 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Test CheckMat\n")); 188 PetscCall(CheckMat(A,B,PETSC_FALSE,"CheckMat")); 189 190 /* test MatDuplicate and MatAXPY */ 191 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Test MatDuplicate and MatAXPY\n")); 192 PetscCall(MatDuplicate(A,MAT_COPY_VALUES,&A2)); 193 PetscCall(CheckMat(A,A2,PETSC_FALSE,"MatDuplicate and MatAXPY")); 194 195 /* test MatConvert */ 196 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Test MatConvert_IS_XAIJ\n")); 197 PetscCall(MatConvert(A2,MATAIJ,MAT_INITIAL_MATRIX,&B2)); 198 PetscCall(CheckMat(B,B2,PETSC_TRUE,"MatConvert_IS_XAIJ MAT_INITIAL_MATRIX")); 199 PetscCall(MatConvert(A2,MATAIJ,MAT_REUSE_MATRIX,&B2)); 200 PetscCall(CheckMat(B,B2,PETSC_TRUE,"MatConvert_IS_XAIJ MAT_REUSE_MATRIX")); 201 PetscCall(MatConvert(A2,MATAIJ,MAT_INPLACE_MATRIX,&A2)); 202 PetscCall(CheckMat(B,A2,PETSC_TRUE,"MatConvert_IS_XAIJ MAT_INPLACE_MATRIX")); 203 PetscCall(MatDestroy(&A2)); 204 PetscCall(MatDestroy(&B2)); 205 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Test MatConvert_XAIJ_IS\n")); 206 PetscCall(MatDuplicate(B,MAT_COPY_VALUES,&B2)); 207 PetscCall(MatConvert(B2,MATIS,MAT_INITIAL_MATRIX,&A2)); 208 PetscCall(CheckMat(A,A2,PETSC_TRUE,"MatConvert_XAIJ_IS MAT_INITIAL_MATRIX")); 209 PetscCall(MatConvert(B2,MATIS,MAT_REUSE_MATRIX,&A2)); 210 PetscCall(CheckMat(A,A2,PETSC_TRUE,"MatConvert_XAIJ_IS MAT_REUSE_MATRIX")); 211 PetscCall(MatConvert(B2,MATIS,MAT_INPLACE_MATRIX,&B2)); 212 PetscCall(CheckMat(A,B2,PETSC_TRUE,"MatConvert_XAIJ_IS MAT_INPLACE_MATRIX")); 213 PetscCall(MatDestroy(&A2)); 214 PetscCall(MatDestroy(&B2)); 215 PetscCall(PetscStrcmp(lmtype,MATSEQAIJ,&isaij)); 216 if (size == 1 && isaij) { /* tests special code paths in MatConvert_IS_XAIJ */ 217 PetscInt ri, ci, rr[3] = {0,1,0}, cr[4] = {1,2,0,1}, rk[3] = {0,2,1}, ck[4] = {1,0,3,2}; 218 ISLocalToGlobalMapping tcmap,trmap; 219 220 for (ri = 0; ri < 2; ri++) { 221 PetscInt *r; 222 223 r = (PetscInt*)(ri == 0 ? rr : rk); 224 for (ci = 0; ci < 2; ci++) { 225 PetscInt *c,rb,cb; 226 227 c = (PetscInt*)(ci == 0 ? cr : ck); 228 for (rb = 1; rb < 4; rb++) { 229 PetscCall(ISCreateBlock(PETSC_COMM_SELF,rb,3,r,PETSC_COPY_VALUES,&is)); 230 PetscCall(ISLocalToGlobalMappingCreateIS(is,&trmap)); 231 PetscCall(ISDestroy(&is)); 232 for (cb = 1; cb < 4; cb++) { 233 Mat T,lT,T2; 234 char testname[256]; 235 236 PetscCall(PetscSNPrintf(testname,sizeof(testname),"MatConvert_IS_XAIJ special case (%" PetscInt_FMT " %" PetscInt_FMT ", bs %" PetscInt_FMT " %" PetscInt_FMT ")",ri,ci,rb,cb)); 237 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Test %s\n",testname)); 238 239 PetscCall(ISCreateBlock(PETSC_COMM_SELF,cb,4,c,PETSC_COPY_VALUES,&is)); 240 PetscCall(ISLocalToGlobalMappingCreateIS(is,&tcmap)); 241 PetscCall(ISDestroy(&is)); 242 243 PetscCall(MatCreate(PETSC_COMM_SELF,&T)); 244 PetscCall(MatSetSizes(T,PETSC_DECIDE,PETSC_DECIDE,rb*3,cb*4)); 245 PetscCall(MatSetType(T,MATIS)); 246 PetscCall(MatSetLocalToGlobalMapping(T,trmap,tcmap)); 247 PetscCall(ISLocalToGlobalMappingDestroy(&tcmap)); 248 PetscCall(MatISGetLocalMat(T,&lT)); 249 PetscCall(MatSetType(lT,MATSEQAIJ)); 250 PetscCall(MatSeqAIJSetPreallocation(lT,cb*4,NULL)); 251 PetscCall(MatSetRandom(lT,NULL)); 252 PetscCall(MatConvert(lT,lmtype,MAT_INPLACE_MATRIX,&lT)); 253 PetscCall(MatISRestoreLocalMat(T,&lT)); 254 PetscCall(MatAssemblyBegin(T,MAT_FINAL_ASSEMBLY)); 255 PetscCall(MatAssemblyEnd(T,MAT_FINAL_ASSEMBLY)); 256 257 PetscCall(MatConvert(T,MATAIJ,MAT_INITIAL_MATRIX,&T2)); 258 PetscCall(CheckMat(T,T2,PETSC_TRUE,"MAT_INITIAL_MATRIX")); 259 PetscCall(MatConvert(T,MATAIJ,MAT_REUSE_MATRIX,&T2)); 260 PetscCall(CheckMat(T,T2,PETSC_TRUE,"MAT_REUSE_MATRIX")); 261 PetscCall(MatDestroy(&T2)); 262 PetscCall(MatDuplicate(T,MAT_COPY_VALUES,&T2)); 263 PetscCall(MatConvert(T2,MATAIJ,MAT_INPLACE_MATRIX,&T2)); 264 PetscCall(CheckMat(T,T2,PETSC_TRUE,"MAT_INPLACE_MATRIX")); 265 PetscCall(MatDestroy(&T)); 266 PetscCall(MatDestroy(&T2)); 267 } 268 PetscCall(ISLocalToGlobalMappingDestroy(&trmap)); 269 } 270 } 271 } 272 } 273 274 /* test MatDiagonalScale */ 275 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Test MatDiagonalScale\n")); 276 PetscCall(MatDuplicate(A,MAT_COPY_VALUES,&A2)); 277 PetscCall(MatDuplicate(B,MAT_COPY_VALUES,&B2)); 278 PetscCall(MatCreateVecs(A,&x,&y)); 279 PetscCall(VecSetRandom(x,NULL)); 280 if (issymmetric) { 281 PetscCall(VecCopy(x,y)); 282 } else { 283 PetscCall(VecSetRandom(y,NULL)); 284 PetscCall(VecScale(y,8.)); 285 } 286 PetscCall(MatDiagonalScale(A2,y,x)); 287 PetscCall(MatDiagonalScale(B2,y,x)); 288 PetscCall(CheckMat(A2,B2,PETSC_FALSE,"MatDiagonalScale")); 289 PetscCall(MatDestroy(&A2)); 290 PetscCall(MatDestroy(&B2)); 291 PetscCall(VecDestroy(&x)); 292 PetscCall(VecDestroy(&y)); 293 294 /* test MatPtAP (A IS and B AIJ) */ 295 if (isaij && m == n) { 296 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Test MatPtAP\n")); 297 PetscCall(MatISStoreL2L(A,PETSC_TRUE)); 298 PetscCall(MatPtAP(A,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&A2)); 299 PetscCall(MatPtAP(B,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&B2)); 300 PetscCall(CheckMat(A2,B2,PETSC_FALSE,"MatPtAP MAT_INITIAL_MATRIX")); 301 PetscCall(MatPtAP(A,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&A2)); 302 PetscCall(CheckMat(A2,B2,PETSC_FALSE,"MatPtAP MAT_REUSE_MATRIX")); 303 PetscCall(MatDestroy(&A2)); 304 PetscCall(MatDestroy(&B2)); 305 } 306 307 /* test MatGetLocalSubMatrix */ 308 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Test MatGetLocalSubMatrix\n")); 309 PetscCall(MatDuplicate(A,MAT_DO_NOT_COPY_VALUES,&A2)); 310 PetscCall(ISCreateStride(PETSC_COMM_SELF,lm/2+lm%2,0,2,&reven)); 311 PetscCall(ISComplement(reven,0,lm,&rodd)); 312 PetscCall(ISCreateStride(PETSC_COMM_SELF,ln/2+ln%2,0,2,&ceven)); 313 PetscCall(ISComplement(ceven,0,ln,&codd)); 314 PetscCall(MatGetLocalSubMatrix(A2,reven,ceven,&Aee)); 315 PetscCall(MatGetLocalSubMatrix(A2,reven,codd,&Aeo)); 316 PetscCall(MatGetLocalSubMatrix(A2,rodd,ceven,&Aoe)); 317 PetscCall(MatGetLocalSubMatrix(A2,rodd,codd,&Aoo)); 318 for (i=0; i<lm; i++) { 319 PetscInt j,je,jo,colse[3], colso[3]; 320 PetscScalar ve[3], vo[3]; 321 PetscScalar v[3]; 322 PetscInt cols[3]; 323 PetscInt row = i/2; 324 325 cols[0] = (i-1+n)%n; 326 cols[1] = i%n; 327 cols[2] = (i+1)%n; 328 v[0] = -1.*(symmetric ? PetscMin(i+1,cols[0]+1) : i+1); 329 v[1] = 2.*(symmetric ? PetscMin(i+1,cols[1]+1) : i+1); 330 v[2] = -1.*(symmetric ? PetscMin(i+1,cols[2]+1) : i+1); 331 PetscCall(ISGlobalToLocalMappingApply(cmap,IS_GTOLM_MASK,3,cols,NULL,cols)); 332 for (j=0,je=0,jo=0;j<3;j++) { 333 if (cols[j]%2) { 334 vo[jo] = v[j]; 335 colso[jo++] = cols[j]/2; 336 } else { 337 ve[je] = v[j]; 338 colse[je++] = cols[j]/2; 339 } 340 } 341 if (i%2) { 342 PetscCall(MatSetValuesLocal(Aoe,1,&row,je,colse,ve,ADD_VALUES)); 343 PetscCall(MatSetValuesLocal(Aoo,1,&row,jo,colso,vo,ADD_VALUES)); 344 } else { 345 PetscCall(MatSetValuesLocal(Aee,1,&row,je,colse,ve,ADD_VALUES)); 346 PetscCall(MatSetValuesLocal(Aeo,1,&row,jo,colso,vo,ADD_VALUES)); 347 } 348 } 349 PetscCall(MatRestoreLocalSubMatrix(A2,reven,ceven,&Aee)); 350 PetscCall(MatRestoreLocalSubMatrix(A2,reven,codd,&Aeo)); 351 PetscCall(MatRestoreLocalSubMatrix(A2,rodd,ceven,&Aoe)); 352 PetscCall(MatRestoreLocalSubMatrix(A2,rodd,codd,&Aoo)); 353 PetscCall(ISDestroy(&reven)); 354 PetscCall(ISDestroy(&ceven)); 355 PetscCall(ISDestroy(&rodd)); 356 PetscCall(ISDestroy(&codd)); 357 PetscCall(MatAssemblyBegin(A2,MAT_FINAL_ASSEMBLY)); 358 PetscCall(MatAssemblyEnd(A2,MAT_FINAL_ASSEMBLY)); 359 PetscCall(MatAXPY(A2,-1.,A,SAME_NONZERO_PATTERN)); 360 PetscCall(CheckMat(A2,NULL,PETSC_FALSE,"MatGetLocalSubMatrix")); 361 PetscCall(MatDestroy(&A2)); 362 363 /* test MatConvert_Nest_IS */ 364 testT = PETSC_FALSE; 365 PetscCall(PetscOptionsGetBool(NULL,NULL,"-test_trans",&testT,NULL)); 366 367 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Test MatConvert_Nest_IS\n")); 368 nr = 2; 369 nc = 2; 370 PetscCall(PetscOptionsGetInt(NULL,NULL,"-nr",&nr,NULL)); 371 PetscCall(PetscOptionsGetInt(NULL,NULL,"-nc",&nc,NULL)); 372 if (testT) { 373 PetscCall(MatGetOwnershipRange(A,&cst,&cen)); 374 PetscCall(MatGetOwnershipRangeColumn(A,&rst,&ren)); 375 } else { 376 PetscCall(MatGetOwnershipRange(A,&rst,&ren)); 377 PetscCall(MatGetOwnershipRangeColumn(A,&cst,&cen)); 378 } 379 PetscCall(PetscMalloc3(nr,&rows,nc,&cols,2*nr*nc,&mats)); 380 for (i=0;i<nr*nc;i++) { 381 if (testT) { 382 PetscCall(MatCreateTranspose(A,&mats[i])); 383 PetscCall(MatTranspose(B,MAT_INITIAL_MATRIX,&mats[i+nr*nc])); 384 } else { 385 PetscCall(MatDuplicate(A,MAT_COPY_VALUES,&mats[i])); 386 PetscCall(MatDuplicate(B,MAT_COPY_VALUES,&mats[i+nr*nc])); 387 } 388 } 389 for (i=0;i<nr;i++) { 390 PetscCall(ISCreateStride(PETSC_COMM_WORLD,ren-rst,i+rst,nr,&rows[i])); 391 } 392 for (i=0;i<nc;i++) { 393 PetscCall(ISCreateStride(PETSC_COMM_WORLD,cen-cst,i+cst,nc,&cols[i])); 394 } 395 PetscCall(MatCreateNest(PETSC_COMM_WORLD,nr,rows,nc,cols,mats,&A2)); 396 PetscCall(MatCreateNest(PETSC_COMM_WORLD,nr,rows,nc,cols,mats+nr*nc,&B2)); 397 for (i=0;i<nr;i++) { 398 PetscCall(ISDestroy(&rows[i])); 399 } 400 for (i=0;i<nc;i++) { 401 PetscCall(ISDestroy(&cols[i])); 402 } 403 for (i=0;i<2*nr*nc;i++) { 404 PetscCall(MatDestroy(&mats[i])); 405 } 406 PetscCall(PetscFree3(rows,cols,mats)); 407 PetscCall(MatConvert(B2,MATAIJ,MAT_INITIAL_MATRIX,&T)); 408 PetscCall(MatDestroy(&B2)); 409 PetscCall(MatConvert(A2,MATIS,MAT_INITIAL_MATRIX,&B2)); 410 PetscCall(CheckMat(B2,T,PETSC_TRUE,"MatConvert_Nest_IS MAT_INITIAL_MATRIX")); 411 PetscCall(MatConvert(A2,MATIS,MAT_REUSE_MATRIX,&B2)); 412 PetscCall(CheckMat(B2,T,PETSC_TRUE,"MatConvert_Nest_IS MAT_REUSE_MATRIX")); 413 PetscCall(MatDestroy(&B2)); 414 PetscCall(MatConvert(A2,MATIS,MAT_INPLACE_MATRIX,&A2)); 415 PetscCall(CheckMat(A2,T,PETSC_TRUE,"MatConvert_Nest_IS MAT_INPLACE_MATRIX")); 416 PetscCall(MatDestroy(&T)); 417 PetscCall(MatDestroy(&A2)); 418 419 /* test MatCreateSubMatrix */ 420 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Test MatCreateSubMatrix\n")); 421 if (rank == 0) { 422 PetscCall(ISCreateStride(PETSC_COMM_WORLD,1,1,1,&is)); 423 PetscCall(ISCreateStride(PETSC_COMM_WORLD,2,0,1,&is2)); 424 } else if (rank == 1) { 425 PetscCall(ISCreateStride(PETSC_COMM_WORLD,1,0,1,&is)); 426 if (n > 3) { 427 PetscCall(ISCreateStride(PETSC_COMM_WORLD,1,3,1,&is2)); 428 } else { 429 PetscCall(ISCreateStride(PETSC_COMM_WORLD,0,0,1,&is2)); 430 } 431 } else if (rank == 2 && n > 4) { 432 PetscCall(ISCreateStride(PETSC_COMM_WORLD,0,0,1,&is)); 433 PetscCall(ISCreateStride(PETSC_COMM_WORLD,n-4,4,1,&is2)); 434 } else { 435 PetscCall(ISCreateStride(PETSC_COMM_WORLD,0,0,1,&is)); 436 PetscCall(ISCreateStride(PETSC_COMM_WORLD,0,0,1,&is2)); 437 } 438 PetscCall(MatCreateSubMatrix(A,is,is,MAT_INITIAL_MATRIX,&A2)); 439 PetscCall(MatCreateSubMatrix(B,is,is,MAT_INITIAL_MATRIX,&B2)); 440 PetscCall(CheckMat(A2,B2,PETSC_TRUE,"first MatCreateSubMatrix")); 441 442 PetscCall(MatCreateSubMatrix(A,is,is,MAT_REUSE_MATRIX,&A2)); 443 PetscCall(MatCreateSubMatrix(B,is,is,MAT_REUSE_MATRIX,&B2)); 444 PetscCall(CheckMat(A2,B2,PETSC_FALSE,"reuse MatCreateSubMatrix")); 445 PetscCall(MatDestroy(&A2)); 446 PetscCall(MatDestroy(&B2)); 447 448 if (!issymmetric) { 449 PetscCall(MatCreateSubMatrix(A,is,is2,MAT_INITIAL_MATRIX,&A2)); 450 PetscCall(MatCreateSubMatrix(B,is,is2,MAT_INITIAL_MATRIX,&B2)); 451 PetscCall(MatCreateSubMatrix(A,is,is2,MAT_REUSE_MATRIX,&A2)); 452 PetscCall(MatCreateSubMatrix(B,is,is2,MAT_REUSE_MATRIX,&B2)); 453 PetscCall(CheckMat(A2,B2,PETSC_FALSE,"second MatCreateSubMatrix")); 454 } 455 456 PetscCall(MatDestroy(&A2)); 457 PetscCall(MatDestroy(&B2)); 458 PetscCall(ISDestroy(&is)); 459 PetscCall(ISDestroy(&is2)); 460 461 /* test MatCreateSubMatrices */ 462 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Test MatCreateSubMatrices\n")); 463 PetscCall(MatGetLayouts(A,&rlayout,&clayout)); 464 PetscCall(PetscLayoutGetRanges(rlayout,&rrange)); 465 PetscCall(PetscLayoutGetRanges(clayout,&crange)); 466 lrank = (size + rank - 1)%size; 467 rrank = (rank + 1)%size; 468 PetscCall(ISCreateStride(PETSC_COMM_SELF,(rrange[lrank+1] - rrange[lrank]),rrange[lrank],1,&irow[0])); 469 PetscCall(ISCreateStride(PETSC_COMM_SELF,(crange[rrank+1] - crange[rrank]),crange[rrank],1,&icol[0])); 470 PetscCall(ISCreateStride(PETSC_COMM_SELF,(rrange[rrank+1] - rrange[rrank]),rrange[rrank],1,&irow[1])); 471 PetscCall(ISCreateStride(PETSC_COMM_SELF,(crange[lrank+1] - crange[lrank]),crange[lrank],1,&icol[1])); 472 PetscCall(MatCreateSubMatrices(A,2,irow,icol,MAT_INITIAL_MATRIX,&Asub)); 473 PetscCall(MatCreateSubMatrices(B,2,irow,icol,MAT_INITIAL_MATRIX,&Bsub)); 474 PetscCall(CheckMat(Asub[0],Bsub[0],PETSC_FALSE,"MatCreateSubMatrices[0]")); 475 PetscCall(CheckMat(Asub[1],Bsub[1],PETSC_FALSE,"MatCreateSubMatrices[1]")); 476 PetscCall(MatCreateSubMatrices(A,2,irow,icol,MAT_REUSE_MATRIX,&Asub)); 477 PetscCall(MatCreateSubMatrices(B,2,irow,icol,MAT_REUSE_MATRIX,&Bsub)); 478 PetscCall(CheckMat(Asub[0],Bsub[0],PETSC_FALSE,"MatCreateSubMatrices[0]")); 479 PetscCall(CheckMat(Asub[1],Bsub[1],PETSC_FALSE,"MatCreateSubMatrices[1]")); 480 PetscCall(MatDestroySubMatrices(2,&Asub)); 481 PetscCall(MatDestroySubMatrices(2,&Bsub)); 482 PetscCall(ISDestroy(&irow[0])); 483 PetscCall(ISDestroy(&irow[1])); 484 PetscCall(ISDestroy(&icol[0])); 485 PetscCall(ISDestroy(&icol[1])); 486 487 /* Create an IS required by MatZeroRows(): just rank zero provides the rows to be eliminated */ 488 if (size > 1) { 489 if (rank == 0) { 490 PetscInt st,len; 491 492 st = (m+1)/2; 493 len = PetscMin(m/2,PetscMax(m-(m+1)/2-1,0)); 494 PetscCall(ISCreateStride(PETSC_COMM_WORLD,len,st,1,&is)); 495 } else { 496 PetscCall(ISCreateStride(PETSC_COMM_WORLD,0,0,1,&is)); 497 } 498 } else { 499 PetscCall(ISCreateStride(PETSC_COMM_WORLD,1,0,1,&is)); 500 } 501 502 if (squaretest) { /* tests for square matrices only, with same maps for rows and columns */ 503 PetscInt *idx0, *idx1, n0, n1; 504 IS Ais[2], Bis[2]; 505 506 /* test MatDiagonalSet */ 507 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Test MatDiagonalSet\n")); 508 PetscCall(MatDuplicate(A,MAT_COPY_VALUES,&A2)); 509 PetscCall(MatDuplicate(B,MAT_COPY_VALUES,&B2)); 510 PetscCall(MatCreateVecs(A,NULL,&x)); 511 PetscCall(VecSetRandom(x,NULL)); 512 PetscCall(MatDiagonalSet(A2,x,INSERT_VALUES)); 513 PetscCall(MatDiagonalSet(B2,x,INSERT_VALUES)); 514 PetscCall(CheckMat(A2,B2,PETSC_FALSE,"MatDiagonalSet")); 515 PetscCall(VecDestroy(&x)); 516 PetscCall(MatDestroy(&A2)); 517 PetscCall(MatDestroy(&B2)); 518 519 /* test MatShift (MatShift_IS internally uses MatDiagonalSet_IS with ADD_VALUES) */ 520 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Test MatShift\n")); 521 PetscCall(MatDuplicate(A,MAT_COPY_VALUES,&A2)); 522 PetscCall(MatDuplicate(B,MAT_COPY_VALUES,&B2)); 523 PetscCall(MatShift(A2,2.0)); 524 PetscCall(MatShift(B2,2.0)); 525 PetscCall(CheckMat(A2,B2,PETSC_FALSE,"MatShift")); 526 PetscCall(MatDestroy(&A2)); 527 PetscCall(MatDestroy(&B2)); 528 529 /* nonzero diag value is supported for square matrices only */ 530 PetscCall(TestMatZeroRows(A,B,PETSC_TRUE,is,diag)); 531 532 /* test MatIncreaseOverlap */ 533 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Test MatIncreaseOverlap\n")); 534 PetscCall(MatGetOwnershipRange(A,&rst,&ren)); 535 n0 = (ren - rst)/2; 536 n1 = (ren - rst)/3; 537 PetscCall(PetscMalloc1(n0,&idx0)); 538 PetscCall(PetscMalloc1(n1,&idx1)); 539 for (i = 0; i < n0; i++) idx0[i] = ren - i - 1; 540 for (i = 0; i < n1; i++) idx1[i] = rst + i; 541 PetscCall(ISCreateGeneral(PETSC_COMM_WORLD,n0,idx0,PETSC_OWN_POINTER,&Ais[0])); 542 PetscCall(ISCreateGeneral(PETSC_COMM_WORLD,n1,idx1,PETSC_OWN_POINTER,&Ais[1])); 543 PetscCall(ISCreateGeneral(PETSC_COMM_WORLD,n0,idx0,PETSC_COPY_VALUES,&Bis[0])); 544 PetscCall(ISCreateGeneral(PETSC_COMM_WORLD,n1,idx1,PETSC_COPY_VALUES,&Bis[1])); 545 PetscCall(MatIncreaseOverlap(A,2,Ais,3)); 546 PetscCall(MatIncreaseOverlap(B,2,Bis,3)); 547 /* Non deterministic output! */ 548 PetscCall(ISSort(Ais[0])); 549 PetscCall(ISSort(Ais[1])); 550 PetscCall(ISSort(Bis[0])); 551 PetscCall(ISSort(Bis[1])); 552 PetscCall(ISView(Ais[0],NULL)); 553 PetscCall(ISView(Bis[0],NULL)); 554 PetscCall(ISView(Ais[1],NULL)); 555 PetscCall(ISView(Bis[1],NULL)); 556 PetscCall(MatCreateSubMatrices(A,2,Ais,Ais,MAT_INITIAL_MATRIX,&Asub)); 557 PetscCall(MatCreateSubMatrices(B,2,Bis,Ais,MAT_INITIAL_MATRIX,&Bsub)); 558 PetscCall(CheckMat(Asub[0],Bsub[0],PETSC_FALSE,"MatIncreaseOverlap[0]")); 559 PetscCall(CheckMat(Asub[1],Bsub[1],PETSC_FALSE,"MatIncreaseOverlap[1]")); 560 PetscCall(MatDestroySubMatrices(2,&Asub)); 561 PetscCall(MatDestroySubMatrices(2,&Bsub)); 562 PetscCall(ISDestroy(&Ais[0])); 563 PetscCall(ISDestroy(&Ais[1])); 564 PetscCall(ISDestroy(&Bis[0])); 565 PetscCall(ISDestroy(&Bis[1])); 566 567 } 568 PetscCall(TestMatZeroRows(A,B,squaretest,is,0.0)); 569 PetscCall(ISDestroy(&is)); 570 571 /* test MatTranspose */ 572 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Test MatTranspose\n")); 573 PetscCall(MatTranspose(A,MAT_INITIAL_MATRIX,&A2)); 574 PetscCall(MatTranspose(B,MAT_INITIAL_MATRIX,&B2)); 575 PetscCall(CheckMat(A2,B2,PETSC_FALSE,"initial matrix MatTranspose")); 576 577 PetscCall(MatTranspose(A,MAT_REUSE_MATRIX,&A2)); 578 PetscCall(CheckMat(A2,B2,PETSC_FALSE,"reuse matrix (not in place) MatTranspose")); 579 PetscCall(MatDestroy(&A2)); 580 581 PetscCall(MatDuplicate(A,MAT_COPY_VALUES,&A2)); 582 PetscCall(MatTranspose(A2,MAT_INPLACE_MATRIX,&A2)); 583 PetscCall(CheckMat(A2,B2,PETSC_FALSE,"reuse matrix (in place) MatTranspose")); 584 PetscCall(MatDestroy(&A2)); 585 586 PetscCall(MatTranspose(A,MAT_INITIAL_MATRIX,&A2)); 587 PetscCall(CheckMat(A2,B2,PETSC_FALSE,"reuse matrix (different type) MatTranspose")); 588 PetscCall(MatDestroy(&A2)); 589 PetscCall(MatDestroy(&B2)); 590 591 /* test MatISFixLocalEmpty */ 592 if (isaij) { 593 PetscInt r[2]; 594 595 r[0] = 0; 596 r[1] = PetscMin(m,n)-1; 597 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Test MatISFixLocalEmpty\n")); 598 PetscCall(MatDuplicate(A,MAT_COPY_VALUES,&A2)); 599 600 PetscCall(MatISFixLocalEmpty(A2,PETSC_TRUE)); 601 PetscCall(MatAssemblyBegin(A2,MAT_FINAL_ASSEMBLY)); 602 PetscCall(MatAssemblyEnd(A2,MAT_FINAL_ASSEMBLY)); 603 PetscCall(CheckMat(A2,B,PETSC_FALSE,"MatISFixLocalEmpty (null)")); 604 605 PetscCall(MatZeroRows(A2,2,r,0.0,NULL,NULL)); 606 PetscCall(MatViewFromOptions(A2,NULL,"-fixempty_view")); 607 PetscCall(MatDuplicate(B,MAT_COPY_VALUES,&B2)); 608 PetscCall(MatZeroRows(B2,2,r,0.0,NULL,NULL)); 609 PetscCall(MatISFixLocalEmpty(A2,PETSC_TRUE)); 610 PetscCall(MatAssemblyBegin(A2,MAT_FINAL_ASSEMBLY)); 611 PetscCall(MatAssemblyEnd(A2,MAT_FINAL_ASSEMBLY)); 612 PetscCall(MatViewFromOptions(A2,NULL,"-fixempty_view")); 613 PetscCall(CheckMat(A2,B2,PETSC_FALSE,"MatISFixLocalEmpty (rows)")); 614 PetscCall(MatDestroy(&A2)); 615 616 PetscCall(MatDuplicate(A,MAT_COPY_VALUES,&A2)); 617 PetscCall(MatZeroRows(A2,2,r,0.0,NULL,NULL)); 618 PetscCall(MatTranspose(A2,MAT_INPLACE_MATRIX,&A2)); 619 PetscCall(MatTranspose(B2,MAT_INPLACE_MATRIX,&B2)); 620 PetscCall(MatViewFromOptions(A2,NULL,"-fixempty_view")); 621 PetscCall(MatISFixLocalEmpty(A2,PETSC_TRUE)); 622 PetscCall(MatAssemblyBegin(A2,MAT_FINAL_ASSEMBLY)); 623 PetscCall(MatAssemblyEnd(A2,MAT_FINAL_ASSEMBLY)); 624 PetscCall(MatViewFromOptions(A2,NULL,"-fixempty_view")); 625 PetscCall(CheckMat(A2,B2,PETSC_FALSE,"MatISFixLocalEmpty (cols)")); 626 627 PetscCall(MatDestroy(&A2)); 628 PetscCall(MatDestroy(&B2)); 629 630 if (squaretest) { 631 PetscCall(MatDuplicate(A,MAT_COPY_VALUES,&A2)); 632 PetscCall(MatDuplicate(B,MAT_COPY_VALUES,&B2)); 633 PetscCall(MatZeroRowsColumns(A2,2,r,0.0,NULL,NULL)); 634 PetscCall(MatZeroRowsColumns(B2,2,r,0.0,NULL,NULL)); 635 PetscCall(MatViewFromOptions(A2,NULL,"-fixempty_view")); 636 PetscCall(MatISFixLocalEmpty(A2,PETSC_TRUE)); 637 PetscCall(MatAssemblyBegin(A2,MAT_FINAL_ASSEMBLY)); 638 PetscCall(MatAssemblyEnd(A2,MAT_FINAL_ASSEMBLY)); 639 PetscCall(MatViewFromOptions(A2,NULL,"-fixempty_view")); 640 PetscCall(CheckMat(A2,B2,PETSC_FALSE,"MatISFixLocalEmpty (rows+cols)")); 641 PetscCall(MatDestroy(&A2)); 642 PetscCall(MatDestroy(&B2)); 643 } 644 645 } 646 647 /* test MatInvertBlockDiagonal 648 special cases for block-diagonal matrices */ 649 if (m == n) { 650 ISLocalToGlobalMapping map; 651 Mat Abd,Bbd; 652 IS is,bis; 653 const PetscScalar *isbd,*aijbd; 654 PetscScalar *vals; 655 const PetscInt *sts,*idxs; 656 PetscInt *idxs2,diff,perm,nl,bs,st,en,in; 657 PetscBool ok; 658 659 for (diff = 0; diff < 3; diff++) { 660 for (perm = 0; perm < 3; perm++) { 661 for (bs = 1; bs < 4; bs++) { 662 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Test MatInvertBlockDiagonal blockdiag %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n",n,diff,perm,bs)); 663 PetscCall(PetscMalloc1(bs*bs,&vals)); 664 PetscCall(MatGetOwnershipRanges(A,&sts)); 665 switch (diff) { 666 case 1: /* inverted layout by processes */ 667 in = 1; 668 st = sts[size - rank - 1]; 669 en = sts[size - rank]; 670 nl = en - st; 671 break; 672 case 2: /* round-robin layout */ 673 in = size; 674 st = rank; 675 nl = n/size; 676 if (rank < n%size) nl++; 677 break; 678 default: /* same layout */ 679 in = 1; 680 st = sts[rank]; 681 en = sts[rank + 1]; 682 nl = en - st; 683 break; 684 } 685 PetscCall(ISCreateStride(PETSC_COMM_WORLD,nl,st,in,&is)); 686 PetscCall(ISGetLocalSize(is,&nl)); 687 PetscCall(ISGetIndices(is,&idxs)); 688 PetscCall(PetscMalloc1(nl,&idxs2)); 689 for (i=0;i<nl;i++) { 690 switch (perm) { /* invert some of the indices */ 691 case 2: 692 idxs2[i] = rank%2 ? idxs[i] : idxs[nl-i-1]; 693 break; 694 case 1: 695 idxs2[i] = rank%2 ? idxs[nl-i-1] : idxs[i]; 696 break; 697 default: 698 idxs2[i] = idxs[i]; 699 break; 700 } 701 } 702 PetscCall(ISRestoreIndices(is,&idxs)); 703 PetscCall(ISCreateBlock(PETSC_COMM_WORLD,bs,nl,idxs2,PETSC_OWN_POINTER,&bis)); 704 PetscCall(ISLocalToGlobalMappingCreateIS(bis,&map)); 705 PetscCall(MatCreateIS(PETSC_COMM_WORLD,bs,PETSC_DECIDE,PETSC_DECIDE,bs*n,bs*n,map,map,&Abd)); 706 PetscCall(ISLocalToGlobalMappingDestroy(&map)); 707 PetscCall(MatISSetPreallocation(Abd,bs,NULL,0,NULL)); 708 for (i=0;i<nl;i++) { 709 PetscInt b1,b2; 710 711 for (b1=0;b1<bs;b1++) { 712 for (b2=0;b2<bs;b2++) { 713 vals[b1*bs + b2] = i*bs*bs + b1*bs + b2 + 1 + (b1 == b2 ? 1.0 : 0); 714 } 715 } 716 PetscCall(MatSetValuesBlockedLocal(Abd,1,&i,1,&i,vals,INSERT_VALUES)); 717 } 718 PetscCall(MatAssemblyBegin(Abd,MAT_FINAL_ASSEMBLY)); 719 PetscCall(MatAssemblyEnd(Abd,MAT_FINAL_ASSEMBLY)); 720 PetscCall(MatConvert(Abd,MATAIJ,MAT_INITIAL_MATRIX,&Bbd)); 721 PetscCall(MatInvertBlockDiagonal(Abd,&isbd)); 722 PetscCall(MatInvertBlockDiagonal(Bbd,&aijbd)); 723 PetscCall(MatGetLocalSize(Bbd,&nl,NULL)); 724 ok = PETSC_TRUE; 725 for (i=0;i<nl/bs;i++) { 726 PetscInt b1,b2; 727 728 for (b1=0;b1<bs;b1++) { 729 for (b2=0;b2<bs;b2++) { 730 if (PetscAbsScalar(isbd[i*bs*bs+b1*bs + b2]-aijbd[i*bs*bs+b1*bs + b2]) > PETSC_SMALL) ok = PETSC_FALSE; 731 if (!ok) { 732 PetscCall(PetscPrintf(PETSC_COMM_SELF,"[%d] ERROR block %" PetscInt_FMT ", entry %" PetscInt_FMT " %" PetscInt_FMT ": %g %g\n",rank,i,b1,b2,(double)PetscAbsScalar(isbd[i*bs*bs+b1*bs + b2]),(double)PetscAbsScalar(aijbd[i*bs*bs+b1*bs + b2]))); 733 break; 734 } 735 } 736 if (!ok) break; 737 } 738 if (!ok) break; 739 } 740 PetscCall(MatDestroy(&Abd)); 741 PetscCall(MatDestroy(&Bbd)); 742 PetscCall(PetscFree(vals)); 743 PetscCall(ISDestroy(&is)); 744 PetscCall(ISDestroy(&bis)); 745 } 746 } 747 } 748 } 749 750 /* test MatGetDiagonalBlock */ 751 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Test MatGetDiagonalBlock\n")); 752 PetscCall(MatGetDiagonalBlock(A,&A2)); 753 PetscCall(MatGetDiagonalBlock(B,&B2)); 754 PetscCall(CheckMat(A2,B2,PETSC_FALSE,"MatGetDiagonalBlock")); 755 PetscCall(MatScale(A,2.0)); 756 PetscCall(MatScale(B,2.0)); 757 PetscCall(MatGetDiagonalBlock(A,&A2)); 758 PetscCall(MatGetDiagonalBlock(B,&B2)); 759 PetscCall(CheckMat(A2,B2,PETSC_FALSE,"MatGetDiagonalBlock")); 760 761 /* free testing matrices */ 762 PetscCall(ISLocalToGlobalMappingDestroy(&cmap)); 763 PetscCall(ISLocalToGlobalMappingDestroy(&rmap)); 764 PetscCall(MatDestroy(&A)); 765 PetscCall(MatDestroy(&B)); 766 PetscCall(PetscFinalize()); 767 return 0; 768 } 769 770 PetscErrorCode CheckMat(Mat A, Mat B, PetscBool usemult, const char* func) 771 { 772 Mat Bcheck; 773 PetscReal error; 774 775 PetscFunctionBeginUser; 776 if (!usemult && B) { 777 PetscBool hasnorm; 778 779 PetscCall(MatHasOperation(B,MATOP_NORM,&hasnorm)); 780 if (!hasnorm) usemult = PETSC_TRUE; 781 } 782 if (!usemult) { 783 if (B) { 784 MatType Btype; 785 786 PetscCall(MatGetType(B,&Btype)); 787 PetscCall(MatConvert(A,Btype,MAT_INITIAL_MATRIX,&Bcheck)); 788 } else { 789 PetscCall(MatConvert(A,MATAIJ,MAT_INITIAL_MATRIX,&Bcheck)); 790 } 791 if (B) { /* if B is present, subtract it */ 792 PetscCall(MatAXPY(Bcheck,-1.,B,DIFFERENT_NONZERO_PATTERN)); 793 } 794 PetscCall(MatNorm(Bcheck,NORM_INFINITY,&error)); 795 if (error > PETSC_SQRT_MACHINE_EPSILON) { 796 ISLocalToGlobalMapping rl2g,cl2g; 797 798 PetscCall(PetscObjectSetName((PetscObject)Bcheck,"Bcheck")); 799 PetscCall(MatView(Bcheck,NULL)); 800 if (B) { 801 PetscCall(PetscObjectSetName((PetscObject)B,"B")); 802 PetscCall(MatView(B,NULL)); 803 PetscCall(MatDestroy(&Bcheck)); 804 PetscCall(MatConvert(A,MATAIJ,MAT_INITIAL_MATRIX,&Bcheck)); 805 PetscCall(PetscObjectSetName((PetscObject)Bcheck,"Assembled A")); 806 PetscCall(MatView(Bcheck,NULL)); 807 } 808 PetscCall(MatDestroy(&Bcheck)); 809 PetscCall(PetscObjectSetName((PetscObject)A,"A")); 810 PetscCall(MatView(A,NULL)); 811 PetscCall(MatGetLocalToGlobalMapping(A,&rl2g,&cl2g)); 812 if (rl2g) PetscCall(ISLocalToGlobalMappingView(rl2g,NULL)); 813 if (cl2g) PetscCall(ISLocalToGlobalMappingView(cl2g,NULL)); 814 SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"ERROR ON %s: %g",func,(double)error); 815 } 816 PetscCall(MatDestroy(&Bcheck)); 817 } else { 818 PetscBool ok,okt; 819 820 PetscCall(MatMultEqual(A,B,3,&ok)); 821 PetscCall(MatMultTransposeEqual(A,B,3,&okt)); 822 PetscCheck(ok && okt,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"ERROR ON %s: mult ok ? %d, multtranspose ok ? %d",func,ok,okt); 823 } 824 PetscFunctionReturn(0); 825 } 826 827 PetscErrorCode TestMatZeroRows(Mat A, Mat Afull, PetscBool squaretest, IS is, PetscScalar diag) 828 { 829 Mat B,Bcheck,B2 = NULL,lB; 830 Vec x = NULL, b = NULL, b2 = NULL; 831 ISLocalToGlobalMapping l2gr,l2gc; 832 PetscReal error; 833 char diagstr[16]; 834 const PetscInt *idxs; 835 PetscInt rst,ren,i,n,N,d; 836 PetscMPIInt rank; 837 PetscBool miss,haszerorows; 838 839 PetscFunctionBeginUser; 840 if (diag == 0.) { 841 PetscCall(PetscStrcpy(diagstr,"zero")); 842 } else { 843 PetscCall(PetscStrcpy(diagstr,"nonzero")); 844 } 845 PetscCall(ISView(is,NULL)); 846 PetscCall(MatGetLocalToGlobalMapping(A,&l2gr,&l2gc)); 847 /* tests MatDuplicate and MatCopy */ 848 if (diag == 0.) { 849 PetscCall(MatDuplicate(A,MAT_COPY_VALUES,&B)); 850 } else { 851 PetscCall(MatDuplicate(A,MAT_DO_NOT_COPY_VALUES,&B)); 852 PetscCall(MatCopy(A,B,SAME_NONZERO_PATTERN)); 853 } 854 PetscCall(MatISGetLocalMat(B,&lB)); 855 PetscCall(MatHasOperation(lB,MATOP_ZERO_ROWS,&haszerorows)); 856 if (squaretest && haszerorows) { 857 858 PetscCall(MatCreateVecs(B,&x,&b)); 859 PetscCall(MatDuplicate(B,MAT_COPY_VALUES,&B2)); 860 PetscCall(VecSetLocalToGlobalMapping(b,l2gr)); 861 PetscCall(VecSetLocalToGlobalMapping(x,l2gc)); 862 PetscCall(VecSetRandom(x,NULL)); 863 PetscCall(VecSetRandom(b,NULL)); 864 /* mimic b[is] = x[is] */ 865 PetscCall(VecDuplicate(b,&b2)); 866 PetscCall(VecSetLocalToGlobalMapping(b2,l2gr)); 867 PetscCall(VecCopy(b,b2)); 868 PetscCall(ISGetLocalSize(is,&n)); 869 PetscCall(ISGetIndices(is,&idxs)); 870 PetscCall(VecGetSize(x,&N)); 871 for (i=0;i<n;i++) { 872 if (0 <= idxs[i] && idxs[i] < N) { 873 PetscCall(VecSetValue(b2,idxs[i],diag,INSERT_VALUES)); 874 PetscCall(VecSetValue(x,idxs[i],1.,INSERT_VALUES)); 875 } 876 } 877 PetscCall(VecAssemblyBegin(b2)); 878 PetscCall(VecAssemblyEnd(b2)); 879 PetscCall(VecAssemblyBegin(x)); 880 PetscCall(VecAssemblyEnd(x)); 881 PetscCall(ISRestoreIndices(is,&idxs)); 882 /* test ZeroRows on MATIS */ 883 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Test MatZeroRows (diag %s)\n",diagstr)); 884 PetscCall(MatZeroRowsIS(B,is,diag,x,b)); 885 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Test MatZeroRowsColumns (diag %s)\n",diagstr)); 886 PetscCall(MatZeroRowsColumnsIS(B2,is,diag,NULL,NULL)); 887 } else if (haszerorows) { 888 /* test ZeroRows on MATIS */ 889 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Test MatZeroRows (diag %s)\n",diagstr)); 890 PetscCall(MatZeroRowsIS(B,is,diag,NULL,NULL)); 891 b = b2 = x = NULL; 892 } else { 893 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Skipping MatZeroRows (diag %s)\n",diagstr)); 894 b = b2 = x = NULL; 895 } 896 897 if (squaretest && haszerorows) { 898 PetscCall(VecAXPY(b2,-1.,b)); 899 PetscCall(VecNorm(b2,NORM_INFINITY,&error)); 900 PetscCheck(error <= PETSC_SQRT_MACHINE_EPSILON,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"ERROR IN ZEROROWS ON B %g (diag %s)",(double)error,diagstr); 901 } 902 903 /* test MatMissingDiagonal */ 904 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Test MatMissingDiagonal\n")); 905 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank)); 906 PetscCall(MatMissingDiagonal(B,&miss,&d)); 907 PetscCall(MatGetOwnershipRange(B,&rst,&ren)); 908 PetscCall(PetscViewerASCIIPushSynchronized(PETSC_VIEWER_STDOUT_WORLD)); 909 PetscCall(PetscViewerASCIISynchronizedPrintf(PETSC_VIEWER_STDOUT_WORLD, "[%d] [%" PetscInt_FMT ",%" PetscInt_FMT ") Missing %d, row %" PetscInt_FMT " (diag %s)\n",rank,rst,ren,(int)miss,d,diagstr)); 910 PetscCall(PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD)); 911 PetscCall(PetscViewerASCIIPopSynchronized(PETSC_VIEWER_STDOUT_WORLD)); 912 913 PetscCall(VecDestroy(&x)); 914 PetscCall(VecDestroy(&b)); 915 PetscCall(VecDestroy(&b2)); 916 917 /* check the result of ZeroRows with that from MPIAIJ routines 918 assuming that MatConvert_IS_XAIJ and MatZeroRows_MPIAIJ work fine */ 919 if (haszerorows) { 920 PetscCall(MatDuplicate(Afull,MAT_COPY_VALUES,&Bcheck)); 921 PetscCall(MatSetOption(Bcheck,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE)); 922 PetscCall(MatZeroRowsIS(Bcheck,is,diag,NULL,NULL)); 923 PetscCall(CheckMat(B,Bcheck,PETSC_FALSE,"Zerorows")); 924 PetscCall(MatDestroy(&Bcheck)); 925 } 926 PetscCall(MatDestroy(&B)); 927 928 if (B2) { /* test MatZeroRowsColumns */ 929 PetscCall(MatDuplicate(Afull,MAT_COPY_VALUES,&B)); 930 PetscCall(MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE)); 931 PetscCall(MatZeroRowsColumnsIS(B,is,diag,NULL,NULL)); 932 PetscCall(CheckMat(B2,B,PETSC_FALSE,"MatZeroRowsColumns")); 933 PetscCall(MatDestroy(&B)); 934 PetscCall(MatDestroy(&B2)); 935 } 936 PetscFunctionReturn(0); 937 } 938 939 /*TEST 940 941 test: 942 args: -test_trans 943 944 test: 945 suffix: 2 946 nsize: 4 947 args: -matis_convert_local_nest -nr 3 -nc 4 948 949 test: 950 suffix: 3 951 nsize: 5 952 args: -m 11 -n 10 -matis_convert_local_nest -nr 2 -nc 1 953 954 test: 955 suffix: 4 956 nsize: 6 957 args: -m 9 -n 12 -test_trans -nr 2 -nc 7 958 959 test: 960 suffix: 5 961 nsize: 6 962 args: -m 12 -n 12 -test_trans -nr 3 -nc 1 963 964 test: 965 suffix: 6 966 args: -m 12 -n 12 -test_trans -nr 2 -nc 3 -diffmap 967 968 test: 969 suffix: 7 970 args: -m 12 -n 12 -test_trans -nr 2 -nc 3 -diffmap -permmap 971 972 test: 973 suffix: 8 974 args: -m 12 -n 17 -test_trans -nr 2 -nc 3 -permmap 975 976 test: 977 suffix: 9 978 nsize: 5 979 args: -m 12 -n 12 -test_trans -nr 2 -nc 3 -diffmap 980 981 test: 982 suffix: 10 983 nsize: 5 984 args: -m 12 -n 12 -test_trans -nr 2 -nc 3 -diffmap -permmap 985 986 test: 987 suffix: vscat_default 988 nsize: 5 989 args: -m 12 -n 17 -test_trans -nr 2 -nc 3 -permmap 990 output_file: output/ex23_11.out 991 992 test: 993 suffix: 12 994 nsize: 3 995 args: -m 12 -n 12 -symmetric -matis_localmat_type sbaij -test_trans -nr 2 -nc 3 996 997 testset: 998 output_file: output/ex23_13.out 999 nsize: 3 1000 args: -m 12 -n 17 -test_trans -nr 2 -nc 3 -diffmap -permmap 1001 filter: grep -v "type:" 1002 test: 1003 suffix: baij 1004 args: -matis_localmat_type baij 1005 test: 1006 requires: viennacl 1007 suffix: viennacl 1008 args: -matis_localmat_type aijviennacl 1009 test: 1010 requires: cuda 1011 suffix: cusparse 1012 args: -matis_localmat_type aijcusparse 1013 1014 test: 1015 suffix: negrep 1016 nsize: {{1 3}separate output} 1017 args: -m {{5 7}separate output} -n {{5 7}separate output} -test_trans -nr 2 -nc 3 -negmap {{0 1}separate output} -repmap {{0 1}separate output} -permmap -diffmap {{0 1}separate output} 1018 1019 TEST*/ 1020