source/unres/src_MD-M/blas.f

   1 C 10 NOV 94 - MWS - DNRM2: REMOVE FTNCHECK WARNINGS
   2 C 11 JUN 94 - MWS - INCLUDE A COPY OF DGEMV (LEVEL TWO ROUTINE)
   3 C 11 AUG 87 - MWS - SANITIZE FLOATING POINT CONSTANTS IN DNRM2
   4 C 26 MAR 87 - MWS - USE GENERIC SIGN IN DROTG
   5 C 28 NOV 86 - STE - SUPPLY ALL LEVEL ONE BLAS
   6 C  7 JUL 86 - JAB - SANITIZE FLOATING POINT CONSTANTS
   7 C
   8 C BASIC LINEAR ALGEBRA SUBPROGRAMS (BLAS) FROM LINPACK  (LEVEL 1)
   9 C
  10 C   THIS MODULE SHOULD BE COMPILED ONLY IF SPECIALLY CODED
  11 C   VERSIONS OF THESE ROUTINES ARE NOT AVAILABLE ON THE TARGET MACHINE
  12 C
  13 C*MODULE BLAS1   *DECK DASUM
  14       DOUBLE PRECISION FUNCTION DASUM(N,DX,INCX)
  15 C
  16 C     TAKES THE SUM OF THE ABSOLUTE VALUES.
  17 C     JACK DONGARRA, LINPACK, 3/11/78.
  18 C
  19       DOUBLE PRECISION DX(1),DTEMP
  20       INTEGER I,INCX,M,MP1,N,NINCX
  21 C
  22       DASUM = 0.0D+00
  23       DTEMP = 0.0D+00
  24       IF(N.LE.0)RETURN
  25       IF(INCX.EQ.1)GO TO 20
  26 C
  27 C        CODE FOR INCREMENT NOT EQUAL TO 1
  28 C
  29       NINCX = N*INCX
  30       DO 10 I = 1,NINCX,INCX
  31         DTEMP = DTEMP + ABS(DX(I))
  32    10 CONTINUE
  33       DASUM = DTEMP
  34       RETURN
  35 C
  36 C        CODE FOR INCREMENT EQUAL TO 1
  37 C
  38 C
  39 C        CLEAN-UP LOOP
  40 C
  41    20 M = MOD(N,6)
  42       IF( M .EQ. 0 ) GO TO 40
  43       DO 30 I = 1,M
  44         DTEMP = DTEMP + ABS(DX(I))
  45    30 CONTINUE
  46       IF( N .LT. 6 ) GO TO 60
  47    40 MP1 = M + 1
  48       DO 50 I = MP1,N,6
  49         DTEMP = DTEMP + ABS(DX(I)) + ABS(DX(I + 1)) + ABS(DX(I + 2))
  50      *  + ABS(DX(I + 3)) + ABS(DX(I + 4)) + ABS(DX(I + 5))
  51    50 CONTINUE
  52    60 DASUM = DTEMP
  53       RETURN
  54       END
  55 C*MODULE BLAS1   *DECK DAXPY
  56       SUBROUTINE DAXPY(N,DA,DX,INCX,DY,INCY)
  57       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
  58       DIMENSION DX(1),DY(1)
  59 C
  60 C     CONSTANT TIMES A VECTOR PLUS A VECTOR.
  61 C           DY(I) = DY(I) + DA * DX(I)
  62 C     USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE.
  63 C     JACK DONGARRA, LINPACK, 3/11/78.
  64 C
  65       IF(N.LE.0)RETURN
  66       IF (DA .EQ. 0.0D+00) RETURN
  67       IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20
  68 C
  69 C        CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS
  70 C          NOT EQUAL TO 1
  71 C
  72       IX = 1
  73       IY = 1
  74       IF(INCX.LT.0)IX = (-N+1)*INCX + 1
  75       IF(INCY.LT.0)IY = (-N+1)*INCY + 1
  76       DO 10 I = 1,N
  77         DY(IY) = DY(IY) + DA*DX(IX)
  78         IX = IX + INCX
  79         IY = IY + INCY
  80    10 CONTINUE
  81       RETURN
  82 C
  83 C        CODE FOR BOTH INCREMENTS EQUAL TO 1
  84 C
  85 C
  86 C        CLEAN-UP LOOP
  87 C
  88    20 M = MOD(N,4)
  89       IF( M .EQ. 0 ) GO TO 40
  90       DO 30 I = 1,M
  91         DY(I) = DY(I) + DA*DX(I)
  92    30 CONTINUE
  93       IF( N .LT. 4 ) RETURN
  94    40 MP1 = M + 1
  95       DO 50 I = MP1,N,4
  96         DY(I) = DY(I) + DA*DX(I)
  97         DY(I + 1) = DY(I + 1) + DA*DX(I + 1)
  98         DY(I + 2) = DY(I + 2) + DA*DX(I + 2)
  99         DY(I + 3) = DY(I + 3) + DA*DX(I + 3)
 100    50 CONTINUE
 101       RETURN
 102       END
 103 C*MODULE BLAS1   *DECK DCOPY
 104       SUBROUTINE  DCOPY(N,DX,INCX,DY,INCY)
 105       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
 106       DIMENSION DX(*),DY(*)
 107 C
 108 C     COPIES A VECTOR.
 109 C           DY(I) <== DX(I)
 110 C     USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE.
 111 C     JACK DONGARRA, LINPACK, 3/11/78.
 112 C
 113       IF(N.LE.0)RETURN
 114       IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20
 115 C
 116 C        CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS
 117 C          NOT EQUAL TO 1
 118 C
 119       IX = 1
 120       IY = 1
 121       IF(INCX.LT.0)IX = (-N+1)*INCX + 1
 122       IF(INCY.LT.0)IY = (-N+1)*INCY + 1
 123       DO 10 I = 1,N
 124         DY(IY) = DX(IX)
 125         IX = IX + INCX
 126         IY = IY + INCY
 127    10 CONTINUE
 128       RETURN
 129 C
 130 C        CODE FOR BOTH INCREMENTS EQUAL TO 1
 131 C
 132 C
 133 C        CLEAN-UP LOOP
 134 C
 135    20 M = MOD(N,7)
 136       IF( M .EQ. 0 ) GO TO 40
 137       DO 30 I = 1,M
 138         DY(I) = DX(I)
 139    30 CONTINUE
 140       IF( N .LT. 7 ) RETURN
 141    40 MP1 = M + 1
 142       DO 50 I = MP1,N,7
 143         DY(I) = DX(I)
 144         DY(I + 1) = DX(I + 1)
 145         DY(I + 2) = DX(I + 2)
 146         DY(I + 3) = DX(I + 3)
 147         DY(I + 4) = DX(I + 4)
 148         DY(I + 5) = DX(I + 5)
 149         DY(I + 6) = DX(I + 6)
 150    50 CONTINUE
 151       RETURN
 152       END
 153 C*MODULE BLAS1   *DECK DDOT
 154       DOUBLE PRECISION FUNCTION DDOT(N,DX,INCX,DY,INCY)
 155       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
 156       DIMENSION DX(1),DY(1)
 157 C
 158 C     FORMS THE DOT PRODUCT OF TWO VECTORS.
 159 C           DOT = DX(I) * DY(I)
 160 C     USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE.
 161 C     JACK DONGARRA, LINPACK, 3/11/78.
 162 C
 163       DDOT = 0.0D+00
 164       DTEMP = 0.0D+00
 165       IF(N.LE.0)RETURN
 166       IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20
 167 C
 168 C        CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS
 169 C          NOT EQUAL TO 1
 170 C
 171       IX = 1
 172       IY = 1
 173       IF(INCX.LT.0)IX = (-N+1)*INCX + 1
 174       IF(INCY.LT.0)IY = (-N+1)*INCY + 1
 175       DO 10 I = 1,N
 176         DTEMP = DTEMP + DX(IX)*DY(IY)
 177         IX = IX + INCX
 178         IY = IY + INCY
 179    10 CONTINUE
 180       DDOT = DTEMP
 181       RETURN
 182 C
 183 C        CODE FOR BOTH INCREMENTS EQUAL TO 1
 184 C
 185 C
 186 C        CLEAN-UP LOOP
 187 C
 188    20 M = MOD(N,5)
 189       IF( M .EQ. 0 ) GO TO 40
 190       DO 30 I = 1,M
 191         DTEMP = DTEMP + DX(I)*DY(I)
 192    30 CONTINUE
 193       IF( N .LT. 5 ) GO TO 60
 194    40 MP1 = M + 1
 195       DO 50 I = MP1,N,5
 196         DTEMP = DTEMP + DX(I)*DY(I) + DX(I + 1)*DY(I + 1) +
 197      *   DX(I + 2)*DY(I + 2) + DX(I + 3)*DY(I + 3) + DX(I + 4)*DY(I + 4)
 198    50 CONTINUE
 199    60 DDOT = DTEMP
 200       RETURN
 201       END
 202 C*MODULE BLAS1   *DECK DNRM2
 203       DOUBLE PRECISION FUNCTION DNRM2 ( N, DX, INCX)
 204       INTEGER          NEXT
 205       DOUBLE PRECISION   DX(1), CUTLO, CUTHI, HITEST, SUM, XMAX,ZERO,ONE
 206       DATA   ZERO, ONE /0.0D+00, 1.0D+00/
 207 C
 208 C     EUCLIDEAN NORM OF THE N-VECTOR STORED IN DX() WITH STORAGE
 209 C     INCREMENT INCX .
 210 C     IF    N .LE. 0 RETURN WITH RESULT = 0.
 211 C     IF N .GE. 1 THEN INCX MUST BE .GE. 1
 212 C
 213 C           C.L.LAWSON, 1978 JAN 08
 214 C
 215 C     FOUR PHASE METHOD     USING TWO BUILT-IN CONSTANTS THAT ARE
 216 C     HOPEFULLY APPLICABLE TO ALL MACHINES.
 217 C         CUTLO = MAXIMUM OF  SQRT(U/EPS)  OVER ALL KNOWN MACHINES.
 218 C         CUTHI = MINIMUM OF  SQRT(V)      OVER ALL KNOWN MACHINES.
 219 C     WHERE
 220 C         EPS = SMALLEST NO. SUCH THAT EPS + 1. .GT. 1.
 221 C         U   = SMALLEST POSITIVE NO.   (UNDERFLOW LIMIT)
 222 C         V   = LARGEST  NO.            (OVERFLOW  LIMIT)
 223 C
 224 C     BRIEF OUTLINE OF ALGORITHM..
 225 C
 226 C     PHASE 1    SCANS ZERO COMPONENTS.
 227 C     MOVE TO PHASE 2 WHEN A COMPONENT IS NONZERO AND .LE. CUTLO
 228 C     MOVE TO PHASE 3 WHEN A COMPONENT IS .GT. CUTLO
 229 C     MOVE TO PHASE 4 WHEN A COMPONENT IS .GE. CUTHI/M
 230 C     WHERE M = N FOR X() REAL AND M = 2*N FOR COMPLEX.
 231 C
 232 C     VALUES FOR CUTLO AND CUTHI..
 233 C     FROM THE ENVIRONMENTAL PARAMETERS LISTED IN THE IMSL CONVERTER
 234 C     DOCUMENT THE LIMITING VALUES ARE AS FOLLOWS..
 235 C     CUTLO, S.P.   U/EPS = 2**(-102) FOR  HONEYWELL.  CLOSE SECONDS ARE
 236 C                   UNIVAC AND DEC AT 2**(-103)
 237 C                   THUS CUTLO = 2**(-51) = 4.44089E-16
 238 C     CUTHI, S.P.   V = 2**127 FOR UNIVAC, HONEYWELL, AND DEC.
 239 C                   THUS CUTHI = 2**(63.5) = 1.30438E19
 240 C     CUTLO, D.P.   U/EPS = 2**(-67) FOR HONEYWELL AND DEC.
 241 C                   THUS CUTLO = 2**(-33.5) = 8.23181D-11
 242 C     CUTHI, D.P.   SAME AS S.P.  CUTHI = 1.30438D+19
 243 C     DATA CUTLO, CUTHI / 8.232D-11,  1.304D+19 /
 244 C     DATA CUTLO, CUTHI / 4.441E-16,  1.304E19 /
 245       DATA CUTLO, CUTHI / 8.232D-11,  1.304D+19 /
 246 C
 247       J=0
 248       IF(N .GT. 0) GO TO 10
 249          DNRM2  = ZERO
 250          GO TO 300
 251 C
 252    10 ASSIGN 30 TO NEXT
 253       SUM = ZERO
 254       NN = N * INCX
 255 C                                                 BEGIN MAIN LOOP
 256       I = 1
 257    20    GO TO NEXT,(30, 50, 70, 110)
 258    30 IF( ABS(DX(I)) .GT. CUTLO) GO TO 85
 259       ASSIGN 50 TO NEXT
 260       XMAX = ZERO
 261 C
 262 C                        PHASE 1.  SUM IS ZERO
 263 C
 264    50 IF( DX(I) .EQ. ZERO) GO TO 200
 265       IF( ABS(DX(I)) .GT. CUTLO) GO TO 85
 266 C
 267 C                                PREPARE FOR PHASE 2.
 268       ASSIGN 70 TO NEXT
 269       GO TO 105
 270 C
 271 C                                PREPARE FOR PHASE 4.
 272 C
 273   100 I = J
 274       ASSIGN 110 TO NEXT
 275       SUM = (SUM / DX(I)) / DX(I)
 276   105 XMAX = ABS(DX(I))
 277       GO TO 115
 278 C
 279 C                   PHASE 2.  SUM IS SMALL.
 280 C                             SCALE TO AVOID DESTRUCTIVE UNDERFLOW.
 281 C
 282    70 IF( ABS(DX(I)) .GT. CUTLO ) GO TO 75
 283 C
 284 C                     COMMON CODE FOR PHASES 2 AND 4.
 285 C                     IN PHASE 4 SUM IS LARGE.  SCALE TO AVOID OVERFLOW.
 286 C
 287   110 IF( ABS(DX(I)) .LE. XMAX ) GO TO 115
 288          SUM = ONE + SUM * (XMAX / DX(I))**2
 289          XMAX = ABS(DX(I))
 290          GO TO 200
 291 C
 292   115 SUM = SUM + (DX(I)/XMAX)**2
 293       GO TO 200
 294 C
 295 C
 296 C                  PREPARE FOR PHASE 3.
 297 C
 298    75 SUM = (SUM * XMAX) * XMAX
 299 C
 300 C
 301 C     FOR REAL OR D.P. SET HITEST = CUTHI/N
 302 C     FOR COMPLEX      SET HITEST = CUTHI/(2*N)
 303 C
 304    85 HITEST = CUTHI/N
 305 C
 306 C                   PHASE 3.  SUM IS MID-RANGE.  NO SCALING.
 307 C
 308       DO 95 J =I,NN,INCX
 309       IF(ABS(DX(J)) .GE. HITEST) GO TO 100
 310    95    SUM = SUM + DX(J)**2
 311       DNRM2 = SQRT( SUM )
 312       GO TO 300
 313 C
 314   200 CONTINUE
 315       I = I + INCX
 316       IF ( I .LE. NN ) GO TO 20
 317 C
 318 C              END OF MAIN LOOP.
 319 C
 320 C              COMPUTE SQUARE ROOT AND ADJUST FOR SCALING.
 321 C
 322       DNRM2 = XMAX * SQRT(SUM)
 323   300 CONTINUE
 324       RETURN
 325       END
 326 C*MODULE BLAS1   *DECK DROT
 327       SUBROUTINE  DROT (N,DX,INCX,DY,INCY,C,S)
 328       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
 329       DIMENSION DX(1),DY(1)
 330 C
 331 C     APPLIES A PLANE ROTATION.
 332 C           DX(I) =  C*DX(I) + S*DY(I)
 333 C           DY(I) = -S*DX(I) + C*DY(I)
 334 C     JACK DONGARRA, LINPACK, 3/11/78.
 335 C
 336       IF(N.LE.0)RETURN
 337       IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20
 338 C
 339 C       CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS NOT EQUAL
 340 C         TO 1
 341 C
 342       IX = 1
 343       IY = 1
 344       IF(INCX.LT.0)IX = (-N+1)*INCX + 1
 345       IF(INCY.LT.0)IY = (-N+1)*INCY + 1
 346       DO 10 I = 1,N
 347         DTEMP = C*DX(IX) + S*DY(IY)
 348         DY(IY) = C*DY(IY) - S*DX(IX)
 349         DX(IX) = DTEMP
 350         IX = IX + INCX
 351         IY = IY + INCY
 352    10 CONTINUE
 353       RETURN
 354 C
 355 C       CODE FOR BOTH INCREMENTS EQUAL TO 1
 356 C
 357    20 DO 30 I = 1,N
 358         DTEMP = C*DX(I) + S*DY(I)
 359         DY(I) = C*DY(I) - S*DX(I)
 360         DX(I) = DTEMP
 361    30 CONTINUE
 362       RETURN
 363       END
 364 C*MODULE BLAS1   *DECK DROTG
 365       SUBROUTINE DROTG(DA,DB,C,S)
 366 C
 367 C     CONSTRUCT GIVENS PLANE ROTATION.
 368 C     JACK DONGARRA, LINPACK, 3/11/78.
 369 C
 370       DOUBLE PRECISION DA,DB,C,S,ROE,SCALE,R,Z
 371       DOUBLE PRECISION ZERO, ONE
 372 C
 373       PARAMETER (ZERO=0.0D+00, ONE=1.0D+00)
 374 C
 375 C-----------------------------------------------------------------------
 376 C
 377 C
 378       ROE = DB
 379       IF( ABS(DA) .GT. ABS(DB) ) ROE = DA
 380       SCALE = ABS(DA) + ABS(DB)
 381       IF( SCALE .NE. ZERO ) GO TO 10
 382          C = ONE
 383          S = ZERO
 384          R = ZERO
 385          GO TO 20
 386 C
 387    10 R = SCALE*SQRT((DA/SCALE)**2 + (DB/SCALE)**2)
 388       R = SIGN(ONE,ROE)*R
 389       C = DA/R
 390       S = DB/R
 391    20 Z = ONE
 392       IF( ABS(DA) .GT. ABS(DB) ) Z = S
 393       IF( ABS(DB) .GE. ABS(DA) .AND. C .NE. ZERO ) Z = ONE/C
 394       DA = R
 395       DB = Z
 396       RETURN
 397       END
 398 C*MODULE BLAS1   *DECK DSCAL
 399       SUBROUTINE  DSCAL(N,DA,DX,INCX)
 400       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
 401       DIMENSION DX(1)
 402 C
 403 C     SCALES A VECTOR BY A CONSTANT.
 404 C           DX(I) = DA * DX(I)
 405 C     USES UNROLLED LOOPS FOR INCREMENT EQUAL TO ONE.
 406 C     JACK DONGARRA, LINPACK, 3/11/78.
 407 C
 408       IF(N.LE.0)RETURN
 409       IF(INCX.EQ.1)GO TO 20
 410 C
 411 C        CODE FOR INCREMENT NOT EQUAL TO 1
 412 C
 413       NINCX = N*INCX
 414       DO 10 I = 1,NINCX,INCX
 415         DX(I) = DA*DX(I)
 416    10 CONTINUE
 417       RETURN
 418 C
 419 C        CODE FOR INCREMENT EQUAL TO 1
 420 C
 421 C
 422 C        CLEAN-UP LOOP
 423 C
 424    20 M = MOD(N,5)
 425       IF( M .EQ. 0 ) GO TO 40
 426       DO 30 I = 1,M
 427         DX(I) = DA*DX(I)
 428    30 CONTINUE
 429       IF( N .LT. 5 ) RETURN
 430    40 MP1 = M + 1
 431       DO 50 I = MP1,N,5
 432         DX(I) = DA*DX(I)
 433         DX(I + 1) = DA*DX(I + 1)
 434         DX(I + 2) = DA*DX(I + 2)
 435         DX(I + 3) = DA*DX(I + 3)
 436         DX(I + 4) = DA*DX(I + 4)
 437    50 CONTINUE
 438       RETURN
 439       END
 440 C*MODULE BLAS1   *DECK DSWAP
 441       SUBROUTINE  DSWAP (N,DX,INCX,DY,INCY)
 442       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
 443       DIMENSION DX(1),DY(1)
 444 C
 445 C     INTERCHANGES TWO VECTORS.
 446 C           DX(I) <==> DY(I)
 447 C     USES UNROLLED LOOPS FOR INCREMENTS EQUAL ONE.
 448 C     JACK DONGARRA, LINPACK, 3/11/78.
 449 C
 450       IF(N.LE.0)RETURN
 451       IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20
 452 C
 453 C       CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS NOT EQUAL
 454 C         TO 1
 455 C
 456       IX = 1
 457       IY = 1
 458       IF(INCX.LT.0)IX = (-N+1)*INCX + 1
 459       IF(INCY.LT.0)IY = (-N+1)*INCY + 1
 460       DO 10 I = 1,N
 461         DTEMP = DX(IX)
 462         DX(IX) = DY(IY)
 463         DY(IY) = DTEMP
 464         IX = IX + INCX
 465         IY = IY + INCY
 466    10 CONTINUE
 467       RETURN
 468 C
 469 C       CODE FOR BOTH INCREMENTS EQUAL TO 1
 470 C
 471 C
 472 C       CLEAN-UP LOOP
 473 C
 474    20 M = MOD(N,3)
 475       IF( M .EQ. 0 ) GO TO 40
 476       DO 30 I = 1,M
 477         DTEMP = DX(I)
 478         DX(I) = DY(I)
 479         DY(I) = DTEMP
 480    30 CONTINUE
 481       IF( N .LT. 3 ) RETURN
 482    40 MP1 = M + 1
 483       DO 50 I = MP1,N,3
 484         DTEMP = DX(I)
 485         DX(I) = DY(I)
 486         DY(I) = DTEMP
 487         DTEMP = DX(I + 1)
 488         DX(I + 1) = DY(I + 1)
 489         DY(I + 1) = DTEMP
 490         DTEMP = DX(I + 2)
 491         DX(I + 2) = DY(I + 2)
 492         DY(I + 2) = DTEMP
 493    50 CONTINUE
 494       RETURN
 495       END
 496 C*MODULE BLAS1   *DECK IDAMAX
 497       INTEGER FUNCTION IDAMAX(N,DX,INCX)
 498       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
 499       DIMENSION DX(1)
 500 C
 501 C     FINDS THE INDEX OF ELEMENT HAVING MAX. ABSOLUTE VALUE.
 502 C     JACK DONGARRA, LINPACK, 3/11/78.
 503 C
 504       IDAMAX = 0
 505       IF( N .LT. 1 ) RETURN
 506       IDAMAX = 1
 507       IF(N.EQ.1)RETURN
 508       IF(INCX.EQ.1)GO TO 20
 509 C
 510 C        CODE FOR INCREMENT NOT EQUAL TO 1
 511 C
 512       IX = 1
 513       RMAX = ABS(DX(1))
 514       IX = IX + INCX
 515       DO 10 I = 2,N
 516          IF(ABS(DX(IX)).LE.RMAX) GO TO 5
 517          IDAMAX = I
 518          RMAX = ABS(DX(IX))
 519     5    IX = IX + INCX
 520    10 CONTINUE
 521       RETURN
 522 C
 523 C        CODE FOR INCREMENT EQUAL TO 1
 524 C
 525    20 RMAX = ABS(DX(1))
 526       DO 30 I = 2,N
 527          IF(ABS(DX(I)).LE.RMAX) GO TO 30
 528          IDAMAX = I
 529          RMAX = ABS(DX(I))
 530    30 CONTINUE
 531       RETURN
 532       END
 533 C*MODULE BLAS    *DECK DGEMV
 534       SUBROUTINE DGEMV(FORMA,M,N,ALPHA,A,LDA,X,INCX,BETA,Y,INCY)
 535       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
 536       CHARACTER*1 FORMA
 537       DIMENSION A(LDA,*),X(*),Y(*)
 538       PARAMETER (ZERO=0.0D+00, ONE=1.0D+00)
 539 C
 540 C        CLONE OF -DGEMV- WRITTEN BY MIKE SCHMIDT
 541 C
 542       LOCY = 1
 543       IF(FORMA.EQ.'T') GO TO 200
 544 C
 545 C                  Y = ALPHA * A * X + BETA * Y
 546 C
 547       IF(ALPHA.EQ.ONE  .AND.  BETA.EQ.ZERO) THEN
 548          DO 110 I=1,M
 549             Y(LOCY) =       DDOT(N,A(I,1),LDA,X,INCX)
 550             LOCY = LOCY+INCY
 551   110    CONTINUE
 552       ELSE
 553          DO 120 I=1,M
 554             Y(LOCY) = ALPHA*DDOT(N,A(I,1),LDA,X,INCX) + BETA*Y(LOCY)
 555             LOCY = LOCY+INCY
 556   120    CONTINUE
 557       END IF
 558       RETURN
 559 C
 560 C                  Y = ALPHA * A-TRANSPOSE * X + BETA * Y
 561 C
 562   200 CONTINUE
 563       IF(ALPHA.EQ.ONE  .AND.  BETA.EQ.ZERO) THEN
 564          DO 210 I=1,N
 565             Y(LOCY) =       DDOT(M,A(1,I),1,X,INCX)
 566             LOCY = LOCY+INCY
 567   210    CONTINUE
 568       ELSE
 569          DO 220 I=1,N
 570             Y(LOCY) = ALPHA*DDOT(M,A(1,I),1,X,INCX) + BETA*Y(LOCY)
 571             LOCY = LOCY+INCY
 572   220    CONTINUE
 573       END IF
 574       RETURN
 575       END