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
main(int argc,char ** argv)25 int main(int argc, char **argv)
26 {
27 PetscMPIInt size, rank;
28 char file[2][PETSC_MAX_PATH_LEN];
29 PetscBool flg;
30 const unsigned int n_dofs = 26559;
31 unsigned int first_local_index;
32 unsigned int last_local_index;
33
34 PetscFunctionBeginUser;
35 PetscCall(PetscInitialize(&argc, &argv, nullptr, help));
36 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
37 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
38 PetscCheck(size <= 2, PETSC_COMM_WORLD, PETSC_ERR_SUP, "This example is for <=2 procs");
39
40 PetscCall(PetscOptionsGetString(NULL, NULL, "-f0", file[0], sizeof(file[0]), &flg));
41 PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_USER, "Must indicate dof indices file for rank 0 with -f0 option");
42 if (size == 2) {
43 PetscCall(PetscOptionsGetString(NULL, NULL, "-f1", file[1], sizeof(file[1]), &flg));
44 PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_USER, "Must indicate dof indices file for rank 1 with -f1 option");
45 }
46
47 if (size == 1) {
48 first_local_index = 0;
49 last_local_index = 26559;
50 } else {
51 if (rank == 0) {
52 first_local_index = 0;
53 last_local_index = 13911;
54 } else {
55 first_local_index = 13911;
56 last_local_index = 26559;
57 }
58 }
59
60 // Read element dof indices from files
61 std::vector<std::vector<std::vector<PetscInt>>> elem_dof_indices(size);
62 for (PetscInt proc_id = 0; proc_id < size; ++proc_id) {
63 std::string line;
64 std::ifstream dof_file(file[proc_id]);
65 PetscCheck(dof_file.good(), PETSC_COMM_SELF, PETSC_ERR_FILE_OPEN, "Could not open file %s", file[proc_id]);
66 while (std::getline(dof_file, line)) {
67 std::vector<PetscInt> dof_indices;
68 std::stringstream sstream(line);
69 std::string token;
70 while (std::getline(sstream, token, ' ')) dof_indices.push_back(std::atoi(token.c_str()));
71 elem_dof_indices[proc_id].push_back(dof_indices);
72 }
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((unsigned int)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