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