C MAIN PROGRAM INTEGER LUNIT C ALLOW 5000 UNDERFLOWS. CALL TRAPS(0,0,5001,0,0) C C OUTPUT UNIT NUMBER C LUNIT = 6 C CALL SQRXX(LUNIT) STOP END SUBROUTINE SQRXX(LUNIT) INTEGER LUNIT REAL RX(20,10),R(10,10),XROW(10),Z(10,4),YROW(10),S(10) REAL SCALE,TINY,HUGE,RHO(2),C(10),SMACH INTEGER N,P,LDX,LDR,LDZ,CASE LOGICAL NOTWRT,OFLOW,UFLOW COMMON SCALE,NOTWRT,OFLOW,UFLOW LDZ = 10 LDX = 20 LDR = 10 NOTWRT = .TRUE. OFLOW = .FALSE. UFLOW = .FALSE. TINY = SMACH(2) HUGE = SMACH(3) SCALE = 1.0E0 DO 60 CASE = 1, 3 GO TO (10, 20, 30), CASE 10 CONTINUE N = 20 P = 10 GO TO 40 20 CONTINUE N = 10 P = 4 GO TO 40 30 CONTINUE N = 10 P = 1 40 CONTINUE CALL SGETRX(RX,LDX,N,P,S) WRITE (LUNIT,130) CASE,N,P IF (NOTWRT) GO TO 50 WRITE (LUNIT,160) CALL SARRAY(RX,LDX,N,P,P,LUNIT) 50 CONTINUE CALL SDRUD1(R,LDR,RX,LDX,N,P,XROW,C,S,LUNIT) CALL SDRUD2(R,LDR,RX,LDX,P,XROW,YROW,Z,LDZ,RHO,C,S,LUNIT) CALL SDRDD(R,LDR,RX,LDX,P,Z,LDZ,XROW,YROW,RHO,C,S,LUNIT) 60 CONTINUE CASE = 4 N = 10 P = 4 WRITE (LUNIT,130) CASE,N,P WRITE (LUNIT,140) CALL SGETRX(RX,LDX,N,P,S) DO 80 J = 1, P DO 70 I = 1, N RX(I,J) = HUGE*RX(I,J) 70 CONTINUE 80 CONTINUE SCALE = HUGE OFLOW = .TRUE. IF (NOTWRT) GO TO 90 WRITE (LUNIT,160) CALL SARRAY(RX,LDX,N,P,P,LUNIT) 90 CONTINUE CALL SDRUD1(R,LDR,RX,LDX,N,P,XROW,C,S,LUNIT) CALL SDRUD2(R,LDR,RX,LDX,P,XROW,YROW,Z,LDZ,RHO,C,S,LUNIT) CALL SDRDD(R,LDR,RX,LDX,P,Z,LDZ,XROW,YROW,RHO,C,S,LUNIT) OFLOW = .FALSE. CASE = 5 N = 10 P = 4 WRITE (LUNIT,130) CASE,N,P WRITE (LUNIT,150) CALL SGETRX(RX,LDX,N,P,S) DO 110 J = 1, P DO 100 I = 1, N RX(I,J) = TINY*RX(I,J) 100 CONTINUE 110 CONTINUE SCALE = TINY UFLOW = .TRUE. IF (NOTWRT) GO TO 120 WRITE (LUNIT,160) CALL SARRAY(RX,LDX,N,P,P,LUNIT) 120 CONTINUE CALL SDRUD1(R,LDR,RX,LDX,N,P,XROW,C,S,LUNIT) CALL SDRUD2(R,LDR,RX,LDX,P,XROW,YROW,Z,LDZ,RHO,C,S,LUNIT) CALL SDRDD(R,LDR,RX,LDX,P,Z,LDZ,XROW,YROW,RHO,C,S,LUNIT) RETURN C 130 FORMAT (11H1 CASE =, I2, 5X, 3HN =, I2, 5X, 3HP =, I2 /////) 140 FORMAT (22H OVERFLOW TEST /////) 150 FORMAT (24H UNDERFLOW TEST /////) 160 FORMAT (5H RX) END SUBROUTINE SARRAY(A,LDA,M,N,NNL,LUNIT) INTEGER LDA,M,N,NNL,LUNIT REAL A(LDA,1) C X C FORTRAN IABS,MIN0 C INTEGER I,J,K,KU,NL NL = IABS(NNL) IF (NNL .LT. 0) GO TO 30 DO 20 I = 1, M WRITE (LUNIT,70) DO 10 K = 1, N, NL KU = MIN0(K+NL-1,N) WRITE (LUNIT,70) (A(I,J), J = K, KU) 10 CONTINUE 20 CONTINUE GO TO 60 30 CONTINUE DO 50 J = 1, N WRITE (LUNIT,70) DO 40 K = 1, M, NL KU = MIN0(K+NL-1,M) WRITE (LUNIT,70) (A(I,J), I = K, KU) 40 CONTINUE 50 CONTINUE 60 CONTINUE RETURN C 70 FORMAT (1H , 4(2E13.6, 4X)) END SUBROUTINE SDRDD(R,LDR,RX,LDX,P,Z,LDZ,XROW,YROW,RHO,C,S,LUNIT) INTEGER LDR,LDX,LDZ,P,LUNIT REAL R(LDR,1),RX(LDX,1),Z(LDZ,1),XROW(1),YROW(1),S(1) REAL RHO(1),SCALE,C(1) LOGICAL NOTWRT,OFLOW,UFLOW COMMON SCALE,NOTWRT,OFLOW,UFLOW INTEGER I,J,JP1,PM1 WRITE (LUNIT,50) CALL SCHDD(R,LDR,P,XROW,Z,LDZ,2,YROW,RHO,C,S,INFO) PM1 = P - 1 IF (PM1 .LT. 1) GO TO 30 DO 20 J = 1, PM1 JP1 = J + 1 DO 10 I = JP1, P R(I,J) = 0.0E0 10 CONTINUE 20 CONTINUE 30 CONTINUE IF (NOTWRT) GO TO 40 WRITE (LUNIT,80) CALL SARRAY(R,LDR,P,P,P,LUNIT) WRITE (LUNIT,60) CALL SARRAY(Z,LDZ,P,2,2,LUNIT) WRITE (LUNIT,70) (RHO(I), I = 1, 2) 40 CONTINUE CALL SMPDD(R,LDR,P,RX,LDX,Z,Z(1,3),LDZ,RHO,LUNIT) RETURN 50 FORMAT ( ///// * 46H STEP THREE DOWNDATING XROW,YROW AND Z, ///) 60 FORMAT ( /// 4H Z) 70 FORMAT ( /// 6H RHO // 1H , 2E13.6) 80 FORMAT (4H R) END SUBROUTINE SDRUD1(R,LDR,RX,LDX,N,P,XROW,C,S,LUNIT) INTEGER N,P,LDR,LDX,LUNIT REAL R(LDR,1),RX(LDX,1),XROW(1),S(1) REAL C(1) LOGICAL NOTWRT,OFLOW REAL RELM,XELM,Y,Z REAL XMAX,RMAX,TEST,T,SCALE,SMACH COMMON SCALE,NOTWRT,OFLOW,UFLOW DO 20 I = 1, P DO 10 J = 1, P R(I,J) = 0.0E0 10 CONTINUE 20 CONTINUE DO 40 I = 1, N DO 30 J = 1, P XROW(J) = RX(I,J) 30 CONTINUE CALL SCHUD(R,LDR,P,XROW,Z,1,0,Y,Y,C,S) 40 CONTINUE WRITE (LUNIT,100) IF (NOTWRT) GO TO 50 WRITE (LUNIT,120) CALL SARRAY(R,LDR,P,P,P,LUNIT) 50 CONTINUE RMAX = 0.0E0 XMAX = 0.0E0 DO 90 I = 1, P DO 80 J = 1, I RELM = 0.0E0 XELM = 0.0E0 NIM = MIN0(I,J) DO 60 K = 1, NIM RELM = RELM + R(K,I)/SCALE*(R(K,J)/SCALE) 60 CONTINUE DO 70 K = 1, N XELM = XELM + RX(K,I)/SCALE*(RX(K,J)/SCALE) 70 CONTINUE T = AMAX1(ABS(XELM),0.0E0) XMAX = AMAX1(XMAX,T) T = AMAX1(ABS(XELM-RELM),0.0E0) RMAX = AMAX1(RMAX,T) 80 CONTINUE 90 CONTINUE TEST = RMAX/XMAX/SMACH(1) WRITE (LUNIT,110) TEST RETURN C 100 FORMAT ( ///// 25H STEP ONE UPDATING X ///) 110 FORMAT ( /// 15H STATISTICS // * 42H RH*R ............................, E10.2) 120 FORMAT (4H R) END SUBROUTINE SDRUD2(R,LDR,RX,LDX,P,XROW,YROW,Z,LDZ,RHO,C,S,LUNIT) INTEGER LDR,LDX,LDZ,P,LUNIT REAL R(LDR,1),RX(LDX,1),XROW(1),YROW(1),Z(LDZ,1),S(1) REAL RHO(1),C(1) REAL SCALE LOGICAL NOTWRT,OFLOW,UFLOW COMMON SCALE,NOTWRT,OFLOW,UFLOW DO 20 I = 1, P DO 10 J = 1, P RX(I,J) = R(I,J) 10 CONTINUE 20 CONTINUE DO 40 I = 1, P Z(I,1) = 0.0E0 DO 30 J = I, P Z(I,1) = Z(I,1) + R(I,J) 30 CONTINUE Z(I,2) = Z(I,1) + SCALE Z(I,3) = Z(I,1) Z(I,4) = Z(I,2) 40 CONTINUE DO 60 I = 1, P XROW(I) = 0.0E0 DO 50 J = 1, I XROW(I) = XROW(I) + R(J,I) 50 CONTINUE 60 CONTINUE YROW(1) = 0.0E0 DO 70 I = 1, P YROW(1) = YROW(1) + XROW(I) 70 CONTINUE YROW(2) = YROW(1) - SCALE RHO(1) = SCALE RHO(2) = SCALE IF (NOTWRT) GO TO 80 WRITE (LUNIT,100) CALL SARRAY(XROW,P,P,1,-P,LUNIT) WRITE (LUNIT,110) CALL SARRAY(Z,LDZ,P,2,2,LUNIT) WRITE (LUNIT,120) CALL SARRAY(YROW,2,2,1,-2,LUNIT) 80 CONTINUE CALL SCHUD(R,LDR,P,XROW,Z,LDZ,2,YROW,RHO,C,S) WRITE (LUNIT,130) IF (NOTWRT) GO TO 90 WRITE (LUNIT,150) CALL SARRAY(R,LDR,P,P,P,LUNIT) WRITE (LUNIT,110) CALL SARRAY(Z,DLZ,P,2,2,LUNIT) WRITE (LUNIT,140) (RHO(I), I = 1, 2) 90 CONTINUE CALL SMPUD(R,LDR,RX,LDX,P,XROW,Z,LDZ,RHO,LUNIT) RETURN C 100 FORMAT ( ///// 7H XROW) 110 FORMAT ( /// 4H Z) 120 FORMAT ( /// 7H YROW) 130 FORMAT ( ///// 41H STEP TWO UPDATING XROW, YROW AND Z ///) 140 FORMAT ( /// 6H RHO // 1H , 2E13.6) 150 FORMAT (4H R) END SUBROUTINE SGETRX(X,LDX,N,P,S) INTEGER N,P,LDX REAL X(LDX,1),S(1) INTEGER PD2,PD2D1,I PD2 = MAX0(P/2,1) PD2D1 = PD2 + 1 DO 10 I = 1, PD2 S(I) = 1.0E0 10 CONTINUE IF (P .LT. PD2D1) GO TO 30 DO 20 I = PD2D1, P S(I) = 0.5E0 20 CONTINUE 30 CONTINUE CALL SXGEN(X,LDX,N,P,S) RETURN END SUBROUTINE SMPDD(R,LDR,P,ROLD,LDRD,Z,ZUP,LDZ,RHO,LUNIT) INTEGER LDR,P,LDRD,LDZ,LUNIT REAL R(LDR,1),ROLD(LDRD,1),Z(LDZ,1),ZUP(LDZ,1) REAL RHO(1),SMACH REAL T INTEGER I,J REAL TR,TZ(2),TRHO(2),MAXOLD,MAXREL,SCALE LOGICAL NOTWRT,OFLOW,UFLOW COMMON SCALE,NOTWRT,OFLOW,UFLOW MAXOLD = 0.0E0 MAXREL = 0.0E0 DO 20 J = 1, P DO 10 I = 1, J T = ROLD(I,J) MAXOLD = AMAX1(MAXOLD,AMAX1(ABS(T),0.0E0)) T = ROLD(I,J) - R(I,J) MAXREL = AMAX1(MAXREL,AMAX1(ABS(T),0.0E0)) 10 CONTINUE 20 CONTINUE TR = MAXREL/MAXOLD/SMACH(1) DO 40 I = 1, 2 MAXOLD = 0.0E0 MAXREL = 0.0E0 DO 30 J = 1, P T = ZUP(J,I) MAXOLD = AMAX1(MAXOLD,AMAX1(ABS(T),0.0E0)) T = ZUP(J,I) - Z(J,I) MAXREL = AMAX1(MAXREL,AMAX1(ABS(T),0.0E0)) 30 CONTINUE TZ(I) = MAXREL/MAXOLD/SMACH(1) TRHO(I) = ABS(RHO(I)/SCALE-1.0E0)/SMACH(1) 40 CONTINUE WRITE (LUNIT,50) TR,TZ,TRHO RETURN C 50 FORMAT ( /// 25H STATSTICS STEP THREE // * 33H R ....................., E10.2 / * 33H Z(1) .................., E10.2 / * 33H Z(2) .................., E10.2 / * 33H RHO(1) ................, E10.2 / * 33H RHO(2) ................, E10.2) END SUBROUTINE SMPUD(R,LDR,RX,LDX,P,XROW,Z,LDZ,RHO,LUNIT) INTEGER LDR,LDX,P,LDZ REAL R(LDR,1),RX(LDX,1),XROW(1),Z(LDZ,1) REAL RHO(1) REAL RELM,XELM,TMP1(10),TMP REAL TEST1,TEST2(2),TEST3,TEST4,MAXR,MAXX,T LOGICAL NOTWRT,OFLOW,UFLOW REAL SCALE,SMACH COMMON SCALE,NOTWRT,OFLOW,UFLOW MAXR = 0.0E0 MAXX = 0.0E0 DO 30 I = 1, P DO 20 J = 1, I RELM = 0.0E0 XELM = 0.0E0 NIM = MIN0(I,J) DO 10 K = 1, NIM RELM = RELM + R(K,I)/SCALE*(R(K,J)/SCALE) XELM = XELM + RX(K,I)/SCALE*(RX(K,J)/SCALE) 10 CONTINUE XELM = XELM + XROW(I)/SCALE*(XROW(J)/SCALE) T = AMAX1(ABS(XELM),0.0E0) MAXX = AMAX1(MAXX,T) T = AMAX1(ABS(XELM-RELM),0.0E0) MAXR = AMAX1(MAXR,T) 20 CONTINUE 30 CONTINUE TEST1 = MAXR/MAXX/SMACH(1) TMP = 0.0E0 DO 50 I = 1, P TMP1(I) = 0.0E0 DO 40 J = I, P TMP1(I) = TMP1(I) + RX(I,J) 40 CONTINUE TMP = TMP + XROW(I) 50 CONTINUE DO 80 K = 1, 2 MAXR = 0.0E0 MAXX = 0.0E0 DO 70 I = 1, P RELM = 0.0E0 XELM = 0.0E0 DO 60 J = 1, I RELM = RELM + R(J,I)/SCALE*(Z(J,K)/SCALE) XELM = XELM + RX(J,I)/SCALE*(TMP1(J)/SCALE) 60 CONTINUE XELM = XELM + XROW(I)/SCALE*(TMP/SCALE) T = AMAX1(ABS(XELM-RELM),0.0E0) MAXR = AMAX1(MAXR,T) T = AMAX1(ABS(XELM),0.0E0) MAXX = AMAX1(MAXX,T) 70 CONTINUE TEST2(K) = MAXR/MAXX/SMACH(1) 80 CONTINUE TEST3 = ABS(RHO(1)/SCALE-1.0E0)/SMACH(1) T = SQRT(2.0E0+FLOAT(P)) TEST4 = ABS(RHO(2)/SCALE-T)/T/SMACH(1) WRITE (LUNIT,90) TEST1,TEST2,TEST3,TEST4 RETURN 90 FORMAT ( /// 24H STATISTICS STEP TWO // * 31H RH*R ................, E10.2 / * 31H RH*Z(1) ............., E10.2 / * 31H RH*Z(2) ............., E10.2 / * 31H RHO(1) .............., E10.2 / * 31H RHO(2) .............., E10.2) END SUBROUTINE SXGEN(X,LDX,N,P,S) INTEGER LDX,N,P REAL X(LDX,1),S(1) INTEGER I,J,M,MP1 REAL FN,FP REAL FAC,RU,T FP = FLOAT(P) FN = FLOAT(N) M = MIN0(N,P) RU = 1.0E0 RU = COS(6.28E0/FLOAT(M+1)) FAC = 1.0E0/FP DO 20 I = 1, M FAC = FAC*RU DO 10 J = 1, P X(I,J) = -2.0E0*S(I)*FAC 10 CONTINUE X(I,I) = X(I,I) + FP*S(I)*FAC 20 CONTINUE IF (M .GE. N) GO TO 50 MP1 = M + 1 DO 40 J = 1, P DO 30 I = MP1, N X(I,J) = 0.0E0 30 CONTINUE 40 CONTINUE 50 CONTINUE DO 80 J = 1, P T = 0.0E0 DO 60 I = 1, N T = T + X(I,J) 60 CONTINUE DO 70 I = 1, N X(I,J) = X(I,J) - 2.0E0*T/FN 70 CONTINUE 80 CONTINUE RETURN END