Mercurial > octave-nkf
diff libcruft/arpack/src/csortc.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/csortc.f Fri Jan 28 14:04:33 2011 -0500 @@ -0,0 +1,322 @@ +c\BeginDoc +c +c\Name: csortc +c +c\Description: +c Sorts the Complex array in X into the order +c specified by WHICH and optionally applies the permutation to the +c Real array Y. +c +c\Usage: +c call csortc +c ( WHICH, APPLY, N, X, Y ) +c +c\Arguments +c WHICH Character*2. (Input) +c 'LM' -> sort X into increasing order of magnitude. +c 'SM' -> sort X into decreasing order of magnitude. +c 'LR' -> sort X with real(X) in increasing algebraic order +c 'SR' -> sort X with real(X) in decreasing algebraic order +c 'LI' -> sort X with imag(X) in increasing algebraic order +c 'SI' -> sort X with imag(X) in decreasing algebraic order +c +c APPLY Logical. (Input) +c APPLY = .TRUE. -> apply the sorted order to array Y. +c APPLY = .FALSE. -> do not apply the sorted order to array Y. +c +c N Integer. (INPUT) +c Size of the arrays. +c +c X Complex array of length N. (INPUT/OUTPUT) +c This is the array to be sorted. +c +c Y Complex array of length N. (INPUT/OUTPUT) +c +c\EndDoc +c +c----------------------------------------------------------------------- +c +c\BeginLib +c +c\Routines called: +c slapy2 LAPACK routine to compute sqrt(x**2+y**2) carefully. +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 Adapted from the sort routine in LANSO. +c +c\SCCS Information: @(#) +c FILE: sortc.F SID: 2.2 DATE OF SID: 4/20/96 RELEASE: 2 +c +c\EndLib +c +c----------------------------------------------------------------------- +c + subroutine csortc (which, apply, n, x, y) +c +c %------------------% +c | Scalar Arguments | +c %------------------% +c + character*2 which + logical apply + integer n +c +c %-----------------% +c | Array Arguments | +c %-----------------% +c + Complex + & x(0:n-1), y(0:n-1) +c +c %---------------% +c | Local Scalars | +c %---------------% +c + integer i, igap, j + Complex + & temp + Real + & temp1, temp2 +c +c %--------------------% +c | External functions | +c %--------------------% +c + Real + & slapy2 +c +c %--------------------% +c | Intrinsic Functions | +c %--------------------% + Intrinsic + & real, aimag +c +c %-----------------------% +c | Executable Statements | +c %-----------------------% +c + igap = n / 2 +c + if (which .eq. 'LM') then +c +c %--------------------------------------------% +c | Sort X into increasing order of magnitude. | +c %--------------------------------------------% +c + 10 continue + if (igap .eq. 0) go to 9000 +c + do 30 i = igap, n-1 + j = i-igap + 20 continue +c + if (j.lt.0) go to 30 +c + temp1 = slapy2(real(x(j)),aimag(x(j))) + temp2 = slapy2(real(x(j+igap)),aimag(x(j+igap))) +c + if (temp1.gt.temp2) then + temp = x(j) + x(j) = x(j+igap) + x(j+igap) = temp +c + if (apply) then + temp = y(j) + y(j) = y(j+igap) + y(j+igap) = temp + end if + else + go to 30 + end if + j = j-igap + go to 20 + 30 continue + igap = igap / 2 + go to 10 +c + else if (which .eq. 'SM') then +c +c %--------------------------------------------% +c | Sort X into decreasing order of magnitude. | +c %--------------------------------------------% +c + 40 continue + if (igap .eq. 0) go to 9000 +c + do 60 i = igap, n-1 + j = i-igap + 50 continue +c + if (j .lt. 0) go to 60 +c + temp1 = slapy2(real(x(j)),aimag(x(j))) + temp2 = slapy2(real(x(j+igap)),aimag(x(j+igap))) +c + if (temp1.lt.temp2) then + temp = x(j) + x(j) = x(j+igap) + x(j+igap) = temp +c + if (apply) then + temp = y(j) + y(j) = y(j+igap) + y(j+igap) = temp + end if + else + go to 60 + endif + j = j-igap + go to 50 + 60 continue + igap = igap / 2 + go to 40 +c + else if (which .eq. 'LR') then +c +c %------------------------------------------------% +c | Sort XREAL into increasing order of algebraic. | +c %------------------------------------------------% +c + 70 continue + if (igap .eq. 0) go to 9000 +c + do 90 i = igap, n-1 + j = i-igap + 80 continue +c + if (j.lt.0) go to 90 +c + if (real(x(j)).gt.real(x(j+igap))) then + temp = x(j) + x(j) = x(j+igap) + x(j+igap) = temp +c + if (apply) then + temp = y(j) + y(j) = y(j+igap) + y(j+igap) = temp + end if + else + go to 90 + endif + j = j-igap + go to 80 + 90 continue + igap = igap / 2 + go to 70 +c + else if (which .eq. 'SR') then +c +c %------------------------------------------------% +c | Sort XREAL into decreasing order of algebraic. | +c %------------------------------------------------% +c + 100 continue + if (igap .eq. 0) go to 9000 + do 120 i = igap, n-1 + j = i-igap + 110 continue +c + if (j.lt.0) go to 120 +c + if (real(x(j)).lt.real(x(j+igap))) then + temp = x(j) + x(j) = x(j+igap) + x(j+igap) = temp +c + if (apply) then + temp = y(j) + y(j) = y(j+igap) + y(j+igap) = temp + end if + else + go to 120 + endif + j = j-igap + go to 110 + 120 continue + igap = igap / 2 + go to 100 +c + else if (which .eq. 'LI') then +c +c %--------------------------------------------% +c | Sort XIMAG into increasing algebraic order | +c %--------------------------------------------% +c + 130 continue + if (igap .eq. 0) go to 9000 + do 150 i = igap, n-1 + j = i-igap + 140 continue +c + if (j.lt.0) go to 150 +c + if (aimag(x(j)).gt.aimag(x(j+igap))) then + temp = x(j) + x(j) = x(j+igap) + x(j+igap) = temp +c + if (apply) then + temp = y(j) + y(j) = y(j+igap) + y(j+igap) = temp + end if + else + go to 150 + endif + j = j-igap + go to 140 + 150 continue + igap = igap / 2 + go to 130 +c + else if (which .eq. 'SI') then +c +c %---------------------------------------------% +c | Sort XIMAG into decreasing algebraic order | +c %---------------------------------------------% +c + 160 continue + if (igap .eq. 0) go to 9000 + do 180 i = igap, n-1 + j = i-igap + 170 continue +c + if (j.lt.0) go to 180 +c + if (aimag(x(j)).lt.aimag(x(j+igap))) then + temp = x(j) + x(j) = x(j+igap) + x(j+igap) = temp +c + if (apply) then + temp = y(j) + y(j) = y(j+igap) + y(j+igap) = temp + end if + else + go to 180 + endif + j = j-igap + go to 170 + 180 continue + igap = igap / 2 + go to 160 + end if +c + 9000 continue + return +c +c %---------------% +c | End of csortc | +c %---------------% +c + end