CMS 3D CMS Logo

BlockSolver.cc
Go to the documentation of this file.
1 
6 int
7 BlockSolver::operator () (const CLHEP::HepMatrix & matrix,
8  const CLHEP::HepVector & vector,
9  CLHEP::HepVector & result)
10 {
11  double threshold = matrix.trace () /
12  static_cast<double> (matrix.num_row ()) ;
13  std::vector<int> holes ;
14  // loop over the matrix rows
15  for (int row = 0 ; row < matrix.num_row () ; ++row)
16  {
17  double sumAbs = 0. ;
18  for (int col = 0 ; col < matrix.num_col () ; ++col)
19  sumAbs += fabs (matrix[row][col]) ;
20  if (sumAbs < threshold) holes.push_back (row) ;
21  } // loop over the matrix rows
22 
23  int dim = matrix.num_col () - holes.size () ;
24 
25  if (holes.empty()) //PG exceptional case!
26  {
27  for (int i = 0 ; i < result.num_row () ; ++i)
28  result[i] = 1. ;
29  }
30  else if (dim > 0)
31  {
32  CLHEP::HepMatrix solution (dim, dim, 0) ;
33  CLHEP::HepVector input (dim, 0) ;
34  shrink (matrix, solution, vector, input, holes) ;
35  CLHEP::HepVector output = solve (solution,input) ;
36  pour (result, output, holes) ;
37  }
38  return holes.size () ;
39 }
40 
41 
42 // ------------------------------------------------------------
43 
44 
45 void
46 BlockSolver::shrink (const CLHEP::HepMatrix & matrix,
47  CLHEP::HepMatrix & solution,
48  const CLHEP::HepVector & result,
49  CLHEP::HepVector & input,
50  const std::vector<int> & where)
51 {
52 
53  int offsetRow = 0 ;
54  std::vector<int>::const_iterator whereRows = where.begin () ;
55  // loop over rows
56  for (int row = 0 ; row < matrix.num_row () ; ++row)
57  {
58  if (row == *whereRows)
59  {
60 // std::cerr << " DEBUG shr hole found " << std::endl ;
61  ++offsetRow ;
62  ++whereRows ;
63  continue ;
64  }
65  input[row-offsetRow] = result[row] ;
66  int offsetCol = 0 ;
67  std::vector<int>::const_iterator whereCols = where.begin () ;
68  // loop over columns
69  for (int col = 0 ; col < matrix.num_col () ; ++col)
70  {
71  if (col == *whereCols)
72  {
73  ++offsetCol ;
74  ++whereCols ;
75  continue ;
76  }
77  solution[row-offsetRow][col-offsetCol] = matrix[row][col] ;
78  }
79  } // loop over rows
80  return ;
81 }
82 
83 
84 // ------------------------------------------------------------
85 
86 
87 void
88 BlockSolver::pour (CLHEP::HepVector & result,
89  const CLHEP::HepVector & output,
90  const std::vector<int> & where)
91 {
92  std::vector<int>::const_iterator whereCols = where.begin () ;
93  int readingIndex = 0 ;
94  //PG loop over the output crystals
95  for (int xtal = 0 ; xtal < result.num_row () ; ++xtal)
96  {
97  if (xtal == *whereCols)
98  {
99  result[xtal] = 1. ;
100  ++whereCols ;
101  }
102  else
103  {
104  result[xtal] = output[readingIndex] ;
105  ++readingIndex ;
106  }
107  } //PG loop over the output crystals
108 
109  return ;
110 }
111 
void shrink(const CLHEP::HepMatrix &matrix, CLHEP::HepMatrix &solution, const CLHEP::HepVector &result, CLHEP::HepVector &input, const std::vector< int > &where)
eliminate empty columns and rows
Definition: BlockSolver.cc:46
static std::string const input
Definition: EdmProvDump.cc:48
return((rh^lh)&mask)
int operator()(const CLHEP::HepMatrix &matrix, const CLHEP::HepVector &vector, CLHEP::HepVector &result)
Definition: BlockSolver.cc:7
col
Definition: cuy.py:1010
void pour(CLHEP::HepVector &result, const CLHEP::HepVector &output, const std::vector< int > &where)
pour results in bigger vector
Definition: BlockSolver.cc:88