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: {
 28:   PetscMPIInt    size, rank;
 29:   char           file[2][PETSC_MAX_PATH_LEN];
 30:   PetscBool      flg;
 31:   const unsigned int n_dofs = 26559;
 32:   unsigned int   first_local_index;
 33:   unsigned int   last_local_index;

 35:   PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
 36:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 37:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
 38:   if (size > 2) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"This example is for <=2 procs");

 40:   PetscOptionsGetString(NULL,NULL,"-f0",file[0],sizeof(file[0]),&flg);
 41:   if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"Must indicate dof indices file for rank 0 with -f0 option");
 42:   if (size == 2) {
 43:     PetscOptionsGetString(NULL,NULL,"-f1",file[1],sizeof(file[1]),&flg);
 44:     if (!flg) SETERRQ(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 SETERRQ1(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<elem_dof_indices.size(); ++i)
 78:   //   {
 79:   //     for (unsigned int j=0; j<elem_dof_indices[i].size(); ++j)
 80:   //       {
 81:   //         for (unsigned int k=0; k<elem_dof_indices[i][j].size(); ++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)
 95:           sparsity[dof_indices[i]].insert(dof_indices[j]);
 96:     }

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

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

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

116:   // for (unsigned int i=0; i<n_oz.size(); ++i)
117:   //   std::cout << n_oz[i] << " ";
118:   // std::cout << std::endl;

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

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

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

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

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

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

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

167: TEST*/