Mercurial > octave-nkf
changeset 2608:bac14003d9bb
[project @ 1997-01-18 00:11:48 by jwe]
line wrap: on
line diff
--- a/README.Linux Fri Jan 17 18:54:19 1997 +0000 +++ b/README.Linux Sat Jan 18 00:13:08 1997 +0000 @@ -24,11 +24,13 @@ * Linux kernel 2.0.6 * gcc/g++ 2.7.2 + * g77 0.5.18 * libg++/libstdc++ 2.7.1.0 * libm 5.0.5 * libc 5.2.18 * libncurses 3.0 * ld.so 1.7.14 + * binutils 2.6 I know from experience that the versions listed above seem to work well together. But if you have a newer version of the kernel, you may @@ -99,6 +101,11 @@ ar cq libieee.a +NOTE: you should fix this problem (either by editing the specs file or +by creating the library) *before* running configure and compiling +Octave. Otherwise, configure may incorrectly determine that your +system doesn't have support for some IEEE math functions. + My system doesn't have g77 -------------------------- @@ -140,4 +147,4 @@ University of Wisconsin-Madison Department of Chemical Engineering -Thu Dec 19 13:07:46 1996 +Wed Jan 15 20:04:54 1997
--- a/libcruft/ChangeLog Fri Jan 17 18:54:19 1997 +0000 +++ b/libcruft/ChangeLog Sat Jan 18 00:13:08 1997 +0000 @@ -1,3 +1,7 @@ +Wed Jan 15 21:04:29 1997 John W. Eaton <jwe@bevo.che.wisc.edu> + + * blas/*.f: Update to latest version from Netlib. + Tue Jan 7 00:17:17 1997 John W. Eaton <jwe@bevo.che.wisc.edu> * Version 2.0.1 released.
--- a/libcruft/blas/dasum.f Fri Jan 17 18:54:19 1997 +0000 +++ b/libcruft/blas/dasum.f Sat Jan 18 00:13:08 1997 +0000 @@ -1,41 +1,43 @@ - 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.0D0 - DTEMP = 0.0D0 - 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 + DABS(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 + DABS(DX(I)) - 30 CONTINUE - IF( N .LT. 6 ) GO TO 60 - 40 MP1 = M + 1 - DO 50 I = MP1,N,6 - DTEMP = DTEMP + DABS(DX(I)) + DABS(DX(I + 1)) + DABS(DX(I + 2)) - * + DABS(DX(I + 3)) + DABS(DX(I + 4)) + DABS(DX(I + 5)) - 50 CONTINUE - 60 DASUM = DTEMP - RETURN - END + double precision function dasum(n,dx,incx) +c +c takes the sum of the absolute values. +c jack dongarra, linpack, 3/11/78. +c modified 3/93 to return if incx .le. 0. +c modified 12/3/93, array(1) declarations changed to array(*) +c + double precision dx(*),dtemp + integer i,incx,m,mp1,n,nincx +c + dasum = 0.0d0 + dtemp = 0.0d0 + if( n.le.0 .or. incx.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 + dabs(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 + dabs(dx(i)) + 30 continue + if( n .lt. 6 ) go to 60 + 40 mp1 = m + 1 + do 50 i = mp1,n,6 + dtemp = dtemp + dabs(dx(i)) + dabs(dx(i + 1)) + dabs(dx(i + 2)) + * + dabs(dx(i + 3)) + dabs(dx(i + 4)) + dabs(dx(i + 5)) + 50 continue + 60 dasum = dtemp + return + end
--- a/libcruft/blas/daxpy.f Fri Jan 17 18:54:19 1997 +0000 +++ b/libcruft/blas/daxpy.f Sat Jan 18 00:13:08 1997 +0000 @@ -1,47 +1,48 @@ - SUBROUTINE DAXPY(N,DA,DX,INCX,DY,INCY) -C -C CONSTANT TIMES A VECTOR PLUS A VECTOR. -C USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE. -C JACK DONGARRA, LINPACK, 3/11/78. -C - DOUBLE PRECISION DX(1),DY(1),DA - INTEGER I,INCX,INCY,IXIY,M,MP1,N -C - IF(N.LE.0)RETURN - IF (DA .EQ. 0.0D0) 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 + subroutine daxpy(n,da,dx,incx,dy,incy) +c +c constant times a vector plus a vector. +c uses unrolled loops for increments equal to one. +c jack dongarra, linpack, 3/11/78. +c modified 12/3/93, array(1) declarations changed to array(*) +c + double precision dx(*),dy(*),da + integer i,incx,incy,ix,iy,m,mp1,n +c + if(n.le.0)return + if (da .eq. 0.0d0) 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
--- a/libcruft/blas/dcopy.f Fri Jan 17 18:54:19 1997 +0000 +++ b/libcruft/blas/dcopy.f Sat Jan 18 00:13:08 1997 +0000 @@ -1,49 +1,50 @@ - SUBROUTINE DCOPY(N,DX,INCX,DY,INCY) -C -C COPIES A VECTOR, X, TO A VECTOR, Y. -C USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE. -C JACK DONGARRA, LINPACK, 3/11/78. -C - DOUBLE PRECISION DX(1),DY(1) - INTEGER I,INCX,INCY,IX,IY,M,MP1,N -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 + subroutine dcopy(n,dx,incx,dy,incy) +c +c copies a vector, x, to a vector, y. +c uses unrolled loops for increments equal to one. +c jack dongarra, linpack, 3/11/78. +c modified 12/3/93, array(1) declarations changed to array(*) +c + double precision dx(*),dy(*) + integer i,incx,incy,ix,iy,m,mp1,n +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
--- a/libcruft/blas/ddot.f Fri Jan 17 18:54:19 1997 +0000 +++ b/libcruft/blas/ddot.f Sat Jan 18 00:13:08 1997 +0000 @@ -1,48 +1,49 @@ - DOUBLE PRECISION FUNCTION DDOT(N,DX,INCX,DY,INCY) -C -C FORMS THE DOT PRODUCT OF TWO VECTORS. -C USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE. -C JACK DONGARRA, LINPACK, 3/11/78. -C - DOUBLE PRECISION DX(1),DY(1),DTEMP - INTEGER I,INCX,INCY,IX,IY,M,MP1,N -C - DDOT = 0.0D0 - DTEMP = 0.0D0 - 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 + double precision function ddot(n,dx,incx,dy,incy) +c +c forms the dot product of two vectors. +c uses unrolled loops for increments equal to one. +c jack dongarra, linpack, 3/11/78. +c modified 12/3/93, array(1) declarations changed to array(*) +c + double precision dx(*),dy(*),dtemp + integer i,incx,incy,ix,iy,m,mp1,n +c + ddot = 0.0d0 + dtemp = 0.0d0 + 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
--- a/libcruft/blas/dgemm.f Fri Jan 17 18:54:19 1997 +0000 +++ b/libcruft/blas/dgemm.f Sat Jan 18 00:13:08 1997 +0000 @@ -1,37 +1,3 @@ -************************************************************************ -* -* File of the DOUBLE PRECISION Level-3 BLAS. -* ========================================== -* -* SUBROUTINE DGEMM ( TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, -* $ BETA, C, LDC ) -* -* SUBROUTINE DSYMM ( SIDE, UPLO, M, N, ALPHA, A, LDA, B, LDB, -* $ BETA, C, LDC ) -* -* SUBROUTINE DSYRK ( UPLO, TRANS, N, K, ALPHA, A, LDA, -* $ BETA, C, LDC ) -* -* SUBROUTINE DSYR2K( UPLO, TRANS, N, K, ALPHA, A, LDA, B, LDB, -* $ BETA, C, LDC ) -* -* SUBROUTINE DTRMM ( SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, -* $ B, LDB ) -* -* SUBROUTINE DTRSM ( SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, -* $ B, LDB ) -* -* See: -* -* Dongarra J. J., Du Croz J. J., Duff I. and Hammarling S. -* A set of Level 3 Basic Linear Algebra Subprograms. Technical -* Memorandum No.88 (Revision 1), Mathematics and Computer Science -* Division, Argonne National Laboratory, 9700 South Cass Avenue, -* Argonne, Illinois 60439. -* -* -************************************************************************ -* SUBROUTINE DGEMM ( TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, $ BETA, C, LDC ) * .. Scalar Arguments ..
--- a/libcruft/blas/dgemv.f Fri Jan 17 18:54:19 1997 +0000 +++ b/libcruft/blas/dgemv.f Sat Jan 18 00:13:08 1997 +0000 @@ -1,63 +1,3 @@ -* -************************************************************************ -* -* File of the DOUBLE PRECISION Level-2 BLAS. -* =========================================== -* -* SUBROUTINE DGEMV ( TRANS, M, N, ALPHA, A, LDA, X, INCX, -* $ BETA, Y, INCY ) -* -* SUBROUTINE DGBMV ( TRANS, M, N, KL, KU, ALPHA, A, LDA, X, INCX, -* $ BETA, Y, INCY ) -* -* SUBROUTINE DSYMV ( UPLO, N, ALPHA, A, LDA, X, INCX, -* $ BETA, Y, INCY ) -* -* SUBROUTINE DSBMV ( UPLO, N, K, ALPHA, A, LDA, X, INCX, -* $ BETA, Y, INCY ) -* -* SUBROUTINE DSPMV ( UPLO, N, ALPHA, AP, X, INCX, BETA, Y, INCY ) -* -* SUBROUTINE DTRMV ( UPLO, TRANS, DIAG, N, A, LDA, X, INCX ) -* -* SUBROUTINE DTBMV ( UPLO, TRANS, DIAG, N, K, A, LDA, X, INCX ) -* -* SUBROUTINE DTPMV ( UPLO, TRANS, DIAG, N, AP, X, INCX ) -* -* SUBROUTINE DTRSV ( UPLO, TRANS, DIAG, N, A, LDA, X, INCX ) -* -* SUBROUTINE DTBSV ( UPLO, TRANS, DIAG, N, K, A, LDA, X, INCX ) -* -* SUBROUTINE DTPSV ( UPLO, TRANS, DIAG, N, AP, X, INCX ) -* -* SUBROUTINE DGER ( M, N, ALPHA, X, INCX, Y, INCY, A, LDA ) -* -* SUBROUTINE DSYR ( UPLO, N, ALPHA, X, INCX, A, LDA ) -* -* SUBROUTINE DSPR ( UPLO, N, ALPHA, X, INCX, AP ) -* -* SUBROUTINE DSYR2 ( UPLO, N, ALPHA, X, INCX, Y, INCY, A, LDA ) -* -* SUBROUTINE DSPR2 ( UPLO, N, ALPHA, X, INCX, Y, INCY, AP ) -* -* See: -* -* Dongarra J. J., Du Croz J. J., Hammarling S. and Hanson R. J.. -* An extended set of Fortran Basic Linear Algebra Subprograms. -* -* Technical Memoranda Nos. 41 (revision 3) and 81, Mathematics -* and Computer Science Division, Argonne National Laboratory, -* 9700 South Cass Avenue, Argonne, Illinois 60439, US. -* -* Or -* -* NAG Technical Reports TR3/87 and TR4/87, Numerical Algorithms -* Group Ltd., NAG Central Office, 256 Banbury Road, Oxford -* OX2 7DE, UK, and Numerical Algorithms Group Inc., 1101 31st -* Street, Suite 100, Downers Grove, Illinois 60515-1263, USA. -* -************************************************************************ -* SUBROUTINE DGEMV ( TRANS, M, N, ALPHA, A, LDA, X, INCX, $ BETA, Y, INCY ) * .. Scalar Arguments ..
--- a/libcruft/blas/dger.f Fri Jan 17 18:54:19 1997 +0000 +++ b/libcruft/blas/dger.f Sat Jan 18 00:13:08 1997 +0000 @@ -1,6 +1,3 @@ -* -************************************************************************ -* SUBROUTINE DGER ( M, N, ALPHA, X, INCX, Y, INCY, A, LDA ) * .. Scalar Arguments .. DOUBLE PRECISION ALPHA
--- a/libcruft/blas/dnrm2.f Fri Jan 17 18:54:19 1997 +0000 +++ b/libcruft/blas/dnrm2.f Sat Jan 18 00:13:08 1997 +0000 @@ -1,122 +1,60 @@ - DOUBLE PRECISION FUNCTION DNRM2 ( N, DX, INCX) - INTEGER NEXT - DOUBLE PRECISION DX(1), CUTLO, CUTHI, HITEST, SUM, XMAX,ZERO,ONE - DATA ZERO, ONE /0.0D0, 1.0D0/ -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 DSQRT(U/EPS) OVER ALL KNOWN MACHINES. -C CUTHI = MINIMUM OF DSQRT(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.30438D19 -C DATA CUTLO, CUTHI / 8.232D-11, 1.304D19 / -C DATA CUTLO, CUTHI / 4.441E-16, 1.304E19 / - DATA CUTLO, CUTHI / 8.232D-11, 1.304D19 / -C - 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( DABS(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( DABS(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 = DABS(DX(I)) - GO TO 115 -C -C PHASE 2. SUM IS SMALL. -C SCALE TO AVOID DESTRUCTIVE UNDERFLOW. -C - 70 IF( DABS(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( DABS(DX(I)) .LE. XMAX ) GO TO 115 - SUM = ONE + SUM * (XMAX / DX(I))**2 - XMAX = DABS(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/FLOAT( N ) -C -C PHASE 3. SUM IS MID-RANGE. NO SCALING. -C - DO 95 J =I,NN,INCX - IF(DABS(DX(J)) .GE. HITEST) GO TO 100 - 95 SUM = SUM + DX(J)**2 - DNRM2 = DSQRT( 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 * DSQRT(SUM) - 300 CONTINUE + DOUBLE PRECISION FUNCTION DNRM2 ( N, X, INCX ) +* .. Scalar Arguments .. + INTEGER INCX, N +* .. Array Arguments .. + DOUBLE PRECISION X( * ) +* .. +* +* DNRM2 returns the euclidean norm of a vector via the function +* name, so that +* +* DNRM2 := sqrt( x'*x ) +* +* +* +* -- This version written on 25-October-1982. +* Modified on 14-October-1993 to inline the call to DLASSQ. +* Sven Hammarling, Nag Ltd. +* +* +* .. Parameters .. + DOUBLE PRECISION ONE , ZERO + PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) +* .. Local Scalars .. + INTEGER IX + DOUBLE PRECISION ABSXI, NORM, SCALE, SSQ +* .. Intrinsic Functions .. + INTRINSIC ABS, SQRT +* .. +* .. Executable Statements .. + IF( N.LT.1 .OR. INCX.LT.1 )THEN + NORM = ZERO + ELSE IF( N.EQ.1 )THEN + NORM = ABS( X( 1 ) ) + ELSE + SCALE = ZERO + SSQ = ONE +* The following loop is equivalent to this call to the LAPACK +* auxiliary routine: +* CALL DLASSQ( N, X, INCX, SCALE, SSQ ) +* + DO 10, IX = 1, 1 + ( N - 1 )*INCX, INCX + IF( X( IX ).NE.ZERO )THEN + ABSXI = ABS( X( IX ) ) + IF( SCALE.LT.ABSXI )THEN + SSQ = ONE + SSQ*( SCALE/ABSXI )**2 + SCALE = ABSXI + ELSE + SSQ = SSQ + ( ABSXI/SCALE )**2 + END IF + END IF + 10 CONTINUE + NORM = SCALE * SQRT( SSQ ) + END IF +* + DNRM2 = NORM RETURN +* +* End of DNRM2. +* END
--- a/libcruft/blas/drot.f Fri Jan 17 18:54:19 1997 +0000 +++ b/libcruft/blas/drot.f Sat Jan 18 00:13:08 1997 +0000 @@ -1,36 +1,37 @@ - SUBROUTINE DROT (N,DX,INCX,DY,INCY,C,S) -C -C APPLIES A PLANE ROTATION. -C JACK DONGARRA, LINPACK, 3/11/78. -C - DOUBLE PRECISION DX(1),DY(1),DTEMP,C,S - INTEGER I,INCX,INCY,IX,IY,N -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 + subroutine drot (n,dx,incx,dy,incy,c,s) +c +c applies a plane rotation. +c jack dongarra, linpack, 3/11/78. +c modified 12/3/93, array(1) declarations changed to array(*) +c + double precision dx(*),dy(*),dtemp,c,s + integer i,incx,incy,ix,iy,n +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
--- a/libcruft/blas/dscal.f Fri Jan 17 18:54:19 1997 +0000 +++ b/libcruft/blas/dscal.f Sat Jan 18 00:13:08 1997 +0000 @@ -1,41 +1,43 @@ - SUBROUTINE DSCAL(N,DA,DX,INCX) -C -C SCALES A VECTOR BY A CONSTANT. -C USES UNROLLED LOOPS FOR INCREMENT EQUAL TO ONE. -C JACK DONGARRA, LINPACK, 3/11/78. -C - DOUBLE PRECISION DA,DX(1) - INTEGER I,INCX,M,MP1,N,NINCX -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 + subroutine dscal(n,da,dx,incx) +c +c scales a vector by a constant. +c uses unrolled loops for increment equal to one. +c jack dongarra, linpack, 3/11/78. +c modified 3/93 to return if incx .le. 0. +c modified 12/3/93, array(1) declarations changed to array(*) +c + double precision da,dx(*) + integer i,incx,m,mp1,n,nincx +c + if( n.le.0 .or. incx.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
--- a/libcruft/blas/dswap.f Fri Jan 17 18:54:19 1997 +0000 +++ b/libcruft/blas/dswap.f Sat Jan 18 00:13:08 1997 +0000 @@ -1,55 +1,56 @@ - SUBROUTINE DSWAP (N,DX,INCX,DY,INCY) -C -C INTERCHANGES TWO VECTORS. -C USES UNROLLED LOOPS FOR INCREMENTS EQUAL ONE. -C JACK DONGARRA, LINPACK, 3/11/78. -C - DOUBLE PRECISION DX(1),DY(1),DTEMP - INTEGER I,INCX,INCY,IX,IY,M,MP1,N -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 + subroutine dswap (n,dx,incx,dy,incy) +c +c interchanges two vectors. +c uses unrolled loops for increments equal one. +c jack dongarra, linpack, 3/11/78. +c modified 12/3/93, array(1) declarations changed to array(*) +c + double precision dx(*),dy(*),dtemp + integer i,incx,incy,ix,iy,m,mp1,n +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
--- a/libcruft/blas/dsyr.f Fri Jan 17 18:54:19 1997 +0000 +++ b/libcruft/blas/dsyr.f Sat Jan 18 00:13:08 1997 +0000 @@ -1,6 +1,3 @@ -* -************************************************************************ -* SUBROUTINE DSYR ( UPLO, N, ALPHA, X, INCX, A, LDA ) * .. Scalar Arguments .. DOUBLE PRECISION ALPHA
--- a/libcruft/blas/dsyrk.f Fri Jan 17 18:54:19 1997 +0000 +++ b/libcruft/blas/dsyrk.f Sat Jan 18 00:13:08 1997 +0000 @@ -1,6 +1,3 @@ -* -************************************************************************ -* SUBROUTINE DSYRK ( UPLO, TRANS, N, K, ALPHA, A, LDA, $ BETA, C, LDC ) * .. Scalar Arguments ..
--- a/libcruft/blas/dtrmm.f Fri Jan 17 18:54:19 1997 +0000 +++ b/libcruft/blas/dtrmm.f Sat Jan 18 00:13:08 1997 +0000 @@ -1,6 +1,3 @@ -* -************************************************************************ -* SUBROUTINE DTRMM ( SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, $ B, LDB ) * .. Scalar Arguments .. @@ -237,7 +234,7 @@ END IF ELSE * -* Form B := alpha*B*A'. +* Form B := alpha*A'*B. * IF( UPPER )THEN DO 110, J = 1, N
--- a/libcruft/blas/dtrmv.f Fri Jan 17 18:54:19 1997 +0000 +++ b/libcruft/blas/dtrmv.f Sat Jan 18 00:13:08 1997 +0000 @@ -1,6 +1,3 @@ -* -************************************************************************ -* SUBROUTINE DTRMV ( UPLO, TRANS, DIAG, N, A, LDA, X, INCX ) * .. Scalar Arguments .. INTEGER INCX, LDA, N
--- a/libcruft/blas/dtrsm.f Fri Jan 17 18:54:19 1997 +0000 +++ b/libcruft/blas/dtrsm.f Sat Jan 18 00:13:08 1997 +0000 @@ -1,6 +1,3 @@ -* -************************************************************************ -* SUBROUTINE DTRSM ( SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, $ B, LDB ) * .. Scalar Arguments ..
--- a/libcruft/blas/dtrsv.f Fri Jan 17 18:54:19 1997 +0000 +++ b/libcruft/blas/dtrsv.f Sat Jan 18 00:13:08 1997 +0000 @@ -1,6 +1,3 @@ -* -************************************************************************ -* SUBROUTINE DTRSV ( UPLO, TRANS, DIAG, N, A, LDA, X, INCX ) * .. Scalar Arguments .. INTEGER INCX, LDA, N
--- a/libcruft/blas/dzasum.f Fri Jan 17 18:54:19 1997 +0000 +++ b/libcruft/blas/dzasum.f Sat Jan 18 00:13:08 1997 +0000 @@ -2,21 +2,21 @@ c c takes the sum of the absolute values. c jack dongarra, 3/11/78. -c modified to correct problem with negative increment, 8/21/90. +c modified 3/93 to return if incx .le. 0. +c modified 12/3/93, array(1) declarations changed to array(*) c - double complex zx(1) + double complex zx(*) double precision stemp,dcabs1 integer i,incx,ix,n c dzasum = 0.0d0 stemp = 0.0d0 - if(n.le.0)return + if( n.le.0 .or. incx.le.0 )return if(incx.eq.1)go to 20 c c code for increment not equal to 1 c ix = 1 - if(incx.lt.0)ix = (-n+1)*incx + 1 do 10 i = 1,n stemp = stemp + dcabs1(zx(ix)) ix = ix + incx
--- a/libcruft/blas/dznrm2.f Fri Jan 17 18:54:19 1997 +0000 +++ b/libcruft/blas/dznrm2.f Sat Jan 18 00:13:08 1997 +0000 @@ -1,138 +1,67 @@ - double precision function dznrm2( n, zx, incx) - logical imag, scale - integer i, incx, ix, n, next - double precision cutlo, cuthi, hitest, sum, xmax, absx, zero, one - double complex zx(1) - double precision dreal,dimag - double complex zdumr,zdumi - dreal(zdumr) = zdumr - dimag(zdumi) = (0.0d0,-1.0d0)*zdumi - data zero, one /0.0d0, 1.0d0/ -c -c unitary norm of the complex n-vector stored in zx() 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 modified to correct problem with negative increment, 8/21/90. -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.30438d19 -c data cutlo, cuthi / 8.232d-11, 1.304d19 / -c data cutlo, cuthi / 4.441e-16, 1.304e19 / - data cutlo, cuthi / 8.232d-11, 1.304d19 / -c - if(n .gt. 0) go to 10 - dznrm2 = zero - go to 300 -c - 10 assign 30 to next - sum = zero - i = 1 - if( incx .lt. 0 )i = (-n+1)*incx + 1 -c begin main loop - do 220 ix = 1,n - absx = dabs(dreal(zx(i))) - imag = .false. - go to next,(30, 50, 70, 90, 110) - 30 if( absx .gt. cutlo) go to 85 - assign 50 to next - scale = .false. -c -c phase 1. sum is zero -c - 50 if( absx .eq. zero) go to 200 - if( absx .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 assign 110 to next - sum = (sum / absx) / absx - 105 scale = .true. - xmax = absx - go to 115 -c -c phase 2. sum is small. -c scale to avoid destructive underflow. -c - 70 if( absx .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( absx .le. xmax ) go to 115 - sum = one + sum * (xmax / absx)**2 - xmax = absx - go to 200 -c - 115 sum = sum + (absx/xmax)**2 - go to 200 -c -c -c prepare for phase 3. -c - 75 sum = (sum * xmax) * xmax -c - 85 assign 90 to next - scale = .false. -c -c for real or d.p. set hitest = cuthi/n -c for complex set hitest = cuthi/(2*n) -c - hitest = cuthi/dble( 2*n ) -c -c phase 3. sum is mid-range. no scaling. -c - 90 if(absx .ge. hitest) go to 100 - sum = sum + absx**2 - 200 continue -c control selection of real and imaginary parts. -c - if(imag) go to 210 - absx = dabs(dimag(zx(i))) - imag = .true. - go to next,( 50, 70, 90, 110 ) -c - 210 continue - i = i + incx - 220 continue -c -c end of main loop. -c compute square root and adjust for scaling. -c - dznrm2 = dsqrt(sum) - if(scale) dznrm2 = dznrm2 * xmax - 300 continue - return - end + DOUBLE PRECISION FUNCTION DZNRM2( N, X, INCX ) +* .. Scalar Arguments .. + INTEGER INCX, N +* .. Array Arguments .. + COMPLEX*16 X( * ) +* .. +* +* DZNRM2 returns the euclidean norm of a vector via the function +* name, so that +* +* DZNRM2 := sqrt( conjg( x' )*x ) +* +* +* +* -- This version written on 25-October-1982. +* Modified on 14-October-1993 to inline the call to ZLASSQ. +* Sven Hammarling, Nag Ltd. +* +* +* .. Parameters .. + DOUBLE PRECISION ONE , ZERO + PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) +* .. Local Scalars .. + INTEGER IX + DOUBLE PRECISION NORM, SCALE, SSQ, TEMP +* .. Intrinsic Functions .. + INTRINSIC ABS, DIMAG, DBLE, SQRT +* .. +* .. Executable Statements .. + IF( N.LT.1 .OR. INCX.LT.1 )THEN + NORM = ZERO + ELSE + SCALE = ZERO + SSQ = ONE +* The following loop is equivalent to this call to the LAPACK +* auxiliary routine: +* CALL ZLASSQ( N, X, INCX, SCALE, SSQ ) +* + DO 10, IX = 1, 1 + ( N - 1 )*INCX, INCX + IF( DBLE( X( IX ) ).NE.ZERO )THEN + TEMP = ABS( DBLE( X( IX ) ) ) + IF( SCALE.LT.TEMP )THEN + SSQ = ONE + SSQ*( SCALE/TEMP )**2 + SCALE = TEMP + ELSE + SSQ = SSQ + ( TEMP/SCALE )**2 + END IF + END IF + IF( DIMAG( X( IX ) ).NE.ZERO )THEN + TEMP = ABS( DIMAG( X( IX ) ) ) + IF( SCALE.LT.TEMP )THEN + SSQ = ONE + SSQ*( SCALE/TEMP )**2 + SCALE = TEMP + ELSE + SSQ = SSQ + ( TEMP/SCALE )**2 + END IF + END IF + 10 CONTINUE + NORM = SCALE * SQRT( SSQ ) + END IF +* + DZNRM2 = NORM + RETURN +* +* End of DZNRM2. +* + END
--- a/libcruft/blas/idamax.f Fri Jan 17 18:54:19 1997 +0000 +++ b/libcruft/blas/idamax.f Sat Jan 18 00:13:08 1997 +0000 @@ -1,37 +1,39 @@ - INTEGER FUNCTION IDAMAX(N,DX,INCX) -C -C FINDS THE INDEX OF ELEMENT HAVING MAX. ABSOLUTE VALUE. -C JACK DONGARRA, LINPACK, 3/11/78. -C - DOUBLE PRECISION DX(1),DMAX - INTEGER I,INCX,IX,N -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 - DMAX = DABS(DX(1)) - IX = IX + INCX - DO 10 I = 2,N - IF(DABS(DX(IX)).LE.DMAX) GO TO 5 - IDAMAX = I - DMAX = DABS(DX(IX)) - 5 IX = IX + INCX - 10 CONTINUE - RETURN -C -C CODE FOR INCREMENT EQUAL TO 1 -C - 20 DMAX = DABS(DX(1)) - DO 30 I = 2,N - IF(DABS(DX(I)).LE.DMAX) GO TO 30 - IDAMAX = I - DMAX = DABS(DX(I)) - 30 CONTINUE - RETURN - END + integer function idamax(n,dx,incx) +c +c finds the index of element having max. absolute value. +c jack dongarra, linpack, 3/11/78. +c modified 3/93 to return if incx .le. 0. +c modified 12/3/93, array(1) declarations changed to array(*) +c + double precision dx(*),dmax + integer i,incx,ix,n +c + idamax = 0 + if( n.lt.1 .or. incx.le.0 ) 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 + dmax = dabs(dx(1)) + ix = ix + incx + do 10 i = 2,n + if(dabs(dx(ix)).le.dmax) go to 5 + idamax = i + dmax = dabs(dx(ix)) + 5 ix = ix + incx + 10 continue + return +c +c code for increment equal to 1 +c + 20 dmax = dabs(dx(1)) + do 30 i = 2,n + if(dabs(dx(i)).le.dmax) go to 30 + idamax = i + dmax = dabs(dx(i)) + 30 continue + return + end
--- a/libcruft/blas/izamax.f Fri Jan 17 18:54:19 1997 +0000 +++ b/libcruft/blas/izamax.f Sat Jan 18 00:13:08 1997 +0000 @@ -2,15 +2,16 @@ c c finds the index of element having max. absolute value. c jack dongarra, 1/15/85. -c modified to correct problem with negative increment, 8/21/90. +c modified 3/93 to return if incx .le. 0. +c modified 12/3/93, array(1) declarations changed to array(*) c - double complex zx(1) + double complex zx(*) double precision smax integer i,incx,ix,n double precision dcabs1 c izamax = 0 - if(n.lt.1)return + if( n.lt.1 .or. incx.le.0 )return izamax = 1 if(n.eq.1)return if(incx.eq.1)go to 20 @@ -18,8 +19,7 @@ c code for increment not equal to 1 c ix = 1 - if(incx.lt.0)ix = (-n+1)*incx + 1 - smax = dcabs1(zx(ix)) + smax = dcabs1(zx(1)) ix = ix + incx do 10 i = 2,n if(dcabs1(zx(ix)).le.smax) go to 5
--- a/libcruft/blas/lsame.f Fri Jan 17 18:54:19 1997 +0000 +++ b/libcruft/blas/lsame.f Sat Jan 18 00:13:08 1997 +0000 @@ -1,89 +1,87 @@ - LOGICAL FUNCTION LSAME ( CA, CB ) + LOGICAL FUNCTION LSAME( CA, CB ) +* +* -- LAPACK auxiliary routine (version 2.0) -- +* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., +* Courant Institute, Argonne National Lab, and Rice University +* January 31, 1994 +* * .. Scalar Arguments .. - CHARACTER*1 CA, CB + CHARACTER CA, CB * .. * * Purpose * ======= * -* LSAME tests if CA is the same letter as CB regardless of case. -* CB is assumed to be an upper case letter. LSAME returns .TRUE. if -* CA is either the same as CB or the equivalent lower case letter. -* -* N.B. This version of the routine is only correct for ASCII code. -* Installers must modify the routine for other character-codes. +* LSAME returns .TRUE. if CA is the same letter as CB regardless of +* case. * -* For EBCDIC systems the constant IOFF must be changed to -64. -* For CDC systems using 6-12 bit representations, the system- -* specific code in comments must be activated. -* -* Parameters -* ========== +* Arguments +* ========= * -* CA - CHARACTER*1 -* CB - CHARACTER*1 -* On entry, CA and CB specify characters to be compared. -* Unchanged on exit. +* CA (input) CHARACTER*1 +* CB (input) CHARACTER*1 +* CA and CB specify the single characters to be compared. * -* -* Auxiliary routine for Level 2 Blas. +* ===================================================================== * -* -- Written on 20-July-1986 -* Richard Hanson, Sandia National Labs. -* Jeremy Du Croz, Nag Central Office. -* -* .. Parameters .. - INTEGER IOFF - PARAMETER ( IOFF=32 ) * .. Intrinsic Functions .. - INTRINSIC ICHAR + INTRINSIC ICHAR +* .. +* .. Local Scalars .. + INTEGER INTA, INTB, ZCODE +* .. * .. Executable Statements .. * * Test if the characters are equal * - LSAME = CA .EQ. CB + LSAME = CA.EQ.CB + IF( LSAME ) + $ RETURN * -* Now test for equivalence +* Now test for equivalence if both characters are alphabetic. +* + ZCODE = ICHAR( 'Z' ) * - IF ( .NOT.LSAME ) THEN - LSAME = ICHAR(CA) - IOFF .EQ. ICHAR(CB) - END IF +* Use 'Z' rather than 'A' so that ASCII can be detected on Prime +* machines, on which ICHAR returns a value with bit 8 set. +* ICHAR('A') on Prime machines returns 193 which is the same as +* ICHAR('A') on an EBCDIC machine. * - RETURN + INTA = ICHAR( CA ) + INTB = ICHAR( CB ) * -* The following comments contain code for CDC systems using 6-12 bit -* representations. + IF( ZCODE.EQ.90 .OR. ZCODE.EQ.122 ) THEN +* +* ASCII is assumed - ZCODE is the ASCII code of either lower or +* upper case 'Z'. * -* .. Parameters .. -* INTEGER ICIRFX -* PARAMETER ( ICIRFX=62 ) -* .. Scalar Arguments .. -* CHARACTER*1 CB -* .. Array Arguments .. -* CHARACTER*1 CA(*) -* .. Local Scalars .. -* INTEGER IVAL -* .. Intrinsic Functions .. -* INTRINSIC ICHAR, CHAR -* .. Executable Statements .. + IF( INTA.GE.97 .AND. INTA.LE.122 ) INTA = INTA - 32 + IF( INTB.GE.97 .AND. INTB.LE.122 ) INTB = INTB - 32 +* + ELSE IF( ZCODE.EQ.233 .OR. ZCODE.EQ.169 ) THEN +* +* EBCDIC is assumed - ZCODE is the EBCDIC code of either lower or +* upper case 'Z'. * -* See if the first character in string CA equals string CB. -* -* LSAME = CA(1) .EQ. CB .AND. CA(1) .NE. CHAR(ICIRFX) -* -* IF (LSAME) RETURN + IF( INTA.GE.129 .AND. INTA.LE.137 .OR. + $ INTA.GE.145 .AND. INTA.LE.153 .OR. + $ INTA.GE.162 .AND. INTA.LE.169 ) INTA = INTA + 64 + IF( INTB.GE.129 .AND. INTB.LE.137 .OR. + $ INTB.GE.145 .AND. INTB.LE.153 .OR. + $ INTB.GE.162 .AND. INTB.LE.169 ) INTB = INTB + 64 * -* The characters are not identical. Now check them for equivalence. -* Look for the 'escape' character, circumflex, followed by the -* letter. + ELSE IF( ZCODE.EQ.218 .OR. ZCODE.EQ.250 ) THEN +* +* ASCII is assumed, on Prime machines - ZCODE is the ASCII code +* plus 128 of either lower or upper case 'Z'. * -* IVAL = ICHAR(CA(2)) -* IF (IVAL.GE.ICHAR('A') .AND. IVAL.LE.ICHAR('Z')) THEN -* LSAME = CA(1) .EQ. CHAR(ICIRFX) .AND. CA(2) .EQ. CB -* END IF + IF( INTA.GE.225 .AND. INTA.LE.250 ) INTA = INTA - 32 + IF( INTB.GE.225 .AND. INTB.LE.250 ) INTB = INTB - 32 + END IF + LSAME = INTA.EQ.INTB * * RETURN * -* End of LSAME. +* End of LSAME * END
--- a/libcruft/blas/sdot.f Fri Jan 17 18:54:19 1997 +0000 +++ b/libcruft/blas/sdot.f Sat Jan 18 00:13:08 1997 +0000 @@ -1,48 +1,49 @@ - REAL FUNCTION SDOT(N,SX,INCX,SY,INCY) -C -C FORMS THE DOT PRODUCT OF TWO VECTORS. -C USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE. -C JACK DONGARRA, LINPACK, 3/11/78. -C - REAL SX(1),SY(1),STEMP - INTEGER I,INCX,INCY,IX,IY,M,MP1,N -C - STEMP = 0.0E0 - SDOT = 0.0E0 - 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 - STEMP = STEMP + SX(IX)*SY(IY) - IX = IX + INCX - IY = IY + INCY - 10 CONTINUE - SDOT = STEMP - 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 - STEMP = STEMP + SX(I)*SY(I) - 30 CONTINUE - IF( N .LT. 5 ) GO TO 60 - 40 MP1 = M + 1 - DO 50 I = MP1,N,5 - STEMP = STEMP + SX(I)*SY(I) + SX(I + 1)*SY(I + 1) + - * SX(I + 2)*SY(I + 2) + SX(I + 3)*SY(I + 3) + SX(I + 4)*SY(I + 4) - 50 CONTINUE - 60 SDOT = STEMP - RETURN - END + real function sdot(n,sx,incx,sy,incy) +c +c forms the dot product of two vectors. +c uses unrolled loops for increments equal to one. +c jack dongarra, linpack, 3/11/78. +c modified 12/3/93, array(1) declarations changed to array(*) +c + real sx(*),sy(*),stemp + integer i,incx,incy,ix,iy,m,mp1,n +c + stemp = 0.0e0 + sdot = 0.0e0 + 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 + stemp = stemp + sx(ix)*sy(iy) + ix = ix + incx + iy = iy + incy + 10 continue + sdot = stemp + 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 + stemp = stemp + sx(i)*sy(i) + 30 continue + if( n .lt. 5 ) go to 60 + 40 mp1 = m + 1 + do 50 i = mp1,n,5 + stemp = stemp + sx(i)*sy(i) + sx(i + 1)*sy(i + 1) + + * sx(i + 2)*sy(i + 2) + sx(i + 3)*sy(i + 3) + sx(i + 4)*sy(i + 4) + 50 continue + 60 sdot = stemp + return + end
--- a/libcruft/blas/xerbla.f Fri Jan 17 18:54:19 1997 +0000 +++ b/libcruft/blas/xerbla.f Sat Jan 18 00:13:08 1997 +0000 @@ -1,45 +1,43 @@ - SUBROUTINE XERBLA ( SRNAME, INFO ) -* .. Scalar Arguments .. + SUBROUTINE XERBLA( SRNAME, INFO ) +* +* -- LAPACK auxiliary routine (preliminary version) -- +* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., +* Courant Institute, Argonne National Lab, and Rice University +* February 29, 1992 +* +* .. Scalar Arguments .. + CHARACTER*6 SRNAME INTEGER INFO - CHARACTER*6 SRNAME * .. * * Purpose * ======= * -* XERBLA is an error handler for the Level 2 BLAS routines. +* XERBLA is an error handler for the LAPACK routines. +* It is called by an LAPACK routine if an input parameter has an +* invalid value. A message is printed and execution stops. * -* It is called by the Level 2 BLAS routines if an input parameter is -* invalid. -* -* Installers should consider modifying the STOP statement in order to +* Installers may consider modifying the STOP statement in order to * call system-specific exception-handling facilities. * -* Parameters -* ========== +* Arguments +* ========= * -* SRNAME - CHARACTER*6. -* On entry, SRNAME specifies the name of the routine which -* called XERBLA. +* SRNAME (input) CHARACTER*6 +* The name of the routine which called XERBLA. * -* INFO - INTEGER. -* On entry, INFO specifies the position of the invalid -* parameter in the parameter-list of the calling routine. +* INFO (input) INTEGER +* The position of the invalid parameter in the parameter list +* of the calling routine. * * -* Auxiliary routine for Level 2 Blas. + WRITE( *, FMT = 9999 )SRNAME, INFO * -* Written on 20-July-1986. -* -* .. Executable Statements .. + STOP * - WRITE (*,99999) SRNAME, INFO -* - CALL XSTOPX (' ') + 9999 FORMAT( ' ** On entry to ', A6, ' parameter number ', I2, ' had ', + $ 'an illegal value' ) * -99999 FORMAT ( ' ** On entry to ', A6, ' parameter number ', I2, - $ ' had an illegal value' ) -* -* End of XERBLA. +* End of XERBLA * END
--- a/libcruft/blas/zaxpy.f Fri Jan 17 18:54:19 1997 +0000 +++ b/libcruft/blas/zaxpy.f Sat Jan 18 00:13:08 1997 +0000 @@ -2,8 +2,9 @@ c c constant times a vector plus a vector. c jack dongarra, 3/11/78. +c modified 12/3/93, array(1) declarations changed to array(*) c - double complex zx(1),zy(1),za + double complex zx(*),zy(*),za integer i,incx,incy,ix,iy,n double precision dcabs1 if(n.le.0)return
--- a/libcruft/blas/zcopy.f Fri Jan 17 18:54:19 1997 +0000 +++ b/libcruft/blas/zcopy.f Sat Jan 18 00:13:08 1997 +0000 @@ -2,8 +2,9 @@ c c copies a vector, x, to a vector, y. c jack dongarra, linpack, 4/11/78. +c modified 12/3/93, array(1) declarations changed to array(*) c - double complex zx(1),zy(1) + double complex zx(*),zy(*) integer i,incx,incy,ix,iy,n c if(n.le.0)return
--- a/libcruft/blas/zdotc.f Fri Jan 17 18:54:19 1997 +0000 +++ b/libcruft/blas/zdotc.f Sat Jan 18 00:13:08 1997 +0000 @@ -2,8 +2,9 @@ c c forms the dot product of a vector. c jack dongarra, 3/11/78. +c modified 12/3/93, array(1) declarations changed to array(*) c - double complex zx(1),zy(1),ztemp + double complex zx(*),zy(*),ztemp integer i,incx,incy,ix,iy,n ztemp = (0.0d0,0.0d0) zdotc = (0.0d0,0.0d0)
--- a/libcruft/blas/zdotu.f Fri Jan 17 18:54:19 1997 +0000 +++ b/libcruft/blas/zdotu.f Sat Jan 18 00:13:08 1997 +0000 @@ -2,8 +2,9 @@ c c forms the dot product of two vectors. c jack dongarra, 3/11/78. +c modified 12/3/93, array(1) declarations changed to array(*) c - double complex zx(1),zy(1),ztemp + double complex zx(*),zy(*),ztemp integer i,incx,incy,ix,iy,n ztemp = (0.0d0,0.0d0) zdotu = (0.0d0,0.0d0)
--- a/libcruft/blas/zdscal.f Fri Jan 17 18:54:19 1997 +0000 +++ b/libcruft/blas/zdscal.f Sat Jan 18 00:13:08 1997 +0000 @@ -2,19 +2,19 @@ c c scales a vector by a constant. c jack dongarra, 3/11/78. -c modified to correct problem with negative increment, 8/21/90. +c modified 3/93 to return if incx .le. 0. +c modified 12/3/93, array(1) declarations changed to array(*) c - double complex zx(1) + double complex zx(*) double precision da integer i,incx,ix,n c - if(n.le.0)return + if( n.le.0 .or. incx.le.0 )return if(incx.eq.1)go to 20 c c code for increment not equal to 1 c ix = 1 - if(incx.lt.0)ix = (-n+1)*incx + 1 do 10 i = 1,n zx(ix) = dcmplx(da,0.0d0)*zx(ix) ix = ix + incx
--- a/libcruft/blas/zherk.f Fri Jan 17 18:54:19 1997 +0000 +++ b/libcruft/blas/zherk.f Sat Jan 18 00:13:08 1997 +0000 @@ -1,9 +1,9 @@ - SUBROUTINE ZHERK ( UPLO, TRANS, N, K, ALPHA, A, LDA, - $ BETA, C, LDC ) + SUBROUTINE ZHERK( UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC ) * .. Scalar Arguments .. - CHARACTER*1 UPLO, TRANS - INTEGER N, K, LDA, LDC + CHARACTER TRANS, UPLO + INTEGER K, LDA, LDC, N DOUBLE PRECISION ALPHA, BETA +* .. * .. Array Arguments .. COMPLEX*16 A( LDA, * ), C( LDC, * ) * .. @@ -61,7 +61,7 @@ * matrix A. K must be at least zero. * Unchanged on exit. * -* ALPHA - DOUBLE PRECISION. +* ALPHA - DOUBLE PRECISION . * On entry, ALPHA specifies the scalar alpha. * Unchanged on exit. * @@ -84,7 +84,7 @@ * On entry, BETA specifies the scalar beta. * Unchanged on exit. * -* C - COMPLEX*16 array of DIMENSION ( LDC, n ). +* C - COMPLEX*16 array of DIMENSION ( LDC, n ). * Before entry with UPLO = 'U' or 'u', the leading n by n * upper triangular part of the array C must contain the upper * triangular part of the hermitian matrix and the strictly @@ -116,28 +116,35 @@ * Jeremy Du Croz, Numerical Algorithms Group Ltd. * Sven Hammarling, Numerical Algorithms Group Ltd. * +* -- Modified 8-Nov-93 to set C(J,J) to DBLE( C(J,J) ) when BETA = 1. +* Ed Anderson, Cray Research Inc. +* * * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME +* .. * .. External Subroutines .. EXTERNAL XERBLA +* .. * .. Intrinsic Functions .. - INTRINSIC DCMPLX, DCONJG, MAX, DBLE + INTRINSIC DBLE, DCMPLX, DCONJG, MAX +* .. * .. Local Scalars .. LOGICAL UPPER INTEGER I, INFO, J, L, NROWA DOUBLE PRECISION RTEMP COMPLEX*16 TEMP +* .. * .. Parameters .. - DOUBLE PRECISION ONE , ZERO - PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) + DOUBLE PRECISION ONE, ZERO + PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) * .. * .. Executable Statements .. * * Test the input parameters. * - IF( LSAME( TRANS, 'N' ) )THEN + IF( LSAME( TRANS, 'N' ) ) THEN NROWA = N ELSE NROWA = K @@ -145,61 +152,59 @@ UPPER = LSAME( UPLO, 'U' ) * INFO = 0 - IF( ( .NOT.UPPER ).AND. - $ ( .NOT.LSAME( UPLO , 'L' ) ) )THEN + IF( ( .NOT.UPPER ) .AND. ( .NOT.LSAME( UPLO, 'L' ) ) ) THEN INFO = 1 - ELSE IF( ( .NOT.LSAME( TRANS, 'N' ) ).AND. - $ ( .NOT.LSAME( TRANS, 'C' ) ) )THEN + ELSE IF( ( .NOT.LSAME( TRANS, 'N' ) ) .AND. + $ ( .NOT.LSAME( TRANS, 'C' ) ) ) THEN INFO = 2 - ELSE IF( N .LT.0 )THEN + ELSE IF( N.LT.0 ) THEN INFO = 3 - ELSE IF( K .LT.0 )THEN + ELSE IF( K.LT.0 ) THEN INFO = 4 - ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN + ELSE IF( LDA.LT.MAX( 1, NROWA ) ) THEN INFO = 7 - ELSE IF( LDC.LT.MAX( 1, N ) )THEN + ELSE IF( LDC.LT.MAX( 1, N ) ) THEN INFO = 10 END IF - IF( INFO.NE.0 )THEN + IF( INFO.NE.0 ) THEN CALL XERBLA( 'ZHERK ', INFO ) RETURN END IF * * Quick return if possible. * - IF( ( N.EQ.0 ).OR. - $ ( ( ( ALPHA.EQ.ZERO ).OR.( K.EQ.0 ) ).AND.( BETA.EQ.ONE ) ) ) - $ RETURN + IF( ( N.EQ.0 ) .OR. ( ( ( ALPHA.EQ.ZERO ) .OR. ( K.EQ.0 ) ) .AND. + $ ( BETA.EQ.ONE ) ) )RETURN * * And when alpha.eq.zero. * - IF( ALPHA.EQ.ZERO )THEN - IF( UPPER )THEN - IF( BETA.EQ.ZERO )THEN - DO 20, J = 1, N - DO 10, I = 1, J + IF( ALPHA.EQ.ZERO ) THEN + IF( UPPER ) THEN + IF( BETA.EQ.ZERO ) THEN + DO 20 J = 1, N + DO 10 I = 1, J C( I, J ) = ZERO 10 CONTINUE 20 CONTINUE ELSE - DO 40, J = 1, N - DO 30, I = 1, J - 1 + DO 40 J = 1, N + DO 30 I = 1, J - 1 C( I, J ) = BETA*C( I, J ) 30 CONTINUE C( J, J ) = BETA*DBLE( C( J, J ) ) 40 CONTINUE END IF ELSE - IF( BETA.EQ.ZERO )THEN - DO 60, J = 1, N - DO 50, I = J, N + IF( BETA.EQ.ZERO ) THEN + DO 60 J = 1, N + DO 50 I = J, N C( I, J ) = ZERO 50 CONTINUE 60 CONTINUE ELSE - DO 80, J = 1, N + DO 80 J = 1, N C( J, J ) = BETA*DBLE( C( J, J ) ) - DO 70, I = J + 1, N + DO 70 I = J + 1, N C( I, J ) = BETA*C( I, J ) 70 CONTINUE 80 CONTINUE @@ -210,51 +215,55 @@ * * Start the operations. * - IF( LSAME( TRANS, 'N' ) )THEN + IF( LSAME( TRANS, 'N' ) ) THEN * * Form C := alpha*A*conjg( A' ) + beta*C. * - IF( UPPER )THEN - DO 130, J = 1, N - IF( BETA.EQ.ZERO )THEN - DO 90, I = 1, J + IF( UPPER ) THEN + DO 130 J = 1, N + IF( BETA.EQ.ZERO ) THEN + DO 90 I = 1, J C( I, J ) = ZERO 90 CONTINUE - ELSE IF( BETA.NE.ONE )THEN - DO 100, I = 1, J - 1 + ELSE IF( BETA.NE.ONE ) THEN + DO 100 I = 1, J - 1 C( I, J ) = BETA*C( I, J ) 100 CONTINUE C( J, J ) = BETA*DBLE( C( J, J ) ) + ELSE + C( J, J ) = DBLE( C( J, J ) ) END IF - DO 120, L = 1, K - IF( A( J, L ).NE.DCMPLX( ZERO ) )THEN + DO 120 L = 1, K + IF( A( J, L ).NE.DCMPLX( ZERO ) ) THEN TEMP = ALPHA*DCONJG( A( J, L ) ) - DO 110, I = 1, J - 1 + DO 110 I = 1, J - 1 C( I, J ) = C( I, J ) + TEMP*A( I, L ) 110 CONTINUE - C( J, J ) = DBLE( C( J, J ) ) + + C( J, J ) = DBLE( C( J, J ) ) + $ DBLE( TEMP*A( I, L ) ) END IF 120 CONTINUE 130 CONTINUE ELSE - DO 180, J = 1, N - IF( BETA.EQ.ZERO )THEN - DO 140, I = J, N + DO 180 J = 1, N + IF( BETA.EQ.ZERO ) THEN + DO 140 I = J, N C( I, J ) = ZERO 140 CONTINUE - ELSE IF( BETA.NE.ONE )THEN + ELSE IF( BETA.NE.ONE ) THEN C( J, J ) = BETA*DBLE( C( J, J ) ) - DO 150, I = J + 1, N + DO 150 I = J + 1, N C( I, J ) = BETA*C( I, J ) 150 CONTINUE + ELSE + C( J, J ) = DBLE( C( J, J ) ) END IF - DO 170, L = 1, K - IF( A( J, L ).NE.DCMPLX( ZERO ) )THEN - TEMP = ALPHA*DCONJG( A( J, L ) ) - C( J, J ) = DBLE( C( J, J ) ) + + DO 170 L = 1, K + IF( A( J, L ).NE.DCMPLX( ZERO ) ) THEN + TEMP = ALPHA*DCONJG( A( J, L ) ) + C( J, J ) = DBLE( C( J, J ) ) + $ DBLE( TEMP*A( J, L ) ) - DO 160, I = J + 1, N + DO 160 I = J + 1, N C( I, J ) = C( I, J ) + TEMP*A( I, L ) 160 CONTINUE END IF @@ -265,46 +274,46 @@ * * Form C := alpha*conjg( A' )*A + beta*C. * - IF( UPPER )THEN - DO 220, J = 1, N - DO 200, I = 1, J - 1 + IF( UPPER ) THEN + DO 220 J = 1, N + DO 200 I = 1, J - 1 TEMP = ZERO - DO 190, L = 1, K + DO 190 L = 1, K TEMP = TEMP + DCONJG( A( L, I ) )*A( L, J ) 190 CONTINUE - IF( BETA.EQ.ZERO )THEN + IF( BETA.EQ.ZERO ) THEN C( I, J ) = ALPHA*TEMP ELSE C( I, J ) = ALPHA*TEMP + BETA*C( I, J ) END IF 200 CONTINUE RTEMP = ZERO - DO 210, L = 1, K + DO 210 L = 1, K RTEMP = RTEMP + DCONJG( A( L, J ) )*A( L, J ) 210 CONTINUE - IF( BETA.EQ.ZERO )THEN + IF( BETA.EQ.ZERO ) THEN C( J, J ) = ALPHA*RTEMP ELSE C( J, J ) = ALPHA*RTEMP + BETA*DBLE( C( J, J ) ) END IF 220 CONTINUE ELSE - DO 260, J = 1, N + DO 260 J = 1, N RTEMP = ZERO - DO 230, L = 1, K + DO 230 L = 1, K RTEMP = RTEMP + DCONJG( A( L, J ) )*A( L, J ) 230 CONTINUE - IF( BETA.EQ.ZERO )THEN + IF( BETA.EQ.ZERO ) THEN C( J, J ) = ALPHA*RTEMP ELSE C( J, J ) = ALPHA*RTEMP + BETA*DBLE( C( J, J ) ) END IF - DO 250, I = J + 1, N + DO 250 I = J + 1, N TEMP = ZERO - DO 240, L = 1, K + DO 240 L = 1, K TEMP = TEMP + DCONJG( A( L, I ) )*A( L, J ) 240 CONTINUE - IF( BETA.EQ.ZERO )THEN + IF( BETA.EQ.ZERO ) THEN C( I, J ) = ALPHA*TEMP ELSE C( I, J ) = ALPHA*TEMP + BETA*C( I, J )
--- a/libcruft/blas/zscal.f Fri Jan 17 18:54:19 1997 +0000 +++ b/libcruft/blas/zscal.f Sat Jan 18 00:13:08 1997 +0000 @@ -2,18 +2,18 @@ c c scales a vector by a constant. c jack dongarra, 3/11/78. -c modified to correct problem with negative increment, 8/21/90. +c modified 3/93 to return if incx .le. 0. +c modified 12/3/93, array(1) declarations changed to array(*) c - double complex za,zx(1) + double complex za,zx(*) integer i,incx,ix,n c - if(n.le.0)return + if( n.le.0 .or. incx.le.0 )return if(incx.eq.1)go to 20 c c code for increment not equal to 1 c ix = 1 - if(incx.lt.0)ix = (-n+1)*incx + 1 do 10 i = 1,n zx(ix) = za*zx(ix) ix = ix + incx
--- a/libcruft/blas/zswap.f Fri Jan 17 18:54:19 1997 +0000 +++ b/libcruft/blas/zswap.f Sat Jan 18 00:13:08 1997 +0000 @@ -2,8 +2,9 @@ c c interchanges two vectors. c jack dongarra, 3/11/78. +c modified 12/3/93, array(1) declarations changed to array(*) c - double complex zx(1),zy(1),ztemp + double complex zx(*),zy(*),ztemp integer i,incx,incy,ix,iy,n c if(n.le.0)return
--- a/libcruft/blas/ztrmm.f Fri Jan 17 18:54:19 1997 +0000 +++ b/libcruft/blas/ztrmm.f Sat Jan 18 00:13:08 1997 +0000 @@ -237,7 +237,7 @@ END IF ELSE * -* Form B := alpha*B*A' or B := alpha*B*conjg( A' ). +* Form B := alpha*A'*B or B := alpha*conjg( A' )*B. * IF( UPPER )THEN DO 120, J = 1, N
--- a/src/ChangeLog Fri Jan 17 18:54:19 1997 +0000 +++ b/src/ChangeLog Sat Jan 18 00:13:08 1997 +0000 @@ -1,3 +1,7 @@ +Wed Jan 8 11:42:44 1997 John W. Eaton <jwe@bevo.che.wisc.edu> + + * log.cc (sqrtm): For complex arg case, compute sqrt, not log. + Tue Jan 7 00:16:41 1997 John W. Eaton <jwe@bevo.che.wisc.edu> * Version 2.0.1 released.