home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 9 / 09.iso / l / l078 / 1.img / HILBD.BAS < prev    next >
Encoding:
BASIC Source File  |  1987-04-02  |  5.1 KB  |  189 lines
  1. '┌───────────────────────────────────────────────────────────────────────────┐
  2. '│                HILBD.BAS                     │
  3. '│                   Version 1.0                     │
  4. '│                                         │
  5. '│    This program demonstrates the increased speed and precision         │
  6. '│    of the TURBO-BASIC compiler:                         │
  7. '│                                         │
  8. '│    --------------------------------------------------             │
  9. '│    From: BASIC Programs for Scientists and Engineers                      │
  10. '│                                         │
  11. '│    Alan R. Miller, Sybex                             │
  12. '│    n x n inverse hilbert matrix                         │
  13. '│    solution is 1 1 1 1 1                             │
  14. '│    double-precision version                             │
  15. '│    --------------------------------------------------                     │
  16. '│                                         │
  17. '│    The program performs simultaneous solution by Gauss-Jordan           │
  18. '│    elimination.                                 │
  19. '│                                         │
  20. '│    In order to run this program do the following:                 │
  21. '│      1. Load Turbo Basic by typing TB at the DOS prompt.              │
  22. '│      2. Load the file HILBD.BAS from the Load option of the File         │
  23. '│         pulldown menu.                             │
  24. '│      3. Select Run from the Main menu                     │
  25. '└───────────────────────────────────────────────────────────────────────────┘
  26.  
  27.  
  28. DEFINT I-N: DEFDBL A-H: DEFDBL O-Z
  29. DIM Z(12), A(12,12), COEF(12), B(12,12)
  30.  
  31. CLS
  32. Max = 12
  33. Nrow = 2
  34. A(1,1) = 1.0D0
  35. PRINT "Inversion of Hilbert Matrices"
  36. DO
  37.   CALL GetData(A(), Z(), Nrow)
  38.   FOR I = 1 TO Nrow
  39.     FOR J = 1 TO Nrow
  40.       B(I,J) = A(I,J)
  41.     NEXT J
  42.   NEXT I
  43.   CALL Gaussj(B(), Z(), Coef(), Nrow, Ierr)
  44.   IF Ierr THEN END                    ' Singular matrix
  45.   CALL WriteData(A(), Z(), Coef(), Nrow)
  46.   Nrow = Nrow+1
  47. LOOP UNTIL Nrow = Max
  48. END ' done
  49.  
  50. SUB GetData(A(2), Z(1), Nrow)
  51.   FOR I = 1 TO Nrow-1
  52.     A(Nrow,I) = 1.0D0 / (Nrow+I-1)
  53.     A(I,Nrow) = A(Nrow,I)
  54.   NEXT I
  55.   A(Nrow,Nrow) = 1.0D0 / (2*Nrow -1)
  56.   FOR I = 1 TO Nrow
  57.     Z(I) = 0
  58.     FOR J = 1 TO Nrow
  59.       Z(I) = Z(I) + A(I,J)
  60.     NEXT J
  61.   NEXT I
  62. END SUB                               ' GetDataput
  63.  
  64. SUB WriteData(A(2), Z(1), Coef(1), N)
  65.   SHARED Determ
  66.   IF N < 6 THEN 'show only smaller sets
  67.     PRINT
  68.     PRINT "         Matrix   Constants"
  69.     FOR I = 1 TO N
  70.       FOR J = 1 TO N
  71.        PRINT USING " #.####^^^^^ "; A(I,J);
  72.       NEXT J
  73.       PRINT USING " = #.####^^^^^"; Z(I)
  74.     NEXT I
  75.   END IF
  76.   PRINT
  77.   PRINT "    Solution for "; N; " equations,";
  78.   PRINT " Determinant =";
  79.   PRINT USING " #.#####^^^^^"; Determ
  80.   FOR I = 1 TO N
  81.     PRINT USING " #.######"; Coef(I);
  82.   NEXT I
  83.   PRINT
  84. END SUB ' WriteData
  85.  
  86. '
  87. 'Gauss-Jordan matrix inversion and solution
  88. '
  89. '  B(N,N) is coeff matrix, becomes inverse
  90. '  Y(N) is original constant vector
  91. '  W(N,Nvec) has constant vectors, becomes solution
  92. '  Determ is determinant
  93. '
  94. SUB Gaussj(B(2), Y(1), Coef(1), N, Ierr)
  95.   SHARED Determ                       ' Determinant available if needed
  96.   DIM W(8,1), Index(8,3)
  97.   False% = 0: True% = NOT False%
  98.   Ierr = False%                       'becomes 1 for singular matrix
  99.   Invers = 1                          'print inverse matrix if zero
  100.   Nvec = 1                            'number of constant vectors
  101.   FOR I = 1 TO N
  102.     W(I,1) = Y(I)
  103.     Index(I,3) = 0
  104.   NEXT I
  105.   Determ = 1.0
  106.   FOR I = 1 TO N
  107.     'search for largest (pivot) element
  108.     Big = 0.0
  109.     FOR J = 1 TO N
  110.     IF (Index(J,3) <> 1) THEN
  111.       FOR K = 1 TO N
  112.         IF (Index(K,3) > 1) THEN Emess
  113.         IF (Index(K,3) <> 1) THEN
  114.           IF (Big < ABS(B(J,K))) THEN
  115.             Irow = J
  116.             Icol = K
  117.             Big = ABS(B(J,K))
  118.           END IF
  119.         END IF
  120.        NEXT K
  121.      END IF
  122.    NEXT J
  123.    Index(Icol,3) = Index(Icol,3) + 1
  124.    Index(I,1) = Irow
  125.    Index(I,2) = Icol
  126.    'interchange rows to put pivot on diagonal
  127.    IF (Irow <> Icol) THEN
  128.      Determ = - Determ
  129.      FOR L = 1 TO N
  130.        SWAP B(Irow,L), B(Icol,L)
  131.      NEXT L
  132.      IF Nvec > 0 THEN
  133.        FOR L = 1 TO Nvec
  134.          SWAP W(Irow,L), W(Icol,L)
  135.        NEXT L
  136.      END IF
  137.    END IF
  138.    'divide pivot row by pivot element
  139.    Pivot = B(Icol, Icol)
  140.    IF Pivot = 0.0 THEN Emess
  141.    Determ = Determ * Pivot
  142.    B(Icol,Icol) = 1.0
  143.    FOR L=1 TO N
  144.      B(Icol,L) = B(Icol,L) / Pivot
  145.    NEXT L
  146.    IF Nvec > 0 THEN
  147.      FOR L = 1 TO Nvec
  148.        W(Icol,L) = W(Icol,L) / Pivot
  149.      NEXT L
  150.    END IF
  151.    'reduce nonpivot rows
  152.    FOR L1 = 1 TO N
  153.     IF (L1 <> Icol) THEN
  154.       T = B(L1, Icol)
  155.       B(L1,Icol) = 0.0
  156.       FOR L = 1 TO N
  157.         B(L1,L) = B(L1,L) - B(Icol,L) * T
  158.       NEXT L
  159.       IF Nvec > 0 THEN
  160.         FOR L = 1 TO Nvec
  161.           W(L1,L) = W(L1,L) - W(Icol,L) * T
  162.         NEXT L
  163.       END IF
  164.     END IF
  165.    NEXT L1
  166.   NEXT I
  167.   'interchange columns
  168.   FOR I = 1 TO N
  169.     L = N - I + 1
  170.     IF (Index(L,1) <> Index(L,2)) THEN
  171.       Irow = Index(L,1)
  172.       Icol = Index(L,2)
  173.       FOR K = 1 TO N
  174.         SWAP B(K,Irow), B(K,Icol)
  175.       NEXT K
  176.     END IF
  177.   NEXT I
  178.   FOR K = 1 TO N
  179.     IF (Index(K,3) <> 1) THEN Emess
  180.   NEXT K
  181.   FOR I = 1 TO N
  182.     Coef(I) = W(I,1)
  183.   NEXT I
  184.   EXIT SUB                            ' normally
  185. Emess:
  186.   Ierr = True%
  187.   PRINT "Matrix singular "
  188. END SUB                               ' Gauss-Jordan
  189.