Ken
unread,May 2, 2012, 6:19:12 PM5/2/12Sign in to reply to author
Sign in to forward
You do not have permission to delete messages in this group
Either email addresses are anonymous for this group or you need the view member email addresses permission to view the original message
to gg95
Thanks. Here it is. I'll likely be embarrased as to how easy the
problem is.
IMPLICIT NONE
INTEGER M, N, I, J, K
DOUBLE PRECISION ALPHA, BETA, TEMP, T, XNORM
DOUBLE PRECISION SM, Y, DA, ONE, ZERO
PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
PARAMETER ( M=5, N=5 )
DOUBLE PRECISION A(M,N), TAU(N)
INTRINSIC SQRT
TAU = ZERO
A(1,1) = 0.0D0
A(2,1) = 408.0D0
A(3,1) = -25.0D0
A(4,1) = 703.0D0
A(5,1) = 237.0D0
A(1,2) = 9.0D0
A(2,2) = 320.0D0
A(3,2) = 556.0D0
A(4,2) = -155.0D0
A(5,2) = -74.0D0
A(1,3) = 296.0D0
A(2,3) = 480.0D0
A(3,3) = 106.0D0
A(4,3) = -431.0D0
A(5,3) = 0.0D0
A(1,4) = A(1,2)
A(2,4) = A(2,2)
A(3,4) = A(3,2)
A(4,4) = A(4,2)
A(5,4) = A(5,2)
A(1,5) = 190.0D0
A(2,5) = -194.0D0
A(3,5) = 229.0D0
A(4,5) = -389.0D0
A(5,5) = 357.0D0
WRITE (*, '(5F18.10)') ((A(I,J),J=1,N),I=1,M)
WRITE (*, *)
DO 70 J = 1, N-1
ALPHA = A(J,J)
SM = ZERO
DO 20 I = J+1, M
SM = SM + A(I,J)**2
20 CONTINUE
XNORM = SQRT(SM)
IF (XNORM.EQ.ZERO) THEN
TAU(J) = ZERO
ELSE
BETA = SQRT( ALPHA**2 + XNORM**2 )
IF (ALPHA.GE.ZERO) THEN
BETA = -BETA
END IF
TAU(J) = (BETA - ALPHA) / BETA
T = -TAU(J)
DA = ONE / (ALPHA - BETA)
DO 30 I = J+1, M
A(I,J) = A(I,J) * DA
30 CONTINUE
A(J,J) = ONE
DO 60 K = J+1, N
Y = ZERO
DO 40 I = J, M
Y = Y + A(I,J) * A(I,K)
40 CONTINUE
IF (Y.NE.ZERO) THEN
TEMP = T * Y
*---different result if write statement commented out ---
WRITE(*,*)
*--------------------------------------------------------
DO 50 I = J, M
A(I,K) = A(I,K) + A(I,J) * TEMP
50 CONTINUE
END IF
60 CONTINUE
A(J,J) = BETA
END IF
70 CONTINUE
WRITE (*, *)
WRITE (*, '(5F18.8)') ((A(I,J),J=1,N),I=1,M)
WRITE (*, *)
WRITE (*, '(5F18.8)') (TAU(J),J=1,N)
END