* Eigenvalues of a Hermitian matrix @ 2001-12-19 13:20 Andrew W Steiner 2001-12-19 13:20 ` Brian Gough 0 siblings, 1 reply; 6+ messages in thread From: Andrew W Steiner @ 2001-12-19 13:20 UTC (permalink / raw) To: gsl-discuss Hello! I want to calculate the derivative of a function of the eigenvalues of a particular Hermitian matrix and so I need the eigenvalues of several matrices whose entries are similar. Is there any way to provide an inital guess to the gsl_eigen_herm routine if one is known? From my understanding of the algorithm, the answer is no. In that case, could anyone suggest an alternative? Thanks, Andrew W. Steiner Nuclear Theory Group Department of Physics and Astronomy State University of New York at Stony Brook http://tonic.physics.sunysb.edu/~asteiner ^ permalink raw reply [flat|nested] 6+ messages in thread
* Re: Eigenvalues of a Hermitian matrix 2001-12-19 13:20 Eigenvalues of a Hermitian matrix Andrew W Steiner @ 2001-12-19 13:20 ` Brian Gough 2001-12-19 13:20 ` Brian Gough 0 siblings, 1 reply; 6+ messages in thread From: Brian Gough @ 2001-12-19 13:20 UTC (permalink / raw) To: Andrew W Steiner; +Cc: gsl-discuss Andrew W Steiner writes: > I want to calculate the derivative of a function of the > eigenvalues of a particular Hermitian matrix and so I need the eigenvalues > of several matrices whose entries are similar. Is there any way to > provide an inital guess to the gsl_eigen_herm routine if one is known? > From my understanding of the algorithm, the answer is no. In that > case, could anyone suggest an alternative? If you premultiply your perturbed matrices by V' (where V is the unitary matrix of unperturbed eigenvectors) they will be close to diagonal. This might improve the computation of the eigenvalues. The same applies if you have only a guess for V. It probably won't make much difference overall, but if it does let me know. regards Brian ^ permalink raw reply [flat|nested] 6+ messages in thread
* Re: Eigenvalues of a Hermitian matrix 2001-12-19 13:20 ` Brian Gough @ 2001-12-19 13:20 ` Brian Gough 0 siblings, 0 replies; 6+ messages in thread From: Brian Gough @ 2001-12-19 13:20 UTC (permalink / raw) To: Andrew W Steiner, gsl-discuss Brian Gough writes: > > If you premultiply your perturbed matrices by V' (where V is the ^^^^^^^^^^^^ Should be "apply a similarity transformation, V' M V", not premultiply, of course. ^ permalink raw reply [flat|nested] 6+ messages in thread
* Eigenvalues of a Hermitian matrix @ 2001-12-19 13:20 Andrew W Steiner 2001-12-19 13:20 ` Brian Gough 0 siblings, 1 reply; 6+ messages in thread From: Andrew W Steiner @ 2001-12-19 13:20 UTC (permalink / raw) To: gsl-discuss Hello! The following code fails to calculate the eigenvalues of the second matrix. (It's only 70 lines, and even that would be shortened if I wouldn't have tried to be so pedantic with my style). Is my usage correct, or is there possibly an error in my compilation of the libraries? I am using Red Hat Linux 6.2 with g++ (version egcs-2.91.66). Could someone possibly double check this for me? The code is supposed to print out and calculate the eigenvalues of two simple matrices but never leaves the second call to gsl_eigen_herm. Thanks, Andrew W. Steiner Nuclear Theory Group Department of Physics and Astronomy State University of New York at Stony Brook http://tonic.physics.sunysb.edu/~asteiner ----------------------------------------------------------------------- #include <iostream> #include <gsl/gsl_math.h> #include <gsl/gsl_eigen.h> #include <gsl/gsl_complex_math.h> #include <gsl/gsl_complex.h> int main(void) { double p; int i,j,k; gsl_complex z, *zp; double real, imag; gsl_matrix_complex *n; gsl_vector *eval; gsl_eigen_herm_workspace *w; zp=new gsl_complex; n=gsl_matrix_complex_alloc(4,4); eval = gsl_vector_alloc (4); w=gsl_eigen_herm_alloc(4); for(i=0;i<4;i++) { for(j=0;j<4;j++) { real=0.0; imag=0.0; if (i==j) real=((double)i); GSL_SET_COMPLEX(zp,real,imag); gsl_matrix_complex_set(n,i,j,*zp); cout << "(" << real << " " << imag << ") "; } cout << endl; } gsl_eigen_herm(n, eval, w); for(k=0;k<4;k++) cout << gsl_vector_get(eval,k) << " "; cout << endl; gsl_eigen_herm_free(w); gsl_matrix_complex_free(n); gsl_vector_free(eval); w=gsl_eigen_herm_alloc(4); n=gsl_matrix_complex_alloc(4,4); eval = gsl_vector_alloc (4); p=1.0; for(i=0;i<4;i++) { for(j=0;j<4;j++) { real=0.0; imag=0.0; if (i==j && j==1) real+=p; if ((i==0 && j==2) || i==2 && j==0) real-=p; if ((i==1 && j==3) || i==3 && j==1) real+=p; GSL_SET_COMPLEX(zp,real,imag); gsl_matrix_complex_set(n,i,j,*zp); cout << "(" << real << " " << imag << ") "; } cout << endl; } gsl_eigen_herm(n, eval, w); for(k=0;k<4;k++) cout << gsl_vector_get(eval,k) << " "; cout << endl; delete zp; return 0; } ^ permalink raw reply [flat|nested] 6+ messages in thread
* Re: Eigenvalues of a Hermitian matrix 2001-12-19 13:20 Andrew W Steiner @ 2001-12-19 13:20 ` Brian Gough 2001-12-19 13:20 ` Brian Gough 0 siblings, 1 reply; 6+ messages in thread From: Brian Gough @ 2001-12-19 13:20 UTC (permalink / raw) To: Andrew W Steiner; +Cc: gsl-discuss Andrew W Steiner writes: > The following code fails to calculate the eigenvalues of the > second matrix. (It's only 70 lines, and even that would be > shortened if I wouldn't have tried to be so pedantic with my > style). Is my usage correct, or is there possibly an error in my > compilation of the libraries? I am using Red Hat Linux 6.2 with g++ > (version egcs-2.91.66). > Could someone possibly double check this for me? The code is > supposed to print out and calculate the eigenvalues of two simple > matrices but never leaves the second call to gsl_eigen_herm. Thanks for the bug report. Your program is valid. The eigenvalue routine was going into an infinite loop. The following patch to linalg/qrstep.c fixes it. I will check it into CVS. After applying the patch, $ ./a.out .... (0 0) (0 0) (-1 0) (0 0) (0 0) (1 0) (0 0) (1 0) (-1 0) (0 0) (0 0) (0 0) (0 0) (1 0) (0 0) (0 0) 1 -1 1.61803 -0.618034 Index: qrstep.c =================================================================== RCS file: /cvs/gsl/gsl/eigen/qrstep.c,v retrieving revision 1.2 diff -u -r1.2 qrstep.c --- qrstep.c 2001/06/12 12:21:46 1.2 +++ qrstep.c 2001/08/02 17:15:23 @@ -60,18 +60,13 @@ double mu; - if (dt > 0) + if (dt >= 0) { mu = tb - (tab * tab) / (dt + hypot (dt, tab)); } - else if (dt < 0) - { - mu = tb + (tab * tab) / ((-dt) + hypot (dt, tab)); - } else { - /* FIXME: what is the best choice of mu for dt == 0 ? */ - mu = tb; + mu = tb + (tab * tab) / ((-dt) + hypot (dt, tab)); } return mu; ^ permalink raw reply [flat|nested] 6+ messages in thread
* Re: Eigenvalues of a Hermitian matrix 2001-12-19 13:20 ` Brian Gough @ 2001-12-19 13:20 ` Brian Gough 0 siblings, 0 replies; 6+ messages in thread From: Brian Gough @ 2001-12-19 13:20 UTC (permalink / raw) To: Andrew W Steiner, gsl-discuss Brian Gough writes: > The following patch to linalg/qrstep.c fixes it. I will check it into > CVS. There was also a bug in the accumulation of the eigenvectors. I've added this to the test suite. Index: symmv.c =================================================================== RCS file: /cvs/gsl/gsl/eigen/symmv.c,v retrieving revision 1.2 retrieving revision 1.3 diff -u -r1.2 -r1.3 --- symmv.c 2001/07/01 21:48:38 1.2 +++ symmv.c 2001/08/02 19:10:53 1.3 @@ -192,10 +192,10 @@ for (k = 0; k < N; k++) { - double qki = gsl_matrix_get (evec, k, i); - double qkj = gsl_matrix_get (evec, k, i + 1); - gsl_matrix_set (evec, k, i, qki * c - qkj * s); - gsl_matrix_set (evec, k, i + 1, qki * s + qkj * c); + double qki = gsl_matrix_get (evec, k, a + i); + double qkj = gsl_matrix_get (evec, k, a + i + 1); + gsl_matrix_set (evec, k, a + i, qki * c - qkj * s); + gsl_matrix_set (evec, k, a + i + 1, qki * s + qkj * c); } } Index: hermv.c =================================================================== RCS file: /cvs/gsl/gsl/eigen/hermv.c,v retrieving revision 1.3 retrieving revision 1.4 diff -u -r1.3 -r1.4 --- hermv.c 2001/07/01 21:48:38 1.3 +++ hermv.c 2001/08/02 19:10:53 1.4 @@ -215,8 +215,8 @@ for (k = 0; k < N; k++) { - gsl_complex qki = gsl_matrix_complex_get (evec, k, i); - gsl_complex qkj = gsl_matrix_complex_get (evec, k, i + 1); + gsl_complex qki = gsl_matrix_complex_get (evec, k, a + i); + gsl_complex qkj = gsl_matrix_complex_get (evec, k, a + i + 1); /* qki <= qki * c - qkj * s */ /* qkj <= qki * s + qkj * c */ gsl_complex x1 = gsl_complex_mul_real(qki, c); @@ -228,8 +228,8 @@ gsl_complex qqki = gsl_complex_add(x1, y1); gsl_complex qqkj = gsl_complex_add(x2, y2); - gsl_matrix_complex_set (evec, k, i, qqki); - gsl_matrix_complex_set (evec, k, i + 1, qqkj); + gsl_matrix_complex_set (evec, k, a + i, qqki); + gsl_matrix_complex_set (evec, k, a + i + 1, qqkj); } } ^ permalink raw reply [flat|nested] 6+ messages in thread
end of thread, other threads:[~2001-12-19 13:20 UTC | newest] Thread overview: 6+ messages (download: mbox.gz / follow: Atom feed) -- links below jump to the message on this page -- 2001-12-19 13:20 Eigenvalues of a Hermitian matrix Andrew W Steiner 2001-12-19 13:20 ` Brian Gough 2001-12-19 13:20 ` Brian Gough -- strict thread matches above, loose matches on Subject: below -- 2001-12-19 13:20 Andrew W Steiner 2001-12-19 13:20 ` Brian Gough 2001-12-19 13:20 ` Brian Gough
This is a public inbox, see mirroring instructions for how to clone and mirror all data and code used for this inbox; as well as URLs for read-only IMAP folder(s) and NNTP newsgroup(s).