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