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