dimension my(12500000),mm(2,20000) character*80 cnfile integer(1) :: a(0:25000000) c c Fixed-format Fortran90 program, c to impute values for saturated pixels (= 65535), c using Principal Components model, c and output rescaled scan image c c input array (NB this assumes a particular style of TIFF image and c may need modifying for other TIFF styles) c write (6,*) 'input: name of tiff input file' read (5,'(a80)') cnfile write (6,*) 'input: file length in bytes (max 25000000)' read (5,*) nbyte nbytez=nbyte-1 call sbrdtf(a,nbytez,my,nc,nr,cnfile) c c input spot centres c write (6,*) 'input: name of file with spot centres' read (5,'(a80)') cnfile write (6,*) 'input: number of spots (max 20000)' read (5,*) nk open(2,file=cnfile) read (2,*) mm(:,1:nk) c c the program assumes kth spot centre is c col location = mm(1,k)/10 + 1 c row location = mm(2,k)/10 + 1 c write (6,*) 'input: spot radius (max 12)' read (5,*) nrad c c impute saturated values c call sbfill(my,nc,nr,mm,nk,nrad) c c output array (NB this assumes a particular style of TIFF image and c may need modifying for other TIFF styles) c write (6,*) 'input: name of tiff output file' read (5,'(a80)') cnfile call sbout(a,nbytez,my,nc,nr,cnfile) c stop end ccccccccc*ccccccccc*ccccccccc*ccccccccc*ccccccccc*ccccccccc*ccccccccc*cc subroutine sbrdtf(a,nbytez,m,nc,nr,cnfile) dimension m(12500000) character*80 cnfile integer(1) :: a(0:nbytez) c c read microarray-formatted tiff file c =================================== c c reads nc * nr data into m(12500000) from file "cnfile" c nbyte=nbytez+1 open(unit=2,form='unformatted',recl=nbyte,access='direct', * status='old',file=cnfile) read(2,rec=1) a c n=((itrans(a(7))*256+itrans(a(6)))*256+itrans(a(5)))*256 * +itrans(a(4)) write(6,'(10i4)') (itrans(a(i)),i=n+10,n+23) nc=itrans(a(n+11))*256+itrans(a(n+10)) nr=itrans(a(n+23))*256+itrans(a(n+22)) c j=8 do 1 i=1,nc*nr m(i)=itrans(a(j+1))*256+itrans(a(j)) 1 j=j+2 return end function itrans(b) integer(1):: b integer it it=b if(b.lt.0) it=it+256 itrans=it return end ccccccccc*ccccccccc*ccccccccc*ccccccccc*ccccccccc*ccccccccc*ccccccccc*cc subroutine sbout(a,nbytez,my,nc,nr,cnfile) dimension my(nc,nr) character*80 cnfile integer(1) :: a(0:nbytez) c c output tiff file c mymax=maxval(my) ymult=(65535-0.5)/mymax k=8 do i=1,nr do j=1,nc ky=nint(my(j,i)*ymult) kk=ky/256 a(k)=ky-kk*256 a(k+1)=kk k=k+2 enddo enddo c nbyte=nbytez+1 open(unit=2,form='unformatted',recl=nbyte,access='direct', * file=cnfile) write(2,rec=1) a c return end ccccccccc*ccccccccc*ccccccccc*ccccccccc*ccccccccc*ccccccccc*ccccccccc*cc subroutine sbfill(my,nc,nr,mm,nk,nrad) dimension my(nc,nr),mm(2,nk),xm(625),v(625,625), * va(625),vb(625,625),mymax(20000),ms(20000), * x(100),xx(100,100),xxlamb(100,100),b(100),yss(10), * mloc1(1),w1(625),w2(625) c nl=(2*nrad+1)**2 nvv=625 npc=min(100,nl/2) nxx=100 nlamb=10 ylim=65535-0.5 c c rank spots by max value c do k=1,nk ic=mm(2,k)/10+1 jc=mm(1,k)/10+1 mymax(k)=maxval(my(jc-nrad:jc+nrad,ic-nrad:ic+nrad)) ms(k)=k enddo call isort(mymax,ms,nk,2) ccc write (6,*) mymax(1),mymax(nk) c c mean spot shape c xm=0 do kk=1,nk if (mymax(kk).lt.ylim) then ns=kk+1 k=ms(kk) ic=mm(2,k)/10+1 jc=mm(1,k)/10+1 l=0 do i=-nrad,nrad do j=-nrad,nrad l=l+1 xm(l)=xm(l)+my(jc+j,ic+i) enddo enddo endif enddo xm=xm/(ns-1) c c variance matrix c v=0 do kk=1,ns-1 k=ms(kk) ic=mm(2,k)/10+1 jc=mm(1,k)/10+1 l=0 do i=-nrad,nrad do j=-nrad,nrad l=l+1 ll=0 do ii=-nrad,nrad do jj=-nrad,nrad ll=ll+1 v(ll,l)=v(ll,l)+(my(jc+j,ic+i)-xm(l))* * (my(jc+jj,ic+ii)-xm(ll)) enddo enddo enddo enddo enddo v=v/(ns-2) c c eigenvectors and eigenvalues c call rs(nvv,nl,v,va,1,vb,w1,w2,ierr) c do n=1,npc nn=nl+1-n vaa=sqrt(va(nn)) do l=1,nl vb(l,n)=vb(l,nn)*vaa enddo enddo c c choose xlamb c nxval=max(100,nk-ns+1) ylima=mymax(ns-nxval-1) ccc write (6,*) ns,nxval,ylima yss=0 do kk=ns-nxval,ns-1 k=ms(kk) x=0 xx=0 ic=mm(2,k)/10+1 jc=mm(1,k)/10+1 l=0 do i=-nrad,nrad do j=-nrad,nrad l=l+1 if (my(jc+j,ic+i).le.ylima) then do n=1,npc x(n)=x(n)+(my(jc+j,ic+i)-xm(l))*vb(l,n) do nn=1,npc xx(n,nn)=xx(n,nn)+vb(l,n)*vb(l,nn) enddo enddo endif enddo enddo c do ilamb=1,nlamb xlamb=10.0**(ilamb-2) xxlamb=xx do n=1,npc xxlamb(n,n)=xxlamb(n,n)+xlamb enddo b=x call spoir(xxlamb,nxx,npc,b,1,ind,v) c ys=0 l=0 do i=-nrad,nrad do j=-nrad,nrad l=l+1 if (my(jc+j,ic+i).gt.ylima) then yy=xm(l) do n=1,npc yy=yy+b(n)*vb(l,n) enddo yy=max(yy,ylima) ys=ys+my(jc+j,ic+i)-yy endif enddo enddo yss(ilamb)=yss(ilamb)+(ys/nl)**2 enddo enddo ccc write (6,200) sqrt(yss/nxval) 200 format (10f8.1) c mloc1=minloc(yss) c c impute censored spots c yst=0 do kk=ns,nk k=ms(kk) x=0 xx=0 ic=mm(2,k)/10+1 jc=mm(1,k)/10+1 l=0 do i=-nrad,nrad do j=-nrad,nrad l=l+1 if (my(jc+j,ic+i).le.ylim) then do n=1,npc x(n)=x(n)+(my(jc+j,ic+i)-xm(l))*vb(l,n) do nn=1,npc xx(n,nn)=xx(n,nn)+vb(l,n)*vb(l,nn) enddo enddo endif enddo enddo c xlamb=10.0**(mloc1(1)-2) xxlamb=xx do n=1,npc xxlamb(n,n)=xxlamb(n,n)+xlamb enddo b=x call spoir(xxlamb,nxx,npc,b,1,ind,v) c ys=0 l=0 do i=-nrad,nrad do j=-nrad,nrad l=l+1 if (my(jc+j,ic+i).gt.ylim) then yy=xm(l) do n=1,npc yy=yy+b(n)*vb(l,n) enddo yy=max(yy,ylim) ys=ys+my(jc+j,ic+i)-yy my(jc+j,ic+i)=nint(yy) endif enddo enddo yst=yst+(ys/nl)**2 enddo c ccc write (6,200) xlamb,sqrt(yst/(nk-ns+1)) return end ccccccccc*ccccccccc*ccccccccc*ccccccccc*ccccccccc*ccccccccc*ccccccccc*cc SUBROUTINE ISORT(X,Y,N,KFLAG) C***BEGIN PROLOGUE ISORT C***DATE WRITTEN 761118 (YYMMDD) C***REVISION DATE 820801 (YYMMDD) C***CATEGORY NO. N6A2A C***KEYWORDS QUICKSORT,SINGLETON QUICKSORT,SORT,SORTING C***AUTHOR JONES, R. E., (SNLA) C KAHANER, D. K., (NBS) C WISNIEWSKI, J. A., (SNLA) C***PURPOSE ISORT sorts integer array X and optionally makes the same C interchanges in integer array Y. The array X may be C sorted in increasing order or decreasing order. A C slightly modified QUICKSORT algorithm is used. C***DESCRIPTION C C Written by Rondall E Jones C Modified by John A. Wisniewski to use the Singleton QUICKSORT C algorithm. Date 18 November 1976. C C Further modified by David K. Kahaner C NATIONAL BUREAU OF STANDARDS C August, 1981 C C Abstract C ISORT sorts integer array X and optionally makes the same C interchanges in integer array Y. The array X may be sorted in C INCREASING order or DECREASING order. A slightly modified C QUICKSORT algorithm is used. C C Reference C Singleton,R.C., Algorithm 347, An Efficient Algorithm For C Sorting With Minimal Storage, CACM,12(3),1969,185-7. C C Description of Parameters C X - integer array of values to be sorted C Y - integer array to be (optionally) carried along C N - number of values in integer array X to be sorted C KFLAG - control parameter C = 2 means sort X in INCREASING order and carry Y along. C = 1 means sort X in INCREASING order (ignoring Y) C =-1 means sort X in DECREASING order (ignoring Y) C =-2 means sort X in DECREASING order and carry Y along. C***REFERENCES SINGLETON, R. C., ALGORITHM 347, AN EFFICIENT C ALGORITHM FOR SORTING WITH MINIMAL STORAGE, CACM, C VOL. 12, NO. 3, 1969, PP. 185-187. C***ROUTINES CALLED XERROR C***END PROLOGUE ISORT DIMENSION IL(21),IU(21) INTEGER X(N),Y(N),T,TT,TY,TTY C***FIRST EXECUTABLE STATEMENT ISORT NN = N IF (NN.GE.1) GO TO 10 CALL XERROR ( 'ISORT- THE NUMBER OF VALUES TO BE SORTED WAS NOT PO 1SITIVE.',58,1,1) RETURN 10 KK = IABS(KFLAG) IF ((KK.EQ.1).OR.(KK.EQ.2)) GO TO 15 CALL XERROR ( 'ISORT- THE SORT CONTROL PARAMETER, K, WAS NOT 2, 1, 1 -1, OR -2.',62,2,1) RETURN C C ALTER ARRAY X TO GET DECREASING ORDER IF NEEDED C 15 IF (KFLAG.GE.1) GO TO 30 DO 20 I=1,NN 20 X(I) = -X(I) 30 GO TO (100,200),KK C C SORT X ONLY C 100 CONTINUE M=1 I=1 J=NN R=.375 110 IF (I .EQ. J) GO TO 155 115 IF (R .GT. .5898437) GO TO 120 R=R+3.90625E-2 GO TO 125 120 R=R-.21875 125 K=I C SELECT A CENTRAL ELEMENT OF THE C ARRAY AND SAVE IT IN LOCATION T IJ = I + IFIX (FLOAT (J-I) * R) T=X(IJ) C IF FIRST ELEMENT OF ARRAY IS GREATER C THAN T, INTERCHANGE WITH T IF (X(I) .LE. T) GO TO 130 X(IJ)=X(I) X(I)=T T=X(IJ) 130 L=J C IF LAST ELEMENT OF ARRAY IS LESS THAN C T, INTERCHANGE WITH T IF (X(J) .GE. T) GO TO 140 X(IJ)=X(J) X(J)=T T=X(IJ) C IF FIRST ELEMENT OF ARRAY IS GREATER C THAN T, INTERCHANGE WITH T IF (X(I) .LE. T) GO TO 140 X(IJ)=X(I) X(I)=T T=X(IJ) GO TO 140 135 TT=X(L) X(L)=X(K) X(K)=TT C FIND AN ELEMENT IN THE SECOND HALF OF C THE ARRAY WHICH IS SMALLER THAN T 140 L=L-1 IF (X(L) .GT. T) GO TO 140 C FIND AN ELEMENT IN THE FIRST HALF OF C THE ARRAY WHICH IS GREATER THAN T 145 K=K+1 IF (X(K) .LT. T) GO TO 145 C INTERCHANGE THESE ELEMENTS IF (K .LE. L) GO TO 135 C SAVE UPPER AND LOWER SUBSCRIPTS OF C THE ARRAY YET TO BE SORTED IF (L-I .LE. J-K) GO TO 150 IL(M)=I IU(M)=L I=K M=M+1 GO TO 160 150 IL(M)=K IU(M)=J J=L M=M+1 GO TO 160 C BEGIN AGAIN ON ANOTHER PORTION OF C THE UNSORTED ARRAY 155 M=M-1 IF (M .EQ. 0) GO TO 300 I=IL(M) J=IU(M) 160 IF (J-I .GE. 1) GO TO 125 IF (I .EQ. 1) GO TO 110 I=I-1 165 I=I+1 IF (I .EQ. J) GO TO 155 T=X(I+1) IF (X(I) .LE. T) GO TO 165 K=I 170 X(K+1)=X(K) K=K-1 IF (T .LT. X(K)) GO TO 170 X(K+1)=T GO TO 165 C C SORT X AND CARRY Y ALONG C 200 CONTINUE M=1 I=1 J=NN R=.375 210 IF (I .EQ. J) GO TO 255 215 IF (R .GT. .5898437) GO TO 220 R=R+3.90625E-2 GO TO 225 220 R=R-.21875 225 K=I C SELECT A CENTRAL ELEMENT OF THE C ARRAY AND SAVE IT IN LOCATION T IJ = I + IFIX (FLOAT (J-I) *R) T=X(IJ) TY= Y(IJ) C IF FIRST ELEMENT OF ARRAY IS GREATER C THAN T, INTERCHANGE WITH T IF (X(I) .LE. T) GO TO 230 X(IJ)=X(I) X(I)=T T=X(IJ) Y(IJ)= Y(I) Y(I)=TY TY= Y(IJ) 230 L=J C IF LAST ELEMENT OF ARRAY IS LESS THAN C T, INTERCHANGE WITH T IF (X(J) .GE. T) GO TO 240 X(IJ)=X(J) X(J)=T T=X(IJ) Y(IJ)= Y(J) Y(J)=TY TY= Y(IJ) C IF FIRST ELEMENT OF ARRAY IS GREATER C THAN T, INTERCHANGE WITH T IF (X(I) .LE. T) GO TO 240 X(IJ)=X(I) X(I)=T T=X(IJ) Y(IJ)= Y(I) Y(I)=TY TY= Y(IJ) GO TO 240 235 TT=X(L) X(L)=X(K) X(K)=TT TTY= Y(L) Y(L)= Y(K) Y(K)=TTY C FIND AN ELEMENT IN THE SECOND HALF OF C THE ARRAY WHICH IS SMALLER THAN T 240 L=L-1 IF (X(L) .GT. T) GO TO 240 C FIND AN ELEMENT IN THE FIRST HALF OF C THE ARRAY WHICH IS GREATER THAN T 245 K=K+1 IF (X(K) .LT. T) GO TO 245 C INTERCHANGE THESE ELEMENTS IF (K .LE. L) GO TO 235 C SAVE UPPER AND LOWER SUBSCRIPTS OF C THE ARRAY YET TO BE SORTED IF (L-I .LE. J-K) GO TO 250 IL(M)=I IU(M)=L I=K M=M+1 GO TO 260 250 IL(M)=K IU(M)=J J=L M=M+1 GO TO 260 C BEGIN AGAIN ON ANOTHER PORTION OF C THE UNSORTED ARRAY 255 M=M-1 IF (M .EQ. 0) GO TO 300 I=IL(M) J=IU(M) 260 IF (J-I .GE. 1) GO TO 225 IF (I .EQ. 1) GO TO 210 I=I-1 265 I=I+1 IF (I .EQ. J) GO TO 255 T=X(I+1) TY= Y(I+1) IF (X(I) .LE. T) GO TO 265 K=I 270 X(K+1)=X(K) Y(K+1)= Y(K) K=K-1 IF (T .LT. X(K)) GO TO 270 X(K+1)=T Y(K+1)=TY GO TO 265 C C CLEAN UP C 300 IF (KFLAG.GE.1) RETURN DO 310 I=1,NN 310 X(I) = -X(I) RETURN END SUBROUTINE RS(NM,N,A,W,MATZ,Z,FV1,FV2,IERR) C***BEGIN PROLOGUE RS C***DATE WRITTEN 760101 (YYMMDD) C***REVISION DATE 830518 (YYMMDD) C***CATEGORY NO. D4A1 C***KEYWORDS EIGENVALUES,EIGENVECTORS,EISPACK C***AUTHOR SMITH, B. T., ET AL. C***PURPOSE Computes eigenvalues and, optionally, eigenvectors of C real symmetric matrix. C***DESCRIPTION C C This subroutine calls the recommended sequence of C subroutines from the eigensystem subroutine package (EISPACK) C to find the eigenvalues and eigenvectors (if desired) C of a REAL SYMMETRIC matrix. C C On Input C C NM must be set to the row dimension of the two-dimensional C array parameters as declared in the calling program C dimension statement. C C N is the order of the matrix A. C C A contains the real symmetric matrix. C C MATZ is an integer variable set equal to zero if C only eigenvalues are desired. Otherwise it is set to C any non-zero integer for both eigenvalues and eigenvectors. C C On Output C C W contains the eigenvalues in ascending order. C C Z contains the eigenvectors if MATZ is not zero. C C IERR is an integer output variable set equal to an C error completion code described in section 2B of the C documentation. The normal completion code is zero. C C FV1 and FV2 are temporary storage arrays. C C Questions and comments should be directed to B. S. Garbow, C APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY C ------------------------------------------------------------------ C***REFERENCES B. T. SMITH, J. M. BOYLE, J. J. DONGARRA, B. S. GARBOW, C Y. IKEBE, V. C. KLEMA, C. B. MOLER, *MATRIX EIGEN- C SYSTEM ROUTINES - EISPACK GUIDE*, SPRINGER-VERLAG, C 1976. C***ROUTINES CALLED TQL2,TQLRAT,TRED1,TRED2 C***END PROLOGUE RS C INTEGER N,NM,IERR,MATZ REAL A(NM,N),W(N),Z(NM,N),FV1(N),FV2(N) C C***FIRST EXECUTABLE STATEMENT RS IF (N .LE. NM) GO TO 10 IERR = 10 * N GO TO 50 C 10 IF (MATZ .NE. 0) GO TO 20 C .......... FIND EIGENVALUES ONLY .......... CALL TRED1(NM,N,A,W,FV1,FV2) CALL TQLRAT(N,W,FV2,IERR) GO TO 50 C .......... FIND BOTH EIGENVALUES AND EIGENVECTORS .......... 20 CALL TRED2(NM,N,A,W,FV1,Z) CALL TQL2(NM,N,W,FV1,Z,IERR) 50 RETURN END SUBROUTINE TQL2(NM,N,D,E,Z,IERR) C***BEGIN PROLOGUE TQL2 C***DATE WRITTEN 760101 (YYMMDD) C***REVISION DATE 830518 (YYMMDD) C***CATEGORY NO. D4A5,D4C2A C***KEYWORDS EIGENVALUES,EIGENVECTORS,EISPACK C***AUTHOR SMITH, B. T., ET AL. C***PURPOSE Compute eigenvalues and eigenvectors of symmetric C tridiagonal matrix. C***DESCRIPTION C C This subroutine is a translation of the ALGOL procedure TQL2, C NUM. MATH. 11, 293-306(1968) by Bowdler, Martin, Reinsch, and C Wilkinson. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 227-240(1971). C C This subroutine finds the eigenvalues and eigenvectors C of a SYMMETRIC TRIDIAGONAL matrix by the QL method. C The eigenvectors of a FULL SYMMETRIC matrix can also C be found if TRED2 has been used to reduce this C full matrix to tridiagonal form. C C On Input C C NM must be set to the row dimension of two-dimensional C array parameters as declared in the calling program C dimension statement. C C N is the order of the matrix. C C D contains the diagonal elements of the input matrix. C C E contains the subdiagonal elements of the input matrix C in its last N-1 positions. E(1) is arbitrary. C C Z contains the transformation matrix produced in the C reduction by TRED2, if performed. If the eigenvectors C of the tridiagonal matrix are desired, Z must contain C the identity matrix. C C On Output C C D contains the eigenvalues in ascending order. If an C error exit is made, the eigenvalues are correct but C unordered for indices 1,2,...,IERR-1. C C E has been destroyed. C C Z contains orthonormal eigenvectors of the symmetric C tridiagonal (or full) matrix. If an error exit is made, C Z contains the eigenvectors associated with the stored C eigenvalues. C C IERR is set to C Zero for normal return, C J if the J-th eigenvalue has not been C determined after 30 iterations. C C Calls PYTHAG(A,B) for sqrt(A**2 + B**2). C C Questions and comments should be directed to B. S. Garbow, C APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY C ------------------------------------------------------------------ C***REFERENCES B. T. SMITH, J. M. BOYLE, J. J. DONGARRA, B. S. GARBOW, C Y. IKEBE, V. C. KLEMA, C. B. MOLER, *MATRIX EIGEN- C SYSTEM ROUTINES - EISPACK GUIDE*, SPRINGER-VERLAG, C 1976. C***ROUTINES CALLED PYTHAG C***END PROLOGUE TQL2 C INTEGER I,J,K,L,M,N,II,L1,L2,NM,MML,IERR REAL D(N),E(N),Z(NM,N) REAL B,C,C2,C3,DL1,EL1,F,G,H,P,R,S,S2 REAL PYTHAG C C***FIRST EXECUTABLE STATEMENT TQL2 IERR = 0 IF (N .EQ. 1) GO TO 1001 C DO 100 I = 2, N 100 E(I-1) = E(I) C F = 0.0E0 B = 0.0E0 E(N) = 0.0E0 C DO 240 L = 1, N J = 0 H = ABS(D(L)) + ABS(E(L)) IF (B .LT. H) B = H C .......... LOOK FOR SMALL SUB-DIAGONAL ELEMENT .......... DO 110 M = L, N IF (B + ABS(E(M)) .EQ. B) GO TO 120 C .......... E(N) IS ALWAYS ZERO, SO THERE IS NO EXIT C THROUGH THE BOTTOM OF THE LOOP .......... 110 CONTINUE C 120 IF (M .EQ. L) GO TO 220 130 IF (J .EQ. 30) GO TO 1000 J = J + 1 C .......... FORM SHIFT .......... L1 = L + 1 L2 = L1 + 1 G = D(L) P = (D(L1) - G) / (2.0E0 * E(L)) R = PYTHAG(P,1.0E0) D(L) = E(L) / (P + SIGN(R,P)) D(L1) = E(L) * (P + SIGN(R,P)) DL1 = D(L1) H = G - D(L) IF (L2 .GT. N) GO TO 145 C DO 140 I = L2, N 140 D(I) = D(I) - H C 145 F = F + H C .......... QL TRANSFORMATION .......... P = D(M) C = 1.0E0 C2 = C EL1 = E(L1) S = 0.0E0 MML = M - L C .......... FOR I=M-1 STEP -1 UNTIL L DO -- .......... DO 200 II = 1, MML C3 = C2 C2 = C S2 = S I = M - II G = C * E(I) H = C * P IF (ABS(P) .LT. ABS(E(I))) GO TO 150 C = E(I) / P R = SQRT(C*C+1.0E0) E(I+1) = S * P * R S = C / R C = 1.0E0 / R GO TO 160 150 C = P / E(I) R = SQRT(C*C+1.0E0) E(I+1) = S * E(I) * R S = 1.0E0 / R C = C * S 160 P = C * D(I) - S * G D(I+1) = H + S * (C * G + S * D(I)) C .......... FORM VECTOR .......... DO 180 K = 1, N H = Z(K,I+1) Z(K,I+1) = S * Z(K,I) + C * H Z(K,I) = C * Z(K,I) - S * H 180 CONTINUE C 200 CONTINUE C P = -S * S2 * C3 * EL1 * E(L) / DL1 E(L) = S * P D(L) = C * P IF (B + ABS(E(L)) .GT. B) GO TO 130 220 D(L) = D(L) + F 240 CONTINUE C .......... ORDER EIGENVALUES AND EIGENVECTORS .......... DO 300 II = 2, N I = II - 1 K = I P = D(I) C DO 260 J = II, N IF (D(J) .GE. P) GO TO 260 K = J P = D(J) 260 CONTINUE C IF (K .EQ. I) GO TO 300 D(K) = D(I) D(I) = P C DO 280 J = 1, N P = Z(J,I) Z(J,I) = Z(J,K) Z(J,K) = P 280 CONTINUE C 300 CONTINUE C GO TO 1001 C .......... SET ERROR -- NO CONVERGENCE TO AN C EIGENVALUE AFTER 30 ITERATIONS .......... 1000 IERR = L 1001 RETURN END SUBROUTINE TQLRAT(N,D,E2,IERR) C***BEGIN PROLOGUE TQLRAT C***DATE WRITTEN 760101 (YYMMDD) C***REVISION DATE 830518 (YYMMDD) C***CATEGORY NO. D4A5,D4C2A C***KEYWORDS EIGENVALUES,EIGENVECTORS,EISPACK C***AUTHOR SMITH, B. T., ET AL. C***PURPOSE Computes eigenvalues of symmetric tridiagonal matrix C a rational variant of the QL method. C***DESCRIPTION C C This subroutine is a translation of the ALGOL procedure TQLRAT, C ALGORITHM 464, COMM. ACM 16, 689(1973) by Reinsch. C C This subroutine finds the eigenvalues of a SYMMETRIC C TRIDIAGONAL matrix by the rational QL method. C C On Input C C N is the order of the matrix. C C D contains the diagonal elements of the input matrix. C C E2 contains the squares of the subdiagonal elements of the C input matrix in its last N-1 positions. E2(1) is arbitrary. C C On Output C C D contains the eigenvalues in ascending order. If an C error exit is made, the eigenvalues are correct and C ordered for indices 1,2,...IERR-1, but may not be C the smallest eigenvalues. C C E2 has been destroyed. C C IERR is set to C Zero for normal return, C J if the J-th eigenvalue has not been C determined after 30 iterations. C C Calls PYTHAG(A,B) for sqrt(A**2 + B**2). C C Questions and comments should be directed to B. S. Garbow, C APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY C ------------------------------------------------------------------ C***REFERENCES B. T. SMITH, J. M. BOYLE, J. J. DONGARRA, B. S. GARBOW, C Y. IKEBE, V. C. KLEMA, C. B. MOLER, *MATRIX EIGEN- C SYSTEM ROUTINES - EISPACK GUIDE*, SPRINGER-VERLAG, C 1976. C***ROUTINES CALLED PYTHAG C***END PROLOGUE TQLRAT C INTEGER I,J,L,M,N,II,L1,MML,IERR REAL D(N),E2(N) REAL B,C,F,G,H,P,R,S,MACHEP REAL PYTHAG C DATA MACHEP/1.0E0/ C***FIRST EXECUTABLE STATEMENT TQLRAT IF (MACHEP .NE. 1.0E0) GO TO 10 C C --- This code fails to compute MACHEP correctly on IBM machines. --- C --- Replaced by call to R1MACH on 15 Jun 94 by Ron Boisvert. --- C C 05 MACHEP = 0.5E0*MACHEP C IF (1.0E0 + MACHEP .GT. 1.0E0) GO TO 05 C MACHEP = 2.0E0*MACHEP C MACHEP = R1MACH(4) C 10 IERR = 0 IF (N .EQ. 1) GO TO 1001 C DO 100 I = 2, N 100 E2(I-1) = E2(I) C F = 0.0E0 B = 0.0E0 E2(N) = 0.0E0 C DO 290 L = 1, N J = 0 H = MACHEP * (ABS(D(L)) + SQRT(E2(L))) IF (B .GT. H) GO TO 105 B = H C = B * B C .......... LOOK FOR SMALL SQUARED SUB-DIAGONAL ELEMENT .......... 105 DO 110 M = L, N IF (E2(M) .LE. C) GO TO 120 C .......... E2(N) IS ALWAYS ZERO, SO THERE IS NO EXIT C THROUGH THE BOTTOM OF THE LOOP .......... 110 CONTINUE C 120 IF (M .EQ. L) GO TO 210 130 IF (J .EQ. 30) GO TO 1000 J = J + 1 C .......... FORM SHIFT .......... L1 = L + 1 S = SQRT(E2(L)) G = D(L) P = (D(L1) - G) / (2.0E0 * S) R = PYTHAG(P,1.0E0) D(L) = S / (P + SIGN(R,P)) H = G - D(L) C DO 140 I = L1, N 140 D(I) = D(I) - H C F = F + H C .......... RATIONAL QL TRANSFORMATION .......... G = D(M) IF (G .EQ. 0.0E0) G = B H = G S = 0.0E0 MML = M - L C .......... FOR I=M-1 STEP -1 UNTIL L DO -- .......... DO 200 II = 1, MML I = M - II P = G * H R = P + E2(I) E2(I+1) = S * R S = E2(I) / R D(I+1) = H + S * (H + D(I)) G = D(I) - E2(I) / G IF (G .EQ. 0.0E0) G = B H = G * P / R 200 CONTINUE C E2(L) = S * G D(L) = H C .......... GUARD AGAINST UNDERFLOW IN CONVERGENCE TEST .......... IF (H .EQ. 0.0E0) GO TO 210 IF (ABS(E2(L)) .LE. ABS(C/H)) GO TO 210 E2(L) = H * E2(L) IF (E2(L) .NE. 0.0E0) GO TO 130 210 P = D(L) + F C .......... ORDER EIGENVALUES .......... IF (L .EQ. 1) GO TO 250 C .......... FOR I=L STEP -1 UNTIL 2 DO -- .......... DO 230 II = 2, L I = L + 2 - II IF (P .GE. D(I-1)) GO TO 270 D(I) = D(I-1) 230 CONTINUE C 250 I = 1 270 D(I) = P 290 CONTINUE C GO TO 1001 C .......... SET ERROR -- NO CONVERGENCE TO AN C EIGENVALUE AFTER 30 ITERATIONS .......... 1000 IERR = L 1001 RETURN END SUBROUTINE TRED1(NM,N,A,D,E,E2) C***BEGIN PROLOGUE TRED1 C***DATE WRITTEN 760101 (YYMMDD) C***REVISION DATE 830518 (YYMMDD) C***CATEGORY NO. D4C1B1 C***KEYWORDS EIGENVALUES,EIGENVECTORS,EISPACK C***AUTHOR SMITH, B. T., ET AL. C***PURPOSE Reduce real symmetric matrix to symmetric tridiagonal C matrix using orthogonal similarity transformations. C***DESCRIPTION C C This subroutine is a translation of the ALGOL procedure TRED1, C NUM. MATH. 11, 181-195(1968) by Martin, Reinsch, and Wilkinson. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971). C C This subroutine reduces a REAL SYMMETRIC matrix C to a symmetric tridiagonal matrix using C orthogonal similarity transformations. C C On Input C C NM must be set to the row dimension of two-dimensional C array parameters as declared in the calling program C dimension statement. C C N is the order of the matrix. C C A contains the real symmetric input matrix. Only the C lower triangle of the matrix need be supplied. C C On Output C C A contains information about the orthogonal trans- C formations used in the reduction in its strict lower C triangle. The full upper triangle of A is unaltered. C C D contains the diagonal elements of the tridiagonal matrix. C C E contains the subdiagonal elements of the tridiagonal C matrix in its last N-1 positions. E(1) is set to zero. C C E2 contains the squares of the corresponding elements of E. C E2 may coincide with E if the squares are not needed. C C Questions and comments should be directed to B. S. Garbow, C APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY C ------------------------------------------------------------------ C***REFERENCES B. T. SMITH, J. M. BOYLE, J. J. DONGARRA, B. S. GARBOW, C Y. IKEBE, V. C. KLEMA, C. B. MOLER, *MATRIX EIGEN- C SYSTEM ROUTINES - EISPACK GUIDE*, SPRINGER-VERLAG, C 1976. C***ROUTINES CALLED (NONE) C***END PROLOGUE TRED1 C INTEGER I,J,K,L,N,II,NM,JP1 REAL A(NM,N),D(N),E(N),E2(N) REAL F,G,H,SCALE C C***FIRST EXECUTABLE STATEMENT TRED1 DO 100 I = 1, N 100 D(I) = A(I,I) C .......... FOR I=N STEP -1 UNTIL 1 DO -- .......... DO 300 II = 1, N I = N + 1 - II L = I - 1 H = 0.0E0 SCALE = 0.0E0 IF (L .LT. 1) GO TO 130 C .......... SCALE ROW (ALGOL TOL THEN NOT NEEDED) .......... DO 120 K = 1, L 120 SCALE = SCALE + ABS(A(I,K)) C IF (SCALE .NE. 0.0E0) GO TO 140 130 E(I) = 0.0E0 E2(I) = 0.0E0 GO TO 290 C 140 DO 150 K = 1, L A(I,K) = A(I,K) / SCALE H = H + A(I,K) * A(I,K) 150 CONTINUE C E2(I) = SCALE * SCALE * H F = A(I,L) G = -SIGN(SQRT(H),F) E(I) = SCALE * G H = H - F * G A(I,L) = F - G IF (L .EQ. 1) GO TO 270 F = 0.0E0 C DO 240 J = 1, L G = 0.0E0 C .......... FORM ELEMENT OF A*U .......... DO 180 K = 1, J 180 G = G + A(J,K) * A(I,K) C JP1 = J + 1 IF (L .LT. JP1) GO TO 220 C DO 200 K = JP1, L 200 G = G + A(K,J) * A(I,K) C .......... FORM ELEMENT OF P .......... 220 E(J) = G / H F = F + E(J) * A(I,J) 240 CONTINUE C H = F / (H + H) C .......... FORM REDUCED A .......... DO 260 J = 1, L F = A(I,J) G = E(J) - H * F E(J) = G C DO 260 K = 1, J A(J,K) = A(J,K) - F * E(K) - G * A(I,K) 260 CONTINUE C 270 DO 280 K = 1, L 280 A(I,K) = SCALE * A(I,K) C 290 H = D(I) D(I) = A(I,I) A(I,I) = H 300 CONTINUE C RETURN END SUBROUTINE TRED2(NM,N,A,D,E,Z) C***BEGIN PROLOGUE TRED2 C***DATE WRITTEN 760101 (YYMMDD) C***REVISION DATE 830518 (YYMMDD) C***CATEGORY NO. D4C1B1 C***KEYWORDS EIGENVALUES,EIGENVECTORS,EISPACK C***AUTHOR SMITH, B. T., ET AL. C***PURPOSE Reduce real symmetric matrix to symmetric tridiagonal C matrix using and accumulating orthogonal transformation C***DESCRIPTION C C This subroutine is a translation of the ALGOL procedure TRED2, C NUM. MATH. 11, 181-195(1968) by Martin, Reinsch, and Wilkinson. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971). C C This subroutine reduces a REAL SYMMETRIC matrix to a C symmetric tridiagonal matrix using and accumulating C orthogonal similarity transformations. C C On Input C C NM must be set to the row dimension of two-dimensional C array parameters as declared in the calling program C dimension statement. C C N is the order of the matrix. C C A contains the real symmetric input matrix. Only the C lower triangle of the matrix need be supplied. C C On Output C C D contains the diagonal elements of the tridiagonal matrix. C C E contains the subdiagonal elements of the tridiagonal C matrix in its last N-1 positions. E(1) is set to zero. C C Z contains the orthogonal transformation matrix C produced in the reduction. C C A and Z may coincide. If distinct, A is unaltered. C C Questions and comments should be directed to B. S. Garbow, C APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY C ------------------------------------------------------------------ C***REFERENCES B. T. SMITH, J. M. BOYLE, J. J. DONGARRA, B. S. GARBOW, C Y. IKEBE, V. C. KLEMA, C. B. MOLER, *MATRIX EIGEN- C SYSTEM ROUTINES - EISPACK GUIDE*, SPRINGER-VERLAG, C 1976. C***ROUTINES CALLED (NONE) C***END PROLOGUE TRED2 C INTEGER I,J,K,L,N,II,NM,JP1 REAL A(NM,N),D(N),E(N),Z(NM,N) REAL F,G,H,HH,SCALE C C***FIRST EXECUTABLE STATEMENT TRED2 DO 100 I = 1, N C DO 100 J = 1, I Z(I,J) = A(I,J) 100 CONTINUE C IF (N .EQ. 1) GO TO 320 C .......... FOR I=N STEP -1 UNTIL 2 DO -- .......... DO 300 II = 2, N I = N + 2 - II L = I - 1 H = 0.0E0 SCALE = 0.0E0 IF (L .LT. 2) GO TO 130 C .......... SCALE ROW (ALGOL TOL THEN NOT NEEDED) .......... DO 120 K = 1, L 120 SCALE = SCALE + ABS(Z(I,K)) C IF (SCALE .NE. 0.0E0) GO TO 140 130 E(I) = Z(I,L) GO TO 290 C 140 DO 150 K = 1, L Z(I,K) = Z(I,K) / SCALE H = H + Z(I,K) * Z(I,K) 150 CONTINUE C F = Z(I,L) G = -SIGN(SQRT(H),F) E(I) = SCALE * G H = H - F * G Z(I,L) = F - G F = 0.0E0 C DO 240 J = 1, L Z(J,I) = Z(I,J) / H G = 0.0E0 C .......... FORM ELEMENT OF A*U .......... DO 180 K = 1, J 180 G = G + Z(J,K) * Z(I,K) C JP1 = J + 1 IF (L .LT. JP1) GO TO 220 C DO 200 K = JP1, L 200 G = G + Z(K,J) * Z(I,K) C .......... FORM ELEMENT OF P .......... 220 E(J) = G / H F = F + E(J) * Z(I,J) 240 CONTINUE C HH = F / (H + H) C .......... FORM REDUCED A .......... DO 260 J = 1, L F = Z(I,J) G = E(J) - HH * F E(J) = G C DO 260 K = 1, J Z(J,K) = Z(J,K) - F * E(K) - G * Z(I,K) 260 CONTINUE C 290 D(I) = H 300 CONTINUE C 320 D(1) = 0.0E0 E(1) = 0.0E0 C .......... ACCUMULATION OF TRANSFORMATION MATRICES .......... DO 500 I = 1, N L = I - 1 IF (D(I) .EQ. 0.0E0) GO TO 380 C DO 360 J = 1, L G = 0.0E0 C DO 340 K = 1, L 340 G = G + Z(I,K) * Z(K,J) C DO 360 K = 1, L Z(K,J) = Z(K,J) - G * Z(K,I) 360 CONTINUE C 380 D(I) = Z(I,I) Z(I,I) = 1.0E0 IF (L .LT. 1) GO TO 500 C DO 400 J = 1, L Z(I,J) = 0.0E0 Z(J,I) = 0.0E0 400 CONTINUE C 500 CONTINUE C RETURN END REAL FUNCTION R1MACH(I) C***BEGIN PROLOGUE R1MACH C***DATE WRITTEN 790101 (YYMMDD) C***REVISION DATE 860825 (YYMMDD) C***CATEGORY NO. R1 C***KEYWORDS MACHINE CONSTANTS C***AUTHOR FOX, P. A., (BELL LABS) C HALL, A. D., (BELL LABS) C SCHRYER, N. L., (BELL LABS) C***PURPOSE Returns single precision machine dependent constants C***DESCRIPTION C C This is the CMLIB version of R1MACH, the real machine C constants subroutine originally developed for the PORT library. C C R1MACH can be used to obtain machine-dependent parameters C for the local machine environment. It is a function C subroutine with one (input) argument, and can be called C as follows, for example C C A = R1MACH(I) C C where I=1,...,5. The (output) value of A above is C determined by the (input) value of I. The results for C various values of I are discussed below. C C Single-Precision Machine Constants C R1MACH(1) = B**(EMIN-1), the smallest positive magnitude. C R1MACH(2) = B**EMAX*(1 - B**(-T)), the largest magnitude. C R1MACH(3) = B**(-T), the smallest relative spacing. C R1MACH(4) = B**(1-T), the largest relative spacing. C R1MACH(5) = LOG10(B) C***REFERENCES FOX, P.A., HALL, A.D., SCHRYER, N.L, *FRAMEWORK FOR C A PORTABLE LIBRARY*, ACM TRANSACTIONS ON MATHE- C MATICAL SOFTWARE, VOL. 4, NO. 2, JUNE 1978, C PP. 177-188. C***ROUTINES CALLED XERROR C***END PROLOGUE R1MACH C INTEGER SMALL(2) INTEGER LARGE(2) INTEGER RIGHT(2) INTEGER DIVER(2) INTEGER LOG10(2) C REAL RMACH(5) C EQUIVALENCE (RMACH(1),SMALL(1)) EQUIVALENCE (RMACH(2),LARGE(1)) EQUIVALENCE (RMACH(3),RIGHT(1)) EQUIVALENCE (RMACH(4),DIVER(1)) EQUIVALENCE (RMACH(5),LOG10(1)) C C C MACHINE CONSTANTS FOR IEEE ARITHMETIC MACHINES, SUCH AS THE AT&T C 3B SERIES, MOTOROLA 68000 BASED MACHINES (E.G. SUN 3 AND AT&T C PC 7300), SUN SPARCSTATIONS, SILICON GRAPHCS WORKSTATIONS, HP C 9000 WORKSTATIONS, IBM RS/6000 WORKSTATIONS, AND 8087 BASED C MICROS (E.G. IBM PC AND AT&T 6300). C C === MACHINE = ATT.3B C === MACHINE = ATT.6300 C === MACHINE = ATT.7300 C === MACHINE = HP.9000 C === MACHINE = IBM.PC C === MACHINE = IBM.RS6000 C === MACHINE = IEEE.MOST-SIG-BYTE-FIRST C === MACHINE = IEEE.LEAST-SIG-BYTE-FIRST C === MACHINE = SGI C === MACHINE = SUN C === MACHINE = 68000 C === MACHINE = 8087 C DATA SMALL(1) / 8388608 / C DATA LARGE(1) / 2139095039 / C DATA RIGHT(1) / 864026624 / C DATA DIVER(1) / 872415232 / C DATA LOG10(1) / 1050288283 / C C MACHINE CONSTANTS FOR SUN WORKSTATIONS. f77 WITH -r8 OPTION. C MACHINE CONSTANTS FOR IBM RS/6000 WORKSTATIONS WITH -qautodbl=dblpad. C C === MACHINE = IBM.RS6000.XLF-WITH-AUTODBL-OPTION C === MACHINE = SUN.F77-WITH-R8-OPTION C === MACHINE = SGI.ORIGIN.F77-WITH-R8-D16-OPTIONS C DATA RMACH(1) / 2.2250738585072E-307 / C DATA RMACH(2) / 1.7976931348623E308 / C DATA RMACH(3) / 1.1102230246252E-16 / C DATA RMACH(4) / 2.2204460492503E-16 / C DATA RMACH(5) / 0.30102999566398 / C C MACHINE CONSTANTS FOR AMDAHL MACHINES. C C === MACHINE = AMDAHL C DATA SMALL(1) / 1048576 / C DATA LARGE(1) / 2147483647 / C DATA RIGHT(1) / 990904320 / C DATA DIVER(1) / 1007681536 / C DATA LOG10(1) / 1091781651 / C C MACHINE CONSTANTS FOR THE BURROUGHS 1700 SYSTEM. C C === MACHINE = BURROUGHS.1700 C DATA RMACH(1) / Z400800000 / C DATA RMACH(2) / Z5FFFFFFFF / C DATA RMACH(3) / Z4E9800000 / C DATA RMACH(4) / Z4EA800000 / C DATA RMACH(5) / Z500E730E8 / C C MACHINE CONSTANTS FOR THE BURROUGHS 5700/6700/7700 SYSTEMS. C C === MACHINE = BURROUGHS.5700 C === MACHINE = BURROUGHS.6700 C === MACHINE = BURROUGHS.7700 C DATA RMACH(1) / O1771000000000000 / C DATA RMACH(2) / O0777777777777777 / C DATA RMACH(3) / O1311000000000000 / C DATA RMACH(4) / O1301000000000000 / C DATA RMACH(5) / O1157163034761675 / C C MACHINE CONSTANTS FOR THE CONVEX C1, C2, C3 SERIES (NATIVE MODE) C C === MACHINE = CONVEX C DATA RMACH(1) / 2.9387360E-39 / C DATA RMACH(2) / 1.7014117E+38 / C DATA RMACH(3) / 5.9604645E-08 / C DATA RMACH(4) / 1.1920929E-07 / C DATA RMACH(5) / 3.0102999E-01 / C C MACHINE CONSTANTS FOR THE CONVEX C1, C2, C3 (NATIVE MODE) C WITH -P8 OPTION C C === MACHINE = CONVEX.P8 C DATA RMACH(1) / 5.562684646268007E-309 / C DATA RMACH(2) / 8.988465674311577E+307 / C DATA RMACH(3) / 1.110223024625157E-016 / C DATA RMACH(4) / 2.220446049250313E-016 / C DATA RMACH(5) / 3.010299956639812E-001 / C C MACHINE CONSTANTS FOR THE CONVEX C1, C2, C3 (IEEE MODE) C C === MACHINE = CONVEX.IEEE C DATA RMACH(1) / 1.1754945E-38 / C DATA RMACH(2) / 3.4028234E+38 / C DATA RMACH(3) / 5.9604645E-08 / C DATA RMACH(4) / 1.1920929E-07 / C DATA RMACH(5) / 3.0102999E-01 / C C MACHINE CONSTANTS FOR THE CONVEX C1, C2, C3 (IEEE MODE) C WITH -P8 OPTION C C === MACHINE = CONVEX.IEEE.P8 C DATA RMACH(1) / 2.225073858507202E-308 / C DATA RMACH(2) / 1.797693134862315E+308 / C DATA RMACH(3) / 1.110223024625157E-016 / C DATA RMACH(4) / 2.220446049250313E-016 / C DATA RMACH(5) / 3.010299956639812E-001 / C C MACHINE CONSTANTS FOR THE CYBER 170/180 SERIES USING NOS (FTN5). C C === MACHINE = CYBER.170.NOS C === MACHINE = CYBER.180.NOS C DATA RMACH(1) / O"00014000000000000000" / C DATA RMACH(2) / O"37767777777777777777" / C DATA RMACH(3) / O"16404000000000000000" / C DATA RMACH(4) / O"16414000000000000000" / C DATA RMACH(5) / O"17164642023241175720" / C C MACHINE CONSTANTS FOR THE CDC 180 SERIES USING NOS/VE C C === MACHINE = CYBER.180.NOS/VE C DATA RMACH(1) / Z"3001800000000000" / C DATA RMACH(2) / Z"4FFEFFFFFFFFFFFE" / C DATA RMACH(3) / Z"3FD2800000000000" / C DATA RMACH(4) / Z"3FD3800000000000" / C DATA RMACH(5) / Z"3FFF9A209A84FBCF" / C C MACHINE CONSTANTS FOR THE CYBER 205 C C === MACHINE = CYBER.205 C DATA RMACH(1) / X'9000400000000000' / C DATA RMACH(2) / X'6FFF7FFFFFFFFFFF' / C DATA RMACH(3) / X'FFA3400000000000' / C DATA RMACH(4) / X'FFA4400000000000' / C DATA RMACH(5) / X'FFD04D104D427DE8' / C C MACHINE CONSTANTS FOR THE CDC 6000/7000 SERIES. C C === MACHINE = CDC.6000 C === MACHINE = CDC.7000 C DATA RMACH(1) / 00014000000000000000B / C DATA RMACH(2) / 37767777777777777777B / C DATA RMACH(3) / 16404000000000000000B / C DATA RMACH(4) / 16414000000000000000B / C DATA RMACH(5) / 17164642023241175720B / C C MACHINE CONSTANTS FOR THE CRAY 1, XMP, 2, AND 3. C C === MACHINE = CRAY.46-BIT-INTEGER C === MACHINE = CRAY.64-BIT-INTEGER C DATA RMACH(1) / 200034000000000000000B / C DATA RMACH(2) / 577767777777777777776B / C DATA RMACH(3) / 377224000000000000000B / C DATA RMACH(4) / 377234000000000000000B / C DATA RMACH(5) / 377774642023241175720B / C C MACHINE CONSTANTS FOR THE DATA GENERAL ECLIPSE S/200 C C NOTE - IT MAY BE APPROPRIATE TO INCLUDE THE FOLLOWING LINE - C STATIC RMACH(5) C C === MACHINE = DATA_GENERAL.ECLIPSE.S/200 C DATA SMALL/20K,0/,LARGE/77777K,177777K/ C DATA RIGHT/35420K,0/,DIVER/36020K,0/ C DATA LOG10/40423K,42023K/ C C ELXSI 6400 C C === MACHINE = ELSXI.6400 C DATA SMALL(1) / '00800000'X / C DATA LARGE(1) / '7F7FFFFF'X / C DATA RIGHT(1) / '33800000'X / C DATA DIVER(1) / '34000000'X / C DATA LOG10(1) / '3E9A209B'X / C C MACHINE CONSTANTS FOR THE HARRIS 220 C MACHINE CONSTANTS FOR THE HARRIS SLASH 6 AND SLASH 7 C C === MACHINE = HARRIS.220 C === MACHINE = HARRIS.SLASH6 C === MACHINE = HARRIS.SLASH7 C DATA SMALL(1),SMALL(2) / '20000000, '00000201 / C DATA LARGE(1),LARGE(2) / '37777777, '00000177 / C DATA RIGHT(1),RIGHT(2) / '20000000, '00000352 / C DATA DIVER(1),DIVER(2) / '20000000, '00000353 / C DATA LOG10(1),LOG10(2) / '23210115, '00000377 / C C MACHINE CONSTANTS FOR THE HONEYWELL 600/6000 SERIES. C MACHINE CONSTANTS FOR THE HONEYWELL DPS 8/70 SERIES. C C === MACHINE = HONEYWELL.600/6000 C === MACHINE = HONEYWELL.DPS.8/70 C DATA RMACH(1) / O402400000000 / C DATA RMACH(2) / O376777777777 / C DATA RMACH(3) / O714400000000 / C DATA RMACH(4) / O716400000000 / C DATA RMACH(5) / O776464202324 / C C MACHINE CONSTANTS FOR THE HP 2100 C 3 WORD DOUBLE PRECISION WITH FTN4 C C === MACHINE = HP.2100.3_WORD_DP C DATA SMALL(1), SMALL(2) / 40000B, 1 / C DATA LARGE(1), LARGE(2) / 77777B, 177776B / C DATA RIGHT(1), RIGHT(2) / 40000B, 325B / C DATA DIVER(1), DIVER(2) / 40000B, 327B / C DATA LOG10(1), LOG10(2) / 46420B, 46777B / C C MACHINE CONSTANTS FOR THE HP 2100 C 4 WORD DOUBLE PRECISION WITH FTN4 C C === MACHINE = HP.2100.4_WORD_DP C DATA SMALL(1), SMALL(2) / 40000B, 1 / C DATA LARGE91), LARGE(2) / 77777B, 177776B / C DATA RIGHT(1), RIGHT(2) / 40000B, 325B / C DATA DIVER(1), DIVER(2) / 40000B, 327B / C DATA LOG10(1), LOG10(2) / 46420B, 46777B / C C MACHINE CONSTANTS FOR THE IBM 360/370 SERIES, C THE XEROX SIGMA 5/7/9, THE SEL SYSTEMS 85/86 AND C THE INTERDATA 3230 AND INTERDATA 7/32. C C === MACHINE = IBM.360 C === MACHINE = IBM.370 C === MACHINE = XEROX.SIGMA.5 C === MACHINE = XEROX.SIGMA.7 C === MACHINE = XEROX.SIGMA.9 C === MACHINE = SEL.85 C === MACHINE = SEL.86 C === MACHINE = INTERDATA.3230 C === MACHINE = INTERDATA.7/32 C DATA RMACH(1) / Z00100000 / C DATA RMACH(2) / Z7FFFFFFF / C DATA RMACH(3) / Z3B100000 / C DATA RMACH(4) / Z3C100000 / C DATA RMACH(5) / Z41134413 / C C MACHINE CONSTANTS FOR THE INTERDATA 8/32 C WITH THE UNIX SYSTEM FORTRAN 77 COMPILER. C C FOR THE INTERDATA FORTRAN VII COMPILER REPLACE C THE Z'S SPECIFYING HEX CONSTANTS WITH Y'S. C C === MACHINE = INTERDATA.8/32.UNIX C DATA RMACH(1) / Z'00100000' / C DATA RMACH(2) / Z'7EFFFFFF' / C DATA RMACH(3) / Z'3B100000' / C DATA RMACH(4) / Z'3C100000' / C DATA RMACH(5) / Z'41134413' / C C MACHINE CONSTANTS FOR THE PDP-10 (KA OR KI PROCESSOR). C C === MACHINE = PDP-10.KA C === MACHINE = PDP-10.KI C DATA RMACH(1) / "000400000000 / C DATA RMACH(2) / "377777777777 / C DATA RMACH(3) / "146400000000 / C DATA RMACH(4) / "147400000000 / C DATA RMACH(5) / "177464202324 / C C MACHINE CONSTANTS FOR PDP-11 FORTRAN SUPPORTING C 32-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL). C C === MACHINE = PDP-11.32-BIT C DATA SMALL(1) / 8388608 / C DATA LARGE(1) / 2147483647 / C DATA RIGHT(1) / 880803840 / C DATA DIVER(1) / 889192448 / C DATA LOG10(1) / 1067065499 / C C DATA RMACH(1) / O00040000000 / C DATA RMACH(2) / O17777777777 / C DATA RMACH(3) / O06440000000 / C DATA RMACH(4) / O06500000000 / C DATA RMACH(5) / O07746420233 / C C MACHINE CONSTANTS FOR PDP-11 FORTRAN SUPPORTING C 16-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL). C C === MACHINE = PDP-11.16-BIT C DATA SMALL(1),SMALL(2) / 128, 0 / C DATA LARGE(1),LARGE(2) / 32767, -1 / C DATA RIGHT(1),RIGHT(2) / 13440, 0 / C DATA DIVER(1),DIVER(2) / 13568, 0 / C DATA LOG10(1),LOG10(2) / 16282, 8347 / C C DATA SMALL(1),SMALL(2) / O000200, O000000 / C DATA LARGE(1),LARGE(2) / O077777, O177777 / C DATA RIGHT(1),RIGHT(2) / O032200, O000000 / C DATA DIVER(1),DIVER(2) / O032400, O000000 / C DATA LOG10(1),LOG10(2) / O037632, O020233 / C C MACHINE CONSTANTS FOR THE SEQUENT BALANCE 8000. C C === MACHINE = SEQUENT.BALANCE.8000 C DATA SMALL(1) / $00800000 / C DATA LARGE(1) / $7F7FFFFF / C DATA RIGHT(1) / $33800000 / C DATA DIVER(1) / $34000000 / C DATA LOG10(1) / $3E9A209B / C C MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES. C C === MACHINE = UNIVAC.1100 C DATA RMACH(1) / O000400000000 / C DATA RMACH(2) / O377777777777 / C DATA RMACH(3) / O146400000000 / C DATA RMACH(4) / O147400000000 / C DATA RMACH(5) / O177464202324 / C C MACHINE CONSTANTS FOR THE VAX 11/780 C (EXPRESSED IN INTEGER AND HEXADECIMAL) C *** THE INTEGER FORMAT SHOULD BE OK FOR UNIX SYSTEMS*** C C === MACHINE = VAX.11/780 C DATA SMALL(1) / 128 / C DATA LARGE(1) / -32769 / C DATA RIGHT(1) / 13440 / C DATA DIVER(1) / 13568 / C DATA LOG10(1) / 547045274 / C C ***THE HEX FORMAT BELOW MAY NOT BE SUITABLE FOR UNIX SYSTEMS*** C C DATA SMALL(1) / Z00000080 / C DATA LARGE(1) / ZFFFF7FFF / C DATA RIGHT(1) / Z00003480 / C DATA DIVER(1) / Z00003500 / C DATA LOG10(1) / Z209B3F9A / C C C***FIRST EXECUTABLE STATEMENT R1MACH C R1MACH = RMACH(I) RETURN C END REAL FUNCTION PYTHAG(A,B) C***BEGIN PROLOGUE PYTHAG C***REFER TO EISDOC C C Finds sqrt(A**2+B**2) without overflow or destructive underflow C***ROUTINES CALLED (NONE) C***END PROLOGUE PYTHAG REAL A,B C REAL P,Q,R,S,T C***FIRST EXECUTABLE STATEMENT PYTHAG P = AMAX1(ABS(A),ABS(B)) Q = AMIN1(ABS(A),ABS(B)) IF (Q .EQ. 0.0E0) GO TO 20 10 CONTINUE R = (Q/P)**2 T = 4.0E0 + R IF (T .EQ. 4.0E0) GO TO 20 S = R/T P = P + 2.0E0*P*S Q = Q*S GO TO 10 20 PYTHAG = P RETURN END SUBROUTINE SPOIR(A,LDA,N,V,ITASK,IND,WORK) C***BEGIN PROLOGUE SPOIR C***DATE WRITTEN 800528 (YYMMDD) C***REVISION DATE 820801 (YYMMDD) C***REVISION HISTORY (YYMMDD) C 000330 Modified array declarations. (JEC) C***CATEGORY NO. D2B1B C***KEYWORDS LINEAR EQUATIONS,POSITIVE DEFINITE,SYMMETRIC C***AUTHOR VOORHEES, E., (LANL) C***PURPOSE SPOIR solves a POSITIVE DEFINITE SYMMETRIC C real NXN system of linear equations. Itera- C tive refinement is used to obtain an error C estimate. C***DESCRIPTION C C Subroutine SPOIR solves a real positive definite symmetric C NxN system of single precision linear equations using LINPACK C subroutines SPOFA and SPOSL. One pass of iterative refine- C ment is used only to obtain an estimate of the accuracy. That C is, if A is an NxN real positive definite symmetric matrix C and if X and B are real N-vectors, then SPOIR solves the C equation C C A*X=B. C C The matrix A is first factored into upper and lower C triangular matrices R and R-TRANSPOSE. These C factors are used to calculate the solution, X. C Then the residual vector is found and used C to calculate an estimate of the relative error, IND. C IND estimates the accuracy of the solution only when the C input matrix and the right hand side are represented C exactly in the computer and does not take into account C any errors in the input data. C C If the equation A*X=B is to be solved for more than one vector C B, the factoring of A does not need to be performed again and C the option to only solve (ITASK .GT. 1) will be faster for C the succeeding solutions. In this case, the contents of A, C LDA, N, and WORK must not have been altered by the user C following factorization (ITASK=1). IND will not be changed C by SPOIR in this case. C C Argument Description *** C A REAL(LDA,N) C the doubly subscripted array with dimension (LDA,N) C which contains the coefficient matrix. Only the C upper triangle, including the diagonal, of the C coefficient matrix need be entered. A is not C altered by the routine. C LDA INTEGER C the leading dimension of the array A. LDA must be great- C er than or equal to N. (Terminal error message IND=-1) C N INTEGER C the order of the matrix A. N must be greater than C or equal to one. (Terminal error message IND=-2) C V REAL(N) C on entry, the singly subscripted array(vector) of di- C mension N which contains the right hand side B of a C system of simultaneous linear equations A*X=B. C on return, V contains the solution vector, X . C ITASK INTEGER C If ITASK = 1, the matrix A is factored and then the C linear equation is solved. C If ITASK .GT. 1, the equation is solved using the existing C factored matrix A (stored in WORK). C If ITASK .LT. 1, then terminal terminal error IND=-3 is C printed. C IND INTEGER C GT. 0 IND is a rough estimate of the number of digits C of accuracy in the solution, X. IND=75 means C that the solution vector X is zero. C LT. 0 See error message corresponding to IND below. C WORK REAL(N*(N+1)) C a singly subscripted array of dimension at least N*(N+1). C C Error Messages Printed *** C C IND=-1 terminal N is greater than LDA. C IND=-2 terminal N is less than one. C IND=-3 terminal ITASK is less than one. C IND=-4 Terminal The matrix A is computationally singular C or is not positive definite. C A solution has not been computed. C IND=-10 warning The solution has no apparent significance. C The solution may be inaccurate or the matrix C A may be poorly scaled. C C Note- The above terminal(*fatal*) error messages are C designed to be handled by XERRWV in which C LEVEL=1 (recoverable) and IFLAG=2 . LEVEL=0 C for warning error messages from XERROR. Unless C the user provides otherwise, an error message C will be printed followed by an abort. C***REFERENCES SUBROUTINE SPOIR WAS DEVELOPED BY GROUP C-3, LOS ALAMOS C SCIENTIFIC LABORATORY, LOS ALAMOS, NM 87545. C THE LINPACK SUBROUTINES USED BY SPOIR ARE DESCRIBED IN C DETAIL IN THE *LINPACK USERS GUIDE* PUBLISHED BY C THE SOCIETY FOR INDUSTRIAL AND APPLIED MATHEMATICS C (SIAM) DATED 1979. C***ROUTINES CALLED DSDOT,R1MACH,SASUM,SCOPY,SPOFA,SPOSL,XERROR,XERRWV C***END PROLOGUE SPOIR C INTEGER LDA,N,ITASK,IND,INFO,J REAL A(LDA,N),V(N),WORK(N,*),SASUM,XNORM,DNORM,R1MACH DOUBLE PRECISION DSDOT C***FIRST EXECUTABLE STATEMENT SPOIR IF (LDA.LT.N) GO TO 101 IF (N.LE.0) GO TO 102 IF (ITASK.LT.1) GO TO 103 IF (ITASK.GT.1) GO TO 20 C MOVE MATRIX A TO WORK DO 10 J=1,N CALL SCOPY(N,A(1,J),1,WORK(1,J),1) 10 CONTINUE C C FACTOR MATRIX A INTO R CALL SPOFA(WORK,N,N,INFO) C C CHECK FOR SINGULAR OR NOT POS.DEF. MATRIX IF (INFO.NE.0) GO TO 104 C C SOLVE AFTER FACTORING 20 CONTINUE C MOVE VECTOR B TO WORK CALL SCOPY(N,V(1),1,WORK(1,N+1),1) CALL SPOSL(WORK,N,N,V) C C FORM NORM OF X0 XNORM=SASUM(N,V(1),1) IF(XNORM.NE.0.0) GO TO 30 IND=75 RETURN 30 CONTINUE C C COMPUTE RESIDUAL DO 40 J=1,N WORK(J,N+1)=-DBLE(WORK(J,N+1)) 1 +DSDOT(J-1,A(1,J),1,V(1),1) 2 +DSDOT(N-J+1,A(J,J),LDA,V(J),1) 40 CONTINUE C C SOLVE A*DELTA=R CALL SPOSL(WORK,N,N,WORK(1,N+1)) C C FORM NORM OF DELTA DNORM=SASUM(N,WORK(1,N+1),1) C C COMPUTE IND (ESTIMATE OF NO. OF SIGNIFICANT DIGITS) IND=-INT(ALOG10(AMAX1(R1MACH(4),DNORM/XNORM))) C C CHECK FOR IND GREATER THAN ZERO IF (IND.GT.0) RETURN IND=-10 CALL XERROR( 'SPOIR ERROR (IND=-10) -- SOLUTION MAY HAVE NO SIGNIF 1ICANCE',58,-10,0) RETURN C C IF LDA.LT.N, IND=-1, TERMINAL XERRWV MESSAGE 101 IND=-1 CALL XERRWV( 'SPOIR ERROR (IND=-1) -- LDA=I1 IS LESS THAN N=I2', 148,-1,1,2,LDA,N,0,0,0) RETURN C C IF N.LT.1, IND=-2, TERMINAL XERRWV MESSAGE 102 IND=-2 CALL XERRWV( 'SPOIR ERROR (IND=-2) -- N=I1 IS LESS THAN 1', 143,-2,1,1,N,0,0,0,0) RETURN C C IF ITASK.LT.1, IND=-3, TERMINAL XERRWV MESSAGE 103 IND=-3 CALL XERRWV( 'SPOIR ERROR (IND=-3) -- ITASK=I1 IS LESS THAN 1', 147,-3,1,1,ITASK,0,0,0,0) RETURN C C IF SINGULAR MATRIX, IND=-4, TERMINAL XERRWV MESSAGE 104 IND=-4 CALL XERRWV( 'SPOIR ERROR (IND=-4) SINGULAR OR NOT POS DEF--NO SO 1LUTION',58,-4,1,0,0,0,0,0,0) RETURN C END SUBROUTINE SPOFA(A,LDA,N,INFO) C***BEGIN PROLOGUE SPOFA C***DATE WRITTEN 780814 (YYMMDD) C***REVISION DATE 820801 (YYMMDD) C***REVISION HISTORY (YYMMDD) C 000330 Modified array declarations. (JEC) C***CATEGORY NO. D2B1B C***KEYWORDS FACTOR,LINEAR ALGEBRA,LINPACK,MATRIX,POSITIVE DEFINITE C***AUTHOR MOLER, C. B., (U. OF NEW MEXICO) C***PURPOSE Factors a real SYMMETRIC POSITIVE DEFINITE matrix. C***DESCRIPTION C C SPOFA factors a real symmetric positive definite matrix. C C SPOFA is usually called by SPOCO, but it can be called C directly with a saving in time if RCOND is not needed. C (Time for SPOCO) = (1 + 18/N)*(Time for SPOFA) . C C On Entry C C A REAL(LDA, N) C the symmetric matrix to be factored. Only the C diagonal and upper triangle are used. C C LDA INTEGER C the leading dimension of the array A . C C N INTEGER C the order of the matrix A . C C On Return C C A an upper triangular matrix R so that A = TRANS(R)*R C where TRANS(R) is the transpose. C The strict lower triangle is unaltered. C If INFO .NE. 0 , the factorization is not complete. C C INFO INTEGER C = 0 for normal return. C = K signals an error condition. The leading minor C of order K is not positive definite. C C LINPACK. This version dated 08/14/78 . C Cleve Moler, University of New Mexico, Argonne National Lab. C C Subroutines and Functions C C BLAS SDOT C Fortran SQRT C***REFERENCES DONGARRA J.J., BUNCH J.R., MOLER C.B., STEWART G.W., C *LINPACK USERS GUIDE*, SIAM, 1979. C***ROUTINES CALLED SDOT C***END PROLOGUE SPOFA INTEGER LDA,N,INFO REAL A(LDA,*) C REAL SDOT,T REAL S INTEGER J,JM1,K C BEGIN BLOCK WITH ...EXITS TO 40 C C C***FIRST EXECUTABLE STATEMENT SPOFA DO 30 J = 1, N INFO = J S = 0.0E0 JM1 = J - 1 IF (JM1 .LT. 1) GO TO 20 DO 10 K = 1, JM1 T = A(K,J) - SDOT(K-1,A(1,K),1,A(1,J),1) T = T/A(K,K) A(K,J) = T S = S + T*T 10 CONTINUE 20 CONTINUE S = A(J,J) - S C ......EXIT IF (S .LE. 0.0E0) GO TO 40 A(J,J) = SQRT(S) 30 CONTINUE INFO = 0 40 CONTINUE RETURN END SUBROUTINE SPOSL(A,LDA,N,B) C***BEGIN PROLOGUE SPOSL C***DATE WRITTEN 780814 (YYMMDD) C***REVISION DATE 820801 (YYMMDD) C***REVISION HISTORY (YYMMDD) C 000330 Modified array declarations. (JEC) C***CATEGORY NO. D2B1B C***KEYWORDS LINEAR ALGEBRA,LINPACK,MATRIX,POSITIVE DEFINITE,SOLVE C***AUTHOR MOLER, C. B., (U. OF NEW MEXICO) C***PURPOSE Solves the real SYMMETRIC POSITIVE DEFINITE system A*X=B C using the factors computed by SPOCO or SPOFA. C***DESCRIPTION C C SPOSL solves the real symmetric positive definite system C A * X = B C using the factors computed by SPOCO or SPOFA. C C On Entry C C A REAL(LDA, N) C the output from SPOCO or SPOFA. C C LDA INTEGER C the leading dimension of the array A . C C N INTEGER C the order of the matrix A . C C B REAL(N) C the right hand side vector. C C On Return C C B the solution vector X . C C Error Condition C C A division by zero will occur if the input factor contains C a zero on the diagonal. Technically, this indicates C singularity, but it is usually caused by improper subroutine C arguments. It will not occur if the subroutines are called C correctly and INFO .EQ. 0 . C C To compute INVERSE(A) * C where C is a matrix C with P columns C CALL SPOCO(A,LDA,N,RCOND,Z,INFO) C IF (RCOND is too small .OR. INFO .NE. 0) GO TO ... C DO 10 J = 1, P C CALL SPOSL(A,LDA,N,C(1,J)) C 10 CONTINUE C C LINPACK. This version dated 08/14/78 . C Cleve Moler, University of New Mexico, Argonne National Lab. C C Subroutines and Functions C C BLAS SAXPY,SDOT C***REFERENCES DONGARRA J.J., BUNCH J.R., MOLER C.B., STEWART G.W., C *LINPACK USERS GUIDE*, SIAM, 1979. C***ROUTINES CALLED SAXPY,SDOT C***END PROLOGUE SPOSL INTEGER LDA,N REAL A(LDA,*),B(*) C REAL SDOT,T INTEGER K,KB C C SOLVE TRANS(R)*Y = B C C***FIRST EXECUTABLE STATEMENT SPOSL DO 10 K = 1, N T = SDOT(K-1,A(1,K),1,B(1),1) B(K) = (B(K) - T)/A(K,K) 10 CONTINUE C C SOLVE R*X = Y C DO 20 KB = 1, N K = N + 1 - KB B(K) = B(K)/A(K,K) T = -B(K) CALL SAXPY(K-1,T,A(1,K),1,B(1),1) 20 CONTINUE RETURN END SUBROUTINE XERROR(MESSG,NMESSG,NERR,LEVEL) C***BEGIN PROLOGUE XERROR C***DATE WRITTEN 790801 (YYMMDD) C***REVISION DATE 820801 (YYMMDD) C***CATEGORY NO. R3C C***KEYWORDS ERROR,XERROR PACKAGE C***AUTHOR JONES, R. E., (SNLA) C***PURPOSE Processes an error (diagnostic) message. C***DESCRIPTION C Abstract C XERROR processes a diagnostic message, in a manner C determined by the value of LEVEL and the current value C of the library error control flag, KONTRL. C (See subroutine XSETF for details.) C C Description of Parameters C --Input-- C MESSG - the Hollerith message to be processed, containing C no more than 72 characters. C NMESSG- the actual number of characters in MESSG. C NERR - the error number associated with this message. C NERR must not be zero. C LEVEL - error category. C =2 means this is an unconditionally fatal error. C =1 means this is a recoverable error. (I.e., it is C non-fatal if XSETF has been appropriately called.) C =0 means this is a warning message only. C =-1 means this is a warning message which is to be C printed at most once, regardless of how many C times this call is executed. C C Examples C CALL XERROR('SMOOTH -- NUM WAS ZERO.',23,1,2) C CALL XERROR('INTEG -- LESS THAN FULL ACCURACY ACHIEVED.', C 43,2,1) C CALL XERROR('ROOTER -- ACTUAL ZERO OF F FOUND BEFORE INTERVAL F C 1ULLY COLLAPSED.',65,3,0) C CALL XERROR('EXP -- UNDERFLOWS BEING SET TO ZERO.',39,1,-1) C C Latest revision --- 19 MAR 1980 C Written by Ron Jones, with SLATEC Common Math Library Subcommittee C***REFERENCES JONES R.E., KAHANER D.K., "XERROR, THE SLATEC ERROR- C HANDLING PACKAGE", SAND82-0800, SANDIA LABORATORIES, C 1982. C***ROUTINES CALLED XERRWV C***END PROLOGUE XERROR CHARACTER*(*) MESSG C***FIRST EXECUTABLE STATEMENT XERROR CALL XERRWV(MESSG,NMESSG,NERR,LEVEL,0,0,0,0,0.,0.) RETURN END SUBROUTINE XERRWV(MESSG,NMESSG,NERR,LEVEL,NI,I1,I2,NR,R1,R2) C***BEGIN PROLOGUE XERRWV C***DATE WRITTEN 800319 (YYMMDD) C***REVISION DATE 820801 (YYMMDD) C***CATEGORY NO. R3C C***KEYWORDS ERROR,XERROR PACKAGE C***AUTHOR JONES, R. E., (SNLA) C***PURPOSE Processes error message allowing 2 integer and two real C values to be included in the message. C***DESCRIPTION C Abstract C XERRWV processes a diagnostic message, in a manner C determined by the value of LEVEL and the current value C of the library error control flag, KONTRL. C (See subroutine XSETF for details.) C In addition, up to two integer values and two real C values may be printed along with the message. C C Description of Parameters C --Input-- C MESSG - the Hollerith message to be processed. C NMESSG- the actual number of characters in MESSG. C NERR - the error number associated with this message. C NERR must not be zero. C LEVEL - error category. C =2 means this is an unconditionally fatal error. C =1 means this is a recoverable error. (I.e., it is C non-fatal if XSETF has been appropriately called.) C =0 means this is a warning message only. C =-1 means this is a warning message which is to be C printed at most once, regardless of how many C times this call is executed. C NI - number of integer values to be printed. (0 to 2) C I1 - first integer value. C I2 - second integer value. C NR - number of real values to be printed. (0 to 2) C R1 - first real value. C R2 - second real value. C C Examples C CALL XERRWV('SMOOTH -- NUM (=I1) WAS ZERO.',29,1,2, C 1 1,NUM,0,0,0.,0.) C CALL XERRWV('QUADXY -- REQUESTED ERROR (R1) LESS THAN MINIMUM ( C 1R2).,54,77,1,0,0,0,2,ERRREQ,ERRMIN) C C Latest revision --- 19 MAR 1980 C Written by Ron Jones, with SLATEC Common Math Library Subcommittee C***REFERENCES JONES R.E., KAHANER D.K., "XERROR, THE SLATEC ERROR- C HANDLING PACKAGE", SAND82-0800, SANDIA LABORATORIES, C 1982. C***ROUTINES CALLED FDUMP,I1MACH,J4SAVE,XERABT,XERCTL,XERPRT,XERSAV, C XGETUA C***END PROLOGUE XERRWV CHARACTER*(*) MESSG CHARACTER*20 LFIRST CHARACTER*37 FORM DIMENSION LUN(5) C GET FLAGS C***FIRST EXECUTABLE STATEMENT XERRWV LKNTRL = J4SAVE(2,0,.FALSE.) MAXMES = J4SAVE(4,0,.FALSE.) C CHECK FOR VALID INPUT IF ((NMESSG.GT.0).AND.(NERR.NE.0).AND. 1 (LEVEL.GE.(-1)).AND.(LEVEL.LE.2)) GO TO 10 IF (LKNTRL.GT.0) CALL XERPRT('FATAL ERROR IN...',17) CALL XERPRT('XERROR -- INVALID INPUT',23) IF (LKNTRL.GT.0) CALL FDUMP IF (LKNTRL.GT.0) CALL XERPRT('JOB ABORT DUE TO FATAL ERROR.', 1 29) IF (LKNTRL.GT.0) CALL XERSAV(' ',0,0,0,KDUMMY) CALL XERABT('XERROR -- INVALID INPUT',23) RETURN 10 CONTINUE C RECORD MESSAGE JUNK = J4SAVE(1,NERR,.TRUE.) CALL XERSAV(MESSG,NMESSG,NERR,LEVEL,KOUNT) C LET USER OVERRIDE LFIRST = MESSG LMESSG = NMESSG LERR = NERR LLEVEL = LEVEL CALL XERCTL(LFIRST,LMESSG,LERR,LLEVEL,LKNTRL) C RESET TO ORIGINAL VALUES LMESSG = NMESSG LERR = NERR LLEVEL = LEVEL LKNTRL = MAX0(-2,MIN0(2,LKNTRL)) MKNTRL = IABS(LKNTRL) C DECIDE WHETHER TO PRINT MESSAGE IF ((LLEVEL.LT.2).AND.(LKNTRL.EQ.0)) GO TO 100 IF (((LLEVEL.EQ.(-1)).AND.(KOUNT.GT.MIN0(1,MAXMES))) 1.OR.((LLEVEL.EQ.0) .AND.(KOUNT.GT.MAXMES)) 2.OR.((LLEVEL.EQ.1) .AND.(KOUNT.GT.MAXMES).AND.(MKNTRL.EQ.1)) 3.OR.((LLEVEL.EQ.2) .AND.(KOUNT.GT.MAX0(1,MAXMES)))) GO TO 100 IF (LKNTRL.LE.0) GO TO 20 CALL XERPRT(' ',1) C INTRODUCTION IF (LLEVEL.EQ.(-1)) CALL XERPRT 1('WARNING MESSAGE...THIS MESSAGE WILL ONLY BE PRINTED ONCE.',57) IF (LLEVEL.EQ.0) CALL XERPRT('WARNING IN...',13) IF (LLEVEL.EQ.1) CALL XERPRT 1 ('RECOVERABLE ERROR IN...',23) IF (LLEVEL.EQ.2) CALL XERPRT('FATAL ERROR IN...',17) 20 CONTINUE C MESSAGE CALL XERPRT(MESSG,LMESSG) CALL XGETUA(LUN,NUNIT) ISIZEI = LOG10(FLOAT(I1MACH(9))) + 1.0 ISIZEF = LOG10(FLOAT(I1MACH(10))**I1MACH(11)) + 1.0 DO 50 KUNIT=1,NUNIT IUNIT = LUN(KUNIT) IF (IUNIT.EQ.0) IUNIT = I1MACH(4) DO 22 I=1,MIN(NI,2) WRITE (FORM,21) I,ISIZEI 21 FORMAT ('(11X,21HIN ABOVE MESSAGE, I',I1,'=,I',I2,') ') IF (I.EQ.1) WRITE (IUNIT,FORM) I1 IF (I.EQ.2) WRITE (IUNIT,FORM) I2 22 CONTINUE DO 24 I=1,MIN(NR,2) WRITE (FORM,23) I,ISIZEF+10,ISIZEF 23 FORMAT ('(11X,21HIN ABOVE MESSAGE, R',I1,'=,E', 1 I2,'.',I2,')') IF (I.EQ.1) WRITE (IUNIT,FORM) R1 IF (I.EQ.2) WRITE (IUNIT,FORM) R2 24 CONTINUE IF (LKNTRL.LE.0) GO TO 40 C ERROR NUMBER WRITE (IUNIT,30) LERR 30 FORMAT (15H ERROR NUMBER =,I10) 40 CONTINUE 50 CONTINUE C TRACE-BACK IF (LKNTRL.GT.0) CALL FDUMP 100 CONTINUE IFATAL = 0 IF ((LLEVEL.EQ.2).OR.((LLEVEL.EQ.1).AND.(MKNTRL.EQ.2))) 1IFATAL = 1 C QUIT HERE IF MESSAGE IS NOT FATAL IF (IFATAL.LE.0) RETURN IF ((LKNTRL.LE.0).OR.(KOUNT.GT.MAX0(1,MAXMES))) GO TO 120 C PRINT REASON FOR ABORT IF (LLEVEL.EQ.1) CALL XERPRT 1 ('JOB ABORT DUE TO UNRECOVERED ERROR.',35) IF (LLEVEL.EQ.2) CALL XERPRT 1 ('JOB ABORT DUE TO FATAL ERROR.',29) C PRINT ERROR SUMMARY CALL XERSAV(' ',-1,0,0,KDUMMY) 120 CONTINUE C ABORT IF ((LLEVEL.EQ.2).AND.(KOUNT.GT.MAX0(1,MAXMES))) LMESSG = 0 CALL XERABT(MESSG,LMESSG) RETURN END SUBROUTINE XERSAV(MESSG,NMESSG,NERR,LEVEL,ICOUNT) C***BEGIN PROLOGUE XERSAV C***DATE WRITTEN 800319 (YYMMDD) C***REVISION DATE 820801 (YYMMDD) C***CATEGORY NO. Z C***KEYWORDS ERROR,XERROR PACKAGE C***AUTHOR JONES, R. E., (SNLA) C***PURPOSE Records that an error occurred. C***DESCRIPTION C Abstract C Record that this error occurred. C C Description of Parameters C --Input-- C MESSG, NMESSG, NERR, LEVEL are as in XERROR, C except that when NMESSG=0 the tables will be C dumped and cleared, and when NMESSG is less than zero the C tables will be dumped and not cleared. C --Output-- C ICOUNT will be the number of times this message has C been seen, or zero if the table has overflowed and C does not contain this message specifically. C When NMESSG=0, ICOUNT will not be altered. C C Written by Ron Jones, with SLATEC Common Math Library Subcommittee C Latest revision --- 19 Mar 1980 C***REFERENCES JONES R.E., KAHANER D.K., "XERROR, THE SLATEC ERROR- C HANDLING PACKAGE", SAND82-0800, SANDIA LABORATORIES, C 1982. C***ROUTINES CALLED I1MACH,S88FMT,XGETUA C***END PROLOGUE XERSAV INTEGER LUN(5) CHARACTER*(*) MESSG CHARACTER*20 MESTAB(10),MES DIMENSION NERTAB(10),LEVTAB(10),KOUNT(10) SAVE MESTAB,NERTAB,LEVTAB,KOUNT,KOUNTX C NEXT TWO DATA STATEMENTS ARE NECESSARY TO PROVIDE A BLANK C ERROR TABLE INITIALLY DATA KOUNT(1),KOUNT(2),KOUNT(3),KOUNT(4),KOUNT(5), 1 KOUNT(6),KOUNT(7),KOUNT(8),KOUNT(9),KOUNT(10) 2 /0,0,0,0,0,0,0,0,0,0/ DATA KOUNTX/0/ C***FIRST EXECUTABLE STATEMENT XERSAV IF (NMESSG.GT.0) GO TO 80 C DUMP THE TABLE IF (KOUNT(1).EQ.0) RETURN C PRINT TO EACH UNIT CALL XGETUA(LUN,NUNIT) DO 60 KUNIT=1,NUNIT IUNIT = LUN(KUNIT) IF (IUNIT.EQ.0) IUNIT = I1MACH(4) C PRINT TABLE HEADER WRITE (IUNIT,10) 10 FORMAT (32H0 ERROR MESSAGE SUMMARY/ 1 51H MESSAGE START NERR LEVEL COUNT) C PRINT BODY OF TABLE DO 20 I=1,10 IF (KOUNT(I).EQ.0) GO TO 30 WRITE (IUNIT,15) MESTAB(I),NERTAB(I),LEVTAB(I),KOUNT(I) 15 FORMAT (1X,A20,3I10) 20 CONTINUE 30 CONTINUE C PRINT NUMBER OF OTHER ERRORS IF (KOUNTX.NE.0) WRITE (IUNIT,40) KOUNTX 40 FORMAT (41H0OTHER ERRORS NOT INDIVIDUALLY TABULATED=,I10) WRITE (IUNIT,50) 50 FORMAT (1X) 60 CONTINUE IF (NMESSG.LT.0) RETURN C CLEAR THE ERROR TABLES DO 70 I=1,10 70 KOUNT(I) = 0 KOUNTX = 0 RETURN 80 CONTINUE C PROCESS A MESSAGE... C SEARCH FOR THIS MESSG, OR ELSE AN EMPTY SLOT FOR THIS MESSG, C OR ELSE DETERMINE THAT THE ERROR TABLE IS FULL. MES = MESSG DO 90 I=1,10 II = I IF (KOUNT(I).EQ.0) GO TO 110 IF (MES.NE.MESTAB(I)) GO TO 90 IF (NERR.NE.NERTAB(I)) GO TO 90 IF (LEVEL.NE.LEVTAB(I)) GO TO 90 GO TO 100 90 CONTINUE C THREE POSSIBLE CASES... C TABLE IS FULL KOUNTX = KOUNTX+1 ICOUNT = 1 RETURN C MESSAGE FOUND IN TABLE 100 KOUNT(II) = KOUNT(II) + 1 ICOUNT = KOUNT(II) RETURN C EMPTY SLOT FOUND FOR NEW MESSAGE 110 MESTAB(II) = MES NERTAB(II) = NERR LEVTAB(II) = LEVEL KOUNT(II) = 1 ICOUNT = 1 RETURN END SUBROUTINE XGETUA(IUNITA,N) C***BEGIN PROLOGUE XGETUA C***DATE WRITTEN 790801 (YYMMDD) C***REVISION DATE 820801 (YYMMDD) C***CATEGORY NO. R3C C***KEYWORDS ERROR,XERROR PACKAGE C***AUTHOR JONES, R. E., (SNLA) C***PURPOSE Returns unit number(s) to which error messages are being C sent. C***DESCRIPTION C Abstract C XGETUA may be called to determine the unit number or numbers C to which error messages are being sent. C These unit numbers may have been set by a call to XSETUN, C or a call to XSETUA, or may be a default value. C C Description of Parameters C --Output-- C IUNIT - an array of one to five unit numbers, depending C on the value of N. A value of zero refers to the C default unit, as defined by the I1MACH machine C constant routine. Only IUNIT(1),...,IUNIT(N) are C defined by XGETUA. The values of IUNIT(N+1),..., C IUNIT(5) are not defined (for N .LT. 5) or altered C in any way by XGETUA. C N - the number of units to which copies of the C error messages are being sent. N will be in the C range from 1 to 5. C C Latest revision --- 19 MAR 1980 C Written by Ron Jones, with SLATEC Common Math Library Subcommittee C***REFERENCES JONES R.E., KAHANER D.K., "XERROR, THE SLATEC ERROR- C HANDLING PACKAGE", SAND82-0800, SANDIA LABORATORIES, C 1982. C***ROUTINES CALLED J4SAVE C***END PROLOGUE XGETUA DIMENSION IUNITA(5) C***FIRST EXECUTABLE STATEMENT XGETUA N = J4SAVE(5,0,.FALSE.) DO 30 I=1,N INDEX = I+4 IF (I.EQ.1) INDEX = 3 IUNITA(I) = J4SAVE(INDEX,0,.FALSE.) 30 CONTINUE RETURN END INTEGER FUNCTION I1MACH(I) C***BEGIN PROLOGUE I1MACH C***DATE WRITTEN 750101 (YYMMDD) C***REVISION DATE 840405 (YYMMDD) C***CATEGORY NO. R1 C***KEYWORDS MACHINE CONSTANTS C***AUTHOR FOX, P. A., (BELL LABS) C HALL, A. D., (BELL LABS) C SCHRYER, N. L., (BELL LABS) C***PURPOSE Returns integer machine dependent constants C***DESCRIPTION C C This is the CMLIB version of I1MACH, the integer machine C constants subroutine originally developed for the PORT library. C C I1MACH can be used to obtain machine-dependent parameters C for the local machine environment. It is a function C subroutine with one (input) argument, and can be called C as follows, for example C C K = I1MACH(I) C C where I=1,...,16. The (output) value of K above is C determined by the (input) value of I. The results for C various values of I are discussed below. C C I/O unit numbers. C I1MACH( 1) = the standard input unit. C I1MACH( 2) = the standard output unit. C I1MACH( 3) = the standard punch unit. C I1MACH( 4) = the standard error message unit. C C Words. C I1MACH( 5) = the number of bits per integer storage unit. C I1MACH( 6) = the number of characters per integer storage unit. C C Integers. C assume integers are represented in the S-digit, base-A form C C sign ( X(S-1)*A**(S-1) + ... + X(1)*A + X(0) ) C C where 0 .LE. X(I) .LT. A for I=0,...,S-1. C I1MACH( 7) = A, the base. C I1MACH( 8) = S, the number of base-A digits. C I1MACH( 9) = A**S - 1, the largest magnitude. C C Floating-Point Numbers. C Assume floating-point numbers are represented in the T-digit, C base-B form C sign (B**E)*( (X(1)/B) + ... + (X(T)/B**T) ) C C where 0 .LE. X(I) .LT. B for I=1,...,T, C 0 .LT. X(1), and EMIN .LE. E .LE. EMAX. C I1MACH(10) = B, the base. C C Single-Precision C I1MACH(11) = T, the number of base-B digits. C I1MACH(12) = EMIN, the smallest exponent E. C I1MACH(13) = EMAX, the largest exponent E. C C Double-Precision C I1MACH(14) = T, the number of base-B digits. C I1MACH(15) = EMIN, the smallest exponent E. C I1MACH(16) = EMAX, the largest exponent E. C C To alter this function for a particular environment, C the desired set of DATA statements should be activated by C removing the C from column 1. Also, the values of C I1MACH(1) - I1MACH(4) should be checked for consistency C with the local operating system. C***REFERENCES FOX P.A., HALL A.D., SCHRYER N.L.,*FRAMEWORK FOR A C PORTABLE LIBRARY*, ACM TRANSACTIONS ON MATHEMATICAL C SOFTWARE, VOL. 4, NO. 2, JUNE 1978, PP. 177-188. C***ROUTINES CALLED (NONE) C***END PROLOGUE I1MACH C INTEGER IMACH(16),OUTPUT EQUIVALENCE (IMACH(4),OUTPUT) C C MACHINE CONSTANTS FOR IEEE ARITHMETIC MACHINES, SUCH AS THE AT&T C 3B SERIES, MOTOROLA 68000 BASED MACHINES (E.G. SUN 3 AND AT&T C PC 7300), SUN SPARCSTATIONS, SILICON GRAPHCS WORKSTATIONS, HP C 9000 WORKSTATIONS, IBM RS/6000 WORKSTATIONS, AND 8087 BASED C MICROS (E.G. IBM PC AND AT&T 6300). C C === MACHINE = ATT.3B C === MACHINE = ATT.6300 C === MACHINE = ATT.7300 C === MACHINE = HP.9000 C === MACHINE = IBM.PC C === MACHINE = IBM.RS6000 C === MACHINE = IEEE.MOST-SIG-BYTE-FIRST C === MACHINE = IEEE.LEAST-SIG-BYTE-FIRST C === MACHINE = SGI C === MACHINE = SUN C === MACHINE = 68000 C === MACHINE = 8087 C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 7 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 32 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 31 / C DATA IMACH( 9) / 2147483647 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 24 / C DATA IMACH(12) / -125 / C DATA IMACH(13) / 128 / C DATA IMACH(14) / 53 / C DATA IMACH(15) / -1021 / C DATA IMACH(16) / 1024 / C C MACHINE CONSTANTS FOR SUN WORKSTATIONS. f77 WITH -r8 OPTION. C C === MACHINE = SUN.F77-WITH-R8-OPTION C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 6 / C DATA IMACH( 4) / 0 / C DATA IMACH( 5) / 32 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 31 / C DATA IMACH( 9) / 2147483647 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 53 / C DATA IMACH(12) / -1021 / C DATA IMACH(13) / 1024 / C DATA IMACH(14) / 113 / C DATA IMACH(15) / -16382 / C DATA IMACH(16) / 16384 / C C MACHINE CONSTANTS FOR SGI Origin 2000 with -r8 -d16 options. C C === MACHINE = SGI.ORIGIN.F77-WITH-R8-D16-OPTIONS C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 6 / C DATA IMACH( 4) / 0 / C DATA IMACH( 5) / 32 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 31 / C DATA IMACH( 9) / 2147483647 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 53 / C DATA IMACH(12) / -1022 / C DATA IMACH(13) / 1024 / C DATA IMACH(14) / 107 / C DATA IMACH(15) / -916 / C DATA IMACH(16) / 1025 / C C MACHINE CONSTANTS FOR IBM RS/6000 WORKSTATIONS WITH -qautodbl=dblpad. C C === MACHINE = IBM.RS6000.XLF-WITH-AUTODBL-OPTION C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 6 / C DATA IMACH( 4) / 0 / C DATA IMACH( 5) / 32 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 31 / C DATA IMACH( 9) / 2147483647 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 53 / C DATA IMACH(12) / -1021 / C DATA IMACH(13) / 1024 / C DATA IMACH(14) / 114 / C DATA IMACH(15) / -1021 / C DATA IMACH(16) / 1024 / C C MACHINE CONSTANTS FOR AMDAHL MACHINES. C C === MACHINE = AMDAHL C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 7 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 32 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 31 / C DATA IMACH( 9) / 2147483647 / C DATA IMACH(10) / 16 / C DATA IMACH(11) / 6 / C DATA IMACH(12) / -64 / C DATA IMACH(13) / 63 / C DATA IMACH(14) / 14 / C DATA IMACH(15) / -64 / C DATA IMACH(16) / 63 / C C MACHINE CONSTANTS FOR THE BURROUGHS 1700 SYSTEM. C C === MACHINE = BURROUGHS.1700 C DATA IMACH( 1) / 7 / C DATA IMACH( 2) / 2 / C DATA IMACH( 3) / 2 / C DATA IMACH( 4) / 2 / C DATA IMACH( 5) / 36 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 33 / C DATA IMACH( 9) / Z1FFFFFFFF / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 24 / C DATA IMACH(12) / -256 / C DATA IMACH(13) / 255 / C DATA IMACH(14) / 60 / C DATA IMACH(15) / -256 / C DATA IMACH(16) / 255 / C C MACHINE CONSTANTS FOR THE BURROUGHS 5700 SYSTEM. C C === MACHINE = BURROUGHS.5700 C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 7 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 48 / C DATA IMACH( 6) / 6 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 39 / C DATA IMACH( 9) / O0007777777777777 / C DATA IMACH(10) / 8 / C DATA IMACH(11) / 13 / C DATA IMACH(12) / -50 / C DATA IMACH(13) / 76 / C DATA IMACH(14) / 26 / C DATA IMACH(15) / -50 / C DATA IMACH(16) / 76 / C C MACHINE CONSTANTS FOR THE BURROUGHS 6700/7700 SYSTEMS. C C === MACHINE = BURROUGHS.6700 C === MACHINE = BURROUGHS.7700 C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 7 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 48 / C DATA IMACH( 6) / 6 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 39 / C DATA IMACH( 9) / O0007777777777777 / C DATA IMACH(10) / 8 / C DATA IMACH(11) / 13 / C DATA IMACH(12) / -50 / C DATA IMACH(13) / 76 / C DATA IMACH(14) / 26 / C DATA IMACH(15) / -32754 / C DATA IMACH(16) / 32780 / C C MACHINE CONSTANTS FOR THE CONVEX C1, C2, C3 SERIES C C === MACHINE = CONVEX C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 0 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 32 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 31 / C DATA IMACH( 9) / 2147483647 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 24 / C DATA IMACH(12) / -127 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 53 / C DATA IMACH(15) / -1023 / C DATA IMACH(16) / 1023 / C C MACHINE CONSTANTS FOR THE CONVEX C1, C2, C3 SERIES (NATIVE MODE) C WITH -P8 OPTION C C === MACHINE = CONVEX.P8 C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 0 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 64 / C DATA IMACH( 6) / 8 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 63 / C DATA IMACH( 9) / 9223372036854775807 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 53 / C DATA IMACH(12) / -1023 / C DATA IMACH(13) / 1023 / C DATA IMACH(14) / 112 / C DATA IMACH(15) / -16383 / C DATA IMACH(16) / 16383 / C C MACHINE CONSTANTS FOR THE CONVEX C1, C2, C3 SERIES (IEEE MODE) C C === MACHINE = CONVEX.IEEE C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 0 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 32 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 31 / C DATA IMACH( 9) / 2147483647 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 24 / C DATA IMACH(12) / -125 / C DATA IMACH(13) / 128 / C DATA IMACH(14) / 53 / C DATA IMACH(15) / -1021 / C DATA IMACH(16) / 1024 / C C MACHINE CONSTANTS FOR THE CONVEX C1, C2, C3 SERIES(IEEE MODE) C WITH -P8 OPTION C C === MACHINE = CONVEX.IEEE.P8 C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 0 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 32 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 31 / C DATA IMACH( 9) / 2147483647 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 53 / C DATA IMACH(12) / -1021 / C DATA IMACH(13) / 1024 / C DATA IMACH(14) / 53 / C DATA IMACH(15) / -1021 / C DATA IMACH(16) / 1024 / C C MACHINE CONSTANTS FOR THE CYBER 170/180 SERIES USING NOS (FTN5). C C === MACHINE = CYBER.170.NOS C === MACHINE = CYBER.180.NOS C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 7 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 60 / C DATA IMACH( 6) / 10 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 48 / C DATA IMACH( 9) / O"00007777777777777777" / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 48 / C DATA IMACH(12) / -974 / C DATA IMACH(13) / 1070 / C DATA IMACH(14) / 96 / C DATA IMACH(15) / -927 / C DATA IMACH(16) / 1070 / C C MACHINE CONSTANTS FOR THE CDC 180 SERIES USING NOS/VE C C === MACHINE = CYBER.180.NOS/VE C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 7 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 64 / C DATA IMACH( 6) / 8 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 63 / C DATA IMACH( 9) / 9223372036854775807 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 47 / C DATA IMACH(12) / -4095 / C DATA IMACH(13) / 4094 / C DATA IMACH(14) / 94 / C DATA IMACH(15) / -4095 / C DATA IMACH(16) / 4094 / C C MACHINE CONSTANTS FOR THE CYBER 205 C C === MACHINE = CYBER.205 C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 7 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 64 / C DATA IMACH( 6) / 8 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 47 / C DATA IMACH( 9) / X'00007FFFFFFFFFFF' / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 47 / C DATA IMACH(12) / -28625 / C DATA IMACH(13) / 28718 / C DATA IMACH(14) / 94 / C DATA IMACH(15) / -28625 / C DATA IMACH(16) / 28718 / C C MACHINE CONSTANTS FOR THE CDC 6000/7000 SERIES. C C === MACHINE = CDC.6000 C === MACHINE = CDC.7000 C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 7 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 60 / C DATA IMACH( 6) / 10 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 48 / C DATA IMACH( 9) / 00007777777777777777B / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 48 / C DATA IMACH(12) / -974 / C DATA IMACH(13) / 1070 / C DATA IMACH(14) / 96 / C DATA IMACH(15) / -927 / C DATA IMACH(16) / 1070 / C C MACHINE CONSTANTS FOR THE CRAY 1, XMP, 2, AND 3. C USING THE 46 BIT INTEGER COMPILER OPTION C C === MACHINE = CRAY.46-BIT-INTEGER C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 102 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 64 / C DATA IMACH( 6) / 8 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 46 / C DATA IMACH( 9) / 777777777777777777777B / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 47 / C DATA IMACH(12) / -8189 / C DATA IMACH(13) / 8190 / C DATA IMACH(14) / 94 / C DATA IMACH(15) / -8099 / C DATA IMACH(16) / 8190 / C C MACHINE CONSTANTS FOR THE CRAY 1, XMP, 2, AND 3. C USING THE 64 BIT INTEGER COMPILER OPTION C C === MACHINE = CRAY.64-BIT-INTEGER C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 102 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 64 / C DATA IMACH( 6) / 8 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 63 / C DATA IMACH( 9) / 777777777777777777777B / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 47 / C DATA IMACH(12) / -8189 / C DATA IMACH(13) / 8190 / C DATA IMACH(14) / 94 / C DATA IMACH(15) / -8099 / C DATA IMACH(16) / 8190 /C C MACHINE CONSTANTS FOR THE DATA GENERAL ECLIPSE S/200 C C === MACHINE = DATA_GENERAL.ECLIPSE.S/200 C DATA IMACH( 1) / 11 / C DATA IMACH( 2) / 12 / C DATA IMACH( 3) / 8 / C DATA IMACH( 4) / 10 / C DATA IMACH( 5) / 16 / C DATA IMACH( 6) / 2 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 15 / C DATA IMACH( 9) /32767 / C DATA IMACH(10) / 16 / C DATA IMACH(11) / 6 / C DATA IMACH(12) / -64 / C DATA IMACH(13) / 63 / C DATA IMACH(14) / 14 / C DATA IMACH(15) / -64 / C DATA IMACH(16) / 63 / C C ELXSI 6400 C C === MACHINE = ELSXI.6400 C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 6 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 32 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 32 / C DATA IMACH( 9) / 2147483647 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 24 / C DATA IMACH(12) / -126 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 53 / C DATA IMACH(15) / -1022 / C DATA IMACH(16) / 1023 / C C MACHINE CONSTANTS FOR THE HARRIS 220 C MACHINE CONSTANTS FOR THE HARRIS SLASH 6 AND SLASH 7 C C === MACHINE = HARRIS.220 C === MACHINE = HARRIS.SLASH6 C === MACHINE = HARRIS.SLASH7 C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 0 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 24 / C DATA IMACH( 6) / 3 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 23 / C DATA IMACH( 9) / 8388607 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 23 / C DATA IMACH(12) / -127 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 38 / C DATA IMACH(15) / -127 / C DATA IMACH(16) / 127 / C C MACHINE CONSTANTS FOR THE HONEYWELL 600/6000 SERIES. C MACHINE CONSTANTS FOR THE HONEYWELL DPS 8/70 SERIES. C C === MACHINE = HONEYWELL.600/6000 C === MACHINE = HONEYWELL.DPS.8/70 C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 43 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 36 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 35 / C DATA IMACH( 9) / O377777777777 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 27 / C DATA IMACH(12) / -127 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 63 / C DATA IMACH(15) / -127 / C DATA IMACH(16) / 127 / C C MACHINE CONSTANTS FOR THE HP 2100 C 3 WORD DOUBLE PRECISION OPTION WITH FTN4 C C === MACHINE = HP.2100.3_WORD_DP C DATA IMACH(1) / 5/ C DATA IMACH(2) / 6 / C DATA IMACH(3) / 4 / C DATA IMACH(4) / 1 / C DATA IMACH(5) / 16 / C DATA IMACH(6) / 2 / C DATA IMACH(7) / 2 / C DATA IMACH(8) / 15 / C DATA IMACH(9) / 32767 / C DATA IMACH(10)/ 2 / C DATA IMACH(11)/ 23 / C DATA IMACH(12)/ -128 / C DATA IMACH(13)/ 127 / C DATA IMACH(14)/ 39 / C DATA IMACH(15)/ -128 / C DATA IMACH(16)/ 127 / C C MACHINE CONSTANTS FOR THE HP 2100 C 4 WORD DOUBLE PRECISION OPTION WITH FTN4 C C === MACHINE = HP.2100.4_WORD_DP C DATA IMACH(1) / 5 / C DATA IMACH(2) / 6 / C DATA IMACH(3) / 4 / C DATA IMACH(4) / 1 / C DATA IMACH(5) / 16 / C DATA IMACH(6) / 2 / C DATA IMACH(7) / 2 / C DATA IMACH(8) / 15 / C DATA IMACH(9) / 32767 / C DATA IMACH(10)/ 2 / C DATA IMACH(11)/ 23 / C DATA IMACH(12)/ -128 / C DATA IMACH(13)/ 127 / C DATA IMACH(14)/ 55 / C DATA IMACH(15)/ -128 / C DATA IMACH(16)/ 127 / C C MACHINE CONSTANTS FOR THE IBM 360/370 SERIES, C THE XEROX SIGMA 5/7/9 AND THE SEL SYSTEMS 85/86 AND C THE INTERDATA 3230 AND INTERDATA 7/32. C C === MACHINE = IBM.360 C === MACHINE = IBM.370 C === MACHINE = XEROX.SIGMA.5 C === MACHINE = XEROX.SIGMA.7 C === MACHINE = XEROX.SIGMA.9 C === MACHINE = SEL.85 C === MACHINE = SEL.86 C === MACHINE = INTERDATA.3230 C === MACHINE = INTERDATA.7/32 C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 7 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 32 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 31 / C DATA IMACH( 9) / Z7FFFFFFF / C DATA IMACH(10) / 16 / C DATA IMACH(11) / 6 / C DATA IMACH(12) / -64 / C DATA IMACH(13) / 63 / C DATA IMACH(14) / 14 / C DATA IMACH(15) / -64 / C DATA IMACH(16) / 63 / C C MACHINE CONSTANTS FOR THE INTERDATA 8/32 C WITH THE UNIX SYSTEM FORTRAN 77 COMPILER. C C FOR THE INTERDATA FORTRAN VII COMPILER REPLACE C THE Z'S SPECIFYING HEX CONSTANTS WITH Y'S. C C === MACHINE = INTERDATA.8/32.UNIX C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 6 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 32 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 31 / C DATA IMACH( 9) / Z'7FFFFFFF' / C DATA IMACH(10) / 16 / C DATA IMACH(11) / 6 / C DATA IMACH(12) / -64 / C DATA IMACH(13) / 62 / C DATA IMACH(14) / 14 / C DATA IMACH(15) / -64 / C DATA IMACH(16) / 62 / C C MACHINE CONSTANTS FOR THE PDP-10 (KA PROCESSOR). C C === MACHINE = PDP-10.KA C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 7 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 36 / C DATA IMACH( 6) / 5 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 35 / C DATA IMACH( 9) / "377777777777 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 27 / C DATA IMACH(12) / -128 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 54 / C DATA IMACH(15) / -101 / C DATA IMACH(16) / 127 / C C MACHINE CONSTANTS FOR THE PDP-10 (KI PROCESSOR). C C === MACHINE = PDP-10.KI C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 7 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 36 / C DATA IMACH( 6) / 5 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 35 / C DATA IMACH( 9) / "377777777777 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 27 / C DATA IMACH(12) / -128 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 62 / C DATA IMACH(15) / -128 / C DATA IMACH(16) / 127 / C C MACHINE CONSTANTS FOR PDP-11 FORTRANS SUPPORTING C 32-BIT INTEGER ARITHMETIC. C C === MACHINE = PDP-11.32-BIT C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 7 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 32 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 31 / C DATA IMACH( 9) / 2147483647 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 24 / C DATA IMACH(12) / -127 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 56 / C DATA IMACH(15) / -127 / C DATA IMACH(16) / 127 / C C MACHINE CONSTANTS FOR PDP-11 FORTRANS SUPPORTING C 16-BIT INTEGER ARITHMETIC. C C === MACHINE = PDP-11.16-BIT C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 7 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 16 / C DATA IMACH( 6) / 2 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 15 / C DATA IMACH( 9) / 32767 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 24 / C DATA IMACH(12) / -127 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 56 / C DATA IMACH(15) / -127 / C DATA IMACH(16) / 127 / C C MACHINE CONSTANTS FOR THE SEQUENT BALANCE 8000. C C === MACHINE = SEQUENT.BALANCE.8000 C DATA IMACH( 1) / 0 / C DATA IMACH( 2) / 0 / C DATA IMACH( 3) / 7 / C DATA IMACH( 4) / 0 / C DATA IMACH( 5) / 32 / C DATA IMACH( 6) / 1 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 31 / C DATA IMACH( 9) / 2147483647 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 24 / C DATA IMACH(12) / -125 / C DATA IMACH(13) / 128 / C DATA IMACH(14) / 53 / C DATA IMACH(15) / -1021 / C DATA IMACH(16) / 1024 / C C MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES. FTN COMPILER C C === MACHINE = UNIVAC.1100 C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 1 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 36 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 35 / C DATA IMACH( 9) / O377777777777 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 27 / C DATA IMACH(12) / -128 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 60 / C DATA IMACH(15) /-1024 / C DATA IMACH(16) / 1023 / C C MACHINE CONSTANTS FOR THE VAX 11/780 C C === MACHINE = VAX.11/780 C DATA IMACH(1) / 5 / C DATA IMACH(2) / 6 / C DATA IMACH(3) / 5 / C DATA IMACH(4) / 6 / C DATA IMACH(5) / 32 / C DATA IMACH(6) / 4 / C DATA IMACH(7) / 2 / C DATA IMACH(8) / 31 / C DATA IMACH(9) /2147483647 / C DATA IMACH(10)/ 2 / C DATA IMACH(11)/ 24 / C DATA IMACH(12)/ -127 / C DATA IMACH(13)/ 127 / C DATA IMACH(14)/ 56 / C DATA IMACH(15)/ -127 / C DATA IMACH(16)/ 127 / C C C***FIRST EXECUTABLE STATEMENT I1MACH C I1MACH=IMACH(I) RETURN C END DOUBLE PRECISION FUNCTION DSDOT(N,SX,INCX,SY,INCY) C***BEGIN PROLOGUE DSDOT C***DATE WRITTEN 791001 (YYMMDD) C***REVISION DATE 820801 (YYMMDD) C***REVISION HISTORY (YYMMDD) C 000330 Modified array declarations. (JEC) C C***CATEGORY NO. D1A4 C***KEYWORDS BLAS,DOUBLE PRECISION,INNER PRODUCT,LINEAR ALGEBRA,VECTOR C***AUTHOR LAWSON, C. L., (JPL) C HANSON, R. J., (SNLA) C KINCAID, D. R., (U. OF TEXAS) C KROGH, F. T., (JPL) C***PURPOSE D.P inner product of s.p. vectors C***DESCRIPTION C C B L A S Subprogram C Description of Parameters C C --Input-- C N number of elements in input vector(s) C SX single precision vector with N elements C INCX storage spacing between elements of SX C SY single precision vector with N elements C INCY storage spacing between elements of SY C C --Output-- C DSDOT double precision dot product (zero if N.LE.0) C C Returns D.P. dot product accumulated in D.P., for S.P. SX and SY C DSDOT = sum for I = 0 to N-1 of SX(LX+I*INCX) * SY(LY+I*INCY), C where LX = 1 if INCX .GE. 0, else LX = (-INCX)*N, and LY is C defined in a similar way using INCY. C***REFERENCES LAWSON C.L., HANSON R.J., KINCAID D.R., KROGH F.T., C *BASIC LINEAR ALGEBRA SUBPROGRAMS FOR FORTRAN USAGE*, C ALGORITHM NO. 539, TRANSACTIONS ON MATHEMATICAL C SOFTWARE, VOLUME 5, NUMBER 3, SEPTEMBER 1979, 308-323 C***ROUTINES CALLED (NONE) C***END PROLOGUE DSDOT C REAL SX(*),SY(*) C***FIRST EXECUTABLE STATEMENT DSDOT DSDOT = 0.D0 IF(N .LE. 0)RETURN IF(INCX.EQ.INCY.AND.INCX.GT.0) GO TO 20 KX = 1 KY = 1 IF(INCX.LT.0) KX = 1+(1-N)*INCX IF(INCY.LT.0) KY = 1+(1-N)*INCY DO 10 I = 1,N DSDOT = DSDOT + DBLE(SX(KX))*DBLE(SY(KY)) KX = KX + INCX KY = KY + INCY 10 CONTINUE RETURN 20 CONTINUE NS = N*INCX DO 30 I=1,NS,INCX DSDOT = DSDOT + DBLE(SX(I))*DBLE(SY(I)) 30 CONTINUE RETURN END REAL FUNCTION SASUM(N,SX,INCX) C***BEGIN PROLOGUE SASUM C***DATE WRITTEN 791001 (YYMMDD) C***REVISION DATE 820801 (YYMMDD) C***REVISION HISTORY (YYMMDD) C 000330 Modified array declarations. (JEC) C C***CATEGORY NO. D1A3A C***KEYWORDS ADD,BLAS,LINEAR ALGEBRA,MAGNITUDE,SUM,VECTOR C***AUTHOR LAWSON, C. L., (JPL) C HANSON, R. J., (SNLA) C KINCAID, D. R., (U. OF TEXAS) C KROGH, F. T., (JPL) C***PURPOSE Sum of magnitudes of s.p vector components C***DESCRIPTION C C B L A S Subprogram C Description of Parameters C C --Input-- C N number of elements in input vector(S) C SX single precision vector with N elements C INCX storage spacing between elements of SX C C --Output-- C SASUM single precision result (zero if N .LE. 0) C C Returns sum of magnitudes of single precision SX. C SASUM = sum from 0 to N-1 of ABS(SX(1+I*INCX)) C***REFERENCES LAWSON C.L., HANSON R.J., KINCAID D.R., KROGH F.T., C *BASIC LINEAR ALGEBRA SUBPROGRAMS FOR FORTRAN USAGE*, C ALGORITHM NO. 539, TRANSACTIONS ON MATHEMATICAL C SOFTWARE, VOLUME 5, NUMBER 3, SEPTEMBER 1979, 308-323 C***ROUTINES CALLED (NONE) C***END PROLOGUE SASUM C REAL SX(*) C***FIRST EXECUTABLE STATEMENT SASUM SASUM = 0.0E0 IF(N.LE.0)RETURN IF(INCX.EQ.1)GOTO 20 C C CODE FOR INCREMENTS NOT EQUAL TO 1. C NS = N*INCX DO 10 I=1,NS,INCX SASUM = SASUM + ABS(SX(I)) 10 CONTINUE RETURN C C CODE FOR INCREMENTS EQUAL TO 1. C C C CLEAN-UP LOOP SO REMAINING VECTOR LENGTH IS A MULTIPLE OF 6. C 20 M = MOD(N,6) IF( M .EQ. 0 ) GO TO 40 DO 30 I = 1,M SASUM = SASUM + ABS(SX(I)) 30 CONTINUE IF( N .LT. 6 ) RETURN 40 MP1 = M + 1 DO 50 I = MP1,N,6 SASUM = SASUM + ABS(SX(I)) + ABS(SX(I + 1)) + ABS(SX(I + 2)) 1 + ABS(SX(I + 3)) + ABS(SX(I + 4)) + ABS(SX(I + 5)) 50 CONTINUE RETURN END SUBROUTINE SAXPY(N,SA,SX,INCX,SY,INCY) C***BEGIN PROLOGUE SAXPY C***DATE WRITTEN 791001 (YYMMDD) C***REVISION DATE 820801 (YYMMDD) C***REVISION HISTORY (YYMMDD) C 000330 Modified array declarations. (JEC) C C***CATEGORY NO. D1A7 C***KEYWORDS BLAS,LINEAR ALGEBRA,TRIAD,VECTOR C***AUTHOR LAWSON, C. L., (JPL) C HANSON, R. J., (SNLA) C KINCAID, D. R., (U. OF TEXAS) C KROGH, F. T., (JPL) C***PURPOSE S.P. computation y = a*x + y C***DESCRIPTION C C B L A S Subprogram C Description of Parameters C C --Input-- C N number of elements in input vector(s) C SA single precision scalar multiplier C SX single precision vector with N elements C INCX storage spacing between elements of SX C SY single precision vector with N elements C INCY storage spacing between elements of SY C C --Output-- C SY single precision result (unchanged if N .LE. 0) C C Overwrite single precision SY with single precision SA*SX +SY. C For I = 0 to N-1, replace SY(LY+I*INCY) with SA*SX(LX+I*INCX) + C SY(LY+I*INCY), where LX = 1 if INCX .GE. 0, else LX = (-INCX)*N C and LY is defined in a similar way using INCY. C***REFERENCES LAWSON C.L., HANSON R.J., KINCAID D.R., KROGH F.T., C *BASIC LINEAR ALGEBRA SUBPROGRAMS FOR FORTRAN USAGE*, C ALGORITHM NO. 539, TRANSACTIONS ON MATHEMATICAL C SOFTWARE, VOLUME 5, NUMBER 3, SEPTEMBER 1979, 308-323 C***ROUTINES CALLED (NONE) C***END PROLOGUE SAXPY C REAL SX(*),SY(*),SA C***FIRST EXECUTABLE STATEMENT SAXPY IF(N.LE.0.OR.SA.EQ.0.E0) RETURN IF(INCX.EQ.INCY) IF(INCX-1) 5,20,60 5 CONTINUE C C CODE FOR NONEQUAL OR NONPOSITIVE INCREMENTS. 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 SY(IY) = SY(IY) + SA*SX(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 SO REMAINING VECTOR LENGTH IS A MULTIPLE OF 4. C 20 M = MOD(N,4) IF( M .EQ. 0 ) GO TO 40 DO 30 I = 1,M SY(I) = SY(I) + SA*SX(I) 30 CONTINUE IF( N .LT. 4 ) RETURN 40 MP1 = M + 1 DO 50 I = MP1,N,4 SY(I) = SY(I) + SA*SX(I) SY(I + 1) = SY(I + 1) + SA*SX(I + 1) SY(I + 2) = SY(I + 2) + SA*SX(I + 2) SY(I + 3) = SY(I + 3) + SA*SX(I + 3) 50 CONTINUE RETURN C C CODE FOR EQUAL, POSITIVE, NONUNIT INCREMENTS. C 60 CONTINUE NS = N*INCX DO 70 I=1,NS,INCX SY(I) = SA*SX(I) + SY(I) 70 CONTINUE RETURN END SUBROUTINE SCOPY(N,SX,INCX,SY,INCY) C***BEGIN PROLOGUE SCOPY C***DATE WRITTEN 791001 (YYMMDD) C***REVISION DATE 820801 (YYMMDD) C***REVISION HISTORY (YYMMDD) C 000330 Modified array declarations. (JEC) C C***CATEGORY NO. D1A5 C***KEYWORDS BLAS,COPY,LINEAR ALGEBRA,VECTOR C***AUTHOR LAWSON, C. L., (JPL) C HANSON, R. J., (SNLA) C KINCAID, D. R., (U. OF TEXAS) C KROGH, F. T., (JPL) C***PURPOSE Copy s.p. vector y = x C***DESCRIPTION C C B L A S Subprogram C Description of Parameters C C --Input-- C N number of elements in input vector(s) C SX single precision vector with N elements C INCX storage spacing between elements of SX C SY single precision vector with N elements C INCY storage spacing between elements of SY C C --Output-- C SY copy of vector SX (unchanged if N .LE. 0) C C Copy single precision SX to single precision SY. C For I = 0 to N-1, copy SX(LX+I*INCX) to SY(LY+I*INCY), C where LX = 1 if INCX .GE. 0, else LX = (-INCX)*N, and LY is C defined in a similar way using INCY. C***REFERENCES LAWSON C.L., HANSON R.J., KINCAID D.R., KROGH F.T., C *BASIC LINEAR ALGEBRA SUBPROGRAMS FOR FORTRAN USAGE*, C ALGORITHM NO. 539, TRANSACTIONS ON MATHEMATICAL C SOFTWARE, VOLUME 5, NUMBER 3, SEPTEMBER 1979, 308-323 C***ROUTINES CALLED (NONE) C***END PROLOGUE SCOPY C REAL SX(*),SY(*) C***FIRST EXECUTABLE STATEMENT SCOPY IF(N.LE.0)RETURN IF(INCX.EQ.INCY) IF(INCX-1) 5,20,60 5 CONTINUE C C CODE FOR UNEQUAL OR NONPOSITIVE INCREMENTS. 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 SY(IY) = SX(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 SO REMAINING VECTOR LENGTH IS A MULTIPLE OF 7. C 20 M = MOD(N,7) IF( M .EQ. 0 ) GO TO 40 DO 30 I = 1,M SY(I) = SX(I) 30 CONTINUE IF( N .LT. 7 ) RETURN 40 MP1 = M + 1 DO 50 I = MP1,N,7 SY(I) = SX(I) SY(I + 1) = SX(I + 1) SY(I + 2) = SX(I + 2) SY(I + 3) = SX(I + 3) SY(I + 4) = SX(I + 4) SY(I + 5) = SX(I + 5) SY(I + 6) = SX(I + 6) 50 CONTINUE RETURN C C CODE FOR EQUAL, POSITIVE, NONUNIT INCREMENTS. C 60 CONTINUE NS = N*INCX DO 70 I=1,NS,INCX SY(I) = SX(I) 70 CONTINUE RETURN END REAL FUNCTION SDOT(N,SX,INCX,SY,INCY) C***BEGIN PROLOGUE SDOT C***DATE WRITTEN 791001 (YYMMDD) C***REVISION DATE 820801 (YYMMDD) C***REVISION HISTORY (YYMMDD) C 000330 Modified array declarations. (JEC) C C***CATEGORY NO. D1A4 C***KEYWORDS BLAS,INNER PRODUCT,LINEAR ALGEBRA,VECTOR C***AUTHOR LAWSON, C. L., (JPL) C HANSON, R. J., (SNLA) C KINCAID, D. R., (U. OF TEXAS) C KROGH, F. T., (JPL) C***PURPOSE S.P. inner product of s.p. vectors C***DESCRIPTION C C B L A S Subprogram C Description of Parameters C C --Input-- C N number of elements in input vector(s) C SX single precision vector with N elements C INCX storage spacing between elements of SX C SY single precision vector with N elements C INCY storage spacing between elements of SY C C --Output-- C SDOT single precision dot product (zero if N .LE. 0) C C Returns the dot product of single precision SX and SY. C SDOT = sum for I = 0 to N-1 of SX(LX+I*INCX) * SY(LY+I*INCY), C where LX = 1 if INCX .GE. 0, else LX = (-INCX)*N, and LY is C defined in a similar way using INCY. C***REFERENCES LAWSON C.L., HANSON R.J., KINCAID D.R., KROGH F.T., C *BASIC LINEAR ALGEBRA SUBPROGRAMS FOR FORTRAN USAGE*, C ALGORITHM NO. 539, TRANSACTIONS ON MATHEMATICAL C SOFTWARE, VOLUME 5, NUMBER 3, SEPTEMBER 1979, 308-323 C***ROUTINES CALLED (NONE) C***END PROLOGUE SDOT C REAL SX(*),SY(*) C***FIRST EXECUTABLE STATEMENT SDOT SDOT = 0.0E0 IF(N.LE.0)RETURN IF(INCX.EQ.INCY) IF(INCX-1)5,20,60 5 CONTINUE C C CODE FOR UNEQUAL INCREMENTS OR NONPOSITIVE INCREMENTS. 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 SDOT = SDOT + SX(IX)*SY(IY) IX = IX + INCX IY = IY + INCY 10 CONTINUE RETURN C C CODE FOR BOTH INCREMENTS EQUAL TO 1 C C C CLEAN-UP LOOP SO REMAINING VECTOR LENGTH IS A MULTIPLE OF 5. C 20 M = MOD(N,5) IF( M .EQ. 0 ) GO TO 40 DO 30 I = 1,M SDOT = SDOT + SX(I)*SY(I) 30 CONTINUE IF( N .LT. 5 ) RETURN 40 MP1 = M + 1 DO 50 I = MP1,N,5 SDOT = SDOT + SX(I)*SY(I) + SX(I + 1)*SY(I + 1) + 1 SX(I + 2)*SY(I + 2) + SX(I + 3)*SY(I + 3) + SX(I + 4)*SY(I + 4) 50 CONTINUE RETURN C C CODE FOR POSITIVE EQUAL INCREMENTS .NE.1. C 60 CONTINUE NS=N*INCX DO 70 I=1,NS,INCX SDOT = SDOT + SX(I)*SY(I) 70 CONTINUE RETURN END SUBROUTINE FDUMP C***BEGIN PROLOGUE FDUMP C***DATE WRITTEN 790801 (YYMMDD) C***REVISION DATE 820801 (YYMMDD) C***CATEGORY NO. Z C***KEYWORDS ERROR,XERROR PACKAGE C***AUTHOR JONES, R. E., (SNLA) C***PURPOSE Symbolic dump (should be locally written). C***DESCRIPTION C ***Note*** Machine Dependent Routine C FDUMP is intended to be replaced by a locally written C version which produces a symbolic dump. Failing this, C it should be replaced by a version which prints the C subprogram nesting list. Note that this dump must be C printed on each of up to five files, as indicated by the C XGETUA routine. See XSETUA and XGETUA for details. C C Written by Ron Jones, with SLATEC Common Math Library Subcommittee C Latest revision --- 23 May 1979 C***ROUTINES CALLED (NONE) C***END PROLOGUE FDUMP C***FIRST EXECUTABLE STATEMENT FDUMP RETURN END FUNCTION J4SAVE(IWHICH,IVALUE,ISET) C***BEGIN PROLOGUE J4SAVE C***REFER TO XERROR C Abstract C J4SAVE saves and recalls several global variables needed C by the library error handling routines. C C Description of Parameters C --Input-- C IWHICH - Index of item desired. C = 1 Refers to current error number. C = 2 Refers to current error control flag. C = 3 Refers to current unit number to which error C messages are to be sent. (0 means use standard.) C = 4 Refers to the maximum number of times any C message is to be printed (as set by XERMAX). C = 5 Refers to the total number of units to which C each error message is to be written. C = 6 Refers to the 2nd unit for error messages C = 7 Refers to the 3rd unit for error messages C = 8 Refers to the 4th unit for error messages C = 9 Refers to the 5th unit for error messages C IVALUE - The value to be set for the IWHICH-th parameter, C if ISET is .TRUE. . C ISET - If ISET=.TRUE., the IWHICH-th parameter will BE C given the value, IVALUE. If ISET=.FALSE., the C IWHICH-th parameter will be unchanged, and IVALUE C is a dummy parameter. C --Output-- C The (old) value of the IWHICH-th parameter will be returned C in the function value, J4SAVE. C C Written by Ron Jones, with SLATEC Common Math Library Subcommittee C Adapted from Bell Laboratories PORT Library Error Handler C Latest revision --- 23 MAY 1979 C***REFERENCES JONES R.E., KAHANER D.K., "XERROR, THE SLATEC ERROR- C HANDLING PACKAGE", SAND82-0800, SANDIA LABORATORIES, C 1982. C***ROUTINES CALLED (NONE) C***END PROLOGUE J4SAVE LOGICAL ISET INTEGER IPARAM(9) SAVE IPARAM DATA IPARAM(1),IPARAM(2),IPARAM(3),IPARAM(4)/0,2,0,10/ DATA IPARAM(5)/1/ DATA IPARAM(6),IPARAM(7),IPARAM(8),IPARAM(9)/0,0,0,0/ C***FIRST EXECUTABLE STATEMENT J4SAVE J4SAVE = IPARAM(IWHICH) IF (ISET) IPARAM(IWHICH) = IVALUE RETURN END SUBROUTINE XERABT(MESSG,NMESSG) C***BEGIN PROLOGUE XERABT C***DATE WRITTEN 790801 (YYMMDD) C***REVISION DATE 820801 (YYMMDD) C***CATEGORY NO. R3C C***KEYWORDS ERROR,XERROR PACKAGE C***AUTHOR JONES, R. E., (SNLA) C***PURPOSE Aborts program execution and prints error message. C***DESCRIPTION C Abstract C ***Note*** machine dependent routine C XERABT aborts the execution of the program. C The error message causing the abort is given in the calling C sequence, in case one needs it for printing on a dayfile, C for example. C C Description of Parameters C MESSG and NMESSG are as in XERROR, except that NMESSG may C be zero, in which case no message is being supplied. C C Written by Ron Jones, with SLATEC Common Math Library Subcommittee C Latest revision --- 19 MAR 1980 C***REFERENCES JONES R.E., KAHANER D.K., "XERROR, THE SLATEC ERROR- C HANDLING PACKAGE", SAND82-0800, SANDIA LABORATORIES, C 1982. C***ROUTINES CALLED (NONE) C***END PROLOGUE XERABT CHARACTER*(*) MESSG C***FIRST EXECUTABLE STATEMENT XERABT STOP END SUBROUTINE XERCTL(MESSG1,NMESSG,NERR,LEVEL,KONTRL) C***BEGIN PROLOGUE XERCTL C***DATE WRITTEN 790801 (YYMMDD) C***REVISION DATE 820801 (YYMMDD) C***CATEGORY NO. R3C C***KEYWORDS ERROR,XERROR PACKAGE C***AUTHOR JONES, R. E., (SNLA) C***PURPOSE Allows user control over handling of individual errors. C***DESCRIPTION C Abstract C Allows user control over handling of individual errors. C Just after each message is recorded, but before it is C processed any further (i.e., before it is printed or C a decision to abort is made), a call is made to XERCTL. C If the user has provided his own version of XERCTL, he C can then override the value of KONTROL used in processing C this message by redefining its value. C KONTRL may be set to any value from -2 to 2. C The meanings for KONTRL are the same as in XSETF, except C that the value of KONTRL changes only for this message. C If KONTRL is set to a value outside the range from -2 to 2, C it will be moved back into that range. C C Description of Parameters C C --Input-- C MESSG1 - the first word (only) of the error message. C NMESSG - same as in the call to XERROR or XERRWV. C NERR - same as in the call to XERROR or XERRWV. C LEVEL - same as in the call to XERROR or XERRWV. C KONTRL - the current value of the control flag as set C by a call to XSETF. C C Abstract C XSETF sets the error control flag value to KONTRL. C (KONTRL is an input parameter only.) C The following table shows how each message is treated, C depending on the values of KONTRL and LEVEL. (See XERROR C for description of LEVEL.) C C If KONTRL is zero or negative, no information other than the C message itself (including numeric values, if any) will be C printed. If KONTRL is positive, introductory messages, C trace-backs, etc., will be printed in addition to the message. C C IABS(KONTRL) C LEVEL 0 1 2 C value C 2 fatal fatal fatal C C 1 not printed printed fatal C C 0 not printed printed printed C C -1 not printed printed printed C only only C once once C --Output-- C KONTRL - the new value of KONTRL. If KONTRL is not C defined, it will remain at its original value. C This changed value of control affects only C the current occurrence of the current message. C***REFERENCES JONES R.E., KAHANER D.K., "XERROR, THE SLATEC ERROR- C HANDLING PACKAGE", SAND82-0800, SANDIA LABORATORIES, C 1982. C***ROUTINES CALLED (NONE) C***END PROLOGUE XERCTL CHARACTER*20 MESSG1 C***FIRST EXECUTABLE STATEMENT XERCTL KONTRL=0 RETURN END SUBROUTINE XERPRT(MESSG,NMESSG) C***BEGIN PROLOGUE XERPRT C***DATE WRITTEN 790801 (YYMMDD) C***REVISION DATE 820801 (YYMMDD) C***CATEGORY NO. Z C***KEYWORDS ERROR,XERROR PACKAGE C***AUTHOR JONES, R. E., (SNLA) C***PURPOSE Prints error messages. C***DESCRIPTION C Abstract C Print the Hollerith message in MESSG, of length NMESSG, C on each file indicated by XGETUA. C Latest revision --- 19 MAR 1980 C***REFERENCES JONES R.E., KAHANER D.K., "XERROR, THE SLATEC ERROR- C HANDLING PACKAGE", SAND82-0800, SANDIA LABORATORIES, C 1982. C***ROUTINES CALLED I1MACH,S88FMT,XGETUA C***END PROLOGUE XERPRT INTEGER LUN(5) CHARACTER*(*) MESSG C OBTAIN UNIT NUMBERS AND WRITE LINE TO EACH UNIT C***FIRST EXECUTABLE STATEMENT XERPRT CALL XGETUA(LUN,NUNIT) LENMES = LEN(MESSG) DO 20 KUNIT=1,NUNIT IUNIT = LUN(KUNIT) IF (IUNIT.EQ.0) IUNIT = I1MACH(4) DO 10 ICHAR=1,LENMES,72 LAST = MIN0(ICHAR+71 , LENMES) WRITE (IUNIT,'(1X,A)') MESSG(ICHAR:LAST) 10 CONTINUE 20 CONTINUE RETURN END