Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
VMatrix.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file VMatrix.cc
1 /*
2  * VMatrix.cpp
3  *
4  * Created on: Feb 15, 2012
5  * Author: kleinwrt
6  */
7 
30 #include "VMatrix.h"
31 
33 namespace gbl {
34 
35 /*********** simple Matrix based on std::vector<double> **********/
36 
37 VMatrix::VMatrix(const unsigned int nRows, const unsigned int nCols) :
38  numRows(nRows), numCols(nCols), theVec(nRows * nCols) {
39 }
40 
41 VMatrix::VMatrix(const VMatrix &aMatrix) :
42  numRows(aMatrix.numRows), numCols(aMatrix.numCols), theVec(
43  aMatrix.theVec) {
44 
45 }
46 
47 VMatrix::~VMatrix() {
48 }
49 
51 
55 void VMatrix::resize(const unsigned int nRows, const unsigned int nCols) {
56  numRows = nRows;
57  numCols = nCols;
58  theVec.resize(nRows * nCols);
59 }
60 
62 
65 VMatrix VMatrix::transpose() const {
66  VMatrix aResult(numCols, numRows);
67  for (unsigned int i = 0; i < numRows; ++i) {
68  for (unsigned int j = 0; j < numCols; ++j) {
69  aResult(j, i) = theVec[numCols * i + j];
70  }
71  }
72  return aResult;
73 }
74 
76 
79 unsigned int VMatrix::getNumRows() const {
80  return numRows;
81 }
82 
84 
87 unsigned int VMatrix::getNumCols() const {
88  return numCols;
89 }
90 
92 void VMatrix::print() const {
93  std::cout << " VMatrix: " << numRows << "*" << numCols << std::endl;
94  for (unsigned int i = 0; i < numRows; ++i) {
95  for (unsigned int j = 0; j < numCols; ++j) {
96  if (j % 5 == 0) {
97  std::cout << std::endl << std::setw(4) << i << ","
98  << std::setw(4) << j << "-" << std::setw(4)
99  << std::min(j + 4, numCols) << ":";
100  }
101  std::cout << std::setw(13) << theVec[numCols * i + j];
102  }
103  }
104  std::cout << std::endl << std::endl;
105 }
106 
108 VVector VMatrix::operator*(const VVector &aVector) const {
109  VVector aResult(numRows);
110  for (unsigned int i = 0; i < this->numRows; ++i) {
111  double sum = 0.0;
112  for (unsigned int j = 0; j < this->numCols; ++j) {
113  sum += theVec[numCols * i + j] * aVector(j);
114  }
115  aResult(i) = sum;
116  }
117  return aResult;
118 }
119 
121 VMatrix VMatrix::operator*(const VMatrix &aMatrix) const {
122 
123  VMatrix aResult(numRows, aMatrix.numCols);
124  for (unsigned int i = 0; i < numRows; ++i) {
125  for (unsigned int j = 0; j < aMatrix.numCols; ++j) {
126  double sum = 0.0;
127  for (unsigned int k = 0; k < numCols; ++k) {
128  sum += theVec[numCols * i + k] * aMatrix(k, j);
129  }
130  aResult(i, j) = sum;
131  }
132  }
133  return aResult;
134 }
135 
137 VMatrix VMatrix::operator+(const VMatrix &aMatrix) const {
138  VMatrix aResult(numRows, numCols);
139  for (unsigned int i = 0; i < numRows; ++i) {
140  for (unsigned int j = 0; j < numCols; ++j) {
141  aResult(i, j) = theVec[numCols * i + j] + aMatrix(i, j);
142  }
143  }
144  return aResult;
145 }
146 
148 VMatrix &VMatrix::operator=(const VMatrix &aMatrix) {
149  if (this != &aMatrix) { // Gracefully handle self assignment
150  numRows = aMatrix.getNumRows();
151  numCols = aMatrix.getNumCols();
152  theVec.resize(numRows * numCols);
153  for (unsigned int i = 0; i < numRows; ++i) {
154  for (unsigned int j = 0; j < numCols; ++j) {
155  theVec[numCols * i + j] = aMatrix(i, j);
156  }
157  }
158  }
159  return *this;
160 }
161 
162 /*********** simple symmetric Matrix based on std::vector<double> **********/
163 
164 VSymMatrix::VSymMatrix(const unsigned int nRows) :
165  numRows(nRows), theVec((nRows * nRows + nRows) / 2) {
166 }
167 
168 VSymMatrix::~VSymMatrix() {
169 }
170 
172 
175 void VSymMatrix::resize(const unsigned int nRows) {
176  numRows = nRows;
177  theVec.resize((nRows * nRows + nRows) / 2);
178 }
179 
181 
184 unsigned int VSymMatrix::getNumRows() const {
185  return numRows;
186 }
187 
189 void VSymMatrix::print() const {
190  std::cout << " VSymMatrix: " << numRows << "*" << numRows << std::endl;
191  for (unsigned int i = 0; i < numRows; ++i) {
192  for (unsigned int j = 0; j <= i; ++j) {
193  if (j % 5 == 0) {
194  std::cout << std::endl << std::setw(4) << i << ","
195  << std::setw(4) << j << "-" << std::setw(4)
196  << std::min(j + 4, i) << ":";
197  }
198  std::cout << std::setw(13) << theVec[(i * i + i) / 2 + j];
199  }
200  }
201  std::cout << std::endl << std::endl;
202 }
203 
205 VSymMatrix VSymMatrix::operator-(const VMatrix &aMatrix) const {
206  VSymMatrix aResult(numRows);
207  for (unsigned int i = 0; i < numRows; ++i) {
208  for (unsigned int j = 0; j <= i; ++j) {
209  aResult(i, j) = theVec[(i * i + i) / 2 + j] - aMatrix(i, j);
210  }
211  }
212  return aResult;
213 }
214 
216 VVector VSymMatrix::operator*(const VVector &aVector) const {
217  VVector aResult(numRows);
218  for (unsigned int i = 0; i < numRows; ++i) {
219  aResult(i) = theVec[(i * i + i) / 2 + i] * aVector(i);
220  for (unsigned int j = 0; j < i; ++j) {
221  aResult(j) += theVec[(i * i + i) / 2 + j] * aVector(i);
222  aResult(i) += theVec[(i * i + i) / 2 + j] * aVector(j);
223  }
224  }
225  return aResult;
226 }
227 
229 VMatrix VSymMatrix::operator*(const VMatrix &aMatrix) const {
230  unsigned int nCol = aMatrix.getNumCols();
231  VMatrix aResult(numRows, nCol);
232  for (unsigned int l = 0; l < nCol; ++l) {
233  for (unsigned int i = 0; i < numRows; ++i) {
234  aResult(i, l) = theVec[(i * i + i) / 2 + i] * aMatrix(i, l);
235  for (unsigned int j = 0; j < i; ++j) {
236  aResult(j, l) += theVec[(i * i + i) / 2 + j] * aMatrix(i, l);
237  aResult(i, l) += theVec[(i * i + i) / 2 + j] * aMatrix(j, l);
238  }
239  }
240  }
241  return aResult;
242 }
243 
244 /*********** simple Vector based on std::vector<double> **********/
245 
246 VVector::VVector(const unsigned int nRows) :
247  numRows(nRows), theVec(nRows) {
248 }
249 
250 VVector::VVector(const VVector &aVector) :
251  numRows(aVector.numRows), theVec(aVector.theVec) {
252 
253 }
254 
255 VVector::~VVector() {
256 }
257 
259 
262 void VVector::resize(const unsigned int nRows) {
263  numRows = nRows;
264  theVec.resize(nRows);
265 }
266 
268 
273 VVector VVector::getVec(unsigned int len, unsigned int start) const {
274  VVector aResult(len);
275  std::memcpy(&aResult.theVec[0], &theVec[start], sizeof(double) * len);
276  return aResult;
277 }
278 
280 
284 void VVector::putVec(const VVector &aVector, unsigned int start) {
285  std::memcpy(&theVec[start], &aVector.theVec[0],
286  sizeof(double) * aVector.numRows);
287 }
288 
290 
293 unsigned int VVector::getNumRows() const {
294  return numRows;
295 }
296 
298 void VVector::print() const {
299  std::cout << " VVector: " << numRows << std::endl;
300  for (unsigned int i = 0; i < numRows; ++i) {
301 
302  if (i % 5 == 0) {
303  std::cout << std::endl << std::setw(4) << i << "-" << std::setw(4)
304  << std::min(i + 4, numRows) << ":";
305  }
306  std::cout << std::setw(13) << theVec[i];
307  }
308  std::cout << std::endl << std::endl;
309 }
310 
312 VVector VVector::operator-(const VVector &aVector) const {
313  VVector aResult(numRows);
314  for (unsigned int i = 0; i < numRows; ++i) {
315  aResult(i) = theVec[i] - aVector(i);
316  }
317  return aResult;
318 }
319 
321 VVector &VVector::operator=(const VVector &aVector) {
322  if (this != &aVector) { // Gracefully handle self assignment
323  numRows = aVector.getNumRows();
324  theVec.resize(numRows);
325  for (unsigned int i = 0; i < numRows; ++i) {
326  theVec[i] = aVector(i);
327  }
328  }
329  return *this;
330 }
331 
332 /*============================================================================
333  from mpnum.F (MillePede-II by V. Blobel, Univ. Hamburg)
334  ============================================================================*/
336 
348 unsigned int VSymMatrix::invert() {
349 
350  const double eps = 1.0E-10;
351  unsigned int aSize = numRows;
352  std::vector<int> next(aSize);
353  std::vector<double> diag(aSize);
354  int nSize = aSize;
355 
356  int first = 1;
357  for (int i = 1; i <= nSize; ++i) {
358  next[i - 1] = i + 1; // set "next" pointer
359  diag[i - 1] = fabs(theVec[(i * i + i) / 2 - 1]); // save abs of diagonal elements
360  }
361  next[aSize - 1] = -1; // end flag
362 
363  unsigned int nrank = 0;
364  for (int i = 1; i <= nSize; ++i) { // start of loop
365  int k = 0;
366  double vkk = 0.0;
367 
368  int j = first;
369  int previous = 0;
370  int last = previous;
371  // look for pivot
372  while (j > 0) {
373  int jj = (j * j + j) / 2 - 1;
374  if (fabs(theVec[jj]) > std::max(fabs(vkk), eps * diag[j - 1])) {
375  vkk = theVec[jj];
376  k = j;
377  last = previous;
378  }
379  previous = j;
380  j = next[j - 1];
381  }
382  // pivot found
383  if (k > 0) {
384  int kk = (k * k + k) / 2 - 1;
385  if (last <= 0) {
386  first = next[k - 1];
387  } else {
388  next[last - 1] = next[k - 1];
389  }
390  next[k - 1] = 0; // index is used, reset
391  nrank++; // increase rank and ...
392 
393  vkk = 1.0 / vkk;
394  theVec[kk] = -vkk;
395  int jk = kk - k;
396  int jl = -1;
397  for (int m = 1; m <= nSize; ++m) { // elimination
398  if (m == k) {
399  jk = kk;
400  jl += m;
401  } else {
402  if (m < k) {
403  ++jk;
404  } else {
405  jk += m - 1;
406  }
407 
408  double vjk = theVec[jk];
409  theVec[jk] = vkk * vjk;
410  int lk = kk - k;
411  if (m >= k) {
412  for (int l = 1; l <= k - 1; ++l) {
413  ++jl;
414  ++lk;
415  theVec[jl] -= theVec[lk] * vjk;
416  }
417  ++jl;
418  lk = kk;
419  for (int l = k + 1; l <= m; ++l) {
420  ++jl;
421  lk += l - 1;
422  theVec[jl] -= theVec[lk] * vjk;
423  }
424  } else {
425  for (int l = 1; l <= m; ++l) {
426  ++jl;
427  ++lk;
428  theVec[jl] -= theVec[lk] * vjk;
429  }
430  }
431  }
432  }
433  } else {
434  for (int m = 1; m <= nSize; ++m) {
435  if (next[m - 1] >= 0) {
436  int kk = (m * m - m) / 2 - 1;
437  for (int n = 1; n <= nSize; ++n) {
438  if (next[n - 1] >= 0) {
439  theVec[kk + n] = 0.0; // clear matrix row/col
440  }
441  }
442  }
443  }
444  throw 1; // singular
445  }
446  }
447  for (int ij = 0; ij < (nSize * nSize + nSize) / 2; ++ij) {
448  theVec[ij] = -theVec[ij]; // finally reverse sign of all matrix elements
449  }
450  return nrank;
451 }
452 }