toon-members
[Top][All Lists]

## Re: [Toon-members] NonSymEigen for TooN, and (SVD full....)

 From: E. Rosten Subject: Re: [Toon-members] NonSymEigen for TooN, and (SVD full....) Date: Thu, 28 Jan 2010 20:50:56 +0000 (GMT) User-agent: Alpine 2.00 (LSU 1167 2008-08-23)

```On Thu, 28 Jan 2010, Gabriel Nützi wrote:

```
```Hello guys

```
I have adapted the SymEigen for Non-Symetric Matrix, because there is no Complex Support in TooN we have just Imaginary and real parts for the EigenValues and Eigenvectors
```
```
In principle, TooN does support Vector<std::complex<double> >. I think some minor modification will be required, but not much.
```
```
One particular problem is that LAPACK seems to want to return the real and imaginary parts of the eigenvalues in different arrays. I find this a little odd seeing as FORTRAN supports complex types.
```
```
To make that work with TooN, we would need a new Allocator type which stores two pointers, and upon indexing constructs and returns a complex<double>. This would require that Vector also be parameterized on the data pointer. Perhaps the underlying allocator type could provide this type. Such parameterization would allow a Vector to be constructed off a pair<double*,double*> rather than a complex<double>*.
```
```
An intermediate solution would be to simply copy the two Vectors in to a single Vector<complex>. I can't see that this would be a problem since copying is a fast O(N) operation and eigendecomposition is a rather slow O(N^3) operation.
```
```
Is returning a Vector<complex<double> > better than returning two vectors?
```
```
Making that change would allow one to be able to decouple the underlying datatype from the datatype returned by operator[]. This would allow one to take the conjugate transpose of a matrix cleanly.
```
Thoughts?

```
--> NonSymEigen.h (it's not finished (I leave the documentation and add ons for this class to other guys because I am not the pro here) but would be nice to have that in the later version of TooN :-). It has been tested and it works.... :-) uses Lapack dgeev_ and sgeev_
```
```
We can leave out all the helper functions until someone needs them. Better to have missing features than to put in untested code.
```

```
```TestCode :

Matrix<4,4> M;
...
NonSymEigen<> myEig(M);

Matrix<>a= myEig.get_evector();
Vector<>b= myEig.get_evaluesIm();
Vector<>c= myEig.get_evaluesRe();

====================================

```
--> The SVD.h file contains an advanced version for full or not full decomposition, so we dont have to Pad with zeros....
```
```
This seems to be directly supported by LAPACK, so it is probably more efficient this way than padding the matrix. If it is not more efficient, then the easiest solution is to have a separate FullSVD class which wraps SVD and pads the matrix.
```
```
```--> Dont know what you think about that but its handy :-)
```
```
```
I think it would be better to specify it using a boolean template parameter. This would ensure that the returned matrices are statically sized. This could be made compatible with existing code by renaming the SVD class, and making two wrapper classes, SVD and FullSVD.
```
That aside, I have a couple of qureies:

```
Even with full_compute=true, it appears that get_vertical() depends on the input matrix. Given that, it appeats that get_U and get_VT are missing some if(full_compute) code paths.
```
-Ed

```
```Thanks for the support and maintaining the class :-)

Gabriel

```
```
--
(You can't go wrong with psycho-rats.)(http://mi.eng.cam.ac.uk/~er258)

/d{def}def/f{/Times s selectfont}d/s{11}d/r{roll}d f 2/m{moveto}d -1
r 230 350 m 0 1 179{ 1 index show 88 rotate 4 mul 0 rmoveto}for/s 12
d f pop 235 420 translate 0 0 moveto 1 2 scale show showpage```