~ubuntu-branches/ubuntu/saucy/python-scipy/saucy

« back to all changes in this revision

Viewing changes to Lib/sandbox/arpack/ARPACK/SRC/dngets.f

  • Committer: Bazaar Package Importer
  • Author(s): Ondrej Certik
  • Date: 2008-06-16 22:58:01 UTC
  • mfrom: (2.1.24 intrepid)
  • Revision ID: james.westby@ubuntu.com-20080616225801-irdhrpcwiocfbcmt
Tags: 0.6.0-12
* The description updated to match the current SciPy (Closes: #489149).
* Standards-Version bumped to 3.8.0 (no action needed)
* Build-Depends: netcdf-dev changed to libnetcdf-dev

Show diffs side-by-side

added added

removed removed

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