C$TEST LYAD C TO RUN AS A MAIN PROGRAM REMOVE NEXT LINE SUBROUTINE LYAD C*********************************************************************** C C TEST OF THE PORT PROGRAMS DSYCE AND FRIENDS C C*********************************************************************** C MAIN PROGRAM INTEGER IWRITE,I1MACH C ALLOW 5000 UNDERFLOWS. C C SET OUTPUT UNIT NUMBER C IWRITE = I1MACH(2) C CALL SPOTS(IWRITE) STOP END SUBROUTINE SPOTS(IWRITE) C IWRITE IS THE OUTPUT UNIT NUMBER C C TESTS C DSYCE;DSYFBS;DSYML;DSYLE;DBPCE;DBPFS;DBPBS;DBPML;DBPLE C C C SUBROUTINES AND FUNCTIONS C C PORT DSYCE,DSYFBS,DSYML,DSYLE,DBPCE,DBPFS,DBPBS,DBPML,DBPLE C EXTERNAL SPOXX,D1MACH C BLAS DAXPY,DDOT,DASUM C FORTRAN DABS,DMAX1,FLOAT,MAX0 C C INTERNAL VARIABLES C INTEGER IPVT(15),IPVTS(15) INTEGER I,IQ(10),I1,J,JB,INDEX INTEGER K,KASE,KB,KBFAIL,KNPD,KOUNT,KPFAIL INTEGER KSUSP(10),LDA,IWRITE,M,N,NPRINT REAL Q(10),QS(10) DOUBLE PRECISION APSAVE(120),AB(15,15),AINV(15,15),ASAVE(15,15) DOUBLE PRECISION AP(120),B(15),DDOT,X(15),XB(15),XEXACT(15) DOUBLE PRECISION ABSAVE(15,15),DENOM DOUBLE PRECISION XP(15),T,Z(15) DOUBLE PRECISION ANORM,AINORM,COND,COND1,AIBNO,EBNORM DOUBLE PRECISION EN,ENORM,EPS,FNORM,RCOND,RCONDB DOUBLE PRECISION RCONDP,RNORM,DASUM,D1MACH,XNORM DOUBLE PRECISION AINVB(15,15),AIS,SC LOGICAL KBF,KPF C LDA = 15 C C WRITE MATRIX AND SOLUTIONS IF N .LE. NPRINT C NPRINT = 3 C WRITE (IWRITE,560) WRITE (IWRITE,1000) C DO 10 I = 1,10 KSUSP(I) = 0 10 CONTINUE KNPD = 0 KPFAIL = 0 KBFAIL = 0 C C SET EPS TO ROUNDING UNIT FOR DOUBLE PRECISION ARITHMETIC C EPS = D1MACH(4) WRITE (IWRITE,570) EPS WRITE (IWRITE,550) C CALL ENTER(1) C START MAIN LOOP C KASE = 1 20 CONTINUE C C GENERATE TEST MATRIX C CALL SPOXX(ASAVE,LDA,N,KASE,IWRITE) C C N = 0 SIGNALS NO MORE TEST MATRICES C C ...EXIT IF (N .LE. 0) GO TO 540 ANORM = 0.0D0 DO 30 J = 1, N ANORM = DMAX1(ANORM,DASUM(N,ASAVE(1,J),1)) 30 CONTINUE WRITE (IWRITE,720) ANORM C IF (N .GT. NPRINT) GO TO 50 WRITE (IWRITE,550) DO 40 I = 1, N WRITE (IWRITE,760) (ASAVE(I,J), J = 1, N) 40 CONTINUE WRITE (IWRITE,550) 50 CONTINUE C C GENERATE EXACT SOLUTION C XEXACT(1) = 1.0D0 IF (N .GE. 2) XEXACT(2) = 0.0D0 IF (N .LE. 2) GO TO 70 DO 60 I = 3, N XEXACT(I) = -XEXACT(I-2) 60 CONTINUE 70 CONTINUE C C C PUT INTO PACKED FORM K = 0 DO 130 J = 1, N DO 120 I = J,N K = K + 1 AP(K) = ASAVE(I,J) APSAVE(K)=AP(K) 120 CONTINUE 130 CONTINUE CALL DSYCE(N,AP,IPVTS,RCONDP) IF (NERROR(IERR).NE.0) CALL ERROFF C FACTOR BAND FORM AND COMPARE C KBF = .FALSE. M = 0 DO 200 J = 1, N DO 190 I = 1, J IF (ASAVE(I,J) .NE. 0.0D0) M = MAX0(M,J-I) 190 CONTINUE 200 CONTINUE C DO 220 J = 1, N I1=MIN0(N,J+M) DO 210 I = J, I1 K =I-J+1 AB(K,J) = ASAVE(I,J) ABSAVE(K,J)=AB(K,J) 210 CONTINUE 220 CONTINUE WRITE (IWRITE,840) M CALL DBPCE(N,M+1,AB,LDA,RCONDB) IF (NERROR(IERR).EQ.0) GO TO 230 CALL ERROFF IF ((IERR).LT.10+N)GO TO 226 WRITE(IWRITE,860)IERR INDEX = IERR - N - 10 WRITE(IWRITE,229)AB(1,INDEX) 229 FORMAT( 19H OFFENDING DIAGONAL,D17.7) GO TO 530 226 WRITE(IWRITE,580) WRITE (IWRITE,930) RCONDP,RCONDB GO TO 530 230 CONTINUE WRITE(IWRITE,930)RCONDP,RCONDB C C COMPUTE INVERSE AND COND1 = TRUE CONDITION C DO 290 J = 1, N DO 280 I = 1, N AINV(I,J)=0.D0 AINVB(I,J)=0.D0 280 CONTINUE AINV(J,J)=1.D0 AINVB(J,J)=1.D0 290 CONTINUE CALL DBPFS(N,M+1,AB,LDA,AINVB,LDA,N) CALL DBPBS(N,M+1,AB,LDA,AINVB,LDA,N) CALL DSYFBS(N,AP,AINV,LDA,N,IPVTS) AINORM=0.D0 DO 310 J=1,N AIS=DASUM(N,AINV(1,J),1) AINORM = DMAX1(AINORM,AIS) 310 CONTINUE COND1 = ANORM*AINORM WRITE (IWRITE,600) COND1 C C GENERATE RIGHT HAND SIDE FOR BOTH SYMMETRIC AND BAND C CALL DSYML(N,APSAVE,XEXACT,B) CALL MOVEFD(N,B,X) CALL DBPML(N,M+1,ABSAVE,LDA,XEXACT,XB) C SOLVE A*X = B C CALL DSYLE(N,APSAVE,X,N,1) IF (NERROR(IRE).NE.0) CALL ERROFF CALL DBPLE(N,M+1,ABSAVE,LDA,XB,N,1) IF (IRE+NERROR(IRB).EQ.0) GO TO 311 IF (IRB.NE.0) CALL ERROFF WRITE(IWRITE,580) GO TO 530 311 CONTINUE C IF (N .GT. NPRINT) GO TO 330 WRITE (IWRITE,610) DO 320 I = 1, N WRITE (IWRITE,790) X(I), XB(I) 320 CONTINUE WRITE (IWRITE,550) 330 CONTINUE C C C COMPUTE ERRORS AND RESIDUALS C E = X - XEXACT C EB = XB - XEXACT C R = B - A*X C F = A - TRANS(R)*R C XNORM = DASUM(N,X,1) ENORM = 0.0D0 EBNORM = 0.E0 DO 460 J = 1, N ENORM = ENORM + DABS(X(J)-XEXACT(J)) EBNORM = EBNORM + DABS(XB(J)-XEXACT(J)) T = -X(J) CALL DAXPY(N,T,ASAVE(1,J),1,B,1) 460 CONTINUE RNORM = DASUM(N,B,1) C C A*INV(A) - I C AINORM = 0.0D0 AIBNO = 0.E0 DO 490 J = 1, N DO 470 I = 1, N B(I) = 0.0D0 XB(I) = 0.0D0 470 CONTINUE DO 480 K = 1, N T = AINV(K,J) CALL DAXPY(N,T,ASAVE(1,K),1,B,1) T=AINVB(K,J) CALL DAXPY(N,T,ASAVE(1,K),1,XB,1) 480 CONTINUE B(J) = B(J) - 1.0D0 XB(J) = XB(J) - 1.0D0 AINORM = DMAX1(AINORM,DASUM(N,B,1)) AIBNO = DMAX1(AIBNO,DASUM(N,XB,1)) 490 CONTINUE FNORM = 0.0D0 ML=M+1 NP1=N+1 MLP1=ML+1 DO 495 J=1,N NUMAX= MIN0(ML,NP1-J) JM1=J-1 IEND=JM1+NUMAX DO 491 I=J,IEND B(I)=-ASAVE(I,J) 491 CONTINUE KBEGIN=MAX0(1,J-M) L=MIN0(J,ML) NUM=MLP1-L IF (JM1.LT.KBEGIN) GO TO 493 DO 492 K=KBEGIN,JM1 SC=AB(L,K)*AB(1,K) CALL DAXPY(NUM,SC,AB(L,K),1,B(J),1) L=L-1 NUM=MIN0(NUM+1,NUMAX) 492 CONTINUE 493 SC=AB(1,J) CALL DAXPY(NUM-1,SC,AB(2,J),1,B(J+1),1) B(J)=B(J)+AB(1,J) FNORM=DMAX1(FNORM,DASUM(NUMAX,B(J),1)) 495 CONTINUE C WRITE (IWRITE,620) ENORM, EBNORM WRITE (IWRITE,630) RNORM WRITE (IWRITE,730) FNORM WRITE (IWRITE,740) AINORM,AIBNO C C COMPUTE TEST RATIOS C Q(1) = RCONDP/COND1 Q(2) = RCONDB/COND1 Q(3) = COND1/RCONDP Q(4) = COND1/RCONDB Q(5) = ENORM/(EPS*RCONDP*XNORM) Q(6) = EBNORM/(EPS*RCONDP*XNORM) DENOM=DMAX1(1.0D2*D1MACH(1),EPS*ANORM*XNORM) Q(7)= RNORM/DENOM DENOM=DMAX1(1.0D2*D1MACH(1),EPS*ANORM) Q(8)=FNORM/DENOM Q(9) = AINORM/(EPS*RCONDP) Q(10) = AIBNO/(EPS*RCONDP) WRITE (IWRITE,550) WRITE (IWRITE,640) WRITE (IWRITE,550) WRITE (IWRITE,690) WRITE (IWRITE,700) WRITE (IWRITE,710) WRITE (IWRITE,550) WRITE (IWRITE,750) (Q(I), I = 1, 10) WRITE (IWRITE,550) C C LOOK FOR SUSPICIOUS RATIOS C QS(1) = 1.0D0 + 4.0D0*EPS QS(2) = QS(1) QS(3) = 10.0D0 QS(4) =QS(3) EN = DBLE(FLOAT(N)) IF (N .EQ. 1) EN = 2.0D0 DO 500 I=5,10 QS(I) = EN 500 CONTINUE KOUNT = 0 DO 520 I = 1, 10 IQ(I) = 0 IF (Q(I) .LE. QS(I)) GO TO 510 IQ(I) = 1 KSUSP(I) = KSUSP(I) + 1 KOUNT = KOUNT + 1 510 CONTINUE 520 CONTINUE IF (KOUNT .EQ. 0) WRITE (IWRITE,980) IF (KOUNT .NE. 0) WRITE (IWRITE,990) (IQ(I), I = 1,10) WRITE (IWRITE,550) 530 CONTINUE C WRITE (IWRITE,650) KASE = KASE + 1 GO TO 20 540 CONTINUE C C FINISH MAIN LOOP C C SUMMARY C WRITE (IWRITE,660) KASE = KASE - 1 WRITE (IWRITE,670) KASE WRITE (IWRITE,680) KSUSP WRITE (IWRITE,910) RETURN C C MOST FORMATS, ALSO SOME IN SPOXX C 550 FORMAT (1H ) 560 FORMAT(24H1PORT TESTER,DSY** DBP**) 570 FORMAT ( / 14H EPSILON =, 1PD15.5) 580 FORMAT ( / 16H MAYBE SINGULAR. /) 600 FORMAT (14H ACTUAL COND =, 1PD15.5) 610 FORMAT ( / 4H X =) 620 FORMAT (14H ERROR NORM =, 2(1PD15.5)) 630 FORMAT (14H RESID NORM =, 1P1D15.5) 640 FORMAT (26H TEST RATIOS.. E = EPSILON) 650 FORMAT ( / 14H ************* /) 660 FORMAT (8H1SUMMARY) 670 FORMAT (18H NUMBER OF TESTS =, I4) 680 FORMAT ( / 30H NUMBER OF SUSPICIOUS RATIOS =, 10I4) 690 FORMAT( 42H COND COND(B) ACTUAL ACTUAL ERROR , 1 40HERROR(B) RESID A-RT*R A*AI-I A*AI-I(B)) 700 FORMAT (10(8H -----)) 710 FORMAT(42H ACTUAL ACTUAL COND COND(B) E*COND*X, 1 40H E*COND*X E*A*X E*A E*COND E*COND ) 720 FORMAT (14H NORM(A) =, 1PD15.5) 730 FORMAT (14H NORM(A-RT*R)=, 1PD15.5) 740 FORMAT (14H NORM(A*AI-I)=, 2(1PD15.5)) 750 FORMAT (10(1X, F7.2)) 760 FORMAT(1H ,3D15.5) 780 FORMAT (2D18.6) 790 FORMAT (2D18.6) 830 FORMAT ( / 28H BAND ROUTINES DO NOT AGREE,) 840 FORMAT (5H M =, I2) 860 FORMAT (30H NOT POSITIVE DEFINITE, INFO =, I2) 910 FORMAT ( / 12H END OF TEST) 930 FORMAT (8H RCOND =, 3(1PD15.5)) 980 FORMAT (21H NO SUSPICIOUS RATIOS) 990 FORMAT (I8, 5I10 / 7X, 28H1 INDICATES SUSPICIOUS RATIO) 1000 FORMAT (29H THIS VERSION DATED 09/21/78.) END SUBROUTINE SPOXX(A,LDA,N,KASE,IWRITE) C C GENERATES DOUBLE PRECISION POSITIVE DEFINITE TEST MATRICES C C EXTERNAL D1MACH C FORTRAN DABS,FLOAT,IABS,MAX0,MIN0 INTEGER LDA,N,KASE,IWRITE INTEGER I,J DOUBLE PRECISION A(LDA,1) DOUBLE PRECISION T DOUBLE PRECISION TINY,HUGE,D1MACH C GO TO (10, 10, 10, 50, 50, 70, 70, 70, 120, 160, 200, 240, 290, * 340), KASE C C KASE 1, 2 AND 3 C 10 CONTINUE N = 5*KASE WRITE (IWRITE,20) KASE,N 20 FORMAT (5H KASE, I3, 3X, 16HHILBERT / 4H N =, I4) T = 1.0D0 T = DSIGN(1.0D0,T) DO 40 J = 1, N DO 30 I = 1, N A(I,J) = T**(I - J)/DBLE(FLOAT(I+J-1)) C FOR DOUBLE PRECISION MATRICES, A(I,J) = 1.0/FLOAT(I+J-1) 30 CONTINUE 40 CONTINUE GO TO 350 C C KASE 4 AND 5 C 50 CONTINUE N = 1 WRITE (IWRITE,60) KASE,N 60 FORMAT (5H KASE, I3, 3X, 16HMONOELEMENTAL / 4H N =, I4) IF (KASE .EQ. 4) A(1,1) = 3.0D0 IF (KASE .EQ. 5) A(1,1) = 0.0D0 GO TO 350 C C KASE 6, 7 AND 8 C 70 CONTINUE N = 15 IF (KASE .NE. 8) WRITE (IWRITE,80) KASE,N 80 FORMAT (5H KASE, I3, 3X, 16HTRIDIAGONAL / 4H N =, I4) IF (KASE .EQ. 8) WRITE (IWRITE,90) KASE,N 90 FORMAT (5H KASE, I3, 3X, 16HDIAGONAL / 4H N =, I4) T = 1.0D0 IF (KASE .EQ. 7) T = 2.0D0 IF (KASE .EQ. 8) T = 0.0D0 DO 110 J = 1, N DO 100 I = 1, J A(I,J) = 0.0D0 IF (I .EQ. J) A(I,I) = 4.0D0 IF (I .EQ. J - 1) A(I,J) = T A(J,I) = A(I,J) 100 CONTINUE 110 CONTINUE GO TO 350 C C KASE 9 C 120 CONTINUE N = 5 WRITE (IWRITE,130) KASE,N 130 FORMAT (5H KASE, I3, 3X, 16HPENTADIAGONAL / 4H N =, I4) DO 150 J = 1, N DO 140 I = 1, N A(I,J) = 0.0D0 IF (IABS(I-J) .LE. 2) * A(I,J) = (5.0D0 - DBLE(FLOAT(IABS(I-J))))**(10 - I - J) 140 CONTINUE 150 CONTINUE GO TO 350 C C KASE 10 C 160 CONTINUE N = 6 WRITE (IWRITE,170) KASE,N 170 FORMAT (5H KASE, I3, 3X, 16HTRIDIAG INVERSE / 4H N =, I4) DO 190 J = 1, N DO 180 I = 1, J A(I,J) = DBLE(FLOAT(N+1-J)) A(J,I) = A(I,J) 180 CONTINUE 190 CONTINUE GO TO 350 C C KASE 11 C 200 CONTINUE N = 15 WRITE (IWRITE,210) KASE,N 210 FORMAT (5H KASE, I3, 3X, 16HTEST COND / 4H N =, I4) DO 230 J = 1, N DO 220 I = 1, N IF (I .EQ. J) A(I,J) = DBLE(FLOAT(I)) IF (I .GT. J) A(I,J) = DBLE(FLOAT(J-2)) IF (I .LT. J) A(I,J) = DBLE(FLOAT(I-2)) 220 CONTINUE 230 CONTINUE GO TO 350 C C KASE 12 C 240 CONTINUE N = 5 WRITE (IWRITE,250) KASE,N 250 FORMAT (5H KASE, I3, 3X, 16HNEAR UNDERFLOW / 4H N =, I4) TINY = D1MACH(1)*DBLE(FLOAT(N*N*100)) WRITE (IWRITE,260) TINY 260 FORMAT (14H TINY =, 1PD15.5) DO 280 I = 1, N DO 270 J = 1, N A(I,J) = TINY*DBLE(FLOAT(MIN0(I,J)))/DBLE(FLOAT(MAX0(I,J))) 270 CONTINUE 280 CONTINUE GO TO 350 C C KASE 13 C 290 CONTINUE N = 5 WRITE (IWRITE,300) KASE,N 300 FORMAT (5H KASE, I3, 3X, 16HNEAR OVERFLOW / 4H N =, I4) HUGE = D1MACH(2)/DBLE(FLOAT(N*N)) WRITE (IWRITE,310) HUGE 310 FORMAT (14H HUGE =, 1PD15.5) DO 330 I = 1, N DO 320 J = 1, N A(I,J) = HUGE*DBLE(FLOAT(MIN0(I,J)))/DBLE(FLOAT(MAX0(I,J))) 320 CONTINUE 330 CONTINUE GO TO 350 C 340 CONTINUE N = 0 350 CONTINUE RETURN C END