1 2 static char help[] = "Test interface of Elemental. \n\n"; 3 4 #include <petscmat.h> 5 6 int main(int argc, char **args) 7 { 8 Mat C, Caij; 9 PetscInt i, j, m = 5, n, nrows, ncols; 10 const PetscInt *rows, *cols; 11 IS isrows, iscols; 12 PetscBool flg, Test_MatMatMult = PETSC_FALSE, mats_view = PETSC_FALSE; 13 PetscScalar *v; 14 PetscMPIInt rank, size; 15 16 PetscFunctionBeginUser; 17 PetscCall(PetscInitialize(&argc, &args, (char *)0, help)); 18 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank)); 19 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 20 PetscCall(PetscOptionsHasName(NULL, NULL, "-mats_view", &mats_view)); 21 22 /* Get local block or element size*/ 23 PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL)); 24 n = m; 25 PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL)); 26 27 PetscCall(MatCreate(PETSC_COMM_WORLD, &C)); 28 PetscCall(MatSetSizes(C, m, n, PETSC_DECIDE, PETSC_DECIDE)); 29 PetscCall(MatSetType(C, MATELEMENTAL)); 30 PetscCall(MatSetFromOptions(C)); 31 PetscCall(MatSetUp(C)); 32 33 PetscCall(PetscOptionsHasName(NULL, NULL, "-row_oriented", &flg)); 34 if (flg) PetscCall(MatSetOption(C, MAT_ROW_ORIENTED, PETSC_TRUE)); 35 PetscCall(MatGetOwnershipIS(C, &isrows, &iscols)); 36 PetscCall(PetscOptionsHasName(NULL, NULL, "-Cexp_view_ownership", &flg)); 37 if (flg) { /* View ownership of explicit C */ 38 IS tmp; 39 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Ownership of explicit C:\n")); 40 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Row index set:\n")); 41 PetscCall(ISOnComm(isrows, PETSC_COMM_WORLD, PETSC_USE_POINTER, &tmp)); 42 PetscCall(ISView(tmp, PETSC_VIEWER_STDOUT_WORLD)); 43 PetscCall(ISDestroy(&tmp)); 44 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Column index set:\n")); 45 PetscCall(ISOnComm(iscols, PETSC_COMM_WORLD, PETSC_USE_POINTER, &tmp)); 46 PetscCall(ISView(tmp, PETSC_VIEWER_STDOUT_WORLD)); 47 PetscCall(ISDestroy(&tmp)); 48 } 49 50 /* Set local matrix entries */ 51 PetscCall(ISGetLocalSize(isrows, &nrows)); 52 PetscCall(ISGetIndices(isrows, &rows)); 53 PetscCall(ISGetLocalSize(iscols, &ncols)); 54 PetscCall(ISGetIndices(iscols, &cols)); 55 PetscCall(PetscMalloc1(nrows * ncols, &v)); 56 for (i = 0; i < nrows; i++) { 57 for (j = 0; j < ncols; j++) { 58 /*v[i*ncols+j] = (PetscReal)(rank);*/ 59 v[i * ncols + j] = (PetscReal)(rank * 10000 + 100 * rows[i] + cols[j]); 60 } 61 } 62 PetscCall(MatSetValues(C, nrows, rows, ncols, cols, v, INSERT_VALUES)); 63 PetscCall(ISRestoreIndices(isrows, &rows)); 64 PetscCall(ISRestoreIndices(iscols, &cols)); 65 PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY)); 66 PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY)); 67 68 /* Test MatView() */ 69 if (mats_view) PetscCall(MatView(C, PETSC_VIEWER_STDOUT_WORLD)); 70 71 /* Test MatMissingDiagonal() */ 72 PetscCall(MatMissingDiagonal(C, &flg, NULL)); 73 PetscCheck(!flg, PETSC_COMM_SELF, PETSC_ERR_ARG_NOTSAMETYPE, "MatMissingDiagonal() did not return false!"); 74 75 /* Set unowned matrix entries - add subdiagonals and diagonals from proc[0] */ 76 if (rank == 0) { 77 PetscInt M, N, cols[2]; 78 PetscCall(MatGetSize(C, &M, &N)); 79 for (i = 0; i < M; i++) { 80 cols[0] = i; 81 v[0] = i + 0.5; 82 cols[1] = i - 1; 83 v[1] = 0.5; 84 if (i) { 85 PetscCall(MatSetValues(C, 1, &i, 2, cols, v, ADD_VALUES)); 86 } else { 87 PetscCall(MatSetValues(C, 1, &i, 1, &i, v, ADD_VALUES)); 88 } 89 } 90 } 91 PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY)); 92 PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY)); 93 94 /* Test MatMult() */ 95 PetscCall(MatComputeOperator(C, MATAIJ, &Caij)); 96 PetscCall(MatMultEqual(C, Caij, 5, &flg)); 97 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_ARG_NOTSAMETYPE, "C != Caij. MatMultEqual() fails"); 98 PetscCall(MatMultTransposeEqual(C, Caij, 5, &flg)); 99 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_ARG_NOTSAMETYPE, "C != Caij. MatMultTransposeEqual() fails"); 100 101 /* Test MatMultAdd() and MatMultTransposeAddEqual() */ 102 PetscCall(MatMultAddEqual(C, Caij, 5, &flg)); 103 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_ARG_NOTSAMETYPE, "C != Caij. MatMultAddEqual() fails"); 104 PetscCall(MatMultTransposeAddEqual(C, Caij, 5, &flg)); 105 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_ARG_NOTSAMETYPE, "C != Caij. MatMultTransposeAddEqual() fails"); 106 107 /* Test MatMatMult() */ 108 PetscCall(PetscOptionsHasName(NULL, NULL, "-test_matmatmult", &Test_MatMatMult)); 109 if (Test_MatMatMult) { 110 Mat CCelem, CCaij; 111 PetscCall(MatMatMult(C, C, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &CCelem)); 112 PetscCall(MatMatMult(Caij, Caij, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &CCaij)); 113 PetscCall(MatMultEqual(CCelem, CCaij, 5, &flg)); 114 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_ARG_NOTSAMETYPE, "CCelem != CCaij. MatMatMult() fails"); 115 PetscCall(MatDestroy(&CCaij)); 116 PetscCall(MatDestroy(&CCelem)); 117 } 118 119 PetscCall(PetscFree(v)); 120 PetscCall(ISDestroy(&isrows)); 121 PetscCall(ISDestroy(&iscols)); 122 PetscCall(MatDestroy(&C)); 123 PetscCall(MatDestroy(&Caij)); 124 PetscCall(PetscFinalize()); 125 return 0; 126 } 127 128 /*TEST 129 130 test: 131 nsize: 2 132 requires: elemental 133 args: -mat_type elemental -m 2 -n 3 134 135 test: 136 suffix: 2 137 nsize: 6 138 requires: elemental 139 args: -mat_type elemental -m 2 -n 2 140 141 test: 142 suffix: 3 143 nsize: 6 144 requires: elemental 145 args: -mat_type elemental -m 2 -n 2 -test_matmatmult 146 output_file: output/ex38_2.out 147 148 TEST*/ 149