wham and cluster_wham Adam's new constr_dist multichain
[unres.git] / 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