Mercurial > octave
diff libcruft/arpack/src/ssgets.f @ 12274:9f5d2ef078e8 release-3-4-x
import ARPACK sources to libcruft from Debian package libarpack2 2.1+parpack96.dfsg-3+b1
author | John W. Eaton <jwe@octave.org> |
---|---|
date | Fri, 28 Jan 2011 14:04:33 -0500 |
parents | |
children |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/arpack/src/ssgets.f Fri Jan 28 14:04:33 2011 -0500 @@ -0,0 +1,219 @@ +c----------------------------------------------------------------------- +c\BeginDoc +c +c\Name: ssgets +c +c\Description: +c Given the eigenvalues of the symmetric tridiagonal matrix H, +c computes the NP shifts AMU that are zeros of the polynomial of +c degree NP which filters out components of the unwanted eigenvectors +c corresponding to the AMU's based on some given criteria. +c +c NOTE: This is called even in the case of user specified shifts in +c order to sort the eigenvalues, and error bounds of H for later use. +c +c\Usage: +c call ssgets +c ( ISHIFT, WHICH, KEV, NP, RITZ, BOUNDS, SHIFTS ) +c +c\Arguments +c ISHIFT Integer. (INPUT) +c Method for selecting the implicit shifts at each iteration. +c ISHIFT = 0: user specified shifts +c ISHIFT = 1: exact shift with respect to the matrix H. +c +c WHICH Character*2. (INPUT) +c Shift selection criteria. +c 'LM' -> KEV eigenvalues of largest magnitude are retained. +c 'SM' -> KEV eigenvalues of smallest magnitude are retained. +c 'LA' -> KEV eigenvalues of largest value are retained. +c 'SA' -> KEV eigenvalues of smallest value are retained. +c 'BE' -> KEV eigenvalues, half from each end of the spectrum. +c If KEV is odd, compute one more from the high end. +c +c KEV Integer. (INPUT) +c KEV+NP is the size of the matrix H. +c +c NP Integer. (INPUT) +c Number of implicit shifts to be computed. +c +c RITZ Real array of length KEV+NP. (INPUT/OUTPUT) +c On INPUT, RITZ contains the eigenvalues of H. +c On OUTPUT, RITZ are sorted so that the unwanted eigenvalues +c are in the first NP locations and the wanted part is in +c the last KEV locations. When exact shifts are selected, the +c unwanted part corresponds to the shifts to be applied. +c +c BOUNDS Real array of length KEV+NP. (INPUT/OUTPUT) +c Error bounds corresponding to the ordering in RITZ. +c +c SHIFTS Real array of length NP. (INPUT/OUTPUT) +c On INPUT: contains the user specified shifts if ISHIFT = 0. +c On OUTPUT: contains the shifts sorted into decreasing order +c of magnitude with respect to the Ritz estimates contained in +c BOUNDS. If ISHIFT = 0, SHIFTS is not modified on exit. +c +c\EndDoc +c +c----------------------------------------------------------------------- +c +c\BeginLib +c +c\Local variables: +c xxxxxx real +c +c\Routines called: +c ssortr ARPACK utility sorting routine. +c ivout ARPACK utility routine that prints integers. +c arscnd ARPACK utility routine for timing. +c svout ARPACK utility routine that prints vectors. +c scopy Level 1 BLAS that copies one vector to another. +c sswap Level 1 BLAS that swaps the contents of two vectors. +c +c\Author +c Danny Sorensen Phuong Vu +c Richard Lehoucq CRPC / Rice University +c Dept. of Computational & Houston, Texas +c Applied Mathematics +c Rice University +c Houston, Texas +c +c\Revision history: +c xx/xx/93: Version ' 2.1' +c +c\SCCS Information: @(#) +c FILE: sgets.F SID: 2.4 DATE OF SID: 4/19/96 RELEASE: 2 +c +c\Remarks +c +c\EndLib +c +c----------------------------------------------------------------------- +c + subroutine ssgets ( ishift, which, kev, np, ritz, bounds, shifts ) +c +c %----------------------------------------------------% +c | Include files for debugging and timing information | +c %----------------------------------------------------% +c + include 'debug.h' + include 'stat.h' +c +c %------------------% +c | Scalar Arguments | +c %------------------% +c + character*2 which + integer ishift, kev, np +c +c %-----------------% +c | Array Arguments | +c %-----------------% +c + Real + & bounds(kev+np), ritz(kev+np), shifts(np) +c +c %------------% +c | Parameters | +c %------------% +c + Real + & one, zero + parameter (one = 1.0E+0, zero = 0.0E+0) +c +c %---------------% +c | Local Scalars | +c %---------------% +c + integer kevd2, msglvl +c +c %----------------------% +c | External Subroutines | +c %----------------------% +c + external sswap, scopy, ssortr, arscnd +c +c %---------------------% +c | Intrinsic Functions | +c %---------------------% +c + intrinsic max, min +c +c %-----------------------% +c | Executable Statements | +c %-----------------------% +c +c %-------------------------------% +c | Initialize timing statistics | +c | & message level for debugging | +c %-------------------------------% +c + call arscnd (t0) + msglvl = msgets +c + if (which .eq. 'BE') then +c +c %-----------------------------------------------------% +c | Both ends of the spectrum are requested. | +c | Sort the eigenvalues into algebraically increasing | +c | order first then swap high end of the spectrum next | +c | to low end in appropriate locations. | +c | NOTE: when np < floor(kev/2) be careful not to swap | +c | overlapping locations. | +c %-----------------------------------------------------% +c + call ssortr ('LA', .true., kev+np, ritz, bounds) + kevd2 = kev / 2 + if ( kev .gt. 1 ) then + call sswap ( min(kevd2,np), ritz, 1, + & ritz( max(kevd2,np)+1 ), 1) + call sswap ( min(kevd2,np), bounds, 1, + & bounds( max(kevd2,np)+1 ), 1) + end if +c + else +c +c %----------------------------------------------------% +c | LM, SM, LA, SA case. | +c | Sort the eigenvalues of H into the desired order | +c | and apply the resulting order to BOUNDS. | +c | The eigenvalues are sorted so that the wanted part | +c | are always in the last KEV locations. | +c %----------------------------------------------------% +c + call ssortr (which, .true., kev+np, ritz, bounds) + end if +c + if (ishift .eq. 1 .and. np .gt. 0) then +c +c %-------------------------------------------------------% +c | Sort the unwanted Ritz values used as shifts so that | +c | the ones with largest Ritz estimates are first. | +c | This will tend to minimize the effects of the | +c | forward instability of the iteration when the shifts | +c | are applied in subroutine ssapps. | +c %-------------------------------------------------------% +c + call ssortr ('SM', .true., np, bounds, ritz) + call scopy (np, ritz, 1, shifts, 1) + end if +c + call arscnd (t1) + tsgets = tsgets + (t1 - t0) +c + if (msglvl .gt. 0) then + call ivout (logfil, 1, kev, ndigit, '_sgets: KEV is') + call ivout (logfil, 1, np, ndigit, '_sgets: NP is') + call svout (logfil, kev+np, ritz, ndigit, + & '_sgets: Eigenvalues of current H matrix') + call svout (logfil, kev+np, bounds, ndigit, + & '_sgets: Associated Ritz estimates') + end if +c + return +c +c %---------------% +c | End of ssgets | +c %---------------% +c + end