cgv
Loading...
Searching...
No Matches
eig.h
1#pragma once
2
3#include "functions.h"
4#include "mat.h"
5#include "diag_mat.h"
6
7#include <limits>
8#include <algorithm>
9
10
11#include "vec.h"
12#include <complex>
13
14namespace cgv {
15namespace math {
16
17
18//if x is an eigen vector of A the rayleigh quotient returns the corresponding eigenvalue
19template<typename T>
20T rayleigh_quotient(mat<T>& A,vec<T>&x)
21{
22 return dot(x,A*x)/dot(x,x);
23}
24
25
26
27
28
29template <typename T>
30void rot(mat<T> &a, const T s, const T tau, const int i,
31
32 const int j, const int k, const int l)
33
34 {
35
36 T g=a(i,j);
37
38 T h=a(k,l);
39
40 a(i,j)=g-s*(h+g*tau);
41
42 a(k,l)=h+s*(g-h*tau);
43
44 }
45
46
47
48template <typename T>
49void eigsrt(diag_mat<T> &d, mat<T> &v)
50
51{
52
53 unsigned k;
54
55 unsigned n=d.size();
56
57 for (unsigned i=0;i<n-1;i++)
58
59 {
60
61 T p=d(k=i);
62
63 for (unsigned j=i;j<n;j++)
64
65 if (d(j) >= p) p=d(k=j);
66
67 if (k != i)
68
69 {
70
71 d(k)=d(i);
72 d(i)=p;
73
74
75 for (unsigned j=0;j<n;j++)
76 {
77
78 p=v(j,i);
79 v(j,i)=v(j,k);
80 v(j,k)=p;
81
82 }
83 }
84 }
85}
86
91template <typename T>
92bool eig_sym(const mat<T> &a, mat<T> &v, diag_mat<T> &d,bool ordering=true, unsigned maxiter=50)
93{
94 mat<T>aa = a;
95 unsigned n =aa.nrows();
96
97 v.identity(n);
98
99 d.resize(n);
100
101 unsigned nrot=0;
102
103 const T eps = std::numeric_limits<T>::epsilon();
104
105 T tresh,theta,tau,t,sm,s,h,g,c;
106
107 vec<T> b(n),z;
108 z.zeros(n);
109 unsigned ip,iq;
110
111 for(unsigned i = 0; i < n; i++)
112 d(i)=b(i)=aa(i,i);
113
114 for (unsigned i=1;i<=maxiter;i++)
115
116 {
117
118 sm=(T)0.0;
119
120
121
122 for (ip=0;ip<n-1;ip++)
123
124 {
125
126 for (iq=ip+1;iq<n;iq++)
127
128 sm += std::abs(aa(ip,iq));
129
130 }
131
132 if (sm == (T)0.0)
133
134 {
135
136 if(ordering)
137
138 eigsrt(d,v);
139
140 return true;
141
142 }
143
144 if (i < 4)
145
146 tresh=(T)(0.2*sm/(n*n));
147
148 else
149
150 tresh=(T)0.0;
151
152 for (ip=0;ip<n-1;ip++)
153
154 {
155
156 for (iq=ip+1;iq<n;iq++)
157
158 {
159
160 g=((T)100.0)*std::abs(aa(ip,iq));
161
162 if (i > 4 && g <= eps*std::abs(d(ip)) && g <= eps*std::abs(d(iq)))
163
164 aa(ip,iq)=(T)0.0;
165
166 else if (std::abs(aa(ip,iq)) > tresh)
167
168 {
169
170 h=d(iq)-d(ip);
171
172 if (g <= eps*std::abs(h))
173
174 t=(aa(ip,iq))/h;
175
176 else {
177
178 theta=(T)(0.5*h/(aa(ip,iq)));
179
180 t=(T)(1.0/(std::abs(theta)+sqrt(1.0+theta*theta)));
181
182 if (theta < 0.0) t = -t;
183
184 }
185
186 c=(T)(1.0/sqrt(1+t*t));
187
188 s=t*c;
189
190 tau=(T)(s/(1.0+c));
191
192 h=t*aa(ip,iq);
193
194 z(ip) -= h;
195
196 z(iq) += h;
197
198 d(ip) -= h;
199
200 d(iq) += h;
201
202 aa(ip,iq)=(T)0.0;
203
204 for (unsigned j=0;j<ip;j++)
205
206 rot(aa,s,tau,j,ip,j,iq);
207
208 for (unsigned j=ip+1;j<iq;j++)
209
210 rot(aa,s,tau,ip,j,j,iq);
211
212 for (unsigned j=iq+1;j<n;j++)
213
214 rot(aa,s,tau,ip,j,iq,j);
215
216 for (unsigned j=0;j<n;j++)
217
218 rot(v,s,tau,j,ip,j,iq);
219
220 ++nrot;
221
222 }
223
224 }
225
226 }
227
228 for (ip=0;ip<n;ip++) {
229
230 b(ip) += z(ip);
231
232 d(ip)=b(ip);
233
234 z(ip)=(T)0.0;
235
236 }
237
238 }
239
240 return false;
241 //Too many iterations in routine jacobi
242
243}
244
245
246
247template <typename T>
248struct Unsymmeig
249{
250 int n;
251 cgv::math::mat<T> a,zz;
252 cgv::math::diag_mat<std::complex<T> > wri;
253 cgv::math::vec<T> scale;
254 cgv::math::vec<int> perm;
255 bool yesvecs,hessen;
256
257 Unsymmeig(const cgv::math::mat<T>&aa, bool yesvec=true, bool hessenb=false) :
258 n(aa.nrows()), a(aa), zz(n,n,0.0), wri(n), scale(n), perm(n),
259 yesvecs(yesvec), hessen(hessenb)
260 {
261 scale.ones();
262 balance();
263 if (!hessen) elmhes();
264 if (yesvecs) {
265 for (int i=0;i<n;i++)
266 zz(i,i)=1.0;
267 if (!hessen) eltran();
268 hqr2();
269 balbak();
270 sortvecs();
271 } else {
272 hqr();
273 sort();
274 }
275 }
276 void balance();
277 void elmhes();
278 void eltran();
279 void hqr();
280 void hqr2();
281 void balbak();
282 void sort();
283 void sortvecs();
284
285
286};
287template <typename T>
288T SIGN(const T &a, const T &b)
289 {return (T)(b >= 0 ? (a >= 0 ? a : -a) : (a >= 0 ? -a : a));}
290template <typename T>
291void Unsymmeig<T>::balance()
292{
293 const T RADIX = std::numeric_limits<T>::radix;
294 bool done=false;
295 T sqrdx=RADIX*RADIX;
296 while (!done) {
297 done=true;
298 for (int i=0;i<n;i++) {
299 T r=0.0,c=0.0;
300 for (int j=0;j<n;j++)
301 if (j != i) {
302 c += std::abs(a(j,i));
303 r += std::abs(a(i,j));
304 }
305 if (c != 0.0 && r != 0.0) {
306 T g=r/RADIX;
307 T f=1.0;
308 T s=c+r;
309 while (c<g) {
310 f *= RADIX;
311 c *= sqrdx;
312 }
313 g=r*RADIX;
314 while (c>g) {
315 f /= RADIX;
316 c /= sqrdx;
317 }
318 if ((c+r)/f < 0.95*s)
319 {
320 done=false;
321 g=1.0/f;
322 scale[i] *= f;
323 for (int j=0;j<n;j++) a(i,j) *= g;
324 for (int j=0;j<n;j++) a(j,i) *= f;
325 }
326 }
327 }
328 }
329}
330template <typename T>
331void Unsymmeig<T>::balbak()
332{
333 for (int i=0;i<n;i++)
334 for (int j=0;j<n;j++)
335 zz(i,j) *= scale[i];
336}
337
338template <typename T>
339void Unsymmeig<T>::elmhes()
340{
341 for (int m=1;m<n-1;m++) {
342 T x=0.0;
343 int i=m;
344 for (int j=m;j<n;j++)
345 {
346 if (std::abs(a(j,m-1)) > std::abs(x))
347 {
348 x=a(j,m-1);
349 i=j;
350 }
351 }
352 perm[m]=i;
353 if (i != m)
354 {
355 for (int j=m-1;j<n;j++) std::swap(a(i,j),a(m,j));
356 for (int j=0;j<n;j++) std::swap(a(j,i),a(j,m));
357 }
358 if (x != 0.0)
359 {
360 for (i=m+1;i<n;i++)
361 {
362 T y=a(i,m-1);
363 if (y != 0.0)
364 {
365 y /= x;
366 a(i,m-1)=y;
367 for (int j=m;j<n;j++) a(i,j) -= y*a(m,j);
368 for (int j=0;j<n;j++) a(j,m) += y*a(j,i);
369 }
370 }
371 }
372 }
373}
374
375template <typename T>
376void Unsymmeig<T>::eltran()
377{
378 for (int mp=n-2;mp>0;mp--)
379 {
380 for (int k=mp+1;k<n;k++)
381 zz(k,mp)=a(k,mp-1);
382 int i=perm[mp];
383 if (i != mp) {
384 for (int j=mp;j<n;j++) {
385 zz(mp,j)=zz(i,j);
386 zz(i,j)=0.0;
387 }
388 zz(i,mp)=1.0;
389 }
390 }
391}
392
393template <typename T>
394void Unsymmeig<T>::hqr()
395{
396 int nn,m,l,k,j,its,i,mmin;
397 T z,y,x,w,v,u,t,s,r,q,p,anorm=0.0;
398
399 const T EPS=std::numeric_limits<T>::epsilon();
400 for (i=0;i<n;i++)
401 for (j=std::max(i-1,0);j<n;j++)
402 anorm += std::abs(a(i,j));
403 nn=n-1;
404 t=0.0;
405 while (nn >= 0) {
406 its=0;
407 do {
408 for (l=nn;l>0;l--) {
409 s=std::abs(a(l-1,l-1))+std::abs(a(l,l));
410 if (s == 0.0) s=anorm;
411 if (std::abs(a(l,l-1)) <= EPS*s)
412 {
413 a(l,l-1) = 0.0;
414 break;
415 }
416 }
417 x=a(nn,nn);
418 if (l == nn) {
419 wri[nn--]=x+t;
420 } else {
421 y=a(nn-1,nn-1);
422 w=a(nn,nn-1)*a(nn-1,nn);
423 if (l == nn-1)
424 {
425 p=0.5*(y-x);
426 q=p*p+w;
427 z=sqrt(std::abs(q));
428 x += t;
429 if (q >= 0.0) {
430 z=p+SIGN(z,p);
431 wri[nn-1]=wri[nn]=x+z;
432 if (z != 0.0) wri[nn]=x-w/z;
433 } else {
434 wri[nn]=std::complex<T>(x+p,-z);
435 wri[nn-1]=conj(wri[nn]);
436 }
437 nn -= 2;
438 } else {
439 if (its == 30) throw("Too many iterations in hqr");
440 if (its == 10 || its == 20) {
441 t += x;
442 for (i=0;i<nn+1;i++) a(i,i) -= x;
443 s=std::abs(a(nn,nn-1))+std::abs(a(nn-1,nn-2));
444 y=x=0.75*s;
445 w = -0.4375*s*s;
446 }
447 ++its;
448 for (m=nn-2;m>=l;m--) {
449 z=a(m,m);
450 r=x-z;
451 s=y-z;
452 p=(r*s-w)/a(m+1,m)+a(m,m+1);
453 q=a(m+1,m+1)-z-r-s;
454 r=a(m+2,m+1);
455 s=std::abs(p)+std::abs(q)+std::abs(r);
456 p /= s;
457 q /= s;
458 r /= s;
459 if (m == l) break;
460 u=std::abs(a(m,m-1))*(std::abs(q)+std::abs(r));
461 v=std::abs(p)*(std::abs(a(m-1,m-1))+std::abs(z)+std::abs(a(m+1,m+1)));
462 if (u <= EPS*v) break;
463 }
464 for (i=m;i<nn-1;i++) {
465 a(i+2,i)=0.0;
466 if (i != m) a(i+2,i-1)=0.0;
467 }
468 for (k=m;k<nn;k++) {
469 if (k != m) {
470 p=a(k,k-1);
471 q=a(k+1,k-1);
472 r=0.0;
473 if (k+1 != nn) r=a(k+2,k-1);
474 if ((x=std::abs(p)+std::abs(q)+std::abs(r)) != 0.0) {
475 p /= x;
476 q /= x;
477 r /= x;
478 }
479 }
480 if ((s=SIGN(sqrt(p*p+q*q+r*r),p)) != 0.0) {
481 if (k == m) {
482 if (l != m)
483 a(k,k-1) = -a(k,k-1);
484 } else
485 a(k,k-1) = -s*x;
486 p += s;
487 x=p/s;
488 y=q/s;
489 z=r/s;
490 q /= p;
491 r /= p;
492 for (j=k;j<nn+1;j++)
493 {
494 p=a(k,j)+q*a(k+1,j);
495 if (k+1 != nn) {
496 p += r*a(k+2,j);
497 a(k+2,j) -= p*z;
498 }
499 a(k+1,j) -= p*y;
500 a(k,j) -= p*x;
501 }
502 mmin = nn < k+3 ? nn : k+3;
503 for (i=l;i<mmin+1;i++)
504 {
505 p=x*a(i,k)+y*a(i,k+1);
506 if (k+1 != nn)
507 {
508 p += z*a(i,k+2);
509 a(i,k+2) -= p*r;
510 }
511 a(i,k+1) -= p*q;
512 a(i,k) -= p;
513 }
514 }
515 }
516 }
517 }
518 } while (l+1 < nn);
519 }
520}
521
522template <typename T>
523void Unsymmeig<T>::hqr2()
524{
525 int nn,m,l,k,j,its,i,mmin,na;
526 T z,y,x,w,v,u,t,s,r,q,p,anorm=0.0,ra,sa,vr,vi;
527
528 const T EPS=std::numeric_limits<T>::epsilon();
529 for (i=0;i<n;i++)
530 for (j=std::max(i-1,0);j<n;j++)
531 anorm += std::abs(a(i,j));
532 nn=n-1;
533 t=0.0;
534 while (nn >= 0) {
535 its=0;
536 do {
537 for (l=nn;l>0;l--) {
538 s=std::abs(a(l-1,l-1))+std::abs(a(l,l));
539 if (s == 0.0) s=anorm;
540 if (std::abs(a(l,l-1)) <= EPS*s) {
541 a(l,l-1) = 0.0;
542 break;
543 }
544 }
545 x=a(nn,nn);
546 if (l == nn) {
547 wri(nn)=a(nn,nn)=x+t;
548 nn--;
549 } else {
550 y=a(nn-1,nn-1);
551 w=a(nn,nn-1)*a(nn-1,nn);
552 if (l == nn-1) {
553 p=0.5*(y-x);
554 q=p*p+w;
555 z=sqrt(std::abs(q));
556 x += t;
557 a(nn,nn)=x;
558 a(nn-1,nn-1)=y+t;
559 if (q >= 0.0) {
560 z=p+SIGN(z,p);
561 wri[nn-1]=wri[nn]=x+z;
562 if (z != 0.0) wri[nn]=x-w/z;
563 x=a(nn,nn-1);
564 s=std::abs(x)+std::abs(z);
565 p=x/s;
566 q=z/s;
567 r=sqrt(p*p+q*q);
568 p /= r;
569 q /= r;
570 for (j=nn-1;j<n;j++) {
571 z=a(nn-1,j);
572 a(nn-1,j)=q*z+p*a(nn,j);
573 a(nn,j)=q*a(nn,j)-p*z;
574 }
575 for (i=0;i<=nn;i++) {
576 z=a(i,nn-1);
577 a(i,nn-1)=q*z+p*a(i,nn);
578 a(i,nn)=q*a(i,nn)-p*z;
579 }
580 for (i=0;i<n;i++) {
581 z=zz(i,nn-1);
582 zz(i,nn-1)=q*z+p*zz(i,nn);
583 zz(i,nn)=q*zz(i,nn)-p*z;
584 }
585 } else {
586 wri[nn]=std::complex<T>(x+p,-z);
587 wri[nn-1]=std::conj(wri[nn]);
588 }
589 nn -= 2;
590 } else {
591 if (its == 30) throw("Too many iterations in hqr");
592 if (its == 10 || its == 20)
593 {
594 t += x;
595 for (i=0;i<nn+1;i++) a(i,i) -= x;
596 s=std::abs(a(nn,nn-1))+std::abs(a(nn-1,nn-2));
597 y=x=0.75*s;
598 w = -0.4375*s*s;
599 }
600 ++its;
601 for (m=nn-2;m>=l;m--)
602 {
603 z=a(m,m);
604 r=x-z;
605 s=y-z;
606 p=(r*s-w)/a(m+1,m)+a(m,m+1);
607 q=a(m+1,m+1)-z-r-s;
608 r=a(m+2,m+1);
609 s=std::abs(p)+std::abs(q)+std::abs(r);
610 p /= s;
611 q /= s;
612 r /= s;
613 if (m == l) break;
614 u=std::abs(a(m,m-1))*(std::abs(q)+std::abs(r));
615 v=std::abs(p)*(std::abs(a(m-1,m-1))+std::abs(z)+std::abs(a(m+1,m+1)));
616 if (u <= EPS*v) break;
617 }
618 for (i=m;i<nn-1;i++) {
619 a(i+2,i)=0.0;
620 if (i != m) a(i+2,i-1)=0.0;
621 }
622 for (k=m;k<nn;k++) {
623 if (k != m) {
624 p=a(k,k-1);
625 q=a(k+1,k-1);
626 r=0.0;
627 if (k+1 != nn) r=a(k+2,k-1);
628 if ((x=std::abs(p)+std::abs(q)+std::abs(r)) != 0.0) {
629 p /= x;
630 q /= x;
631 r /= x;
632 }
633 }
634 if ((s=SIGN(sqrt(p*p+q*q+r*r),p)) != 0.0) {
635 if (k == m) {
636 if (l != m)
637 a(k,k-1) = -a(k,k-1);
638 } else
639 a(k,k-1) = -s*x;
640 p += s;
641 x=p/s;
642 y=q/s;
643 z=r/s;
644 q /= p;
645 r /= p;
646 for (j=k;j<n;j++) {
647 p=a(k,j)+q*a(k+1,j);
648 if (k+1 != nn) {
649 p += r*a(k+2,j);
650 a(k+2,j) -= p*z;
651 }
652 a(k+1,j) -= p*y;
653 a(k,j) -= p*x;
654 }
655 mmin = nn < k+3 ? nn : k+3;
656 for (i=0;i<mmin+1;i++) {
657 p=x*a(i,k)+y*a(i,k+1);
658 if (k+1 != nn) {
659 p += z*a(i,k+2);
660 a(i,k+2) -= p*r;
661 }
662 a(i,k+1) -= p*q;
663 a(i,k) -= p;
664 }
665 for (i=0; i<n; i++) {
666 p=x*zz(i,k)+y*zz(i,k+1);
667 if (k+1 != nn) {
668 p += z*zz(i,k+2);
669 zz(i,k+2) -= p*r;
670 }
671 zz(i,k+1) -= p*q;
672 zz(i,k) -= p;
673 }
674 }
675 }
676 }
677 }
678 } while (l+1 < nn);
679 }
680 if (anorm != 0.0) {
681 for (nn=n-1;nn>=0;nn--) {
682 p=real(wri[nn]);
683 q=imag(wri[nn]);
684 na=nn-1;
685 if (q == 0.0) {
686 m=nn;
687 a(nn,nn)=1.0;
688 for (i=nn-1;i>=0;i--) {
689 w=a(i,i)-p;
690 r=0.0;
691 for (j=m;j<=nn;j++)
692 r += a(i,j)*a(j,nn);
693 if (imag(wri[i]) < 0.0) {
694 z=w;
695 s=r;
696 } else {
697 m=i;
698
699 if (imag(wri[i]) == 0.0) {
700 t=w;
701 if (t == 0.0)
702 t=EPS*anorm;
703 a(i,nn)=-r/t;
704 } else {
705 x=a(i,i+1);
706 y=a(i+1,i);
707 q=cgv::math::sqr(real(wri[i])-p)+cgv::math::sqr(imag(wri[i]));
708 t=(x*s-z*r)/q;
709 a(i,nn)=t;
710 if (std::abs(x) > std::abs(z))
711 a(i+1,nn)=(-r-w*t)/x;
712 else
713 a(i+1,nn)=(-s-y*t)/z;
714 }
715 t=std::abs(a(i,nn));
716 if (EPS*t*t > 1)
717 for (j=i;j<=nn;j++)
718 a(j,nn) /= t;
719 }
720 }
721 } else if (q < 0.0) {
722 m=na;
723 if (std::abs(a(nn,na)) > std::abs(a(na,nn))) {
724 a(na,na)=q/a(nn,na);
725 a(na,nn)=-(a(nn,nn)-p)/a(nn,na);
726 } else {
727 std::complex<T> temp=std::complex<T>(0.0,-a(na,nn))/std::complex<T>(a(na,na)-p,q);
728 a(na,na)=real(temp);
729 a(na,nn)=imag(temp);
730 }
731 a(nn,na)=0.0;
732 a(nn,nn)=1.0;
733 for (i=nn-2;i>=0;i--) {
734 w=a(i,i)-p;
735 ra=sa=0.0;
736 for (j=m;j<=nn;j++) {
737 ra += a(i,j)*a(j,na);
738 sa += a(i,j)*a(j,nn);
739 }
740 if (imag(wri[i]) < 0.0) {
741 z=w;
742 r=ra;
743 s=sa;
744 } else {
745 m=i;
746 if (imag(wri[i]) == 0.0) {
747 std::complex<T> temp = std::complex<T>(-ra,-sa)/std::complex<T>(w,q);
748 a(i,na)=real(temp);
749 a(i,nn)=imag(temp);
750 } else {
751 x=a(i,i+1);
752 y=a(i+1,i);
753 vr=cgv::math::sqr(real(wri(i))-p)+cgv::math::sqr(imag(wri(i)))-q*q;
754 vi=2.0*q*(real(wri(i))-p);
755 if (vr == 0.0 && vi == 0.0)
756 vr=EPS*anorm*(std::abs(w)+std::abs(q)+std::abs(x)+std::abs(y)+std::abs(z));
757 std::complex<T> temp=std::complex<T>(x*r-z*ra+q*sa,x*s-z*sa-q*ra)/
758 std::complex<T>(vr,vi);
759 a(i,na)=real(temp);
760 a(i,nn)=imag(temp);
761 if (std::abs(x) > std::abs(z)+std::abs(q)) {
762 a(i+1,na)=(-ra-w*a(i,na)+q*a(i,nn))/x;
763 a(i+1,nn)=(-sa-w*a(i,nn)-q*a(i,na))/x;
764 } else {
765 std::complex<T> temp=std::complex<T>(-r-y*a(i,na),-s-y*a(i,nn))/
766 std::complex<T>(z,q);
767 a(i+1,na)=real(temp);
768 a(i+1,nn)=imag(temp);
769 }
770 }
771 }
772 t=std::max(std::abs(a(i,na)),std::abs(a(i,nn)));
773 if (EPS*t*t > 1)
774 for (j=i;j<=nn;j++) {
775 a(j,na) /= t;
776 a(j,nn) /= t;
777 }
778 }
779 }
780 }
781 for (j=n-1;j>=0;j--)
782 for (i=0;i<n;i++) {
783 z=0.0;
784 for (k=0;k<=j;k++)
785 z += zz(i,k)*a(k,j);
786 zz(i,j)=z;
787 }
788 }
789}
790template <typename T>
791void Unsymmeig<T>::sort()
792{
793 int i;
794 for (int j=1;j<n;j++)
795 {
796 std::complex<T> x=wri[j];
797 for (i=j-1;i>=0;i--) {
798 if (std::real(wri[i]) >= real(x)) break;
799 wri[i+1]=wri[i];
800 }
801 wri[i+1]=x;
802 }
803}
804
805template <typename T>
806void Unsymmeig<T>::sortvecs()
807{
808 int i;
809 cgv::math::vec<T> temp(n);
810 for (int j=1;j<n;j++) {
811 std::complex<T> x=wri[j];
812 for (int k=0;k<n;k++)
813 temp[k]=zz(k,j);
814 for (i=j-1;i>=0;i--) {
815 if (real(wri[i]) >= std::real(x)) break;
816 wri[i+1]=wri[i];
817 for (int k=0;k<n;k++)
818 zz(k,i+1)=zz(k,i);
819 }
820 wri[i+1]=x;
821 for (int k=0;k<n;k++)
822 zz(k,i+1)=temp[k];
823 }
824
825}
826
827
828
829
830
831
832//compute eigenvalues of a which is a real unsymmetric matrix
833//if a is still in hessenberg form set hessenb to true (default is false)
834template <typename T>
835void eig_unsym(const cgv::math::mat<T> &a,cgv::math::diag_mat<std::complex<T> >& eigvals, bool hessenb=false)
836{
837 Unsymmeig<T> eigsolver((cgv::math::mat<T>)a,false,hessenb);
838 eigvals = eigsolver.wri;
839
840}
841
842
843//compute eigenvectors and eigenvalues of matrix a which is a real unsymmetric matrix
844//if a is still in hessenberg form set hessenb to true (default is false)
845//eigenvectors are not normalized
846template <typename T>
847void eig_unsym(const cgv::math::mat<T>& a,cgv::math::mat<T>& eigvecs,cgv::math::diag_mat<std::complex<T> >& eigvals,bool normalize=true, bool hessenb=false)
848{
849 cgv::math::mat<T> A=a;
850 Unsymmeig<T> eigsolver(A,true,hessenb);
851 eigvals = eigsolver.wri;
852 if(normalize)
853 {
854 for(unsigned i = 0; i < eigsolver.zz.ncols();i++)
855 eigsolver.zz.set_col(i,cgv::math::normalize(eigsolver.zz.col(i)));
856 }
857 eigvecs=eigsolver.zz;
858
859}
860
861}}
cgv::math::mat
A matrix type (full column major storage) The matrix can be loaded directly into OpenGL without need ...
Definition mat.h:208
cgv::math::mat::nrows
unsigned nrows() const
number of rows
Definition mat.h:540
cgv::math::vec
A column vector class.
Definition vec.h:28
cgv::math::vec::ones
void ones()
fill the vector with ones
Definition vec.h:476
cgv
the cgv namespace
Definition print.h:11
vr
the vr namespace for virtual reality support
Definition gl_vr_display.cxx:5
cgv::math::diag_mat
A diagonal matrix type which internally stores the elements on the main diagonal in a vector.
Definition diag_mat.h:16