NESSi  v1.0.2 The NonEquilibrium Systems Simulation Library
cntr_elements.hpp
Go to the documentation of this file.
1 #ifndef CNTR_ELEMENTS_H
2 #define CNTR_ELEMENTS_H
3
4 #include "eigen_map.hpp"
5 #include "linalg.hpp"
6
7 namespace cntr {
8
9 // a matrix-matrix multiplication, like gemm in blas
10
11 #define LARGESIZE (-1) // Fall back to dynamic size
12 #define BLAS_SIZE 4 // use blas for larger matrices
13 // #define USE_BLAS
14 /* #######################################################################################
15 #
16 # general matrix operations ...
17 #
18 ########################################################################################*/
19 #define CPLX std::complex<T>
20 #define VALUE_TYPE T
21 template<typename T,int DIM1,int DIM2>
23 inline void element_set_zero(int size1,int size2,CPLX * z){
24  memset(z,0,sizeof(CPLX)*size1*size2);
25 }
27 template<typename T,int DIM1,int DIM2>
28 inline void element_set(int size1,int size2,CPLX *z,CPLX *z1) {
29  memcpy(z,z1,sizeof(CPLX)*size1*size2);
30 }
32 template<typename T,int DIM1,int DIM2>
33 inline void element_swap(int size1,int size2,CPLX *z,CPLX *z1) {
34  int nm=size1*size2,i;
35  CPLX x;
36  for(i=0;i<nm;i++){x=z[i];z[i]=z1[i];z1[i]=x;}
37 }
38 // hermitian conjugate
40 template<typename T,int DIM1,int DIM2>
41 inline void element_conj(int size1,int size2,CPLX * z,CPLX *z1){
42  // z1 is of size [ size2 x size1 ]
43  int l,m;
44  CPLX temp;
45  for(l=0;l<size1;l++){
46  for(m=0;m<size2;m++){
47  temp=z1[m*size1+l];
48  z[l*size2 + m ] = CPLX(temp.real(),-temp.imag());
49  }
50  }
51 }
53 template<typename T,int DIM1,int DIM2>
54 inline void element_conj(int size1,int size2,CPLX * z){
55  CPLX *z1 = new CPLX [size1*size2];
56  element_set<T,DIM2,DIM1>(size2,size1,z1,z);
57  element_conj<T,DIM1,DIM2>(size1,size2,z,z1);
58  delete [] z1;
59 }
61 template<typename T,int DIM1,int DIM2>
62 inline void element_minusconj(int size1,int size2,CPLX * z,CPLX *z1){
63  // z1 is of size [ size2 x size1 ]
64  int l,m;
65  CPLX temp;
66  for(l=0;l<size1;l++){
67  for(m=0;m<size2;m++){
68  temp=z1[m*size1+l];
69  z[l*size2 + m ] = CPLX(-temp.real(),temp.imag());
70  }
71  }
72 }
74 template<typename T,int DIM1,int DIM2>
75 inline void element_minusconj(int size1,int size2,CPLX * z){
76  CPLX *z1 = new CPLX [size1*size2];
77  element_set<T,DIM2,DIM1>(size2,size1,z1,z);
78  element_minusconj<T,DIM1,DIM2>(size1,size2,z,z1);
79  delete [] z1;
80 }
82 template<typename T,int DIM1,int DIM2>
83 inline T element_norm2(int size1,int size2,CPLX * z){
84  int i,nm=size1*size2;
85  T r=0.0;
86  for(i=0;i<nm;i++) r += z[i].real()*z[i].real()+z[i].imag()*z[i].imag();
87  return sqrt(r);
88 }
89 // scalar multiplication
91 template<typename T,int DIM1,int DIM2>
92 inline void element_smul(int size1,int size2,CPLX *z,T z0) {
93  int nm=size1*size2,i;
94  for(i=0;i<nm;i++) z[i] *= z0;
95 }
97 template<typename T,int DIM1,int DIM2>
98 inline void element_smul(int size1,int size2,CPLX *z,CPLX z0) {
99  int nm=size1*size2,i;
100  for(i=0;i<nm;i++) z[i] *= z0;
101 }
102 // multiplication
104 template<typename T,int DIM1,int DIM2,int DIM3>
105 inline void element_mult(int size1,int size2,int size3,CPLX *z,CPLX *z1,CPLX *z2){
106  element_map<DIM1, DIM2>(size1, size2, z) =
107  element_map<DIM1, DIM3>(size1, size3, z1) *
108  element_map<DIM3, DIM2>(size3, size2, z2);
109 }
110 // increment with product and other
112 template<typename T,int DIM1,int DIM2>
113 inline void element_incr(int size1,int size2,CPLX *z,CPLX *z1) {
114  if (z != z1) {
115  element_map<DIM1, DIM2>(size1, size2, z).noalias() +=
116  element_map<DIM1, DIM2>(size1, size2, z1);
117  } else {
118  element_map<DIM1, DIM2>(size1, size2, z) +=
119  element_map<DIM1, DIM2>(size1, size2, z1);
120  }
121 }
123 template<typename T,int DIM1,int DIM2>
124 inline void element_incr(int size1,int size2,CPLX *z,CPLX alpha,CPLX *z1) {
125  if (z != z1) {
126  element_map<DIM1, DIM2>(size1, size2, z).noalias() +=
127  alpha * element_map<DIM1, DIM2>(size1, size2, z1);
128  } else {
129  element_map<DIM1, DIM2>(size1, size2, z) +=
130  alpha * element_map<DIM1, DIM2>(size1, size2, z1);
131  }
132 }
134 template<typename T,int DIM1,int DIM2,int DIM3>
135 inline void element_incr(int size1,int size2,int size3,CPLX *z,CPLX *z1,CPLX *z2){
136  if (z != z1 && z != z2) {
137  element_map<DIM1, DIM2>(size1, size2, z).noalias() +=
138  element_map<DIM1, DIM3>(size1, size3, z1) *
139  element_map<DIM3, DIM2>(size3, size2, z2);
140  } else {
141  element_map<DIM1, DIM2>(size1, size2, z) +=
142  element_map<DIM1, DIM3>(size1, size3, z1) *
143  element_map<DIM3, DIM2>(size3, size2, z2);
144  }
145 }
147 template<typename T,int DIM1,int DIM2,int DIM3>
148 inline void element_incr(int size1,int size2,int size3,CPLX *z,CPLX alpha,CPLX *z1,CPLX *z2){
149  if (z != z1 && z != z2) {
150  element_map<DIM1, DIM2>(size1, size2, z).noalias() +=
151  alpha * (element_map<DIM1, DIM3>(size1, size3, z1) *
152  element_map<DIM3, DIM2>(size3, size2, z2));
153  } else {
154  element_map<DIM1, DIM2>(size1, size2, z) +=
155  alpha * (element_map<DIM1, DIM3>(size1, size3, z1) *
156  element_map<DIM3, DIM2>(size3, size2, z2));
157  }
158 }
159 // specialization for DIM=1, T=double
160 #undef CPLX
161 #define CPLX std::complex<double>
162 #undef VALUE_TYPE
163 #define VALUE_TYPE double
164 template<> inline void element_set_zero<VALUE_TYPE,1,1>(int size1,int size2,CPLX * z){
166  *z=0.0;
167 }
169 template<> inline void element_set<VALUE_TYPE,1,1>(int size1,int size2,CPLX *z,CPLX *z1) {
170  *z=*z1;
171 }
173 template<> inline void element_swap<VALUE_TYPE,1,1>(int size1,int size2,CPLX *z,CPLX *z1) {
174  CPLX x;
175  x=*z;*z=*z1;*z1=x;
176 }
178 template<> inline void element_conj<VALUE_TYPE,1,1>(int size1,int size2,CPLX * z,CPLX *z1){
179  CPLX temp=*z1;
180  //temp.imag()=-temp.imag();
181  temp = std::conj( temp );
182  *z=temp;
183 }
185 template<> inline void element_conj<VALUE_TYPE,1,1>(int size1,int size2,CPLX * z){
186  CPLX temp=*z;
187  //temp.imag()=-temp.imag();
188  temp = std::conj( temp );
189  *z=temp;
190 }
192 template<> inline void element_minusconj<VALUE_TYPE,1,1>(int size1,int size2,CPLX * z,CPLX *z1){
193  CPLX temp=*z1;
194  //temp.real()=-temp.real();
195  temp = -std::conj( temp );
196  *z=temp;
197 }
199 template<> inline void element_minusconj<VALUE_TYPE,1,1>(int size1,int size2,CPLX * z){
200  CPLX temp=*z;
201  //temp.real()=-temp.real();
202  temp = -std::conj( temp );
203  *z=temp;
204 }
206 template<> inline VALUE_TYPE element_norm2<VALUE_TYPE,1,1>(int size1,int size2,CPLX * z){
207  CPLX x=*z;
208  return sqrt(x.real()*x.real()+x.imag()*x.imag());
209 }
211 template<> inline void element_smul<VALUE_TYPE,1,1>(int size1,int size2,CPLX *z,VALUE_TYPE alpha) {
212  *z=*z*alpha;
213 }
215 template<> inline void element_smul<VALUE_TYPE,1,1>(int size1,int size2,CPLX *z,CPLX alpha){
216  *z=*z*alpha;
217 }
219 template<> inline void element_mult<VALUE_TYPE,1,1,1>(int size1,int size2,int size3,CPLX *z,CPLX *z1,CPLX *z2){
220  *z=(*z1)*(*z2);
221 }
223 template<> inline void element_incr<VALUE_TYPE,1,1>(int size1,int size2,CPLX *z,CPLX *z1) {
224  *z=*z+*z1;
225 }
227 template<> inline void element_incr<VALUE_TYPE,1,1>(int size1,int size2,CPLX *z,CPLX alpha,CPLX *z1) {
228  *z=*z+alpha*(*z1);
229 }
231 template<> inline void element_incr<VALUE_TYPE,1,1,1>(int size1,int size2,int size3,CPLX *z,CPLX *z1,CPLX *z2){
232  *z=*z+(*z1)*(*z2);
233 }
235 template<> inline void element_incr<VALUE_TYPE,1,1,1>(int size1,int size2,int size3,CPLX *z,CPLX alpha,CPLX *z1,CPLX *z2){
236  *z=*z+alpha*(*z1)*(*z2);
237 }
238 // specialization for DIM=1, T=float
239 #undef CPLX
240 #define CPLX std::complex<float>
241 #undef VALUE_TYPE
242 #define VALUE_TYPE float
243 template<> inline void element_set_zero<VALUE_TYPE,1,1>(int size1,int size2,CPLX * z){
245  *z=0.0;
246 }
248 template<> inline void element_set<VALUE_TYPE,1,1>(int size1,int size2,CPLX *z,CPLX *z1) {
249  *z=*z1;
250 }
252 template<> inline void element_swap<VALUE_TYPE,1,1>(int size1,int size2,CPLX *z,CPLX *z1) {
253  CPLX x;
254  x=*z;*z=*z1;*z1=x;
255 }
257 template<> inline void element_conj<VALUE_TYPE,1,1>(int size1,int size2,CPLX * z,CPLX *z1){
258  CPLX temp=*z1;
259  //temp.imag()=-temp.imag();
260  temp = std::conj( temp );
261  *z=temp;
262 }
264 template<> inline void element_conj<VALUE_TYPE,1,1>(int size1,int size2,CPLX * z){
265  CPLX temp=*z;
266  //temp.imag()=-temp.imag();
267  temp = std::conj( temp );
268  *z=temp;
269 }
271 template<> inline void element_minusconj<VALUE_TYPE,1,1>(int size1,int size2,CPLX * z,CPLX *z1){
272  CPLX temp=*z1;
273  //temp.real()=-temp.real();
274  temp = -std::conj( temp );
275  *z=temp;
276 }
278 template<> inline void element_minusconj<VALUE_TYPE,1,1>(int size1,int size2,CPLX * z){
279  CPLX temp=*z;
280  //temp.real()=-temp.real();
281  temp = -std::conj( temp );
282  *z=temp;
283 }
285 template<> inline VALUE_TYPE element_norm2<VALUE_TYPE,1,1>(int size1,int size2,CPLX * z){
286  CPLX x=*z;
287  return sqrt(x.real()*x.real()+x.imag()*x.imag());
288 }
290 template<> inline void element_smul<VALUE_TYPE,1,1>(int size1,int size2,CPLX *z,VALUE_TYPE alpha) {
291  *z=*z*alpha;
292 }
294 template<> inline void element_smul<VALUE_TYPE,1,1>(int size1,int size2,CPLX *z,CPLX alpha) {
295  *z=*z*alpha;
296 }
298 template<> inline void element_mult<VALUE_TYPE,1,1,1>(int size1,int size2,int size3,CPLX *z,CPLX *z1,CPLX *z2){
299  *z=(*z1)*(*z2);
300 }
302 template<> inline void element_incr<VALUE_TYPE,1,1>(int size1,int size2,CPLX *z,CPLX *z1) {
303  *z=*z+*z1;
304 }
306 template<> inline void element_incr<VALUE_TYPE,1,1>(int size1,int size2,CPLX *z,CPLX alpha,CPLX *z1) {
307  *z=*z+alpha*(*z1);
308 }
310 template<> inline void element_incr<VALUE_TYPE,1,1,1>(int size1,int size2,int size3,CPLX *z,CPLX *z1,CPLX *z2){
311  *z=*z+(*z1)*(*z2);
312 }
314 template<> inline void element_incr<VALUE_TYPE,1,1,1>(int size1,int size2,int size3,CPLX *z,CPLX alpha,CPLX *z1,CPLX *z2){
315  *z=*z+alpha*(*z1)*(*z2);
316 }
317
318 /* #######################################################################################
319 #
320 # square matrix operations ... partly derived from the general ones,
321 # from which they inherit the specialization
322 #
323 ########################################################################################*/
324 #undef CPLX
325 #define CPLX std::complex<T>
326 #undef VALUE_TYPE
327 #define VALUE_TYPE T
328 template<typename T,int DIM>
330 inline void element_set_zero(int size1,CPLX * z){
331  element_set_zero<T,DIM,DIM>(size1,size1,z);
332 }
334 template<typename T,int DIM>
335 inline void element_set(int size1,CPLX *z,CPLX z0) {
336  element_set_zero<T,DIM>(size1,z);
337  for(int i=0;i<size1;i++) z[i*size1+i]=z0;
338 }
340 template<typename T,int DIM>
341 inline void element_set(int size1,CPLX *z,CPLX *z1) {
342  element_set<T,DIM,DIM>(size1,size1,z,z1);
343 }
345 template<typename T,int DIM>
346 inline void element_swap(int size1,CPLX *z,CPLX *z1) {
347  element_swap<T,DIM,DIM>(size1,size1,z,z1);
348 }
350 template<typename T,int DIM>
351 inline CPLX element_trace(int size1,CPLX * z){
352  CPLX result=0.0;
353  for(int i=0;i<size1;i++) result += z[i*size1+i];
354  return result;
355 }
356 // hermitian conjugate
358 template<typename T,int DIM>
359 inline void element_conj(int size1,CPLX * z,CPLX *z1){
360  element_conj<T,DIM,DIM>(size1,size1,z,z1);
361 }
363 template<typename T,int DIM>
364 inline void element_conj(int size1,CPLX * z){
365  element_conj<T,DIM,DIM>(size1,size1,z);
366 }
368 template<typename T,int DIM>
369 inline void element_minusconj(int size1,CPLX * z,CPLX *z1){
370  element_minusconj<T,DIM,DIM>(size1,size1,z,z1);
371 }
373 template<typename T,int DIM>
374 inline void element_minusconj(int size1,CPLX * z){
375  element_minusconj<T,DIM,DIM>(size1,size1,z);
376 }
378 template<typename T,int DIM>
379 inline T element_norm2(int size1,CPLX * z){
380  return element_norm2<T,DIM,DIM>(size1,size1,z);
381 }
382 // scalar multiplication
384 template<typename T,int DIM>
385 inline void element_smul(int size1,CPLX *z,T z0) {
386  element_smul<T,DIM,DIM>(size1,size1,z,z0);
387 }
389 template<typename T,int DIM>
390 inline void element_smul(int size1,CPLX *z,CPLX z0) {
391  element_smul<T,DIM,DIM>(size1,size1,z,z0);
392 }
393 // multiplication
395 template<typename T,int DIM>
396 inline void element_mult(int size1,CPLX *z,CPLX *z1,CPLX *z2){
397  element_mult<T,DIM,DIM,DIM>(size1,size1,size1,z,z1,z2);
398 }
399 // increment with product etc.
401 template<typename T,int DIM>
402 inline void element_incr(int size1,CPLX *z,CPLX z0) {
403  for(int i=0;i<size1;i++) z[i*size1+i]+=z0;
404 }
406 template<typename T,int DIM>
407 inline void element_incr(int size1,CPLX *z,CPLX *z1) {
408  element_incr<T,DIM,DIM>(size1,size1,z,z1);
409 }
411 template<typename T,int DIM>
412 inline void element_incr(int size1,CPLX *z,CPLX alpha,CPLX *z1) {
413  element_incr<T,DIM,DIM>(size1,size1,z,alpha,z1);
414 }
416 template<typename T,int DIM>
417 inline void element_incr(int size1,CPLX *z,CPLX *z1,CPLX *z2){
418  element_incr<T,DIM,DIM,DIM>(size1,size1,size1,z,z1,z2);
419 }
421 template<typename T,int DIM>
422 inline void element_incr(int size1,CPLX *z,CPLX alpha,CPLX *z1,CPLX *z2){
423  element_incr<T,DIM,DIM,DIM>(size1,size1,size1,z,alpha,z1,z2);
424 }
425 // matrix inverse , use double precision because i am lazy to look
426 // for the float inverse routine
428 template<typename T,int DIM>
429 inline void element_inverse(int size1,CPLX *z,CPLX *z1) {
430  int i,nm=size1*size1;
431  std::complex<double> *zd = new std::complex<double> [size1*size1];
432  std::complex<double> *z1d = new std::complex<double> [size1*size1];
433  for(i=0;i<nm;i++) z1d[i]= (std::complex<double>) z1[i];
434  linalg::cplx_matrix_inverse(z1d,zd,size1);
435  for(i=0;i<nm;i++) z[i]= (CPLX)zd[i];
436  delete [] zd;
437  delete [] z1d;
438 }
439 // solution of linear equation M_ij X_j = Q_i (right) and X_j M_ij =Q_i (left)
440 // make this more efficient at a later time, allowing for a reuse of the matrix M (?)
442 template<typename T,int DIM>
443 inline void element_linsolve_right(int size1,int n,CPLX *X,CPLX *M,CPLX *Q){
444  // this uses currently always double precision
445  int size=size1,llen=size*n,l,s2=size*size;
446  int p,q,m;
447  std::complex<double> *mtemp = new std::complex<double> [llen*llen];
448  std::complex<double> *qtemp = new std::complex<double> [n*s2];
449  std::complex<double> *xtemp = new std::complex<double> [n*s2];
450  // set up the matrix - reshuffle elements line by line
451  for(l=0;l<n;l++){
452  for(m=0;m<n;m++){
453  for(p=0;p<size;p++){
454  for(q=0;q<size;q++){
455  mtemp[ (l*size + p)*llen + m*size+q] = (std::complex<double>) M[ (l*n+m)*s2+p*size+q];
456  }
457  }
458  }
459  }
460  for(l=0;l<n;l++){
461  for(p=0;p<size;p++){
462  for(q=0;q<size;q++){
463  qtemp[ q*n*size + l*size+p] = (std::complex<double>) Q[ l*s2+p*size+q];
464  }
465  }
466  }
467  linalg::cplx_sq_solve_many(mtemp,qtemp,xtemp,llen,size);
468  for(l=0;l<n;l++){
469  for(p=0;p<size;p++){
470  for(q=0;q<size;q++){
471  X[ l*s2+p*size+q]= (CPLX) xtemp[ q*n*size + l*size+p];
472  }
473  }
474  }
475  delete [] xtemp;
476  delete [] qtemp;
477  delete [] mtemp;
478 }
480 template<typename T,int DIM>
481 inline void element_linsolve_right(int size1,CPLX *X,CPLX *M,CPLX *Q){
482  element_linsolve_right<T,DIM>(size1,1,X,M,Q);
483 }
485 template<typename T,int DIM>
486 inline void element_linsolve_left(int size1,int n,CPLX *X,CPLX *M,CPLX *Q){
487  int l,m,element_size=size1*size1;
488  // transpose all elements
489  for(l=0;l<n;l++) for(m=0;m<n;m++) element_conj<T,DIM>(size1, M + (l*n+m)*element_size);
490  for(l=0;l<n;l++) element_conj<T,DIM>(size1, Q + l*element_size);
491  element_linsolve_right<T,DIM>(size1,n,X,M,Q);
492  for(l=0;l<n;l++) element_conj<T,DIM>(size1, X + l*element_size);
493 }
495 template<typename T,int DIM>
496 inline void element_linsolve_left(int size1,CPLX *X,CPLX *M,CPLX *Q){
497  element_linsolve_left<T,DIM>(size1,1,X,M,Q);
498 }
499
500 // specialization of square matrix funtions which are not derived from general
501 // for DIM=1,T=double
502 #undef CPLX
503 #define CPLX std::complex<double>
504 #undef VALUE_TYPE
505 #define VALUE_TYPE double
506 template<> inline void element_set<VALUE_TYPE,1>(int size1,CPLX *z,CPLX z0) {
508  *z=z0;
509 }
511 template<> inline CPLX element_trace<VALUE_TYPE,1>(int size1,CPLX * z){
512  return *z;
513 }
515 template<> inline void element_incr<VALUE_TYPE,1>(int size1,CPLX *z,CPLX z0) {
516  *z=*z+z0;
517 }
519 template<> inline void element_linsolve_right<VALUE_TYPE,1>(int size1,CPLX *X,CPLX *M,CPLX *Q) {
520  *X = *Q / (*M);
521 }
523 template<> inline void element_linsolve_left<VALUE_TYPE,1>(int size1,CPLX *X,CPLX *M,CPLX *Q) {
524  *X = *Q / (*M);
525 }
527 template<> inline void element_inverse<VALUE_TYPE,1>(int size1,CPLX *z,CPLX *z1) {
528  *z= 1.0/(*z1);
529 }
530 // specialization of square matrix funtions which are not derived from general
531 // for DIM=1,T=float
532 #undef CPLX
533 #define CPLX std::complex<float>
534 #undef VALUE_TYPE
535 #define VALUE_TYPE float
536 template<> inline void element_set<VALUE_TYPE,1>(int size1,CPLX *z,CPLX z0) {
538  *z=z0;
539 }
541 template<> inline CPLX element_trace<VALUE_TYPE,1>(int size1,CPLX * z){
542  return *z;
543 }
545 template<> inline void element_incr<VALUE_TYPE,1>(int size1,CPLX *z,CPLX z0) {
546  *z=*z+z0;
547 }
549 template<> inline void element_linsolve_right<VALUE_TYPE,1>(int size1,CPLX *X,CPLX *M,CPLX *Q) {
550  *X = *Q / (*M);
551 }
553 template<> inline void element_linsolve_left<VALUE_TYPE,1>(int size1,CPLX *X,CPLX *M,CPLX *Q) {
554  *X = *Q / (*M);
555 }
557 template<> inline void element_inverse<VALUE_TYPE,1>(int size1,CPLX *z,CPLX *z1) {
558  *z= CPLX(1.0,0.0)/(*z1);
559 }
560
561
562 #undef CPLX
563 #undef VALUE_TYPE
564
565 } // namespace cntr
566
567 #endif // CNTR_ELEMENTS_H
void cplx_matrix_inverse(void *a, void *x, int n)
Evaluate the inverse matrix of a complex matrix .
void cplx_sq_solve_many(void *a, void *b, void *x, int dim, int d)
Solve a linear equation .