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