C$TEST PRAD
C TO RUN AS A MAIN PROGRAM REMOVE NEXT LINE
      SUBROUTINE PRAD
C***********************************************************************
C
C  TEST OF THE PORT DOUBLE PRECISION SPARSE MATRIX PACKAGE
C
C***********************************************************************
C  THIS IS THE TESTER FOR LINDA KAUFMAN'S SPARSE
C  MATRIX PACKAGE - DOUBLE PRECISION
C
C  MARCH 20, 1981
C
C     MAIN PROGRAM
      INTEGER IWRITE, I1MACH
C     ALLOW 5000 UNDERFLOWS.
C
C     SET OUTPUT UNIT NUMBER
C
      IWRITE = I1MACH(2)
C
      CALL SGETS(IWRITE)
      STOP
      END
      SUBROUTINE SGETS(IWRITE)
C     IWRITE IS THE OUTPUT UNIT NUMBER
C
C     TESTS
C        GECE,GEFS,GEBS,GEML,SPMOR,SPFOR,SPMCE,SPFCE,SPMLE,SPFLE,SPMSL
C        DSPFSL,DSPMNF,DSPFNF,SPFSF,DSPMSF,DSPMIN
C
C
C     SUBROUTINES AND FUNCTIONS
C
C     PORT GECE,GEFS,GEBS,GEML,SPMOR,SPFOR,SPMCE,SPFCE,SPMLE,SPFLE
C     PORT DSPMSL,DSPFSL,DSPMNF,DSPFNF,DSPMS,SPFSF.DSPMIN
C     PORT UTILITIES ERROFF,ENTER,LEAVE,R1MACH
C     EXTERNAL SGEXX
C     BLAS DAXPY,DDOT,SSCAL,DASUM
C     FORTRAN ABS,DMAX1,DFLOAT,MAX0,MIN0
C
C     INTERNAL VARIABLES
C
      INTEGER IIA(101), JJA(2500)
      INTEGER JSPSAV(2500),IWORK(2800),IL(101)
      INTEGER MR(100),MC(100),FR(100)
      INTEGER FC(100),MIC(100),JA(2500),IA(101)
      INTEGER I,IPVT(100),IPVTB(100),IQ(8),I1,I2,J
      INTEGER K,KASE,KB,KBFAIL,KOUNT,KP1,KSING,KSUSP(8)
      INTEGER L,LDA,LDAB,IWRITE,M,ML,MU,N,NM1,NPRINT
      DOUBLE PRECISION A(100,100),AS(2500),AINV(100,100),Z(100)
      DOUBLE PRECISION ASAVE(100,100),UL(2500)
      DOUBLE PRECISION XEXACT(100,2),B(100,2),BM(100,2)
      DOUBLE PRECISION BM2(100,2),BF(100,2), D1MACH, DFLOAT
      DOUBLE PRECISION BF2(100,2)
      DOUBLE PRECISION AFNORM,AMNORM,ANRM,CONDF1,CONDF,CONDM1
      DOUBLE PRECISION CONDM,GROWTH,AA(2500)
      DOUBLE PRECISION SPSAV(2500)
      DOUBLE PRECISION AINORM,ANORM,SMACH,COND,COND1,EN,ENORM,EPS, EPS1
      DOUBLE PRECISION ETNORM,FNI,FNORM,ONEPX,RCOND,RCONDB,RNORM
      DOUBLE PRECISION RTNORM,Q(10),QS(10),DASUM,XNORM
      DOUBLE PRECISION EB(2),XN(2),EBM(2),EBF(2),EBM2(2),EBF2(2)
      LOGICAL KBF
C
      COMMON /SPA/ IIA,JJA,AA
      COMMON /CSTAK/ D(12000)
C
      EXTERNAL IIROW,AROW
C
      CALL ISTKIN(12000,3)
C
      LDA = 100
      IAMAX = 2500
      MAXUL = 2500
      IWMAX = 2500
      CALL ENTER(1)
C
C     WRITE MATRIX AND SOLUTIONS IF  N .LE. NPRINT
C
      NPRINT =3
C
      WRITE (IWRITE,480)
      WRITE (IWRITE,770)
C
      DO 10 I = 1, 8
         KSUSP(I) = 0
 10   CONTINUE
      KSING = 0
      KBFAIL = 0
C
C     SET EPS TO ROUNDING UNIT
C
      EPS =D1MACH(4)
      WRITE (IWRITE,490) EPS
      WRITE (IWRITE,470)
C
C     START MAIN LOOP
C
      KASE = 1
 20   CONTINUE
C
C        GENERATE TEST MATRIX
C
         CALL SGEXX(A,LDA,N,KASE,IWRITE)
C
C        N = 0 SIGNALS NO MORE TEST MATRICES
C
C     ...EXIT
         IF (N .LE. 0) GO TO 460
         ANORM = 0.0D0
         DO 30 J = 1, N
            ANORM = DMAX1(ANORM,DASUM(N,A(1,J),1))
 30      CONTINUE
         WRITE (IWRITE,650) ANORM
C
         IF (N .GT. NPRINT) GO TO 50
            WRITE (IWRITE,470)
            DO 40 I = 1, N
               WRITE (IWRITE,690) (A(I,J), J = 1, N)
 40         CONTINUE
            WRITE (IWRITE,470)
 50      CONTINUE
C
C        GENERATE EXACT SOLUTION
C
         XEXACT(1,1) = 1.0D0
         XEXACT(1,2) = 1.0D0
         XEXACT(2,2) = 1.0D0
         IF (N .GE. 2) XEXACT(2,1) = 0.0D0
         IF (N .LE. 2) GO TO 70
            DO 60 I = 3, N
               XEXACT(I,1) = -XEXACT(I-2,1)
                XEXACT(I,2) = 1.0D0
 60         CONTINUE
 70      CONTINUE
C
C        SAVE MATRIX AND GENERATE R.H.S.
C
         DO 90 I = 1, N
            DO 80 J = 1, N
               ASAVE(I,J) = A(I,J)
 80         CONTINUE
 90      CONTINUE
         CALL DGEML(N,A,LDA,XEXACT,B)
         CALL DGEML(N,A,LDA,XEXACT(1,2),B(1,2))
         NN=LDA*2
         CALL MOVEFD(NN,B,BM)
         CALL MOVEFD(NN,B,BF)
         CALL MOVEFD(NN,B,BM2)
         CALL MOVEFD(NN,B,BF2)
      CALL DGETSP(N,A, LDA, IA,JA,AS,MAXA)
      CALL DGETSP(N,A,LDA,IIA,JJA,AA,MAXA)
         DO 100 I=1,MAXA
            SPSAV(I)=AS(I)
            JSPSAV(I)=JA(I)
 100     CONTINUE
         NP1=N+1
C
C CHECK ORDERING SUBROUTINES
C
         CALL SPMOR(N,IA,JA,MC,MIC)
         CALL SPFOR(N,IIROW,FC)
         IBAD=0
         DO 120 I=1,N
            MR(I)=MC(I)
            FR(I)=FC(I)
            IF (MC(I).EQ.FC(I)) GO TO 120
            IBAD=IBAD+1
            WRITE(IWRITE,110)I,MC(I),FC(I)
 110         FORMAT(19H ORDER DISAGREEMENT,3I5)
 120    CONTINUE
        IF (IBAD.EQ.0)WRITE(IWRITE,130)
 130     FORMAT(28H NO DISAGREEMENT IN ORDERING)
C
C        FACTOR AND ESTIMATE CONDITION
C
         CALL DGECE(N,A,LDA,IPVT,RCOND)
         IF (NERROR(INE).NE.0) CALL  ERROFF
C
C
C        FACTOR SPARSE FORM AND COMPARE
C
        CALL DSPMCE(N,MR,MC,AS,IA,JA,IAMAX,IL,ISIZE,CONDM,Z)
        IF(NERROR(INM).NE.0) CALL ERROFF
         CALL DSPFCE(N,FR,FC,AROW,IWORK,UL,MAXUL,ISIZE,CONDF,Z)
         IF (NERROR(INF).NE.0) CALL ERROFF
        WRITE(IWRITE,140)
 140     FORMAT(21H CONDITION ESTIMATION)
        WRITE(IWRITE,150)
 150     FORMAT(42H    GENERAL        SPARSE M       SPARSE F)
        WRITE(IWRITE,160)RCOND,CONDM,CONDF
 160    FORMAT(1H ,3D15.5)
C
C
C           TEST FOR SINGULARITY
C
          IF (INE+INM+INF.EQ.0)GO TO 190
            WRITE(IWRITE,170)INE,INM,INF
 170         FORMAT(13H INE, INM,INF,3I5)
 180          FORMAT(25H SPARSE ROUTINES DISAGREE,2I4)
               WRITE (IWRITE,500)
               KSING = KSING + 1
            GO TO 440
 190        CONTINUE
C
C              COMPUTE INVERSE AND COND1 = TRUE CONDITION
C
               DO 210 J = 1, N
                  DO 200 I = 1, N
                     AINV(I,J) = 0.D0
 200              CONTINUE
               AINV(J,J)=1.D0
 210           CONTINUE
               CALL DGEFS(N,A,LDA,AINV,LDA,N,IPVT)
               CALL DGEBS(N,A,LDA,AINV,LDA,N)
               AINORM=ANRM(N,ASAVE,LDA,AINV,LDA, ANORM, COND1)
               WRITE (IWRITE,520) COND1
               DO 230 I=1,N
                  DO 220 J=1,N
                     AINV(I,J)=0.0D0
 220               CONTINUE
                  AINV(I,I)=1.0D0
 230           CONTINUE
               CALL DSPMSL(N,MR,MC,IA,JA,AS,IL,AINV,LDA,N)
               AMNORM=ANRM(N,ASAVE,LDA,AINV,LDA,ANORM,CONDM1)
               WRITE(IWRITE,260)CONDM1
               DO 250 I=1,N
                   DO 240 J=1,N
                      AINV(I,J)=0.0D0
 240               CONTINUE
                   AINV(I,I)=1.0D0
 250           CONTINUE
               CALL DSPFSL(N,FR,FC,IWORK,UL,AINV,LDA,N)
               AFNORM=ANRM(N,ASAVE,LDA,AINV,LDA,ANORM,CONDF1)
               WRITE(IWRITE,270)CONDF1
 260           FORMAT(19H COND FROM AINV*A M,D15.5)
 270           FORMAT(19H COND FROM AINV*A F,D15.5)
C
C              SOLVE  A*X = B
C
               CALL DGEFS(N,A,LDA,B,LDA,2,IPVT)
               CALL DGEBS(N,A,LDA,B,LDA,2)
C
C              MORE SPAR COMPARE
C
               DO 280 I=1,N
                  MC(I)=MR(I)
                  FC(I)=FR(I)
 280           CONTINUE
               CALL SPMSF(N,MR,MIC,IA,JSPSAV,IWORK,IWMAX,IFILL)
               DO 290 J=1,N
                  I=MR(J)
                  NUM=IA(I+1)-IA(I)
                 IB=IA(I)
                   CALL DSPMIN(N,MIC,IWORK,J,SPSAV(IB),JSPSAV(IB),NUM,
     1          J,UL)
 290           CONTINUE
              EPS1=EPS*ANORM
               CALL DSPMNF(N,IWORK,UL,EPS1,GROWTH)
               IF (NERROR(IMNF).EQ.0) GO TO 300
                 CALL ERROFF
                 GO TO 310
 300           CONTINUE
               CALL DSPSOL(N,MR,MC,IWORK,UL,BM,LDA,2)
 310           CONTINUE
               CALL SPFSF(N,FR,FC,IIROW,IWORK,IWMAX,IFILL)
               CALL DSPFNF(N,FR,FC,AROW,IWORK,UL,GROWTH,EPS1)
               IF (NERROR(IFNF).EQ.0) GO TO 320
                  CALL ERROFF
                  GO TO 330
 320           CONTINUE
               CALL DSPSOL(N,FR,FC,IWORK,UL,BF,LDA,2)
 330           CONTINUE
               CALL DSPMLE(N,.TRUE.,IA,JSPSAV,SPSAV,ISIZE,BM2,LDA,2)
               IF(NERROR(IMLE).NE.0)CALL ERROFF
               CALL DSPFLE(N,.TRUE.,AROW,ISIZE,BF2,LDA,2)
                IF (NERROR(IFLE).NE.0) CALL ERROFF
C COMPUTE ERRORS IN SOLUTION WITH 2 RIGHT HAND SIDES
               IF(IFLE+IMLE+IMNF+IFNF.EQ.0)GO TO 350
                 WRITE(IWRITE,340)IMNF,IFNF,IMLE,IFLE
 340    FORMAT(42H NUMERICALLY SINGULAR- IMNF,IFNF,IMLE,IFLE,4I5)
                GO TO 440
 350           CONTINUE
               DO 360 I=1,2
                  XN(I)=DASUM(N,B(1,I),1)
                  CALL DAXPY(N,-1.0D0,XEXACT(1,I),1,B(1,I),1)
                  CALL DAXPY(N,-1.0D0,XEXACT(1,I),1,BM(1,I),1)
                  CALL DAXPY(N,-1.0D0,XEXACT(1,I),1,BF(1,I),1)
                  CALL DAXPY(N,-1.0D0,XEXACT(1,I),1,BM2(1,I),1)
                  CALL DAXPY(N,-1.0D0,XEXACT(1,I),1,BF2(1,I),1)
                  EB(I)=DASUM(N,B(1,I),1)
                  EBM(I)=DASUM(N,BM(1,I),1)
                  EBF(I)=DASUM(N,BF(1,I),1)
                  EBM2(I)=DASUM(N,BM2(1,I),1)
                  EBF2(I)=DASUM(N,BF2(1,I),1)
 360              CONTINUE
                  WRITE(IWRITE,370)
 370              FORMAT(29H ERRORS IN SOLUTION -ABSOLUTE)
                  WRITE(IWRITE,380)
 380           FORMAT(15X,43H GENERAL SPMN      SPFN       SPMLE   SPFLE
     *)
                  DO 400 I=1,2
                 WRITE(IWRITE,390)I,EB(I),EBM(I),EBF(I),EBM2(I),EBF2(I)
 390      FORMAT(8H PROBLEM,I3,5D15.5)
 400              CONTINUE
C
C
               Q(1) = RCOND/COND1
               Q(2) = CONDM/COND1
               Q(3) = COND1/RCOND
               Q(4)=COND1/CONDM
               Q(5) = EB(1)/(EPS*RCOND*XN(1))
               Q(6) = EBM(1)/(EPS*RCOND*XN(1))
               Q(7) = EBM2(1)/(EPS*RCOND*XN(1))
               Q(8) = AINORM/(EPS*RCOND)
               Q(9)=AMNORM/(EPS*RCOND)
               WRITE (IWRITE,470)
               WRITE (IWRITE,560)
               WRITE (IWRITE,470)
               WRITE (IWRITE,620)
               WRITE (IWRITE,630)
               WRITE (IWRITE,640)
               WRITE (IWRITE,470)
               WRITE (IWRITE,680) (Q(I), I = 1, 9)
               WRITE (IWRITE,470)
C
C              LOOK FOR SUSPICIOUS RATIOS
C
               QS(1) = 1.0D0 + 4.0D0*EPS
               QS(2)=QS(1)
               QS(4) = 10.0D0
               QS(3)=10.0D0
               EN = DFLOAT(N)
               IF (N .EQ. 1) EN = 2.0D0
               DO 410 I = 3, 10
                  QS(I) = EN
 410           CONTINUE
               KOUNT = 0
               DO 430 I = 1, 9
                  IQ(I) = 0
                  IF (Q(I) .LE. QS(I)) GO TO 420
                     IQ(I) = 1
                     KSUSP(I) = KSUSP(I) + 1
                     KOUNT = KOUNT + 1
 420              CONTINUE
 430           CONTINUE
               IF (KOUNT .EQ. 0) WRITE (IWRITE,750)
               IF (KOUNT .NE. 0) WRITE (IWRITE,760) (IQ(I), I = 1, 9)
               WRITE (IWRITE,470)
 440        CONTINUE
 450     CONTINUE
C
         WRITE (IWRITE,570)
         KASE = KASE + 1
      GO TO 20
 460  CONTINUE
C
C     FINISH MAIN LOOP
C
C     SUMMARY
C
      WRITE (IWRITE,580)
      KASE = KASE - 1
      WRITE (IWRITE,590) KASE
      WRITE (IWRITE,600) KSING
      WRITE (IWRITE,610) KSUSP
      WRITE (IWRITE,730)
      CALL LEAVE
      RETURN
C
C     MOST FORMATS, ALSO SOME IN SGEXX
C
 470  FORMAT (1H )
 480  FORMAT (29H1  PORT  TESTER, DGE**, DSP**)
 490  FORMAT ( / 14H EPSILON     =, 1PD13.5)
 500  FORMAT ( / 19H EXACT SINGULARITY. /)
 510  FORMAT ( / 16H MAYBE SINGULAR. /)
 520  FORMAT (14H ACTUAL COND =, 1PD13.5)
 530  FORMAT(/14H X AND XBAND =)
 540  FORMAT (14H ERROR NORMS =, 2(1PD13.5))
 550  FORMAT (14H RESID NORMS =, 2(1PD13.5))
 560  FORMAT (26H TEST RATIOS.. E = EPSILON)
 570  FORMAT ( / 14H ************* /)
 580  FORMAT (8H1SUMMARY)
 590  FORMAT (18H NUMBER OF TESTS =, I4)
 600  FORMAT (30H NUMBER OF SINGULAR MATRICES =, I4)
 610  FORMAT (30H NUMBER OF SUSPICIOUS RATIOS =, 10I4)
 620  FORMAT(40H    COND  COND(M)  ACTUAL  ACTUAL  ERROR,
     1       36H ERROR(M) ERROR(M2) A*AI-I A*AI-I(M))
 630  FORMAT (9(8H   -----))
 640  FORMAT(42H    ACTUAL ACTUAL   COND  COND(M) E*COND*X,
     1       31H E*CONDX E*COND*X E*COND E*COND)
 650  FORMAT (14H NORM(A)     =, 1PD13.5)
 660  FORMAT (14H NORM(A - LU)=, 1PD13.5)
 670  FORMAT (14H NORM(A*AI-I)=, 2(1PD13.5))
 680  FORMAT (10(1X, F7.2))
 690  FORMAT (1H , 6G11.4)
 700  FORMAT (14H 1/COND      =, 1PD13.5)
 710  FORMAT (2G14.6)
 720  FORMAT (5H ML =, I2, 6H  MU =, I2)
 730  FORMAT ( / 12H END OF TEST)
 740  FORMAT(7H COND =,1PD13.5,13H COND(BAND) =,1PD13.5 /)
 750  FORMAT (21H NO SUSPICIOUS RATIOS)
 760  FORMAT (I8, 9I10 / 7X, 28H1 INDICATES SUSPICIOUS RATIO)
 770  FORMAT (29H THIS VERSION DATED 03/11/78.)
      END
      SUBROUTINE SGEXX(A,LDA,N,KASE,IWRITE)
C
C     GENERATES DOUBLE PRECISION GENERAL TEST MATRICES
C
C     EXTERNAL SMACH
C     FORTRAN DFLOAT,MAX0
      INTEGER I,J,IRAND
      INTEGER LDA,N,KASE,IWRITE
      DOUBLE PRECISION A(LDA,1), D1MACH, DFLOAT
      DOUBLE PRECISION T1,T2
      DOUBLE PRECISION SMACH,HUGE,TINY
C
      GO TO (10, 10, 10, 60, 60, 80, 80, 80, 120, 160, 200, 250, 290,
     *       330, 370, 420, 470, 530, 600), KASE
C
C     KASE 1, 2 AND 3
C
 10   CONTINUE
         N = 3*KASE
         WRITE (IWRITE,20) KASE,N
 20      FORMAT (5H KASE, I3, 3X, 16HHILBERT SLICE    / 4H N =, I4)
         DO 50 J = 1, N
            DO 40 I = 1, N
               A(I,J) = 0.0D0
               IF (I .GT. J + 2) GO TO 30
               IF (I .LT. J - 3) GO TO 30
                  A(I,J) = 1.0D0/DFLOAT(I+J-1)
 30            CONTINUE
 40         CONTINUE
 50      CONTINUE
      GO TO 610
C
C     KASE 4 AND 5
C
 60   CONTINUE
         N = 1
         WRITE (IWRITE,70) KASE,N
 70      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 610
C
C     KASE 6, 7 AND 8
C
 80   CONTINUE
         N = 15
         WRITE (IWRITE,90) KASE,N
 90      FORMAT (5H KASE, I3, 3X, 16HTRIDIAGONAL      / 4H N =, I4)
         T1 = 1.0D0
         T2 = 1.0D0
         IF (KASE .EQ. 7) T1 = 100.0D0
         IF (KASE .EQ. 8) T2 = 100.0D0
         DO 110 I = 1, N
            DO 100 J = 1, N
               A(I,J) = 0.0D0
               IF (I .EQ. J) A(I,I) = 4.0D0
               IF (I .EQ. J - 1) A(I,J) = T1
               IF (I .EQ. J + 1) A(I,J) = T2
 100        CONTINUE
 110     CONTINUE
      GO TO 610
C
C     KASE 9
C
 120  CONTINUE
         N = 5
         WRITE (IWRITE,130) KASE,N
 130     FORMAT (5H KASE, I3, 3X, 16HRANK ONE         / 4H N =, I4)
         DO 150 I = 1, N
            DO 140 J = 1, N
               A(I,J) = 10.0D0**(I - J)
 140        CONTINUE
 150     CONTINUE
      GO TO 610
C
C     KASE 10
C
 160  CONTINUE
         N = 4
         WRITE (IWRITE,170) KASE,N
 170     FORMAT (5H KASE, I3, 3X, 16HZERO COLUMN      / 4H N =, I4)
         DO 190 I = 1, N
            DO 180 J = 1, N
               T1 = DFLOAT(J-3)
               T2 = DFLOAT(I)
               A(I,J) = T1/T2
 180        CONTINUE
 190     CONTINUE
      GO TO 610
C
C     KASE 11
C
 200  CONTINUE
         N = 5
         WRITE (IWRITE,210) KASE,N
 210     FORMAT (5H KASE, I3, 3X, 16HTEST COND        / 4H N =, I4)
         WRITE(IWRITE,220)(A(I,I),I=1,N)
 220       FORMAT(15H DIGAG FROM GEN,5E15.5)
         DO 240 I = 1, N
            DO 230 J = 1, N
               IF (I .EQ. J) A(I,J) = DFLOAT(I)
               IF (I .GT. J) A(I,J) = DFLOAT(J-2)
               IF (I .LT. J) A(I,J) = DFLOAT(I-2)
 230        CONTINUE
 240     CONTINUE
      GO TO 610
C
C     KASE 12
C
 250  CONTINUE
         N = 3
         WRITE (IWRITE,260) KASE,N
 260     FORMAT (5H KASE, I3, 3X, 16HIDENTITY         / 4H N =, I4)
         DO 280 I = 1, N
            DO 270 J = 1, N
               IF (I .EQ. J) A(I,I) = 1.0D0
               IF (I .NE. J) A(I,J) = 0.0D0
 270        CONTINUE
 280     CONTINUE
      GO TO 610
C
C     KASE 13
C
 290  CONTINUE
         N = 6
         WRITE (IWRITE,300) KASE,N
 300     FORMAT (5H KASE, I3, 3X, 16HUPPER TRIANGULAR / 4H N =, I4)
         DO 320 I = 1, N
            DO 310 J = 1, N
               IF (I .GT. J) A(I,J) = 0.0D0
               IF (I .LE. J) A(I,J) = DFLOAT(J-I+1)
 310        CONTINUE
 320     CONTINUE
      GO TO 610
C
C     KASE 14
C
 330  CONTINUE
         N = 6
         WRITE (IWRITE,340) KASE,N
 340     FORMAT (5H KASE, I3, 3X, 16HLOWER TRIANGULAR / 4H N =, I4)
         DO 360 I = 1, N
            DO 350 J = 1, N
               IF (I .LT. J) A(I,J) = 0.0D0
               IF (I .GE. J) A(I,J) = DFLOAT(I-J+1)
 350        CONTINUE
 360     CONTINUE
      GO TO 610
C
C     KASE 15
C
 370  CONTINUE
         N = 5
         WRITE (IWRITE,380) KASE,N
 380     FORMAT (5H KASE, I3, 3X, 16HNEAR UNDERFLOW   / 4H N =, I4)
         TINY =D1MACH(1)*DFLOAT(N*N)*1000.0
         WRITE (IWRITE,390) TINY
 390     FORMAT (14H TINY        =, 1PD13.5)
         DO 410 I = 1, N
            DO 400 J = 1, N
               A(I,J) = TINY*DFLOAT(J)/DFLOAT(MAX0(I,J))
 400        CONTINUE
 410     CONTINUE
      GO TO 610
C
C     KASE 16
C
 420  CONTINUE
         N = 5
         WRITE (IWRITE,430) KASE,N
 430     FORMAT (5H KASE, I3, 3X, 16HNEAR OVERFLOW    / 4H N =, I4)
         HUGE =D1MACH(2)/DFLOAT(12000*N*N)
         WRITE (IWRITE,440) HUGE
 440     FORMAT (14H HUGE        =, 1PD13.5)
         DO 460 I = 1, N
            DO 450 J = 1, N
               A(I,J) = HUGE*(DFLOAT(J)/DFLOAT(MAX0(I,J)))
 450        CONTINUE
 460     CONTINUE
      GO TO 610
C
 470  CONTINUE
C
C THIS IS RANDOM SPARSE ROUTINE
C
        WRITE(IWRITE,480)KASE,N
 480    FORMAT(5H KASE, I3,14H RANDOM SPARSE,/4H N =, I4)
        DO 500 I=1,N
           DO 490 J=1,N
              A(I,J)=0.0D0
 490       CONTINUE
 500     CONTINUE
        K=N/2
         DO 520 I=1,N
            J=0
            IRAND = UNI(0)
 510        JJ=J+K*IRAND+1
             IF (J.LT.I.AND.JJ.GT.I) JJ=I
             J=JJ
             IF (J.GT.N) GO TO 520
               A(I,J)=UNI(0)
               GO TO 510
 520     CONTINUE
         GO TO 610
C
C QUEING PROBLEM
C
 530    CONTINUE
        N=100
        WRITE(IWRITE,540)KASE,N
 540    FORMAT(5H KASE,I3, 9H QUEUEING,/4H N =, I4)
        DO 560 I=1,N
           DO 550 J=1,N
              A(I,J)=0.0D0
 550       CONTINUE
 560    CONTINUE
        M=SQRT(FLOAT(N))
        L=0
       A1=M
        DO 580 I=1,M
           DO 570 J=1,M
              L=L+1
               A(L,L)=-2.0*A1-DFLOAT(I+J-2)
              IF (J.GT.1)A(L,L-1)=A1
              IF (J.LT.M)A(L,L+1)=J
              IF (I.GT.1)A(L,L-M)=A1
              IF (I.LT.M)A(L,L+M)=I
 570        CONTINUE
 580    CONTINUE
        DO 590 I=1,N
 590         A(N,I)=1.0
          GO TO 610
 600       CONTINUE
         N = 0
 610  CONTINUE
      RETURN
C
      END
      SUBROUTINE AROW(I, ROW, JCOL, NUM)
      INTEGER JCOL(200), IA(101), JA(2500), JPIB
      DOUBLE PRECISION ROW(200),A(2500)
      COMMON /SPA/IA,JA,A
       NUM=IA(I+1)-IA(I)
      IB=IA(I)-1
      DO 10 J=1,NUM
          JPIB = J+IB
          ROW(J)=A(JPIB)
          JCOL(J)=JA(JPIB)
 10   CONTINUE
      RETURN
      END
      SUBROUTINE IIROW(I,  JCOL, NUM)
      INTEGER JCOL (200), IA(101), JA(2500), JPIB
      DOUBLE PRECISION A(2500)
      COMMON /SPA/IA,JA,A
      NUM=IA(I+1)-IA(I)
      IB=IA(I)-1
      DO 10 J=1,NUM
          JPIB = J+IB
          JCOL(J)=JA(JPIB)
 10   CONTINUE
      RETURN
      END
      DOUBLE PRECISION FUNCTION ANRM(N,AS,LDA1,AINV,LDA,ANORM,COND1)
      DOUBLE PRECISION AS(LDA,N),AINV(LDA1,N)
      DOUBLE PRECISION AIN,ANORM,COND1,DASUM,T
C
C     THIS SUBROUTINE COMPUTES THE CONDITION NUMBER FROM
C     A AND AINVERSE
C
C     IT ALSO COMPUTES THE ERROR A*AINV-I AND STORES IT IN ANRM
      DOUBLE PRECISION B(200)
      ANRM=0.0D0
      DO 30 J =1,N
         DO 10 I=1,N
            B(I)=0.0D0
 10      CONTINUE
         DO 20 K=1,N
            T=AINV(K,J)
            CALL DAXPY(N,T,AS(1,K),1,B,1)
 20       CONTINUE
         B(J)=B(J)-1.0D0
         ANRM=DMAX1(ANRM,DASUM(N,B,1))
 30   CONTINUE
      AIN=0.0D0
      DO 40 I=1,N
         AIN=DMAX1(AIN,DASUM(N,AINV(1,I),1))
 40   CONTINUE
      COND1=ANORM*AIN
      RETURN
      END
         SUBROUTINE DGETSP(N,A,LDA,IA,JA,AS,L)
         DOUBLE PRECISION A(LDA,N),AS(1)
         INTEGER JA(1),IA(1)
         L=0
         IA(1)=1
         DO 30 I=1,N
            DO 10 J=1,N
               IF (A(I,J).EQ.0.0D0) GO TO 10
               L=L+1
               AS(L)=A(I,J)
               JA(L)=J
 10         CONTINUE
           IF (L.GE.IA(I))GO TO 20
               L=L+1
               AS(L)=0.0
               JA(L)=1
 20        CONTINUE
            IA(I+1)=L+1
 30     CONTINUE
        RETURN
        END

.