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