1c-----------------------------------------------------------------------
2c\BeginDoc
3c
4c\Name: ssgets
5c
6c\Description:
7c  Given the eigenvalues of the symmetric tridiagonal matrix H,
8c  computes the NP shifts AMU that are zeros of the polynomial of
9c  degree NP which filters out components of the unwanted eigenvectors
10c  corresponding to the AMU's based on some given criteria.
11c
12c  NOTE: This is called even in the case of user specified shifts in
13c  order to sort the eigenvalues, and error bounds of H for later use.
14c
15c\Usage:
16c  call ssgets
17c     ( ISHIFT, WHICH, KEV, NP, RITZ, BOUNDS, SHIFTS )
18c
19c\Arguments
20c  ISHIFT  Integer.  (INPUT)
21c          Method for selecting the implicit shifts at each iteration.
22c          ISHIFT = 0: user specified shifts
23c          ISHIFT = 1: exact shift with respect to the matrix H.
24c
25c  WHICH   Character*2.  (INPUT)
26c          Shift selection criteria.
27c          'LM' -> KEV eigenvalues of largest magnitude are retained.
28c          'SM' -> KEV eigenvalues of smallest magnitude are retained.
29c          'LA' -> KEV eigenvalues of largest value are retained.
30c          'SA' -> KEV eigenvalues of smallest value are retained.
31c          'BE' -> KEV eigenvalues, half from each end of the spectrum.
32c                  If KEV is odd, compute one more from the high end.
33c
34c  KEV      Integer.  (INPUT)
35c          KEV+NP is the size of the matrix H.
36c
37c  NP      Integer.  (INPUT)
38c          Number of implicit shifts to be computed.
39c
40c  RITZ    Real array of length KEV+NP.  (INPUT/OUTPUT)
41c          On INPUT, RITZ contains the eigenvalues of H.
42c          On OUTPUT, RITZ are sorted so that the unwanted eigenvalues
43c          are in the first NP locations and the wanted part is in
44c          the last KEV locations.  When exact shifts are selected, the
45c          unwanted part corresponds to the shifts to be applied.
46c
47c  BOUNDS  Real array of length KEV+NP.  (INPUT/OUTPUT)
48c          Error bounds corresponding to the ordering in RITZ.
49c
50c  SHIFTS  Real array of length NP.  (INPUT/OUTPUT)
51c          On INPUT:  contains the user specified shifts if ISHIFT = 0.
52c          On OUTPUT: contains the shifts sorted into decreasing order
53c          of magnitude with respect to the Ritz estimates contained in
54c          BOUNDS. If ISHIFT = 0, SHIFTS is not modified on exit.
55c
56c\EndDoc
57c
58c-----------------------------------------------------------------------
59c
60c\BeginLib
61c
62c\Local variables:
63c     xxxxxx  real
64c
65c\Routines called:
66c     ssortr  ARPACK utility sorting routine.
67c     ivout   ARPACK utility routine that prints integers.
68c     arscnd  ARPACK utility routine for timing.
69c     svout   ARPACK utility routine that prints vectors.
70c     scopy   Level 1 BLAS that copies one vector to another.
71c     sswap   Level 1 BLAS that swaps the contents of two vectors.
72c
73c\Author
74c     Danny Sorensen               Phuong Vu
75c     Richard Lehoucq              CRPC / Rice University
76c     Dept. of Computational &     Houston, Texas
77c     Applied Mathematics
78c     Rice University
79c     Houston, Texas
80c
81c\Revision history:
82c     xx/xx/93: Version ' 2.1'
83c
84c\SCCS Information: @(#)
85c FILE: sgets.F   SID: 2.4   DATE OF SID: 4/19/96   RELEASE: 2
86c
87c\Remarks
88c
89c\EndLib
90c
91c-----------------------------------------------------------------------
92c
93      subroutine ssgets ( ishift, which, kev, np, ritz, bounds, shifts )
94c
95c     %----------------------------------------------------%
96c     | Include files for debugging and timing information |
97c     %----------------------------------------------------%
98c
99      include   'debug.h'
100      include   'stat.h'
101c
102c     %------------------%
103c     | Scalar Arguments |
104c     %------------------%
105c
106      character*2 which
107      integer    ishift, kev, np
108c
109c     %-----------------%
110c     | Array Arguments |
111c     %-----------------%
112c
113      Real
114     &           bounds(kev+np), ritz(kev+np), shifts(np)
115c
116c     %------------%
117c     | Parameters |
118c     %------------%
119c
120      Real
121     &           one, zero
122      parameter (one = 1.0E+0, zero = 0.0E+0)
123c
124c     %---------------%
125c     | Local Scalars |
126c     %---------------%
127c
128      integer    kevd2, msglvl
129c
130c     %----------------------%
131c     | External Subroutines |
132c     %----------------------%
133c
134      external   sswap, scopy, ssortr, arscnd
135c
136c     %---------------------%
137c     | Intrinsic Functions |
138c     %---------------------%
139c
140      intrinsic    max, min
141c
142c     %-----------------------%
143c     | Executable Statements |
144c     %-----------------------%
145c
146c     %-------------------------------%
147c     | Initialize timing statistics  |
148c     | & message level for debugging |
149c     %-------------------------------%
150c
151      call arscnd (t0)
152      msglvl = msgets
153c
154      if (which .eq. 'BE') then
155c
156c        %-----------------------------------------------------%
157c        | Both ends of the spectrum are requested.            |
158c        | Sort the eigenvalues into algebraically increasing  |
159c        | order first then swap high end of the spectrum next |
160c        | to low end in appropriate locations.                |
161c        | NOTE: when np < floor(kev/2) be careful not to swap |
162c        | overlapping locations.                              |
163c        %-----------------------------------------------------%
164c
165         call ssortr ('LA', .true., kev+np, ritz, bounds)
166         kevd2 = kev / 2
167         if ( kev .gt. 1 ) then
168            call sswap ( min(kevd2,np), ritz, 1,
169     &                   ritz( max(kevd2,np)+1 ), 1)
170            call sswap ( min(kevd2,np), bounds, 1,
171     &                   bounds( max(kevd2,np)+1 ), 1)
172         end if
173c
174      else
175c
176c        %----------------------------------------------------%
177c        | LM, SM, LA, SA case.                               |
178c        | Sort the eigenvalues of H into the desired order   |
179c        | and apply the resulting order to BOUNDS.           |
180c        | The eigenvalues are sorted so that the wanted part |
181c        | are always in the last KEV locations.               |
182c        %----------------------------------------------------%
183c
184         call ssortr (which, .true., kev+np, ritz, bounds)
185      end if
186c
187      if (ishift .eq. 1 .and. np .gt. 0) then
188c
189c        %-------------------------------------------------------%
190c        | Sort the unwanted Ritz values used as shifts so that  |
191c        | the ones with largest Ritz estimates are first.       |
192c        | This will tend to minimize the effects of the         |
193c        | forward instability of the iteration when the shifts  |
194c        | are applied in subroutine ssapps.                     |
195c        %-------------------------------------------------------%
196c
197         call ssortr ('SM', .true., np, bounds, ritz)
198         call scopy (np, ritz, 1, shifts, 1)
199      end if
200c
201      call arscnd (t1)
202      tsgets = tsgets + (t1 - t0)
203c
204      if (msglvl .gt. 0) then
205         call ivout (logfil, 1, [kev], ndigit, '_sgets: KEV is')
206         call ivout (logfil, 1, [np], ndigit, '_sgets: NP is')
207         call svout (logfil, kev+np, ritz, ndigit,
208     &        '_sgets: Eigenvalues of current H matrix')
209         call svout (logfil, kev+np, bounds, ndigit,
210     &        '_sgets: Associated Ritz estimates')
211      end if
212c
213      return
214c
215c     %---------------%
216c     | End of ssgets |
217c     %---------------%
218c
219      end
220