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