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, (PetscInt)info.nz_allocated, 139 (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, (PetscInt)info.nz_allocated, (PetscInt)info.nz_unneeded, 143 (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, (PetscInt)info.nz_allocated, (PetscInt)info.nz_unneeded, 146 (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++) PetscCall(ISCreateStride(PETSC_COMM_WORLD, ren - rst, i + rst, nr, &rows[i])); 391 for (i = 0; i < nc; i++) PetscCall(ISCreateStride(PETSC_COMM_WORLD, cen - cst, i + cst, nc, &cols[i])); 392 PetscCall(MatCreateNest(PETSC_COMM_WORLD, nr, rows, nc, cols, mats, &A2)); 393 PetscCall(MatCreateNest(PETSC_COMM_WORLD, nr, rows, nc, cols, mats + nr * nc, &B2)); 394 for (i = 0; i < nr; i++) PetscCall(ISDestroy(&rows[i])); 395 for (i = 0; i < nc; i++) PetscCall(ISDestroy(&cols[i])); 396 for (i = 0; i < 2 * nr * nc; i++) PetscCall(MatDestroy(&mats[i])); 397 PetscCall(PetscFree3(rows, cols, mats)); 398 PetscCall(MatConvert(B2, MATAIJ, MAT_INITIAL_MATRIX, &T)); 399 PetscCall(MatDestroy(&B2)); 400 PetscCall(MatConvert(A2, MATIS, MAT_INITIAL_MATRIX, &B2)); 401 PetscCall(CheckMat(B2, T, PETSC_TRUE, "MatConvert_Nest_IS MAT_INITIAL_MATRIX")); 402 PetscCall(MatConvert(A2, MATIS, MAT_REUSE_MATRIX, &B2)); 403 PetscCall(CheckMat(B2, T, PETSC_TRUE, "MatConvert_Nest_IS MAT_REUSE_MATRIX")); 404 PetscCall(MatDestroy(&B2)); 405 PetscCall(MatConvert(A2, MATIS, MAT_INPLACE_MATRIX, &A2)); 406 PetscCall(CheckMat(A2, T, PETSC_TRUE, "MatConvert_Nest_IS MAT_INPLACE_MATRIX")); 407 PetscCall(MatDestroy(&T)); 408 PetscCall(MatDestroy(&A2)); 409 410 /* test MatCreateSubMatrix */ 411 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatCreateSubMatrix\n")); 412 if (rank == 0) { 413 PetscCall(ISCreateStride(PETSC_COMM_WORLD, 1, 1, 1, &is)); 414 PetscCall(ISCreateStride(PETSC_COMM_WORLD, 2, 0, 1, &is2)); 415 } else if (rank == 1) { 416 PetscCall(ISCreateStride(PETSC_COMM_WORLD, 1, 0, 1, &is)); 417 if (n > 3) { 418 PetscCall(ISCreateStride(PETSC_COMM_WORLD, 1, 3, 1, &is2)); 419 } else { 420 PetscCall(ISCreateStride(PETSC_COMM_WORLD, 0, 0, 1, &is2)); 421 } 422 } else if (rank == 2 && n > 4) { 423 PetscCall(ISCreateStride(PETSC_COMM_WORLD, 0, 0, 1, &is)); 424 PetscCall(ISCreateStride(PETSC_COMM_WORLD, n - 4, 4, 1, &is2)); 425 } else { 426 PetscCall(ISCreateStride(PETSC_COMM_WORLD, 0, 0, 1, &is)); 427 PetscCall(ISCreateStride(PETSC_COMM_WORLD, 0, 0, 1, &is2)); 428 } 429 PetscCall(MatCreateSubMatrix(A, is, is, MAT_INITIAL_MATRIX, &A2)); 430 PetscCall(MatCreateSubMatrix(B, is, is, MAT_INITIAL_MATRIX, &B2)); 431 PetscCall(CheckMat(A2, B2, PETSC_TRUE, "first MatCreateSubMatrix")); 432 433 PetscCall(MatCreateSubMatrix(A, is, is, MAT_REUSE_MATRIX, &A2)); 434 PetscCall(MatCreateSubMatrix(B, is, is, MAT_REUSE_MATRIX, &B2)); 435 PetscCall(CheckMat(A2, B2, PETSC_FALSE, "reuse MatCreateSubMatrix")); 436 PetscCall(MatDestroy(&A2)); 437 PetscCall(MatDestroy(&B2)); 438 439 if (!issymmetric) { 440 PetscCall(MatCreateSubMatrix(A, is, is2, MAT_INITIAL_MATRIX, &A2)); 441 PetscCall(MatCreateSubMatrix(B, is, is2, MAT_INITIAL_MATRIX, &B2)); 442 PetscCall(MatCreateSubMatrix(A, is, is2, MAT_REUSE_MATRIX, &A2)); 443 PetscCall(MatCreateSubMatrix(B, is, is2, MAT_REUSE_MATRIX, &B2)); 444 PetscCall(CheckMat(A2, B2, PETSC_FALSE, "second MatCreateSubMatrix")); 445 } 446 447 PetscCall(MatDestroy(&A2)); 448 PetscCall(MatDestroy(&B2)); 449 PetscCall(ISDestroy(&is)); 450 PetscCall(ISDestroy(&is2)); 451 452 /* test MatCreateSubMatrices */ 453 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatCreateSubMatrices\n")); 454 PetscCall(MatGetLayouts(A, &rlayout, &clayout)); 455 PetscCall(PetscLayoutGetRanges(rlayout, &rrange)); 456 PetscCall(PetscLayoutGetRanges(clayout, &crange)); 457 lrank = (size + rank - 1) % size; 458 rrank = (rank + 1) % size; 459 PetscCall(ISCreateStride(PETSC_COMM_SELF, (rrange[lrank + 1] - rrange[lrank]), rrange[lrank], 1, &irow[0])); 460 PetscCall(ISCreateStride(PETSC_COMM_SELF, (crange[rrank + 1] - crange[rrank]), crange[rrank], 1, &icol[0])); 461 PetscCall(ISCreateStride(PETSC_COMM_SELF, (rrange[rrank + 1] - rrange[rrank]), rrange[rrank], 1, &irow[1])); 462 PetscCall(ISCreateStride(PETSC_COMM_SELF, (crange[lrank + 1] - crange[lrank]), crange[lrank], 1, &icol[1])); 463 PetscCall(MatCreateSubMatrices(A, 2, irow, icol, MAT_INITIAL_MATRIX, &Asub)); 464 PetscCall(MatCreateSubMatrices(B, 2, irow, icol, MAT_INITIAL_MATRIX, &Bsub)); 465 PetscCall(CheckMat(Asub[0], Bsub[0], PETSC_FALSE, "MatCreateSubMatrices[0]")); 466 PetscCall(CheckMat(Asub[1], Bsub[1], PETSC_FALSE, "MatCreateSubMatrices[1]")); 467 PetscCall(MatCreateSubMatrices(A, 2, irow, icol, MAT_REUSE_MATRIX, &Asub)); 468 PetscCall(MatCreateSubMatrices(B, 2, irow, icol, MAT_REUSE_MATRIX, &Bsub)); 469 PetscCall(CheckMat(Asub[0], Bsub[0], PETSC_FALSE, "MatCreateSubMatrices[0]")); 470 PetscCall(CheckMat(Asub[1], Bsub[1], PETSC_FALSE, "MatCreateSubMatrices[1]")); 471 PetscCall(MatDestroySubMatrices(2, &Asub)); 472 PetscCall(MatDestroySubMatrices(2, &Bsub)); 473 PetscCall(ISDestroy(&irow[0])); 474 PetscCall(ISDestroy(&irow[1])); 475 PetscCall(ISDestroy(&icol[0])); 476 PetscCall(ISDestroy(&icol[1])); 477 478 /* Create an IS required by MatZeroRows(): just rank zero provides the rows to be eliminated */ 479 if (size > 1) { 480 if (rank == 0) { 481 PetscInt st, len; 482 483 st = (m + 1) / 2; 484 len = PetscMin(m / 2, PetscMax(m - (m + 1) / 2 - 1, 0)); 485 PetscCall(ISCreateStride(PETSC_COMM_WORLD, len, st, 1, &is)); 486 } else { 487 PetscCall(ISCreateStride(PETSC_COMM_WORLD, 0, 0, 1, &is)); 488 } 489 } else { 490 PetscCall(ISCreateStride(PETSC_COMM_WORLD, 1, 0, 1, &is)); 491 } 492 493 if (squaretest) { /* tests for square matrices only, with same maps for rows and columns */ 494 PetscInt *idx0, *idx1, n0, n1; 495 IS Ais[2], Bis[2]; 496 497 /* test MatDiagonalSet */ 498 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatDiagonalSet\n")); 499 PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &A2)); 500 PetscCall(MatDuplicate(B, MAT_COPY_VALUES, &B2)); 501 PetscCall(MatCreateVecs(A, NULL, &x)); 502 PetscCall(VecSetRandom(x, NULL)); 503 PetscCall(MatDiagonalSet(A2, x, INSERT_VALUES)); 504 PetscCall(MatDiagonalSet(B2, x, INSERT_VALUES)); 505 PetscCall(CheckMat(A2, B2, PETSC_FALSE, "MatDiagonalSet")); 506 PetscCall(VecDestroy(&x)); 507 PetscCall(MatDestroy(&A2)); 508 PetscCall(MatDestroy(&B2)); 509 510 /* test MatShift (MatShift_IS internally uses MatDiagonalSet_IS with ADD_VALUES) */ 511 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatShift\n")); 512 PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &A2)); 513 PetscCall(MatDuplicate(B, MAT_COPY_VALUES, &B2)); 514 PetscCall(MatShift(A2, 2.0)); 515 PetscCall(MatShift(B2, 2.0)); 516 PetscCall(CheckMat(A2, B2, PETSC_FALSE, "MatShift")); 517 PetscCall(MatDestroy(&A2)); 518 PetscCall(MatDestroy(&B2)); 519 520 /* nonzero diag value is supported for square matrices only */ 521 PetscCall(TestMatZeroRows(A, B, PETSC_TRUE, is, diag)); 522 523 /* test MatIncreaseOverlap */ 524 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatIncreaseOverlap\n")); 525 PetscCall(MatGetOwnershipRange(A, &rst, &ren)); 526 n0 = (ren - rst) / 2; 527 n1 = (ren - rst) / 3; 528 PetscCall(PetscMalloc1(n0, &idx0)); 529 PetscCall(PetscMalloc1(n1, &idx1)); 530 for (i = 0; i < n0; i++) idx0[i] = ren - i - 1; 531 for (i = 0; i < n1; i++) idx1[i] = rst + i; 532 PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, n0, idx0, PETSC_OWN_POINTER, &Ais[0])); 533 PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, n1, idx1, PETSC_OWN_POINTER, &Ais[1])); 534 PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, n0, idx0, PETSC_COPY_VALUES, &Bis[0])); 535 PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, n1, idx1, PETSC_COPY_VALUES, &Bis[1])); 536 PetscCall(MatIncreaseOverlap(A, 2, Ais, 3)); 537 PetscCall(MatIncreaseOverlap(B, 2, Bis, 3)); 538 /* Non deterministic output! */ 539 PetscCall(ISSort(Ais[0])); 540 PetscCall(ISSort(Ais[1])); 541 PetscCall(ISSort(Bis[0])); 542 PetscCall(ISSort(Bis[1])); 543 PetscCall(ISView(Ais[0], NULL)); 544 PetscCall(ISView(Bis[0], NULL)); 545 PetscCall(ISView(Ais[1], NULL)); 546 PetscCall(ISView(Bis[1], NULL)); 547 PetscCall(MatCreateSubMatrices(A, 2, Ais, Ais, MAT_INITIAL_MATRIX, &Asub)); 548 PetscCall(MatCreateSubMatrices(B, 2, Bis, Ais, MAT_INITIAL_MATRIX, &Bsub)); 549 PetscCall(CheckMat(Asub[0], Bsub[0], PETSC_FALSE, "MatIncreaseOverlap[0]")); 550 PetscCall(CheckMat(Asub[1], Bsub[1], PETSC_FALSE, "MatIncreaseOverlap[1]")); 551 PetscCall(MatDestroySubMatrices(2, &Asub)); 552 PetscCall(MatDestroySubMatrices(2, &Bsub)); 553 PetscCall(ISDestroy(&Ais[0])); 554 PetscCall(ISDestroy(&Ais[1])); 555 PetscCall(ISDestroy(&Bis[0])); 556 PetscCall(ISDestroy(&Bis[1])); 557 } 558 PetscCall(TestMatZeroRows(A, B, squaretest, is, 0.0)); 559 PetscCall(ISDestroy(&is)); 560 561 /* test MatTranspose */ 562 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatTranspose\n")); 563 PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &A2)); 564 PetscCall(MatTranspose(B, MAT_INITIAL_MATRIX, &B2)); 565 PetscCall(CheckMat(A2, B2, PETSC_FALSE, "initial matrix MatTranspose")); 566 567 PetscCall(MatTranspose(A, MAT_REUSE_MATRIX, &A2)); 568 PetscCall(CheckMat(A2, B2, PETSC_FALSE, "reuse matrix (not in place) MatTranspose")); 569 PetscCall(MatDestroy(&A2)); 570 571 PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &A2)); 572 PetscCall(MatTranspose(A2, MAT_INPLACE_MATRIX, &A2)); 573 PetscCall(CheckMat(A2, B2, PETSC_FALSE, "reuse matrix (in place) MatTranspose")); 574 PetscCall(MatDestroy(&A2)); 575 576 PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &A2)); 577 PetscCall(CheckMat(A2, B2, PETSC_FALSE, "reuse matrix (different type) MatTranspose")); 578 PetscCall(MatDestroy(&A2)); 579 PetscCall(MatDestroy(&B2)); 580 581 /* test MatISFixLocalEmpty */ 582 if (isaij) { 583 PetscInt r[2]; 584 585 r[0] = 0; 586 r[1] = PetscMin(m, n) - 1; 587 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatISFixLocalEmpty\n")); 588 PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &A2)); 589 590 PetscCall(MatISFixLocalEmpty(A2, PETSC_TRUE)); 591 PetscCall(MatAssemblyBegin(A2, MAT_FINAL_ASSEMBLY)); 592 PetscCall(MatAssemblyEnd(A2, MAT_FINAL_ASSEMBLY)); 593 PetscCall(CheckMat(A2, B, PETSC_FALSE, "MatISFixLocalEmpty (null)")); 594 595 PetscCall(MatZeroRows(A2, 2, r, 0.0, NULL, NULL)); 596 PetscCall(MatViewFromOptions(A2, NULL, "-fixempty_view")); 597 PetscCall(MatDuplicate(B, MAT_COPY_VALUES, &B2)); 598 PetscCall(MatZeroRows(B2, 2, r, 0.0, NULL, NULL)); 599 PetscCall(MatISFixLocalEmpty(A2, PETSC_TRUE)); 600 PetscCall(MatAssemblyBegin(A2, MAT_FINAL_ASSEMBLY)); 601 PetscCall(MatAssemblyEnd(A2, MAT_FINAL_ASSEMBLY)); 602 PetscCall(MatViewFromOptions(A2, NULL, "-fixempty_view")); 603 PetscCall(CheckMat(A2, B2, PETSC_FALSE, "MatISFixLocalEmpty (rows)")); 604 PetscCall(MatDestroy(&A2)); 605 606 PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &A2)); 607 PetscCall(MatZeroRows(A2, 2, r, 0.0, NULL, NULL)); 608 PetscCall(MatTranspose(A2, MAT_INPLACE_MATRIX, &A2)); 609 PetscCall(MatTranspose(B2, MAT_INPLACE_MATRIX, &B2)); 610 PetscCall(MatViewFromOptions(A2, NULL, "-fixempty_view")); 611 PetscCall(MatISFixLocalEmpty(A2, PETSC_TRUE)); 612 PetscCall(MatAssemblyBegin(A2, MAT_FINAL_ASSEMBLY)); 613 PetscCall(MatAssemblyEnd(A2, MAT_FINAL_ASSEMBLY)); 614 PetscCall(MatViewFromOptions(A2, NULL, "-fixempty_view")); 615 PetscCall(CheckMat(A2, B2, PETSC_FALSE, "MatISFixLocalEmpty (cols)")); 616 617 PetscCall(MatDestroy(&A2)); 618 PetscCall(MatDestroy(&B2)); 619 620 if (squaretest) { 621 PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &A2)); 622 PetscCall(MatDuplicate(B, MAT_COPY_VALUES, &B2)); 623 PetscCall(MatZeroRowsColumns(A2, 2, r, 0.0, NULL, NULL)); 624 PetscCall(MatZeroRowsColumns(B2, 2, r, 0.0, NULL, NULL)); 625 PetscCall(MatViewFromOptions(A2, NULL, "-fixempty_view")); 626 PetscCall(MatISFixLocalEmpty(A2, PETSC_TRUE)); 627 PetscCall(MatAssemblyBegin(A2, MAT_FINAL_ASSEMBLY)); 628 PetscCall(MatAssemblyEnd(A2, MAT_FINAL_ASSEMBLY)); 629 PetscCall(MatViewFromOptions(A2, NULL, "-fixempty_view")); 630 PetscCall(CheckMat(A2, B2, PETSC_FALSE, "MatISFixLocalEmpty (rows+cols)")); 631 PetscCall(MatDestroy(&A2)); 632 PetscCall(MatDestroy(&B2)); 633 } 634 } 635 636 /* test MatInvertBlockDiagonal 637 special cases for block-diagonal matrices */ 638 if (m == n) { 639 ISLocalToGlobalMapping map; 640 Mat Abd, Bbd; 641 IS is, bis; 642 const PetscScalar *isbd, *aijbd; 643 PetscScalar *vals; 644 const PetscInt *sts, *idxs; 645 PetscInt *idxs2, diff, perm, nl, bs, st, en, in; 646 PetscBool ok; 647 648 for (diff = 0; diff < 3; diff++) { 649 for (perm = 0; perm < 3; perm++) { 650 for (bs = 1; bs < 4; bs++) { 651 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatInvertBlockDiagonal blockdiag %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n", n, diff, perm, bs)); 652 PetscCall(PetscMalloc1(bs * bs, &vals)); 653 PetscCall(MatGetOwnershipRanges(A, &sts)); 654 switch (diff) { 655 case 1: /* inverted layout by processes */ 656 in = 1; 657 st = sts[size - rank - 1]; 658 en = sts[size - rank]; 659 nl = en - st; 660 break; 661 case 2: /* round-robin layout */ 662 in = size; 663 st = rank; 664 nl = n / size; 665 if (rank < n % size) nl++; 666 break; 667 default: /* same layout */ 668 in = 1; 669 st = sts[rank]; 670 en = sts[rank + 1]; 671 nl = en - st; 672 break; 673 } 674 PetscCall(ISCreateStride(PETSC_COMM_WORLD, nl, st, in, &is)); 675 PetscCall(ISGetLocalSize(is, &nl)); 676 PetscCall(ISGetIndices(is, &idxs)); 677 PetscCall(PetscMalloc1(nl, &idxs2)); 678 for (i = 0; i < nl; i++) { 679 switch (perm) { /* invert some of the indices */ 680 case 2: 681 idxs2[i] = rank % 2 ? idxs[i] : idxs[nl - i - 1]; 682 break; 683 case 1: 684 idxs2[i] = rank % 2 ? idxs[nl - i - 1] : idxs[i]; 685 break; 686 default: 687 idxs2[i] = idxs[i]; 688 break; 689 } 690 } 691 PetscCall(ISRestoreIndices(is, &idxs)); 692 PetscCall(ISCreateBlock(PETSC_COMM_WORLD, bs, nl, idxs2, PETSC_OWN_POINTER, &bis)); 693 PetscCall(ISLocalToGlobalMappingCreateIS(bis, &map)); 694 PetscCall(MatCreateIS(PETSC_COMM_WORLD, bs, PETSC_DECIDE, PETSC_DECIDE, bs * n, bs * n, map, map, &Abd)); 695 PetscCall(ISLocalToGlobalMappingDestroy(&map)); 696 PetscCall(MatISSetPreallocation(Abd, bs, NULL, 0, NULL)); 697 for (i = 0; i < nl; i++) { 698 PetscInt b1, b2; 699 700 for (b1 = 0; b1 < bs; b1++) { 701 for (b2 = 0; b2 < bs; b2++) vals[b1 * bs + b2] = i * bs * bs + b1 * bs + b2 + 1 + (b1 == b2 ? 1.0 : 0); 702 } 703 PetscCall(MatSetValuesBlockedLocal(Abd, 1, &i, 1, &i, vals, INSERT_VALUES)); 704 } 705 PetscCall(MatAssemblyBegin(Abd, MAT_FINAL_ASSEMBLY)); 706 PetscCall(MatAssemblyEnd(Abd, MAT_FINAL_ASSEMBLY)); 707 PetscCall(MatConvert(Abd, MATAIJ, MAT_INITIAL_MATRIX, &Bbd)); 708 PetscCall(MatInvertBlockDiagonal(Abd, &isbd)); 709 PetscCall(MatInvertBlockDiagonal(Bbd, &aijbd)); 710 PetscCall(MatGetLocalSize(Bbd, &nl, NULL)); 711 ok = PETSC_TRUE; 712 for (i = 0; i < nl / bs; i++) { 713 PetscInt b1, b2; 714 715 for (b1 = 0; b1 < bs; b1++) { 716 for (b2 = 0; b2 < bs; b2++) { 717 if (PetscAbsScalar(isbd[i * bs * bs + b1 * bs + b2] - aijbd[i * bs * bs + b1 * bs + b2]) > PETSC_SMALL) ok = PETSC_FALSE; 718 if (!ok) { 719 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]))); 720 break; 721 } 722 } 723 if (!ok) break; 724 } 725 if (!ok) break; 726 } 727 PetscCall(MatDestroy(&Abd)); 728 PetscCall(MatDestroy(&Bbd)); 729 PetscCall(PetscFree(vals)); 730 PetscCall(ISDestroy(&is)); 731 PetscCall(ISDestroy(&bis)); 732 } 733 } 734 } 735 } 736 737 /* test MatGetDiagonalBlock */ 738 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatGetDiagonalBlock\n")); 739 PetscCall(MatGetDiagonalBlock(A, &A2)); 740 PetscCall(MatGetDiagonalBlock(B, &B2)); 741 PetscCall(CheckMat(A2, B2, PETSC_FALSE, "MatGetDiagonalBlock")); 742 PetscCall(MatScale(A, 2.0)); 743 PetscCall(MatScale(B, 2.0)); 744 PetscCall(MatGetDiagonalBlock(A, &A2)); 745 PetscCall(MatGetDiagonalBlock(B, &B2)); 746 PetscCall(CheckMat(A2, B2, PETSC_FALSE, "MatGetDiagonalBlock")); 747 748 /* free testing matrices */ 749 PetscCall(ISLocalToGlobalMappingDestroy(&cmap)); 750 PetscCall(ISLocalToGlobalMappingDestroy(&rmap)); 751 PetscCall(MatDestroy(&A)); 752 PetscCall(MatDestroy(&B)); 753 PetscCall(PetscFinalize()); 754 return 0; 755 } 756 757 PetscErrorCode CheckMat(Mat A, Mat B, PetscBool usemult, const char *func) 758 { 759 Mat Bcheck; 760 PetscReal error; 761 762 PetscFunctionBeginUser; 763 if (!usemult && B) { 764 PetscBool hasnorm; 765 766 PetscCall(MatHasOperation(B, MATOP_NORM, &hasnorm)); 767 if (!hasnorm) usemult = PETSC_TRUE; 768 } 769 if (!usemult) { 770 if (B) { 771 MatType Btype; 772 773 PetscCall(MatGetType(B, &Btype)); 774 PetscCall(MatConvert(A, Btype, MAT_INITIAL_MATRIX, &Bcheck)); 775 } else { 776 PetscCall(MatConvert(A, MATAIJ, MAT_INITIAL_MATRIX, &Bcheck)); 777 } 778 if (B) { /* if B is present, subtract it */ 779 PetscCall(MatAXPY(Bcheck, -1., B, DIFFERENT_NONZERO_PATTERN)); 780 } 781 PetscCall(MatNorm(Bcheck, NORM_INFINITY, &error)); 782 if (error > PETSC_SQRT_MACHINE_EPSILON) { 783 ISLocalToGlobalMapping rl2g, cl2g; 784 785 PetscCall(PetscObjectSetName((PetscObject)Bcheck, "Bcheck")); 786 PetscCall(MatView(Bcheck, NULL)); 787 if (B) { 788 PetscCall(PetscObjectSetName((PetscObject)B, "B")); 789 PetscCall(MatView(B, NULL)); 790 PetscCall(MatDestroy(&Bcheck)); 791 PetscCall(MatConvert(A, MATAIJ, MAT_INITIAL_MATRIX, &Bcheck)); 792 PetscCall(PetscObjectSetName((PetscObject)Bcheck, "Assembled A")); 793 PetscCall(MatView(Bcheck, NULL)); 794 } 795 PetscCall(MatDestroy(&Bcheck)); 796 PetscCall(PetscObjectSetName((PetscObject)A, "A")); 797 PetscCall(MatView(A, NULL)); 798 PetscCall(MatGetLocalToGlobalMapping(A, &rl2g, &cl2g)); 799 if (rl2g) PetscCall(ISLocalToGlobalMappingView(rl2g, NULL)); 800 if (cl2g) PetscCall(ISLocalToGlobalMappingView(cl2g, NULL)); 801 SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "ERROR ON %s: %g", func, (double)error); 802 } 803 PetscCall(MatDestroy(&Bcheck)); 804 } else { 805 PetscBool ok, okt; 806 807 PetscCall(MatMultEqual(A, B, 3, &ok)); 808 PetscCall(MatMultTransposeEqual(A, B, 3, &okt)); 809 PetscCheck(ok && okt, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "ERROR ON %s: mult ok ? %d, multtranspose ok ? %d", func, ok, okt); 810 } 811 PetscFunctionReturn(0); 812 } 813 814 PetscErrorCode TestMatZeroRows(Mat A, Mat Afull, PetscBool squaretest, IS is, PetscScalar diag) 815 { 816 Mat B, Bcheck, B2 = NULL, lB; 817 Vec x = NULL, b = NULL, b2 = NULL; 818 ISLocalToGlobalMapping l2gr, l2gc; 819 PetscReal error; 820 char diagstr[16]; 821 const PetscInt *idxs; 822 PetscInt rst, ren, i, n, N, d; 823 PetscMPIInt rank; 824 PetscBool miss, haszerorows; 825 826 PetscFunctionBeginUser; 827 if (diag == 0.) { 828 PetscCall(PetscStrcpy(diagstr, "zero")); 829 } else { 830 PetscCall(PetscStrcpy(diagstr, "nonzero")); 831 } 832 PetscCall(ISView(is, NULL)); 833 PetscCall(MatGetLocalToGlobalMapping(A, &l2gr, &l2gc)); 834 /* tests MatDuplicate and MatCopy */ 835 if (diag == 0.) { 836 PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &B)); 837 } else { 838 PetscCall(MatDuplicate(A, MAT_DO_NOT_COPY_VALUES, &B)); 839 PetscCall(MatCopy(A, B, SAME_NONZERO_PATTERN)); 840 } 841 PetscCall(MatISGetLocalMat(B, &lB)); 842 PetscCall(MatHasOperation(lB, MATOP_ZERO_ROWS, &haszerorows)); 843 if (squaretest && haszerorows) { 844 PetscCall(MatCreateVecs(B, &x, &b)); 845 PetscCall(MatDuplicate(B, MAT_COPY_VALUES, &B2)); 846 PetscCall(VecSetLocalToGlobalMapping(b, l2gr)); 847 PetscCall(VecSetLocalToGlobalMapping(x, l2gc)); 848 PetscCall(VecSetRandom(x, NULL)); 849 PetscCall(VecSetRandom(b, NULL)); 850 /* mimic b[is] = x[is] */ 851 PetscCall(VecDuplicate(b, &b2)); 852 PetscCall(VecSetLocalToGlobalMapping(b2, l2gr)); 853 PetscCall(VecCopy(b, b2)); 854 PetscCall(ISGetLocalSize(is, &n)); 855 PetscCall(ISGetIndices(is, &idxs)); 856 PetscCall(VecGetSize(x, &N)); 857 for (i = 0; i < n; i++) { 858 if (0 <= idxs[i] && idxs[i] < N) { 859 PetscCall(VecSetValue(b2, idxs[i], diag, INSERT_VALUES)); 860 PetscCall(VecSetValue(x, idxs[i], 1., INSERT_VALUES)); 861 } 862 } 863 PetscCall(VecAssemblyBegin(b2)); 864 PetscCall(VecAssemblyEnd(b2)); 865 PetscCall(VecAssemblyBegin(x)); 866 PetscCall(VecAssemblyEnd(x)); 867 PetscCall(ISRestoreIndices(is, &idxs)); 868 /* test ZeroRows on MATIS */ 869 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatZeroRows (diag %s)\n", diagstr)); 870 PetscCall(MatZeroRowsIS(B, is, diag, x, b)); 871 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatZeroRowsColumns (diag %s)\n", diagstr)); 872 PetscCall(MatZeroRowsColumnsIS(B2, is, diag, NULL, NULL)); 873 } else if (haszerorows) { 874 /* test ZeroRows on MATIS */ 875 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatZeroRows (diag %s)\n", diagstr)); 876 PetscCall(MatZeroRowsIS(B, is, diag, NULL, NULL)); 877 b = b2 = x = NULL; 878 } else { 879 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Skipping MatZeroRows (diag %s)\n", diagstr)); 880 b = b2 = x = NULL; 881 } 882 883 if (squaretest && haszerorows) { 884 PetscCall(VecAXPY(b2, -1., b)); 885 PetscCall(VecNorm(b2, NORM_INFINITY, &error)); 886 PetscCheck(error <= PETSC_SQRT_MACHINE_EPSILON, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "ERROR IN ZEROROWS ON B %g (diag %s)", (double)error, diagstr); 887 } 888 889 /* test MatMissingDiagonal */ 890 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatMissingDiagonal\n")); 891 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank)); 892 PetscCall(MatMissingDiagonal(B, &miss, &d)); 893 PetscCall(MatGetOwnershipRange(B, &rst, &ren)); 894 PetscCall(PetscViewerASCIIPushSynchronized(PETSC_VIEWER_STDOUT_WORLD)); 895 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)); 896 PetscCall(PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD)); 897 PetscCall(PetscViewerASCIIPopSynchronized(PETSC_VIEWER_STDOUT_WORLD)); 898 899 PetscCall(VecDestroy(&x)); 900 PetscCall(VecDestroy(&b)); 901 PetscCall(VecDestroy(&b2)); 902 903 /* check the result of ZeroRows with that from MPIAIJ routines 904 assuming that MatConvert_IS_XAIJ and MatZeroRows_MPIAIJ work fine */ 905 if (haszerorows) { 906 PetscCall(MatDuplicate(Afull, MAT_COPY_VALUES, &Bcheck)); 907 PetscCall(MatSetOption(Bcheck, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE)); 908 PetscCall(MatZeroRowsIS(Bcheck, is, diag, NULL, NULL)); 909 PetscCall(CheckMat(B, Bcheck, PETSC_FALSE, "Zerorows")); 910 PetscCall(MatDestroy(&Bcheck)); 911 } 912 PetscCall(MatDestroy(&B)); 913 914 if (B2) { /* test MatZeroRowsColumns */ 915 PetscCall(MatDuplicate(Afull, MAT_COPY_VALUES, &B)); 916 PetscCall(MatSetOption(B, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE)); 917 PetscCall(MatZeroRowsColumnsIS(B, is, diag, NULL, NULL)); 918 PetscCall(CheckMat(B2, B, PETSC_FALSE, "MatZeroRowsColumns")); 919 PetscCall(MatDestroy(&B)); 920 PetscCall(MatDestroy(&B2)); 921 } 922 PetscFunctionReturn(0); 923 } 924 925 /*TEST 926 927 test: 928 args: -test_trans 929 930 test: 931 suffix: 2 932 nsize: 4 933 args: -matis_convert_local_nest -nr 3 -nc 4 934 935 test: 936 suffix: 3 937 nsize: 5 938 args: -m 11 -n 10 -matis_convert_local_nest -nr 2 -nc 1 939 940 test: 941 suffix: 4 942 nsize: 6 943 args: -m 9 -n 12 -test_trans -nr 2 -nc 7 944 945 test: 946 suffix: 5 947 nsize: 6 948 args: -m 12 -n 12 -test_trans -nr 3 -nc 1 949 950 test: 951 suffix: 6 952 args: -m 12 -n 12 -test_trans -nr 2 -nc 3 -diffmap 953 954 test: 955 suffix: 7 956 args: -m 12 -n 12 -test_trans -nr 2 -nc 3 -diffmap -permmap 957 958 test: 959 suffix: 8 960 args: -m 12 -n 17 -test_trans -nr 2 -nc 3 -permmap 961 962 test: 963 suffix: 9 964 nsize: 5 965 args: -m 12 -n 12 -test_trans -nr 2 -nc 3 -diffmap 966 967 test: 968 suffix: 10 969 nsize: 5 970 args: -m 12 -n 12 -test_trans -nr 2 -nc 3 -diffmap -permmap 971 972 test: 973 suffix: vscat_default 974 nsize: 5 975 args: -m 12 -n 17 -test_trans -nr 2 -nc 3 -permmap 976 output_file: output/ex23_11.out 977 978 test: 979 suffix: 12 980 nsize: 3 981 args: -m 12 -n 12 -symmetric -matis_localmat_type sbaij -test_trans -nr 2 -nc 3 982 983 testset: 984 output_file: output/ex23_13.out 985 nsize: 3 986 args: -m 12 -n 17 -test_trans -nr 2 -nc 3 -diffmap -permmap 987 filter: grep -v "type:" 988 test: 989 suffix: baij 990 args: -matis_localmat_type baij 991 test: 992 requires: viennacl 993 suffix: viennacl 994 args: -matis_localmat_type aijviennacl 995 test: 996 requires: cuda 997 suffix: cusparse 998 args: -matis_localmat_type aijcusparse 999 1000 test: 1001 suffix: negrep 1002 nsize: {{1 3}separate output} 1003 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} 1004 1005 TEST*/ 1006