X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD-NEWSC%2Fblas.f;fp=source%2Funres%2Fsrc_MD-NEWSC%2Fblas.f;h=142d82153839bdd94918381126daf55b7ca41ef8;hb=d101c97dea752458d76055fdbae49c26fff03c1f;hp=0000000000000000000000000000000000000000;hpb=325eda160c9ad2982501e091ca40606a29043712;p=unres.git diff --git a/source/unres/src_MD-NEWSC/blas.f b/source/unres/src_MD-NEWSC/blas.f new file mode 100644 index 0000000..142d821 --- /dev/null +++ b/source/unres/src_MD-NEWSC/blas.f @@ -0,0 +1,575 @@ +C 10 NOV 94 - MWS - DNRM2: REMOVE FTNCHECK WARNINGS +C 11 JUN 94 - MWS - INCLUDE A COPY OF DGEMV (LEVEL TWO ROUTINE) +C 11 AUG 87 - MWS - SANITIZE FLOATING POINT CONSTANTS IN DNRM2 +C 26 MAR 87 - MWS - USE GENERIC SIGN IN DROTG +C 28 NOV 86 - STE - SUPPLY ALL LEVEL ONE BLAS +C 7 JUL 86 - JAB - SANITIZE FLOATING POINT CONSTANTS +C +C BASIC LINEAR ALGEBRA SUBPROGRAMS (BLAS) FROM LINPACK (LEVEL 1) +C +C THIS MODULE SHOULD BE COMPILED ONLY IF SPECIALLY CODED +C VERSIONS OF THESE ROUTINES ARE NOT AVAILABLE ON THE TARGET MACHINE +C +C*MODULE BLAS1 *DECK DASUM + DOUBLE PRECISION FUNCTION DASUM(N,DX,INCX) +C +C TAKES THE SUM OF THE ABSOLUTE VALUES. +C JACK DONGARRA, LINPACK, 3/11/78. +C + DOUBLE PRECISION DX(1),DTEMP + INTEGER I,INCX,M,MP1,N,NINCX +C + DASUM = 0.0D+00 + DTEMP = 0.0D+00 + IF(N.LE.0)RETURN + IF(INCX.EQ.1)GO TO 20 +C +C CODE FOR INCREMENT NOT EQUAL TO 1 +C + NINCX = N*INCX + DO 10 I = 1,NINCX,INCX + DTEMP = DTEMP + ABS(DX(I)) + 10 CONTINUE + DASUM = DTEMP + RETURN +C +C CODE FOR INCREMENT EQUAL TO 1 +C +C +C CLEAN-UP LOOP +C + 20 M = MOD(N,6) + IF( M .EQ. 0 ) GO TO 40 + DO 30 I = 1,M + DTEMP = DTEMP + ABS(DX(I)) + 30 CONTINUE + IF( N .LT. 6 ) GO TO 60 + 40 MP1 = M + 1 + DO 50 I = MP1,N,6 + DTEMP = DTEMP + ABS(DX(I)) + ABS(DX(I + 1)) + ABS(DX(I + 2)) + * + ABS(DX(I + 3)) + ABS(DX(I + 4)) + ABS(DX(I + 5)) + 50 CONTINUE + 60 DASUM = DTEMP + RETURN + END +C*MODULE BLAS1 *DECK DAXPY + SUBROUTINE DAXPY(N,DA,DX,INCX,DY,INCY) + IMPLICIT DOUBLE PRECISION(A-H,O-Z) + DIMENSION DX(1),DY(1) +C +C CONSTANT TIMES A VECTOR PLUS A VECTOR. +C DY(I) = DY(I) + DA * DX(I) +C USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE. +C JACK DONGARRA, LINPACK, 3/11/78. +C + IF(N.LE.0)RETURN + IF (DA .EQ. 0.0D+00) RETURN + IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20 +C +C CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS +C NOT EQUAL TO 1 +C + IX = 1 + IY = 1 + IF(INCX.LT.0)IX = (-N+1)*INCX + 1 + IF(INCY.LT.0)IY = (-N+1)*INCY + 1 + DO 10 I = 1,N + DY(IY) = DY(IY) + DA*DX(IX) + IX = IX + INCX + IY = IY + INCY + 10 CONTINUE + RETURN +C +C CODE FOR BOTH INCREMENTS EQUAL TO 1 +C +C +C CLEAN-UP LOOP +C + 20 M = MOD(N,4) + IF( M .EQ. 0 ) GO TO 40 + DO 30 I = 1,M + DY(I) = DY(I) + DA*DX(I) + 30 CONTINUE + IF( N .LT. 4 ) RETURN + 40 MP1 = M + 1 + DO 50 I = MP1,N,4 + DY(I) = DY(I) + DA*DX(I) + DY(I + 1) = DY(I + 1) + DA*DX(I + 1) + DY(I + 2) = DY(I + 2) + DA*DX(I + 2) + DY(I + 3) = DY(I + 3) + DA*DX(I + 3) + 50 CONTINUE + RETURN + END +C*MODULE BLAS1 *DECK DCOPY + SUBROUTINE DCOPY(N,DX,INCX,DY,INCY) + IMPLICIT DOUBLE PRECISION(A-H,O-Z) + DIMENSION DX(*),DY(*) +C +C COPIES A VECTOR. +C DY(I) <== DX(I) +C USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE. +C JACK DONGARRA, LINPACK, 3/11/78. +C + IF(N.LE.0)RETURN + IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20 +C +C CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS +C NOT EQUAL TO 1 +C + IX = 1 + IY = 1 + IF(INCX.LT.0)IX = (-N+1)*INCX + 1 + IF(INCY.LT.0)IY = (-N+1)*INCY + 1 + DO 10 I = 1,N + DY(IY) = DX(IX) + IX = IX + INCX + IY = IY + INCY + 10 CONTINUE + RETURN +C +C CODE FOR BOTH INCREMENTS EQUAL TO 1 +C +C +C CLEAN-UP LOOP +C + 20 M = MOD(N,7) + IF( M .EQ. 0 ) GO TO 40 + DO 30 I = 1,M + DY(I) = DX(I) + 30 CONTINUE + IF( N .LT. 7 ) RETURN + 40 MP1 = M + 1 + DO 50 I = MP1,N,7 + DY(I) = DX(I) + DY(I + 1) = DX(I + 1) + DY(I + 2) = DX(I + 2) + DY(I + 3) = DX(I + 3) + DY(I + 4) = DX(I + 4) + DY(I + 5) = DX(I + 5) + DY(I + 6) = DX(I + 6) + 50 CONTINUE + RETURN + END +C*MODULE BLAS1 *DECK DDOT + DOUBLE PRECISION FUNCTION DDOT(N,DX,INCX,DY,INCY) + IMPLICIT DOUBLE PRECISION(A-H,O-Z) + DIMENSION DX(1),DY(1) +C +C FORMS THE DOT PRODUCT OF TWO VECTORS. +C DOT = DX(I) * DY(I) +C USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE. +C JACK DONGARRA, LINPACK, 3/11/78. +C + DDOT = 0.0D+00 + DTEMP = 0.0D+00 + IF(N.LE.0)RETURN + IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20 +C +C CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS +C NOT EQUAL TO 1 +C + IX = 1 + IY = 1 + IF(INCX.LT.0)IX = (-N+1)*INCX + 1 + IF(INCY.LT.0)IY = (-N+1)*INCY + 1 + DO 10 I = 1,N + DTEMP = DTEMP + DX(IX)*DY(IY) + IX = IX + INCX + IY = IY + INCY + 10 CONTINUE + DDOT = DTEMP + RETURN +C +C CODE FOR BOTH INCREMENTS EQUAL TO 1 +C +C +C CLEAN-UP LOOP +C + 20 M = MOD(N,5) + IF( M .EQ. 0 ) GO TO 40 + DO 30 I = 1,M + DTEMP = DTEMP + DX(I)*DY(I) + 30 CONTINUE + IF( N .LT. 5 ) GO TO 60 + 40 MP1 = M + 1 + DO 50 I = MP1,N,5 + DTEMP = DTEMP + DX(I)*DY(I) + DX(I + 1)*DY(I + 1) + + * DX(I + 2)*DY(I + 2) + DX(I + 3)*DY(I + 3) + DX(I + 4)*DY(I + 4) + 50 CONTINUE + 60 DDOT = DTEMP + RETURN + END +C*MODULE BLAS1 *DECK DNRM2 + DOUBLE PRECISION FUNCTION DNRM2 ( N, DX, INCX) + INTEGER NEXT + DOUBLE PRECISION DX(1), CUTLO, CUTHI, HITEST, SUM, XMAX,ZERO,ONE + DATA ZERO, ONE /0.0D+00, 1.0D+00/ +C +C EUCLIDEAN NORM OF THE N-VECTOR STORED IN DX() WITH STORAGE +C INCREMENT INCX . +C IF N .LE. 0 RETURN WITH RESULT = 0. +C IF N .GE. 1 THEN INCX MUST BE .GE. 1 +C +C C.L.LAWSON, 1978 JAN 08 +C +C FOUR PHASE METHOD USING TWO BUILT-IN CONSTANTS THAT ARE +C HOPEFULLY APPLICABLE TO ALL MACHINES. +C CUTLO = MAXIMUM OF SQRT(U/EPS) OVER ALL KNOWN MACHINES. +C CUTHI = MINIMUM OF SQRT(V) OVER ALL KNOWN MACHINES. +C WHERE +C EPS = SMALLEST NO. SUCH THAT EPS + 1. .GT. 1. +C U = SMALLEST POSITIVE NO. (UNDERFLOW LIMIT) +C V = LARGEST NO. (OVERFLOW LIMIT) +C +C BRIEF OUTLINE OF ALGORITHM.. +C +C PHASE 1 SCANS ZERO COMPONENTS. +C MOVE TO PHASE 2 WHEN A COMPONENT IS NONZERO AND .LE. CUTLO +C MOVE TO PHASE 3 WHEN A COMPONENT IS .GT. CUTLO +C MOVE TO PHASE 4 WHEN A COMPONENT IS .GE. CUTHI/M +C WHERE M = N FOR X() REAL AND M = 2*N FOR COMPLEX. +C +C VALUES FOR CUTLO AND CUTHI.. +C FROM THE ENVIRONMENTAL PARAMETERS LISTED IN THE IMSL CONVERTER +C DOCUMENT THE LIMITING VALUES ARE AS FOLLOWS.. +C CUTLO, S.P. U/EPS = 2**(-102) FOR HONEYWELL. CLOSE SECONDS ARE +C UNIVAC AND DEC AT 2**(-103) +C THUS CUTLO = 2**(-51) = 4.44089E-16 +C CUTHI, S.P. V = 2**127 FOR UNIVAC, HONEYWELL, AND DEC. +C THUS CUTHI = 2**(63.5) = 1.30438E19 +C CUTLO, D.P. U/EPS = 2**(-67) FOR HONEYWELL AND DEC. +C THUS CUTLO = 2**(-33.5) = 8.23181D-11 +C CUTHI, D.P. SAME AS S.P. CUTHI = 1.30438D+19 +C DATA CUTLO, CUTHI / 8.232D-11, 1.304D+19 / +C DATA CUTLO, CUTHI / 4.441E-16, 1.304E19 / + DATA CUTLO, CUTHI / 8.232D-11, 1.304D+19 / +C + J=0 + IF(N .GT. 0) GO TO 10 + DNRM2 = ZERO + GO TO 300 +C + 10 ASSIGN 30 TO NEXT + SUM = ZERO + NN = N * INCX +C BEGIN MAIN LOOP + I = 1 + 20 GO TO NEXT,(30, 50, 70, 110) + 30 IF( ABS(DX(I)) .GT. CUTLO) GO TO 85 + ASSIGN 50 TO NEXT + XMAX = ZERO +C +C PHASE 1. SUM IS ZERO +C + 50 IF( DX(I) .EQ. ZERO) GO TO 200 + IF( ABS(DX(I)) .GT. CUTLO) GO TO 85 +C +C PREPARE FOR PHASE 2. + ASSIGN 70 TO NEXT + GO TO 105 +C +C PREPARE FOR PHASE 4. +C + 100 I = J + ASSIGN 110 TO NEXT + SUM = (SUM / DX(I)) / DX(I) + 105 XMAX = ABS(DX(I)) + GO TO 115 +C +C PHASE 2. SUM IS SMALL. +C SCALE TO AVOID DESTRUCTIVE UNDERFLOW. +C + 70 IF( ABS(DX(I)) .GT. CUTLO ) GO TO 75 +C +C COMMON CODE FOR PHASES 2 AND 4. +C IN PHASE 4 SUM IS LARGE. SCALE TO AVOID OVERFLOW. +C + 110 IF( ABS(DX(I)) .LE. XMAX ) GO TO 115 + SUM = ONE + SUM * (XMAX / DX(I))**2 + XMAX = ABS(DX(I)) + GO TO 200 +C + 115 SUM = SUM + (DX(I)/XMAX)**2 + GO TO 200 +C +C +C PREPARE FOR PHASE 3. +C + 75 SUM = (SUM * XMAX) * XMAX +C +C +C FOR REAL OR D.P. SET HITEST = CUTHI/N +C FOR COMPLEX SET HITEST = CUTHI/(2*N) +C + 85 HITEST = CUTHI/N +C +C PHASE 3. SUM IS MID-RANGE. NO SCALING. +C + DO 95 J =I,NN,INCX + IF(ABS(DX(J)) .GE. HITEST) GO TO 100 + 95 SUM = SUM + DX(J)**2 + DNRM2 = SQRT( SUM ) + GO TO 300 +C + 200 CONTINUE + I = I + INCX + IF ( I .LE. NN ) GO TO 20 +C +C END OF MAIN LOOP. +C +C COMPUTE SQUARE ROOT AND ADJUST FOR SCALING. +C + DNRM2 = XMAX * SQRT(SUM) + 300 CONTINUE + RETURN + END +C*MODULE BLAS1 *DECK DROT + SUBROUTINE DROT (N,DX,INCX,DY,INCY,C,S) + IMPLICIT DOUBLE PRECISION(A-H,O-Z) + DIMENSION DX(1),DY(1) +C +C APPLIES A PLANE ROTATION. +C DX(I) = C*DX(I) + S*DY(I) +C DY(I) = -S*DX(I) + C*DY(I) +C JACK DONGARRA, LINPACK, 3/11/78. +C + IF(N.LE.0)RETURN + IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20 +C +C CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS NOT EQUAL +C TO 1 +C + IX = 1 + IY = 1 + IF(INCX.LT.0)IX = (-N+1)*INCX + 1 + IF(INCY.LT.0)IY = (-N+1)*INCY + 1 + DO 10 I = 1,N + DTEMP = C*DX(IX) + S*DY(IY) + DY(IY) = C*DY(IY) - S*DX(IX) + DX(IX) = DTEMP + IX = IX + INCX + IY = IY + INCY + 10 CONTINUE + RETURN +C +C CODE FOR BOTH INCREMENTS EQUAL TO 1 +C + 20 DO 30 I = 1,N + DTEMP = C*DX(I) + S*DY(I) + DY(I) = C*DY(I) - S*DX(I) + DX(I) = DTEMP + 30 CONTINUE + RETURN + END +C*MODULE BLAS1 *DECK DROTG + SUBROUTINE DROTG(DA,DB,C,S) +C +C CONSTRUCT GIVENS PLANE ROTATION. +C JACK DONGARRA, LINPACK, 3/11/78. +C + DOUBLE PRECISION DA,DB,C,S,ROE,SCALE,R,Z + DOUBLE PRECISION ZERO, ONE +C + PARAMETER (ZERO=0.0D+00, ONE=1.0D+00) +C +C----------------------------------------------------------------------- +C +C + ROE = DB + IF( ABS(DA) .GT. ABS(DB) ) ROE = DA + SCALE = ABS(DA) + ABS(DB) + IF( SCALE .NE. ZERO ) GO TO 10 + C = ONE + S = ZERO + R = ZERO + GO TO 20 +C + 10 R = SCALE*SQRT((DA/SCALE)**2 + (DB/SCALE)**2) + R = SIGN(ONE,ROE)*R + C = DA/R + S = DB/R + 20 Z = ONE + IF( ABS(DA) .GT. ABS(DB) ) Z = S + IF( ABS(DB) .GE. ABS(DA) .AND. C .NE. ZERO ) Z = ONE/C + DA = R + DB = Z + RETURN + END +C*MODULE BLAS1 *DECK DSCAL + SUBROUTINE DSCAL(N,DA,DX,INCX) + IMPLICIT DOUBLE PRECISION(A-H,O-Z) + DIMENSION DX(1) +C +C SCALES A VECTOR BY A CONSTANT. +C DX(I) = DA * DX(I) +C USES UNROLLED LOOPS FOR INCREMENT EQUAL TO ONE. +C JACK DONGARRA, LINPACK, 3/11/78. +C + IF(N.LE.0)RETURN + IF(INCX.EQ.1)GO TO 20 +C +C CODE FOR INCREMENT NOT EQUAL TO 1 +C + NINCX = N*INCX + DO 10 I = 1,NINCX,INCX + DX(I) = DA*DX(I) + 10 CONTINUE + RETURN +C +C CODE FOR INCREMENT EQUAL TO 1 +C +C +C CLEAN-UP LOOP +C + 20 M = MOD(N,5) + IF( M .EQ. 0 ) GO TO 40 + DO 30 I = 1,M + DX(I) = DA*DX(I) + 30 CONTINUE + IF( N .LT. 5 ) RETURN + 40 MP1 = M + 1 + DO 50 I = MP1,N,5 + DX(I) = DA*DX(I) + DX(I + 1) = DA*DX(I + 1) + DX(I + 2) = DA*DX(I + 2) + DX(I + 3) = DA*DX(I + 3) + DX(I + 4) = DA*DX(I + 4) + 50 CONTINUE + RETURN + END +C*MODULE BLAS1 *DECK DSWAP + SUBROUTINE DSWAP (N,DX,INCX,DY,INCY) + IMPLICIT DOUBLE PRECISION(A-H,O-Z) + DIMENSION DX(1),DY(1) +C +C INTERCHANGES TWO VECTORS. +C DX(I) <==> DY(I) +C USES UNROLLED LOOPS FOR INCREMENTS EQUAL ONE. +C JACK DONGARRA, LINPACK, 3/11/78. +C + IF(N.LE.0)RETURN + IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20 +C +C CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS NOT EQUAL +C TO 1 +C + IX = 1 + IY = 1 + IF(INCX.LT.0)IX = (-N+1)*INCX + 1 + IF(INCY.LT.0)IY = (-N+1)*INCY + 1 + DO 10 I = 1,N + DTEMP = DX(IX) + DX(IX) = DY(IY) + DY(IY) = DTEMP + IX = IX + INCX + IY = IY + INCY + 10 CONTINUE + RETURN +C +C CODE FOR BOTH INCREMENTS EQUAL TO 1 +C +C +C CLEAN-UP LOOP +C + 20 M = MOD(N,3) + IF( M .EQ. 0 ) GO TO 40 + DO 30 I = 1,M + DTEMP = DX(I) + DX(I) = DY(I) + DY(I) = DTEMP + 30 CONTINUE + IF( N .LT. 3 ) RETURN + 40 MP1 = M + 1 + DO 50 I = MP1,N,3 + DTEMP = DX(I) + DX(I) = DY(I) + DY(I) = DTEMP + DTEMP = DX(I + 1) + DX(I + 1) = DY(I + 1) + DY(I + 1) = DTEMP + DTEMP = DX(I + 2) + DX(I + 2) = DY(I + 2) + DY(I + 2) = DTEMP + 50 CONTINUE + RETURN + END +C*MODULE BLAS1 *DECK IDAMAX + INTEGER FUNCTION IDAMAX(N,DX,INCX) + IMPLICIT DOUBLE PRECISION(A-H,O-Z) + DIMENSION DX(1) +C +C FINDS THE INDEX OF ELEMENT HAVING MAX. ABSOLUTE VALUE. +C JACK DONGARRA, LINPACK, 3/11/78. +C + IDAMAX = 0 + IF( N .LT. 1 ) RETURN + IDAMAX = 1 + IF(N.EQ.1)RETURN + IF(INCX.EQ.1)GO TO 20 +C +C CODE FOR INCREMENT NOT EQUAL TO 1 +C + IX = 1 + RMAX = ABS(DX(1)) + IX = IX + INCX + DO 10 I = 2,N + IF(ABS(DX(IX)).LE.RMAX) GO TO 5 + IDAMAX = I + RMAX = ABS(DX(IX)) + 5 IX = IX + INCX + 10 CONTINUE + RETURN +C +C CODE FOR INCREMENT EQUAL TO 1 +C + 20 RMAX = ABS(DX(1)) + DO 30 I = 2,N + IF(ABS(DX(I)).LE.RMAX) GO TO 30 + IDAMAX = I + RMAX = ABS(DX(I)) + 30 CONTINUE + RETURN + END +C*MODULE BLAS *DECK DGEMV + SUBROUTINE DGEMV(FORMA,M,N,ALPHA,A,LDA,X,INCX,BETA,Y,INCY) + IMPLICIT DOUBLE PRECISION(A-H,O-Z) + CHARACTER*1 FORMA + DIMENSION A(LDA,*),X(*),Y(*) + PARAMETER (ZERO=0.0D+00, ONE=1.0D+00) +C +C CLONE OF -DGEMV- WRITTEN BY MIKE SCHMIDT +C + LOCY = 1 + IF(FORMA.EQ.'T') GO TO 200 +C +C Y = ALPHA * A * X + BETA * Y +C + IF(ALPHA.EQ.ONE .AND. BETA.EQ.ZERO) THEN + DO 110 I=1,M + Y(LOCY) = DDOT(N,A(I,1),LDA,X,INCX) + LOCY = LOCY+INCY + 110 CONTINUE + ELSE + DO 120 I=1,M + Y(LOCY) = ALPHA*DDOT(N,A(I,1),LDA,X,INCX) + BETA*Y(LOCY) + LOCY = LOCY+INCY + 120 CONTINUE + END IF + RETURN +C +C Y = ALPHA * A-TRANSPOSE * X + BETA * Y +C + 200 CONTINUE + IF(ALPHA.EQ.ONE .AND. BETA.EQ.ZERO) THEN + DO 210 I=1,N + Y(LOCY) = DDOT(M,A(1,I),1,X,INCX) + LOCY = LOCY+INCY + 210 CONTINUE + ELSE + DO 220 I=1,N + Y(LOCY) = ALPHA*DDOT(M,A(1,I),1,X,INCX) + BETA*Y(LOCY) + LOCY = LOCY+INCY + 220 CONTINUE + END IF + RETURN + END