real a(100,10),b(100,10),ab(400,10),bb(400,10) real y(100),w(100),yy(100) real e(400),d(400),dd(400) real add(100,10) real ee(400),e2(400) 1000 read(5,10001)n,k km1=k-1 10001 format(2i3) if (n.eq.0)go to 999 write(6,10002)n,k 10002 format(" n",i5,"k=",i5) c do 20 i=1,n c do 10 j=1,n c a(i,j)=00.0 c b(i,j)=0.0 c10 continue c20 continue c do 30 i=1,n c a(i,i)=(i+1)*(i+1) c b(i,i)=(i+1)*(i+1) c do 25 j=1,km1 c if (i.le.j) go to 25 c a(i,i-j)=0.5 c a(i-j,i)=a(i,i-j) c b(i-j,i)=1.0 c b(i,i-j)=b(i-j,i) c25 continue c30 continue c do 40 i=1,n c do 35 j=1,k c if (j-k+i.le.0) go to 35 c ab(i,j)=a(i,j-k+i) c bb(i,j)=b(i,j-k+i) c35 continue c40 continue nm=400 c t1=xectim(0) c call reduc(nm,n,a,b,y,ierr) c call tred1(nm,n,a,dd,ee,e2) c t1=xectim(0)-t1 c write(6,673)t1 c673 format(" time for reduc and tred1",e15.5) c t2=xectim(0) c call tql1(n,dd,ee,ierr) c t2=xectim(0)-t2 c write(6,674)t2 c674 format(" time for tql1",e15.5) c ia=10 c write(6,84)(dd(i),i=1,n) c84 format(" eigenvalues",5e15.4) do 424 i=1,n ab(i,k)=(i+1)*(i+1) bb(i,k)=ab(i,k) do 423 j=1,km1 ab(i,j)=.5 bb(i,j)=.5 423 continue 424 continue c t3=xectim(0) call bndchl(nm,n,k,bb) call redbnd(nm,n,k,k,ab,bb) c t3=xectim(0)-t3 write(6,675)t3 675 format(" time for redub",e15.5) c t5=xectim(0) call bandr2(nm,n,k,ab,d,e,e2) c t5=xectim(0)-t5 write(6,676)t5 676 format(" time for bandr",e15.5) go to 1000 c call tql1(n,d,e,ierr) c do 80 i=1,n c write(6,81)d(i),dd(i) c80 continue 81 format(" eigenvalues",2e15.5) 999 stop go to 1000 end SUBROUTINE REDBND(NMAX,N,MBA,MBB,A,B) C A AND B ARE BAND MATRICES. A IS SYMMETRIC AND B IS LOWER- C TRIANGULAR. AT EXIT FROM THIS SUBROUTINE A CONTAINS C, A C SYMMETRIC BAND MATRIX ORTHOGONALLY SIMILAR TO B-INV*A*B'INV , C WHERE B-INV IS THE INVERSE OF B AND B'INV IS THE INVERSE OF C THE TRANSPOSE OF B. C C DESCRIPTION OF PARAMETERS IN THE CALLING SEQUENCE C NMAX AN INTEGER VARIABLE CONTAINING THE FIRST DIMENSION C OF THE ARRAYS A AND B. C C N THE ORDER OF THE MATRICES A AND B. C C MBA THE BANDWIDTH OF A. USING ORDINARY MATRIX NOTATION C A=0 WHENEVER \I-J\ >=MBA. C C C MBB THE BANDWIDTH OF B. (SEE DESCRIPTION OF MBA.). C C A A DOUBLY DIMENSIONED ARRAY IN WHICH THE MATRIX A IS C STORED SUCH THAT THE MATRIX ELEMENT FROM ROW K AND C COLUMN L IS STORED IN THE (K,MBA-K+L) POSITION IN THE C OF THE ARRAY A. THE ARRAY A CONTAINS ONLY THE LOWER C TRIANGLE OF THE MATRIX A. C C B A DOUBLY DIMENSIONED ARRAY IN WHICH THE MATRIX B IS C STORED SUCH THAT THE MATRIX ELEMENT FROM ROW K AND C COLUMN L IS STORED IN THE (K,MBB-K+L) POSITION IN THE C ARRAY B. THE ARRAY B CONTAINS ONLY THE LOWER TRIANGLE C OF THE MATRIX B. C C OUTPUT: THE MATRIX C IS OVERWRITTEN IN THE ARRAY A. C THE MATRIX B IS UNCHANGED. C ARRAY ELEMENTS (A(1,J),J=1,MBA-1) ARE USED AS SCRATCH SPACE. C C SUBROUTINES CALLED: SQUASH,CHASE,ROTATE,ROT. C C RESTRICTIONS: MBB<= MBA. INTEGER NMAX,N,MBB,MBA REAL A(NMAX,MBA),B(NMAX,MBB) REAL AJJ,AXK,BJJ,AXJ MA=MBA-1 MB=MBB-1 BJJ=B(1,MBB) DO 1 K=2,MBA 1 A(K,MBA-K+1)=A(K,MBA-K+1)/BJJ A(1,MBA)=A(1,MBA)/(BJJ*BJJ) DO 22 J=2,N BJJ=B(J,MBB) C:::::::UPDATE ROW(J) OF A. C LET X=(0,0,...,B,...,B,0,...,0) C THEN (AX)=SUM FROM L=MAX(J-MB,K-MA) TO MIN(J-1,K+MA) C OF A(*B. C NOTE: ONLY COMPONENTS J-MA-MB,...,J+MA-1 OF AX ARE NON-ZERO. C A <--(A - (AX))/B FOR J]=K. C A <-- (A -2*(AX) + X*AX)/B**2. :::::::::::: C AJJ=0. C:::::::COMPUTE THE J-TH COMPONENT OF AX, AX. ::::::::::::::::::::: AXJ=0. KK1=MAX0(1,J-MB) KK2=J-1 DO 2 KK=KK1,KK2 2 AXJ=AXJ+ A(J,MBA-J+KK)*B(J,MBB-J+KK) C:::::::NOTE: WE ARE NOW FREE TO MODIFY ROW AND COLUMN J OF A.::::::: LENU=0 C IF(J-MA.LE.1) GO TO 9 C:::::::COMPUTE THE FILL, A AND SAVE IN A(1,K-K1+1) FOR K < J-MA.:: K1=MAX0(1,J-MA-MB) K2= J-MA-1 LENU=K2-K1+1 DO 8 K=K1,K2 C:::::::::::COMPUTE THE K-TH COMPONENT OF AX,AX.::::::::::::::::::: AXK=0. KK1=MAX0(1,K-MA,J-MB) IF(KK1.GT.K) GO TO 6 DO 5 KK=KK1,K 5 AXK=AXK + A(K,MBA-K+KK)*B(J,MBB-J+KK) 6 KK1=MAX0(K+1,J-MB) KK2=K+MA DO 7 KK=KK1,KK2 7 AXK=AXK + A(KK,MBA-KK+K)*B(J,MBB-J+KK) C:::::::::::COMPUTE X*AX.::::::::::::::::::::::::::::::::::::::::::::: IF(J-K .LE. MB) AJJ=AJJ+B(J,MBB-J+K)*AXK 8 A(1,K-K1+1)= -AXK/BJJ C C:::::::UPDATE A FOR J-MA<= K < J. :::::::::::::::::::::::::::::: 9 K1=MAX0(1,J-MA) K2=J-1 DO 17 K=K1,K2 C:::::::::COMPUTE THE K-TH COMPONENT OF AX, AX.:::::::::::::::::::: AXK=0 KK1=MAX0(J-MB,K-MA,1) IF(KK1.GT.K) GO TO 13 DO 12 KK=KK1,K 12 AXK =AXK+ A(K,MBA-K+KK)*B(J,MBB-J+KK) 13 KK1=MAX0(K+1,J-MB) KK2=MIN0(J-1,K+MA) IF(KK1.GT.KK2) GO TO 16 DO 15 KK=KK1,KK2 15 AXK=AXK+ A(KK,MBA-KK+K)*B(J,MBB-J+KK) C:::::::::::COMPUTE X*AX.::::::::::::::::::::::::::::::::::::::::::::: 16 IF(J-K .LE. MB) AJJ=AJJ+B(J,MBB-J+K)*AXK 17 A(J,MBA-J+K)=(A(J,MBA-J+K) - AXK)/BJJ C C:::::::UPDATE A.:::::::::::::::::::::::::::::::::::::::::::::::: A(J,MBA)= (A(J,MBA) - 2*AXJ + AJJ)/(BJJ*BJJ) C C:::::::UPDATE A FOR J < K < J-MA.::::::::::::::::::::::::::::::: K1=J+1 K2=MIN0(N,J+MA-1) IF(K1.GT.K2) GO TO 21 DO 20 K=K1,K2 C:::::::::::COMPUTE THE K-TH COMPONENT OF AX, AX.:::::::::::::::::: AXK=0 KK1=MAX0(1,K-MA,J-MB) KK2=J-1 DO 19 KK=KK1,KK2 19 AXK=AXK+ A(K,MBA-K+KK)*B(J,MBB-J+KK) 20 A(K,MBA-K+J)= (A(K,MBA-K+J) - AXK)/BJJ C C:::::::UPDATE A. :::::::::::::::::::::::::::::::::::::::::::: 21 IF(J+MA.LE.N) A(J+MA,1)= A(J+MA,1)/BJJ C C:::::::TAKE CARE OF FILL, IF ANY. ::::::::::::::::::::::::::::::::::: IF(LENU.GT.0) CALL SQUASH(NMAX,N,MBA,A,J,LENU) 22 CONTINUE RETURN END SUBROUTINE CHASE(NMAX,N,MB,A,J,GG) C PROGRAM TO ELIMINATE FILL-IN FROM A BAND SYMMETRIC MATRIX USING C ROTATIONS IN PLANES (K,K+1)FOR K <= J. THE FILL-IN,WHOSE VALUE C IS IN GG,IS AT A. C IT IS ASSUMED THAT J+MB <= N. C OPS = ORDER(J). C DESCRIPTION OF PARAMETERS IN THE CALLING SEQUENCE: C NMAX THE ROW DIMENSION OF THE ARRAY A. C C N THE ORDER OF THE MATRIX A. C C MB THE BANDWIDTH OF THE MATRIX A. USING ORDINARY MATRIX C NOTATION, A = 0 WHENEVER \I-J\ >= MB. C C A A DOUBLY DIMENSIONED ARRAY IN WHICH THE MATRIX A IS C STORED SUCH THAT THE MATRIX ELEMENT FROM ROW K AND C COLUMN L IS STORED IN THE (K,MBA-K+L) POSITION IN THE C ARRAY A. THE ARRAY A CONTAINS ONLY THE LOWER TRIANGLE C OF THE MATRIX A. C SUBROUTINES CALLED: ROTATE,ROT INTEGER NMAX,MB,M,J,K,KK REAL A(NMAX,MB),G,GG,C,S,R,T REAL ABS,SQRT M=MB-1 G=GG C:::::FILL IS CHASED FROM COLUMN J BACKWARD IN STEPS OF M OFF THE C UPPER LEFT CORNER OF A.::::::::::::::::::::::::::::::::::::::::: DO 9 KK=1,J,M C::::::::THIS LOOP WILL ALWAYS TERMINATE VIA ONE OF THE EXTRA- C ORDINARY EXITS, I.E. KK WILL NEVER EXCEED J.::::::::::::::::: K=J-KK+1 IF(G.EQ.0.) GO TO 10 CALL ROT(G,A(K+MB,1),C,S,R) A(K+MB,1)=R CALL ROTATE(NMAX,N,MB,A,K,C,S) C::::::::EXIT IF(K-M.LT.1) ::::::::::::::::::::::::::::::::::::::::::: IF(K-M.LT.1) GO TO 10 C::::::::COMPUTE FILL AT A CAUSED BY A.::::::::::::::: G= -S*A(K,1) 9 A(K,1)=C*A(K,1) 10 CONTINUE RETURN END SUBROUTINE SQUASH(NMAX,N,MB,A,J,LENU) C PROGRAM TO ELIMINATE FILL-IN IN A BAND SYMMETRIC MATRIX USING C ROTATIONS IN PLANES (K,K+1) FOR K<=J-MB. THE FILL-IN, WHOSE C VALUES ARE IN A(1,1),...,A(1,LENU), ARE AT MATRIX POSITIONS C A, FOR K= J-MB-LENU+1,...,J-MB. C IT IS ASSUMED THAT: 1) MB <= J <=N, AND 0 = 0 WHENEVER \I-J\ >= MB. C C A A DOUBLY DIMENSIONED ARRAY IN WHICH THE MATRIX A IS C STORED SUCH THAT THE MATRIX ELEMENT FROM ROW K AND C COLUMN L IS STORED IN THE (K,MB-K+L) POSITION IN THE C ARRAY A. THE ARRAY A CONTAINS ONLY THE LOWER TRIANGLE C OF THE MATRIX A. C THE VALUES OF THE FILL ELEMENTS ARE ASSUMED TO BE STORED C IN ELEMENTS (1,1),...,(1,LENU) OF THE ARRAY A. C C SUBROUTINES CALLED: CHASE,ROTATE,ROT INTEGER NMAX,MB,M,J,LENU,KFIRST,KLAST REAL A(NMAX,MB),GG,C,S,R REAL ABS,SQRT M=MB-1 KFIRST=J-M-LENU KLAST=J-MB IF(LENU.GT.1) GO TO 2 C::::::::ONLY ONE ELEMENT TO CHASE. :::::::::::::::::::::::::::::::::: CALL CHASE(NMAX,N,MB,A,KLAST,A(1,1)) RETURN C:::::::::::::MOVE FILL-IN IN COLUMN J TO FORM EXTRA SUPER-DIAGONAL. ::: 2 K2=KLAST-1 DO 15 K=KFIRST,K2 KK=K-KFIRST+1 CALL ROT(A(1,KK) ,A(1,KK+1),C,S,R) A(1,KK+1)=R A(1,KK)= S*A(K+MB,1) A(K+MB,1)= C*A(K+MB,1) CALL ROTATE(NMAX,N,MB,A,K,C,S) IF(K.LE.M) GO TO 15 C::::::::::::::ELIMINATE EXTRA FILL. ::::::::::::::::::::::::::::::::: 11 KM=K-M GG=-S*A(KM+M,1) A(KM+M,1)= C*A(KM+M,1) CALL CHASE(NMAX,N,MB,A,KM,GG) 15 CONTINUE C C::::::::LOOP TO ELIMINATE EXTRA SUPER-DIAGONAL, A, K=KFIRST TO C KLAST AND WHOSE VALUES ARE STORE IN A(1,1) TO A(1,LENU). :::::: 17 IF(KLAST.LT.KFIRST) GO TO 37 CALL ROT(A(1,LENU),A(KLAST+MB,1),C,S,R) A(KLAST+MB,1)=R CALL ROTATE(NMAX,N,MB,A,KLAST,C,S) IF(KLAST.LE.M) GO TO 33 C:::::::::::::::::COMPUTE FILL AT A, C CAUSED BY A.::::::::::::::::::::::::: GG=-S*A(KLAST,1) A(KLAST,1)= C*A(KLAST,1) K=KLAST-M C:::::::::::::::::TAKE CARE OF THE FILL. ::::::::::::::::::::::::::::: IF(K.GE.KFIRST-1) GO TO 27 C:::::::::::::::::::NEW FILL IS SEPARARTED FROM OLD. ::::::::::::::::: CALL CHASE(NMAX,N,MB,A,K,GG) GO TO 33 C:::::::::::::::::::NEW FILL IS NOT SEPARATED FROM OLD SO LEAP-FROG C BACK.::::::::: 27 DO 29 KKK=2,LENU KU=LENU-KKK+2 29 A(1,KU)= A(1,KU-1) A(1,1)=GG KFIRST=KFIRST-1 33 CONTINUE KLAST=KLAST-1 LENU=KLAST-KFIRST+1 GO TO 17 37 CONTINUE RETURN END SUBROUTINE ROTATE(NMAX,N,MB,A,K,C,S) C APPLY A ROTATION IN PLANE(K,K+1) TO A BAND-SYMMETRIC MATRIX. C ONLY ROWS AND COLUMNS BETWEEN I1 AND I2 ARE CHANGED.. CHANGES C TO THE MATRIX CAUSED BY THE ROTATION BUT OCCURRING OUT SIDE C THESE LIMITS ARE CALCULATED ELSEWHERE. THESE CHANGES INCLUDE C FILL-IN OUTSIDE THE BAND STRUCTURE OF THE MATRIX. C DESCRIPTION OF PARAMETERS IN CALLING SEQUENCE: C NMAX THE ROW DIMENSION OF THE ARRAY A. C C N THE ORDER OF THE MATRIX A. C C MB THE BANDWIDTH OF THE MATRIX A. USING ORDINARY MATRIX C NOTATION, A = 0 WHENEVER \I-J\ >= MB. C C A A DOUBLY DIMENSIONED ARRAY IN WHICH THE MATRIX A IS C STORED SUCH THAT THE MATRIX ELEMENT FROM ROW K AND C COLUMN L IS STORED IN THE (K,MB-K+L) POSITION IN THE C ARRAY A. THE ARRAY A CONTAINS ONLY THE LOWER TRIANGLE C OF THE MATRIX A. C OPS.= (I2 -I1-1)(4M + 4A) + 10M + 7A <= 2(M-1)(4M+2A)+10M+7A INTEGER NMAX,N,MB,M,I1,I2,IR,IM,IM2 REAL A(NMAX,MB),C,S,R,G,H,T,AK,AK1,AKK,AKK1,C2,S2,SC M=MB-1 I1=MAX0(1,K-M+1) I2=MIN0(N,K+M) IF(I1.GE.K) GO TO 10 KM1=K-1 DO 8 IR=I1,KM1 IM=MB-K+IR AK=A(K,IM) AK1=A(K+1,IM-1) A(K,IM)= C*AK + S*AK1 8 A(K+1,IM-1)= -S*AK + C*AK1 10 CONTINUE K2=K+2 IF(K2.GT.I2) GO TO 18 DO 16 IR=K2,I2 IM=MB-IR+K AK= A(IR,IM) AK1=A(IR,IM+1) A(IR,IM)= C*AK + S*AK1 16 A(IR,IM+1)=-S*AK + C*AK1 18 CONTINUE C2=C*C S2=S*S SC=S*C AKK=A(K,MB) AK1=A(K+1,MB) AKK1=A(K+1,M) A(K,MB)= C2*AKK + 2*SC*AKK1 + S2*AK1 A(K+1,MB) = C2*AK1 - 2*SC*AKK1 + S2*AKK A(K+1,M)= (C2-S2)*AKK1 + SC*(AK1-AKK) RETURN END SUBROUTINE ROT(G,H,C,S,R) C FIND THE PARAMETERS OF A PLANE ROTATION WHICH MAPS THE VECTOR C (G,H) TO (0,R)#&W8ERE R = SQRT(G*G+R*R). C AND S ARE CHOSEN SO C THAT C*G + S*H = 0.0. REAL C,S,R,H,G,T REAL SQRT,ABS C THERE ARE FOUR CASES: IF (G.EQ.0) GO TO 1 IF(H.EQ.0.0) GO TO 2 IF (ABS(G).GT.ABS(H)) GO TO 3 C SECANT FORMULA T=-G/H R = SQRT(1.D0 + T*T) C= 1.D0/R S=C*T R=R*H RETURN C TRIVIAL CASE 1 C=1.D0 S=0.D0 R=H RETURN C INTERCHANGE COORDINATES. 2 C=0.D0 S=1.D0 R=G RETURN C COSECANT 3 T=-H/G R=SQRT(1.D0+T*T) S=1.D0/R C=S*T R=-R*G RETURN END SUBROUTINE BNDCHL(NM,N,MB,B) C SUBROUTINE TO FIND THE CHOLESKY DECOMPOSITION OF A REAL , C SYMMETRIC, POSITIVE-DEFINITE BAND-MATRIX B. C DESCRIPTION OF PARAMETERS: C NM THE ROW DIMENSION OF THE ARRAY B. C C N ORDER OF THE MATRX B. C C MB HALF-BANDWIDTH OF THE MATRIX B. THE MATRIX ELEMENT C (I,J) OF B IS ZERO WHEN \I-J\ >= MB. C C B A DOUBLY DIMENSIONED ARRAY IN WHICH THE MATRIX B IS C STORED SUCH THAT THE MATRIX ELEMENT FROM ROW K AND C COLUMN L IS STORED IN THE (K,MB-K+L) POSITION IN THE C ARRAY B. ONLY THE LOWER TRIANGLE OF THE MATRIX IS C STORED IN THE ARRAY. C C OUTPUT; THE MATRIX L, THE CHOLESKY FACTOR OF B IS OVER-WRITTEN IN C THE ARRAY. B. C IF THE MATRIX B IS NOT POSITIVE DEFINITE, (POSSIBLY C BECAUSE OF ROUNDING ERRORS), ARRAY ELEMENT B(N,MB) IS SET C TO ZERO. REAL B(NM,MB),SUM,D,SQRT M=MB-1 DO 4 J=1,N JPM=MIN0(J+M,N) JM1=J-1 DO 4 I=J,JPM IMM=MAX0(I-M,1) SUM=B(I,MB-I+J) IF(IMM.GT.JM1) GO TO 2 DO 1 K=IMM,JM1 1 SUM=SUM- B(I,MB-I+K)*B(J,MB-J+K) 2 IF(I.NE.J) GO TO 3 C *******ERROR EXIT**************** IF(SUM.LE.0.) GO TO 99 D=SQRT(SUM) B(I,MB)=D GO TO 4 3 CONTINUE B(I,MB-I+J)=SUM/D 4 CONTINUE RETURN C **ERROR RETURN******************* 99 B(N,MB)=0. RETURN END subroutine bandr2(nm, n, mb, a, d, c, s) integer n, mb, nm real a(nm, mb), d(n), c(n), s(n) integer jxbeg, jlast, jcolp1, i, j, l integer m, j1, jl, ml, jx, lm1 integer nm2, mp1, min0, jxm1, gbeg, lbeg integer gend, iend, lend, jcol, itim, lpml integer lpjx, l2end, l3end real r1mach, e, g, r, u, cc real ce, ch, al, dq, sc, sg real sh, ss, sc2, dia, abs, tiny real sqrt, amin1, amax1 c c this is a vectorized bandr routine c c mb is number of nonzero diagonals below diagonal+diagonal c ml = mb-1 mp1 = mb+1 tiny = 10.*r1mach(1) gbeg = mb+2 do 100 i=1,n d(i)=0.0 100 continue gend = 1 m = mb nm2 = n-2 do 14 i = 1, nm2 jlast = min0(ml, n-i+2) do 13 j = 1, jlast gbeg = gbeg-1 jl = 1 lbeg = gbeg lend = gend iend = 2 if (j .eq. 1) iend = 1 c call check(n,a,nm,d,mb) do 12 itim = 1, iend if (lbeg .gt. lend) goto 9 c write(6,93)i,j,itim,lbeg,lend 93 format(" i,j,itim,lbeg,lend",5i5) c write(6,94)(d(k),k=1,n) 94 format(" d",5e15.5) cdir$ ivdep do 1 l = lbeg, lend, m lm1 = l-1 al = a(lm1, jl) r = amin1(abs(d(l)), abs(al)) sc = amax1(abs(d(l)), abs(al)) sc2 = amax1(sc, tiny) r = r/sc2 dq = sqrt(r*r+1.) dia = sc2*dq cc = (al+tiny)/dia ss = d(l)/dia d(l)=0.0 a(l-1, jl) = sc*dq c(l) = cc s(l) = ss e = a(lm1, m) g = a(l, m) ce = cc*e sg = ss*g sh = ss*a(l, ml) ch = cc*a(l, ml) a(lm1, m) = cc*ce+2.*ss*ch+ss*sg a(l, m) = e+g-a(lm1, m) a(l, ml) = ss*(cc*(g-e)-sh)+cc*ch 1 continue c c do row transformations jxbeg = 2 if (itim .eq. 2) jxbeg = j+1 if (ml .lt. jxbeg) goto 4 do 3 jx = jxbeg, ml do 2 l = lbeg, lend, m u = a(l-1, jx) a(l-1, jx) = c(l)*u+s(l)*a(l, jx-1) a(l, jx-1) = (-s(l))*u+c(l)*a(l, jx-1) 2 continue 3 continue c c do column transformations c 4 jxbeg = 2 do 7 jx = jxbeg, ml jcol = m-jx jcolp1=jcol+1 jxm1 = jx-1 l2end = min0(lend, n-jxm1) if (l2end .lt. lbeg) goto 6 do 5 l = lbeg, l2end, m lpjx = l+jxm1 u = a(lpjx, jcol) a(lpjx, jcol) = c(l)*u+s(l)*a(lpjx, jcolp1) a(lpjx, jcolp1) = (-s(l))*u+c(l)*a(lpjx, 1 jcolp1) 5 continue 6 continue 7 continue c c creation of new elements l3end = min0(lend, n-ml) if (l3end .lt. lbeg) goto 9 do 8 l = lbeg, l3end, m lpml = l+ml d(lpml) = a(lpml, 1)*s(l) a(lpml, 1) = a(lpml, 1)*c(l) 8 continue 9 if (itim .ne. 1 .or. j .le. 1) goto 10 lbeg = i-j+mp1 lend = lbeg jl = j d(lbeg) = a(lbeg, j-1) a(lbeg,j-1)=0.0 10 continue 12 continue if (j .eq. 1) gbeg = gbeg+m gend= min0(n,gend+ml) 13 continue 14 continue do 15 i = 2, n c(i) = a(i, ml) s(i) = c(i)*c(i) d(i) = a(i, m) 15 continue d(1) = a(1, m) c(1) = 0 s(1) = 0. return end real function r1mach(i) c c single-precision machine constants c c r1mach(1) = b**(emin-1), the smallest positive magnitude. c c r1mach(2) = b**emax*(1 - b**(-t)), the largest magnitude. c c r1mach(3) = b**(-t), the smallest relative spacing. c c r1mach(4) = b**(1-t), the largest relative spacing. c c r1mach(5) = log10(b) 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. c c where possible, octal or hexadecimal constants have been used c to specify the constants exactly which has in some cases c required the use of equivalent integer arrays. 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 machine constants for the burroughs 1700 system. c 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 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 cdc 6000/7000 series. c 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 c c data rmach(1) / 200004000000000000000b / c data rmach(2) / 577777777777777777777b / c data rmach(3) / 377214000000000000000b / c data rmach(4) / 377224000000000000000b / 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 card - c static rmach(5) c c data small/20k,0/,large/77777k,177777k/ c data right/35420k,0/,diver/36020k,0/ c data log10/40423k,42023k/ c c machine constants for the harris 220 c 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 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 ibm 360/370 series, c the xerox sigma 5/7/9 and the sel systems 85/86. c 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 pdp-10 (ka or ki processor). c 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's supporting c 32-bit integers (expressed in integer and octal). c 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's supporting c 16-bit integers (expressed in integer and octal). c 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 univac 1100 series. c 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 c data small(1) / 128 / data large(1) / -32769 / data right(1) / 13440 / data diver(1) / 13568 / data log10(1) / 547045274 / c if (i .lt. 1 .or. i .gt. 5) c 1 call seterr(24hr1mach - i out of bounds,24,1,2) c r1mach = rmach(i) return c end real a(100,10),b(100,10),ab(400,10),bb(400,10) real y(100),w(100),yy(100) real e(400),d(400),dd(400) real add(100,10) real ee(400),e2(400) 1000 read(5,10001)n,k km1=k-1 10001 format(2i3) if (n.eq.0)go to 999 write(6,10002)n,k 10002 format(" n",i5,"k=",i5) c do 20 i=1,n c do 10 j=1,n c a(i,j)=00.0 c b(i,j)=0.0 c10 continue c20 continue c do 30 i=1,n c a(i,i)=(i+1)*(i+1) c b(i,i)=(i+1)*(i+1) c do 25 j=1,km1 c if (i.le.j) go to 25 c a(i,i-j)=0.5 c a(i-j,i)=a(i,i-j) c b(i-j,i)=1.0 c b(i,i-j)=b(i-j,i) c25 continue c30 continue c do 40 i=1,n c do 35 j=1,k c if (j-k+i.le.0) go to 35 c ab(i,j)=a(i,j-k+i) c bb(i,j)=b(i,j-k+i) c35 continue c40 continue nm=400 c t1=xectim(0) c call reduc(nm,n,a,b,y,ierr) c call tred1(nm,n,a,dd,ee,e2) c t1=xectim(0)-t1 c write(6,673)t1 c673 format(" time for reduc and tred1",e15.5) c t2=xectim(0) c call tql1(n,dd,ee,ierr) c t2=xectim(0)-t2 c write(6,674)t2 c674 format(" time for tql1",e15.5) c ia=10 c write(6,84)(dd(i),i=1,n) c84 format(" eigenvalues",5e15.4) do 424 i=1,n ab(i,k)=(i+1)*(i+1) bb(i,k)=ab(i,k) do 423 j=1,km1 ab(i,j)=.5 bb(i,j)=.5 423 continue 424 continue c t3=xectim(0) call bndchl(nm,n,k,bb) call redbnd(nm,n,k,k,ab,bb) c t3=xectim(0)-t3 write(6,675)t3 675 format(" time for redub",e15.5) c t5=xectim(0) call bandr2(nm,n,k,ab,d,e,e2) c t5=xectim(0)-t5 write(6,676)t5 676 format(" time for bandr",e15.5) go to 1000 c call tql1(n,d,e,ierr) c do 80 i=1,n c write(6,81)d(i),dd(i) c80 continue 81 format(" eigenvalues",2e15.5) 999 stop go to 1000 end SUBROUTINE REDBND(NMAX,N,MBA,MBB,A,B) C A AND B ARE BAND MATRICES. A IS SYMMETRIC AND B IS LOWER- C TRIANGULAR. AT EXIT FROM THIS SUBROUTINE A CONTAINS C, A C SYMMETRIC BAND MATRIX ORTHOGONALLY SIMILAR TO B-INV*A*B'INV , C WHERE B-INV IS THE INVERSE OF B AND B'INV IS THE INVERSE OF C THE TRANSPOSE OF B. C C DESCRIPTION OF PARAMETERS IN THE CALLING SEQUENCE C NMAX AN INTEGER VARIABLE CONTAINING THE FIRST DIMENSION C OF THE ARRAYS A AND B. C C N THE ORDER OF THE MATRICES A AND B. C C MBA THE BANDWIDTH OF A. USING ORDINARY MATRIX NOTATION C A=0 WHENEVER \I-J\ >=MBA. C C C MBB THE BANDWIDTH OF B. (SEE DESCRIPTION OF MBA.). C C A A DOUBLY DIMENSIONED ARRAY IN WHICH THE MATRIX A IS C STORED SUCH THAT THE MATRIX ELEMENT FROM ROW K AND C COLUMN L IS STORED IN THE (K,MBA-K+L) POSITION IN THE C OF THE ARRAY A. THE ARRAY A CONTAINS ONLY THE LOWER C TRIANGLE OF THE MATRIX A. C C B A DOUBLY DIMENSIONED ARRAY IN WHICH THE MATRIX B IS C STORED SUCH THAT THE MATRIX ELEMENT FROM ROW K AND C COLUMN L IS STORED IN THE (K,MBB-K+L) POSITION IN THE C ARRAY B. THE ARRAY B CONTAINS ONLY THE LOWER TRIANGLE C OF THE MATRIX B. C C OUTPUT: THE MATRIX C IS OVERWRITTEN IN THE ARRAY A. C THE MATRIX B IS UNCHANGED. C ARRAY ELEMENTS (A(1,J),J=1,MBA-1) ARE USED AS SCRATCH SPACE. C C SUBROUTINES CALLED: SQUASH,CHASE,ROTATE,ROT. C C RESTRICTIONS: MBB<= MBA. INTEGER NMAX,N,MBB,MBA REAL A(NMAX,MBA),B(NMAX,MBB) REAL AJJ,AXK,BJJ,AXJ MA=MBA-1 MB=MBB-1 BJJ=B(1,MBB) DO 1 K=2,MBA 1 A(K,MBA-K+1)=A(K,MBA-K+1)/BJJ A(1,MBA)=A(1,MBA)/(BJJ*BJJ) DO 22 J=2,N BJJ=B(J,MBB) C:::::::UPDATE ROW(J) OF A. C LET X=(0,0,...,B,...,B,0,...,0) C THEN (AX)=SUM FROM L=MAX(J-MB,K-MA) TO MIN(J-1,K+MA) C OF A(*B. C NOTE: ONLY COMPONENTS J-MA-MB,...,J+MA-1 OF AX ARE NON-ZERO. C A <--(A - (AX))/B FOR J]=K. C A <-- (A -2*(AX) + X*AX)/B**2. :::::::::::: C AJJ=0. C:::::::COMPUTE THE J-TH COMPONENT OF AX, AX. ::::::::::::::::::::: AXJ=0. KK1=MAX0(1,J-MB) KK2=J-1 DO 2 KK=KK1,KK2 2 AXJ=AXJ+ A(J,MBA-J+KK)*B(J,MBB-J+KK) C:::::::NOTE: WE ARE NOW FREE TO MODIFY ROW AND COLUMN J OF A.::::::: LENU=0 C IF(J-MA.LE.1) GO TO 9 C:::::::COMPUTE THE FILL, A AND SAVE IN A(1,K-K1+1) FOR K < J-MA.:: K1=MAX0(1,J-MA-MB) K2= J-MA-1 LENU=K2-K1+1 DO 8 K=K1,K2 C:::::::::::COMPUTE THE K-TH COMPONENT OF AX,AX.::::::::::::::::::: AXK=0. KK1=MAX0(1,K-MA,J-MB) IF(KK1.GT.K) GO TO 6 DO 5 KK=KK1,K 5 AXK=AXK + A(K,MBA-K+KK)*B(J,MBB-J+KK) 6 KK1=MAX0(K+1,J-MB) KK2=K+MA DO 7 KK=KK1,KK2 7 AXK=AXK + A(KK,MBA-KK+K)*B(J,MBB-J+KK) C:::::::::::COMPUTE X*AX.::::::::::::::::::::::::::::::::::::::::::::: IF(J-K .LE. MB) AJJ=AJJ+B(J,MBB-J+K)*AXK 8 A(1,K-K1+1)= -AXK/BJJ C C:::::::UPDATE A FOR J-MA<= K < J. :::::::::::::::::::::::::::::: 9 K1=MAX0(1,J-MA) K2=J-1 DO 17 K=K1,K2 C:::::::::COMPUTE THE K-TH COMPONENT OF AX, AX.:::::::::::::::::::: AXK=0 KK1=MAX0(J-MB,K-MA,1) IF(KK1.GT.K) GO TO 13 DO 12 KK=KK1,K 12 AXK =AXK+ A(K,MBA-K+KK)*B(J,MBB-J+KK) 13 KK1=MAX0(K+1,J-MB) KK2=MIN0(J-1,K+MA) IF(KK1.GT.KK2) GO TO 16 DO 15 KK=KK1,KK2 15 AXK=AXK+ A(KK,MBA-KK+K)*B(J,MBB-J+KK) C:::::::::::COMPUTE X*AX.::::::::::::::::::::::::::::::::::::::::::::: 16 IF(J-K .LE. MB) AJJ=AJJ+B(J,MBB-J+K)*AXK 17 A(J,MBA-J+K)=(A(J,MBA-J+K) - AXK)/BJJ C C:::::::UPDATE A.:::::::::::::::::::::::::::::::::::::::::::::::: A(J,MBA)= (A(J,MBA) - 2*AXJ + AJJ)/(BJJ*BJJ) C C:::::::UPDATE A FOR J < K < J-MA.::::::::::::::::::::::::::::::: K1=J+1 K2=MIN0(N,J+MA-1) IF(K1.GT.K2) GO TO 21 DO 20 K=K1,K2 C:::::::::::COMPUTE THE K-TH COMPONENT OF AX, AX.:::::::::::::::::: AXK=0 KK1=MAX0(1,K-MA,J-MB) KK2=J-1 DO 19 KK=KK1,KK2 19 AXK=AXK+ A(K,MBA-K+KK)*B(J,MBB-J+KK) 20 A(K,MBA-K+J)= (A(K,MBA-K+J) - AXK)/BJJ C C:::::::UPDATE A. :::::::::::::::::::::::::::::::::::::::::::: 21 IF(J+MA.LE.N) A(J+MA,1)= A(J+MA,1)/BJJ C C:::::::TAKE CARE OF FILL, IF ANY. ::::::::::::::::::::::::::::::::::: IF(LENU.GT.0) CALL SQUASH(NMAX,N,MBA,A,J,LENU) 22 CONTINUE RETURN END SUBROUTINE CHASE(NMAX,N,MB,A,J,GG) C PROGRAM TO ELIMINATE FILL-IN FROM A BAND SYMMETRIC MATRIX USING C ROTATIONS IN PLANES (K,K+1)FOR K <= J. THE FILL-IN,WHOSE VALUE C IS IN GG,IS AT A. C IT IS ASSUMED THAT J+MB <= N. C OPS = ORDER(J). C DESCRIPTION OF PARAMETERS IN THE CALLING SEQUENCE: C NMAX THE ROW DIMENSION OF THE ARRAY A. C C N THE ORDER OF THE MATRIX A. C C MB THE BANDWIDTH OF THE MATRIX A. USING ORDINARY MATRIX C NOTATION, A = 0 WHENEVER \I-J\ >= MB. C C A A DOUBLY DIMENSIONED ARRAY IN WHICH THE MATRIX A IS C STORED SUCH THAT THE MATRIX ELEMENT FROM ROW K AND C COLUMN L IS STORED IN THE (K,MBA-K+L) POSITION IN THE C ARRAY A. THE ARRAY A CONTAINS ONLY THE LOWER TRIANGLE C OF THE MATRIX A. C SUBROUTINES CALLED: ROTATE,ROT INTEGER NMAX,MB,M,J,K,KK REAL A(NMAX,MB),G,GG,C,S,R,T REAL ABS,SQRT M=MB-1 G=GG C:::::FILL IS CHASED FROM COLUMN J BACKWARD IN STEPS OF M OFF THE C UPPER LEFT CORNER OF A.::::::::::::::::::::::::::::::::::::::::: DO 9 KK=1,J,M C::::::::THIS LOOP WILL ALWAYS TERMINATE VIA ONE OF THE EXTRA- C ORDINARY EXITS, I.E. KK WILL NEVER EXCEED J.::::::::::::::::: K=J-KK+1 IF(G.EQ.0.) GO TO 10 CALL ROT(G,A(K+MB,1),C,S,R) A(K+MB,1)=R CALL ROTATE(NMAX,N,MB,A,K,C,S) C::::::::EXIT IF(K-M.LT.1) ::::::::::::::::::::::::::::::::::::::::::: IF(K-M.LT.1) GO TO 10 C::::::::COMPUTE FILL AT A CAUSED BY A.::::::::::::::: G= -S*A(K,1) 9 A(K,1)=C*A(K,1) 10 CONTINUE RETURN END SUBROUTINE SQUASH(NMAX,N,MB,A,J,LENU) C PROGRAM TO ELIMINATE FILL-IN IN A BAND SYMMETRIC MATRIX USING C ROTATIONS IN PLANES (K,K+1) FOR K<=J-MB. THE FILL-IN, WHOSE C VALUES ARE IN A(1,1),...,A(1,LENU), ARE AT MATRIX POSITIONS C A, FOR K= J-MB-LENU+1,...,J-MB. C IT IS ASSUMED THAT: 1) MB <= J <=N, AND 0 = 0 WHENEVER \I-J\ >= MB. C C A A DOUBLY DIMENSIONED ARRAY IN WHICH THE MATRIX A IS C STORED SUCH THAT THE MATRIX ELEMENT FROM ROW K AND C COLUMN L IS STORED IN THE (K,MB-K+L) POSITION IN THE C ARRAY A. THE ARRAY A CONTAINS ONLY THE LOWER TRIANGLE C OF THE MATRIX A. C THE VALUES OF THE FILL ELEMENTS ARE ASSUMED TO BE STORED C IN ELEMENTS (1,1),...,(1,LENU) OF THE ARRAY A. C C SUBROUTINES CALLED: CHASE,ROTATE,ROT INTEGER NMAX,MB,M,J,LENU,KFIRST,KLAST REAL A(NMAX,MB),GG,C,S,R REAL ABS,SQRT M=MB-1 KFIRST=J-M-LENU KLAST=J-MB IF(LENU.GT.1) GO TO 2 C::::::::ONLY ONE ELEMENT TO CHASE. :::::::::::::::::::::::::::::::::: CALL CHASE(NMAX,N,MB,A,KLAST,A(1,1)) RETURN C:::::::::::::MOVE FILL-IN IN COLUMN J TO FORM EXTRA SUPER-DIAGONAL. ::: 2 K2=KLAST-1 DO 15 K=KFIRST,K2 KK=K-KFIRST+1 CALL ROT(A(1,KK) ,A(1,KK+1),C,S,R) A(1,KK+1)=R A(1,KK)= S*A(K+MB,1) A(K+MB,1)= C*A(K+MB,1) CALL ROTATE(NMAX,N,MB,A,K,C,S) IF(K.LE.M) GO TO 15 C::::::::::::::ELIMINATE EXTRA FILL. ::::::::::::::::::::::::::::::::: 11 KM=K-M GG=-S*A(KM+M,1) A(KM+M,1)= C*A(KM+M,1) CALL CHASE(NMAX,N,MB,A,KM,GG) 15 CONTINUE C C::::::::LOOP TO ELIMINATE EXTRA SUPER-DIAGONAL, A, K=KFIRST TO C KLAST AND WHOSE VALUES ARE STORE IN A(1,1) TO A(1,LENU). :::::: 17 IF(KLAST.LT.KFIRST) GO TO 37 CALL ROT(A(1,LENU),A(KLAST+MB,1),C,S,R) A(KLAST+MB,1)=R CALL ROTATE(NMAX,N,MB,A,KLAST,C,S) IF(KLAST.LE.M) GO TO 33 C:::::::::::::::::COMPUTE FILL AT A, C CAUSED BY A.::::::::::::::::::::::::: GG=-S*A(KLAST,1) A(KLAST,1)= C*A(KLAST,1) K=KLAST-M C:::::::::::::::::TAKE CARE OF THE FILL. ::::::::::::::::::::::::::::: IF(K.GE.KFIRST-1) GO TO 27 C:::::::::::::::::::NEW FILL IS SEPARARTED FROM OLD. ::::::::::::::::: CALL CHASE(NMAX,N,MB,A,K,GG) GO TO 33 C:::::::::::::::::::NEW FILL IS NOT SEPARATED FROM OLD SO LEAP-FROG C BACK.::::::::: 27 DO 29 KKK=2,LENU KU=LENU-KKK+2 29 A(1,KU)= A(1,KU-1) A(1,1)=GG KFIRST=KFIRST-1 33 CONTINUE KLAST=KLAST-1 LENU=KLAST-KFIRST+1 GO TO 17 37 CONTINUE RETURN END SUBROUTINE ROTATE(NMAX,N,MB,A,K,C,S) C APPLY A ROTATION IN PLANE(K,K+1) TO A BAND-SYMMETRIC MATRIX. C ONLY ROWS AND COLUMNS BETWEEN I1 AND I2 ARE CHANGED.. CHANGES C TO THE MATRIX CAUSED BY THE ROTATION BUT OCCURRING OUT SIDE C THESE LIMITS ARE CALCULATED ELSEWHERE. THESE CHANGES INCLUDE C FILL-IN OUTSIDE THE BAND STRUCTURE OF THE MATRIX. C DESCRIPTION OF PARAMETERS IN CALLING SEQUENCE: C NMAX THE ROW DIMENSION OF THE ARRAY A. C C N THE ORDER OF THE MATRIX A. C C MB THE BANDWIDTH OF THE MATRIX A. USING ORDINARY MATRIX C NOTATION, A = 0 WHENEVER \I-J\ >= MB. C C A A DOUBLY DIMENSIONED ARRAY IN WHICH THE MATRIX A IS C STORED SUCH THAT THE MATRIX ELEMENT FROM ROW K AND C COLUMN L IS STORED IN THE (K,MB-K+L) POSITION IN THE C ARRAY A. THE ARRAY A CONTAINS ONLY THE LOWER TRIANGLE C OF THE MATRIX A. C OPS.= (I2 -I1-1)(4M + 4A) + 10M + 7A <= 2(M-1)(4M+2A)+10M+7A INTEGER NMAX,N,MB,M,I1,I2,IR,IM,IM2 REAL A(NMAX,MB),C,S,R,G,H,T,AK,AK1,AKK,AKK1,C2,S2,SC M=MB-1 I1=MAX0(1,K-M+1) I2=MIN0(N,K+M) IF(I1.GE.K) GO TO 10 KM1=K-1 DO 8 IR=I1,KM1 IM=MB-K+IR AK=A(K,IM) AK1=A(K+1,IM-1) A(K,IM)= C*AK + S*AK1 8 A(K+1,IM-1)= -S*AK + C*AK1 10 CONTINUE K2=K+2 IF(K2.GT.I2) GO TO 18 DO 16 IR=K2,I2 IM=MB-IR+K AK= A(IR,IM) AK1=A(IR,IM+1) A(IR,IM)= C*AK + S*AK1 16 A(IR,IM+1)=-S*AK + C*AK1 18 CONTINUE C2=C*C S2=S*S SC=S*C AKK=A(K,MB) AK1=A(K+1,MB) AKK1=A(K+1,M) A(K,MB)= C2*AKK + 2*SC*AKK1 + S2*AK1 A(K+1,MB) = C2*AK1 - 2*SC*AKK1 + S2*AKK A(K+1,M)= (C2-S2)*AKK1 + SC*(AK1-AKK) RETURN END SUBROUTINE ROT(G,H,C,S,R) C FIND THE PARAMETERS OF A PLANE ROTATION WHICH MAPS THE VECTOR C (G,H) TO (0,R)#&W8ERE R = SQRT(G*G+R*R). C AND S ARE CHOSEN SO C THAT C*G + S*H = 0.0. REAL C,S,R,H,G,T REAL SQRT,ABS C THERE ARE FOUR CASES: IF (G.EQ.0) GO TO 1 IF(H.EQ.0.0) GO TO 2 IF (ABS(G).GT.ABS(H)) GO TO 3 C SECANT FORMULA T=-G/H R = SQRT(1.D0 + T*T) C= 1.D0/R S=C*T R=R*H RETURN C TRIVIAL CASE 1 C=1.D0 S=0.D0 R=H RETURN C INTERCHANGE COORDINATES. 2 C=0.D0 S=1.D0 R=G RETURN C COSECANT 3 T=-H/G R=SQRT(1.D0+T*T) S=1.D0/R C=S*T R=-R*G RETURN END SUBROUTINE BNDCHL(NM,N,MB,B) C SUBROUTINE TO FIND THE CHOLESKY DECOMPOSITION OF A REAL , C SYMMETRIC, POSITIVE-DEFINITE BAND-MATRIX B. C DESCRIPTION OF PARAMETERS: C NM THE ROW DIMENSION OF THE ARRAY B. C C N ORDER OF THE MATRX B. C C MB HALF-BANDWIDTH OF THE MATRIX B. THE MATRIX ELEMENT C (I,J) OF B IS ZERO WHEN \I-J\ >= MB. C C B A DOUBLY DIMENSIONED ARRAY IN WHICH THE MATRIX B IS C STORED SUCH THAT THE MATRIX ELEMENT FROM ROW K AND C COLUMN L IS STORED IN THE (K,MB-K+L) POSITION IN THE C ARRAY B. ONLY THE LOWER TRIANGLE OF THE MATRIX IS C STORED IN THE ARRAY. C C OUTPUT; THE MATRIX L, THE CHOLESKY FACTOR OF B IS OVER-WRITTEN IN C THE ARRAY. B. C IF THE MATRIX B IS NOT POSITIVE DEFINITE, (POSSIBLY C BECAUSE OF ROUNDING ERRORS), ARRAY ELEMENT B(N,MB) IS SET C TO ZERO. REAL B(NM,MB),SUM,D,SQRT M=MB-1 DO 4 J=1,N JPM=MIN0(J+M,N) JM1=J-1 DO 4 I=J,JPM IMM=MAX0(I-M,1) SUM=B(I,MB-I+J) IF(IMM.GT.JM1) GO TO 2 DO 1 K=IMM,JM1 1 SUM=SUM- B(I,MB-I+K)*B(J,MB-J+K) 2 IF(I.NE.J) GO TO 3 C *******ERROR EXIT**************** IF(SUM.LE.0.) GO TO 99 D=SQRT(SUM) B(I,MB)=D GO TO 4 3 CONTINUE B(I,MB-I+J)=SUM/D 4 CONTINUE RETURN C **ERROR RETURN******************* 99 B(N,MB)=0. RETURN END subroutine bandr2(nm, n, mb, a, d, c, s) integer n, mb, nm real a(nm, mb), d(n), c(n), s(n) integer jxbeg, jlast, jcolp1, i, j, l integer m, j1, jl, ml, jx, lm1 integer nm2, mp1, min0, jxm1, gbeg, lbeg integer gend, iend, lend, jcol, itim, lpml integer lpjx, l2end, l3end real r1mach, e, g, r, u, cc real ce, ch, al, dq, sc, sg real sh, ss, sc2, dia, abs, tiny real sqrt, amin1, amax1 c c this is a vectorized bandr routine c c mb is number of nonzero diagonals below diagonal+diagonal c ml = mb-1 mp1 = mb+1 tiny = 10.*r1mach(1) gbeg = mb+2 do 100 i=1,n d(i)=0.0 100 continue gend = 1 m = mb nm2 = n-2 do 14 i = 1, nm2 jlast = min0(ml, n-i+2) do 13 j = 1, jlast gbeg = gbeg-1 jl = 1 lbeg = gbeg lend = gend iend = 2 if (j .eq. 1) iend = 1 c call check(n,a,nm,d,mb) do 12 itim = 1, iend if (lbeg .gt. lend) goto 9 c write(6,93)i,j,itim,lbeg,lend 93 format(" i,j,itim,lbeg,lend",5i5) c write(6,94)(d(k),k=1,n) 94 format(" d",5e15.5) cdir$ ivdep do 1 l = lbeg, lend, m lm1 = l-1 al = a(lm1, jl) r = amin1(abs(d(l)), abs(al)) sc = amax1(abs(d(l)), abs(al)) sc2 = amax1(sc, tiny) r = r/sc2 dq = sqrt(r*r+1.) dia = sc2*dq cc = (al+tiny)/dia ss = d(l)/dia d(l)=0.0 a(l-1, jl) = sc*dq c(l) = cc s(l) = ss e = a(lm1, m) g = a(l, m) ce = cc*e sg = ss*g sh = ss*a(l, ml) ch = cc*a(l, ml) a(lm1, m) = cc*ce+2.*ss*ch+ss*sg a(l, m) = e+g-a(lm1, m) a(l, ml) = ss*(cc*(g-e)-sh)+cc*ch 1 continue c c do row transformations jxbeg = 2 if (itim .eq. 2) jxbeg = j+1 if (ml .lt. jxbeg) goto 4 do 3 jx = jxbeg, ml do 2 l = lbeg, lend, m u = a(l-1, jx) a(l-1, jx) = c(l)*u+s(l)*a(l, jx-1) a(l, jx-1) = (-s(l))*u+c(l)*a(l, jx-1) 2 continue 3 continue c c do column transformations c 4 jxbeg = 2 do 7 jx = jxbeg, ml jcol = m-jx jcolp1=jcol+1 jxm1 = jx-1 l2end = min0(lend, n-jxm1) if (l2end .lt. lbeg) goto 6 do 5 l = lbeg, l2end, m lpjx = l+jxm1 u = a(lpjx, jcol) a(lpjx, jcol) = c(l)*u+s(l)*a(lpjx, jcolp1) a(lpjx, jcolp1) = (-s(l))*u+c(l)*a(lpjx, 1 jcolp1) 5 continue 6 continue 7 continue c c creation of new elements l3end = min0(lend, n-ml) if (l3end .lt. lbeg) goto 9 do 8 l = lbeg, l3end, m lpml = l+ml d(lpml) = a(lpml, 1)*s(l) a(lpml, 1) = a(lpml, 1)*c(l) 8 continue 9 if (itim .ne. 1 .or. j .le. 1) goto 10 lbeg = i-j+mp1 lend = lbeg jl = j d(lbeg) = a(lbeg, j-1) a(lbeg,j-1)=0.0 10 continue 12 continue if (j .eq. 1) gbeg = gbeg+m gend= min0(n,gend+ml) 13 continue 14 continue do 15 i = 2, n c(i) = a(i, ml) s(i) = c(i)*c(i) d(i) = a(i, m) 15 continue d(1) = a(1, m) c(1) = 0 s(1) = 0. return end real function r1mach(i) c c single-precision machine constants c c r1mach(1) = b**(emin-1), the smallest positive magnitude. c c r1mach(2) = b**emax*(1 - b**(-t)), the largest magnitude. c c r1mach(3) = b**(-t), the smallest relative spacing. c c r1mach(4) = b**(1-t), the largest relative spacing. c c r1mach(5) = log10(b) 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. c c where possible, octal or hexadecimal constants have been used c to specify the constants exactly which has in some cases c required the use of equivalent integer arrays. 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 machine constants for the burroughs 1700 system. c 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 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 cdc 6000/7000 series. c 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 c c data rmach(1) / 200004000000000000000b / c data rmach(2) / 577777777777777777777b / c data rmach(3) / 377214000000000000000b / c data rmach(4) / 377224000000000000000b / 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 card - c static rmach(5) c c data small/20k,0/,large/77777k,177777k/ c data right/35420k,0/,diver/36020k,0/ c data log10/40423k,42023k/ c c machine constants for the harris 220 c 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 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 ibm 360/370 series, c the xerox sigma 5/7/9 and the sel systems 85/86. c 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 pdp-10 (ka or ki processor). c 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's supporting c 32-bit integers (expressed in integer and octal). c 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's supporting c 16-bit integers (expressed in integer and octal). c 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 univac 1100 series. c 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 c data small(1) / 128 / data large(1) / -32769 / data right(1) / 13440 / data diver(1) / 13568 / data log10(1) / 547045274 / c if (i .lt. 1 .or. i .gt. 5) c 1 call seterr(24hr1mach - i out of bounds,24,1,2) c r1mach = rmach(i) return c end .