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