1 static char help[] = "A test for MatAssembly that heavily relies on PetscSortIntWithArrayPair\n"; 2 3 /* 4 The characteristic of the array (about 4M in length) to sort in this test is that it has 5 many duplicated values that already clustered together (around 95 duplicates per unique integer). 6 7 It was gotten from a petsc performance bug report from the Moose project. One can use 8 it for future performance study. 9 10 Contributed-by: Fande Kong <fdkong.jd@gmail.com>, John Peterson <jwpeterson@gmail.com> 11 */ 12 13 // PETSc includes 14 #include <petscmat.h> 15 16 // C++ includes 17 #include <iostream> 18 #include <fstream> 19 #include <sstream> 20 #include <algorithm> 21 #include <vector> 22 #include <string> 23 #include <set> 24 25 int main(int argc, char **argv) { 26 PetscMPIInt size, rank; 27 char file[2][PETSC_MAX_PATH_LEN]; 28 PetscBool flg; 29 const unsigned int n_dofs = 26559; 30 unsigned int first_local_index; 31 unsigned int last_local_index; 32 33 PetscFunctionBeginUser; 34 PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 35 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 36 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank)); 37 PetscCheck(size <= 2, PETSC_COMM_WORLD, PETSC_ERR_SUP, "This example is for <=2 procs"); 38 39 PetscCall(PetscOptionsGetString(NULL, NULL, "-f0", file[0], sizeof(file[0]), &flg)); 40 PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_USER, "Must indicate dof indices file for rank 0 with -f0 option"); 41 if (size == 2) { 42 PetscCall(PetscOptionsGetString(NULL, NULL, "-f1", file[1], sizeof(file[1]), &flg)); 43 PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_USER, "Must indicate dof indices file for rank 1 with -f1 option"); 44 } 45 46 if (size == 1) { 47 first_local_index = 0; 48 last_local_index = 26559; 49 } else { 50 if (rank == 0) { 51 first_local_index = 0; 52 last_local_index = 13911; 53 } else { 54 first_local_index = 13911; 55 last_local_index = 26559; 56 } 57 } 58 59 // Read element dof indices from files 60 std::vector<std::vector<std::vector<PetscInt>>> elem_dof_indices(size); 61 for (PetscInt proc_id = 0; proc_id < size; ++proc_id) { 62 std::string line; 63 std::ifstream dof_file(file[proc_id]); 64 if (dof_file.good()) { 65 while (std::getline(dof_file, line)) { 66 std::vector<PetscInt> dof_indices; 67 std::stringstream sstream(line); 68 std::string token; 69 while (std::getline(sstream, token, ' ')) { dof_indices.push_back(std::atoi(token.c_str())); } 70 elem_dof_indices[proc_id].push_back(dof_indices); 71 } 72 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_OPEN, "Could not open file %s", file[proc_id]); 73 } 74 75 // Debugging: Verify we read in elem_dof_indices correctly 76 // for (unsigned int i=0; i<elem_dof_indices.size(); ++i) 77 // { 78 // for (unsigned int j=0; j<elem_dof_indices[i].size(); ++j) 79 // { 80 // for (unsigned int k=0; k<elem_dof_indices[i][j].size(); ++k) 81 // std::cout << elem_dof_indices[i][j][k] << " "; 82 // std::cout << std::endl; 83 // } 84 // std::cout << std::endl; 85 // } 86 87 // Set up the (global) sparsity pattern 88 std::vector<std::set<unsigned int>> sparsity(n_dofs); 89 for (PetscInt proc_id = 0; proc_id < size; ++proc_id) 90 for (unsigned int k = 0; k < elem_dof_indices[proc_id].size(); k++) { 91 std::vector<PetscInt> &dof_indices = elem_dof_indices[proc_id][k]; 92 for (unsigned int i = 0; i < dof_indices.size(); ++i) 93 for (unsigned int j = 0; j < dof_indices.size(); ++j) sparsity[dof_indices[i]].insert(dof_indices[j]); 94 } 95 96 // Determine the local nonzeros on this processor 97 const unsigned int n_local_dofs = last_local_index - first_local_index; 98 std::vector<PetscInt> n_nz(n_local_dofs); 99 std::vector<PetscInt> n_oz(n_local_dofs); 100 101 for (unsigned int i = 0; i < n_local_dofs; ++i) { 102 for (std::set<unsigned int>::iterator iter = sparsity[i + first_local_index].begin(); iter != sparsity[i + first_local_index].end(); iter++) { 103 unsigned int dof = *iter; 104 if ((dof >= first_local_index) && (dof < last_local_index)) n_nz[i]++; 105 else n_oz[i]++; 106 } 107 } 108 109 // Debugging: print number of on/off diagonal nonzeros 110 // for (unsigned int i=0; i<n_nz.size(); ++i) 111 // std::cout << n_nz[i] << " "; 112 // std::cout << std::endl; 113 114 // for (unsigned int i=0; i<n_oz.size(); ++i) 115 // std::cout << n_oz[i] << " "; 116 // std::cout << std::endl; 117 118 // Compute and print max number of on- and off-diagonal nonzeros. 119 PetscInt n_nz_max = *std::max_element(n_nz.begin(), n_nz.end()); 120 PetscInt n_oz_max = *std::max_element(n_oz.begin(), n_oz.end()); 121 122 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "Max on-diagonal non-zeros: = %" PetscInt_FMT "\n", n_nz_max)); 123 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "Max off-diagonal non-zeros: = %" PetscInt_FMT "\n", n_oz_max)); 124 PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT)); 125 126 // Initialize the matrix similar to what we do in the PetscMatrix 127 // ctor and init() routines. 128 Mat mat; 129 PetscCall(MatCreate(PETSC_COMM_WORLD, &mat)); 130 PetscCall(MatSetSizes(mat, n_local_dofs, n_local_dofs, n_dofs, n_dofs)); 131 PetscCall(MatSetBlockSize(mat, 1)); 132 PetscCall(MatSetType(mat, MATAIJ)); // Automatically chooses seqaij or mpiaij 133 PetscCall(MatSeqAIJSetPreallocation(mat, 0, &n_nz[0])); 134 PetscCall(MatMPIAIJSetPreallocation(mat, 0, &n_nz[0], 0, &n_oz[0])); 135 PetscCall(MatSetOption(mat, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE)); 136 137 // Local "element" loop 138 for (unsigned int k = 0; k < elem_dof_indices[rank].size(); k++) { 139 std::vector<PetscInt> &dof_indices = elem_dof_indices[rank][k]; 140 // DenseMatrix< Number > zero_mat( dof_indices.size(), dof_indices.size()); 141 // B.add_matrix( zero_mat, dof_indices); 142 std::vector<PetscScalar> ones(dof_indices.size() * dof_indices.size(), 1.); 143 PetscCall(MatSetValues(mat, dof_indices.size(), &dof_indices[0], dof_indices.size(), &dof_indices[0], &ones[0], ADD_VALUES)); 144 } 145 146 // Matrix assembly 147 PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY)); 148 PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY)); 149 150 // Clean up 151 PetscCall(MatDestroy(&mat)); 152 PetscCall(PetscFinalize()); 153 return 0; 154 } 155 /*TEST 156 build: 157 requires: !defined(PETSC_HAVE_SUN_CXX) 158 159 test: 160 nsize: 2 161 args: -f0 ${DATAFILESPATH}/meshes/moose_dof_indices_np_2_proc_0.txt -f1 ${DATAFILESPATH}/meshes/moose_dof_indices_np_2_proc_1.txt 162 # Skip the test for Sun C++ compiler because of its warnings/errors in petscmat.h 163 requires: datafilespath 164 165 TEST*/ 166