home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD1.mdf / fortran / library / ssp / matinsl / rslmc.for < prev    next >
Encoding:
Text File  |  1985-11-29  |  4.3 KB  |  151 lines
  1. C
  2. C     ..................................................................
  3. C
  4. C        SUBROUTINE RSLMC
  5. C
  6. C        PURPOSE
  7. C           SOLUTION OF A SYSTEM OF LINEAR EQUATIONS AX=B
  8. C
  9. C        USAGE
  10. C           CALL RSLMC(A,AF,B,X,N,EPSI,IER,IA,V,PER)
  11. C
  12. C        DESCRIPTION OF PARAMETERS
  13. C           A      INPUT MATRIX
  14. C           AF     ARRAY OF THE FACTORIZATION OF THE ORIGINAL MATRIX
  15. C           B      RIGHT HAND SIDE VECTOR
  16. C           X      VECTOR CONTAINING THE SOLUTION ON RETURN
  17. C           N      ORDER OF THE SYSTEM
  18. C           EPSI   RELATIVE PRECISION INDICATOR(REQUIRED INPUT)
  19. C           IER    ERROR INDICATOR
  20. C                     =0 IF EACH COMPONENT OF X MEETS THE PRECISION EPSI
  21. C                     =1 IF ONLY THE NORM OF X MEETS THIS PRECISION
  22. C                     =2 IF THE PRECISION IN THE NORM OF THE COMPUTED
  23. C                        SOLUTION IS LOWER THAN EPSI
  24. C                     =3 IF THE SOLUTION OBTAINED HAS NO MEANING AT ALL
  25. C                     =4 IF A DIAGONAL TERM OF THE UPPER TRIANGULAR
  26. C                        FACTOR IS ZERO
  27. C           IA     SIZE OF THE FIRST DIMENSION ASSIGNED TO THE ARRAY A
  28. C                  IN THE CALLING PROGRAM WHEN THE MATRIX IS IN DOUBLE
  29. C                  SUBSCRIPTED DATA STORAGE MODE.  IA=N WHEN THE MATRIX
  30. C                  IS IN SSP VECTOR STORAGE MODE.
  31. C           V      WORKING STORAGE VECTOR
  32. C                  DIMENSION OF V MUST BE GREATER THAN OR EQUAL TO N
  33. C           PER    VECTOR WHERE PERMUTATIONS OF ROWS OF THE MATRIX ARE
  34. C                  STORED
  35. C                  DIMENSION OF PER MUST BE GREATER THAN OR EQUAL TO N
  36. C
  37. C        REMARKS
  38. C           THE MATRIX OF THE SYSTEM MAY BE FACTORIZED BY THE SUBROUTINE
  39. C           FACTR IN THE ARRAY AF PRIOR TO ENTRY TO THIS SUBROUTINE.
  40. C           THE LOWER TRIANGULAR FACTOR MUST HAVE AN UNIT DIAGONAL.
  41. C           EPSI IS MODIFIED WHEN IER=2
  42. C
  43. C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
  44. C           NONE
  45. C
  46. C        METHOD
  47. C           A TRIAL SOLUTION IS FIRST COMPUTED.  THEN CORRECTIONS ARE
  48. C           CALCULATED FROM RESIDUAL VECTORS.
  49. C
  50. C        REFERENCES
  51. C           J. H. WILKINSON - THE ALGEBRAIC EIGENVALUE PROBLEM -
  52. C           CLARENDON PRESS OXFORD, 1965.  H. J. BOWDLER, R. S. MARTIN,
  53. C           G. PETERS, AND J. H. WILKINSON.  'SOLUTION OF REAL AND
  54. C           COMPLEX SYSTEMS OF LINEAR EQUATIONS', NUMERISCHE MATHEMATIK,
  55. C           VOL. 8, NO. 3, 1966, 217-234.
  56. C
  57. C     ..................................................................
  58. C
  59.       SUBROUTINE RSLMC (A,AF,B,X,N,EPSI,IER,IA,V,PER)
  60.       DIMENSION A(1),AF(1),B(1),X(1),V(1),PER(1)
  61.       DOUBLE PRECISION DP
  62. C
  63. C        INITIALIZATION
  64. C
  65.       D0=0.
  66.       IER=0
  67.       ITE=0
  68.       DO 10 I=1,N
  69.       V(I)=B(I)
  70.    10 X(I)=0.
  71.    20 ITE=ITE+1
  72. C
  73. C        THE PERMUTATIONS OF ROWS OF A ARE APPLIED TO V
  74. C
  75.       DO 30 I=1,N
  76.       K=PER(I)
  77.       IF (K-I)25,30,25
  78.    25 D1=V(K)
  79.       V(K)=V(I)
  80.       V(I)=D1
  81.    30 CONTINUE
  82. C
  83. C        SOLUTION OF THE LOWER TRIANGULAR SYSTEM
  84. C
  85.       DO 50 I=2,N
  86.       IM1=I-1
  87.       DP=V(I)
  88.       IK=I
  89.       DO 40 K=1,IM1
  90.       DP=DP-1.D0*AF(IK)*V(K)
  91.    40 IK=IK+IA
  92.    50 V(I)=DP
  93. C
  94. C        SOLUTION OF THE UPPER TRIANGULAR SYSTEM
  95. C
  96.       IF(AF(IK)) 58,54,58
  97.    54 IER=4
  98.       GO TO 82
  99.    58 V(N)=DP/AF(IK)
  100.       DO 70 I=2,N
  101.       IM1=N-I+1
  102.       INF=IM1+1
  103.       DP=V(IM1)
  104.       IK=(IM1-1)*IA+IM1
  105.       D1=AF(IK)
  106.       DO 60 K=INF,N
  107.       IK=IK+IA
  108.    60 DP=DP-1.D0*AF(IK)*V(K)
  109.    70 V(IM1)=DP/D1
  110. C
  111. C        TEST OF PRECISION
  112. C
  113.       D1=0.
  114.       D2=0.
  115.       KLE=0
  116.       DO 80 I=1,N
  117.       D1=D1+ABS(V(I))
  118.       D2=D2+ABS(X(I))
  119.       IF (ABS(V(I))-EPSI*ABS(X(I))) 80,80,75
  120.    75 KLE=1
  121.    80 CONTINUE
  122.       IF (KLE)140,82,85
  123.    82 RETURN
  124.    85 IF (ITE-1)140,90,87
  125. C
  126. C        ITERATIONS ARE STOPPED WHEN THE NORM OF THE CORRECTION IS MORE
  127. C        THAN HALF OF THE ONE OF THE FORMER
  128. C
  129.    87 IF (D0-2.*D1)120,90,90
  130.    90 DO 95 I=1,N
  131.    95 X(I)=X(I)+V(I)
  132.       DO 110 I=1,N
  133.       DP=B(I)
  134.       IK=I
  135.       DO 100 K=1,N
  136.       DP=DP-1.D0*A(IK)*X(K)
  137.   100 IK=IK+IA
  138.   110 V(I)=DP
  139.       D0=D1
  140.       GO TO 20
  141.   120 IF(ITE-2)140,140,125
  142.   125 IF (D1-EPSI*D2)127,127,130
  143.   127 IER=1
  144.       RETURN
  145.   130 IER=2
  146.       EPSI=D1/D2
  147.       RETURN
  148.   140 IER=3
  149.       RETURN
  150.       END
  151.