GeneralBrokenLines  trunk_rev46
 All Classes Files Functions Variables Typedefs Pages
VMatrix.cpp
Go to the documentation of this file.
1 /*
2  * VMatrix.cpp
3  *
4  * Created on: Feb 15, 2012
5  * Author: kleinwrt
6  */
7 
8 #include "VMatrix.h"
9 
10 /*********** simple Matrix based on std::vector<double> **********/
11 
12 VMatrix::VMatrix(const unsigned int nRows, const unsigned int nCols) :
13  numRows(nRows), numCols(nCols), theVec(nRows * nCols) {
14 }
15 
16 VMatrix::VMatrix(const VMatrix &aMatrix) :
17  numRows(aMatrix.numRows), numCols(aMatrix.numCols), theVec(
18  aMatrix.theVec) {
19 
20 }
21 
22 VMatrix::~VMatrix() {
23 }
24 
26 
30 void VMatrix::resize(const unsigned int nRows, const unsigned int nCols) {
31  numRows = nRows;
32  numCols = nCols;
33  theVec.resize(nRows * nCols);
34 }
35 
37 
40 VMatrix VMatrix::transpose() const {
41  VMatrix aResult(numCols, numRows);
42  for (unsigned int i = 0; i < numRows; ++i) {
43  for (unsigned int j = 0; j < numCols; ++j) {
44  aResult(j, i) = theVec[numCols * i + j];
45  }
46  }
47  return aResult;
48 }
49 
51 
54 unsigned int VMatrix::getNumRows() const {
55  return numRows;
56 }
57 
59 
62 unsigned int VMatrix::getNumCols() const {
63  return numCols;
64 }
65 
67 void VMatrix::print() const {
68  std::cout << " VMatrix: " << numRows << "*" << numCols << std::endl;
69  for (unsigned int i = 0; i < numRows; ++i) {
70  for (unsigned int j = 0; j < numCols; ++j) {
71  if (j % 5 == 0) {
72  std::cout << std::endl << std::setw(4) << i << ","
73  << std::setw(4) << j << "-" << std::setw(4)
74  << std::min(j + 4, numCols) << ":";
75  }
76  std::cout << std::setw(13) << theVec[numCols * i + j];
77  }
78  }
79  std::cout << std::endl << std::endl;
80 }
81 
83 VVector VMatrix::operator*(const VVector &aVector) const {
84  VVector aResult(numRows);
85  for (unsigned int i = 0; i < this->numRows; ++i) {
86  double sum = 0.0;
87  for (unsigned int j = 0; j < this->numCols; ++j) {
88  sum += theVec[numCols * i + j] * aVector(j);
89  }
90  aResult(i) = sum;
91  }
92  return aResult;
93 }
94 
96 VMatrix VMatrix::operator*(const VMatrix &aMatrix) const {
97 
98  VMatrix aResult(numRows, aMatrix.numCols);
99  for (unsigned int i = 0; i < numRows; ++i) {
100  for (unsigned int j = 0; j < aMatrix.numCols; ++j) {
101  double sum = 0.0;
102  for (unsigned int k = 0; k < numCols; ++k) {
103  sum += theVec[numCols * i + k] * aMatrix(k, j);
104  }
105  aResult(i, j) = sum;
106  }
107  }
108  return aResult;
109 }
110 
112 VMatrix VMatrix::operator+(const VMatrix &aMatrix) const {
113  VMatrix aResult(numRows, numCols);
114  for (unsigned int i = 0; i < numRows; ++i) {
115  for (unsigned int j = 0; j < numCols; ++j) {
116  aResult(i, j) = theVec[numCols * i + j] + aMatrix(i, j);
117  }
118  }
119  return aResult;
120 }
121 
122 /*********** simple symmetric Matrix based on std::vector<double> **********/
123 
124 VSymMatrix::VSymMatrix(const unsigned int nRows) :
125  numRows(nRows), theVec((nRows * nRows + nRows) / 2) {
126 }
127 
128 VSymMatrix::~VSymMatrix() {
129 }
130 
132 
135 void VSymMatrix::resize(const unsigned int nRows) {
136  numRows = nRows;
137  theVec.resize((nRows * nRows + nRows) / 2);
138 }
139 
141 
144 unsigned int VSymMatrix::getNumRows() const {
145  return numRows;
146 }
147 
149 void VSymMatrix::print() const {
150  std::cout << " VSymMatrix: " << numRows << "*" << numRows << std::endl;
151  for (unsigned int i = 0; i < numRows; ++i) {
152  for (unsigned int j = 0; j <= i; ++j) {
153  if (j % 5 == 0) {
154  std::cout << std::endl << std::setw(4) << i << ","
155  << std::setw(4) << j << "-" << std::setw(4)
156  << std::min(j + 4, i) << ":";
157  }
158  std::cout << std::setw(13) << theVec[(i * i + i) / 2 + j];
159  }
160  }
161  std::cout << std::endl << std::endl;
162 }
163 
165 VSymMatrix VSymMatrix::operator-(const VMatrix &aMatrix) const {
166  VSymMatrix aResult(numRows);
167  for (unsigned int i = 0; i < numRows; ++i) {
168  for (unsigned int j = 0; j <= i; ++j) {
169  aResult(i, j) = theVec[(i * i + i) / 2 + j] - aMatrix(i, j);
170  }
171  }
172  return aResult;
173 }
174 
176 VVector VSymMatrix::operator*(const VVector &aVector) const {
177  VVector aResult(numRows);
178  for (unsigned int i = 0; i < numRows; ++i) {
179  aResult(i) = theVec[(i * i + i) / 2 + i] * aVector(i);
180  for (unsigned int j = 0; j < i; ++j) {
181  aResult(j) += theVec[(i * i + i) / 2 + j] * aVector(i);
182  aResult(i) += theVec[(i * i + i) / 2 + j] * aVector(j);
183  }
184  }
185  return aResult;
186 }
187 
189 VMatrix VSymMatrix::operator*(const VMatrix &aMatrix) const {
190  unsigned int nCol = aMatrix.getNumCols();
191  VMatrix aResult(numRows, nCol);
192  for (unsigned int l = 0; l < nCol; ++l) {
193  for (unsigned int i = 0; i < numRows; ++i) {
194  aResult(i, l) = theVec[(i * i + i) / 2 + i] * aMatrix(i, l);
195  for (unsigned int j = 0; j < i; ++j) {
196  aResult(j, l) += theVec[(i * i + i) / 2 + j] * aMatrix(i, l);
197  aResult(i, l) += theVec[(i * i + i) / 2 + j] * aMatrix(j, l);
198  }
199  }
200  }
201  return aResult;
202 }
203 
204 /*********** simple Vector based on std::vector<double> **********/
205 
206 VVector::VVector(const unsigned int nRows) :
207  numRows(nRows), theVec(nRows) {
208 }
209 
210 VVector::VVector(const VVector &aVector) :
211  numRows(aVector.numRows), theVec(aVector.theVec) {
212 
213 }
214 
215 VVector::~VVector() {
216 }
217 
219 
222 void VVector::resize(const unsigned int nRows) {
223  numRows = nRows;
224  theVec.resize(nRows);
225 }
226 
228 
233 VVector VVector::getVec(unsigned int len, unsigned int start) const {
234  VVector aResult(len);
235  std::memcpy(&aResult.theVec[0], &theVec[start], sizeof(double) * len);
236  return aResult;
237 }
238 
240 
244 void VVector::putVec(const VVector &aVector, unsigned int start) {
245  std::memcpy(&theVec[start], &aVector.theVec[0],
246  sizeof(double) * aVector.numRows);
247 }
248 
250 
253 unsigned int VVector::getNumRows() const {
254  return numRows;
255 }
256 
258 void VVector::print() const {
259  std::cout << " VVector: " << numRows << std::endl;
260  for (unsigned int i = 0; i < numRows; ++i) {
261 
262  if (i % 5 == 0) {
263  std::cout << std::endl << std::setw(4) << i << "-" << std::setw(4)
264  << std::min(i + 4, numRows) << ":";
265  }
266  std::cout << std::setw(13) << theVec[i];
267  }
268  std::cout << std::endl << std::endl;
269 }
270 
272 VVector VVector::operator-(const VVector &aVector) const {
273  VVector aResult(numRows);
274  for (unsigned int i = 0; i < numRows; ++i) {
275  aResult(i) = theVec[i] - aVector(i);
276  }
277  return aResult;
278 }
279 
280 /*============================================================================
281  from mpnum.F (MillePede-II by V. Blobel, Univ. Hamburg)
282  ============================================================================*/
284 
296 unsigned int VSymMatrix::invert() {
297 
298  const double eps = 1.0E-10;
299  unsigned int aSize = numRows;
300  std::vector<int> next(aSize);
301  std::vector<double> diag(aSize);
302  int nSize = aSize;
303 
304  int first = 1;
305  for (int i = 1; i <= nSize; ++i) {
306  next[i - 1] = i + 1; // set "next" pointer
307  diag[i - 1] = fabs(theVec[(i * i + i) / 2 - 1]); // save abs of diagonal elements
308  }
309  next[aSize - 1] = -1; // end flag
310 
311  unsigned int nrank = 0;
312  for (int i = 1; i <= nSize; ++i) { // start of loop
313  int k = 0;
314  double vkk = 0.0;
315 
316  int j = first;
317  int previous = 0;
318  int last = previous;
319  // look for pivot
320  while (j > 0) {
321  int jj = (j * j + j) / 2 - 1;
322  if (fabs(theVec[jj]) > std::max(fabs(vkk), eps * diag[j - 1])) {
323  vkk = theVec[jj];
324  k = j;
325  last = previous;
326  }
327  previous = j;
328  j = next[j - 1];
329  }
330  // pivot found
331  if (k > 0) {
332  int kk = (k * k + k) / 2 - 1;
333  if (last <= 0) {
334  first = next[k - 1];
335  } else {
336  next[last - 1] = next[k - 1];
337  }
338  next[k - 1] = 0; // index is used, reset
339  nrank++; // increase rank and ...
340 
341  vkk = 1.0 / vkk;
342  theVec[kk] = -vkk;
343  int jk = kk - k;
344  int jl = -1;
345  for (int j = 1; j <= nSize; ++j) { // elimination
346  if (j == k) {
347  jk = kk;
348  jl += j;
349  } else {
350  if (j < k) {
351  ++jk;
352  } else {
353  jk += j - 1;
354  }
355 
356  double vjk = theVec[jk];
357  theVec[jk] = vkk * vjk;
358  int lk = kk - k;
359  if (j >= k) {
360  for (int l = 1; l <= k - 1; ++l) {
361  ++jl;
362  ++lk;
363  theVec[jl] -= theVec[lk] * vjk;
364  }
365  ++jl;
366  lk = kk;
367  for (int l = k + 1; l <= j; ++l) {
368  ++jl;
369  lk += l - 1;
370  theVec[jl] -= theVec[lk] * vjk;
371  }
372  } else {
373  for (int l = 1; l <= j; ++l) {
374  ++jl;
375  ++lk;
376  theVec[jl] -= theVec[lk] * vjk;
377  }
378  }
379  }
380  }
381  } else {
382  for (int k = 1; k <= nSize; ++k) {
383  if (next[k - 1] >= 0) {
384  int kk = (k * k - k) / 2 - 1;
385  for (int j = 1; j <= nSize; ++j) {
386  if (next[j - 1] >= 0) {
387  theVec[kk + j] = 0.0; // clear matrix row/col
388  }
389  }
390  }
391  }
392  throw 1; // singular
393  }
394  }
395  for (int ij = 0; ij < (nSize * nSize + nSize) / 2; ++ij) {
396  theVec[ij] = -theVec[ij]; // finally reverse sign of all matrix elements
397  }
398  return nrank;
399 }
400 
Generated on Thu Sep 20 2012 15:45:38 for GeneralBrokenLines by   doxygen 1.8.1