1 2 static char help[] = "Tests the use of interface functions for MATIS matrices.\n\ 3 This example tests: MatZeroRows(), MatZeroRowsLocal(), MatView(), MatDuplicate(),\n\ 4 MatCopy(), MatCreateSubMatrix(), MatGetLocalSubMatrix(), MatAXPY(), MatShift()\n\ 5 MatDiagonalSet(), MatTranspose() and MatPtAP(). It also tests some\n\ 6 conversion routines.\n"; 7 8 #include <petscmat.h> 9 10 PetscErrorCode TestMatZeroRows(Mat,Mat,PetscBool,IS,PetscScalar); 11 PetscErrorCode CheckMat(Mat,Mat,PetscBool,const char*); 12 13 int main(int argc,char **args) 14 { 15 Mat A,B,A2,B2,T; 16 Mat Aee,Aeo,Aoe,Aoo; 17 Mat *mats; 18 Vec x,y; 19 MatInfo info; 20 ISLocalToGlobalMapping cmap,rmap; 21 IS is,is2,reven,rodd,ceven,codd; 22 IS *rows,*cols; 23 MatType lmtype; 24 PetscScalar diag = 2.; 25 PetscInt n,m,i,lm,ln; 26 PetscInt rst,ren,cst,cen,nr,nc; 27 PetscMPIInt rank,size; 28 PetscBool testT,squaretest,isaij; 29 PetscBool permute = PETSC_FALSE, negmap = PETSC_FALSE, repmap = PETSC_FALSE; 30 PetscBool diffmap = PETSC_TRUE, symmetric = PETSC_FALSE, issymmetric; 31 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 /* Create an IS required by MatZeroRows(): just rank zero provides the rows to be eliminated */ 463 if (size > 1) { 464 if (rank == 0) { 465 PetscInt st,len; 466 467 st = (m+1)/2; 468 len = PetscMin(m/2,PetscMax(m-(m+1)/2-1,0)); 469 PetscCall(ISCreateStride(PETSC_COMM_WORLD,len,st,1,&is)); 470 } else { 471 PetscCall(ISCreateStride(PETSC_COMM_WORLD,0,0,1,&is)); 472 } 473 } else { 474 PetscCall(ISCreateStride(PETSC_COMM_WORLD,1,0,1,&is)); 475 } 476 477 if (squaretest) { /* tests for square matrices only, with same maps for rows and columns */ 478 /* test MatDiagonalSet */ 479 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Test MatDiagonalSet\n")); 480 PetscCall(MatDuplicate(A,MAT_COPY_VALUES,&A2)); 481 PetscCall(MatDuplicate(B,MAT_COPY_VALUES,&B2)); 482 PetscCall(MatCreateVecs(A,NULL,&x)); 483 PetscCall(VecSetRandom(x,NULL)); 484 PetscCall(MatDiagonalSet(A2,x,INSERT_VALUES)); 485 PetscCall(MatDiagonalSet(B2,x,INSERT_VALUES)); 486 PetscCall(CheckMat(A2,B2,PETSC_FALSE,"MatDiagonalSet")); 487 PetscCall(VecDestroy(&x)); 488 PetscCall(MatDestroy(&A2)); 489 PetscCall(MatDestroy(&B2)); 490 491 /* test MatShift (MatShift_IS internally uses MatDiagonalSet_IS with ADD_VALUES) */ 492 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Test MatShift\n")); 493 PetscCall(MatDuplicate(A,MAT_COPY_VALUES,&A2)); 494 PetscCall(MatDuplicate(B,MAT_COPY_VALUES,&B2)); 495 PetscCall(MatShift(A2,2.0)); 496 PetscCall(MatShift(B2,2.0)); 497 PetscCall(CheckMat(A2,B2,PETSC_FALSE,"MatShift")); 498 PetscCall(MatDestroy(&A2)); 499 PetscCall(MatDestroy(&B2)); 500 501 /* nonzero diag value is supported for square matrices only */ 502 PetscCall(TestMatZeroRows(A,B,PETSC_TRUE,is,diag)); 503 } 504 PetscCall(TestMatZeroRows(A,B,squaretest,is,0.0)); 505 PetscCall(ISDestroy(&is)); 506 507 /* test MatTranspose */ 508 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Test MatTranspose\n")); 509 PetscCall(MatTranspose(A,MAT_INITIAL_MATRIX,&A2)); 510 PetscCall(MatTranspose(B,MAT_INITIAL_MATRIX,&B2)); 511 PetscCall(CheckMat(A2,B2,PETSC_FALSE,"initial matrix MatTranspose")); 512 513 PetscCall(MatTranspose(A,MAT_REUSE_MATRIX,&A2)); 514 PetscCall(CheckMat(A2,B2,PETSC_FALSE,"reuse matrix (not in place) MatTranspose")); 515 PetscCall(MatDestroy(&A2)); 516 517 PetscCall(MatDuplicate(A,MAT_COPY_VALUES,&A2)); 518 PetscCall(MatTranspose(A2,MAT_INPLACE_MATRIX,&A2)); 519 PetscCall(CheckMat(A2,B2,PETSC_FALSE,"reuse matrix (in place) MatTranspose")); 520 PetscCall(MatDestroy(&A2)); 521 522 PetscCall(MatTranspose(A,MAT_INITIAL_MATRIX,&A2)); 523 PetscCall(CheckMat(A2,B2,PETSC_FALSE,"reuse matrix (different type) MatTranspose")); 524 PetscCall(MatDestroy(&A2)); 525 PetscCall(MatDestroy(&B2)); 526 527 /* test MatISFixLocalEmpty */ 528 if (isaij) { 529 PetscInt r[2]; 530 531 r[0] = 0; 532 r[1] = PetscMin(m,n)-1; 533 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Test MatISFixLocalEmpty\n")); 534 PetscCall(MatDuplicate(A,MAT_COPY_VALUES,&A2)); 535 536 PetscCall(MatISFixLocalEmpty(A2,PETSC_TRUE)); 537 PetscCall(MatAssemblyBegin(A2,MAT_FINAL_ASSEMBLY)); 538 PetscCall(MatAssemblyEnd(A2,MAT_FINAL_ASSEMBLY)); 539 PetscCall(CheckMat(A2,B,PETSC_FALSE,"MatISFixLocalEmpty (null)")); 540 541 PetscCall(MatZeroRows(A2,2,r,0.0,NULL,NULL)); 542 PetscCall(MatViewFromOptions(A2,NULL,"-fixempty_view")); 543 PetscCall(MatDuplicate(B,MAT_COPY_VALUES,&B2)); 544 PetscCall(MatZeroRows(B2,2,r,0.0,NULL,NULL)); 545 PetscCall(MatISFixLocalEmpty(A2,PETSC_TRUE)); 546 PetscCall(MatAssemblyBegin(A2,MAT_FINAL_ASSEMBLY)); 547 PetscCall(MatAssemblyEnd(A2,MAT_FINAL_ASSEMBLY)); 548 PetscCall(MatViewFromOptions(A2,NULL,"-fixempty_view")); 549 PetscCall(CheckMat(A2,B2,PETSC_FALSE,"MatISFixLocalEmpty (rows)")); 550 PetscCall(MatDestroy(&A2)); 551 552 PetscCall(MatDuplicate(A,MAT_COPY_VALUES,&A2)); 553 PetscCall(MatZeroRows(A2,2,r,0.0,NULL,NULL)); 554 PetscCall(MatTranspose(A2,MAT_INPLACE_MATRIX,&A2)); 555 PetscCall(MatTranspose(B2,MAT_INPLACE_MATRIX,&B2)); 556 PetscCall(MatViewFromOptions(A2,NULL,"-fixempty_view")); 557 PetscCall(MatISFixLocalEmpty(A2,PETSC_TRUE)); 558 PetscCall(MatAssemblyBegin(A2,MAT_FINAL_ASSEMBLY)); 559 PetscCall(MatAssemblyEnd(A2,MAT_FINAL_ASSEMBLY)); 560 PetscCall(MatViewFromOptions(A2,NULL,"-fixempty_view")); 561 PetscCall(CheckMat(A2,B2,PETSC_FALSE,"MatISFixLocalEmpty (cols)")); 562 563 PetscCall(MatDestroy(&A2)); 564 PetscCall(MatDestroy(&B2)); 565 566 if (squaretest) { 567 PetscCall(MatDuplicate(A,MAT_COPY_VALUES,&A2)); 568 PetscCall(MatDuplicate(B,MAT_COPY_VALUES,&B2)); 569 PetscCall(MatZeroRowsColumns(A2,2,r,0.0,NULL,NULL)); 570 PetscCall(MatZeroRowsColumns(B2,2,r,0.0,NULL,NULL)); 571 PetscCall(MatViewFromOptions(A2,NULL,"-fixempty_view")); 572 PetscCall(MatISFixLocalEmpty(A2,PETSC_TRUE)); 573 PetscCall(MatAssemblyBegin(A2,MAT_FINAL_ASSEMBLY)); 574 PetscCall(MatAssemblyEnd(A2,MAT_FINAL_ASSEMBLY)); 575 PetscCall(MatViewFromOptions(A2,NULL,"-fixempty_view")); 576 PetscCall(CheckMat(A2,B2,PETSC_FALSE,"MatISFixLocalEmpty (rows+cols)")); 577 PetscCall(MatDestroy(&A2)); 578 PetscCall(MatDestroy(&B2)); 579 } 580 581 } 582 583 /* test MatInvertBlockDiagonal 584 special cases for block-diagonal matrices */ 585 if (m == n) { 586 ISLocalToGlobalMapping map; 587 Mat Abd,Bbd; 588 IS is,bis; 589 const PetscScalar *isbd,*aijbd; 590 PetscScalar *vals; 591 const PetscInt *sts,*idxs; 592 PetscInt *idxs2,diff,perm,nl,bs,st,en,in; 593 PetscBool ok; 594 595 for (diff = 0; diff < 3; diff++) { 596 for (perm = 0; perm < 3; perm++) { 597 for (bs = 1; bs < 4; bs++) { 598 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Test MatInvertBlockDiagonal blockdiag %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n",n,diff,perm,bs)); 599 PetscCall(PetscMalloc1(bs*bs,&vals)); 600 PetscCall(MatGetOwnershipRanges(A,&sts)); 601 switch (diff) { 602 case 1: /* inverted layout by processes */ 603 in = 1; 604 st = sts[size - rank - 1]; 605 en = sts[size - rank]; 606 nl = en - st; 607 break; 608 case 2: /* round-robin layout */ 609 in = size; 610 st = rank; 611 nl = n/size; 612 if (rank < n%size) nl++; 613 break; 614 default: /* same layout */ 615 in = 1; 616 st = sts[rank]; 617 en = sts[rank + 1]; 618 nl = en - st; 619 break; 620 } 621 PetscCall(ISCreateStride(PETSC_COMM_WORLD,nl,st,in,&is)); 622 PetscCall(ISGetLocalSize(is,&nl)); 623 PetscCall(ISGetIndices(is,&idxs)); 624 PetscCall(PetscMalloc1(nl,&idxs2)); 625 for (i=0;i<nl;i++) { 626 switch (perm) { /* invert some of the indices */ 627 case 2: 628 idxs2[i] = rank%2 ? idxs[i] : idxs[nl-i-1]; 629 break; 630 case 1: 631 idxs2[i] = rank%2 ? idxs[nl-i-1] : idxs[i]; 632 break; 633 default: 634 idxs2[i] = idxs[i]; 635 break; 636 } 637 } 638 PetscCall(ISRestoreIndices(is,&idxs)); 639 PetscCall(ISCreateBlock(PETSC_COMM_WORLD,bs,nl,idxs2,PETSC_OWN_POINTER,&bis)); 640 PetscCall(ISLocalToGlobalMappingCreateIS(bis,&map)); 641 PetscCall(MatCreateIS(PETSC_COMM_WORLD,bs,PETSC_DECIDE,PETSC_DECIDE,bs*n,bs*n,map,map,&Abd)); 642 PetscCall(ISLocalToGlobalMappingDestroy(&map)); 643 PetscCall(MatISSetPreallocation(Abd,bs,NULL,0,NULL)); 644 for (i=0;i<nl;i++) { 645 PetscInt b1,b2; 646 647 for (b1=0;b1<bs;b1++) { 648 for (b2=0;b2<bs;b2++) { 649 vals[b1*bs + b2] = i*bs*bs + b1*bs + b2 + 1 + (b1 == b2 ? 1.0 : 0); 650 } 651 } 652 PetscCall(MatSetValuesBlockedLocal(Abd,1,&i,1,&i,vals,INSERT_VALUES)); 653 } 654 PetscCall(MatAssemblyBegin(Abd,MAT_FINAL_ASSEMBLY)); 655 PetscCall(MatAssemblyEnd(Abd,MAT_FINAL_ASSEMBLY)); 656 PetscCall(MatConvert(Abd,MATAIJ,MAT_INITIAL_MATRIX,&Bbd)); 657 PetscCall(MatInvertBlockDiagonal(Abd,&isbd)); 658 PetscCall(MatInvertBlockDiagonal(Bbd,&aijbd)); 659 PetscCall(MatGetLocalSize(Bbd,&nl,NULL)); 660 ok = PETSC_TRUE; 661 for (i=0;i<nl/bs;i++) { 662 PetscInt b1,b2; 663 664 for (b1=0;b1<bs;b1++) { 665 for (b2=0;b2<bs;b2++) { 666 if (PetscAbsScalar(isbd[i*bs*bs+b1*bs + b2]-aijbd[i*bs*bs+b1*bs + b2]) > PETSC_SMALL) ok = PETSC_FALSE; 667 if (!ok) { 668 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]))); 669 break; 670 } 671 } 672 if (!ok) break; 673 } 674 if (!ok) break; 675 } 676 PetscCall(MatDestroy(&Abd)); 677 PetscCall(MatDestroy(&Bbd)); 678 PetscCall(PetscFree(vals)); 679 PetscCall(ISDestroy(&is)); 680 PetscCall(ISDestroy(&bis)); 681 } 682 } 683 } 684 } 685 /* free testing matrices */ 686 PetscCall(ISLocalToGlobalMappingDestroy(&cmap)); 687 PetscCall(ISLocalToGlobalMappingDestroy(&rmap)); 688 PetscCall(MatDestroy(&A)); 689 PetscCall(MatDestroy(&B)); 690 PetscCall(PetscFinalize()); 691 return 0; 692 } 693 694 PetscErrorCode CheckMat(Mat A, Mat B, PetscBool usemult, const char* func) 695 { 696 Mat Bcheck; 697 PetscReal error; 698 699 PetscFunctionBeginUser; 700 if (!usemult && B) { 701 PetscBool hasnorm; 702 703 PetscCall(MatHasOperation(B,MATOP_NORM,&hasnorm)); 704 if (!hasnorm) usemult = PETSC_TRUE; 705 } 706 if (!usemult) { 707 if (B) { 708 MatType Btype; 709 710 PetscCall(MatGetType(B,&Btype)); 711 PetscCall(MatConvert(A,Btype,MAT_INITIAL_MATRIX,&Bcheck)); 712 } else { 713 PetscCall(MatConvert(A,MATAIJ,MAT_INITIAL_MATRIX,&Bcheck)); 714 } 715 if (B) { /* if B is present, subtract it */ 716 PetscCall(MatAXPY(Bcheck,-1.,B,DIFFERENT_NONZERO_PATTERN)); 717 } 718 PetscCall(MatNorm(Bcheck,NORM_INFINITY,&error)); 719 if (error > PETSC_SQRT_MACHINE_EPSILON) { 720 ISLocalToGlobalMapping rl2g,cl2g; 721 722 PetscCall(PetscObjectSetName((PetscObject)Bcheck,"Assembled Bcheck")); 723 PetscCall(MatView(Bcheck,NULL)); 724 if (B) { 725 PetscCall(PetscObjectSetName((PetscObject)B,"Assembled AIJ")); 726 PetscCall(MatView(B,NULL)); 727 PetscCall(MatDestroy(&Bcheck)); 728 PetscCall(MatConvert(A,MATAIJ,MAT_INITIAL_MATRIX,&Bcheck)); 729 PetscCall(PetscObjectSetName((PetscObject)Bcheck,"Assembled IS")); 730 PetscCall(MatView(Bcheck,NULL)); 731 } 732 PetscCall(MatDestroy(&Bcheck)); 733 PetscCall(PetscObjectSetName((PetscObject)A,"MatIS")); 734 PetscCall(MatView(A,NULL)); 735 PetscCall(MatGetLocalToGlobalMapping(A,&rl2g,&cl2g)); 736 PetscCall(ISLocalToGlobalMappingView(rl2g,NULL)); 737 PetscCall(ISLocalToGlobalMappingView(cl2g,NULL)); 738 SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"ERROR ON %s: %g",func,(double)error); 739 } 740 PetscCall(MatDestroy(&Bcheck)); 741 } else { 742 PetscBool ok,okt; 743 744 PetscCall(MatMultEqual(A,B,3,&ok)); 745 PetscCall(MatMultTransposeEqual(A,B,3,&okt)); 746 PetscCheck(ok && okt,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"ERROR ON %s: mult ok ? %d, multtranspose ok ? %d",func,ok,okt); 747 } 748 PetscFunctionReturn(0); 749 } 750 751 PetscErrorCode TestMatZeroRows(Mat A, Mat Afull, PetscBool squaretest, IS is, PetscScalar diag) 752 { 753 Mat B,Bcheck,B2 = NULL,lB; 754 Vec x = NULL, b = NULL, b2 = NULL; 755 ISLocalToGlobalMapping l2gr,l2gc; 756 PetscReal error; 757 char diagstr[16]; 758 const PetscInt *idxs; 759 PetscInt rst,ren,i,n,N,d; 760 PetscMPIInt rank; 761 PetscBool miss,haszerorows; 762 763 PetscFunctionBeginUser; 764 if (diag == 0.) { 765 PetscCall(PetscStrcpy(diagstr,"zero")); 766 } else { 767 PetscCall(PetscStrcpy(diagstr,"nonzero")); 768 } 769 PetscCall(ISView(is,NULL)); 770 PetscCall(MatGetLocalToGlobalMapping(A,&l2gr,&l2gc)); 771 /* tests MatDuplicate and MatCopy */ 772 if (diag == 0.) { 773 PetscCall(MatDuplicate(A,MAT_COPY_VALUES,&B)); 774 } else { 775 PetscCall(MatDuplicate(A,MAT_DO_NOT_COPY_VALUES,&B)); 776 PetscCall(MatCopy(A,B,SAME_NONZERO_PATTERN)); 777 } 778 PetscCall(MatISGetLocalMat(B,&lB)); 779 PetscCall(MatHasOperation(lB,MATOP_ZERO_ROWS,&haszerorows)); 780 if (squaretest && haszerorows) { 781 782 PetscCall(MatCreateVecs(B,&x,&b)); 783 PetscCall(MatDuplicate(B,MAT_COPY_VALUES,&B2)); 784 PetscCall(VecSetLocalToGlobalMapping(b,l2gr)); 785 PetscCall(VecSetLocalToGlobalMapping(x,l2gc)); 786 PetscCall(VecSetRandom(x,NULL)); 787 PetscCall(VecSetRandom(b,NULL)); 788 /* mimic b[is] = x[is] */ 789 PetscCall(VecDuplicate(b,&b2)); 790 PetscCall(VecSetLocalToGlobalMapping(b2,l2gr)); 791 PetscCall(VecCopy(b,b2)); 792 PetscCall(ISGetLocalSize(is,&n)); 793 PetscCall(ISGetIndices(is,&idxs)); 794 PetscCall(VecGetSize(x,&N)); 795 for (i=0;i<n;i++) { 796 if (0 <= idxs[i] && idxs[i] < N) { 797 PetscCall(VecSetValue(b2,idxs[i],diag,INSERT_VALUES)); 798 PetscCall(VecSetValue(x,idxs[i],1.,INSERT_VALUES)); 799 } 800 } 801 PetscCall(VecAssemblyBegin(b2)); 802 PetscCall(VecAssemblyEnd(b2)); 803 PetscCall(VecAssemblyBegin(x)); 804 PetscCall(VecAssemblyEnd(x)); 805 PetscCall(ISRestoreIndices(is,&idxs)); 806 /* test ZeroRows on MATIS */ 807 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Test MatZeroRows (diag %s)\n",diagstr)); 808 PetscCall(MatZeroRowsIS(B,is,diag,x,b)); 809 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Test MatZeroRowsColumns (diag %s)\n",diagstr)); 810 PetscCall(MatZeroRowsColumnsIS(B2,is,diag,NULL,NULL)); 811 } else if (haszerorows) { 812 /* test ZeroRows on MATIS */ 813 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Test MatZeroRows (diag %s)\n",diagstr)); 814 PetscCall(MatZeroRowsIS(B,is,diag,NULL,NULL)); 815 b = b2 = x = NULL; 816 } else { 817 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Skipping MatZeroRows (diag %s)\n",diagstr)); 818 b = b2 = x = NULL; 819 } 820 821 if (squaretest && haszerorows) { 822 PetscCall(VecAXPY(b2,-1.,b)); 823 PetscCall(VecNorm(b2,NORM_INFINITY,&error)); 824 PetscCheck(error <= PETSC_SQRT_MACHINE_EPSILON,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"ERROR IN ZEROROWS ON B %g (diag %s)",(double)error,diagstr); 825 } 826 827 /* test MatMissingDiagonal */ 828 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Test MatMissingDiagonal\n")); 829 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank)); 830 PetscCall(MatMissingDiagonal(B,&miss,&d)); 831 PetscCall(MatGetOwnershipRange(B,&rst,&ren)); 832 PetscCall(PetscViewerASCIIPushSynchronized(PETSC_VIEWER_STDOUT_WORLD)); 833 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)); 834 PetscCall(PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD)); 835 PetscCall(PetscViewerASCIIPopSynchronized(PETSC_VIEWER_STDOUT_WORLD)); 836 837 PetscCall(VecDestroy(&x)); 838 PetscCall(VecDestroy(&b)); 839 PetscCall(VecDestroy(&b2)); 840 841 /* check the result of ZeroRows with that from MPIAIJ routines 842 assuming that MatConvert_IS_XAIJ and MatZeroRows_MPIAIJ work fine */ 843 if (haszerorows) { 844 PetscCall(MatDuplicate(Afull,MAT_COPY_VALUES,&Bcheck)); 845 PetscCall(MatSetOption(Bcheck,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE)); 846 PetscCall(MatZeroRowsIS(Bcheck,is,diag,NULL,NULL)); 847 PetscCall(CheckMat(B,Bcheck,PETSC_FALSE,"Zerorows")); 848 PetscCall(MatDestroy(&Bcheck)); 849 } 850 PetscCall(MatDestroy(&B)); 851 852 if (B2) { /* test MatZeroRowsColumns */ 853 PetscCall(MatDuplicate(Afull,MAT_COPY_VALUES,&B)); 854 PetscCall(MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE)); 855 PetscCall(MatZeroRowsColumnsIS(B,is,diag,NULL,NULL)); 856 PetscCall(CheckMat(B2,B,PETSC_FALSE,"MatZeroRowsColumns")); 857 PetscCall(MatDestroy(&B)); 858 PetscCall(MatDestroy(&B2)); 859 } 860 PetscFunctionReturn(0); 861 } 862 863 /*TEST 864 865 test: 866 args: -test_trans 867 868 test: 869 suffix: 2 870 nsize: 4 871 args: -matis_convert_local_nest -nr 3 -nc 4 872 873 test: 874 suffix: 3 875 nsize: 5 876 args: -m 11 -n 10 -matis_convert_local_nest -nr 2 -nc 1 877 878 test: 879 suffix: 4 880 nsize: 6 881 args: -m 9 -n 12 -test_trans -nr 2 -nc 7 882 883 test: 884 suffix: 5 885 nsize: 6 886 args: -m 12 -n 12 -test_trans -nr 3 -nc 1 887 888 test: 889 suffix: 6 890 args: -m 12 -n 12 -test_trans -nr 2 -nc 3 -diffmap 891 892 test: 893 suffix: 7 894 args: -m 12 -n 12 -test_trans -nr 2 -nc 3 -diffmap -permmap 895 896 test: 897 suffix: 8 898 args: -m 12 -n 17 -test_trans -nr 2 -nc 3 -permmap 899 900 test: 901 suffix: 9 902 nsize: 5 903 args: -m 12 -n 12 -test_trans -nr 2 -nc 3 -diffmap 904 905 test: 906 suffix: 10 907 nsize: 5 908 args: -m 12 -n 12 -test_trans -nr 2 -nc 3 -diffmap -permmap 909 910 test: 911 suffix: vscat_default 912 nsize: 5 913 args: -m 12 -n 17 -test_trans -nr 2 -nc 3 -permmap 914 output_file: output/ex23_11.out 915 916 test: 917 suffix: 12 918 nsize: 3 919 args: -m 12 -n 12 -symmetric -matis_localmat_type sbaij -test_trans -nr 2 -nc 3 920 921 testset: 922 output_file: output/ex23_13.out 923 nsize: 3 924 args: -m 12 -n 17 -test_trans -nr 2 -nc 3 -diffmap -permmap 925 filter: grep -v "type:" 926 test: 927 suffix: baij 928 args: -matis_localmat_type baij 929 test: 930 requires: viennacl 931 suffix: viennacl 932 args: -matis_localmat_type aijviennacl 933 test: 934 requires: cuda 935 suffix: cusparse 936 args: -matis_localmat_type aijcusparse 937 938 test: 939 suffix: negrep 940 nsize: {{1 3}separate output} 941 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} 942 943 TEST*/ 944