Free Electron
SparseMatrix3x3.h
Go to the documentation of this file.
1 /* Copyright (C) 2003-2021 Free Electron Organization
2  Any use of this software requires a license. If a valid license
3  was not distributed with this file, visit freeelectron.org. */
4 
5 /** @file */
6 
7 #ifndef __SparseMatrix3x3_h__
8 #define __SparseMatrix3x3_h__
9 
10 namespace fe
11 {
12 namespace ext
13 {
14 
15 template <typename T>
16 class FE_DL_EXPORT SparseMatrix3x3 : public Component,
17  public CastableAs< SparseMatrix3x3<T> >
18 {
19  public:
20 typedef T t_real;
21 typedef Vector<3, t_real> t_v3;
22 typedef Matrix<3,3,t_real> t_matrix;
23 
24 virtual ~SparseMatrix3x3(void) {}
25 virtual void multiply(std::vector<t_v3> &a_b,
26  const std::vector<t_v3> &a_x) const = 0;
27 virtual void scale(t_real a_scale) = 0;
28 };
29 
30 /// Vector Map Rows (vector of row maps)
31 class FE_DL_EXPORT FullVMR : public SparseMatrix3x3<Real>
32 {
33  public:
34  FullVMR(void);
35 virtual ~FullVMR(void);
36 
37 virtual void multiply(std::vector<SpatialVector> &a_b,
38  const std::vector<SpatialVector> &a_x) const;
39 virtual void scale(Real a_scale);
40 
41  typedef std::map<unsigned int, SpatialMatrix> t_row;
42 
43  void clear(void);
44  void setRows(unsigned int a_count);
45  unsigned int rows(void);
46  t_row &row(unsigned int a_index);
47  void zero(void);
48 
49  SpatialMatrix &block(unsigned int a_i, unsigned int a_j);
50 
51  protected:
52  std::vector<t_row> m_rows;
53 
54 };
55 
56 /// Upper Triangular Vector Map Rows (vector of row maps)
57 class FE_DL_EXPORT UpperTriangularVMR : public FullVMR
58 {
59  public:
60  UpperTriangularVMR(void) {}
61 virtual ~UpperTriangularVMR(void) {}
62 
63 virtual void multiply(std::vector<SpatialVector> &a_b,
64  const std::vector<SpatialVector> &a_x) const;
65  private:
66 
67 };
68 
69 template <typename T>
70 class FE_DL_EXPORT CRS
71 {
72  public:
73  CRS(void){}
74 virtual ~CRS(void){}
75 
76  std::vector<T> m_values;
77  std::vector<int> m_colind;
78  std::vector<int> m_rowptr;
79  unsigned int m_r;
80 };
81 
82 template <typename T>
83 class FE_DL_EXPORT FullBCRS : public SparseMatrix3x3<T>
84 {
85  public:
86 typedef T t_real;
87 typedef Vector<3, t_real> t_v3;
88 typedef Matrix<3,3,t_real> t_matrix;
89  FullBCRS(void);
90 virtual ~FullBCRS(void);
91 
92 
93 virtual void multiply(std::vector<t_v3> &a_b,
94  const std::vector<t_v3> &a_x) const;
95 virtual void scale(t_real a_scale);
96 
97  void clear(void);
98  void startrow(void);
99  void done(void);
100  t_matrix &insert(unsigned int a_columnIndex,
101  const t_matrix &a_block);
102 
103  void dump(void) const;
104 
105  void write(const String &a_filename);
106  void read(const String &a_filename);
107 
108  void writeDense(const String &a_filename);
109  void readDense(const String &a_filename, unsigned int a_M, unsigned int a_N);
110 
111  void writeStructure(const String &a_filename) const;
112 
113  void sparse(const sp<FullBCRS> &a_other);
114  void project(const sp<FullBCRS> &a_other);
115 
116  void write(CRS<T> &aCRS) const;
117  void read(const CRS<T> &aCRS);
118 
119 
120  public:
121  std::vector<t_matrix> &blocks(void) { return m_blocks; }
122  std::vector<unsigned int> &colind(void) { return m_colind; }
123  std::vector<unsigned int> &rowptr(void) { return m_rowptr; }
124 
125  protected:
126  std::vector<t_matrix> m_blocks;
127  std::vector<unsigned int> m_colind;
128  std::vector<unsigned int> m_rowptr;
129  unsigned int m_r;
130  t_v3 m_zeroVector;
131 };
132 
133 template <typename T>
134 class FE_DL_EXPORT UpperTriangularBCRS : public FullBCRS<T>
135 {
136  public:
137 typedef T t_real;
138 typedef Vector<3, t_real> t_v3;
139 typedef Matrix<3,3,t_real> t_matrix;
140  UpperTriangularBCRS(void){}
141 virtual ~UpperTriangularBCRS(void){}
142 
143 virtual void multiply(std::vector<t_v3> &a_b,
144  const std::vector<t_v3> &a_x) const;
145 
146  void incompleteSqrt(void);
147  void backSub(std::vector<t_v3> &a_x, const std::vector<t_v3> &a_b);
148  void transposeForeSub(std::vector<t_v3> &a_x, const std::vector<t_v3> &a_b);
149  void backSolve(std::vector<t_v3> &a_x, const std::vector<t_v3> &a_b);
150 
151 };
152 
153 
154 
155 template <typename T>
156 inline FullBCRS<T>::FullBCRS(void)
157 {
158  m_r = 0;
159  m_zeroVector = t_v3(0.0,0.0,0.0);
160 }
161 
162 template <typename T>
163 inline FullBCRS<T>::~FullBCRS(void)
164 {
165 }
166 
167 template <typename T>
168 inline void FullBCRS<T>::clear(void)
169 {
170  m_blocks.clear();
171  m_colind.clear();
172  m_rowptr.clear();
173  m_r = 0;
174 }
175 
176 template <typename T>
177 inline void FullBCRS<T>::startrow(void)
178 {
179  m_rowptr.push_back(m_r);
180 }
181 
182 template <typename T>
183 inline void FullBCRS<T>::done(void)
184 {
185  m_rowptr.push_back(m_r);
186 }
187 
188 template <typename T>
189 inline Matrix<3,3,T> &FullBCRS<T>::insert(unsigned int a_columnIndex, const t_matrix &a_block)
190 {
191  m_blocks.push_back(a_block);
192  m_colind.push_back(a_columnIndex);
193  m_r++;
194  return m_blocks.back();
195 }
196 
197 template <typename T>
198 inline void FullBCRS<T>::multiply(std::vector<t_v3> &a_b, const std::vector<t_v3> &a_x) const
199 {
200  unsigned int n = m_rowptr.size() - 1;
201  for(unsigned int i = 0; i < n; i++)
202  {
203  a_b[i] = m_zeroVector;
204  for(unsigned int k = m_rowptr[i]; k < m_rowptr[i+1]; k++)
205  {
206  a_b[i] += fe::multiply(m_blocks[k], a_x[m_colind[k]]);
207  }
208  }
209 }
210 
211 template <typename T>
212 inline void FullBCRS<T>::dump(void) const
213 {
214  t_matrix offdiagonal_sum;
215  setAll(offdiagonal_sum, 0.0);
216  unsigned int n = m_rowptr.size() - 1;
217  for(unsigned int i = 0; i < n; i++)
218  {
219  for(unsigned int k = m_rowptr[i]; k < m_rowptr[i+1]; k++)
220  {
221  if(m_colind[k] != i) // diagonal
222  {
223  offdiagonal_sum = offdiagonal_sum + m_blocks[k];
224  }
225  }
226  }
227 
228  String s;
229  s = print(offdiagonal_sum);
230  feLogGroup("BCRS", "off diagonal\n%s\n", s.c_str());
231 }
232 
233 template <typename T>
234 inline void FullBCRS<T>::writeStructure(const String &a_filename) const
235 {
236  FILE *fp = fopen(a_filename.c_str(), "w");
237  if(!fp)
238  {
239  feLog("WARNING: FullBCRS::writeStructure failed\n");
240  return;
241  }
242  t_matrix zero;
243  setAll(zero, 0.0);
244  unsigned int n = m_rowptr.size() - 1;
245  for(unsigned int i = 0; i < n; i++)
246  {
247  unsigned int j = 0;
248  for(unsigned int k = m_rowptr[i]; k < m_rowptr[i+1]; k++)
249  {
250  while(m_colind[k] > j)
251  {
252  fe_fprintf(fp, "-");
253  j++;
254  }
255  if(m_blocks[k] == zero)
256  {
257  fe_fprintf(fp, "0");
258  }
259  else
260  {
261  fe_fprintf(fp, "*");
262  }
263  j++;
264  }
265  while(j<n)
266  {
267  fe_fprintf(fp, "-");
268  j++;
269  }
270  fe_fprintf(fp, "\n");
271  }
272  fclose(fp);
273 }
274 
275 template <typename T>
276 inline void FullBCRS<T>::write(const String &a_filename)
277 {
278  FILE *fp = fopen(a_filename.c_str(), "w");
279  if(!fp)
280  {
281  feLog("WARNING: FullBCRS::write failed\n");
282  return;
283  }
284  fe_fprintf(fp, "%d\n", (int)m_blocks.size());
285  for(unsigned int i = 0; i < m_blocks.size(); i++)
286  {
287  fe_fprintf(fp, "%g %g %g\n", m_blocks[i](0,0), m_blocks[i](0,1), m_blocks[i](0,2));
288  fe_fprintf(fp, "%g %g %g\n", m_blocks[i](1,0), m_blocks[i](1,1), m_blocks[i](1,2));
289  fe_fprintf(fp, "%g %g %g\n", m_blocks[i](2,0), m_blocks[i](2,1), m_blocks[i](2,2));
290  }
291  fe_fprintf(fp, "%d\n", (int)m_rowptr.size());
292  for(unsigned int i = 0; i < m_rowptr.size(); i++)
293  {
294  fe_fprintf(fp, "%d ", m_rowptr[i]);
295  }
296  fe_fprintf(fp, "\n");
297 
298  fe_fprintf(fp, "%d\n", (int)m_colind.size());
299  for(unsigned int i = 0; i < m_colind.size(); i++)
300  {
301  fe_fprintf(fp, "%d ", m_colind[i]);
302  }
303  fe_fprintf(fp, "\n");
304 
305  fe_fprintf(fp, "%d\n", m_r);
306  fclose(fp);
307 }
308 
309 template <typename T>
310 inline void FullBCRS<T>::read(const String &a_filename)
311 {
312  std::ifstream infile(a_filename.c_str());
313  int m;
314  infile >> m;
315  m_blocks.resize(m);
316  for(int i = 0; i < m; i++)
317  {
318  infile >> m_blocks[i](0,0);
319  infile >> m_blocks[i](0,1);
320  infile >> m_blocks[i](0,2);
321  infile >> m_blocks[i](1,0);
322  infile >> m_blocks[i](1,1);
323  infile >> m_blocks[i](1,2);
324  infile >> m_blocks[i](2,0);
325  infile >> m_blocks[i](2,1);
326  infile >> m_blocks[i](2,2);
327  }
328 
329  infile >> m;
330  m_rowptr.resize(m);
331  for(int i = 0; i < m; i++)
332  {
333  infile >> m_rowptr[i];
334  }
335 
336  infile >> m;
337  m_colind.resize(m);
338  for(int i = 0; i < m; i++)
339  {
340  infile >> m_colind[i];
341  }
342 
343  infile >> m_r;
344 
345  infile.close();
346 }
347 
348 template <typename T>
349 inline void FullBCRS<T>::writeDense(const String &a_filename)
350 {
351  FILE *fp = fopen(a_filename.c_str(), "w");
352  if(!fp)
353  {
354  feLog("WARNING: FullBCRS::write failed\n");
355  return;
356  }
357  unsigned int n = m_rowptr.size() - 1;
358 
359  std::vector< std::vector<double> > dense_matrix;
360  dense_matrix.resize(n*3);
361  for(unsigned int i = 0; i < n*3; i++)
362  {
363  dense_matrix[i].resize(n*3);
364  }
365 
366 #define FMT " % 10.8e % 10.8e % 10.8e"
367 #if 0
368  for(unsigned int i = 0; i < n; i++)
369  {
370  for(unsigned int k = m_rowptr[i]; k < m_rowptr[i+1]; k++)
371  {
372  unsigned int j = m_colind[k];
373 
374  for(unsigned int ii = 0; ii < 3; ii++)
375  {
376  for(unsigned int jj = 0; jj < 3; jj++)
377  {
378  unsigned int m = i*3 + ii;
379  unsigned int n = j*3 + jj;
380  dense_matrix[m][n] = m_blocks[k](ii,jj);
381  dense_matrix[n][m] = dense_matrix[m][n];
382  }
383  }
384  }
385  }
386 
387  for(unsigned int i = 0; i < n*3; i++)
388  {
389  for(unsigned int j = 0; j < n*3; j++)
390  {
391  fe_fprintf(fp, "%+10.8e ", dense_matrix[i][j]);
392  }
393  fe_fprintf(fp, "\n");
394  }
395 #endif
396 
397 #if 1
398  for(unsigned int i = 0; i < n; i++)
399  {
400  for(unsigned int r = 0; r < 3; r++)
401  {
402  unsigned int m = 0;
403  for(unsigned int k = m_rowptr[i]; k < m_rowptr[i+1]; k++)
404  {
405  unsigned int c = m_colind[k];
406  while(m < c)
407  {
408  m++;
409  fe_fprintf(fp, FMT, 0.0, 0.0, 0.0);
410  }
411  m++;
412  fe_fprintf(fp, FMT, m_blocks[k](r,0), m_blocks[k](r,1), m_blocks[k](r,2));
413  }
414  while(m<n)
415  {
416  m++;
417  fe_fprintf(fp, FMT, 0.0, 0.0, 0.0);
418  }
419  fe_fprintf(fp, "\n");
420  }
421  }
422 #endif
423 
424  fclose(fp);
425 }
426 
427 template <typename T>
428 inline void FullBCRS<T>::readDense(const String &a_filename, unsigned int a_M, unsigned int a_N)
429 {
430  m_blocks.clear();
431  m_colind.clear();
432  m_rowptr.clear();
433  for(unsigned int i = 0; i < a_M; i++)
434  {
435  m_rowptr.push_back(m_colind.size());
436  for(unsigned int j = 0; j < a_N; j++)
437  {
438  if(j>=i)
439  {
440  m_colind.push_back(j);
441  }
442  }
443  }
444  m_rowptr.push_back(m_colind.size());
445  m_r = m_colind.size();
446  m_blocks.resize(m_r);
447 
448  std::ifstream infile(a_filename.c_str());
449  unsigned int k = 0;
450  for(unsigned int i = 0; i < a_M; i++)
451  {
452  unsigned int kk = k;
453  for(unsigned int ii = 0; ii < 3; ii++)
454  {
455  k = kk;
456  for(unsigned int j = 0; j < a_N; j++)
457  {
458  for(unsigned int jj = 0; jj < 3; jj++)
459  {
460  if(j>=i)
461  {
462  infile >> m_blocks[k](ii,jj);
463  }
464  else
465  {
466  double dbl;
467  infile >> dbl;
468  }
469  }
470  if(j>=i){k++;}
471  }
472  }
473  }
474 
475  infile.close();
476 }
477 
478 template <typename T>
479 inline void FullBCRS<T>::write(CRS<T> &aCRS) const
480 {
481  unsigned int n = m_rowptr.size() - 1;
482 
483  aCRS.m_values.resize(m_blocks.size() * 9);
484  aCRS.m_colind.resize(m_blocks.size() * 9);
485  aCRS.m_rowptr.resize(n*3 + 1);
486 
487  unsigned int ii = 0;
488  for(unsigned int i = 0; i < n; i++) // block rows
489  {
490  aCRS.m_rowptr[i] = ii;
491  for(unsigned int r = 0; r < 3; r++) // unit rows
492  {
493  unsigned int m = i*3 + r;
494  for(unsigned int k = m_rowptr[i]; k < m_rowptr[i+1]; k++) // blocks
495  {
496  unsigned int c = m_colind[k];
497  aCRS.m_colind[ii + 0] = c*3 + 0;
498  aCRS.m_colind[ii + 1] = c*3 + 1;
499  aCRS.m_colind[ii + 2] = c*3 + 2;
500 
501  aCRS.m_values[ii + 0] = m_blocks[k](r, 0);
502  aCRS.m_values[ii + 1] = m_blocks[k](r, 1);
503  aCRS.m_values[ii + 2] = m_blocks[k](r, 2);
504 
505  ii += 3;
506  }
507  }
508  }
509  aCRS.m_rowptr[n] = ii;
510  assert(ii == aCRS.m_values.size());
511 
512 
513 #if 0
514  for(unsigned int i = 0; i < n; i++)
515  {
516  for(unsigned int r = 0; r < 3; r++)
517  {
518  unsigned int m = 0;
519  for(unsigned int k = m_rowptr[i]; k < m_rowptr[i+1]; k++)
520  {
521  unsigned int c = m_colind[k];
522  while(m < c)
523  {
524  m++;
525  fe_fprintf(fp, FMT, 0.0, 0.0, 0.0);
526  }
527  m++;
528  fe_fprintf(fp, FMT, m_blocks[k](r,0), m_blocks[k](r,1), m_blocks[k](r,2));
529  }
530  while(m<n)
531  {
532  m++;
533  fe_fprintf(fp, FMT, 0.0, 0.0, 0.0);
534  }
535  fe_fprintf(fp, "\n");
536  }
537  }
538 #endif
539 
540 }
541 
542 template <typename T>
543 inline void FullBCRS<T>::read(const CRS<T> &aCRS)
544 {
545 }
546 
547 
548 template <typename T>
549 inline void FullBCRS<T>::scale(t_real a_scale)
550 {
551  unsigned int n = m_blocks.size();
552  for(unsigned int i = 0; i < n; i++)
553  {
554  m_blocks[i] *= a_scale;
555  }
556 }
557 
558 template <typename T>
559 inline void UpperTriangularBCRS<T>::multiply(std::vector<t_v3> &a_b, const std::vector<t_v3> &a_x) const
560 {
561  unsigned int n = FullBCRS<T>::m_rowptr.size() - 1;
562  for(unsigned int i = 0; i < n; i++)
563  {
564  a_b[i] = FullBCRS<T>::m_zeroVector;
565  }
566  for(unsigned int i = 0; i < n; i++)
567  {
568  unsigned int k = FullBCRS<T>::m_rowptr[i];
569  a_b[i] += fe::multiply(FullBCRS<T>::m_blocks[k], a_x[i]);
570  for(k=k+1; k < FullBCRS<T>::m_rowptr[i+1]; k++)
571  {
572  a_b[i] += fe::multiply(FullBCRS<T>::m_blocks[k], a_x[FullBCRS<T>::m_colind[k]]);
573  a_b[FullBCRS<T>::m_colind[k]] += fe::transposeMultiply(FullBCRS<T>::m_blocks[k], a_x[i]);
574  }
575  }
576 }
577 
578 template <typename T>
579 inline void just_upper(Matrix<3,3,T> &a_m)
580 {
581  a_m(1,0) = 0.0;
582  a_m(2,0) = 0.0;
583  a_m(2,1) = 0.0;
584 }
585 
586 
587 template <typename T>
588 inline void FullBCRS<T>::sparse(const sp<FullBCRS> &a_other)
589 {
590  m_blocks.clear();
591  m_colind.clear();
592  m_rowptr.clear();
593  m_r = 0;
594  t_matrix zero;
595  setAll(zero, 0.0);
596  unsigned int n = a_other->m_rowptr.size()-1;
597  for(unsigned int i = 0; i < n; i++)
598  {
599  startrow();
600  int d = a_other->m_rowptr[i];
601  for(int k=d; k < a_other->m_rowptr[i+1]; k++)
602  {
603  int c = a_other->m_colind[k];
604  if(!(zero == a_other->m_blocks[k]))
605  {
606  insert(c, a_other->m_blocks[k]);
607  }
608  }
609  }
610  done();
611 }
612 
613 template <typename T>
614 inline void FullBCRS<T>::project(const sp<FullBCRS> &a_other)
615 {
616  unsigned int n = a_other->m_rowptr.size()-1;
617  for(unsigned int i = 0; i < n; i++)
618  {
619  int d = a_other->m_rowptr[i];
620  int d_self = m_rowptr[i];
621  for(int k=d; k < a_other->m_rowptr[i+1]; k++)
622  {
623  int c = a_other->m_colind[k];
624  for(int k_self=d_self; k_self < m_rowptr[i+1]; k_self++)
625  {
626  int c_self = m_colind[k_self];
627  if(c_self == c)
628  {
629  m_blocks[k_self] = a_other->m_blocks[k];
630  }
631  }
632  }
633  }
634 }
635 
636 template <typename T>
637 inline void UpperTriangularBCRS<T>::backSub(std::vector<t_v3> &a_x, const std::vector<t_v3> &a_b)
638 {
639 #if 0
640 t_matrix zero;
641 setAll(zero, 0.0);
642 #endif
643  unsigned int n = FullBCRS<T>::m_rowptr.size()-1;
644 //fe_fprintf(stderr, "root %d %d\n", n-1, m_colind[m_rowptr[n-1]]);
645  fe::backSub(FullBCRS<T>::m_blocks[FullBCRS<T>::m_rowptr[n-1]], a_x[n-1], a_b[n-1]);
646  t_v3 tmp;
647 
648  for(int i = n-2; i >= 0; i--)
649  {
650 //fe_fprintf(stderr, "bs i %d\n", i);
651  a_x[i] = a_b[i];
652  int d = FullBCRS<T>::m_rowptr[i];
653  for(unsigned int k=d+1; k < FullBCRS<T>::m_rowptr[i+1]; k++)
654  {
655  int c = FullBCRS<T>::m_colind[k];
656  fe::multiply(tmp, FullBCRS<T>::m_blocks[k], a_x[c]);
657  a_x[i] -= tmp;
658 #if 0
659 if(! (zero == m_blocks[k]))
660 {
661 fe_fprintf(stderr, "bs c %d -- %d <- T %d\n", c, i, k);
662 fprint(stderr, m_blocks[k]);
663 fe_fprintf(stderr, "\n");
664 }
665 #endif
666  //a_x[i] -= fe::transposeMultiply(m_blocks[k], a_x[c]);
667  }
668  fe::backSub(FullBCRS<T>::m_blocks[d], a_x[i], a_x[i]);
669 //fprint(stderr, m_blocks[d]);
670 //fe_fprintf(stderr, "\n");
671  }
672 }
673 
674 // operates on upper tri as if it is a transpose of the upper and is lower
675 template <typename T>
676 inline void UpperTriangularBCRS<T>::transposeForeSub(std::vector<t_v3> &a_x, const std::vector<t_v3> &a_b)
677 {
678 //t_matrix zero;
679 //setAll(zero, 0.0);
680  unsigned int n = FullBCRS<T>::m_rowptr.size()-1;
681  std::vector<unsigned int> ptr = FullBCRS<T>::m_rowptr;
682  fe::transposeForeSub(FullBCRS<T>::m_blocks[FullBCRS<T>::m_rowptr[0]], a_x[0], a_b[0]);
683  t_v3 tmp;
684 
685  for(unsigned int i = 1; i < n; i++)
686  {
687 //fe_fprintf(stderr, "fs i %d\n", i);
688  a_x[i] = a_b[i];
689  for(unsigned int m = 0; m < i; m++)
690  {
691  unsigned int c = FullBCRS<T>::m_colind[ptr[m]];
692  while(c < i && ptr[m] < FullBCRS<T>::m_rowptr[m+1]-1)
693  {
694  ptr[m]++;
695  c = FullBCRS<T>::m_colind[ptr[m]];
696  }
697  if(c == i)
698  {
699 //if(! (zero == m_blocks[ptr[m]]))
700 //{
701 //fe_fprintf(stderr, "fs m %d -- i %d <- T ptr[m] %d m_rowptr[m+1] %d\n", m, i, ptr[m], m_rowptr[m+1]);
702 //fprint(stderr, m_blocks[ptr[m]]);
703 //fe_fprintf(stderr, "\n");
704 //}
705  //fe::multiply(tmp, m_blocks[ptr[m]], a_x[m]);
706  //a_x[i] -= tmp;
707  a_x[i] -= fe::transposeMultiply(FullBCRS<T>::m_blocks[ptr[m]], a_x[m]);
708  }
709  }
710 //fe_fprintf(stderr, "fs m_rowptr[%d] %d col %d\n", i, m_rowptr[i], m_colind[m_rowptr[i]]);
711  fe::transposeForeSub(FullBCRS<T>::m_blocks[FullBCRS<T>::m_rowptr[i]], a_x[i], a_x[i]);
712 //fprint(stderr, m_blocks[m_rowptr[i]]);
713 //fe_fprintf(stderr, "\n");
714  }
715 }
716 
717 template <typename T>
718 inline void UpperTriangularBCRS<T>::backSolve(std::vector<t_v3> &a_x, const std::vector<t_v3> &a_b)
719 {
720  std::vector<t_v3> y(a_x.size());
721 
722  transposeForeSub(y, a_b);
723  backSub(a_x, y);
724 }
725 
726 #if 0
727 inline bool nanCheck(const t_matrix &aM, const char *aMsg)
728 {
729  for(unsigned int i = 0; i < 3; i++)
730  {
731  for(unsigned int j = 0; j < 3; j++)
732  {
733  if(isnan(aM(i,j)))
734  {
735  fe_fprintf(stderr, "M: %s\n", aMsg);
736  exit(101);
737  return false;
738  }
739  }
740  }
741  return true;
742 }
743 
744 inline bool nanCheck(const t_v3 &aV, const char *aMsg)
745 {
746  for(unsigned int i = 0; i < 3; i++)
747  {
748  if(isnan(aV[i]))
749  {
750  fe_fprintf(stderr, "V: %s\n", aMsg);
751  exit(101);
752  return false;
753  }
754  }
755  return true;
756 }
757 #endif
758 
759 template <typename T>
760 inline void UpperTriangularBCRS<T>::incompleteSqrt(void)
761 {
762  unsigned int d;
763  unsigned int n = FullBCRS<T>::m_rowptr.size()-1;
764 //fe_fprintf(stderr, "%d rows\n", n);
765  for(unsigned int k = 0; k < n-1; k++)
766  {
767  d = FullBCRS<T>::m_rowptr[k];
768  t_matrix z;
769  just_upper(FullBCRS<T>::m_blocks[d]);
770 //fe_fprintf(stderr, "WA %d attempt (row %d)\n", d, k);
771  //t_matrix tz = m_blocks[d] + transpose(m_blocks[d]);
772  //fprint(stderr, m_blocks[d]); fe_fprintf(stderr, "\n");
773  //nanCheck(m_blocks[d], "pre squareroot");
774  squareroot(FullBCRS<T>::m_blocks[d], FullBCRS<T>::m_blocks[d]);
775  //nanCheck(m_blocks[d], "post squareroot");
776  z = FullBCRS<T>::m_blocks[d];
777  //fprint(stderr, m_blocks[d]); fe_fprintf(stderr, "\n");
778 
779  just_upper(z);
780 
781  //t_v3 sv(0.2,0.5,0.7);
782  //t_v3 x, y;
783 
784  t_matrix inv_z;
785  invert(inv_z, z);
786 
787 #if 0
788 //fprint(stderr, tz); fe_fprintf(stderr, "\n");
789  t_matrix inv_tz;
790  invert(inv_tz, tz);
791  invert(tz, inv_tz);
792 //fprint(stderr, tz); fe_fprintf(stderr, "\n");
793 //fprint(stderr, inv_tz); fe_fprintf(stderr, "\n");
794  fe::multiply(x, inv_tz, sv);
795  //fe_fprintf(stderr, "via inv %f %f %f\n", x[0], x[1], x[2]);
796  t_v3 x_sv;
797  fe::multiply(x_sv, tz, x);
798  //fe_fprintf(stderr, "via inv back to sv %f %f %f\n", x_sv[0], x_sv[1], x_sv[2]);
799 
800  //fe::backSub(z, x, sv);
801  squareroot(z, tz);
802  t_matrix re_tz;
803  fe::multiply(re_tz, transpose(z), z);
804 //fprint(stderr, re_tz); fe_fprintf(stderr, "\n");
805  //setTranspose(z);
806  fe::transposeForeSub(z, y, sv);
807  fe::backSub(z, x, y);
808  //fe::transposeForeSub(z, x, y);
809  //fe_fprintf(stderr, "via bs %f %f %f\n", x[0], x[1], x[2]);
810  fe::multiply(x_sv, tz, x);
811  //fe_fprintf(stderr, "via bs back to sv %f %f %f\n", x_sv[0], x_sv[1], x_sv[2]);
812 //fe_fprintf(stderr, "-----------------------------------------------\n");
813 #endif
814 
815  for(unsigned int i=d+1; i < FullBCRS<T>::m_rowptr[k+1]; i++)
816  {
817 //fe_fprintf(stderr, "WB %d\n", i);
818  //fe::multiply(m_blocks[i], inv_z, m_blocks[i]);
819  //fprint(stderr, m_blocks[i]); fe_fprintf(stderr, "\n");
820  //fprint(stderr, inv_z); fe_fprintf(stderr, "\n");
821 
822  //nanCheck(m_blocks[i], "pre multiply A");
823  //nanCheck(inv_z, "pre multiply inv_z");
824  fe::multiply(z, FullBCRS<T>::m_blocks[i], inv_z);
825  setTranspose(z);
826 
827  //fe::multiply(m_blocks[i], inv_z, m_blocks[i]);
828  //setTranspose(m_blocks[i]);
829 
830  //fprint(stderr, m_blocks[i]); fe_fprintf(stderr, "\n");
831  //fprint(stderr, z); fe_fprintf(stderr, "\n");
832  FullBCRS<T>::m_blocks[i] = z;
833  }
834 
835  for(unsigned int i=d+1; i < FullBCRS<T>::m_rowptr[k+1]; i++)
836  {
837 //fe_fprintf(stderr, "i %d\n", i);
838  z = transpose(FullBCRS<T>::m_blocks[i]);
839  //z = m_blocks[i];
840  unsigned int h = FullBCRS<T>::m_colind[i];
841  unsigned int g = i;
842 
843  for(unsigned int j = FullBCRS<T>::m_rowptr[h]; j < FullBCRS<T>::m_rowptr[h+1]; j++)
844  {
845 //fe_fprintf(stderr, "j %d (%d < %d) (%d <= %d)\n", j, g, m_rowptr[k+1], m_colind[g+1], m_colind[j]);
846  for(; g < FullBCRS<T>::m_rowptr[k+1] && FullBCRS<T>::m_colind[g] <= FullBCRS<T>::m_colind[j]; g++)
847  {
848 //fe_fprintf(stderr, " g (%d < %d) (%d <= %d)\n", g, m_rowptr[k+1], m_colind[g+1], m_colind[j]);
849  if (FullBCRS<T>::m_colind[g] == FullBCRS<T>::m_colind[j])
850  {
851 //fe_fprintf(stderr, "WC[%d,%d] %d h: %d col: %d\n", i, g, j, h, m_colind[g]);
852  t_matrix t;
853  //nanCheck(m_blocks[i], "pre multiply B");
854  fe::multiply(t, z, FullBCRS<T>::m_blocks[g]);
855  //nanCheck(m_blocks[i], "post multiply BA");
856  //nanCheck(z, "post multiply BA");
857  //nanCheck(t, "post multiply BA");
858  setTranspose(t);
859  subtract(FullBCRS<T>::m_blocks[j], FullBCRS<T>::m_blocks[j], t);
860  //fprint(stderr, m_blocks[j]); fe_fprintf(stderr, "\n");
861  }
862  }
863  }
864  }
865  }
866  d = FullBCRS<T>::m_rowptr[n-1];
867 //fe_fprintf(stderr, "==------------------------ pre %d %d\n", d, m_blocks.size());
868  t_matrix tmp;
869  just_upper(FullBCRS<T>::m_blocks[d]);
870  //fprint(stderr, m_blocks[d]);
871  //fe_fprintf(stderr, "\n");
872  squareroot(tmp, FullBCRS<T>::m_blocks[d]);
873 //fe_fprintf(stderr, "==------------------------\n");
874  //fprint(stderr, tmp);
875  //fe_fprintf(stderr, "\n");
876  FullBCRS<T>::m_blocks[d] = tmp;
877 //fe_fprintf(stderr, "==------------------------ post\n");
878 }
879 
880 
881 inline FullVMR::FullVMR(void)
882 {
883 }
884 
885 inline FullVMR::~FullVMR(void)
886 {
887 }
888 
889 inline void UpperTriangularVMR::multiply(std::vector<t_v3> &a_b, const std::vector<t_v3> &a_x) const
890 {
891  unsigned int n = a_b.size();
892  for(unsigned int i = 0; i < n; i++)
893  {
894  a_b[i] = zeroVector;
895  }
896 
897  unsigned int r,c;
898  n = m_rows.size();
899  for(unsigned int i = 0; i < n; i++)
900  {
901  t_row::const_iterator i_row = m_rows[i].begin();
902  r = i;
903  const SpatialMatrix &block = (*i_row).second;
904  assert((*i_row).first == r);
905  a_b[r] += fe::multiply(block, a_x[r]);
906  i_row++;
907 
908  for(;i_row != m_rows[i].end(); i_row++)
909  {
910  const SpatialMatrix &block = (*i_row).second;
911  c = (*i_row).first;
912 
913  a_b[r] += fe::multiply(block, a_x[c]);
914  a_b[c] += fe::transposeMultiply(block, a_x[r]);
915  }
916  }
917 }
918 
919 inline void FullVMR::multiply(std::vector<t_v3> &a_b, const std::vector<t_v3> &a_x) const
920 {
921  unsigned int r,c;
922  unsigned int n = m_rows.size();
923  for(unsigned int i = 0; i < n; i++)
924  {
925  a_b[i] = zeroVector;
926  r = i;
927  for(t_row::const_iterator i_row = m_rows[i].begin();
928  i_row != m_rows[i].end(); i_row++)
929  {
930  const SpatialMatrix &block = (*i_row).second;
931  c = (*i_row).first;
932  a_b[r] += fe::multiply(block, a_x[c]);
933  }
934  }
935 }
936 
937 inline void FullVMR::zero(void)
938 {
939  unsigned int n = m_rows.size();
940  for(unsigned int i = 0; i < n; i++)
941  {
942  for(t_row::iterator i_row = m_rows[i].begin();
943  i_row != m_rows[i].end(); i_row++)
944  {
945  SpatialMatrix &block = (*i_row).second;
946  setAll(block, 0.0);
947  }
948  }
949 }
950 
951 inline void FullVMR::scale(t_real a_scale)
952 {
953  unsigned int n = m_rows.size();
954  for(unsigned int i = 0; i < n; i++)
955  {
956  for(t_row::iterator i_row = m_rows[i].begin();
957  i_row != m_rows[i].end(); i_row++)
958  {
959  SpatialMatrix &block = (*i_row).second;
960  fe::multiply(block, a_scale);
961  }
962  }
963 }
964 
965 inline void FullVMR::clear(void)
966 {
967  m_rows.clear();
968 }
969 
970 inline void FullVMR::setRows(unsigned int a_count)
971 {
972  m_rows.resize(a_count);
973 }
974 
975 inline unsigned int FullVMR::rows(void)
976 {
977  return m_rows.size();
978 }
979 
980 inline UpperTriangularVMR::t_row &FullVMR::row(unsigned int a_index)
981 {
982  return m_rows[a_index];
983 }
984 
985 inline SpatialMatrix &FullVMR::block(unsigned int a_i, unsigned int a_j)
986 {
987  //assert(a_i < m_rows.size());
988  t_row::iterator i_row = m_rows[a_i].find(a_j);
989  if(i_row != m_rows[a_i].end())
990  {
991  return i_row->second;
992  }
993  SpatialMatrix &m = m_rows[a_i][a_j];
994  setAll(m, 0.0);
995  return m;
996 }
997 
998 
999 } /* namespace ext */
1000 } /* namespace fe */
1001 
1002 #endif /* __SparseMatrix3x3_h__ */
1003 
fe::ext::FullVMR
Vector Map Rows (vector of row maps)
Definition: SparseMatrix3x3.h:31
fe::String::c_str
const FESTRING_I8 * c_str(void) const
Return the contents of the 8-bit buffer cast as signed bytes.
Definition: String.h:352
fe
kernel
Definition: namespace.dox:3
fe::invert
Matrix< 4, 4, T > & invert(Matrix< 4, 4, T > &a_inverted, const Matrix< 4, 4, T > &a_matrix)
4x4 full matrix inversion
Definition: Matrix.h:1033
fe::ext::UpperTriangularVMR
Upper Triangular Vector Map Rows (vector of row maps)
Definition: SparseMatrix3x3.h:57
fe::String
Automatically reference-counted string container.
Definition: String.h:128
fe::project
void project(Vector< 4, T > &a_r, const Matrix< 4, 4, T > &a_m, const Vector< 4, T > &a_v)
project vector through matrix. divides by transformed w
Definition: Matrix.h:1182
fe::sp
Intrusive Smart Pointer.
Definition: src/core/ptr.h:53
fe::squareroot
void squareroot(Matrix< M, N, T > &a_U, const Matrix< M, N, T > &a_A)
Square root of a matrix.
Definition: Matrix.h:567
Generated by   doxygen 1.8.13