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