C      ALGORITHM 663, COLLECTED ALGORITHMS FROM ACM.
C      THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE,
C      VOL. 14, NO. 2, P.177.
      SUBROUTINE BLASDOC
C                                                                      C
C -------------------------------------------------------------------- C
C                                                                      C
C            Documentation of vectorized BLAS for the Cyber 205        C
C                                                                      C
C -------------------------------------------------------------------- C
C Language : FORTRAN 200                                               C
C ========                                                             C
C -------------------------------------------------------------------- C
C                                                                      C
C Purpose :                                                            C
C =======                                                              C
C The BLAS, or, Basic Linear Algebra Subprograms, is a set of          C
C subprograms, performing operations on vectors. The BLAS[1]           C
C describes a package of 38 Fortran-callable subprograms. This         C
C package contains only the REAL and COMPLEX subprograms performing    C
C operations on vectors. This means that no DOUBLE PRECISION and no    C
C Givens transformation subprograms have been implemented (see under   C
C "Double Precision" and "Givens transformation").                     C
C According to the LINPACK[2], the subprogram CSROT has been added.    C
C                                                                      C
C For more information on this vectorized BLAS, see [3] and [4].       C
C                                                                      C
C -------------------------------------------------------------------- C
C                                                                      C
C References :                                                         C
C ==========                                                           C
C 1. C.L. Lawson, R.J. Hanson, D.R. Kincaid, and F.T. Krogh (1979).    C
C    Basic Linear Algebra Subprograms for Fortran Usage, ACM Trans-    C
C    actions on Mathematical Software, 5, 308-323.                     C
C 2. J.J. Dongarra, J.R. Bunch, C.B. Moler, and G.W. Stewart (1979).   C
C    LINPACK user's guide, SIAM publications.                          C
C 3. M. Louter-Nool (1987). Basic Linear Algebra Subprograms (BLAS) on C
C    the CDC Cyber 205, Parallel Computing, 4, 143-165.                C
C 4. M. Louter-Nool. Translation of Algorithm 539: Basic Linear        C
C    Algebra Subprograms for Fortran Usage in FORTRAN 200 for the      C
C    Cyber 205, Report NM-R8702, CWI, Amsterdam, 1987. (Submitted      C
C    to ACM TOMS)                                                      C
C 5. CDC Cyber 200 FORTRAN reference manual, version 1,                C
C    publ. number 60480200.                                            C
C 6. CDC Cyber 200 ASSEMBLER reference manual, version 2,              C
C    publ. number 60485010.                                            C
C 7. CDC VSOS reference manual, version 2,                             C
C    publ. number 6049410.                                             C
C                                                                      C
C -------------------------------------------------------------------- C
C                                                                      C
C The increment parameters :                                           C
C ========================                                             C
C The parameters INCX and INCY allow the BLAS to operate on vectors    C
C whose elements are non-contiguous in memory. Both positive and       C
C negative values are permitted. Only positive values of INCX are      C
C allowed for operations that have a single vector argument.           C
C If INCX > 0, then the components x(i), i = 1,...,N, of a vector x    C
C are stored in                                                        C
C                                                                      C
C     X(1 + (i - 1) * INCX).                                           C
C                                                                      C
C If INCX < 0, then x(i) is stored in                                  C
C                                                                      C
C     X(1 - (N - i) * INCX).                                           C
C                                                                      C
C The original BLAS permits an INCX value of 0. Here, this option was  C
C not implemented. An error message will be printed.                   C
C                                                                      C
C -------------------------------------------------------------------- C
C                                                                      C
C Parameters :                                                         C
C ==========                                                           C
C The parameters of BLAS have the following type :                     C
C                                                                      C
C     INTEGER   N, INCX, INCY                                          C
C     REAL      SX(NX), SY(NY), SA, SC, SS                             C
C     COMPLEX   CX(NX), CY(NY), CA                                     C
C                                                                      C
C The dimensions NX of SX (CX), and NY of SY (CY) must satisfy         C
C                                                                      C
C      NX >= 1 + IABS (INCX) * (N - 1).                                C
C      NY >= 1 + IABS (INCY) * (N - 1).                                C
C                                                                      C
C -------------------------------------------------------------------- C
C                                                                      C
C Type declarations for function names :                               C
C ====================================                                 C
C The function names are of the following type :                       C
C                                                                      C
C     INTEGER   ISAMAX, ICAMAX                                         C
C     REAL      SASUM, SCASUM, SNRM2, SCNRM2, SDOT                     C
C     COMPLEX   CDOTC, CDOTU                                           C
C                                                                      C
C -------------------------------------------------------------------- C
C                                                                      C
C The subprograms :                                                    C
C ===============                                                      C
C This package contains the following subprograms :                    C
C                                                                      C
C 1.  IMAX = ISAMAX (N, SX, INCX)                                      C
C 2.  IMAX = ICAMAX (N, CX, INCX)                                      C
C     finding the largest component of a vector                        C
C                                                                      C
C 3.  SW = SNRM2  (N, SX, INCX)                                        C
C 4.  SW = SCNRM2 (N, CX, INCX)                                        C
C     computing the Euclidian length or L2 norm of a vector            C
C                                                                      C
C 5.  SW = SASUM  (N, SX, INCX)                                        C
C 6.  SW = SCASUM (N, CX, INCX)                                        C
C     computing the sum of magnitudes of vector components             C
C                                                                      C
C 7.  SW = SDOT  (N, SX, INCX, SY, INCY)                               C
C 8.  SW = CDOTU (N, SX, INCX, SY, INCY)                               C
C 9.  SW = CDOTC (N, SX, INCX, SY, INCY)                               C
C     computing the dot products                                       C
C                                                                      C
C 10. CALL SAXPY (N, SA, SX, INCX, SY, INCY)                           C
C 11. CALL CAXPY (N, CA, CX, INCX, CY, INCY)                           C
C     performing the elementary vector operation: y := a.x + y         C
C                                                                      C
C 12. CALL SSCAL  (N, SA, SX, INCX)                                    C
C 13. CALL CSCAL  (N, CA, CX, INCX)                                    C
C 14. CALL CSSCAL (N, SA, CX, INCX)                                    C
C     performing the vector scaling: x := a.x                          C
C                                                                      C
C 15. CALL SCOPY (N, SX, INCX, SY, INCY)                               C
C 16. CALL CCOPY (N, CX, INCX, CY, INCY)                               C
C     copying a vector x to a vector y:  y := x                        C
C                                                                      C
C 17. CALL SSWAP (N, SX, INCX, SY, INCY)                               C
C 18. CALL CSWAP (N, CX, INCX, CY, INCY)                               C
C     interchanging the vectors x and y: x :=: y                       C
C                                                                      C
C 19. CALL SROT  (N, SX, INCX, SY, INCY, SC, SS)                       C
C 20. CALL CSROT (N, CX, INCX, CY, INCY, SC, SS)                       C
C     applying a plane rotation                                        C
C                                                                      C
C Note :                                                               C
C ----                                                                 C
C For both SROT and CSROT two implementations are available, namely,   C
C SROT - the original SROT3 from [3, 4] - and SROT4, and CSROT         C
C - the original CSROT3 from [3, 4] - and CSROT4, respectively.        C
C The _ROT implementations perform the operations in 3 vector          C
C operations, and the _ROT4 in 4 vector operations. By omitting the    C
C suffix 4 the _ROT4 routines can be loaded instead of the _ROT3's.    C
C The numerical results may differ slightly. For more information      C
C see [3].
C                                                                      C
C -------------------------------------------------------------------- C
C                                                                      C
C General comments :                                                   C
C ================                                                     C
C For the REAL vectors the value N should not exceed the number of     C
C data elements on a LARGE PAGE (65535). For COMPLEX vectors the       C
C length is restricted to 32767. If larger values of N are given,      C
C an error message will be printed.                                    C
C                                                                      C
C When operating on non-contiguously stored vectors, a great part      C
C of the execution time will be spent on data movements. So, if        C
C a non-contiguous vector is frequently used, it is wise to copy such  C
C a vector to a contiguously stored vector, e.g., by means of a        C
C BLAS _COPY routine.                                                  C
C                                                                      C
C When both vectors would have a negative increment value, all         C
C subprograms act as if both were positive. This implies that          C
C the process is NOT carried out in reverse storage order. This saves  C
C a lot of probably superfluous reverse operations. If one really      C
C wants to operate in reverse order, we recommend the use of a         C
C _COPY operation, like                                                C
C                                                                      C
C     ASSIGN SXD, .DYN.N                                               C
C     ASSIGN SYD, .DYN.N                                               C
C     CALL SCOPY (N, SX, INCX, SXD, 1)                                 C
C     CALL SCOPY (N, SY, INCY, SYD, 1)                                 C
C     ......                                                           C
C     CALL       (N, SXD, 1, SYD, 1)                                   C
C     ......                                                           C
C     CALL SCOPY (N, SXD, INCX, SX, 1)                                 C
C     CALL SCOPY (N, SYD, INCY, SY, 1)                                 C
C                                                                      C
C                                                                      C
C Note: For one positive and one negative increment value the order    C
C       of processing depends on the increment values and the          C
C       operation involved. In all subprograms the data movements      C
C       were kept to a minimum. For more information see [4].          C
C                                                                      C
C -------------------------------------------------------------------- C
C                                                                      C
C Implementation :                                                     C
C ==============                                                       C
C The subprograms are coded in CONTROL DATA FORTRAN 200, suitable for  C
C Cyber 200 vector machines. This language is a superset of the ANSI   C
C FORTRAN 77. In order to obtain an optimal code vector intrinsic      C
C functions and SPECIAL CALLS are used. A special call statement       C
C directly generates machine language instructions.                    C
C                                                                      C
C Vector intrinsic functions used (see [5]) :                          C
C -----------------------------------------                            C
C     i = Q8SMAXI  (v)           find the largest element of vector v  C
C     u = Q8VGATHP (v,i,n;u)     gather n elements of v with period i  C
C                                into vector u                         C
C     u = Q8VMERG  (v1,v2,cv;u)  merge the vectors v1 and v2 into u    C
C                                under control of bit vector cv        C
C     u = Q8VREV   (v;u)         store v in reverse order into u       C
C     u = Q8VSCATP (v,i,n;u)     inverse operation of Q8VGATHP         C
C                                                                      C
C Special calls used (see [6]) :                                       C
C ----------------------------                                         C
C     CALL Q8ADDL  (Rf,Sf,Tf)        add lower, full word              C
C     CALL Q8ADDN  (Rf,Sf,Tf)        add normalized, full word         C
C     CALL Q8ADDNS (G,X,A,Y,B,Z,C)   add normalized, sparse vector     C
C     CALL Q8ADDNV (G,X,A,Y,B,Z,C)   add normalized, vector            C
C     CALL Q8ADDX  (Rf,Sf,Tf)        add index, full word              C
C     CALL Q8DOTV  (G,X,A,Y,B,Z,C)   dot product vector                C
C     CALL Q8MAX   (G,X,A, ,B,Z,C)   vector maximum                    C
C     CALL Q8PACK  (Rf,Sf,Tf)        pack, full word                   C
C     CALL Q8SUM   (G,X,A, , ,Z,C)   vector sum                        C
C     CALL Q8VREV  (G,X,A, , ,Z,C)   transmit vector reversed to       C
C                                    vector                            C
C                                                                      C
C -------------------------------------------------------------------- C
C                                                                      C
C Error messages :                                                     C
C ==============                                                       C
C Error messages will be printed in the following cases :              C
C    1. a) for REAL    vectors when N > 65535                          C
C       b) for COMPLEX vectors when N > 32767                          C
C    2. a) for vectors with one vector argument  when INCX < 1         C
C       b) for vectors with two vector arguments when INCX = 0         C
C                                                                      C
C -------------------------------------------------------------------- C
C                                                                      C
C Storage :                                                            C
C =======                                                              C
C In each subprogram some BIT vector of length 65535 are declared      C
C and initialized by means of a DATA statement.                        C
C All other auxiliary vectors storage is dynamically allocated.        C
C                                                                      C
C -------------------------------------------------------------------- C
C                                                                      C
C Execution time :                                                     C
C ==============                                                       C
C In [4] a complete summary of the expected execution times for all    C
C BLAS subprograms as a function of N is given. In general, the        C
C execution time depends on p, the number of pipes of your Cyber 205.  C
C In case of more pipes the execution time is accordingly reduced.     C
C                                                                      C
C -------------------------------------------------------------------- C
C                                                                      C
C Double precision :                                                   C
C ================                                                     C
C The double precision subprograms have not been vectorized. The       C
C Cyber 205 employs 64-bits full words and 32-bits half words.         C
C The execution time can be increased by a factor of 20 when           C
C operating in double precision instead of single precision. If you    C
C do need double precision results, then use the scalar CDC version    C
C of the BLAS, which is also available via TOMS.                       C
C                                                                      C
C -------------------------------------------------------------------- C
C                                                                      C
C Givens transformation subprograms :                                  C
C =================================                                    C
C All these subprograms don't use vectors. The codes can be obtained   C
C via TOMS.                                                            C
C                                                                      C
C -------------------------------------------------------------------- C
C                                                                      C
C How to use the vectorized BLAS :                                     C
C ==============================                                       C
C The subprograms can be translated with the FORTRAN 200 compiler,     C
C e.g., as follows :                                                   C
C                                                                      C
C      FTN200,OPT,SC,UNS,LO.                                           C
C                                                                      C
C This means :                                                         C
C      The object code is written to the file BINARY/16.               C
C      OPT causes that all optimizations are performed.                C
C      SC interprets names as machine instruction references.          C
C      UNS allows unsafe optimization.                                 C
C      LO delivers the source code and the cross-reference list.       C
C                                                                      C
C For more information about the FTN200 control statement, see [5].    C
C                                                                      C
C Now the binary code can be added to a library at the Cyber 205.      C
C The next control statements are sufficient to create a library       C
C with the BLAS subprograms named OLDLIB.                              C
C                                                                      C
C     OLE,I=BINARY,N=OLDLIB.                                           C
C     DEFINE,OLDLIB.                                                   C
C                                                                      C
C If you want to add the binary code to your existing library, e.g.,   C
C called OLDLIB, then you can use :                                    C
C                                                                      C
C     ATTACH,OLDLIB.                                                   C
C     OLE,I=BINARY,OLDLIB,N=NEWLIB.                                    C
C     DEFINE,NEWLIB.                                                   C
C     PURGE,OLDLIB.                                                    C
C     RETURN,OLDLIB.                                                   C
C     SWITCH,NEWLIB,OLDLIB.                                            C
C                                                                      C
C Your library OLDLIB is replaced by a new library with the same name, C
C including the BLAS subprograms with the same name.                   C
C                                                                      C
C Since FORTRAN 200 is a superset of ANSI FORTRAN 77, any correct      C
C FORTRAN 77 program can be compiled and executed successfully, on     C
C the Cyber 205. Use of the BLAS subprograms can improve the           C
C readability, portability and efficiency of your programs.            C
C When using these vectorized BLAS a great part of the execution       C
C can be performed at vector speed. If N is sufficiently large (N>50)  C
C the execution time will decrease extremely as compared with          C
C scalar execution.                                                    C
C                                                                      C
C An example of the control statements needed for execution on the     C
C Cyber 205.                                                           C
C                                                                      C
C     ATTACH,OLDLIB.                                                   C
C     FTN200,I=PROGRAM.                                                C
C     LOAD,LIB=OLDLIB,GRLPALL= ,L=O.                                   C
C     GO.                                                              C
C                                                                      C
C The file PROGRAM is translated by the FTN200 compiler. The used      C
C subprograms included in the library OLDLIB are loaded, and the       C
C binary code is executed. The parameter GRLPALL indicates that all    C
C modules, data base and common blocks are grouped and mapped on       C
C large page boundary.                                                 C
C                                                                      C
C For more information about the control statements, see [7].          C
C                                                                      C
C -------------------------------------------------------------------- C
C                                                                      C
      RETURN
      END
      FUNCTION ISAMAX (N, X, INCX)
      INTEGER          ISAMAX, N, INCX
      REAL             X(*)
C
C MARGREET LOUTER-NOOL, CWI, 1987 OCTOBER 9
C
C PURPOSE
C =======
C
C ISAMAX CALCULATES THE SMALLEST INTEGER I, SUCH THAT
C
C    ABS (X ) = MAX [ ABS (X )  : J = 1, ..., N ]
C          I                J
C
C WHERE X IS A REAL VECTOR, AND INCX >= 1
C
C ===================================================
C
      BIT        REPXD,
     +           REP102 (65536) /32768 * B'10'/
      DESCRIPTOR REPXD
      INTEGER    IMAX, IMAX2, LASTX
      INTEGER    MAXVL /65535/
      REAL       SMAX, SMAX2
      REAL       XD
      DESCRIPTOR XD
C
C
      IF (N .LE. 0) THEN
         ISAMAX = 0
         RETURN
      ELSE IF (N .GT. 65535) THEN
         STOP 'BLAS ERROR (ISAMAX): N > 65535'
      ELSE IF (INCX .LT. 1) THEN
         STOP 'BLAS ERROR (ISAMAX): INCX < 1'
      END IF
C
C
      IF (INCX .EQ. 1) THEN
C
C        ............
C        . INCX = 1 .
C        ............
C
         ASSIGN XD, X(1;N)
C
         CALL Q8MAX (X'05',,XD,,IMAX,,)
C
         ISAMAX = IMAX + 1
C
         RETURN
C
      ELSE IF (INCX .GT. 2) THEN
C
C        ............
C        . INCX > 2 .
C        ............
C
         ASSIGN XD, .DYN.N
         XD = Q8VGATHP (X(1; 1), INCX, N; XD)
C
         CALL Q8MAX (X'05',,XD,,IMAX,,)
C
         ISAMAX = IMAX + 1
C
         RETURN
C
      ELSE
C
C        ..................................
C        . INCX = 2 => USE CONTROL VECTOR .
C        ..................................
C
         LASTX = (N-1) * INCX + 1
C
         IF (LASTX.LE.MAXVL) THEN
C
            ASSIGN XD, X(1;LASTX)
            ASSIGN REPXD, REP102(1;LASTX)
            CALL Q8MAX (X'05',,XD,,IMAX,REPXD,)
            ISAMAX = IMAX / 2 + 1
            RETURN
C
         ELSE
C
            ASSIGN XD, X(1;MAXVL)
            ASSIGN REPXD, REP102(1;MAXVL)
            CALL Q8MAX (X'05',,XD,,IMAX,REPXD,SMAX)
C
            LASTX = LASTX - MAXVL
            ASSIGN XD, X(MAXVL+1;LASTX)
            ASSIGN REPXD, REP102(2;LASTX)
            CALL Q8MAX (X'05',,XD,,IMAX2,,SMAX2)
C
            IF (SMAX2.GT.SMAX) THEN
               ISAMAX = (IMAX2 + MAXVL) / 2 + 1
            ELSE
               ISAMAX = IMAX / 2 + 1
            END IF
C
            RETURN
C
         END IF
C
      END IF
C
      RETURN
      END
      FUNCTION SASUM (N, X, INCX)
      INTEGER         N, INCX
      REAL            SASUM, X(*)
C
C MARGREET LOUTER-NOOL, CWI, 1987 OCTOBER 9
C
C PURPOSE
C =======
C
C SASUM CALCULATES THE SUM OF THE ABSOLUTE VALUES OF THE VECTOR X.
C
C    SASUM = SUM [ ABS (X ) ]  FOR I = 1, ..., N
C                        I
C
C WHERE X IS A REAL VECTOR, AND INCX >= 1
C
C ===================================================
C
      BIT        REPXD,
     +           REP102 (65536) /32768 * B'10'/
      DESCRIPTOR REPXD
      INTEGER    LASTX, NDIV2, NMIN1, NT64
      REAL       SUM
      REAL       ASUMD, XD, X1D, X2D
      DESCRIPTOR ASUMD, XD, X1D, X2D
C
C
      IF (N .LE. 0) THEN
         SASUM = 0.
         RETURN
      ELSE IF (N .EQ. 1) THEN
         SASUM = ABS (X(1))
         RETURN
      ELSE IF (N .GT. 65535) THEN
         STOP 'BLAS ERROR (SASUM): N > 65535'
      ELSE IF (INCX .LT. 1) THEN
         STOP 'BLAS ERROR (SASUM): INCX < 1'
      END IF
C
C
      LASTX = (N-1) * INCX + 1
C
      IF (INCX .NE. 2) THEN
C
         IF (INCX .EQ. 1) THEN
C
C           ............
C           . INCX = 1 .
C           ............
C
            ASSIGN XD, X(1;N)
C
         ELSE
C
C           ............
C           . INCX > 2 .
C           ............
C
            ASSIGN XD, .DYN.N
            XD = Q8VGATHP (X(1; 1), INCX, N; XD)
C
         END IF
C
         NDIV2 = N / 2
C
C THE VECTORS REFERENCED BY X1D AND X2D SHARE THE SAME STORAGE;
C X1D REFERENCES THE FIRST HALF OF XD AND X2D REFENCES THE
C SECOND HALF OF XD;
C
         NT64 = NDIV2 * 64
         CALL Q8PACK (NDIV2, XD, X1D)
         CALL Q8ADDX (X1D  ,NT64,X2D)
C
C THE VECTOR ASUMD BECOMES THE (ABSOLUTE) SUM OF THE VECTORS
C X1D AND X2D
C
         ASSIGN ASUMD, .DYN.NDIV2
         CALL Q8ADDNV (X'05',,X1D,,X2D,,ASUMD)
C
C SUM BECOMES THE SUM OF ALL ELEMENTS OF ASUMD;
C NOTE ASUMD IS OF LENGTH N / 2
C REGISTERS 4 AND 5 CONTAIN THE UNNORMALIZED DOUBLE PRECISION RESULT
C
         CALL Q8SUM  (X'00',,ASUMD,,,,4)
         CALL Q8ADDN (4,0,4)
         CALL Q8ADDN (4,5,SUM)
C
      ELSE
C
C        ..................................
C        . INCX = 2 => USE CONTROL VECTOR .
C        ..................................
C
         ASSIGN XD, X(1;LASTX)
C
         NMIN1 = N - 1
         ASSIGN REPXD, REP102 (1; NMIN1)
C
C THE VECTORS REFERENCED BY X1D AND X2D SHARE THE SAME STORAGE;
C X1D REFERENCES THE FIRST HALF OF XD AND X2D REFENCES THE
C SECOND HALF OF XD;
C
C IF N EVEN => 1 0 1 .... 0 1 + 1 0 1 .... 0 1
C IF N ODD  => 1 0 1 .... 0   + 1 0 1 .... 0
C
         IF (MOD (N,2) .NE. 0) THEN
            NT64 = NMIN1 * 64
         ELSE
            NT64 = N     * 64
         END IF
C
         CALL Q8PACK (NMIN1, XD, X1D)
         CALL Q8ADDX (X1D, NT64, X2D)
C
C THE VECTOR ASUMD BECOMES THE (ABSOLUTE) SUM OF THE VECTORS
C X1D AND X2D
C
         ASSIGN ASUMD, .DYN.NMIN1
         CALL Q8ADDNV (X'05',,X1D,,X2D,REPXD,ASUMD)
C
C SUM BECOMES THE SUM OF THE ELEMENTS OF ASUMD WITH PERIOD 1010...
C NOTE ASUMD IS OF LENGTH N - 1
C REGISTERS 4 AND 5 CONTAIN THE UNNORMALIZED DOUBLE PRECISION RESULT
C
         CALL Q8SUM  (X'00',,ASUMD,,,REPXD,4)
         CALL Q8ADDN (4,0,4)
         CALL Q8ADDN (4,5,SUM)
C
      END IF
C
C IN CASE OF N BEING ODD, SUM := SUM + ABS (X (LASTX))
C
      IF (MOD (N,2) .NE. 0) THEN
         SUM = SUM + ABS (X (LASTX))
      END IF
C
      SASUM = SUM
C
      RETURN
      END
      FUNCTION SNRM2 (N, X, INCX)
      INTEGER         N, INCX
      REAL            SNRM2, X(*)
C
C MARGREET LOUTER-NOOL, CWI, 1987 OCTOBER 9
C
C PURPOSE
C =======
C
C SNRM2 CALCULATES THE EUCLIDIAN LENGTH OR L -NORM OF A VECTOR
C                                           2
C                          2
C    SNRM2 = SQRT ( SUM [ X ] )   FOR I = 1, ..., N
C                          I
C
C WHERE X IS A REAL VECTOR, AND INCX >= 1
C
C ===================================================
C
      REAL             DD(2), H1, T1
      DOUBLE PRECISION DSUM, D1
      EQUIVALENCE      (DD(1), D1, H1), (DD(2), T1)
C
      BIT        REPXD,
     +           REP102 (65536) /32768 * B'10'/
      DESCRIPTOR REPXD
      INTEGER    LASTX, LEN, MAXVL /65535/, ONE /1/
      REAL       DOTP
      REAL       XD
      DESCRIPTOR XD
C
C
      IF (N .LE. 0) THEN
         SNRM2 = 0.
         RETURN
      ELSE IF (N .GT. 65535) THEN
         STOP 'BLAS ERROR (SNRM2) : N > 65535'
      ELSE IF (INCX .LT. 1) THEN
         STOP 'BLAS ERROR (SNRM2) : INCX < 1'
      END IF
C
C
      IF (INCX .NE. 2) THEN
C
         IF (INCX .EQ. 1) THEN
C
C           ............
C           . INCX = 1 .
C           ............
C
            ASSIGN XD, X(1;N)
C
         ELSE
C
C           ............
C           . INCX > 2 .
C           ............
C
            ASSIGN XD, .DYN.N
            XD = Q8VGATHP (X(1; 1), INCX, N; XD)
C
         END IF
C
C CALCULATE THE DOT PRODUCT DOTP = XD . XD;
C REGISTERS 4 AND 5 CONTAIN THE UNNORMALIZED DOUBLE PRECISION RESULT
C SNRM2 = SQRT (DOTP)
C
         CALL Q8DOTV (X'00',,XD,,XD,,4)
         CALL Q8ADDN (4,0,4)
         CALL Q8ADDN (4,5,DOTP)
C
         SNRM2 = SQRT (DOTP)
         RETURN
C
      ELSE
C
C        ..................................
C        . INCX = 2 => USE CONTROL VECTOR .
C        ..................................
C
         LASTX = (N-1) * INCX + 1
C
         ASSIGN REPXD, REP102 (1; LASTX)
C
C CALCULATE THE DOT PRODUCT DOTP = XD . XD;
C WITH PERIOD 1 0 1 0 ...
C REGISTERS 4 AND 5 CONTAIN THE UNNORMALIZED DOUBLE PRECISION RESULT
C
C IF LASTX > MAXVL (65535) THEN
C THE INTERVAL RESULTS ARE ADDED IN DOUBLE PRECISION
C
         DSUM = 0D0
         DO 1 I = 1, LASTX, MAXVL
C
            LEN = MIN0 (LASTX - I + 1, MAXVL)
C
            ASSIGN XD, X(I;LEN)
            CALL Q8PACK (LEN, REPXD, REPXD)
C
            CALL Q8DOTV (X'00',,XD,,XD,REPXD,4)
            CALL Q8ADDN (4,0,4)
            CALL Q8ADDN (4,5,H1)
            CALL Q8ADDL (4,5,T1)
            DSUM = DSUM + D1
C
            CALL Q8ADDX (REPXD, ONE, REPXD)
C
    1    CONTINUE
C
         SNRM2 = SQRT (SNGL(DSUM))
         RETURN
C
      END IF
C
      RETURN
      END
      FUNCTION SDOT (N, X, INCX, Y, INCY)
      INTEGER        N, INCX, INCY
      REAL           SDOT, X(*), Y(*)
C
C MARGREET LOUTER-NOOL, CWI, 1987 OCTOBER 9
C
C PURPOSE
C =======
C
C SDOT CALCULATES THE DOT PRODUCT OF THE VECTORS X AND Y
C
C    SDOT = SUM [ X . Y ]  FOR I = 1, ..., N
C                  I   I
C
C WHERE X AND Y ARE REAL VECTORS.
C
C ===========================================================
C
      REAL             DD(2), H1, T1
      DOUBLE PRECISION DSUM, D1
      EQUIVALENCE      (DD(1), D1, H1), (DD(2), T1)
C
      BIT        REPXYD,
     +           REP102 (65536) /32768 * B'10'/,
     +           REP103 (65535) /21845 * B'100'/,
     +           REP104 (65544) /16386 * B'1000'/
      DESCRIPTOR REPXYD
      INTEGER    ABSINCX, ABSINCY, FIRSTX, FIRSTY, IXIY
      INTEGER    JSHIFT, LASTXY, LEN, MAXVL /65535/
      REAL       DOTP
      REAL       XD, YD
      DESCRIPTOR XD, YD
C
C
      IF (N .LE. 0) THEN
         SDOT = 0.
         RETURN
      ELSE IF (N .GT. 65535) THEN
         STOP 'BLAS ERROR (SDOT) : N > 65535'
      ELSE IF (INCX .EQ. 0) THEN
         STOP 'BLAS ERROR (SDOT) : INCX = 0'
      ELSE IF (INCY .EQ. 0) THEN
         STOP 'BLAS ERROR (SDOT) : INCY = 0'
      END IF
C
C
      ABSINCX = IABS (INCX)
      ABSINCY = IABS (INCY)
      IXIY    = INCX * INCY
C
      IF (IXIY .EQ. 1) THEN
C
C        .........................
C        . INCX = INCY = 1 OR -1 .
C        .........................
C
         ASSIGN XD, X(1;N)
         ASSIGN YD, Y(1;N)
C
C CALCULATE THE DOT PRODUCT DOTP = XD . YD;
C REGISTERS 4 AND 5 CONTAIN THE UNNORMALIZED DOUBLE PRECISION RESULT
C
         CALL Q8DOTV (X'00',,XD,,YD,,4)
         CALL Q8ADDN (4,0,4)
         CALL Q8ADDN (4,5,DOTP)
C
         SDOT = DOTP
         RETURN
C
      ELSE IF (IXIY .EQ. -1) THEN
C
C        ...........................
C        . INCX = - INCY = 1 OR -1 .
C        ...........................
C
         ASSIGN XD, X(1;N)
         ASSIGN YD, .DYN.N
         YD = Q8VREV (Y(1;N); YD)
C
C CALCULATE THE DOT PRODUCT DOTP = XD . YD;
C REGISTERS 4 AND 5 CONTAIN THE UNNORMALIZED DOUBLE PRECISION RESULT
C
         CALL Q8DOTV (X'00',,XD,,YD,,4)
         CALL Q8ADDN (4,0,4)
         CALL Q8ADDN (4,5,DOTP)
C
         SDOT = DOTP
         RETURN
C
      ELSE IF ((INCX .EQ. INCY) .AND. (ABSINCX .LE. 4)) THEN
C
C        ..............................
C        . 1 < ABS (INCX = INCY) <= 4 .
C        ..............................
C
         LASTXY = (N-1) * ABSINCX + 1
C
         IF (ABSINCX .EQ. 2) THEN
            ASSIGN REPXYD, REP102(1;LASTXY)
         ELSE IF (ABSINCX .EQ. 3) THEN
            ASSIGN REPXYD, REP103(1;LASTXY)
         ELSE
            ASSIGN REPXYD, REP104(1;LASTXY)
         END IF
C
C CALCULATE THE DOT PRODUCT DOTP = XD . YD;
C REGISTERS 4 AND 5 CONTAIN THE UNNORMALIZED DOUBLE PRECISION RESULT
C
C IF LASTX > MAXVL (65535) THEN
C THE INTERVAL RESULTS ARE ADDED IN DOUBLE PRECISION
C
         DSUM = 0D0
         JSHIFT = MOD (MAXVL, ABSINCX)
         DO 1 I = 1, LASTXY, MAXVL
C
            LEN = MIN0 (LASTXY-I+1, MAXVL)
C
            CALL Q8PACK (LEN, REPXYD, REPXYD)
C
            ASSIGN XD, X(I;LEN)
            ASSIGN YD, Y(I;LEN)
C
            CALL Q8DOTV (X'00',,XD,,YD,REPXYD,4)
            CALL Q8ADDN (4,0,4)
            CALL Q8ADDN (4,5,H1)
            CALL Q8ADDL (4,5,T1)
            DSUM = DSUM + D1
C
            CALL Q8ADDX (REPXYD, JSHIFT, REPXYD)
C
    1    CONTINUE
C
         SDOT = SNGL (DSUM)
         RETURN
C
      ELSE
C
         IF (ABSINCX .EQ. 1) THEN
C
C           ...............
C           . ABSINCX = 1 .
C           ...............
C
            ASSIGN XD, X(1;N)
C
            ASSIGN YD, .DYN.N
            IF (IXIY .GT. 0) THEN
C
C              '''''''''''''''''''''''''''''''''''''''''''''''''''
C              ' INCX AND INCY BOTH POSITIVE OR BOTH NEGATIVE => '
C              ' GATHER Y FORWARDS;                              '
C              '''''''''''''''''''''''''''''''''''''''''''''''''''
C
               YD = Q8VGATHP (Y(1; 1), ABSINCY, N; YD)
C
            ELSE
C
C              '''''''''''''''''''''''''''''''''''''''''''''''''''
C              ' INCX AND INCY ARE OF DIFFERENT SIGN  =>         '
C              ' GATHER Y BACKWARDS;                             '
C              '''''''''''''''''''''''''''''''''''''''''''''''''''
C
               FIRSTY = 1 + ABSINCY * (N-1)
               YD = Q8VGATHP (Y(FIRSTY; 1), -ABSINCY, N; YD)
C
            END IF
C
         ELSE IF (ABSINCY .EQ. 1) THEN
C
C           ...............
C           . ABSINCY = 1 .
C           ...............
C
            ASSIGN YD, Y(1;N)
C
            ASSIGN XD, .DYN.N
            IF (IXIY .GT. 0) THEN
C
C              '''''''''''''''''''''''''''''''''''''''''''''''''''
C              ' INCX AND INCY BOTH POSITIVE OR BOTH NEGATIVE => '
C              ' GATHER X FORWARDS;                              '
C              '''''''''''''''''''''''''''''''''''''''''''''''''''
C
               XD = Q8VGATHP (X(1; 1), ABSINCX, N; XD)
C
            ELSE
C
C              '''''''''''''''''''''''''''''''''''''''''''''''''''
C              ' INCX AND INCY ARE OF DIFFERENT SIGN  =>         '
C              ' GATHER X BACKWARDS;                             '
C              '''''''''''''''''''''''''''''''''''''''''''''''''''
C
               FIRSTX = 1 + ABSINCX * (N-1)
               XD = Q8VGATHP (X(FIRSTX; 1), -ABSINCX, N; XD)
C
            END IF
C
         ELSE
C
C           ................................
C           . ABSINCX <> 1 OR ABSINCY <> 1 .
C           ................................
C
            ASSIGN XD, .DYN.N
            XD = Q8VGATHP (X(1; 1), ABSINCX, N; XD)
C
            ASSIGN YD, .DYN.N
C
            IF (IXIY .GT. 0) THEN
C
C              '''''''''''''''''''''''''''''''''''''''''''''''''''
C              ' INCX AND INCY BOTH POSITIVE OR BOTH NEGATIVE => '
C              ' GATHER Y FORWARDS;                              '
C              '''''''''''''''''''''''''''''''''''''''''''''''''''
C
               YD = Q8VGATHP (Y(1; 1), ABSINCY, N; YD)
C
            ELSE
C
C              '''''''''''''''''''''''''''''''''''''''''''''''''''
C              ' INCX AND INCY ARE OF DIFFERENT SIGN  =>         '
C              ' GATHER Y BACKWARDS;                             '
C              '''''''''''''''''''''''''''''''''''''''''''''''''''
C
               FIRSTY = 1 + ABSINCY * (N-1)
               YD = Q8VGATHP (Y(FIRSTY; 1), -ABSINCY, N; YD)
C
            END IF
C
         END IF
C
C CALCULATE THE DOT PRODUCT DOTP = XD . YD;
C REGISTERS 4 AND 5 CONTAIN THE UNNORMALIZED DOUBLE PRECISION RESULT
C
         CALL Q8DOTV (X'00',,XD,,YD,,4)
         CALL Q8ADDN (4,0,4)
         CALL Q8ADDN (4,5,DOTP)
C
         SDOT = DOTP
         RETURN
C
      END IF
C
      RETURN
      END
      SUBROUTINE SAXPY (N, SA, X, INCX, Y, INCY)
      INTEGER           N, INCX, INCY
      REAL              SA, X(*), Y(*)
C
C MARGREET LOUTER-NOOL, CWI, 1987 OCTOBER 9
C
C PURPOSE
C =======
C
C SAXPY PERFORMS THE FOLLOWING OPERATION
C
C    Y  := SA . X  + Y    FOR I = 1, ..., N
C     I          I    I
C
C WHERE X AND Y ARE REAL VECTORS, SA IS A REAL SCALAR.
C
C ====================================================
C
      BIT        REPXYD,
     +           REP102 (65536) /32768 * B'10'/,
     +           REP103 (65535) /21845 * B'100'/,
     +           REP104 (65544) /16386 * B'1000'/,
     +           REP105 (65535) /13107 * B'10000'/
      DESCRIPTOR REPXYD
      INTEGER    ABSINCX, ABSINCY, FIRSTX, FIRSTY, IXIY
      INTEGER    JSHIFT, LASTXY, LEN, MAXVL /65535/
      REAL       ZERO /0.E0/
      REAL       XD, YD
      DESCRIPTOR XD, YD
C
C
      IF (N .LE. 0 .OR. SA .EQ. ZERO) THEN
         RETURN
      ELSE IF (N .GT. 65535) THEN
         STOP 'BLAS ERROR (SAXPY) : N > 65535'
      ELSE IF (INCX .EQ. 0) THEN
         STOP 'BLAS ERROR (SAXPY) : INCX = 0'
      ELSE IF (INCY .EQ. 0) THEN
         STOP 'BLAS ERROR (SAXPY) : INCY = 0'
      END IF
C
C
      ABSINCX = IABS (INCX)
      ABSINCY = IABS (INCY)
      IXIY    = INCX * INCY
C
      IF (IXIY .EQ. 1) THEN
C
C        .........................
C        . INCX = INCY = 1 OR -1 .
C        .........................
C
         ASSIGN XD, X(1;N)
         ASSIGN YD, Y(1;N)
C
         YD = YD + SA * XD
C
         RETURN
C
      ELSE IF (IXIY .EQ. -1) THEN
C
C        ...........................
C        . INCX = - INCY = 1 OR -1 .
C        ...........................
C
         ASSIGN YD, Y(1;N)
C
         ASSIGN XD, .DYN.N
         XD = Q8VREV (X(1;N); XD)
C
         YD = YD + SA * XD
C
         RETURN
C
      ELSE IF ((INCX .EQ. INCY) .AND. (ABSINCY .LE. 5)) THEN
C
C        ..............................
C        . 1 < ABS (INCX = INCY) <= 5 .
C        ..............................
C
         LASTXY = (N-1) * ABSINCY + 1
C
         IF (ABSINCY .EQ. 2) THEN
            ASSIGN REPXYD, REP102(1;LASTXY)
         ELSE IF (ABSINCY .EQ. 3) THEN
            ASSIGN REPXYD, REP103(1;LASTXY)
         ELSE IF (ABSINCY .EQ. 4) THEN
            ASSIGN REPXYD, REP104(1;LASTXY)
         ELSE
            ASSIGN REPXYD, REP105(1;LASTXY)
         END IF
C
C IF LASTXY > MAXVL (65535) THEN
C SEVERAL LINKED TRIADS ARE PERFORMED
C
         JSHIFT = MOD (MAXVL, ABSINCY)
         DO 1 I = 1, LASTXY, MAXVL
C
            LEN = MIN0 (LASTXY-I+1, MAXVL)
C
            CALL Q8PACK (LEN, REPXYD, REPXYD)
C
            ASSIGN XD, X(I;LEN)
            ASSIGN YD, Y(I;LEN)
C
            WHERE (REPXYD)
               YD = YD + SA * XD
            END WHERE
C
            CALL Q8ADDX (REPXYD, JSHIFT, REPXYD)
C
    1    CONTINUE
C
         RETURN
C
      ELSE IF (ABSINCY .EQ. 1) THEN
C
C        ..................
C        . ABS (INCY) = 1 .
C        ..................
C
         ASSIGN YD, Y(1;N)
C
         ASSIGN XD, .DYN.N
         IF (IXIY .GT. 0) THEN
C
C           '''''''''''''''''''''''''''''''''''''''''''''''''''
C           ' INCX AND INCY BOTH POSITIVE OR BOTH NEGATIVE => '
C           ' GATHER X FORWARDS;                              '
C           '''''''''''''''''''''''''''''''''''''''''''''''''''
C
            XD = Q8VGATHP (X(1     ; 1),  ABSINCX, N; XD)
C
         ELSE
C
C           '''''''''''''''''''''''''''''''''''''''''''''''''''
C           ' INCX AND INCY ARE OF DIFFERENT SIGN  =>         '
C           ' GATHER X BACKWARDS;                             '
C           '''''''''''''''''''''''''''''''''''''''''''''''''''
C
            FIRSTX = 1 + ABSINCX * (N-1)
            XD = Q8VGATHP (X(FIRSTX; 1), -ABSINCX, N; XD)
C
         END IF
C
         YD = YD + SA * XD
C
         RETURN
C
      ELSE
C
C        ...............
C        . ABSINCY > 1 .
C        ...............
C
         IF (ABSINCX .EQ. 1) THEN
            ASSIGN XD, X(1;N)
         ELSE
            ASSIGN XD, .DYN.N
            XD = Q8VGATHP (X(1;1), ABSINCX, N; XD)
         END IF
C
         ASSIGN YD, .DYN.N
         IF (IXIY .GT. 0) THEN
C
C           '''''''''''''''''''''''''''''''''''''''''''''''''''
C           ' INCX AND INCY BOTH POSITIVE OR BOTH NEGATIVE => '
C           ' GATHER Y FORWARDS;                              '
C           '''''''''''''''''''''''''''''''''''''''''''''''''''
C
            YD = Q8VGATHP (Y(1; 1), ABSINCY, N; YD)
C
            YD = YD + SA * XD
C
            Y(1; 1) = Q8VSCATP (YD, ABSINCY, N; Y(1; 1))
C
            RETURN
C
         ELSE
C
C           '''''''''''''''''''''''''''''''''''''''''''''''''''
C           ' INCX AND INCY ARE OF DIFFERENT SIGN  =>         '
C           ' GATHER Y BACKWARDS;                             '
C           '''''''''''''''''''''''''''''''''''''''''''''''''''
C
            FIRSTY = 1 + ABSINCY * (N-1)
            YD = Q8VGATHP (Y(FIRSTY; 1), -ABSINCY, N; YD)
C
            YD = YD + SA * XD
C
            Y(FIRSTY; 1) = Q8VSCATP (YD, -ABSINCY, N; Y(FIRSTY; 1))
C
            RETURN
C
         END IF
C
      END IF
C
      RETURN
      END
      SUBROUTINE SSCAL (N, SA, X, INCX)
      INTEGER  N, INCX
      REAL     SA, X(*)
C
C MARGREET LOUTER-NOOL, CWI, 1987 OCTOBER 9
C
C PURPOSE
C =======
C
C SSCAL PERFORMS A VECTOR SCALING OPERATION
C
C    X  := SA . X         FOR I = 1, ..., N
C     I          I
C
C WHERE X IS A REAL VECTOR, SA IS A REAL SCALAR,
C AND INCX >= 1
C
C ====================================================
C
      BIT        REPXD,
     +           REP102 (65536) /32768 * B'10'/,
     +           REP103 (65535) /21845 * B'100'/,
     +           REP104 (65544) /16386 * B'1000'/
      DESCRIPTOR REPXD
      INTEGER    JSHIFT, LASTX, LEN, MAXVL /65535/
      REAL       XD
      DESCRIPTOR XD
      REAL       ONE /1.E0/
C
C
      IF (N .LE. 0 .OR. SA .EQ. ONE) THEN
         RETURN
      ELSE IF (N .GT. 65535) THEN
         STOP 'BLAS ERROR (SSCAL) : N > 65535'
      ELSE IF (INCX .LT. 1) THEN
         STOP 'BLAS ERROR (SSCAL) : INCX < 1'
      END IF
C
C
      IF (INCX .EQ. 1) THEN
C
C        ............
C        . INCX = 1 .
C        ............
C
         X(1;N) = SA * X(1;N)
         RETURN
C
      ELSE IF (INCX .LE. 4) THEN
C
C        .................
C        . 1 < INCX <= 4 .
C        .................
C
         LASTX = (N-1) * INCX + 1
C
         IF (INCX .EQ. 2) THEN
            ASSIGN REPXD, REP102(1;LASTX)
         ELSE IF (INCX .EQ. 3) THEN
            ASSIGN REPXD, REP103(1;LASTX)
         ELSE
            ASSIGN REPXD, REP104(1;LASTX)
         END IF
C
C IF LASTX > MAXVL (65535) THEN
C SEVERAL SCALAR-VECTOR OPERATIONS ARE PERFORMED
C
         JSHIFT = MOD (MAXVL, INCX)
         DO 1 I = 1, LASTX, MAXVL
C
            LEN = MIN0 (LASTX-I+1, MAXVL)
            CALL Q8PACK (LEN, REPXD, REPXD)
            ASSIGN XD, X(I;LEN)
C
            WHERE (REPXD)
               XD = SA * XD
            END WHERE
C
            CALL Q8ADDX (REPXD, JSHIFT, REPXD)
C
    1    CONTINUE
C
         RETURN
C
      ELSE
C
C        ............
C        . INCX > 4 .
C        ............
C
         ASSIGN XD, .DYN.N
         XD = Q8VGATHP (X(1; 1), INCX, N; XD)
C
         XD = SA * XD
C
         X(1; 1) = Q8VSCATP (XD, INCX, N; X(1; 1))
C
         RETURN
C
      END IF
C
      RETURN
      END
      SUBROUTINE SSWAP (N, X, INCX, Y, INCY)
      INTEGER    N, INCX, INCY
      REAL       X(*), Y(*)
C
C MARGREET LOUTER-NOOL, CWI, 1987 OCTOBER 9
C
C PURPOSE
C =======
C
C SSWAP INTERCHANGES THE VECTORS X AND Y
C
C
C    X  :=: Y          FOR I = 1, ..., N
C     I      I
C
C WHERE X AND Y ARE REAL VECTORS.
C
C ======================================
C
      INTEGER    ABSINCX, ABSINCY,
     +           FIRSTX, FIRSTY, IXIY
      REAL       XD, YD
      DESCRIPTOR XD, YD
C
C
      IF (N .LE. 0) THEN
         RETURN
      ELSE IF (N .GT. 65535) THEN
         STOP 'BLAS ERROR (SSWAP) : N > 65535'
      ELSE IF (INCX .EQ. 0) THEN
         STOP 'BLAS ERROR (SSWAP) : INCX = 0'
      ELSE IF (INCY .EQ. 0) THEN
         STOP 'BLAS ERROR (SSWAP) : INCY = 0'
      END IF
C
C
      ABSINCX = IABS (INCX)
      ABSINCY = IABS (INCY)
      IXIY    = INCX * INCY
C
      IF (IXIY .EQ. 1) THEN
C        .........................
C        . INCX = INCY = 1 OR -1 .
C        .........................
C
C STORE X IN XD; COPY Y IN X; COPY XD IN Y
C
         ASSIGN XD, .DYN.N
         XD = X(1;N)
         X(1;N) = Y(1;N)
         Y(1;N) = XD
         RETURN
C
      ELSE IF (IXIY .EQ. -1) THEN
C
C        ..........................
C        . INCX = -INCY = 1 OR -1 .
C        ..........................
C
C STORE REVERSE X IN XD; REVERSE Y IN X; COPY XD IN Y
C
         ASSIGN XD, .DYN.N
         XD = Q8VREV (X(1;N); XD)
         X(1;N) = Q8VREV (Y(1;N); X(1;N))
         Y(1;N) = XD
         RETURN
C
      ELSE IF (ABSINCX .EQ. 1) THEN
C
C        .....................................
C        . ABS (INCX) = 1 AND ABS (INCY) > 1 .
C        .....................................
C
C GATHER Y INTO YD; SCATTER X INTO Y; COPY YD IN X
C
         ASSIGN XD, X(1;N)
C
         ASSIGN YD, .DYN.N
         IF (IXIY .GT. 0) THEN
C
C           '''''''''''''''''''''''''''''''''''''''''''''''''''
C           ' INCX AND INCY BOTH POSITIVE OR BOTH NEGATIVE => '
C           ' GATHER Y FORWARDS;                              '
C           '''''''''''''''''''''''''''''''''''''''''''''''''''
C
            YD = Q8VGATHP (Y(1; 1), ABSINCY, N; YD)
            Y(1; 1) = Q8VSCATP (XD, ABSINCY, N; Y(1; 1))
            XD = YD
            RETURN
C
         ELSE
C
C           '''''''''''''''''''''''''''''''''''''''''''''''''''
C           ' INCX AND INCY ARE OF DIFFERENT SIGN  =>         '
C           ' GATHER Y BACKWARDS;                             '
C           '''''''''''''''''''''''''''''''''''''''''''''''''''
C
            FIRSTY = 1 + ABSINCY * (N-1)
            YD = Q8VGATHP (Y(FIRSTY; 1), -ABSINCY, N; YD)
            Y(FIRSTY; 1) = Q8VSCATP (XD, -ABSINCY, N; Y(FIRSTY; 1))
            XD = YD
            RETURN
C
         END IF
C
      ELSE IF (ABSINCY .EQ. 1) THEN
C
C        .....................................
C        . ABS (INCY) = 1 AND ABS (INCX) > 1 .
C        .....................................
C
C GATHER X INTO XD; SCATTER Y INTO X; COPY XD IN Y
C
         ASSIGN YD, Y(1;N)
C
         ASSIGN XD, .DYN.N
         IF (IXIY .GT. 0) THEN
C
C           '''''''''''''''''''''''''''''''''''''''''''''''''''
C           ' INCX AND INCY BOTH POSITIVE OR BOTH NEGATIVE => '
C           ' GATHER X FORWARDS;                              '
C           '''''''''''''''''''''''''''''''''''''''''''''''''''
C
            XD = Q8VGATHP (X(1; 1), ABSINCX, N; XD)
            X(1; 1) = Q8VSCATP (YD, ABSINCX, N; X(1; 1))
            YD = XD
            RETURN
C
         ELSE
C
C           '''''''''''''''''''''''''''''''''''''''''''''''''''
C           ' INCX AND INCY ARE OF DIFFERENT SIGN  =>         '
C           ' GATHER X BACKWARDS;                             '
C           '''''''''''''''''''''''''''''''''''''''''''''''''''
C
            FIRSTX = 1 + ABSINCX * (N-1)
            XD = Q8VGATHP (X(FIRSTX; 1), -ABSINCX, N; XD)
            X(FIRSTX; 1) = Q8VSCATP (YD, -ABSINCX, N; X(FIRSTX; 1))
            YD = XD
            RETURN
C
         END IF
C
      ELSE
C
C        ...................................
C        . ABS (INCX) > 1 ; ABS (INCY) > 1 .
C        ...................................
C
C GATHER X INTO XD; GATHER Y INTO YD;
C SCATTER XD INTO YD; SCATTER YD INTO Y
C
         ASSIGN XD, .DYN.N
         ASSIGN YD, .DYN.N
         IF (IXIY .GT. 0) THEN
C
C           '''''''''''''''''''''''''''''''''''''''''''''''''''
C           ' INCX AND INCY BOTH POSITIVE OR BOTH NEGATIVE => '
C           ' GATHER X AND Y FORWARDS;                        '
C           '''''''''''''''''''''''''''''''''''''''''''''''''''
C
            XD = Q8VGATHP (X(1; 1), ABSINCX, N; XD)
            YD = Q8VGATHP (Y(1; 1), ABSINCY, N; YD)
C
            X(1; 1) = Q8VSCATP (YD, ABSINCX, N; X(1; 1))
            Y(1; 1) = Q8VSCATP (XD, ABSINCY, N; Y(1; 1))
C
            RETURN
C
         ELSE
C
C           '''''''''''''''''''''''''''''''''''''''''''''''''''
C           ' INCX AND INCY ARE OF DIFFERENT SIGN  =>         '
C           ' GATHER X AND Y BACKWARDS;                       '
C           '''''''''''''''''''''''''''''''''''''''''''''''''''
C
            FIRSTX = 1 + ABSINCX * (N-1)
            FIRSTY = 1 + ABSINCY * (N-1)
            XD = Q8VGATHP (X(FIRSTX; 1), -ABSINCX, N; XD)
            YD = Q8VGATHP (Y(FIRSTY; 1), -ABSINCY, N; YD)
C
            X(1; 1) = Q8VSCATP (YD, ABSINCX, N; X(1; 1))
            Y(1; 1) = Q8VSCATP (XD, ABSINCY, N; Y(1; 1))
C
            RETURN
C
         END IF
C
      END IF
C
      RETURN
      END
      SUBROUTINE SCOPY (N, X, INCX, Y, INCY)
      INTEGER    N, INCX, INCY
      REAL       X(*), Y(*)
C
C MARGREET LOUTER-NOOL, CWI, 1987 OCTOBER 9
C
C PURPOSE
C =======
C
C SCOPY PERFORMS THE FOLLOWING COPY OPERATION
C
C
C    Y  := X           FOR I = 1, ..., N
C     I     I
C
C WHERE X AND Y ARE REAL VECTORS.
C
C ===========================================
C
      BIT        REPXYD,
     +           REP102 (65536) /32768 * B'10'/,
     +           REP103 (65535) /21845 * B'100'/
      DESCRIPTOR REPXYD
      INTEGER    ABSINCX, ABSINCY, FIRSTX, FIRSTY, IXIY
      INTEGER    JSHIFT, LASTXY, LEN, MAXVL /65535/
      REAL       XD
      DESCRIPTOR XD
C
C
      IF (N .LE. 0) THEN
         RETURN
      ELSE IF (N .GT. 65535) THEN
         STOP 'BLAS ERROR (SCOPY) : N > 65535'
      ELSE IF (INCX .EQ. 0) THEN
         STOP 'BLAS ERROR (SCOPY) : INCX = 0'
      ELSE IF (INCY .EQ. 0) THEN
         STOP 'BLAS ERROR (SCOPY) : INCY = 0'
      END IF
C
C
      ABSINCX = IABS (INCX)
      ABSINCY = IABS (INCY)
      IXIY    = INCX * INCY
C
      IF (IXIY .EQ. 1) THEN
C
C        .........................
C        . INCX = INCY = 1 OR -1 .
C        .........................
C
C STORE X IN Y
C
         Y(1;N) = X(1;N)
         RETURN
C
      ELSE IF (IXIY .EQ. -1) THEN
C
C        ..........................
C        . INCX = -INCY = 1 OR -1 .
C        ..........................
C
C STORE REVERSE X IN Y
C
         Y(1;N) = Q8VREV (X(1;N); Y(1;N))
         RETURN
C
      ELSE IF (ABSINCX .EQ. 1) THEN
C
C        .....................................
C        . ABS (INCX) = 1 AND ABS (INCY) > 1 .
C        .....................................
C
C SCATTER X INTO Y
C
         IF (IXIY .GT. 0) THEN
C
C           '''''''''''''''''''''''''''''''''''''''''''''''''''
C           ' INCX AND INCY BOTH POSITIVE OR BOTH NEGATIVE => '
C           ' GATHER X FORWARDS INTO Y                        '
C           '''''''''''''''''''''''''''''''''''''''''''''''''''
C
            Y(1; 1) = Q8VSCATP (X(1; 1), ABSINCY, N; Y(1; 1))
C
         ELSE
C
C           '''''''''''''''''''''''''''''''''''''''''''''''''''
C           ' INCX AND INCY ARE OF DIFFERENT SIGN  =>         '
C           ' GATHER X BACKWARDS INTO Y                       '
C           '''''''''''''''''''''''''''''''''''''''''''''''''''
C
            FIRSTY = 1 + ABSINCY * (N-1)
            Y(FIRSTY; 1) = Q8VSCATP (X(1; 1), -ABSINCY, N; Y(FIRSTY; 1))
C
         END IF
C
         RETURN
C
      ELSE IF (ABSINCY .EQ. 1) THEN
C
C        .....................................
C        . ABS (INCY) = 1 AND ABS (INCX) > 1 .
C        .....................................
C
C GATHER X INTO Y
C
         IF (IXIY .GT. 0) THEN
C
C           '''''''''''''''''''''''''''''''''''''''''''''''''''
C           ' INCX AND INCY BOTH POSITIVE OR BOTH NEGATIVE => '
C           ' GATHER X FORWARDS INTO Y                        '
C           '''''''''''''''''''''''''''''''''''''''''''''''''''
C
            Y(1; 1) = Q8VGATHP (X(1     ; 1),  ABSINCX, N; Y(1; 1))
C
         ELSE
C
C           '''''''''''''''''''''''''''''''''''''''''''''''''''
C           ' INCX AND INCY ARE OF DIFFERENT SIGN  =>         '
C           ' GATHER X BACKWARDS INTO Y                       '
C           '''''''''''''''''''''''''''''''''''''''''''''''''''
C
            FIRSTX = 1 + ABSINCX * (N-1)
            Y(1; 1) = Q8VGATHP (X(FIRSTX; 1), -ABSINCX, N; Y(1; 1))
C
         END IF
         RETURN
C
      ELSE IF ((INCX .EQ. INCY) .AND. (ABSINCX .LE. 3)) THEN
C
C        .............................
C        . 1 < ABS (INCX = INCY) <=3 .
C        .............................
C
         LASTXY = (N-1) * ABSINCX + 1
C
         IF (ABSINCX .EQ. 2) THEN
            ASSIGN REPXYD, REP102(1;LASTXY)
         ELSE
            ASSIGN REPXYD, REP103(1;LASTXY)
         END IF
C
C IF LASTXY > MAXVL (65535) THEN
C SEVERAL COPY OPERATIONS ARE PERFORMED
C
         JSHIFT = MOD (MAXVL, ABSINCX)
         DO 1 I = 1, LASTXY, MAXVL
C
            LEN = MIN0 (LASTXY-I+1, MAXVL)
            CALL Q8PACK (LEN, REPXYD, REPXYD)
C
            WHERE (REPXYD)
               Y(I;LEN) = X(I;LEN)
            END WHERE
C
            CALL Q8ADDX (REPXYD, JSHIFT, REPXYD)
C
    1    CONTINUE
C
         RETURN
C
      ELSE
C
C        ...................................
C        . ABS (INCX) > 1 ; ABS (INCY) > 1 .
C        ...................................
C
C GATHER X INTO XD; SCATTER XD INTO Y
C
         ASSIGN XD, .DYN.N
         XD = Q8VGATHP (X(1; 1), ABSINCX, N; XD)
C
         IF (IXIY .GT. 0) THEN
C
C           '''''''''''''''''''''''''''''''''''''''''''''''''''
C           ' INCX AND INCY BOTH POSITIVE OR BOTH NEGATIVE => '
C           ' SCATTER XD FORWARDS INTO Y                      '
C           '''''''''''''''''''''''''''''''''''''''''''''''''''
C
            Y(1; 1) = Q8VSCATP (XD, ABSINCY, N; Y(1; 1))
C
            RETURN
C
         ELSE
C
C           '''''''''''''''''''''''''''''''''''''''''''''''''''
C           ' INCX AND INCY ARE OF DIFFERENT SIGN  =>         '
C           ' SCATTER XD BACKWARDS INTO Y                     '
C           '''''''''''''''''''''''''''''''''''''''''''''''''''
C
            FIRSTY = 1 - INCY * (N-1)
            Y(FIRSTY; 1) = Q8VSCATP (XD, -ABSINCY, N; Y(FIRSTY; 1))
C
            RETURN
C
         END IF
C
      END IF
C
      RETURN
      END
      SUBROUTINE SROT (N, X, INCX, Y, INCY, SC, SS)
      INTEGER           N, INCX, INCY
      REAL              SC, SS, X(*), Y(*)
C
C MARGREET LOUTER-NOOL, CWI, 1987 OCTOBER 9
C
C PURPOSE
C =======
C
C SROT PERFORMS THE FOLLOWING OPERATION
C
C    X  := SC . X  + SS * Y    FOR I = 1, ..., N
C     I          I         I
C    Y  :=-SS . X  + SC * Y    FOR I = 1, ..., N
C     I          I         I
C
C WHERE X AND Y ARE REAL VECTORS, SC AND SS ARE REAL SCALARS.
C
C ===========================================================
C
      BIT        REPXYD,
     +           REP102 (65536) /32768 * B'10'/,
     +           REP103 (65535) /21845 * B'100'/
      DESCRIPTOR REPXYD
      INTEGER    ABSINCX, ABSINCY, FIRSTX, FIRSTY, IX, IXIY, IY
      INTEGER    JSHIFT, LASTXY, LEN, MAXVL /65535/
      REAL       XD, YD, ZD
      DESCRIPTOR XD, YD, ZD
      REAL       ZERO /0.E0/, ONE /1.E0/
C
C
      IF (N .LE. 0 .OR. (SS .EQ. ZERO .AND. SC. EQ. ONE)) THEN
         RETURN
      ELSE IF (N .GT. 65535) THEN
         STOP 'BLAS ERROR (SROT) : N > 65535'
      ELSE IF (INCX .EQ. 0) THEN
         STOP 'BLAS ERROR (SROT) : INCX = 0'
      ELSE IF (INCY .EQ. 0) THEN
         STOP 'BLAS ERROR (SROT) : INCY = 0'
      END IF
C
C
      ABSINCX = IABS (INCX)
      ABSINCY = IABS (INCY)
      IXIY    = INCX * INCY
C
      IF (IXIY .EQ. 1) THEN
C
C        .........................
C        . INCX = INCY = 1 OR -1 .
C        .........................
C
         ASSIGN XD, X(1;N)
         ASSIGN YD, Y(1;N)
C
C COMPUTE ZD = SS * (XD + YD)
C
         ASSIGN ZD, .DYN.N
         ZD = SS * (XD + YD)
C
C COMPUTE THE NEW XD AND YD
C
         XD = (SC - SS ) * XD + ZD
         YD = (SC + SS ) * YD - ZD
C
         RETURN
C
      ELSE IF (IXIY .EQ. -1) THEN
C
C        ...........................
C        . INCX = - INCY = 1 OR -1 .
C        ...........................
C
         ASSIGN XD, X(1;N)
         ASSIGN YD, .DYN.N
         YD = Q8VREV (Y(1;N); YD)
C
C COMPUTE ZD = SS * (XD + YD)
C
         ASSIGN ZD, .DYN.N
         ZD = SS * (XD + YD)
C
C COMPUTE THE NEW XD AND YD
C
         XD = (SC - SS ) * XD + ZD
         YD = (SC + SS ) * YD - ZD
C
         Y(1;N) = Q8VREV (YD; Y(1;N))
C
         RETURN
C
      ELSE IF ((INCX .EQ. INCY) .AND. (ABSINCX .LE. 3)) THEN
C
C        ..............................
C        . 1 < ABS (INCX = INCY) <= 3 .
C        ..............................
C
         LASTXY = (N-1) * ABSINCX + 1
         LEN = MIN0 (LASTXY, MAXVL)
         ASSIGN ZD, .DYN.LEN
C
         IF (ABSINCX .EQ. 2) THEN
            ASSIGN REPXYD, REP102(1;LASTXY)
         ELSE
            ASSIGN REPXYD, REP103(1;LASTXY)
         END IF
C
C IF LASTXY > MAXVL (65535) THEN
C THE OPERATIONS ARE PERFORMED SEVERAL TIMES
C
         JSHIFT = MOD (MAXVL, ABSINCX)
         DO 1 I = 1, LASTXY, MAXVL
C
            LEN = MIN0 (LASTXY, MAXVL)
            CALL Q8PACK (LEN,REPXYD, REPXYD)
            CALL Q8PACK (LEN, ZD, ZD)
C
            ASSIGN XD, X(I;LEN)
            ASSIGN YD, Y(I;LEN)
C
C COMPUTE ZD = SS * (XD + YD)
C
            ZD = SS * (XD + YD)
C
C COMPUTE THE NEW XD AND YD
C
            WHERE (REPXYD)
               XD = (SC - SS ) * XD + ZD
               YD = (SC + SS ) * YD - ZD
            END WHERE
C
            CALL Q8ADDX (REPXYD, JSHIFT, REPXYD)
C
    1    CONTINUE
C
         RETURN
C
      ELSE
C
C        ..............................
C        . ABSINCX > 1 OR ABSINCY > 1 .
C        ..............................
C
         IF (ABSINCX .EQ. 1) THEN
C
C           ...............
C           . ABSINCX = 1 .
C           ...............
C
            ASSIGN XD, X(1;N)
C
            IF (IXIY .GT. 0) THEN
C
C              '''''''''''''''''''''''''''''''''''''''''''''''''''
C              ' INCX AND INCY BOTH POSITIVE OR BOTH NEGATIVE => '
C              ' GATHER Y FORWARDS;                              '
C              '''''''''''''''''''''''''''''''''''''''''''''''''''
C
               FIRSTY = 1
               IY     = ABSINCY
C
            ELSE
C
C              '''''''''''''''''''''''''''''''''''''''''''''''''''
C              ' INCX AND INCY ARE OF DIFFERENT SIGN  =>         '
C              ' GATHER Y BACKWARDS;                             '
C              '''''''''''''''''''''''''''''''''''''''''''''''''''
C
               FIRSTY = 1 + ABSINCY * (N-1)
               IY     = - ABSINCY
C
            END IF
C
            ASSIGN YD, .DYN.N
            YD = Q8VGATHP (Y (FIRSTY; 1), IY, N; YD)
C
         ELSE IF (ABSINCY .EQ. 1) THEN
C
C           ...............
C           . ABSINCY = 1 .
C           ...............
C
            ASSIGN YD, Y(1;N)
C
            IF (IXIY .GT. 0) THEN
C
C              '''''''''''''''''''''''''''''''''''''''''''''''''''
C              ' INCX AND INCY BOTH POSITIVE OR BOTH NEGATIVE => '
C              ' GATHER X FORWARDS;                              '
C              '''''''''''''''''''''''''''''''''''''''''''''''''''
C
               FIRSTX = 1
               IX     = ABSINCX
C
            ELSE
C
C              '''''''''''''''''''''''''''''''''''''''''''''''''''
C              ' INCX AND INCY ARE OF DIFFERENT SIGN  =>         '
C              ' GATHER X BACKWARDS;                             '
C              '''''''''''''''''''''''''''''''''''''''''''''''''''
C
               FIRSTX = 1 + ABSINCX * (N-1)
               IX     = - ABSINCX
C
            END IF
C
            ASSIGN XD, .DYN.N
            XD = Q8VGATHP (X(FIRSTX; 1), IX, N; XD)
C
         ELSE
C
            IF (INCX .GT. 0) THEN
C
C              '''''''''''''''''''''''''''''''''
C              ' GATHER AND SCATTER X FORWARDS '
C              '''''''''''''''''''''''''''''''''
C
               FIRSTX = 1
C
            ELSE
C
C              ''''''''''''''''''''''''''''''''''
C              ' GATHER AND SCATTER X BACKWARDS '
C              ''''''''''''''''''''''''''''''''''
C
               FIRSTX = 1 - INCX * (N-1)
            END IF
C
            ASSIGN XD, .DYN.N
            XD = Q8VGATHP (X(FIRSTX; 1), INCX, N; XD)
C
            IF (INCY .GT. 0) THEN
C
C              '''''''''''''''''''''''''''''''''
C              ' GATHER AND SCATTER Y FORWARDS '
C              '''''''''''''''''''''''''''''''''
C
               FIRSTY = 1
C
            ELSE
C
C              ''''''''''''''''''''''''''''''''''
C              ' GATHER AND SCATTER Y BACKWARDS '
C              ''''''''''''''''''''''''''''''''''
C
               FIRSTY = 1 - INCY * (N-1)
            END IF
C
            ASSIGN YD, .DYN.N
            YD = Q8VGATHP (Y(FIRSTY; 1), INCY, N; YD)
C
         END IF
C
C COMPUTE ZD = SS * (XD + YD)
C
         ASSIGN ZD, .DYN.N
         ZD = SS * (XD + YD)
C
C COMPUTE THE NEW XD AND YD
C
         XD = (SC - SS ) * XD + ZD
         YD = (SC + SS ) * YD - ZD
C
C IF ABSINCX > 1 AND/OR ABSINCY > 1 THEN
C SCATTER XD INTO X AND/OR YD INTO Y;
C
         IF (ABSINCX .EQ. 1) THEN
            Y(FIRSTY; 1) = Q8VSCATP (YD, IY, N; Y(FIRSTY; 1))
         ELSE IF (ABSINCY .EQ. 1) THEN
            X(FIRSTX; 1) = Q8VSCATP (XD, IX, N; X(FIRSTX; 1))
         ELSE
            X(FIRSTX; 1) = Q8VSCATP (XD, INCX, N; X(FIRSTX; 1))
            Y(FIRSTY; 1) = Q8VSCATP (YD, INCY, N; Y(FIRSTY; 1))
         END IF
C
      END IF
C
      RETURN
      END
      SUBROUTINE SROT4 (N, X, INCX, Y, INCY, SC, SS)
      INTEGER           N, INCX, INCY
      REAL              SC, SS, X(*), Y(*)
C
C MARGREET LOUTER-NOOL, CWI, 1987 OCTOBER 9
C
C PURPOSE
C =======
C
C SROT PERFORMS THE FOLLOWING OPERATION
C
C    X  := SC . X  + SS * Y    FOR I = 1, ..., N
C     I          I         I
C    Y  :=-SS . X  + SC * Y    FOR I = 1, ..., N
C     I          I         I
C
C WHERE X AND Y ARE REAL VECTORS, SC AND SS ARE REAL SCALARS.
C
C ===========================================================
C
      BIT        REPXYD,
     +           REP102 (65536) /32768 * B'10'/
      DESCRIPTOR REPXYD
      INTEGER    ABSINCX, ABSINCY, FIRSTX, FIRSTY, IX, IXIY, IY
      INTEGER    JSHIFT, LASTXY, LEN, MAXVL /65535/
      REAL       XD, YD, VD, WD
      DESCRIPTOR XD, YD, VD, WD
      REAL       ZERO /0.E0/, ONE /1.E0/
C
C
      IF (N .LE. 0 .OR. (SS .EQ. ZERO .AND. SC. EQ. ONE)) THEN
         RETURN
      ELSE IF (N .GT. 65535) THEN
         STOP 'BLAS ERROR (SROT) : N > 65535'
      ELSE IF (INCX .EQ. 0) THEN
         STOP 'BLAS ERROR (SROT) : INCX = 0'
      ELSE IF (INCY .EQ. 0) THEN
         STOP 'BLAS ERROR (SROT) : INCY = 0'
      END IF
C
C
      ABSINCX = IABS (INCX)
      ABSINCY = IABS (INCY)
      IXIY    = INCX * INCY
C
      IF (IXIY .EQ. 1) THEN
C
C        .........................
C        . INCX = INCY = 1 OR -1 .
C        .........................
C
         ASSIGN XD, X(1;N)
         ASSIGN YD, Y(1;N)
C
C COMPUTE VD = SC * XD AND WD = SS * XD
C
         ASSIGN VD, .DYN.N
         ASSIGN WD, .DYN.N
         VD = SC * XD
         WD = SS * XD
C
C COMPUTE THE NEW XD AND YD
C
         XD = VD + SS * YD
         YD = SC * YD - WD
C
         RETURN
C
      ELSE IF (IXIY .EQ. -1) THEN
C
C        ...........................
C        . INCX = - INCY = 1 OR -1 .
C        ...........................
C
         ASSIGN XD, X(1;N)
         ASSIGN YD, .DYN.N
         YD = Q8VREV (Y(1;N); YD)
C
C COMPUTE VD = SC * XD AND WD = SS * XD
C
         ASSIGN VD, .DYN.N
         ASSIGN WD, .DYN.N
         VD = SC * XD
         WD = SS * XD
C
C COMPUTE THE NEW XD AND YD
C
         XD = VD + SS * YD
         YD = SC * YD - WD
C
         Y(1;N) = Q8VREV (YD; Y(1;N))
C
         RETURN
C
      ELSE IF ((INCX .EQ. INCY) .AND. (ABSINCX .LE. 2)) THEN
C
C        .........................
C        . INCX = INCY = 2 OR -2 .
C        .........................
C
         LASTXY = (N-1) * ABSINCX + 1
         LEN = MIN0 (LASTXY, MAXVL)
         ASSIGN VD, .DYN.LEN
         ASSIGN WD, .DYN.LEN
         ASSIGN REPXYD, REP102(1;LEN)
C
C IF LASTXY > MAXVL (65535) THEN
C THE OPERATIONS ARE PERFORMED SEVERAL TIMES
C
         JSHIFT = MOD (MAXVL, ABSINCX)
         DO 1 I = 1, LASTXY, MAXVL
C
            LEN = MIN0 (LASTXY, MAXVL)
            CALL Q8PACK (LEN, REPXYD, REPXYD)
            CALL Q8PACK (LEN, VD, VD)
            CALL Q8PACK (LEN, WD, WD)
C
            ASSIGN XD, X(I;LEN)
            ASSIGN YD, Y(I;LEN)
C
C COMPUTE VD = SC * XD AND WD = SS * XD
C
            VD = SC * XD
            WD = SS * XD
C
C COMPUTE THE NEW XD AND YD
C
            WHERE (REPXYD)
               XD = VD + SS * YD
               YD = SC * YD - WD
            END WHERE
C
            CALL Q8ADDX (REPXYD, JSHIFT, REPXYD)
C
    1    CONTINUE
C
         RETURN
C
      ELSE
C
C        ..............................
C        . ABSINCX > 1 OR ABSINCY > 1 .
C        ..............................
C
         IF (ABSINCX .EQ. 1) THEN
C
C           ...............
C           . ABSINCX = 1 .
C           ...............
C
            ASSIGN XD, X(1;N)
C
            IF (IXIY .GT. 0) THEN
C
C              '''''''''''''''''''''''''''''''''''''''''''''''''''
C              ' INCX AND INCY BOTH POSITIVE OR BOTH NEGATIVE => '
C              ' GATHER Y FORWARDS;                              '
C              '''''''''''''''''''''''''''''''''''''''''''''''''''
C
               FIRSTY = 1
               IY     = ABSINCY
C
            ELSE
C
C              '''''''''''''''''''''''''''''''''''''''''''''''''''
C              ' INCX AND INCY ARE OF DIFFERENT SIGN  =>         '
C              ' GATHER Y BACKWARDS;                             '
C              '''''''''''''''''''''''''''''''''''''''''''''''''''
C
               FIRSTY = 1 + ABSINCY * (N-1)
               IY     = - ABSINCY
C
            END IF
C
            ASSIGN YD, .DYN.N
            YD = Q8VGATHP (Y (FIRSTY; 1), IY, N; YD)
C
         ELSE IF (ABSINCY .EQ. 1) THEN
C
C           ...............
C           . ABSINCY = 1 .
C           ...............
C
            ASSIGN YD, Y(1;N)
C
            IF (IXIY .GT. 0) THEN
C
C              '''''''''''''''''''''''''''''''''''''''''''''''''''
C              ' INCX AND INCY BOTH POSITIVE OR BOTH NEGATIVE => '
C              ' GATHER X FORWARDS;                              '
C              '''''''''''''''''''''''''''''''''''''''''''''''''''
C
               FIRSTX = 1
               IX     = ABSINCX
C
            ELSE
C
C              '''''''''''''''''''''''''''''''''''''''''''''''''''
C              ' INCX AND INCY ARE OF DIFFERENT SIGN  =>         '
C              ' GATHER X BACKWARDS;                             '
C              '''''''''''''''''''''''''''''''''''''''''''''''''''
C
               FIRSTX = 1 + ABSINCX * (N-1)
               IX     = - ABSINCX
C
            END IF
C
            ASSIGN XD, .DYN.N
            XD = Q8VGATHP (X(FIRSTX; 1), IX, N; XD)
C
         ELSE
C
            IF (INCX .GT. 0) THEN
C
C              '''''''''''''''''''''''''''''''''
C              ' GATHER AND SCATTER X FORWARDS '
C              '''''''''''''''''''''''''''''''''
C
               FIRSTX = 1
C
            ELSE
C
C              ''''''''''''''''''''''''''''''''''
C              ' GATHER AND SCATTER X BACKWARDS '
C              ''''''''''''''''''''''''''''''''''
C
               FIRSTX = 1 - INCX * (N-1)
C
            END IF
C
            ASSIGN XD, .DYN.N
            XD = Q8VGATHP (X(FIRSTX; 1), INCX, N; XD)
C
            IF (INCY .GT. 0) THEN
C
C              '''''''''''''''''''''''''''''''''
C              ' GATHER AND SCATTER Y FORWARDS '
C              '''''''''''''''''''''''''''''''''
C
               FIRSTY = 1
C
            ELSE
C
C              ''''''''''''''''''''''''''''''''''
C              ' GATHER AND SCATTER Y BACKWARDS '
C              ''''''''''''''''''''''''''''''''''
C
               FIRSTY = 1 - INCY * (N-1)
C
            END IF
C
            ASSIGN YD, .DYN.N
            YD = Q8VGATHP (Y(FIRSTY; 1), INCY, N; YD)
C
         END IF
C
C COMPUTE VD = SC * XD AND WD = SS * XD
C
         ASSIGN VD, .DYN.N
         ASSIGN WD, .DYN.N
         VD = SC * XD
         WD = SS * XD
C
C COMPUTE THE NEW XD AND YD
C
         XD = VD + SS * YD
         YD = SC * YD - WD
C
C IF ABSINCX <> 1 AND/OR ABSINCY <> 1 THEN
C SCATTER XD INTO X AND/OR YD INTO Y;
C
         IF (ABSINCX .EQ. 1) THEN
            Y(FIRSTY; 1) = Q8VSCATP (YD, IY, N; Y(FIRSTY; 1))
         ELSE IF (ABSINCY .EQ. 1) THEN
            X(FIRSTX; 1) = Q8VSCATP (XD, IX, N; X(FIRSTX; 1))
         ELSE
            X(FIRSTX; 1) = Q8VSCATP (XD, INCX, N; X(FIRSTX; 1))
            Y(FIRSTY; 1) = Q8VSCATP (YD, INCY, N; Y(FIRSTY; 1))
         END IF
C
      END IF
C
      RETURN
      END
      FUNCTION ICAMAX (N, X, INCX)
      INTEGER          ICAMAX, N, INCX
      REAL             X(2, *)
C
C MARGREET LOUTER-NOOL, CWI, 1987 OCTOBER 9
C
C PURPOSE
C =======
C
C ICAMAX CALCULATES THE SMALLEST INTEGER I, SUCH THAT
C
C    ABS (X ) = MAX [ ABS (RE X ) + ABS (IM X )  : J = 1, ... N ]
C          I                   J             J
C
C WHERE X IS A COMPLEX VECTOR, AND INCX >= 1
C
C ===============================================================
C
      BIT         REP101D, REP110D, REP100D,
     +            REP110 (65535) /21845 * B'110'/,
     +            REP100 (65535) /21845 * B'100'/
      DESCRIPTOR  REP101D, REP110D, REP100D
      INTEGER     INCXT2
      INTEGER     JSHIFT, LEN, LENT2M1, LENT3M2, MAXVL /21845/
      REAL        ASUMD, ASUMD1, IMXD, REXD
      DESCRIPTOR  ASUMD, ASUMD1, IMXD, REXD
C
C
      IF (N .LE. 0) THEN
         ICAMAX = 0
         RETURN
      ELSE IF (N .GT. 32767) THEN
         STOP 'BLAS ERROR (ICAMAX) : N > 32767'
      ELSE IF (INCX .LT. 1) THEN
         STOP 'BLAS ERROR (ICAMAX) : INCX < 1'
      END IF
C
C
      IF (INCX .EQ. 1) THEN
C
C        ............
C        . INCX = 1 .
C        ............
C
         ASSIGN ASUMD, .DYN.N
         CALL Q8PACK (N, ASUMD, ASUMD1)
C
C 3N-2 IS LENGTH OF SPARSE VECTOR ADDITION
C IF 3N-2 > 65535 (MAXIMUM VECTOR LENGTH) OR N > 21845
C THEN TWO ADDITIONS ARE PERFORMED
C
         DO 1 I = 1, N, MAXVL
C
            LEN = MIN0 (N - I + 1, MAXVL)
            CALL Q8PACK (LEN, ASUMD1, ASUMD1)
C
            LENT2M1 = LEN * 2 - 1
            LENT3M2 = LEN * 3 - 2
C
            ASSIGN REXD, X(1,I; LENT2M1)
            ASSIGN IMXD, X(2,I; LENT2M1)
C
C MAKE BIT DESCRIPTORS FOR SPARSE ADDITION
C
            ASSIGN REP101D, REP110 (2; LENT3M2)
            ASSIGN REP110D, REP110 (1; LENT3M2)
            ASSIGN REP100D, REP100 (1; LENT3M2)
C
C Q8ADDNS : ASUM  := ABS (RE X ) + ABS (IM X )
C               I             I             I
C
            CALL Q8ADDNS (X'25',REP110D,REXD,REP101D,IMXD,
     +                    REP100D, ASUMD1)
C
            JSHIFT = LEN * 64
            CALL Q8ADDX (ASUMD1, JSHIFT, ASUMD1)
C
    1    CONTINUE
C
      ELSE
C
C        ............
C        . INCX > 1 .
C        ............
C
C REXD := REAL COMPONENTS OF X;
C IMXD := IMAGINARY COMPONENTS OF X;
C
         ASSIGN REXD, .DYN.N
         ASSIGN IMXD, .DYN.N
         INCXT2 = INCX * 2
         REXD = Q8VGATHP (X(1,1; 1), INCXT2, N; REXD)
         IMXD = Q8VGATHP (X(2,1; 1), INCXT2, N; IMXD)
C
C Q8ADDNV : ASUMD := ABS (REXD) + ABS (IMXD)
C
         ASSIGN ASUMD, .DYN.N
         CALL Q8ADDNV (X'05',,REXD,,IMXD,,ASUMD)
C
      END IF
C
      ICAMAX = Q8SMAXI (ASUMD) + 1
C
      RETURN
      END
      FUNCTION SCASUM (N, X, INCX)
      INTEGER          N, INCX
      REAL             SCASUM, X(2, *)
C
C MARGREET LOUTER-NOOL, CWI, 1987 OCTOBER 9
C
C PURPOSE
C =======
C
C SCASUM CALCULATES THE SUM OF THE ABOLUTE VALUES OF THE VECTOR X
C
C    SCASUM = SUM [ ABS (RE X ) + ABS (IM X ) ]  FOR I = 1, ... N
C                      I             I
C
C WHERE X IS A COMPLEX VECTOR, AND INCX >= 1
C
C ===============================================================
C
      INTEGER     INCXT2, JSHIFT, NT2
      REAL        SUM
      REAL        ASUMD, IMXD, REXD, XD
      DESCRIPTOR  ASUMD, IMXD, REXD, XD
C
C
      IF (N .LE. 0) THEN
         SCASUM = 0
         RETURN
      ELSE IF (N .GT. 32767) THEN
         STOP 'BLAS ERROR (SCASUM) : N > 32767'
      ELSE IF (INCX .LT. 1) THEN
         STOP 'BLAS ERROR (SCASUM) : INCX < 1'
      END IF
C
C
      IF (INCX .EQ. 1) THEN
C
C        ............
C        . INCX = 1 .
C        ............
C
C THE VECTORS REFERENCED BY XD, REXD AND IMXD SHARE
C THE SAME STORAGE;
C REXD REFERENCES TO THE FIRST !!!! N ELEMENTS OF X;
C IMXD REFERENCES TO THE LAST  !!!! N ELEMENTS OF X;
C
         NT2    = N * 2
         ASSIGN XD, X(1,1; NT2)
C
         JSHIFT   = N * 64
         CALL Q8PACK (N, XD, REXD)
         CALL Q8ADDX (REXD, JSHIFT, IMXD)
C
      ELSE
C
C        ............
C        . INCX > 1 .
C        ............
C
C REXD := REAL COMPONENTS OF X;
C IMXD := IMAGINARY COMPONENTS OF X;
C
         ASSIGN REXD, .DYN.N
         ASSIGN IMXD, .DYN.N
         INCXT2 = INCX * 2
         REXD = Q8VGATHP (X(1,1; 1), INCXT2, N; REXD)
         IMXD = Q8VGATHP (X(2,1; 1), INCXT2, N; IMXD)
C
      END IF
C
C ASUMD := ABS (REXD) + ABS (IMXD)
C
      ASSIGN ASUMD, .DYN.N
      CALL Q8ADDNV (X'05',,REXD,,IMXD,,ASUMD)
C
C REGISTERS 4 AND 5 CONTAIN THE UNNORMALIZED DOUBLE PRECISION RESULT
C
      CALL Q8SUM   (X'00',,ASUMD,,,,4)
      CALL Q8ADDN  (4,0,4)
      CALL Q8ADDN  (4,5,SUM)
C
      SCASUM = SUM
C
      RETURN
      END
      FUNCTION SCNRM2 (N, X, INCX)
      INTEGER          N, INCX
      REAL             SCNRM2, X(2, *)
C
C MARGREET LOUTER-NOOL, CWI, 1987 OCTOBER 9
C
C PURPOSE
C =======
C
C SCNRM2 CALCULATES THE EUCLIDIAN LENGTH OR L -NORM OF A VECTOR
C                                            2
C                           2
C    SCNRM2 = SQRT ( SUM [ X ] )  FOR I = 1, ..., N
C                           I
C
C WHERE X IS A COMPLEX VECTOR, AND INCX >= 1
C
C =============================================================
C
      INTEGER     INCXT2, JSHIFT, NT2
      REAL        DOTP
      REAL        IMXD, REXD, XD
      DESCRIPTOR  IMXD, REXD, XD
C
C
      IF (N .LE. 0) THEN
         SCNRM2 = 0
         RETURN
      ELSE IF (N .GT. 32767) THEN
         STOP 'BLAS ERROR (SCNRM2) : N > 32767'
      ELSE IF (INCX .LT. 1) THEN
         STOP 'BLAS ERROR (SCNRM2) : INCX < 1'
      END IF
C
C
      NT2    = N * 2
      INCXT2 = INCX * 2
C
      IF (INCX .EQ. 1) THEN
C
C        ............
C        . INCX = 1 .
C        ............
C
         ASSIGN XD, X(1,1; NT2)
C
      ELSE
C
C        ............
C        . INCX > 1 .
C        ............
C
C THE VECTORS REFERENCED BY XD, REXD AND IMXD SHARE THE SAME STORAGE;
C XD BECOMES REXD FOLLOWED BY IMXD;
C
C REXD REFERENCES TO THE REAL ELEMENTS OF X;
C IMXD REFERENCES TO THE IMAGINARY ELEMENTS OF X;
C
         JSHIFT = N * 64
         ASSIGN XD, .DYN.NT2
         CALL Q8PACK (N, XD, REXD)
         CALL Q8ADDX (REXD, JSHIFT, IMXD)
C
         REXD = Q8VGATHP (X(1,1; 1), INCXT2, N; REXD)
         IMXD = Q8VGATHP (X(2,1; 1), INCXT2, N; IMXD)
C
      END IF
C
C CALCULATE THE DOT PRODUCT DOTP = XD . XD;
C REGISTERS 4 AND 5 CONTAIN THE UNNORMALIZED DOUBLE PRECISION RESULT
C
C SCNRM2 = SQRT (DOTP)
C
      CALL Q8DOTV (X'00',,XD,,XD,,4)
      CALL Q8ADDN (4,0,4)
      CALL Q8ADDN (4,5,DOTP)
C
      SCNRM2 = SQRT (DOTP)
C
      RETURN
      END
      FUNCTION CDOTU (N, X, INCX, Y, INCY)
      INTEGER         N, INCX, INCY
      REAL            X(2, *), Y(2, *)
      COMPLEX         CDOTU
C
C MARGREET LOUTER-NOOL, CWI, 1987 OCTOBER 9
C
C PURPOSE
C =======
C
C CDOTU CALCULATES THE DOT PRODUCT OF THE VECTORS X AND Y
C
C    CDOTU = SUM [ X . Y ]   FOR I = 1, ..., N
C                   I   I
C
C WHERE X AND Y ARE COMPLEX VECTORS.
C
C =======================================================
C
      REAL             DD(4), H1, H2, T1, T2
      DOUBLE PRECISION D1, D2
      EQUIVALENCE      (DD(1), D1, H1), (DD(2), T1),
     +                 (DD(3), D2, H2), (DD(4), T2)
C
      BIT         REP10D,
     +            REP10 (65536) /32768 * B'10'/
      DESCRIPTOR  REP10D
      INTEGER     FIRSTX, FIRSTY, INCXT2, INCYT2, IXIY, NT2
      INTEGER     ONE /1/, JSHIFT
      REAL        IMDOT, REDOT
      REAL        IMXD, IMYD, REXD, REYD, XD, YBARD, YD
      DESCRIPTOR  IMXD, IMYD, REXD, REYD, XD, YBARD, YD
C
C
      IF (N .LE. 0) THEN
         CDOTU = CMPLX (0, 0)
         RETURN
      ELSE IF (N .GT. 32767) THEN
         STOP 'BLAS ERROR (CDOTU) : N > 32767'
      ELSE IF (INCX .EQ. 0) THEN
         STOP 'BLAS ERROR (CDOTU) : INCX = 0'
      ELSE IF (INCY .EQ. 0) THEN
         STOP 'BLAS ERROR (CDOTU) : INCY = 0'
      END IF
C
      IXIY    = INCX * INCY
      NT2     = N * 2
C
      IF (IXIY .EQ. 1) THEN
C
C        ..........................
C        . INCX =  INCY = 1 OR -1 .
C        ..........................
C
         ASSIGN XD, X(1,1; NT2)
         ASSIGN YD, Y(1,1; NT2)
C
C RE CDOTU := DOT PRODUKT (RE X, RE Y) - DOT PRODUKT (IM X, IM Y)
C
C Q8DOTV : COMPUTES THE DOUBLE PRECISION DOT PRODUKT OF TWO VECTORS
C REGISTERS 4 AND 5 CONTAIN THE UNNORMALIZED DOUBLE PRECISION RESULT
C
         ASSIGN REP10D, REP10(1;NT2)
         CALL Q8DOTV (X'00',   ,XD,   ,YD,REP10D,4)
         CALL Q8ADDN (4,0,4)
         CALL Q8ADDN (4,5,H1)
         CALL Q8ADDL (4,5,T1)
         CALL Q8DOTV (X'00',ONE,XD,ONE,YD,REP10D,4)
         CALL Q8ADDN (4,0,4)
         CALL Q8ADDN (4,5,H2)
         CALL Q8ADDL (4,5,T2)
         REDOT = SNGL (D1 - D2)
C
C IM CDOTU := DOT PRODUKT (RE X, IM Y) + DOT PRODUKT (IM X, RE Y)
C
         CALL Q8DOTV (X'00',   ,XD,ONE,YD,REP10D,4)
         CALL Q8ADDN (4,0,4)
         CALL Q8ADDN (4,5,H1)
         CALL Q8ADDL (4,5,T1)
         CALL Q8DOTV (X'00',ONE,XD,   ,YD,REP10D,4)
         CALL Q8ADDN (4,0,4)
         CALL Q8ADDN (4,5,H2)
         CALL Q8ADDL (4,5,T2)
         IMDOT = SNGL (D1 + D2)
C
         CDOTU = CMPLX (REDOT,IMDOT)
C
         RETURN
C
      ELSE IF (IXIY .EQ. -1) THEN
C
C        ..........................
C        . INCX = -INCY = 1 OR -1 .
C        ..........................
C
C IF INCX = -INCY THEN XD := REVERSE OF X =>
C XD = IM X  RE X  IM X    RE X    ......
C          N     N     N-1     N-1
C
C THIS YIELDS, ONLY ONE INPRODUCT WITHOUT A CONTROL VECTOR
C FOR THE IMAGINARY PART OF CDOTU WILL SATISFY.
C
         ASSIGN XD, .DYN.NT2
         XD = Q8VREV (X(1,1;NT2); XD)
         ASSIGN YD, Y(1,1; NT2)
C
C RE CDOTU := DOT PRODUKT (RE X, RE Y) - DOT PRODUKT (IM X, IM Y)
C
C Q8DOTV : COMPUTES THE DOUBLE PRECISION DOT PRODUKT OF TWO VECTORS
C REGISTERS 4 AND 5 CONTAIN THE UNNORMALIZED DOUBLE PRECISION RESULT
C
         ASSIGN REP10D, REP10(1;NT2)
         CALL Q8DOTV (X'00',ONE,XD,   ,YD,REP10D,4)
         CALL Q8ADDN (4,0,4)
         CALL Q8ADDN (4,5,H1)
         CALL Q8ADDL (4,5,T1)
         CALL Q8DOTV (X'00',   ,XD,ONE,YD,REP10D,4)
         CALL Q8ADDN (4,0,4)
         CALL Q8ADDN (4,5,H2)
         CALL Q8ADDL (4,5,T2)
         REDOT = SNGL (D1 - D2)
C
C IM CDOTU := DOT PRODUKT (RE X, IM Y) + DOT PRODUKT (IM X, RE Y)
C
         CALL Q8DOTV (X'00',   ,XD,   ,YD,      ,4)
         CALL Q8ADDN (4,0,4)
         CALL Q8ADDN (4,5,IMDOT)
C
         CDOTU = CMPLX (REDOT,IMDOT)
C
         RETURN
C
      ELSE
C
C        ....................................
C        . ABS (INCX) > 1 OR ABS (INCY) > 1 .
C        ....................................
C
C THE VECTORS REFERENCED BY XD, REXD AND IMXD SHARE THE SAME STORAGE;
C XD BECOMES REXD FOLLOWED BY IMXD;
C
C REXD REFERENCES TO THE REAL ELEMENTS OF X;
C IMXD REFERENCES TO THE IMAGINARY ELEMENTS OF X;
C
         JSHIFT = N * 64
         ASSIGN XD, .DYN.NT2
         CALL Q8PACK (N, XD, REXD)
         CALL Q8ADDX (REXD, JSHIFT, IMXD)
C
         IF (INCX .GT. 0) THEN
            FIRSTX = 1
         ELSE
            FIRSTX = 1 - INCX * (N-1)
         END IF
C
         INCXT2  = INCX * 2
         REXD = Q8VGATHP (X(1, FIRSTX; 1), INCXT2, N; REXD)
         IMXD = Q8VGATHP (X(2, FIRSTX; 1), INCXT2, N; IMXD)
C
C THE VECTORS REFERENCED BY YBARD, REYD AND IMYD SHARE THE SAME STORAGE;
C YBARD BECOMES IMYD FOLLOWED BY REYD;
C
C REYD REFERENCES TO THE REAL ELEMENTS OF Y;
C IMYD REFERENCES TO THE IMAGINARY ELEMENTS OF Y;
C
C        JSHIFT = N * 64
         ASSIGN YBARD, .DYN.NT2
         CALL Q8PACK (N, YBARD, IMYD)
         CALL Q8ADDX (IMYD, JSHIFT, REYD)
C
         IF (INCY .GT. 0) THEN
            FIRSTY = 1
         ELSE
            FIRSTY = 1 - INCY * (N-1)
         END IF
C
         INCYT2  = INCY * 2
         REYD = Q8VGATHP (Y(1, FIRSTY; 1), INCYT2, N; REYD)
         IMYD = Q8VGATHP (Y(2, FIRSTY; 1), INCYT2, N; IMYD)
C
C RE CDOTU := DOT PRODUKT (RE X, RE Y) - DOT PRODUKT (IM X, IM Y)
C
C Q8DOTV : COMPUTES THE DOUBLE PRECISION DOT PRODUKT OF TWO VECTORS
C REGISTERS 4 AND 5 CONTAIN THE UNNORMALIZED DOUBLE PRECISION RESULT
C
         CALL Q8DOTV (X'00',,REXD,, REYD,,4)
         CALL Q8ADDN (4,0,4)
         CALL Q8ADDN (4,5,H1)
         CALL Q8ADDL (4,5,T1)
         CALL Q8DOTV (X'00',,IMXD,, IMYD,,4)
         CALL Q8ADDN (4,0,4)
         CALL Q8ADDN (4,5,H2)
         CALL Q8ADDL (4,5,T2)
         REDOT = SNGL (D1 - D2)
C
C IM CDOTU := DOT PRODUKT (RE X, IM Y) + DOT PRODUKT (IM X, RE Y)
C
         CALL Q8DOTV (X'00',,  XD,,YBARD,,4)
         CALL Q8ADDN (4,0,4)
         CALL Q8ADDN (4,5,IMDOT)
C
         CDOTU = CMPLX (REDOT, IMDOT)
C
      END IF
C
      RETURN
      END
      FUNCTION CDOTC (N, X, INCX, Y, INCY)
      INTEGER         N, INCX, INCY
      REAL            X(2, *), Y(2, *)
      COMPLEX         CDOTC
C
C MARGREET LOUTER-NOOL, CWI, 1987 OCTOBER 9
C
C PURPOSE
C =======
C                                                 _
C CDOTC CALCULATES THE DOT PRODUCT OF THE VECTORS X AND Y
C                  _
C    CDOTC = SUM [ X . Y ] FOR I = 1, ..., N )
C                   I   I
C
C WHERE X AND Y ARE COMPLEX VECTORS.
C
C =======================================================
C
      REAL             DD(4), H1, H2, T1, T2
      DOUBLE PRECISION D1, D2
      EQUIVALENCE      (DD(1), D1, H1), (DD(2), T1),
     +                 (DD(3), D2, H2), (DD(4), T2)
C
      BIT         REP10D,
     +            REP10 (65536) /32768 * B'10'/
      DESCRIPTOR  REP10D
      INTEGER     FIRSTX, FIRSTY, INCXT2, INCYT2, IXIY, NT2, NT2M1
      REAL        IMDOT, REDOT
      INTEGER     I64 /64/, JSHIFT, ONE /1/
      REAL        REXD, REYD, IMXD, IMYD, XD, YD
      DESCRIPTOR  REXD, REYD, IMXD, IMYD, XD, YD
C
C
      IF (N .LE. 0) THEN
         CDOTC = CMPLX (0, 0)
         RETURN
      ELSE IF (N .GT. 32767) THEN
         STOP 'BLAS ERROR (CDOTC) : N > 32767'
      ELSE IF (INCX .EQ. 0) THEN
         STOP 'BLAS ERROR (CDOTC) : INCX = 0'
      ELSE IF (INCY .EQ. 0) THEN
         STOP 'BLAS ERROR (CDOTC) : INCY = 0'
      END IF
C
C
      IXIY    = INCX * INCY
      NT2     = N * 2
C
      IF (IABS (IXIY) .EQ. 1) THEN
C
         ASSIGN XD, X(1,1; NT2)
C
         IF (IXIY .EQ. 1) THEN
C
C           ..........................
C           . INCX =  INCY = 1 OR -1 .
C           ..........................
C
            ASSIGN YD, Y(1,1; NT2)
C
         ELSE
C           ........................................
C           . INCX = -INCY = 1 OR -1               .
C           . YD BECOMES THE REVERSE OF Y(1,1;NT2) .
C           ........................................
C
            NT2M1 = NT2 - 1
            ASSIGN YD, .DYN.NT2
            CALL Q8PACK (NT2M1, YD, REYD)
            CALL Q8ADDX (REYD, I64, IMYD)
C
            ASSIGN REP10D, REP10(1;NT2M1)
            CALL Q8VREVV (X'00',,Y(1,1;NT2M1),,,REP10D,REYD)
            CALL Q8VREVV (X'00',,Y(2,1;NT2M1),,,REP10D,IMYD)
C
         END IF
C
C RE CDOTC := DOT PRODUKT (RE X, RE Y) + DOT PRODUCT (IM X, IM Y)
C
C Q8DOTV : COMPUTES THE DOUBLE PRECISION DOT PRODUKT OF TWO VECTORS
C REGISTERS 4 AND 5 CONTAIN THE UNNORMALIZED DOUBLE PRECISION RESULT
C
         ASSIGN REP10D, REP10(1;NT2)
         CALL Q8DOTV (X'00',   ,XD,   ,YD,      ,4)
         CALL Q8ADDN (4,0,4)
         CALL Q8ADDN (4,5,REDOT)
C
C IM CDOTC := DOT PRODUKT (RE X, IM Y) - DOT PRODUCT (IM X, RE Y)
C
         CALL Q8DOTV (X'00',   ,XD,ONE,YD,REP10D,4)
         CALL Q8ADDN (4,0,4)
         CALL Q8ADDN (4,5,H1)
         CALL Q8ADDL (4,5,T1)
         CALL Q8DOTV (X'00',ONE,XD,   ,YD,REP10D,4)
         CALL Q8ADDN (4,0,4)
         CALL Q8ADDN (4,5,H2)
         CALL Q8ADDL (4,5,T2)
         IMDOT = SNGL (D1 - D2)
C
         CDOTC = CMPLX (REDOT, IMDOT)
C
         RETURN
C
      ELSE
C
C        ...........................................
C        . NOT (ABS (INCX) = 1 AND ABS (INCY) = 1) .
C        ...........................................
C
C THE VECTORS REFERENCED BY XD, REXD AND IMXD SHARE THE SAME STORAGE;
C XD BECOMES REXD FOLLOWED BY IMXD;
C
C REXD REFERENCES TO THE REAL ELEMENTS OF X;
C IMXD REFERENCES TO THE IMAGINARY ELEMENTS OF X;
C
         JSHIFT = N * 64
         ASSIGN XD, .DYN.NT2
         CALL Q8PACK (N, XD, REXD)
         CALL Q8ADDX (REXD, JSHIFT, IMXD)
C
         IF (INCX .GT. 0) THEN
            FIRSTX = 1
         ELSE
            FIRSTX = 1 - INCX * (N-1)
         END IF
C
         INCXT2 = INCX * 2
         REXD = Q8VGATHP (X(1, FIRSTX; 1), INCXT2, N; REXD)
         IMXD = Q8VGATHP (X(2, FIRSTX; 1), INCXT2, N; IMXD)
C
C THE VECTORS REFERENCED BY YD, REYD AND IMYD SHARE THE SAME STORAGE;
C YD BECOMES REYD FOLLOWED BY IMYD;
C
C REYD REFERENCES TO THE REAL ELEMENTS OF Y;
C IMYD REFERENCES TO THE IMAGINARY ELEMENTS OF Y;
C
         ASSIGN YD, .DYN.NT2
         CALL Q8PACK (N, YD, REYD)
         CALL Q8ADDX (REYD, JSHIFT, IMYD)
C
         IF (INCY .GT. 0) THEN
            FIRSTY = 1
         ELSE
            FIRSTY = 1 - INCY * (N-1)
         END IF
C
         INCYT2 = INCY * 2
         REYD = Q8VGATHP (Y(1, FIRSTY; 1), INCYT2, N; REYD)
         IMYD = Q8VGATHP (Y(2, FIRSTY; 1), INCYT2, N; IMYD)
C
C RE CDOTC := DOT PRODUKT (RE X, RE Y) + DOT PRODUCT (IM X, IM Y)
C
C Q8DOTV : COMPUTES THE DOUBLE PRECISION DOT PRODUKT OF TWO VECTORS
C REGISTERS 4 AND 5 CONTAIN THE UNNORMALIZED DOUBLE PRECISION RESULT
C
         CALL Q8DOTV (X'00',,  XD,,  YD,,4)
         CALL Q8ADDN (4,0,4)
         CALL Q8ADDN (4,5,REDOT)
C
C IM CDOTC := DOT PRODUKT (RE X, IM Y) - DOT PRODUCT (IM X, RE Y)
C
         CALL Q8DOTV (X'00',,REXD,,IMYD,,4)
         CALL Q8ADDN (4,0,4)
         CALL Q8ADDN (4,5,H1)
         CALL Q8ADDL (4,5,T1)
         CALL Q8DOTV (X'00',,IMXD,,REYD,,4)
         CALL Q8ADDN (4,0,4)
         CALL Q8ADDN (4,5,H2)
         CALL Q8ADDL (4,5,T2)
         IMDOT = SNGL (D1 - D2)
C
         CDOTC = CMPLX (REDOT, IMDOT)
C
         RETURN
C
      END IF
C
      RETURN
      END
      SUBROUTINE CAXPY (N, CA, X, INCX, Y, INCY)
      INTEGER           N, INCX, INCY
      REAL              X(2, *), Y(2, *)
      COMPLEX           CA
C
C MARGREET LOUTER-NOOL, CWI, 1987 OCTOBER 9
C
C PURPOSE
C =======
C
C CAXPY PERFORMS THE FOLLOWING OPERATION
C
C    Y := CA . X  +  Y      FOR I = 1, ..., N
C     I         I     I
C
C WHERE X AND Y ARE COMPLEX VECTORS, CA IS A COMPLEX SCALAR;
C INCX <> 0 AND INCY <> 0.
C
C ==========================================================
C
      BIT         REP10D,
     +            REP10 (65536) /32768 * B'10'/
      DESCRIPTOR  REP10D
      INTEGER     ABSINCX, ABSINCY, FIRSTX, FIRSTY,
     +            INCXT2, INCYT2, IXIY, IXT2,
     +            NT2, NT2M1
      INTEGER     I64 /64/, JSHIFT
      REAL        AI, RE
      REAL        IMXD, IMYD, REXD, REYD, XD, YD
      DESCRIPTOR  IMXD, IMYD, REXD, REYD, XD, YD
      COMPLEX     ZERO /(0.E0, 0.E0)/
C
C
      IF (N .LE. 0 .OR. CA .EQ. ZERO) THEN
         RETURN
      ELSE IF (N .GT. 32767) THEN
         STOP 'BLAS ERROR (CAXPY) : N > 32767'
      ELSE IF (INCX .EQ. 0) THEN
         STOP 'BLAS ERROR (CAXPY) : INCX = 0'
      ELSE IF (INCY .EQ. 0) THEN
         STOP 'BLAS ERROR (CAXPY) : INCY = 0'
      END IF
C
C ------------------------------------------------------------
C CAXPY PERFORMS THE NEXT OPERATIONS :
C
C    YD =   YD  + RE . XD
C
C RE Y  = RE Y  - AI * IM X
C     I       I            I
C IM Y  = IM Y  + AI * RE X
C     I       I            I
C
C ABSINCY = 1   => USE THE MIXED     MODE R  I  R  I   . . .
C                                          1  1  2  2
C ABSINCY > 1   => USE THE SEPERATED MODE R  R  . . . I  I   . . .
C                                          1  2        1  2
C ----------------------------------------------------------------
C
      NT2     = N * 2
      NT2M1   = NT2 - 1
      ABSINCX = IABS (INCX)
      ABSINCY = IABS (INCY)
      IXIY    = INCX * INCY
      INCXT2  = INCX * 2
      INCYT2  = INCY * 2
C
      RE = REAL (CA)
      AI = AIMAG (CA)
C
      IF (ABSINCY .EQ. 1) THEN
C
C        ..................
C        . ABS (INCY) = 1 .
C        ..................
C
         ASSIGN YD, Y(1,1; NT2)
C
         IF (IXIY .EQ. 1) THEN
C
C           .........................
C           . INCX = INCY = 1 OR -1 .
C           .........................
C
            ASSIGN XD, X(1,1; NT2)
C
         ELSE IF (IXIY .EQ. -1) THEN
C
C           ...........................
C           . INCX = - INCY = 1 OR -1 .
C           ...........................
C
C X(1,1;NT2) IS REVERSED INTO XD;
C REXD REFERENCES TO THE FIRST  ELEMENT OF XD;
C IMXD REFERENCES TO THE SECOND ELEMENT OF XD;
C
            ASSIGN XD, .DYN.NT2
            CALL Q8PACK (NT2M1, XD, REXD)
            CALL Q8ADDX (REXD, I64, IMXD)
C
            ASSIGN REP10D, REP10(1;NT2)
            CALL Q8VREVV (X'00',,X(1,1;NT2M1),,,REP10D,REXD)
            CALL Q8VREVV (X'00',,X(2,1;NT2M1),,,REP10D,IMXD)
C
         ELSE
C
            IF (IXIY .GT. 0) THEN
C
C              '''''''''''''''''''''''''''''''''''''''''''''''''''
C              ' INCX AND INCY BOTH POSITIVE OR BOTH NEGATIVE => '
C              ' GATHER REAL AND IMAGINARY PARTS OF X FORWARDS;  '
C              '''''''''''''''''''''''''''''''''''''''''''''''''''
C
               FIRSTX = 1
               IXT2   = IABS (INCXT2)
C
            ELSE
C
C              '''''''''''''''''''''''''''''''''''''''''''''''''''
C              ' INCX AND INCY ARE OF DIFFERENT SIGN =>          '
C              ' GATHER REAL AND IMAGINARY PARTS OF X BACKWARDS; '
C              '''''''''''''''''''''''''''''''''''''''''''''''''''
C
               FIRSTX = 1 + ABSINCX * (N-1)
               IXT2   = - IABS (INCXT2)
C
            END IF
C
C REXD REFERENCES TO THE REAL ELEMENTS OF X;
C IMXD REFERENCES TO THE IMAGINARY ELEMENTS OF X;
C
            ASSIGN REXD, .DYN.N
            ASSIGN IMXD, .DYN.N
C
            REXD = Q8VGATHP (X(1, FIRSTX; 1), IXT2, N; REXD)
            IMXD = Q8VGATHP (X(2, FIRSTX; 1), IXT2, N; IMXD)
C
C ABSINCY = 1 : USE MERGE
C XD REFERENCES TO THE MIXED REAL AND IMAGINARY ELEMENTS OF X;
C
            ASSIGN XD, .DYN.NT2
            ASSIGN REP10D, REP10 (1;NT2)
            XD = Q8VMERG (REXD, IMXD, REP10D; XD)
C
         END IF
C
         YD = YD + RE * XD
C
C THE REFERENCES OF REXD, REYD, IMXD AND IMYD ARE CHANGED !!!
C
C REXD REFERENCES TO THE FIRST/LAST REAL      ELEMENT OF X
C IMXD REFERENCES TO THE FIRST/LAST IMAGINARY ELEMENT OF X
C
C REYD REFERENCES TO THE FIRST REAL      ELEMENT OF THE NEW Y
C IMYD REFERENCES TO THE FIRST IMAGINARY ELEMENT OF THE NEW Y
C
C ALL VECTORS ARE OF LENGTH 2*N-1
C
         CALL Q8PACK (NT2M1, XD, REXD)
         CALL Q8ADDX (REXD, I64, IMXD)
C
         CALL Q8PACK (NT2M1, YD, REYD)
         CALL Q8ADDX (REYD, I64, IMYD)
C
         ASSIGN REP10D, REP10 (1;NT2M1)
         WHERE (REP10D)
            REYD = REYD - AI * IMXD
            IMYD = IMYD + AI * REXD
         END WHERE
C
         RETURN
C
      ELSE
C
C        ..................
C        . ABS (INCY) > 1 .
C        ..................
C
C THE VECTORS REFERENCED BY XD, REXD AND IMXD SHARE THE SAME STORAGE;
C XD BECOMES REXD FOLLOWED BY IMXD;
C
C REXD REFERENCES TO THE REAL ELEMENTS OF X;
C IMXD REFERENCES TO THE IMAGINARY ELEMENTS OF X;
C
         JSHIFT  = N * 64
         ASSIGN XD, .DYN.NT2
         CALL Q8PACK (N, XD, REXD)
         CALL Q8ADDX (REXD, JSHIFT, IMXD)
C
         IF (INCX .GT. 0) THEN
            FIRSTX  = 1
         ELSE
            FIRSTX  = 1 + ABSINCX * (N-1)
         END IF
C
         REXD = Q8VGATHP (X(1, FIRSTX; 1), INCXT2, N; REXD)
         IMXD = Q8VGATHP (X(2, FIRSTX; 1), INCXT2, N; IMXD)
C
C THE VECTORS REFERENCED BY YD, REYD AND IMYD SHARE THE SAME STORAGE;
C YD BECOMES REYD FOLLOWED BY IMYD;
C
C REYD REFERENCES TO THE REAL ELEMENTS OF Y;
C IMYD REFERENCES TO THE IMAGINARY ELEMENTS OF Y;
C
C        JSHIFT  = N * 64
         ASSIGN YD, .DYN.NT2
         CALL Q8PACK (N, YD, REYD)
         CALL Q8ADDX (REYD, JSHIFT, IMYD)
C
         IF (INCY .GT. 0) THEN
            FIRSTY  = 1
         ELSE
            FIRSTY  = 1 + ABSINCY * (N-1)
         END IF
C
         REYD = Q8VGATHP (Y(1, FIRSTY; 1), INCYT2, N; REYD)
         IMYD = Q8VGATHP (Y(2, FIRSTY; 1), INCYT2, N; IMYD)
C
C COMPUTE THE NEW Y
C
         YD = YD + RE * XD
C
         REYD = REYD - AI * IMXD
         IMYD = IMYD + AI * REXD
C
C        ABSINCY > 1 : USE SCATTER
C
         Y(1, FIRSTY; 1) = Q8VSCATP (REYD, INCYT2, N; Y(1, FIRSTY; 1))
         Y(2, FIRSTY; 1) = Q8VSCATP (IMYD, INCYT2, N; Y(2, FIRSTY; 1))
C
      END IF
C
      RETURN
      END
      SUBROUTINE CSCAL (N, CA, X, INCX)
      INTEGER           N, INCX
      REAL              X(2, *)
      COMPLEX           CA
C
C MARGREET LOUTER-NOOL, CWI, 1987 OCTOBER 9
C
C PURPOSE
C =======
C
C CSCAL PERFORMS A VECTOR SCALING OPERATION
C
C    X  := CA . X         FOR I = 1, ..., N
C     I          I
C
C WHERE X IS A COMPLEX VECTOR, CA IS A COMPLEX SCALAR;
C AND INCX >= 1
C
C ====================================================
C
      BIT         REP10D,
     +            REP10 (65536) /32768 * B'10'/
      DESCRIPTOR  REP10D
      INTEGER     INCXT2, NT2, NT2M1
      REAL        AI, RE
      INTEGER     I64 /64/, JSHIFT
      REAL        AITXD, IMAITXD, IMXD, REAITXD, REXD, XD
      DESCRIPTOR  AITXD, IMAITXD, IMXD, REAITXD, REXD, XD
      COMPLEX     ONE /(1.E0, 0.E0)/
C
C
      IF (N .LE. 0 .OR. CA .EQ. ONE) THEN
         RETURN
      ELSE IF (N .GT. 32767) THEN
         STOP 'BLAS ERROR (CSCAL) : N > 32767'
      ELSE IF (INCX .LT. 1) THEN
         STOP 'BLAS ERROR (CSCAL) : INCX < 1'
      END IF
C
C
      NT2    = N * 2
      NT2M1  = NT2 - 1
      INCXT2 = INCX * 2
C
      RE = REAL (CA)
      AI = AIMAG (CA)
C
C ---------------------------------------
C CSCAL PERFORMS THE NEXT OPERATIONS :
C
C AITX  = IM CA . X
C     I            I
C RE X  = RE CA . RE X  - IM AITX
C     I               I          I
C IM X  = RE CA . IM X  + RE AITX
C     I               I          I
C ---------------------------------------
C
      IF (INCX .EQ. 1) THEN
C
C        ............
C        . INCX = 1 .
C        ............
C
         ASSIGN XD, X(1,1; NT2)
         ASSIGN AITXD, .DYN.NT2
         AITXD = AI * XD
C
C THE VECTORS REFERENCED BY XD, REXD AND IMXD SHARE THE SAME STORAGE;
C REXD REFERS TO THE FIRST REAL      ELEMENT OF X;
C IMXD REFERS TO THE FIRST IMAGINARY ELEMENT OF X;
C
C THE VECTORS REFERENCED BY AITXD, REAITXD AND IMAITXD SHARE
C THE SAME STORAGE;
C REAITXD REFERS TO THE FIRST REAL      ELEMENT OF AITX;
C IMAITXD REFERS TO THE FIRST IMAGINARY ELEMENT OF AITX;
C
C ALL VECTORS ARE OF SIZE 2*N-1
C
         CALL Q8PACK (NT2M1, XD, REXD)
         CALL Q8ADDX (REXD, I64, IMXD)
C
         CALL Q8PACK (NT2M1, AITXD, REAITXD)
         CALL Q8ADDX (REAITXD, I64, IMAITXD)
C
         ASSIGN REP10D,  REP10 (1;NT2M1)
         WHERE (REP10D)
            REXD = RE * REXD - IMAITXD
            IMXD = RE * IMXD + REAITXD
         END WHERE
C
         RETURN
C
      ELSE
C
C        ............
C        . INCX > 1 .
C        ............
C
C THE VECTORS REFERENCED BY XD, REXD AND IMXD SHARE THE SAME STORAGE;
C XD BECOMES REXD FOLLOWED BY IMXD;
C
C REXD REFERENCES TO THE REAL ELEMENTS OF X;
C IMXD REFERENCES TO THE IMAGINARY ELEMENTS OF X;
C
         JSHIFT   = N * 64
         ASSIGN XD, .DYN.NT2
         CALL Q8PACK (N, XD, REXD)
         CALL Q8ADDX (REXD, JSHIFT, IMXD)
C
         REXD = Q8VGATHP (X(1,1; 1), INCXT2, N; REXD)
         IMXD = Q8VGATHP (X(2,1; 1), INCXT2, N; IMXD)
C
C THE VECTORS REFERENCED BY AITXD, REAITXD AND IMAITXD SHARE
C THE SAME STORAGE;
C AITXD BECOMES REAITXD FOLLOWED BY IMAITXD;
C
C REAITXD REFERENCES TO THE REAL ELEMENTS OF AITX;
C IMAITXD REFERENCES TO THE IMAGINARY ELEMENTS OF AITX;
C
         ASSIGN AITXD, .DYN.NT2
         AITXD = AI*XD
C
C        JSHIFT   = N * 64
         CALL Q8PACK (N, AITXD, REAITXD)
         CALL Q8ADDX (REAITXD, JSHIFT, IMAITXD)
C
C CALCULATE THE NEW REXD AND IMXD
C
         REXD = RE * REXD - IMAITXD
         IMXD = RE * IMXD + REAITXD
C
C RESTORE REXD AND IMXD IN X
C
         X(1,1; 1) = Q8VSCATP (REXD, INCXT2, N; X(1,1; 1))
         X(2,1; 1) = Q8VSCATP (IMXD, INCXT2, N; X(2,1; 1))
C
      END IF
C
      RETURN
      END
      SUBROUTINE CSSCAL (N, SA, X, INCX)
      INTEGER            N, INCX
      REAL               SA, X(2, *)
C
C MARGREET LOUTER-NOOL, CWI, 1987 OCTOBER 9
C
C PURPOSE
C =======
C
C CSSCAL PERFORMS A VECTOR SCALING OPERATION
C
C    X  := SA . X         FOR I = 1, ..., N
C     I          I
C
C WHERE X IS A COMPLEX VECTOR, SA IS A REAL SCALAR;
C AND INCX >= 1
C
C =================================================
C
      REAL        IMXD, REXD, XD
      DESCRIPTOR  IMXD, REXD, XD
      INTEGER     INCXT2, JSHIFT, NT2
      REAL        ONE /1.E0/
C
C
      IF (N .LE. 0 .OR. SA .EQ. ONE) THEN
         RETURN
      ELSE IF (N .GT. 32767) THEN
         STOP 'BLAS ERROR (CSSCAL) : N > 32767'
      ELSE IF (INCX .LT. 1) THEN
         STOP 'BLAS ERROR (CSSCAL) : INCX < 1'
      END IF
C
C
      NT2    = N * 2
      INCXT2 = INCX * 2
C
      IF (INCX .EQ. 1) THEN
C
C        ............
C        . INCX = 1 .
C        ............
C
         ASSIGN XD, X(1,1; NT2)
C
         XD = SA * XD
C
         RETURN
C
      ELSE
C
C        ............
C        . INCX > 1 .
C        ............
C
C THE VECTORS REFERENCED BY XD, REXD AND IMXD SHARE THE SAME STORAGE;
C XD BECOMES REXD FOLLOWED BY IMXD;
C
C REXD REFERENCES TO THE REAL ELEMENTS OF X;
C IMXD REFERENCES TO THE IMAGINARY ELEMENTS OF X;
C
         JSHIFT = N * 64
         ASSIGN XD, .DYN.NT2
         CALL Q8PACK (N, XD, REXD)
         CALL Q8ADDX (REXD, JSHIFT, IMXD)
C
         REXD = Q8VGATHP (X(1,1; 1), INCXT2, N; REXD)
         IMXD = Q8VGATHP (X(2,1; 1), INCXT2, N; IMXD)
C
C XD = SA * XD
C
         XD = SA * XD
C
C SCATTER REXD AND IMXD INTO X
C
         X(1,1; 1) = Q8VSCATP (REXD, INCXT2, N; X(1,1; 1))
         X(2,1; 1) = Q8VSCATP (IMXD, INCXT2, N; X(2,1; 1))
C
      END IF
C
      RETURN
      END
      SUBROUTINE CSWAP (N, X, INCX, Y, INCY)
      INTEGER           N, INCX, INCY
      REAL              X(2, *), Y(2, *)
C
C MARGREET LOUTER-NOOL, CWI, 1987 OCTOBER 9
C
C PURPOSE
C =======
C
C CSWAP INTERCHANGES THE VECTORS X AND Y
C
C    X :=: Y           FOR I = 1, ..., N
C     I     I
C
C WHERE X AND Y ARE COMPLEX VECTORS, AND INCX <> 0 AND INCY <> 0
C
C ==============================================================
C
      BIT         REP10D,
     +            REP10 (65536) /32768 * B'10'/
      DESCRIPTOR  REP10D
      INTEGER     ABSINCX, ABSINCY, FIRSTX, FIRSTY,
     +            INCXT2, INCYT2, IXIY, IXT2, IYT2,
     +            NT2, NT2M1
      INTEGER     I64 /64/
      REAL        IMXD, IMYD, REXD, REYD, TEMPD
      DESCRIPTOR  IMXD, IMYD, REXD, REYD, TEMPD
C
C
      IF (N .LE. 0) THEN
         RETURN
      ELSE IF (N .GT. 32767) THEN
         STOP 'BLAS ERROR (CSWAP) : N > 32767'
      ELSE IF (INCX .EQ. 0) THEN
         STOP 'BLAS ERROR (CSWAP) : INCX = 0'
      ELSE IF (INCY .EQ. 0) THEN
         STOP 'BLAS ERROR (CSWAP) : INCY = 0'
      END IF
C
C
      IXIY    = INCX * INCY
      NT2     = N * 2
C
      IF (IXIY .EQ. 1) THEN
C
C        .........................
C        . INCX = INCY = 1 OR -1 .
C        .........................
C
         ASSIGN TEMPD, .DYN.NT2
         TEMPD = X(1,1; NT2)
         X(1,1; NT2) = Y(1,1; NT2)
         Y(1,1; NT2) = TEMPD
C
         RETURN
C
      ELSE IF (IXIY .EQ. -1) THEN
C
C        ...........................
C        . INCX = - INCY = 1 OR -1 .
C        ...........................
C
C THE REAL AND IMAGINARY PART OF X ARE REVERSED INTO REXD;
C IMXD REFERENCES TO THE SECOND VECTOR ELEMENT OF REXD;
C THE REAL AND IMAGINARY PART OF Y ARE REVERSED INTO X;
C REXD IS COPIED INTO Y.
C
         NT2M1 = NT2 - 1
         ASSIGN REXD, .DYN.NT2
         CALL Q8PACK (NT2M1, REXD, IMXD)
         CALL Q8ADDX (IMXD, I64, IMXD)
C
         ASSIGN REP10D, REP10(1;NT2M1)
         CALL Q8VREVV (X'00',,X(1,1;NT2M1),,,REP10D,REXD)
         CALL Q8VREVV (X'00',,X(2,1;NT2M1),,,REP10D,IMXD)
         CALL Q8VREVV (X'00',,Y(1,1;NT2M1),,,REP10D,X(1,1;NT2M1))
         CALL Q8VREVV (X'00',,Y(2,1;NT2M1),,,REP10D,X(2,1;NT2M1))
C
         Y(1,1;NT2) = REXD
C
         RETURN
C
      ELSE
C
C        .............................................
C        . NOT ( ABS (INCX) = 1 AND ABS (INCY) = 1 ) .
C        .............................................
C
         ABSINCX = IABS (INCX)
         ABSINCY = IABS (INCY)
         INCXT2  = INCX * 2
         INCYT2  = INCY * 2
C
         IF (IXIY .GT. 0) THEN
C
C           '''''''''''''''''''''''''''''''''''''''''''''''''''
C           ' INCX AND INCY BOTH POSITIVE OR BOTH NEGATIVE => '
C           ' GATHER X FORWARDS INTO XD;                      '
C           ' GATHER Y FORWARDS INTO YD;                      '
C           ' SCATTER XD FORWARDS INTO Y;                     '
C           ' SCATTER YD FORWARDS INTO X;                     '
C           '''''''''''''''''''''''''''''''''''''''''''''''''''
C
            FIRSTX = 1
            FIRSTY = 1
            IXT2   = IABS (INCXT2)
            IYT2   = IABS (INCYT2)
C
         ELSE
C
C           '''''''''''''''''''''''''''''''''''
C           ' INCX AND INCY OF DIFFERENT SIGN '
C           ' GATHER X BACKWARDS INTO XD;     '
C           ' GATHER Y BACKWARDS INTO YD;     '
C           ' SCATTER XD FORWARDS INTO Y;     '
C           ' SCATTER YD FORWARDS INTO X;     '
C           '''''''''''''''''''''''''''''''''''
C
            FIRSTX = 1 + ABSINCX * (N-1)
            FIRSTY = 1 + ABSINCY * (N-1)
            IXT2   = - IABS (INCXT2)
            IYT2   = - IABS (INCYT2)
C
         END IF
C
C REXD REFERENCES TO THE REAL ELEMENTS OF X;
C IMXD REFERENCES TO THE IMAGINARY ELEMENTS OF X;
C
         ASSIGN REXD, .DYN.N
         ASSIGN IMXD, .DYN.N
C
         REXD = Q8VGATHP (X(1,FIRSTX; 1), IXT2, N; REXD)
         IMXD = Q8VGATHP (X(2,FIRSTX; 1), IXT2, N; IMXD)
C
C REYD REFERENCES TO THE REAL ELEMENTS OF Y;
C IMYD REFERENCES TO THE IMAGINARY ELEMENTS OF Y;
C
         ASSIGN REYD, .DYN.N
         ASSIGN IMYD, .DYN.N
C
         REYD = Q8VGATHP (Y(1,FIRSTY; 1), IYT2, N; REYD)
         IMYD = Q8VGATHP (Y(2,FIRSTY; 1), IYT2, N; IMYD)
C
C ABSINCX = 1 (ABSINCY = 1) : USE MERGE
C ABSINCX > 1 (ABSINCY > 1) : USE SCATTER
C
         IF (ABSINCX .EQ. 1) THEN
            ASSIGN REP10D, REP10(1;NT2)
            X(1,1; 1) = Q8VMERG (REYD, IMYD, REP10D; X(1,1; 1))
         ELSE
            IXT2 = IABS (INCXT2)
            X(1,1; 1) = Q8VSCATP (REYD, IXT2, N; X(1,1; 1))
            X(2,1; 1) = Q8VSCATP (IMYD, IXT2, N; X(2,1; 1))
         END IF
C
         IF (ABSINCY .EQ. 1) THEN
            ASSIGN REP10D, REP10(1;NT2)
            Y(1,1; 1) = Q8VMERG (REXD, IMXD, REP10D; Y(1,1; 1))
         ELSE
            IYT2 = IABS (INCYT2)
            Y(1,1; 1) = Q8VSCATP (REXD, IYT2, N; Y(1,1; 1))
            Y(2,1; 1) = Q8VSCATP (IMXD, IYT2, N; Y(2,1; 1))
         END IF
C
      END IF
C
      RETURN
      END
      SUBROUTINE CCOPY (N, X, INCX, Y, INCY)
      INTEGER           N, INCX, INCY
      REAL              X(2, *), Y(2, *)
C
C MARGREET LOUTER-NOOL, CWI, 1987 OCTOBER 9
C
C PURPOSE
C =======
C
C CCOPY COPIES A VECTOR X TO Y
C
C    Y := X           FOR I = 1, ..., N
C     I    I
C
C WHERE X AND Y ARE COMPLEX VECTORS, AND INCX <> 0 AND INCY <> 0
C
C ==============================================================
C
      BIT         REP10D,
     +            REP10 (65536) /32768 * B'10'/
      DESCRIPTOR  REP10D
      INTEGER     FIRSTX, IXIY, IXT2, IYT2, NT2, NT2M1
      REAL        IMXD, REXD
      DESCRIPTOR  IMXD, REXD
C
C
      IF (N .LE. 0) THEN
         RETURN
      ELSE IF (N .GT. 32767) THEN
         STOP 'BLAS ERROR (CCOPY) : N > 32767'
      ELSE IF (INCX .EQ. 0) THEN
         STOP 'BLAS ERROR (CCOPY) : INCX = 0'
      ELSE IF (INCY .EQ. 0) THEN
         STOP 'BLAS ERROR (CCOPY) : INCY = 0'
      END IF
C
C
      IXIY    = INCX * INCY
      NT2     = N * 2
C
      IF (IXIY .EQ. 1) THEN
C
C        .........................
C        . INCX = INCY = 1 OR -1 .
C        .........................
C
         Y(1,1; NT2) = X(1,1; NT2)
C
         RETURN
C
      ELSE IF (IXIY .EQ. -1) THEN
C
C        ...........................
C        . INCX = - INCY = 1 OR -1 .
C        ...........................
C
         NT2M1 = NT2 - 1
         ASSIGN REP10D, REP10(1;NT2M1)
         CALL Q8VREVV (X'00',,X(1,1;NT2M1),,,REP10D,Y(1,1;NT2M1))
         CALL Q8VREVV (X'00',,X(2,1;NT2M1),,,REP10D,Y(2,1;NT2M1))
C
         RETURN
C
      ELSE
C
C        ...............................
C        . NOT (INCX = INCY = 1 OR -1) .
C        ...............................
C
         IF (IXIY .GT. 0) THEN
C
C           '''''''''''''''''''''''''''''''''''''''''''''''''''
C           ' INCX AND INCY BOTH POSITIVE OR BOTH NEGATIVE => '
C           ' GATHER X FORWARDS INTO REXD AND IMXD;           '
C           ' SCATTER REXD AND IMXD FORWARDS INTO Y           '
C           '''''''''''''''''''''''''''''''''''''''''''''''''''
C
            FIRSTX = 1
            IXT2   = IABS (INCX * 2)
C
         ELSE
C
C           ''''''''''''''''''''''''''''''''''''''''''
C           ' INCX AND INCY OF DIFFERENT SIGN =>     '
C           ' GATHER X BACKWARDS INTO REXD AND IMXD  '
C           ' SCATTER REXD AND IMXD FORWARDS INTO Y  '
C           ''''''''''''''''''''''''''''''''''''''''''
C
            FIRSTX = 1 + IABS (INCX) * (N-1)
            IXT2   = - IABS (INCX * 2)
C
         END IF
C
C REXD REFERENCES TO THE REAL ELEMENTS OF X;
C IMXD REFERENCES TO THE IMAGINARY ELEMENTS OF X;
C
         ASSIGN REXD, .DYN.N
         ASSIGN IMXD, .DYN.N
C
         REXD = Q8VGATHP (X(1,FIRSTX; 1), IXT2, N; REXD)
         IMXD = Q8VGATHP (X(2,FIRSTX; 1), IXT2, N; IMXD)
C
C ABSINCY = 1 : USE MERGE
C ABSINCY > 1 : USE SCATTER
C
         IF (IABS (INCY) .EQ. 1) THEN
            ASSIGN REP10D, REP10(1;NT2)
            Y(1,1; 1) = Q8VMERG (REXD, IMXD, REP10D; Y(1,1; 1))
         ELSE
            IYT2 = IABS (INCY*2)
            Y(1,1; 1) = Q8VSCATP (REXD, IYT2, N; Y(1,1; 1))
            Y(2,1; 1) = Q8VSCATP (IMXD, IYT2, N; Y(2,1; 1))
         END IF
C
      END IF
C
      RETURN
      END
      SUBROUTINE CSROT (N, X, INCX, Y, INCY, SC, SS)
      INTEGER            N, INCX, INCY
      REAL               X(2, *), Y(2, *), SC, SS
C
C MARGREET LOUTER-NOOL, CWI, 1987 OCTOBER 9
C
C PURPOSE
C =======
C
C CSROT PERFORMS THE FOLLOWING OPERATION
C
C    X := SC * X  + SS * Y      FOR I = 1, ..., N
C     I         I         I
C    Y :=-SS . X  + SC * Y      FOR I = 1, ..., N
C     I         I         I
C
C WHERE X AND Y ARE COMPLEX VECTORS, SC AND SS ARE REAL SCALARS;
C INCX <> 0 AND INCY <> 0.
C
C ==============================================================
C
      BIT         REP10D,
     +            REP10 (65536) /32768 * B'10'/
      DESCRIPTOR  REP10D
      INTEGER     ABSINCX, ABSINCY, FIRSTX, FIRSTY,
     +            IXIY, IXT2, IYT2,
     +            NT2, NT2M1
      INTEGER     I64 /64/, JSHIFT
      REAL        IMXD, IMYD, REXD, REYD,
     +            XD, YD, ZD
      DESCRIPTOR  IMXD, IMYD, REXD, REYD,
     +            XD, YD, ZD
      REAL       ZERO /0.E0/, ONE /1.E0/
C
C
      IF (N .LE. 0 .OR. (SS .EQ. ZERO. AND. SC .EQ. ONE)) THEN
         RETURN
      ELSE IF (N .GT. 32767) THEN
         STOP 'BLAS ERROR (CSROT) : N > 32767'
      ELSE IF (INCX .EQ. 0) THEN
         STOP 'BLAS ERROR (CSROT) : INCX = 0'
      ELSE IF (INCY .EQ. 0) THEN
         STOP 'BLAS ERROR (CSROT) : INCY = 0'
      END IF
C
C ----------------------------------------------------
C
C CSROT PERFORMS THE NEXT OPERATIONS :
C
C   Z  = SS * (X + Y )
C    I          I   I
C
C   X  = (SC - SS) * X  + Z
C    I                I    I
C   Y  = (SC + SS) * Y  - Z
C    I                I    I
C
C FOR I = 1 , ... , N
C -------------------------------------
C
      IXIY    = INCX * INCY
      NT2     = N * 2
C
      IF (IXIY .EQ. 1) THEN
C
C        .........................
C        . INCX = INCY = 1 OR -1 .
C        .........................
C
         ASSIGN XD, X(1,1; NT2)
         ASSIGN YD, Y(1,1; NT2)
C
      ELSE IF (IXIY .EQ. -1) THEN
C
C        ...........................
C        . INCX = - INCY = 1 OR -1 .
C        ...........................
C
C XD BECOMES THE REVERSE OF X(1,1;NT2)
C
         NT2M1 = NT2 - 1
         ASSIGN XD, .DYN.NT2
         CALL Q8PACK (NT2M1, XD, REXD)
         CALL Q8ADDX (REXD, I64, IMXD)
C
         ASSIGN REP10D, REP10(1;NT2M1)
         CALL Q8VREVV (X'00',,X(1,1;NT2M1),,,REP10D,REXD)
         CALL Q8VREVV (X'00',,X(2,1;NT2M1),,,REP10D,IMXD)
C
         ASSIGN YD, Y(1,1; NT2)
C
      ELSE
C
C THE VECTORS REFERENCED BY XD, REXD AND IMXD SHARE THE SAME STORAGE;
C XD BECOMES REXD FOLLOWED BY IMXD;
C
C REXD REFERENCES TO THE REAL ELEMENTS OF X;
C IMXD REFERENCES TO THE IMAGINARY ELEMENTS OF X;
C
C THE VECTORS REFERENCED BY YD, REYD AND IMYD SHARE THE SAME STORAGE;
C YD BECOMES REYD FOLLOWED BY IMYD;
C
C REYD REFERENCES TO THE REAL ELEMENTS OF Y;
C IMYD REFERENCES TO THE IMAGINARY ELEMENTS OF Y;
C
         JSHIFT    = N * 64
         ASSIGN XD, .DYN.NT2
         CALL Q8PACK (N, XD, REXD)
         CALL Q8ADDX (REXD, JSHIFT, IMXD)
         ASSIGN YD, .DYN.NT2
         CALL Q8PACK (N, YD, REYD)
         CALL Q8ADDX (REYD, JSHIFT, IMYD)
C
         ABSINCX = IABS (INCX)
         ABSINCY = IABS (INCY)
C
         IF (ABSINCX .EQ. 1) THEN
C
C           ..................
C           . ABS (INCX) = 1 .
C           ..................
C           ''''''''''''''''''''''
C           ' GATHER X FORWARDS; '
C           ''''''''''''''''''''''
C
            FIRSTX = 1
            IXT2   = 2
C
            IF (IXIY .GT. 0) THEN
C
C              '''''''''''''''''''''''''''''
C              ' GATHER Y FORWARDS INTO YD '
C              '''''''''''''''''''''''''''''
C
               FIRSTY = 1
               IYT2   = ABSINCY * 2
C
            ELSE
C
C              ''''''''''''''''''''''''''''''
C              ' GATHER Y BACKWARDS INTO YD '
C              ''''''''''''''''''''''''''''''
C
               FIRSTY = 1 + ABSINCY * (N - 1)
               IYT2   = - ABSINCY * 2
C
            END IF
C
         ELSE
C
C           ''''''''''''''''''''''
C           ' GATHER Y FORWARDS; '
C           ''''''''''''''''''''''
C
            FIRSTY = 1
            IYT2   = ABSINCY * 2
C
            IF (IXIY .GT. 0) THEN
C
C              '''''''''''''''''''''''''''''
C              ' GATHER X FORWARDS INTO XD '
C              '''''''''''''''''''''''''''''
C
               FIRSTX = 1
               IXT2   = ABSINCX * 2
            ELSE
C
C              ''''''''''''''''''''''''''''''
C              ' GATHER X BACKWARDS INTO XD '
C              ''''''''''''''''''''''''''''''
C
               FIRSTX = 1 + ABSINCX * (N - 1)
               IXT2   = - ABSINCX * 2
            END IF
C
         END IF
C
         REXD = Q8VGATHP (X(1,FIRSTX; 1), IXT2, N; REXD)
         IMXD = Q8VGATHP (X(2,FIRSTX; 1), IXT2, N; IMXD)
         REYD = Q8VGATHP (Y(1,FIRSTY; 1), IYT2, N; REYD)
         IMYD = Q8VGATHP (Y(2,FIRSTY; 1), IYT2, N; IMYD)
C
      END IF
C
C
C COMPUTE ZD = SS * (XD + YD)
C
      ASSIGN ZD, .DYN.NT2
      ZD = SS * (XD + YD)
C
C COMPUTE THE NEW XD AND YD
C
      XD = (SC - SS) * XD + ZD
      YD = (SC + SS) * YD - ZD
C
      IF (IXIY .EQ. 1) THEN
C
         RETURN
C
      ELSE IF (IXIY .EQ. -1) THEN
C
C REVERSE XD INTO X(1,1;NT2)
C
         CALL Q8PACK (NT2M1, XD, REXD)
         CALL Q8ADDX (REXD, I64, IMXD)
C
         CALL Q8VREVV (X'00',,REXD,,,REP10D,X(1,1;NT2M1))
         CALL Q8VREVV (X'00',,IMXD,,,REP10D,X(2,1;NT2M1))
C
         RETURN
C
      ELSE
C
C        ....................................
C        . ABS (INCX) > 1 OR ABS (INCY) > 1 .
C        ....................................
C
C        ABS (INCX) = 1  (ABS (INCY) = 1) : USE MERGE
C        ABS (INCX) > 1  (ABS (INCY) > 1) : USE SCATTER
C
         IF (ABSINCX .EQ. 1) THEN
            ASSIGN REP10D, REP10(1;NT2)
            X(1,1; 1) = Q8VMERG (REXD, IMXD, REP10D; X(1,1; 1))
         ELSE
            X(1,FIRSTX; 1) = Q8VSCATP (REXD, IXT2, N; X(1,FIRSTX; 1))
            X(2,FIRSTX; 1) = Q8VSCATP (IMXD, IXT2, N; X(2,FIRSTX; 1))
         END IF
C
         IF (ABSINCY .EQ. 1) THEN
            ASSIGN REP10D, REP10(1;NT2)
            Y(1,1; 1) = Q8VMERG (REYD, IMYD, REP10D; Y(1,1; 1))
         ELSE
            Y(1,FIRSTY; 1) = Q8VSCATP (REYD, IYT2, N; Y(1,FIRSTY; 1))
            Y(2,FIRSTY; 1) = Q8VSCATP (IMYD, IYT2, N; Y(2,FIRSTY; 1))
         END IF
C
      END IF
C
      RETURN
      END
      SUBROUTINE CSROT4 (N, X, INCX, Y, INCY, SC, SS)
      INTEGER            N, INCX, INCY
      REAL               X(2, *), Y(2, *), SC, SS
C
C MARGREET LOUTER-NOOL, CWI, 1987 OCTOBER 9
C
C PURPOSE
C =======
C
C CSROT PERFORMS THE FOLLOWING OPERATION
C
C    X := SC * X  + SS * Y      FOR I = 1, ..., N
C     I         I         I
C    Y :=-SS . X  + SC * Y      FOR I = 1, ..., N
C     I         I         I
C
C WHERE X AND Y ARE COMPLEX VECTORS, SC AND SS ARE REAL SCALARS;
C INCX <> 0 AND INCY <> 0.
C
C ==============================================================
C
      BIT         REP10D,
     +            REP10 (65536) /32768 * B'10'/
      DESCRIPTOR  REP10D
      INTEGER     ABSINCX, ABSINCY, FIRSTX, FIRSTY,
     +            IXT2, IYT2,
     +            NT2, NT2M1
      INTEGER     I64 /64/, JSHIFT
      REAL        IMXD, IMYD, REXD, REYD,
     +            VD, WD, XD, YD
      DESCRIPTOR  IMXD, IMYD, REXD, REYD,
     +            VD, WD, XD, YD
      REAL       ZERO /0.E0/, ONE /1.E0/
C
C
      IF (N .LE. 0 .OR. (SS .EQ. ZERO. AND. SC .EQ. ONE)) THEN
         RETURN
      ELSE IF (N .GT. 32767) THEN
         STOP 'BLAS ERROR (CSROT) : N > 32767'
      ELSE IF (INCX .EQ. 0) THEN
         STOP 'BLAS ERROR (CSROT) : INCX = 0'
      ELSE IF (INCY .EQ. 0) THEN
         STOP 'BLAS ERROR (CSROT) : INCY = 0'
      END IF
C
C ----------------------------------------------------
C
C CSROT4 PERFORMS THE NEXT OPERATIONS :
C
C   V  = SC * X
C    I         I
C   W  = SS * X
C    I         I
C
C   X  = V  + SS * Y
C    I              I
C   Y  = SC * Y  - W
C    I         I    I
C
C FOR I = 1 , ... , N
C -------------------------------------
C
      IXIY    = INCX * INCY
      NT2     = N * 2
C
      IF (IXIY .EQ. 1) THEN
C
C        .........................
C        . INCX = INCY = 1 OR -1 .
C        .........................
C
         ASSIGN XD, X(1,1; NT2)
         ASSIGN YD, Y(1,1; NT2)
C
      ELSE IF (IXIY .EQ. -1) THEN
C
C        ...........................
C        . INCX = - INCY = 1 OR -1 .
C        ...........................
C
C XD BECOMES THE REVERSE OF X(1,1;NT2)
C
         NT2M1 = NT2 - 1
         ASSIGN XD, .DYN.NT2
         CALL Q8PACK (NT2M1, XD, REXD)
         CALL Q8ADDX (REXD, I64, IMXD)
C
         ASSIGN REP10D, REP10(1;NT2M1)
         CALL Q8VREVV (X'00',,X(1,1;NT2M1),,,REP10D,REXD)
         CALL Q8VREVV (X'00',,X(2,1;NT2M1),,,REP10D,IMXD)
C
         ASSIGN YD, Y(1,1; NT2)
C
      ELSE
C
C THE VECTORS REFERENCED BY XD, REXD AND IMXD SHARE THE SAME STORAGE;
C XD BECOMES REXD FOLLOWED BY IMXD;
C
C REXD REFERENCES TO THE REAL ELEMENTS OF X;
C IMXD REFERENCES TO THE IMAGINARY ELEMENTS OF X;
C
C THE VECTORS REFERENCED BY YD, REYD AND IMYD SHARE THE SAME STORAGE;
C YD BECOMES REYD FOLLOWED BY IMYD;
C
C REYD REFERENCES TO THE REAL ELEMENTS OF Y;
C IMYD REFERENCES TO THE IMAGINARY ELEMENTS OF Y;
C
         JSHIFT = N * 64
         ASSIGN XD, .DYN.NT2
         CALL Q8PACK (N, XD, REXD)
         CALL Q8ADDX (REXD, JSHIFT, IMXD)
         ASSIGN YD, .DYN.NT2
         CALL Q8PACK (N, YD, REYD)
         CALL Q8ADDX (REYD, JSHIFT, IMYD)
C
         ABSINCX = IABS (INCX)
         ABSINCY = IABS (INCY)
C
         IF (ABSINCX .EQ. 1) THEN
C
C           ..................
C           . ABS (INCX) = 1 .
C           ..................
C           ''''''''''''''''''''''
C           ' GATHER X FORWARDS; '
C           ''''''''''''''''''''''
C
            FIRSTX = 1
            IXT2   = 2
C
            IF (IXIY .GT. 0) THEN
C
C              '''''''''''''''''''''''''''''
C              ' GATHER Y FORWARDS INTO YD '
C              '''''''''''''''''''''''''''''
C
               FIRSTY = 1
               IYT2   = ABSINCY * 2
C
            ELSE
C
C              ''''''''''''''''''''''''''''''
C              ' GATHER Y BACKWARDS INTO YD '
C              ''''''''''''''''''''''''''''''
C
               FIRSTY = 1 + ABSINCY * (N - 1)
               IYT2   = - ABSINCY * 2
C
            END IF
C
         ELSE
C
C           ''''''''''''''''''''''
C           ' GATHER Y FORWARDS; '
C           ''''''''''''''''''''''
C
            FIRSTY = 1
            IYT2   = ABSINCY * 2
C
            IF (IXIY .GT. 0) THEN
C
C              '''''''''''''''''''''''''''''
C              ' GATHER X FORWARDS INTO XD '
C              '''''''''''''''''''''''''''''
C
               FIRSTX = 1
               IXT2   = ABSINCX * 2
C
            ELSE
C
C              ''''''''''''''''''''''''''''''
C              ' GATHER X BACKWARDS INTO XD '
C              ''''''''''''''''''''''''''''''
C
               FIRSTX = 1 + ABSINCX * (N - 1)
               IXT2   = - ABSINCX * 2
C
            END IF
C
         END IF
C
         REXD = Q8VGATHP (X(1,FIRSTX; 1), IXT2, N; REXD)
         IMXD = Q8VGATHP (X(2,FIRSTX; 1), IXT2, N; IMXD)
         REYD = Q8VGATHP (Y(1,FIRSTY; 1), IYT2, N; REYD)
         IMYD = Q8VGATHP (Y(2,FIRSTY; 1), IYT2, N; IMYD)
C
      END IF
C
C
C COMPUTE VD = SC * XD AND WD = SS * XD
C
      ASSIGN VD, .DYN.NT2
      ASSIGN WD, .DYN.NT2
      VD = SC * XD
      WD = SS * XD
C
C COMPUTE THE NEW XD AND YD
C
      XD = VD + SS * YD
      YD = SC * YD - WD
C
      IF (IXIY .EQ. 1) THEN
C
         RETURN
C
      ELSE IF (IXIY .EQ. -1) THEN
C
C REVERSE XD INTO X(1,1;NT2)
C
         CALL Q8PACK (NT2M1, XD, REXD)
         CALL Q8ADDX (REXD, I64, IMXD)
C
         CALL Q8VREVV (X'00',,REXD,,,REP10D,X(1,1;NT2M1))
         CALL Q8VREVV (X'00',,IMXD,,,REP10D,X(2,1;NT2M1))
C
         RETURN
C
      ELSE
C
C        ....................................
C        . ABS (INCX) > 1 OR ABS (INCY) > 1 .
C        ....................................
C
C        ABS (INCX) =  1 (ABS (INCY) = 1) : USE MERGE
C        ABS (INCX) >  1 (ABS (INCY) > 1) : USE SCATTER
C
         IF (ABSINCX .EQ. 1) THEN
            ASSIGN REP10D, REP10(1;NT2)
            X(1,1; 1) = Q8VMERG (REXD, IMXD, REP10D; X(1,1; 1))
         ELSE
            X(1,FIRSTX; 1) = Q8VSCATP (REXD, IXT2, N; X(1,FIRSTX; 1))
            X(2,FIRSTX; 1) = Q8VSCATP (IMXD, IXT2, N; X(2,FIRSTX; 1))
         END IF
C
         IF (ABSINCY .EQ. 1) THEN
            ASSIGN REP10D, REP10(1;NT2)
            Y(1,1; 1) = Q8VMERG (REYD, IMYD, REP10D; Y(1,1; 1))
         ELSE
            Y(1,FIRSTY; 1) = Q8VSCATP (REYD, IYT2, N; Y(1,FIRSTY; 1))
            Y(2,FIRSTY; 1) = Q8VSCATP (IMYD, IYT2, N; Y(2,FIRSTY; 1))
         END IF
C
      END IF
C
      RETURN
      END