home *** CD-ROM | disk | FTP | other *** search
- //------------------------------------------------------------------------------
- //
- // esort
- //
- // Syntax: a=esort(p)
- //
- // This routine sorts the complex continuous eigenvalues in descending
- // order. It sorts based on the real part. The unstable eigenvalues (positive
- // real part) will appear first.
- //
- // The routine passes back the sorted eigenvalues and the cooresponding
- // indices in a list:
- //
- // a.s = sorted eigenvalues
- // a.ndx = index
- //
- // Version JBL 940917
- //------------------------------------------------------------------------------
-
- esort = function(p)
- {
- local(s,ndx,i,k,swap,a,ic,j)
-
- if (p.nr == 1) {
- p=p.';
- }
-
- a=sort(-real(p));
- s=a.val;
- ndx=a.ind';
-
- for (i in 1:p.nc) {
- ic=ndx[;i];
- for (j in 1:p.nr) {
- s[j;i]=p[ndx[j;i];i];
- }
- }
-
- for (i in 1:p.nc) {
- k=1;
- while (k < p.nr) {
- if (imag(s[k;i]) != 0) {
- if (imag(s[k;i]) < 0) {
- s[k:k+1;i]=conj(s[k:k+1;i]);
- swap=ndx[k;i];
- ndx[k;i]=ndx[k+1;i];
- ndx[k+1;i]=swap;
- }
- k=k+2;
- else
- k=k+1;
- }
- }
- }
-
- return << s=s; ndx=ndx >>
- };
-