home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / eispack-1.0-src.tgz / tar.out / contrib / eispack / ex / rsptest.f < prev    next >
Text File  |  1996-09-28  |  15KB  |  434 lines
  1. C
  2. C     THIS DRIVER TESTS  EISPACK  FOR THE CLASS OF REAL SYMMETRIC
  3. C     PACKED MATRICES SUMMARIZING THE FIGURES OF MERIT FOR ALL PATHS.
  4. C
  5. C     THIS DRIVER IS CATALOGUED AS  EISPDRV4(RSPSUMAR).
  6. C
  7. C     THE DIMENSION OF  A  SHOULD BE  NNN  AND THE DIMENSION
  8. C     OF  Z  SHOULD BE  NM  BY  NM.
  9. C     THE DIMENSION OF  W,D,E,E2,IND,RV1,RV2,RV3,RV4,RV5,RV6,
  10. C     W1, AND  W2  SHOULD BE  NM.
  11. C     THE DIMENSION OF  AHOLD  SHOULD BE  NNN.
  12. C     HERE NM = 20, AND  NNN = NM*(NM+1)/2.
  13. C
  14.       DOUBLE PRECISION A(210),Z(20,20),AHOLD(210),W(20),D(20),E(20),
  15.      X        E2( 20),RV1( 20),RV2( 20),RV3( 20),RV4( 20),RV5( 20),
  16.      X        RV6( 20),W1( 20),W2( 20),TCRIT( 8),EPSLON,RESDUL,MAXEIG,
  17.      X        MAXDIF,U,LB,UB,EPS1,DFL
  18.       REAL  XUB,XLB
  19.       INTEGER  IND( 20),IERR( 6),ERROR
  20.       DATA IREAD1/1/,IREADC/5/,IWRITE/6/
  21. C
  22.       OPEN(UNIT=IREAD1,FILE='FILE35')
  23.       OPEN(UNIT=IREADC,FILE='FILE36')
  24.       REWIND IREAD1
  25.       REWIND IREADC
  26. C
  27.       NM = 20
  28.       NNN = (NM*(NM+1))/2
  29.       LCOUNT = 0
  30.       WRITE(IWRITE,1)
  31.     1 FORMAT(1H1,19X,57H EXPLANATION OF COLUMN ENTRIES FOR THE SUMMARY S
  32.      XTATISTICS//1H ,95(1H-)/96H ORDER TQL2   TQLRAT IMTQL2 IMTQL1    LB
  33.      X      UB    M  IMTQLV   TSTURM   BISECT  M1 NO  TRIDIB  /1H ,
  34.      X95(1H-)//48H UNDER 'ORDER' IS THE ORDER OF EACH TEST MATRIX. //
  35.      X95H UNDER 'TQL2   TQLRAT' ARE THREE NUMBERS.  THE FIRST NUMBER, AN
  36.      X INTEGER, IS THE ABSOLUTE SUM OF/
  37.      X61H THE ERROR FLAGS RETURNED SEPARATELY FROM  TQL2  AND  TQLRAT.,
  38.      X34H  THE SECOND NUMBER IS THE MEASURE/
  39.      X62H OF PERFORMANCE BASED UPON THE RESIDUAL COMPUTED FOR THE  TQL2,
  40.      X25H  PATH.  THE THIRD NUMBER        /
  41.      X62H MEASURES THE AGREEMENT OF THE EIGENVALUES FROM THE  TQL2  AND,
  42.      X16H  TQLRAT  PATHS.  //
  43.      X95H UNDER 'IMTQL2 IMTQL1' ARE THREE NUMBERS WITH MEANING LIKE THOS
  44.      XE UNDER  'TQL2   TQLRAT'.       //
  45.      X95H UNDER 'LB' AND 'UB' ARE THE INPUT VARIABLES SPECIFYING THE INT
  46.      XERVAL TO  BISECT  AND  TSTURM.  //
  47.      X61H UNDER 'M' IS THE NUMBER OF EIGENVALUES DETERMINED BY  BISECT,
  48.      X30H  AND  TSTURM  THAT LIE IN THE    /18H INTERVAL (LB,UB).//
  49.      X95H UNDER EACH OF 'IMTQLV', 'TSTURM', 'BISECT', AND 'TRIDIB' ARE T
  50.      XWO NUMBERS.  THE FIRST NUMBER,       )
  51.       WRITE(IWRITE,2)
  52.     2 FORMAT(
  53.      X95H AN INTEGER, IS THE ABSOLUTE SUM OF THE ERROR FLAGS RETURNED FR
  54.      XOM THE RESPECTIVE PATH.         /
  55.      X95H THE SECOND NUMBER IS THE MEASURE OF PERFORMANCE BASED UPON THE
  56.      X RESIDUAL COMPUTED FOR THE PATH.//
  57.      X95H UNDER 'M1' AND 'NO' ARE THE VARIABLES SPECIFYING THE LOWER BOU
  58.      XNDARY INDEX AND THE NUMBER        /
  59.      X28H OF EIGENVALUES TO  TRIDIB.   //
  60.      X62H -1.0  AS THE MEASURE OF PERFORMANCE IS PRINTED IF AN ERROR IN,
  61.      X27H THE CORRESPONDING PATH HAS        /
  62.      X47H PREVENTED THE COMPUTATION OF THE EIGENVECTORS. //
  63.      X64H THE  TQL2    PATH USES THE EISPACK CODES   TRED3-TQL2  -TRBAK3
  64.      X,  /
  65.      X39H AS CALLED FROM DRIVER SUBROUTINE  RSP. /
  66.      X62H THE  TQLRAT  PATH USES THE EISPACK CODES   TRED3-TQLRAT,     /
  67.      X39H AS CALLED FROM DRIVER SUBROUTINE  RSP. /
  68.      X64H THE  IMTQL2  PATH USES THE EISPACK CODES   TRED3-IMTQL2-TRBAK3
  69.      X.  )
  70.       WRITE(IWRITE,3)
  71.     3 FORMAT(
  72.      X62H THE  IMTQL1  PATH USES THE EISPACK CODES   TRED3-IMTQL1.     /
  73.      X63H THE  IMTQLV  PATH USES THE EISPACK CODES   TRED3-IMTQLV-TINVIT
  74.      X ,8H-TRBAK3./
  75.      X64H THE  TSTURM  PATH USES THE EISPACK CODES   TRED3-TSTURM-TRBAK3
  76.      X.  /
  77.      X63H THE  BISECT  PATH USES THE EISPACK CODES   TRED3-BISECT-TINVIT
  78.      X ,8H-TRBAK3. /
  79.      X63H THE  TRIDIB  PATH USES THE EISPACK CODES   TRED3-TRIDIB-TINVIT
  80.      X ,8H-TRBAK3. /)
  81.       WRITE(IWRITE,15)
  82.    15 FORMAT(1X,21HD.P. VERSION 04/15/83 )
  83.     5 FORMAT( 53H1       TABULATION OF THE ERROR FLAG  ERROR  AND THE ,
  84.      X    31HMEASURE OF PERFORMANCE  Y  FOR /5X,
  85.      X    56HTHE  EISPACK  CODES.  THIS RUN DISPLAYS THESE STATISTICS ,
  86.      X    40H FOR REAL SYMMETRIC PACKED MATRICES.        /
  87.      X    55H0ORDER TQL2   TQLRAT IMTQL2 IMTQL1    LB      UB    M   ,
  88.      X    40HIMTQLV   TSTURM   BISECT  M1 NO  TRIDIB )
  89.    10 CALL RMATIN(NNN,N,A,AHOLD,0)
  90.       READ(IREADC,50) MM,LB,UB,M11,NO
  91.    50 FORMAT(I4,2D24.16,2(4X,I4))
  92. C
  93. C     MM,LB,UB,M11,  AND  NO  ARE READ FROM SYSIN AFTER THE MATRIX IS
  94. C     GENERATED.  MM,LB,  AND  UB  SPECIFY TO  BISECT  THE MAXIMUM
  95. C     NUMBER OF EIGENVALUES AND THE BOUNDS FOR THE INTERVAL WHICH IS
  96. C     TO BE SEARCHED.  M11  AND  NO  SPECIFY TO  TRIDIB  THE LOWER
  97. C     BOUNDARY INDEX AND THE NUMBER OF DESIRED EIGENVALUES.
  98. C
  99.       DO  230  ICALL = 1,10
  100.          IF( ICALL .NE. 1 )  CALL  RMATIN(NNN,N,A,AHOLD,1)
  101. C
  102. C     IF  TQLRAT  PATH (LABEL 80) IS TAKEN THEN  TQL2  PATH (LABEL 70)
  103. C     MUST ALSO BE TAKEN IN ORDER THAT THE MEASURE OF PERFORMANCE BE
  104. C     MEANINGFUL.
  105. C     IF  IMTQL1  PATH (LABEL 85) IS TAKEN THEN  IMTQL2  PATH (LABEL 75)
  106. C     MUST ALSO BE TAKEN IN ORDER THAT THE MEASURE OF PERFORMANCE BE
  107. C     MEANINGFUL.
  108. C     IF  TQL2  (IMTQL2)  PATH FAILS, THEN  TQLRAT  (IMTQL1)  PATH IS
  109. C     OMITTED AND PRINTOUT FLAGGED WITH  -1.0.
  110. C
  111.          GO TO  (70,75,80,85,89,90,95,230,110,230),  ICALL
  112. C
  113. C     RSPWZ  USING  TQL2
  114. C     INVOKED FROM DRIVER SUBROUTINE  RSP.
  115. C
  116.    70    ICT = 1
  117.          CALL  RSP(NM,N,NNN,A,W,1,Z,E,E,ERROR)
  118.          IERR(ICT) = ERROR
  119.          M = ERROR - 1
  120.          IF( ERROR .NE. 0 ) GO TO 74
  121.          M = N
  122.          DO 71 I = 1,N
  123.             W1(I) = W(I)
  124.    71    CONTINUE
  125.    74    GO TO  190
  126. C
  127. C     RSPWZ  USING  IMTQL2
  128. C
  129.    75    ICT = 2
  130.          DO 77 I=1,N
  131.            DO 76 J=1,N
  132.    76        Z(I,J)=0.0D0
  133.    77      Z(I,I)=1.0D0
  134.          CALL  TRED3(N,NNN,A,W,E,E)
  135.          CALL  IMTQL2(NM,N,W,E,Z,ERROR)
  136.          IERR(ICT) = ERROR
  137.          M = ERROR - 1
  138.          IF( ERROR .NE. 0 ) GO TO 79
  139.          DO 78 I = 1,N
  140.    78       W2(I) = W(I)
  141.          M = N
  142.    79    CALL  TRBAK3(NM,N,NNN,A,M,Z)
  143.          GO TO  190
  144. C
  145. C     RSPW  USING  TQLRAT
  146. C     INVOKED FROM DRIVER SUBROUTINE  RSP.
  147. C
  148.    80    ICT = 7
  149.          IF( IERR(1) .NE. 0 ) GO TO 200
  150.          CALL  RSP(NM,N,NNN,A,W,0,Z,E,E2,ERROR)
  151.          IERR(1) = ERROR
  152.          IF( ERROR .NE. 0 ) GO TO 200
  153.          MAXEIG = 0.0D0
  154.          MAXDIF = 0.0D0
  155.          DO 81 I = 1,N
  156.             IF( DABS(W(I)) .GT. MAXEIG ) MAXEIG = DABS(W(I))
  157.             U = DABS(W1(I) - W(I))
  158.             IF( U .GT. MAXDIF ) MAXDIF = U
  159.    81    CONTINUE
  160.          IF( MAXEIG .EQ. 0.0D0 ) MAXEIG = 1.0D0
  161.          DFL = 10 * N
  162.          TCRIT(7) = MAXDIF/EPSLON(MAXEIG*DFL)
  163.          GO TO  230
  164. C
  165. C     RSPW  USING  IMTQL1
  166. C
  167.    85    ICT = 8
  168.          IF( IERR(2) .NE. 0 ) GO TO 200
  169.          CALL  TRED3(N,NNN,A,W,E,E)
  170.          CALL  IMTQL1(N,W,E,ERROR)
  171.          IERR(2) = ERROR
  172.          MAXEIG = 0.0D0
  173.          MAXDIF = 0.0D0
  174.          DO 86 I = 1,N
  175.             IF( DABS(W(I)) .GT. MAXEIG ) MAXEIG = DABS(W(I))
  176.             U = DABS(W2(I) - W(I))
  177.             IF( U .GT. MAXDIF ) MAXDIF = U
  178.    86    CONTINUE
  179.          IF( MAXEIG .EQ. 0.0D0 ) MAXEIG = 1.0D0
  180.          DFL = 10 * N
  181.          TCRIT(8) = MAXDIF/EPSLON(MAXEIG*DFL)
  182.          GO TO  230
  183. C
  184. C     RSPW1Z  ( USAGE HERE COMPUTES ALL THE EIGENVECTORS )
  185. C
  186.    89    ICT = 3
  187.          CALL  TRED3(N,NNN,A,D,E,E2)
  188.          CALL  IMTQLV(N,D,E,E2,W,IND,ERROR,RV1)
  189.          IERR(ICT) = ERROR
  190.          M = N
  191.          IF( ERROR .NE. 0 ) M = ERROR - 1
  192.          CALL  TINVIT(NM,N,D,E,E2,M,W,IND,Z,ERROR,RV1,RV2,RV3,RV4,RV6)
  193.          IERR(ICT) = IERR(ICT) + IABS(ERROR)
  194.          CALL  TRBAK3(NM,N,NNN,A,M,Z)
  195.          GO TO 190
  196. C
  197. C     RSP1W1Z  USING  TSTURM
  198. C
  199.    90    ICT = 4
  200.          EPS1 = 0.0D0
  201.          CALL  TRED3(N,NNN,A,D,E,E2)
  202.          CALL  TSTURM(NM,N,EPS1,D,E,E2,LB,UB,MM,M,W,Z,ERROR,
  203.      X                RV1,RV2,RV3,RV4,RV5,RV6)
  204.          IERR(ICT) = ERROR
  205.          XLB = LB
  206.          XUB = UB
  207.          IF( ERROR .EQ. 3*N + 1 ) GO TO 200
  208.          IF( ERROR .GT. 4*N ) M = ERROR - 4*N + 1
  209.          CALL  TRBAK3(NM,N,NNN,A,M,Z)
  210.          GO TO  190
  211. C
  212. C     RSP1W1Z  USING  BISECT  AND  TINVIT
  213. C
  214.    95    ICT = 5
  215.          EPS1 = 0.0D0
  216.          CALL  TRED3(N,NNN,A,D,E,E2)
  217.          CALL  BISECT(N,EPS1,D,E,E2,LB,UB,MM,M,W,IND,ERROR,RV1,RV2)
  218.          IERR(ICT) = ERROR
  219.          MBISCT = M
  220.          XLB = LB
  221.          XUB = UB
  222.          IF( ERROR .NE. 0 ) GO TO 200
  223.          CALL  TINVIT(NM,N,D,E,E2,M,W,IND,Z,ERROR,RV1,RV2,RV3,RV4,RV6)
  224.          IERR(ICT) = IABS(ERROR)
  225.          CALL  TRBAK3(NM,N,NNN,A,M,Z)
  226.          GO TO  190
  227. C
  228. C     RSP1W1Z  USING  TRIDIB  AND  TINVIT
  229. C
  230.   110    ICT = 6
  231.          EPS1 = 0.0D0
  232.          CALL  TRED3(N,NNN,A,D,E,E2)
  233.          CALL  TRIDIB(N,EPS1,D,E,E2,LB,UB,M11,NO,W,IND,ERROR,RV4,RV5)
  234.          IERR(ICT) = ERROR
  235.          IF( ERROR .NE. 0 )  GO TO  200
  236.          M = NO
  237.          CALL  TINVIT(NM,N,D,E,E2,M,W,IND,Z,ERROR,RV1,RV2,RV3,RV4,RV6)
  238.          IERR(ICT) = IABS(ERROR)
  239.          CALL  TRBAK3(NM,N,NNN,A,M,Z)
  240. C
  241.   190    IF( M .EQ. 0 .AND. ERROR .NE. 0 ) GO TO 200
  242.          CALL  RMATIN(NNN,N,A,AHOLD,1)
  243.          CALL  RSPWZR(NM,N,NNN,M,A,W,Z,RV1,RESDUL)
  244.          DFL = 10 * N
  245.          TCRIT(ICT) = RESDUL/EPSLON(DFL)
  246.          GO TO 230
  247.   200    TCRIT(ICT) = -1.0D0
  248.   230    CONTINUE
  249. C
  250.       IF( MOD(LCOUNT,35) .EQ. 0 ) WRITE(IWRITE,5)
  251.       LCOUNT = LCOUNT + 1
  252.       WRITE(IWRITE,240) N,IERR(1),TCRIT(1),TCRIT(7),IERR(2),TCRIT(2),
  253.      X             TCRIT(8),XLB,XUB,MBISCT,(IERR(I),TCRIT(I),I=3,5),
  254.      X             M11,NO,IERR(6),TCRIT(6)
  255.   240 FORMAT(I4,2(I3,2F6.3),2(1PE8.0),I3,3(I3,0PF6.3),3I3,F6.3)
  256.       GO TO  10
  257.       END
  258.       SUBROUTINE RSPWZR(NM,N,NNN,M,A,W,Z,NORM,RESDUL)
  259. C
  260.       DOUBLE PRECISION NORM(M), W(M), A(NNN), Z(NM,M), NORMA, S, SUM,
  261.      X       SUMA, SUMZ, RESDUL, TNORM
  262. C
  263. C     THIS SUBROUTINE FORMS THE 1-NORM OF THE RESIDUAL MATRIX
  264. C     A*Z-Z*DIAG(W)   WHERE  A  IS A REAL SYMMETRIC MATRIX STORED IN
  265. C     PACKED FORM,  W  IS A VECTOR WHICH CONTAINS M EIGENVALUES OF
  266. C     A , AND  Z  IS AN ARRAY WHICH CONTAINS THE CORRESPONDING EIGEN-
  267. C     VECTORS OF  A . ALL NORMS APPEARING IN THE COMMENTS BELOW ARE
  268. C     1-NORMS.
  269. C
  270. C     THIS SUBROUTINE IS CATALOGUED AS EISPDRV4(RSPWZR).
  271. C
  272. C     INPUT.
  273. C
  274. C        NM  IS THE ROW DIMENSION OF TWO-DIMENSIONAL ARRAY PARAMETERS
  275. C           AS DECLARED IN THE CALLING PROGRAM DIMENSION STATEMENT;
  276. C
  277. C        N  IS THE ORDER OF THE MATRIX  A;
  278. C
  279. C        NNN  IS THE DIMENSION OF THE ARRAY PARAMETER  A;
  280. C
  281. C        M  IS THE NUMBER OF EIGENVECTORS FOR WHICH RESIDUALS ARE
  282. C           DESIRED;
  283. C
  284. C        A(N*(N+1)/2) IS A VECTOR WHICH CONTAINS THE ELEMENTS OF THE
  285. C           LOWER TRIANGULAR PART OF THE MATRIX  A  (AS MENTIONED ABOVE)
  286. C           IN PACKED FORM.  BY PACKED FORM, WE MEAN THAT THE FIRST ROW
  287. C           OF THE TRIANGLE IS STORED IN THE FIRST POSITION OF  A , THE
  288. C           SECOND ROW OF THE TRIANGLE IS STORED IN THE NEXT TWO
  289. C           POSITIONS, AND SO FORTH UNTIL WE HAVE THE N-TH ROW STORED
  290. C           IN THE LAST N POSITIONS OF  A;
  291. C
  292. C        W(M)  IS AN ARRAY CONTAINING THE EIGENVALUES OF  A;
  293. C
  294. C        Z(NM,M) IS AN ARRAY WHICH CONTAINS THE EIGENVECTORS OF  A.
  295. C
  296. C
  297. C     OUTPUT.
  298. C
  299. C        Z(NM,M) IS AN ARRAY WHICH CONTAINS THE NORMALIZED
  300. C           APPROXIMATE EIGENVECTORS OF  A.  THE EIGENVECTORS
  301. C           ARE NORMALIZED USING THE 1-NORM IN SUCH A WAY
  302. C           THAT THE FIRST ELEMENT WHOSE MAGNITUDE IS LARGER
  303. C           THAN THE NORM OF THE EIGENVECTOR DIVIDED BY  N  IS
  304. C           POSITIVE;
  305. C
  306. C        NORM(M) IS AN ARRAY SUCH THAT FOR EACH  K
  307. C           NORM(K) = !!A*Z(K)-Z(K)*W(K)!!/(!!A!!*!!Z(K)!!)
  308. C           WHERE  Z(K)  IS THE K-TH EIGENVECTOR;
  309. C
  310. C        RESDUL IS THE REAL NUMBER
  311. C           !!A*Z-Z*DIAG(W)!!/(!!A!!*!!Z!!).
  312. C
  313. C     ----------------------------------------------------------------
  314. C
  315.       RESDUL = 0.0D0
  316.       IF( M .EQ. 0 ) RETURN
  317.       NORMA = 0.0D0
  318.       INCMT = 0
  319. C
  320.       DO 40 I=1,N
  321.          J=I-1
  322.          SUMA = 0.0D0
  323.          ISP = INCMT
  324.          INCMT = INCMT + I
  325.          LSTOP = N+1-I
  326. C
  327.          IF(I .EQ. 1) GO TO 25
  328. C
  329.          DO 20 L=1,J
  330.            L1 = ISP + L
  331.    20      SUMA = SUMA + DABS(A(L1))
  332.    25    ISP = ISP + 1
  333.          DO 30 L=1,LSTOP
  334.             ISP = ISP + J
  335.             SUMA = SUMA + DABS(A(ISP))
  336.    30       J = J+1
  337. C
  338.    40    NORMA = DMAX1(NORMA,SUMA)
  339. C
  340.       IF(NORMA .EQ. 0.0D0) NORMA = 1.0D0
  341. C
  342.       DO 100 I=1,M
  343.          S = 0.0D0
  344.          SUMZ = 0.0D0
  345.          INCMT = 0
  346. C
  347.          DO 65 L=1,N
  348.             SUMZ = SUMZ + DABS(Z(L,I))
  349.             SUM = -W(I)*Z(L,I)
  350.             J = L-1
  351.             ISP = INCMT
  352.             INCMT = INCMT + L
  353.             KSTOP = N+1-L
  354. C
  355.             IF(L .EQ. 1) GO TO 55
  356. C
  357.             DO 50 K=1,J
  358.               K1 = ISP + K
  359.    50         SUM = SUM + A(K1)*Z(K,I)
  360.    55       ISP = ISP + 1
  361.             DO 60 K=1,KSTOP
  362.               ISP =ISP + J
  363.               K2 = K - 1 + L
  364.               SUM =SUM + A(ISP)*Z(K2,I)
  365.    60       J = J + 1
  366.    65     S = S + DABS(SUM)
  367. C
  368.          NORM(I) = SUMZ
  369. C
  370.          IF(SUMZ .EQ. 0.0D0) GO TO 100
  371. C        ..........THIS LOOP WILL NEVER BE COMPLETED SINCE THERE
  372. C                  WILL ALWAYS EXIST AN ELEMENT IN THE VECTOR Z(I)
  373. C                  LARGER THAN !!Z(I)!!/N..........
  374.          DO 70 L=1,N
  375.             IF(DABS(Z(L,I)) .GE. NORM(I)/N) GO TO 80
  376.    70       CONTINUE
  377. C
  378.    80    TNORM = DSIGN(NORM(I),Z(L,I))
  379. C
  380.          DO 90 L=1,N
  381.    90       Z(L,I) = Z(L,I)/TNORM
  382. C
  383.          NORM(I) = S/(NORM(I)*NORMA)
  384.   100    RESDUL = DMAX1(NORM(I),RESDUL)
  385. C
  386.       RETURN
  387.       END
  388.       SUBROUTINE RMATIN(NNM,N,A,AHOLD,INITIL)
  389. C
  390. C     THIS INPUT SUBROUTINE READS A REAL SYMMETRIC MATRIX FROM SYSIN OF
  391. C     ORDER N AND STORES IT IN PACKED FORM. THE INPUT MATRIX IS READ
  392. C     AS A FULL MATRIX  (MM IS EQUAL TO ZERO)  OR AS ROWS OF AN UPPER
  393. C     TRIANGULAR MATRIX  (MM IS NOT EQUAL TO ZERO) .
  394. C     TO GENERATE THE MATRIX  A  INITIALLY,  INITIL  IS TO BE 0.
  395. C     TO REGENERATE THE MATRIX  A  FOR THE PURPOSE OF THE RESIDUAL
  396. C     CALCULATION,  INITIL  IS TO BE  1.
  397. C
  398. C     THIS ROUTINE IS CATALOGUED AS  EISPDRV4(RSPREADI).
  399. C
  400.       DOUBLE PRECISION A(NNM),AHOLD(NNM) 
  401.       INTEGER  IA( 20)
  402.       DATA IREADA/1/,IWRITE/6/
  403. C
  404.       IF( INITIL .EQ. 1 )  GO TO  30
  405.       READ(IREADA,5) N,MM
  406.     5 FORMAT(I6,6X,I6)
  407.       IF( N .EQ. 0 )  GO TO  70
  408.       INCMT=1
  409.       DO 15 L=1,N
  410.          JST = L
  411.          IF( MM .EQ. 0 )  JST = 1
  412.          READ(IREADA,10) (IA(J),J=JST,N) 
  413.    10    FORMAT(6I12)
  414.          J = L - 1
  415.          ISP = INCMT
  416.          INCMT =INCMT + L
  417.          DO 12 K = L,N
  418.            ISP = ISP + J
  419.            A(ISP) = IA(J + 1)
  420.    12      J = J + 1
  421.    15    CONTINUE
  422.       NNN = (N*(N+1))/2
  423.       DO 20 I = 1,NNN
  424.    20   AHOLD(I) = A(I)
  425.       RETURN
  426.    30 NNN = (N*(N+1))/2
  427.       DO 40 I = 1,NNN
  428.    40    A(I) = AHOLD(I)
  429.       RETURN
  430.    70 WRITE(IWRITE,80)
  431.    80 FORMAT(47H0END OF DATA FOR SUBROUTINE  RMATIN(RSPREADI). /1H1)
  432.       STOP
  433.       END
  434.