home *** CD-ROM | disk | FTP | other *** search
/ ProfitPress Mega CDROM2 …eeware (MSDOS)(1992)(Eng) / ProfitPress-MegaCDROM2.B6I / APPS / DATABASE / AMSF20.ZIP / DGELG.FOR < prev    next >
Encoding:
Text File  |  1992-01-06  |  4.6 KB  |  170 lines
  1. C
  2. C     ..................................................................
  3. C
  4. C        SUBROUTINE DGELG
  5. C
  6. C        PURPOSE
  7. C           TO SOLVE A GENERAL SYSTEM OF SIMULTANEOUS LINEAR EQUATIONS.
  8. C
  9. C        USAGE
  10. C           CALL DGELG(R,A,M,N,EPS,IER)
  11. C
  12. C        DESCRIPTION OF PARAMETERS
  13. C           R      - DOUBLE PRECISION M BY N RIGHT HAND SIDE MATRIX
  14. C                    (DESTROYED). ON RETURN R CONTAINS THE SOLUTIONS
  15. C                    OF THE EQUATIONS.
  16. C           A      - DOUBLE PRECISION M BY M COEFFICIENT MATRIX
  17. C                    (DESTROYED).
  18. C           M      - THE NUMBER OF EQUATIONS IN THE SYSTEM.
  19. C           N      - THE NUMBER OF RIGHT HAND SIDE VECTORS.
  20. C           EPS    - SINGLE PRECISION INPUT CONSTANT WHICH IS USED AS
  21. C                    RELATIVE TOLERANCE FOR TEST ON LOSS OF
  22. C                    SIGNIFICANCE.
  23. C           IER    - RESULTING ERROR PARAMETER CODED AS FOLLOWS
  24. C                    IER=0  - NO ERROR,
  25. C                    IER=-1 - NO RESULT BECAUSE OF M LESS THAN 1 OR
  26. C                             PIVOT ELEMENT AT ANY ELIMINATION STEP
  27. C                             EQUAL TO 0,
  28. C                    IER=K  - WARNING DUE TO POSSIBLE LOSS OF SIGNIFI-
  29. C                             CANCE INDICATED AT ELIMINATION STEP K+1,
  30. C                             WHERE PIVOT ELEMENT WAS LESS THAN OR
  31. C                             EQUAL TO THE INTERNAL TOLERANCE EPS TIMES
  32. C                             ABSOLUTELY GREATEST ELEMENT OF MATRIX A.
  33. C
  34. C        REMARKS
  35. C           INPUT MATRICES R AND A ARE ASSUMED TO BE STORED COLUMNWISE
  36. C           IN M*N RESP. M*M SUCCESSIVE STORAGE LOCATIONS. ON RETURN
  37. C           SOLUTION MATRIX R IS STORED COLUMNWISE TOO.
  38. C           THE PROCEDURE GIVES RESULTS IF THE NUMBER OF EQUATIONS M IS
  39. C           GREATER THAN 0 AND PIVOT ELEMENTS AT ALL ELIMINATION STEPS
  40. C           ARE DIFFERENT FROM 0. HOWEVER WARNING IER=K - IF GIVEN -
  41. C           INDICATES POSSIBLE LOSS OF SIGNIFICANCE. IN CASE OF A WELL
  42. C           SCALED MATRIX A AND APPROPRIATE TOLERANCE EPS, IER=K MAY BE
  43. C           INTERPRETED THAT MATRIX A HAS THE RANK K. NO WARNING IS
  44. C           GIVEN IN CASE M=1.
  45. C
  46. C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
  47. C           NONE
  48. C
  49. C        METHOD
  50. C           SOLUTION IS DONE BY MEANS OF GAUSS-ELIMINATION WITH
  51. C           COMPLETE PIVOTING.
  52. C
  53. C     ..................................................................
  54. C
  55.       SUBROUTINE DGELG(R,A,M,N,EPS,IER)
  56. C
  57. C
  58.       DIMENSION A(1),R(1)
  59.       DOUBLE PRECISION R,A,PIV,TB,TOL,PIVI
  60.       IF(M)23,23,1
  61. C
  62. C     SEARCH FOR GREATEST ELEMENT IN MATRIX A
  63.     1 IER=0
  64.       PIV=0.D0
  65.       MM=M*M
  66.       NM=N*M
  67.       DO 3 L=1,MM
  68.       TB=DABS(A(L))
  69.       IF(TB-PIV)3,3,2
  70.     2 PIV=TB
  71.       I=L
  72.     3 CONTINUE
  73.       TOL=EPS*PIV
  74. C     A(I) IS PIVOT ELEMENT. PIV CONTAINS THE ABSOLUTE VALUE OF A(I).
  75. C
  76. C
  77. C     START ELIMINATION LOOP
  78.       LST=1
  79.       DO 17 K=1,M
  80. C
  81. C     TEST ON SINGULARITY
  82.       IF(PIV)23,23,4
  83.     4 IF(IER)7,5,7
  84.     5 IF(PIV-TOL)6,6,7
  85.     6 IER=K-1
  86.     7 PIVI=1.D0/A(I)
  87.       J=(I-1)/M
  88.       I=I-J*M-K
  89.       J=J+1-K
  90. C     I+K IS ROW-INDEX, J+K COLUMN-INDEX OF PIVOT ELEMENT
  91. C
  92. C     PIVOT ROW REDUCTION AND ROW INTERCHANGE IN RIGHT HAND SIDE R
  93.       DO 8 L=K,NM,M
  94.       LL=L+I
  95.       TB=PIVI*R(LL)
  96.       R(LL)=R(L)
  97.     8 R(L)=TB
  98. C
  99. C     IS ELIMINATION TERMINATED
  100.       IF(K-M)9,18,18
  101. C
  102. C     COLUMN INTERCHANGE IN MATRIX A
  103.     9 LEND=LST+M-K
  104.       IF(J)12,12,10
  105.    10 II=J*M
  106.       DO 11 L=LST,LEND
  107.       TB=A(L)
  108.       LL=L+II
  109.       A(L)=A(LL)
  110.    11 A(LL)=TB
  111. C
  112. C     ROW INTERCHANGE AND PIVOT ROW REDUCTION IN MATRIX A
  113.    12 DO 13 L=LST,MM,M
  114.       LL=L+I
  115.       TB=PIVI*A(LL)
  116.       A(LL)=A(L)
  117.    13 A(L)=TB
  118. C
  119. C     SAVE COLUMN INTERCHANGE INFORMATION
  120.       A(LST)=J
  121. C
  122. C     ELEMENT REDUCTION AND NEXT PIVOT SEARCH
  123.       PIV=0.D0
  124.       LST=LST+1
  125.       J=0
  126.       DO 16 II=LST,LEND
  127.       PIVI=-A(II)
  128.       IST=II+M
  129.       J=J+1
  130.       DO 15 L=IST,MM,M
  131.       LL=L-J
  132.       A(L)=A(L)+PIVI*A(LL)
  133.       TB=DABS(A(L))
  134.       IF(TB-PIV)15,15,14
  135.    14 PIV=TB
  136.       I=L
  137.    15 CONTINUE
  138.       DO 16 L=K,NM,M
  139.       LL=L+J
  140.    16 R(LL)=R(LL)+PIVI*R(L)
  141.    17 LST=LST+M
  142. C     END OF ELIMINATION LOOP
  143. C
  144. C
  145. C     BACK SUBSTITUTION AND BACK INTERCHANGE
  146.    18 IF(M-1)23,22,19
  147.    19 IST=MM+M
  148.       LST=M+1
  149.       DO 21 I=2,M
  150.       II=LST-I
  151.       IST=IST-LST
  152.       L=IST-M
  153.       L=A(L)+.5D0
  154.       DO 21 J=II,NM,M
  155.       TB=R(J)
  156.       LL=J
  157.       DO 20 K=IST,MM,M
  158.       LL=LL+1
  159.    20 TB=TB-A(K)*R(LL)
  160.       K=J+L
  161.       R(J)=R(K)
  162.    21 R(K)=TB
  163.    22 RETURN
  164. C
  165. C
  166. C     ERROR RETURN
  167.    23 IER=-1
  168.       RETURN
  169.       END
  170.