Google
 

Trailing-Edge - PDP-10 Archives - decuslib20-02 - decus/20-0026/gelg.ssp
There are 2 other files named gelg.ssp in the archive. Click here to see a list.
C                                                                       GELG  10
C     ..................................................................GELG  20
C                                                                       GELG  30
C        SUBROUTINE GELG                                                GELG  40
C                                                                       GELG  50
C        PURPOSE                                                        GELG  60
C           TO SOLVE A GENERAL SYSTEM OF SIMULTANEOUS LINEAR EQUATIONS. GELG  70
C                                                                       GELG  80
C        USAGE                                                          GELG  90
C           CALL GELG(R,A,M,N,EPS,IER)                                  GELG 100
C                                                                       GELG 110
C        DESCRIPTION OF PARAMETERS                                      GELG 120
C           R      - THE M BY N MATRIX OF RIGHT HAND SIDES.  (DESTROYED)GELG 130
C                    ON RETURN R CONTAINS THE SOLUTION OF THE EQUATIONS.GELG 140
C           A      - THE M BY M COEFFICIENT MATRIX.  (DESTROYED)        GELG 150
C           M      - THE NUMBER OF EQUATIONS IN THE SYSTEM.             GELG 160
C           N      - THE NUMBER OF RIGHT HAND SIDE VECTORS.             GELG 170
C           EPS    - AN INPUT CONSTANT WHICH IS USED AS RELATIVE        GELG 180
C                    TOLERANCE FOR TEST ON LOSS OF SIGNIFICANCE.        GELG 190
C           IER    - RESULTING ERROR PARAMETER CODED AS FOLLOWS         GELG 200
C                    IER=0  - NO ERROR,                                 GELG 210
C                    IER=-1 - NO RESULT BECAUSE OF M LESS THAN 1 OR     GELG 220
C                             PIVOT ELEMENT AT ANY ELIMINATION STEP     GELG 230
C                             EQUAL TO 0,                               GELG 240
C                    IER=K  - WARNING DUE TO POSSIBLE LOSS OF SIGNIFI-  GELG 250
C                             CANCE INDICATED AT ELIMINATION STEP K+1,  GELG 260
C                             WHERE PIVOT ELEMENT WAS LESS THAN OR      GELG 270
C                             EQUAL TO THE INTERNAL TOLERANCE EPS TIMES GELG 280
C                             ABSOLUTELY GREATEST ELEMENT OF MATRIX A.  GELG 290
C                                                                       GELG 300
C        REMARKS                                                        GELG 310
C           INPUT MATRICES R AND A ARE ASSUMED TO BE STORED COLUMNWISE  GELG 320
C           IN M*N RESP. M*M SUCCESSIVE STORAGE LOCATIONS. ON RETURN    GELG 330
C           SOLUTION MATRIX R IS STORED COLUMNWISE TOO.                 GELG 340
C           THE PROCEDURE GIVES RESULTS IF THE NUMBER OF EQUATIONS M IS GELG 350
C           GREATER THAN 0 AND PIVOT ELEMENTS AT ALL ELIMINATION STEPS  GELG 360
C           ARE DIFFERENT FROM 0. HOWEVER WARNING IER=K - IF GIVEN -    GELG 370
C           INDICATES POSSIBLE LOSS OF SIGNIFICANCE. IN CASE OF A WELL  GELG 380
C           SCALED MATRIX A AND APPROPRIATE TOLERANCE EPS, IER=K MAY BE GELG 390
C           INTERPRETED THAT MATRIX A HAS THE RANK K. NO WARNING IS     GELG 400
C           GIVEN IN CASE M=1.                                          GELG 410
C                                                                       GELG 420
C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED                  GELG 430
C           NONE                                                        GELG 440
C                                                                       GELG 450
C        METHOD                                                         GELG 460
C           SOLUTION IS DONE BY MEANS OF GAUSS-ELIMINATION WITH         GELG 470
C           COMPLETE PIVOTING.                                          GELG 480
C                                                                       GELG 490
C     ..................................................................GELG 500
C                                                                       GELG 510
      SUBROUTINE GELG(R,A,M,N,EPS,IER)                                  GELG 520
C                                                                       GELG 530
C                                                                       GELG 540
      DIMENSION A(1),R(1)                                               GELG 550
      IF(M)23,23,1                                                      GELG 560
C                                                                       GELG 570
C     SEARCH FOR GREATEST ELEMENT IN MATRIX A                           GELG 580
    1 IER=0                                                             GELG 590
      PIV=0.                                                            GELG 600
      MM=M*M                                                            GELG 610
      NM=N*M                                                            GELG 620
      DO 3 L=1,MM                                                       GELG 630
      TB=ABS(A(L))                                                      GELG 640
      IF(TB-PIV)3,3,2                                                   GELG 650
    2 PIV=TB                                                            GELG 660
      I=L                                                               GELG 670
    3 CONTINUE                                                          GELG 680
      TOL=EPS*PIV                                                       GELG 690
C     A(I) IS PIVOT ELEMENT. PIV CONTAINS THE ABSOLUTE VALUE OF A(I).   GELG 700
C                                                                       GELG 710
C                                                                       GELG 720
C     START ELIMINATION LOOP                                            GELG 730
      LST=1                                                             GELG 740
      DO 17 K=1,M                                                       GELG 750
C                                                                       GELG 760
C     TEST ON SINGULARITY                                               GELG 770
      IF(PIV)23,23,4                                                    GELG 780
    4 IF(IER)7,5,7                                                      GELG 790
    5 IF(PIV-TOL)6,6,7                                                  GELG 800
    6 IER=K-1                                                           GELG 810
    7 PIVI=1./A(I)                                                      GELG 820
      J=(I-1)/M                                                         GELG 830
      I=I-J*M-K                                                         GELG 840
      J=J+1-K                                                           GELG 850
C     I+K IS ROW-INDEX, J+K COLUMN-INDEX OF PIVOT ELEMENT               GELG 860
C                                                                       GELG 870
C     PIVOT ROW REDUCTION AND ROW INTERCHANGE IN RIGHT HAND SIDE R      GELG 880
      DO 8 L=K,NM,M                                                     GELG 890
      LL=L+I                                                            GELG 900
      TB=PIVI*R(LL)                                                     GELG 910
      R(LL)=R(L)                                                        GELG 920
    8 R(L)=TB                                                           GELG 930
C                                                                       GELG 940
C     IS ELIMINATION TERMINATED                                         GELG 950
      IF(K-M)9,18,18                                                    GELG 960
C                                                                       GELG 970
C     COLUMN INTERCHANGE IN MATRIX A                                    GELG 980
    9 LEND=LST+M-K                                                      GELG 990
      IF(J)12,12,10                                                     GELG1000
   10 II=J*M                                                            GELG1010
      DO 11 L=LST,LEND                                                  GELG1020
      TB=A(L)                                                           GELG1030
      LL=L+II                                                           GELG1040
      A(L)=A(LL)                                                        GELG1050
   11 A(LL)=TB                                                          GELG1060
C                                                                       GELG1070
C     ROW INTERCHANGE AND PIVOT ROW REDUCTION IN MATRIX A               GELG1080
   12 DO 13 L=LST,MM,M                                                  GELG1090
      LL=L+I                                                            GELG1100
      TB=PIVI*A(LL)                                                     GELG1110
      A(LL)=A(L)                                                        GELG1120
   13 A(L)=TB                                                           GELG1130
C                                                                       GELG1140
C     SAVE COLUMN INTERCHANGE INFORMATION                               GELG1150
      A(LST)=J                                                          GELG1160
C                                                                       GELG1170
C     ELEMENT REDUCTION AND NEXT PIVOT SEARCH                           GELG1180
      PIV=0.                                                            GELG1190
      LST=LST+1                                                         GELG1200
      J=0                                                               GELG1210
      DO 16 II=LST,LEND                                                 GELG1220
      PIVI=-A(II)                                                       GELG1230
      IST=II+M                                                          GELG1240
      J=J+1                                                             GELG1250
      DO 15 L=IST,MM,M                                                  GELG1260
      LL=L-J                                                            GELG1270
      A(L)=A(L)+PIVI*A(LL)                                              GELG1280
      TB=ABS(A(L))                                                      GELG1290
      IF(TB-PIV)15,15,14                                                GELG1300
   14 PIV=TB                                                            GELG1310
      I=L                                                               GELG1320
   15 CONTINUE                                                          GELG1330
      DO 16 L=K,NM,M                                                    GELG1340
      LL=L+J                                                            GELG1350
   16 R(LL)=R(LL)+PIVI*R(L)                                             GELG1360
   17 LST=LST+M                                                         GELG1370
C     END OF ELIMINATION LOOP                                           GELG1380
C                                                                       GELG1390
C                                                                       GELG1400
C     BACK SUBSTITUTION AND BACK INTERCHANGE                            GELG1410
   18 IF(M-1)23,22,19                                                   GELG1420
   19 IST=MM+M                                                          GELG1430
      LST=M+1                                                           GELG1440
      DO 21 I=2,M                                                       GELG1450
      II=LST-I                                                          GELG1460
      IST=IST-LST                                                       GELG1470
      L=IST-M                                                           GELG1480
      L=A(L)+.5                                                         GELG1490
      DO 21 J=II,NM,M                                                   GELG1500
      TB=R(J)                                                           GELG1510
      LL=J                                                              GELG1520
      DO 20 K=IST,MM,M                                                  GELG1530
      LL=LL+1                                                           GELG1540
   20 TB=TB-A(K)*R(LL)                                                  GELG1550
      K=J+L                                                             GELG1560
      R(J)=R(K)                                                         GELG1570
   21 R(K)=TB                                                           GELG1580
   22 RETURN                                                            GELG1590
C                                                                       GELG1600
C                                                                       GELG1610
C     ERROR RETURN                                                      GELG1620
   23 IER=-1                                                            GELG1630
      RETURN                                                            GELG1640
      END                                                               GELG1650