"Fossies" - the Fresh Open Source Software Archive

Member "getdp-3.1.0-source/contrib/Arpack/dsortc.f" (31 Jul 2018, 9254 Bytes) of package /linux/privat/getdp-3.1.0-source.tgz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) Fortran 77 source code syntax highlighting (style: standard) with prefixed line numbers. Alternatively you can here view or download the uninterpreted source code file. For more information about "dsortc.f" see the Fossies "Dox" file reference documentation.

    1 c-----------------------------------------------------------------------
    2 c\BeginDoc
    3 c
    4 c\Name: dsortc
    5 c
    6 c\Description:
    7 c  Sorts the complex array in XREAL and XIMAG into the order 
    8 c  specified by WHICH and optionally applies the permutation to the
    9 c  real array Y. It is assumed that if an element of XIMAG is
   10 c  nonzero, then its negative is also an element. In other words,
   11 c  both members of a complex conjugate pair are to be sorted and the
   12 c  pairs are kept adjacent to each other.
   13 c
   14 c\Usage:
   15 c  call dsortc
   16 c     ( WHICH, APPLY, N, XREAL, XIMAG, Y )
   17 c
   18 c\Arguments
   19 c  WHICH   Character*2.  (Input)
   20 c          'LM' -> sort XREAL,XIMAG into increasing order of magnitude.
   21 c          'SM' -> sort XREAL,XIMAG into decreasing order of magnitude.
   22 c          'LR' -> sort XREAL into increasing order of algebraic.
   23 c          'SR' -> sort XREAL into decreasing order of algebraic.
   24 c          'LI' -> sort XIMAG into increasing order of magnitude.
   25 c          'SI' -> sort XIMAG into decreasing order of magnitude.
   26 c          NOTE: If an element of XIMAG is non-zero, then its negative
   27 c                is also an element.
   28 c
   29 c  APPLY   Logical.  (Input)
   30 c          APPLY = .TRUE.  -> apply the sorted order to array Y.
   31 c          APPLY = .FALSE. -> do not apply the sorted order to array Y.
   32 c
   33 c  N       Integer.  (INPUT)
   34 c          Size of the arrays.
   35 c
   36 c  XREAL,  Double precision array of length N.  (INPUT/OUTPUT)
   37 c  XIMAG   Real and imaginary part of the array to be sorted.
   38 c
   39 c  Y       Double precision array of length N.  (INPUT/OUTPUT)
   40 c
   41 c\EndDoc
   42 c
   43 c-----------------------------------------------------------------------
   44 c
   45 c\BeginLib
   46 c
   47 c\Author
   48 c     Danny Sorensen               Phuong Vu
   49 c     Richard Lehoucq              CRPC / Rice University
   50 c     Dept. of Computational &     Houston, Texas
   51 c     Applied Mathematics
   52 c     Rice University           
   53 c     Houston, Texas            
   54 c
   55 c\Revision history:
   56 c     xx/xx/92: Version ' 2.1'
   57 c               Adapted from the sort routine in LANSO.
   58 c
   59 c\SCCS Information: @(#) 
   60 c FILE: sortc.F   SID: 2.3   DATE OF SID: 4/20/96   RELEASE: 2
   61 c
   62 c\EndLib
   63 c
   64 c-----------------------------------------------------------------------
   65 c
   66       subroutine dsortc (which, apply, n, xreal, ximag, y)
   67 c
   68 c     %------------------%
   69 c     | Scalar Arguments |
   70 c     %------------------%
   71 c
   72       character*2 which
   73       logical    apply
   74       integer    n
   75 c
   76 c     %-----------------%
   77 c     | Array Arguments |
   78 c     %-----------------%
   79 c
   80       Double precision     
   81      &           xreal(0:n-1), ximag(0:n-1), y(0:n-1)
   82 c
   83 c     %---------------%
   84 c     | Local Scalars |
   85 c     %---------------%
   86 c
   87       integer    i, igap, j
   88       Double precision     
   89      &           temp, temp1, temp2
   90 c
   91 c     %--------------------%
   92 c     | External Functions |
   93 c     %--------------------%
   94 c
   95       Double precision     
   96      &           dlapy2
   97       external   dlapy2
   98 c
   99 c     %-----------------------%
  100 c     | Executable Statements |
  101 c     %-----------------------%
  102 c
  103       igap = n / 2
  104 c 
  105       if (which .eq. 'LM') then
  106 c
  107 c        %------------------------------------------------------%
  108 c        | Sort XREAL,XIMAG into increasing order of magnitude. |
  109 c        %------------------------------------------------------%
  110 c
  111    10    continue
  112          if (igap .eq. 0) go to 9000
  113 c
  114          do 30 i = igap, n-1
  115             j = i-igap
  116    20       continue
  117 c
  118             if (j.lt.0) go to 30
  119 c
  120             temp1 = dlapy2(xreal(j),ximag(j))
  121             temp2 = dlapy2(xreal(j+igap),ximag(j+igap))
  122 c
  123             if (temp1.gt.temp2) then
  124                 temp = xreal(j)
  125                 xreal(j) = xreal(j+igap)
  126                 xreal(j+igap) = temp
  127 c
  128                 temp = ximag(j)
  129                 ximag(j) = ximag(j+igap)
  130                 ximag(j+igap) = temp
  131 c
  132                 if (apply) then
  133                     temp = y(j)
  134                     y(j) = y(j+igap)
  135                     y(j+igap) = temp
  136                 end if
  137             else
  138                 go to 30
  139             end if
  140             j = j-igap
  141             go to 20
  142    30    continue
  143          igap = igap / 2
  144          go to 10
  145 c
  146       else if (which .eq. 'SM') then
  147 c
  148 c        %------------------------------------------------------%
  149 c        | Sort XREAL,XIMAG into decreasing order of magnitude. |
  150 c        %------------------------------------------------------%
  151 c
  152    40    continue
  153          if (igap .eq. 0) go to 9000
  154 c
  155          do 60 i = igap, n-1
  156             j = i-igap
  157    50       continue
  158 c
  159             if (j .lt. 0) go to 60
  160 c
  161             temp1 = dlapy2(xreal(j),ximag(j))
  162             temp2 = dlapy2(xreal(j+igap),ximag(j+igap))
  163 c
  164             if (temp1.lt.temp2) then
  165                temp = xreal(j)
  166                xreal(j) = xreal(j+igap)
  167                xreal(j+igap) = temp
  168 c
  169                temp = ximag(j)
  170                ximag(j) = ximag(j+igap)
  171                ximag(j+igap) = temp
  172 c 
  173                if (apply) then
  174                   temp = y(j)
  175                   y(j) = y(j+igap)
  176                   y(j+igap) = temp
  177                end if
  178             else
  179                go to 60
  180             endif
  181             j = j-igap
  182             go to 50
  183    60    continue
  184          igap = igap / 2
  185          go to 40
  186 c 
  187       else if (which .eq. 'LR') then
  188 c
  189 c        %------------------------------------------------%
  190 c        | Sort XREAL into increasing order of algebraic. |
  191 c        %------------------------------------------------%
  192 c
  193    70    continue
  194          if (igap .eq. 0) go to 9000
  195 c
  196          do 90 i = igap, n-1
  197             j = i-igap
  198    80       continue
  199 c
  200             if (j.lt.0) go to 90
  201 c
  202             if (xreal(j).gt.xreal(j+igap)) then
  203                temp = xreal(j)
  204                xreal(j) = xreal(j+igap)
  205                xreal(j+igap) = temp
  206 c
  207                temp = ximag(j)
  208                ximag(j) = ximag(j+igap)
  209                ximag(j+igap) = temp
  210 c 
  211                if (apply) then
  212                   temp = y(j)
  213                   y(j) = y(j+igap)
  214                   y(j+igap) = temp
  215                end if
  216             else
  217                go to 90
  218             endif
  219             j = j-igap
  220             go to 80
  221    90    continue
  222          igap = igap / 2
  223          go to 70
  224 c 
  225       else if (which .eq. 'SR') then
  226 c
  227 c        %------------------------------------------------%
  228 c        | Sort XREAL into decreasing order of algebraic. |
  229 c        %------------------------------------------------%
  230 c
  231   100    continue
  232          if (igap .eq. 0) go to 9000
  233          do 120 i = igap, n-1
  234             j = i-igap
  235   110       continue
  236 c
  237             if (j.lt.0) go to 120
  238 c
  239             if (xreal(j).lt.xreal(j+igap)) then
  240                temp = xreal(j)
  241                xreal(j) = xreal(j+igap)
  242                xreal(j+igap) = temp
  243 c
  244                temp = ximag(j)
  245                ximag(j) = ximag(j+igap)
  246                ximag(j+igap) = temp
  247 c 
  248                if (apply) then
  249                   temp = y(j)
  250                   y(j) = y(j+igap)
  251                   y(j+igap) = temp
  252                end if
  253             else
  254                go to 120
  255             endif
  256             j = j-igap
  257             go to 110
  258   120    continue
  259          igap = igap / 2
  260          go to 100
  261 c 
  262       else if (which .eq. 'LI') then
  263 c
  264 c        %------------------------------------------------%
  265 c        | Sort XIMAG into increasing order of magnitude. |
  266 c        %------------------------------------------------%
  267 c
  268   130    continue
  269          if (igap .eq. 0) go to 9000
  270          do 150 i = igap, n-1
  271             j = i-igap
  272   140       continue
  273 c
  274             if (j.lt.0) go to 150
  275 c
  276             if (abs(ximag(j)).gt.abs(ximag(j+igap))) then
  277                temp = xreal(j)
  278                xreal(j) = xreal(j+igap)
  279                xreal(j+igap) = temp
  280 c
  281                temp = ximag(j)
  282                ximag(j) = ximag(j+igap)
  283                ximag(j+igap) = temp
  284 c 
  285                if (apply) then
  286                   temp = y(j)
  287                   y(j) = y(j+igap)
  288                   y(j+igap) = temp
  289                end if
  290             else
  291                go to 150
  292             endif
  293             j = j-igap
  294             go to 140
  295   150    continue
  296          igap = igap / 2
  297          go to 130
  298 c 
  299       else if (which .eq. 'SI') then
  300 c
  301 c        %------------------------------------------------%
  302 c        | Sort XIMAG into decreasing order of magnitude. |
  303 c        %------------------------------------------------%
  304 c
  305   160    continue
  306          if (igap .eq. 0) go to 9000
  307          do 180 i = igap, n-1
  308             j = i-igap
  309   170       continue
  310 c
  311             if (j.lt.0) go to 180
  312 c
  313             if (abs(ximag(j)).lt.abs(ximag(j+igap))) then
  314                temp = xreal(j)
  315                xreal(j) = xreal(j+igap)
  316                xreal(j+igap) = temp
  317 c
  318                temp = ximag(j)
  319                ximag(j) = ximag(j+igap)
  320                ximag(j+igap) = temp
  321 c 
  322                if (apply) then
  323                   temp = y(j)
  324                   y(j) = y(j+igap)
  325                   y(j+igap) = temp
  326                end if
  327             else
  328                go to 180
  329             endif
  330             j = j-igap
  331             go to 170
  332   180    continue
  333          igap = igap / 2
  334          go to 160
  335       end if
  336 c 
  337  9000 continue
  338       return
  339 c
  340 c     %---------------%
  341 c     | End of dsortc |
  342 c     %---------------%
  343 c
  344       end