1c-----------------------------------------------------------------------
2c\BeginDoc
3c
4c\Name: pdngets
5c
6c Message Passing Layer: MPI
7c
8c\Description:
9c  Given the eigenvalues of the upper Hessenberg matrix H,
10c  computes the NP shifts AMU that are zeros of the polynomial of
11c  degree NP which filters out components of the unwanted eigenvectors
12c  corresponding to the AMU's based on some given criteria.
13c
14c  NOTE: call this even in the case of user specified shifts in order
15c  to sort the eigenvalues, and error bounds of H for later use.
16c
17c\Usage:
18c  call pdngets
19c     ( COMM, ISHIFT, WHICH, KEV, NP, RITZR, RITZI, BOUNDS, SHIFTR, SHIFTI )
20c
21c\Arguments
22c  COMM    MPI Communicator for the processor grid.  (INPUT)
23c
24c  ISHIFT  Integer.  (INPUT)
25c          Method for selecting the implicit shifts at each iteration.
26c          ISHIFT = 0: user specified shifts
27c          ISHIFT = 1: exact shift with respect to the matrix H.
28c
29c  WHICH   Character*2.  (INPUT)
30c          Shift selection criteria.
31c          'LM' -> want the KEV eigenvalues of largest magnitude.
32c          'SM' -> want the KEV eigenvalues of smallest magnitude.
33c          'LR' -> want the KEV eigenvalues of largest real part.
34c          'SR' -> want the KEV eigenvalues of smallest real part.
35c          'LI' -> want the KEV eigenvalues of largest imaginary part.
36c          'SI' -> want the KEV eigenvalues of smallest imaginary part.
37c
38c  KEV      Integer.  (INPUT/OUTPUT)
39c           INPUT: KEV+NP is the size of the matrix H.
40c           OUTPUT: Possibly increases KEV by one to keep complex conjugate
41c           pairs together.
42c
43c  NP       Integer.  (INPUT/OUTPUT)
44c           Number of implicit shifts to be computed.
45c           OUTPUT: Possibly decreases NP by one to keep complex conjugate
46c           pairs together.
47c
48c  RITZR,  Double precision array of length KEV+NP.  (INPUT/OUTPUT)
49c  RITZI   On INPUT, RITZR and RITZI contain the real and imaginary
50c          parts of the eigenvalues of H.
51c          On OUTPUT, RITZR and RITZI are sorted so that the unwanted
52c          eigenvalues are in the first NP locations and the wanted
53c          portion is in the last KEV locations.  When exact shifts are
54c          selected, the unwanted part corresponds to the shifts to
55c          be applied. Also, if ISHIFT .eq. 1, the unwanted eigenvalues
56c          are further sorted so that the ones with largest Ritz values
57c          are first.
58c
59c  BOUNDS  Double precision array of length KEV+NP.  (INPUT/OUTPUT)
60c          Error bounds corresponding to the ordering in RITZ.
61c
62c  SHIFTR, SHIFTI  *** USE deprecated as of version 2.1. ***
63c
64c
65c\EndDoc
66c
67c-----------------------------------------------------------------------
68c
69c\BeginLib
70c
71c\Local variables:
72c     xxxxxx  real
73c
74c\Routines called:
75c     dsortc  ARPACK sorting routine.
76c     dcopy   Level 1 BLAS that copies one vector to another .
77c
78c\Author
79c     Danny Sorensen               Phuong Vu
80c     Richard Lehoucq              CRPC / Rice University
81c     Dept. of Computational &     Houston, Texas
82c     Applied Mathematics
83c     Rice University
84c     Houston, Texas
85c
86c\Parallel Modifications
87c     Kristi Maschhoff
88c
89c\Revision history:
90c     Starting Point: Serial Code FILE: ngets.F   SID: 2.2
91c
92c\SCCS Information:
93c FILE: ngets.F   SID: 1.2   DATE OF SID: 2/22/96
94c
95c\Remarks
96c     1. xxxx
97c
98c\EndLib
99c
100c-----------------------------------------------------------------------
101c
102      subroutine pdngets
103     &                 ( comm, ishift, which, kev, np, ritzr, ritzi,
104     &                   bounds, shiftr, shifti )
105c
106c     %--------------------%
107c     | MPI Communicator |
108c     %--------------------%
109c
110      integer   comm
111c
112c     %----------------------------------------------------%
113c     | Include files for debugging and timing information |
114c     %----------------------------------------------------%
115c
116      include   'debug.h'
117      include   'stat.h'
118c
119c     %------------------%
120c     | Scalar Arguments |
121c     %------------------%
122c
123      character*2 which
124      integer    ishift, kev, np
125c
126c     %-----------------%
127c     | Array Arguments |
128c     %-----------------%
129c
130      Double precision
131     &           bounds(kev+np), ritzr(kev+np), ritzi(kev+np),
132     &           shiftr(1), shifti(1)
133c
134c     %------------%
135c     | Parameters |
136c     %------------%
137c
138      Double precision
139     &           one, zero
140      parameter (one = 1.0, zero = 0.0)
141c
142c     %---------------%
143c     | Local Scalars |
144c     %---------------%
145c
146      integer    msglvl
147c
148c     %----------------------%
149c     | External Subroutines |
150c     %----------------------%
151c
152      external   dcopy, dsortc, arscnd
153c
154c     %----------------------%
155c     | Intrinsics Functions |
156c     %----------------------%
157c
158      intrinsic  abs
159c
160c     %-----------------------%
161c     | Executable Statements |
162c     %-----------------------%
163c
164c     %-------------------------------%
165c     | Initialize timing statistics  |
166c     | & message level for debugging |
167c     %-------------------------------%
168c
169      call arscnd (t0)
170      msglvl = mngets
171c
172c     %----------------------------------------------------%
173c     | LM, SM, LR, SR, LI, SI case.                       |
174c     | Sort the eigenvalues of H into the desired order   |
175c     | and apply the resulting order to BOUNDS.           |
176c     | The eigenvalues are sorted so that the wanted part |
177c     | are always in the last KEV locations.              |
178c     | We first do a pre-processing sort in order to keep |
179c     | complex conjugate pairs together                   |
180c     %----------------------------------------------------%
181c
182      if (which .eq. 'LM') then
183         call dsortc ('LR', .true., kev+np, ritzr, ritzi, bounds)
184      else if (which .eq. 'SM') then
185         call dsortc ('SR', .true., kev+np, ritzr, ritzi, bounds)
186      else if (which .eq. 'LR') then
187         call dsortc ('LM', .true., kev+np, ritzr, ritzi, bounds)
188      else if (which .eq. 'SR') then
189         call dsortc ('SM', .true., kev+np, ritzr, ritzi, bounds)
190      else if (which .eq. 'LI') then
191         call dsortc ('LM', .true., kev+np, ritzr, ritzi, bounds)
192      else if (which .eq. 'SI') then
193         call dsortc ('SM', .true., kev+np, ritzr, ritzi, bounds)
194      end if
195c
196      call dsortc (which, .true., kev+np, ritzr, ritzi, bounds)
197c
198c     %-------------------------------------------------------%
199c     | Increase KEV by one if the ( ritzr(np),ritzi(np) )    |
200c     | = ( ritzr(np+1),-ritzi(np+1) ) and ritz(np) .ne. zero |
201c     | Accordingly decrease NP by one. In other words keep   |
202c     | complex conjugate pairs together.                     |
203c     %-------------------------------------------------------%
204c
205      if (       ( ritzr(np+1) - ritzr(np) ) .eq. zero
206     &     .and. ( ritzi(np+1) + ritzi(np) ) .eq. zero ) then
207         np = np - 1
208         kev = kev + 1
209      end if
210c
211      if ( ishift .eq. 1 ) then
212c
213c        %-------------------------------------------------------%
214c        | Sort the unwanted Ritz values used as shifts so that  |
215c        | the ones with largest Ritz estimates are first        |
216c        | This will tend to minimize the effects of the         |
217c        | forward instability of the iteration when they shifts |
218c        | are applied in subroutine pdnapps.                    |
219c        | Be careful and use 'SR' since we want to sort BOUNDS! |
220c        %-------------------------------------------------------%
221c
222         call dsortc ( 'SR', .true., np, bounds, ritzr, ritzi )
223      end if
224c
225      call arscnd (t1)
226      tngets = tngets + (t1 - t0)
227c
228      if (msglvl .gt. 0) then
229         call pivout (comm, logfil, 1, [kev], ndigit, '_ngets: KEV is')
230         call pivout (comm, logfil, 1, [np], ndigit, '_ngets: NP is')
231         call pdvout (comm, logfil, kev+np, ritzr, ndigit,
232     &        '_ngets: Eigenvalues of current H matrix -- real part')
233         call pdvout (comm, logfil, kev+np, ritzi, ndigit,
234     &        '_ngets: Eigenvalues of current H matrix -- imag part')
235         call pdvout (comm, logfil, kev+np, bounds, ndigit,
236     &      '_ngets: Ritz estimates of the current KEV+NP Ritz values')
237      end if
238c
239      return
240c
241c     %----------------%
242c     | End of pdngets |
243c     %----------------%
244c
245      end
246