X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD_DFA%2Fblas.f;fp=source%2Funres%2Fsrc_MD_DFA%2Fblas.f;h=0000000000000000000000000000000000000000;hb=664d495e70d14eed4e97f7b8efd2e107dee2fd4e;hp=142d82153839bdd94918381126daf55b7ca41ef8;hpb=77276d5043cefaf512451c3ad9f669ed22b90d04;p=unres.git diff --git a/source/unres/src_MD_DFA/blas.f b/source/unres/src_MD_DFA/blas.f deleted file mode 100644 index 142d821..0000000 --- a/source/unres/src_MD_DFA/blas.f +++ /dev/null @@ -1,575 +0,0 @@ -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