Actual source code: ex231.cxx

  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:  */

 13: // PETSc includes
 14: #include <petscmat.h>

 16: // C++ includes
 17: #include <iostream>
 18: #include <fstream>
 19: #include <sstream>
 20: #include <algorithm>
 21: #include <vector>
 22: #include <string>
 23: #include <set>

 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;

 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");

 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:   }

 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:   }

 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:     if (dof_file.good()) {
 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:     } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_OPEN, "Could not open file %s", file[proc_id]);
 74:   }

 76:   // Debugging: Verify we read in elem_dof_indices correctly
 77:   // for (unsigned int i=0; i
 78:   //   {
 79:   //     for (unsigned int j=0; j
 80:   //       {
 81:   //         for (unsigned int k=0; k
 82:   //           std::cout << elem_dof_indices[i][j][k] << " ";
 83:   //         std::cout << std::endl;
 84:   //       }
 85:   //     std::cout << std::endl;
 86:   //   }

 88:   // Set up the (global) sparsity pattern
 89:   std::vector<std::set<unsigned int>> sparsity(n_dofs);
 90:   for (PetscInt proc_id = 0; proc_id < size; ++proc_id)
 91:     for (unsigned int k = 0; k < elem_dof_indices[proc_id].size(); k++) {
 92:       std::vector<PetscInt> &dof_indices = elem_dof_indices[proc_id][k];
 93:       for (unsigned int i = 0; i < dof_indices.size(); ++i)
 94:         for (unsigned int j = 0; j < dof_indices.size(); ++j) sparsity[dof_indices[i]].insert(dof_indices[j]);
 95:     }

 97:   // Determine the local nonzeros on this processor
 98:   const unsigned int    n_local_dofs = last_local_index - first_local_index;
 99:   std::vector<PetscInt> n_nz(n_local_dofs);
100:   std::vector<PetscInt> n_oz(n_local_dofs);

102:   for (unsigned int i = 0; i < n_local_dofs; ++i) {
103:     for (std::set<unsigned int>::iterator iter = sparsity[i + first_local_index].begin(); iter != sparsity[i + first_local_index].end(); iter++) {
104:       unsigned int dof = *iter;
105:       if ((dof >= first_local_index) && (dof < last_local_index)) n_nz[i]++;
106:       else n_oz[i]++;
107:     }
108:   }

110:   // Debugging: print number of on/off-diagonal nonzeros
111:   // for (unsigned int i=0; i
112:   //   std::cout << n_nz[i] << " ";
113:   // std::cout << std::endl;

115:   // for (unsigned int i=0; i
116:   //   std::cout << n_oz[i] << " ";
117:   // std::cout << std::endl;

119:   // Compute and print max number of on- and off-diagonal nonzeros.
120:   PetscInt n_nz_max = *std::max_element(n_nz.begin(), n_nz.end());
121:   PetscInt n_oz_max = *std::max_element(n_oz.begin(), n_oz.end());

123:   PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "Max on-diagonal non-zeros: = %" PetscInt_FMT "\n", n_nz_max));
124:   PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "Max off-diagonal non-zeros: = %" PetscInt_FMT "\n", n_oz_max));
125:   PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT));

127:   // Initialize the matrix similar to what we do in the PetscMatrix
128:   // ctor and init() routines.
129:   Mat mat;
130:   PetscCall(MatCreate(PETSC_COMM_WORLD, &mat));
131:   PetscCall(MatSetSizes(mat, n_local_dofs, n_local_dofs, n_dofs, n_dofs));
132:   PetscCall(MatSetBlockSize(mat, 1));
133:   PetscCall(MatSetType(mat, MATAIJ)); // Automatically chooses seqaij or mpiaij
134:   PetscCall(MatSeqAIJSetPreallocation(mat, 0, &n_nz[0]));
135:   PetscCall(MatMPIAIJSetPreallocation(mat, 0, &n_nz[0], 0, &n_oz[0]));
136:   PetscCall(MatSetOption(mat, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE));

138:   // Local "element" loop
139:   for (unsigned int k = 0; k < elem_dof_indices[rank].size(); k++) {
140:     std::vector<PetscInt> &dof_indices = elem_dof_indices[rank][k];
141:     // DenseMatrix< Number >  zero_mat( dof_indices.size(), dof_indices.size());
142:     // B.add_matrix( zero_mat, dof_indices);
143:     std::vector<PetscScalar> ones(dof_indices.size() * dof_indices.size(), 1.);
144:     PetscCall(MatSetValues(mat, dof_indices.size(), &dof_indices[0], dof_indices.size(), &dof_indices[0], &ones[0], ADD_VALUES));
145:   }

147:   // Matrix assembly
148:   PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY));
149:   PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY));

151:   // Clean up
152:   PetscCall(MatDestroy(&mat));
153:   PetscCall(PetscFinalize());
154:   return 0;
155: }
156: /*TEST
157:    build:
158:       requires: !defined(PETSC_HAVE_SUN_CXX)

160:    test:
161:       nsize: 2
162:       args: -f0 ${DATAFILESPATH}/meshes/moose_dof_indices_np_2_proc_0.txt -f1 ${DATAFILESPATH}/meshes/moose_dof_indices_np_2_proc_1.txt
163:       # Skip the test for Sun C++ compiler because of its warnings/errors in petscmat.h
164:       requires: datafilespath

166: TEST*/