Google
 

Trailing-Edge - PDP-10 Archives - decuslib20-02 - decus/20-0026/gelb.ssp
There are 2 other files named gelb.ssp in the archive. Click here to see a list.
C                                                                       GELB  10
C     ..................................................................GELB  20
C                                                                       GELB  30
C        SUBROUTINE GELB                                                GELB  40
C                                                                       GELB  50
C        PURPOSE                                                        GELB  60
C           TO SOLVE A SYSTEM OF SIMULTANEOUS LINEAR EQUATIONS WITH A   GELB  70
C           COEFFICIENT MATRIX OF BAND STRUCTURE.                       GELB  80
C                                                                       GELB  90
C        USAGE                                                          GELB 100
C           CALL GELB(R,A,M,N,MUD,MLD,EPS,IER)                          GELB 110
C                                                                       GELB 120
C        DESCRIPTION OF PARAMETERS                                      GELB 130
C           R      - M BY N RIGHT HAND SIDE MATRIX (DESTROYED).         GELB 140
C                    ON RETURN R CONTAINS THE SOLUTION OF THE EQUATIONS.GELB 150
C           A      - M BY M COEFFICIENT MATRIX WITH BAND STRUCTURE      GELB 160
C                    (DESTROYED).                                       GELB 170
C           M      - THE NUMBER OF EQUATIONS IN THE SYSTEM.             GELB 180
C           N      - THE NUMBER OF RIGHT HAND SIDE VECTORS.             GELB 190
C           MUD    - THE NUMBER OF UPPER CODIAGONALS (THAT MEANS        GELB 200
C                    CODIAGONALS ABOVE MAIN DIAGONAL).                  GELB 210
C           MLD    - THE NUMBER OF LOWER CODIAGONALS (THAT MEANS        GELB 220
C                    CODIAGONALS BELOW MAIN DIAGONAL).                  GELB 230
C           EPS    - AN INPUT CONSTANT WHICH IS USED AS RELATIVE        GELB 240
C                    TOLERANCE FOR TEST ON LOSS OF SIGNIFICANCE.        GELB 250
C           IER    - RESULTING ERROR PARAMETER CODED AS FOLLOWS         GELB 260
C                    IER=0  - NO ERROR,                                 GELB 270
C                    IER=-1 - NO RESULT BECAUSE OF WRONG INPUT PARAME-  GELB 280
C                             TERS M,MUD,MLD OR BECAUSE OF PIVOT ELEMENTGELB 290
C                             AT ANY ELIMINATION STEP EQUAL TO 0,       GELB 300
C                    IER=K  - WARNING DUE TO POSSIBLE LOSS OF SIGNIFI-  GELB 310
C                             CANCE INDICATED AT ELIMINATION STEP K+1,  GELB 320
C                             WHERE PIVOT ELEMENT WAS LESS THAN OR      GELB 330
C                             EQUAL TO THE INTERNAL TOLERANCE EPS TIMES GELB 340
C                             ABSOLUTELY GREATEST ELEMENT OF MATRIX A.  GELB 350
C                                                                       GELB 360
C        REMARKS                                                        GELB 370
C           BAND MATRIX A IS ASSUMED TO BE STORED ROWWISE IN THE FIRST  GELB 380
C           ME SUCCESSIVE STORAGE LOCATIONS OF TOTALLY NEEDED MA        GELB 390
C           STORAGE LOCATIONS, WHERE                                    GELB 400
C             MA=M*MC-ML*(ML+1)/2    AND    ME=MA-MU*(MU+1)/2    WITH   GELB 410
C             MC=MIN(M,1+MUD+MLD),  ML=MC-1-MLD,  MU=MC-1-MUD.          GELB 420
C           RIGHT HAND SIDE MATRIX R IS ASSUMED TO BE STORED COLUMNWISE GELB 430
C           IN N*M SUCCESSIVE STORAGE LOCATIONS. ON RETURN SOLUTION     GELB 440
C           MATRIX R IS STORED COLUMNWISE TOO.                          GELB 450
C           INPUT PARAMETERS M, MUD, MLD SHOULD SATISFY THE FOLLOWING   GELB 460
C           RESTRICTIONS     MUD NOT LESS THAN ZERO                     GELB 470
C                            MLD NOT LESS THAN ZERO                     GELB 480
C                            MUD+MLD NOT GREATER THAN 2*M-2.            GELB 490
C           NO ACTION BESIDES ERROR MESSAGE IER=-1 TAKES PLACE IF THESE GELB 500
C           RESTRICTIONS ARE NOT SATISFIED.                             GELB 510
C           THE PROCEDURE GIVES RESULTS IF THE RESTRICTIONS ON INPUT    GELB 520
C           PARAMETERS ARE SATISFIED AND IF PIVOT ELEMENTS AT ALL       GELB 530
C           ELIMINATION STEPS ARE DIFFERENT FROM 0. HOWEVER WARNING     GELB 540
C           IER=K - IF GIVEN - INDICATES POSSIBLE LOSS OF SIGNIFICANCE. GELB 550
C           IN CASE OF A WELL SCALED MATRIX A AND APPROPRIATE TOLERANCE GELB 560
C           EPS, IER=K MAY BE INTERPRETED THAT MATRIX A HAS THE RANK K. GELB 570
C           NO WARNING IS GIVEN IF MATRIX A HAS NO LOWER CODIAGONAL.    GELB 580
C                                                                       GELB 590
C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED                  GELB 600
C           NONE                                                        GELB 610
C                                                                       GELB 620
C        METHOD                                                         GELB 630
C           SOLUTION IS DONE BY MEANS OF GAUSS ELIMINATION WITH         GELB 640
C           COLUMN PIVOTING ONLY, IN ORDER TO PRESERVE BAND STRUCTURE   GELB 650
C           IN REMAINING COEFFICIENT MATRICES.                          GELB 660
C                                                                       GELB 670
C     ..................................................................GELB 680
C                                                                       GELB 690
      SUBROUTINE GELB(R,A,M,N,MUD,MLD,EPS,IER)                          GELB 700
C                                                                       GELB 710
C                                                                       GELB 720
      DIMENSION R(1),A(1)                                               GELB 730
C                                                                       GELB 740
C     TEST ON WRONG INPUT PARAMETERS                                    GELB 750
      IF(MLD)47,1,1                                                     GELB 760
    1 IF(MUD)47,2,2                                                     GELB 770
    2 MC=1+MLD+MUD                                                      GELB 780
      IF(MC+1-M-M)3,3,47                                                GELB 790
C                                                                       GELB 800
C     PREPARE INTEGER PARAMETERS                                        GELB 810
C        MC=NUMBER OF COLUMNS IN MATRIX A                               GELB 820
C        MU=NUMBER OF ZEROS TO BE INSERTED IN FIRST ROW OF MATRIX A     GELB 830
C        ML=NUMBER OF MISSING ELEMENTS IN LAST ROW OF MATRIX A          GELB 840
C        MR=INDEX OF LAST ROW IN MATRIX A WITH MC ELEMENTS              GELB 850
C        MZ=TOTAL NUMBER OF ZEROS TO BE INSERTED IN MATRIX A            GELB 860
C        MA=TOTAL NUMBER OF STORAGE LOCATIONS NECESSARY FOR MATRIX A    GELB 870
C        NM=NUMBER OF ELEMENTS IN MATRIX R                              GELB 880
    3 IF(MC-M)5,5,4                                                     GELB 890
    4 MC=M                                                              GELB 900
    5 MU=MC-MUD-1                                                       GELB 910
      ML=MC-MLD-1                                                       GELB 920
      MR=M-ML                                                           GELB 930
      MZ=(MU*(MU+1))/2                                                  GELB 940
      MA=M*MC-(ML*(ML+1))/2                                             GELB 950
      NM=N*M                                                            GELB 960
C                                                                       GELB 970
C     MOVE ELEMENTS BACKWARD AND SEARCH FOR ABSOLUTELY GREATEST ELEMENT GELB 980
C     (NOT NECESSARY IN CASE OF A MATRIX WITHOUT LOWER CODIAGONALS)     GELB 990
      IER=0                                                             GELB1000
      PIV=0.                                                            GELB1010
      IF(MLD)14,14,6                                                    GELB1020
    6 JJ=MA                                                             GELB1030
      J=MA-MZ                                                           GELB1040
      KST=J                                                             GELB1050
      DO 9 K=1,KST                                                      GELB1060
      TB=A(J)                                                           GELB1070
      A(JJ)=TB                                                          GELB1080
      TB=ABS(TB)                                                        GELB1090
      IF(TB-PIV)8,8,7                                                   GELB1100
    7 PIV=TB                                                            GELB1110
    8 J=J-1                                                             GELB1120
    9 JJ=JJ-1                                                           GELB1130
C                                                                       GELB1140
C     INSERT ZEROS IN FIRST MU ROWS (NOT NECESSARY IN CASE MZ=0)        GELB1150
      IF(MZ)14,14,10                                                    GELB1160
   10 JJ=1                                                              GELB1170
      J=1+MZ                                                            GELB1180
      IC=1+MUD                                                          GELB1190
      DO 13 I=1,MU                                                      GELB1200
      DO 12 K=1,MC                                                      GELB1210
      A(JJ)=0.                                                          GELB1220
      IF(K-IC)11,11,12                                                  GELB1230
   11 A(JJ)=A(J)                                                        GELB1240
      J=J+1                                                             GELB1250
   12 JJ=JJ+1                                                           GELB1260
   13 IC=IC+1                                                           GELB1270
C                                                                       GELB1280
C     GENERATE TEST VALUE FOR SINGULARITY                               GELB1290
   14 TOL=EPS*PIV                                                       GELB1300
C                                                                       GELB1310
C                                                                       GELB1320
C     START DECOMPOSITION LOOP                                          GELB1330
      KST=1                                                             GELB1340
      IDST=MC                                                           GELB1350
      IC=MC-1                                                           GELB1360
      DO 38 K=1,M                                                       GELB1370
      IF(K-MR-1)16,16,15                                                GELB1380
   15 IDST=IDST-1                                                       GELB1390
   16 ID=IDST                                                           GELB1400
      ILR=K+MLD                                                         GELB1410
      IF(ILR-M)18,18,17                                                 GELB1420
   17 ILR=M                                                             GELB1430
   18 II=KST                                                            GELB1440
C                                                                       GELB1450
C     PIVOT SEARCH IN FIRST COLUMN (ROW INDEXES FROM I=K UP TO I=ILR)   GELB1460
      PIV=0.                                                            GELB1470
      DO 22 I=K,ILR                                                     GELB1480
      TB=ABS(A(II))                                                     GELB1490
      IF(TB-PIV)20,20,19                                                GELB1500
   19 PIV=TB                                                            GELB1510
      J=I                                                               GELB1520
      JJ=II                                                             GELB1530
   20 IF(I-MR)22,22,21                                                  GELB1540
   21 ID=ID-1                                                           GELB1550
   22 II=II+ID                                                          GELB1560
C                                                                       GELB1570
C     TEST ON SINGULARITY                                               GELB1580
      IF(PIV)47,47,23                                                   GELB1590
   23 IF(IER)26,24,26                                                   GELB1600
   24 IF(PIV-TOL)25,25,26                                               GELB1610
   25 IER=K-1                                                           GELB1620
   26 PIV=1./A(JJ)                                                      GELB1630
C                                                                       GELB1640
C     PIVOT ROW REDUCTION AND ROW INTERCHANGE IN RIGHT HAND SIDE R      GELB1650
      ID=J-K                                                            GELB1660
      DO 27 I=K,NM,M                                                    GELB1670
      II=I+ID                                                           GELB1680
      TB=PIV*R(II)                                                      GELB1690
      R(II)=R(I)                                                        GELB1700
   27 R(I)=TB                                                           GELB1710
C                                                                       GELB1720
C     PIVOT ROW REDUCTION AND ROW INTERCHANGE IN COEFFICIENT MATRIX A   GELB1730
      II=KST                                                            GELB1740
      J=JJ+IC                                                           GELB1750
      DO 28 I=JJ,J                                                      GELB1760
      TB=PIV*A(I)                                                       GELB1770
      A(I)=A(II)                                                        GELB1780
      A(II)=TB                                                          GELB1790
   28 II=II+1                                                           GELB1800
C                                                                       GELB1810
C     ELEMENT REDUCTION                                                 GELB1820
      IF(K-ILR)29,34,34                                                 GELB1830
   29 ID=KST                                                            GELB1840
      II=K+1                                                            GELB1850
      MU=KST+1                                                          GELB1860
      MZ=KST+IC                                                         GELB1870
      DO 33 I=II,ILR                                                    GELB1880
C                                                                       GELB1890
C     IN MATRIX A                                                       GELB1900
      ID=ID+MC                                                          GELB1910
      JJ=I-MR-1                                                         GELB1920
      IF(JJ)31,31,30                                                    GELB1930
   30 ID=ID-JJ                                                          GELB1940
   31 PIV=-A(ID)                                                        GELB1950
      J=ID+1                                                            GELB1960
      DO 32 JJ=MU,MZ                                                    GELB1970
      A(J-1)=A(J)+PIV*A(JJ)                                             GELB1980
   32 J=J+1                                                             GELB1990
      A(J-1)=0.                                                         GELB2000
C                                                                       GELB2010
C     IN MATRIX R                                                       GELB2020
      J=K                                                               GELB2030
      DO 33 JJ=I,NM,M                                                   GELB2040
      R(JJ)=R(JJ)+PIV*R(J)                                              GELB2050
   33 J=J+M                                                             GELB2060
   34 KST=KST+MC                                                        GELB2070
      IF(ILR-MR)36,35,35                                                GELB2080
   35 IC=IC-1                                                           GELB2090
   36 ID=K-MR                                                           GELB2100
      IF(ID)38,38,37                                                    GELB2110
   37 KST=KST-ID                                                        GELB2120
   38 CONTINUE                                                          GELB2130
C     END OF DECOMPOSITION LOOP                                         GELB2140
C                                                                       GELB2150
C                                                                       GELB2160
C     BACK SUBSTITUTION                                                 GELB2170
      IF(MC-1)46,46,39                                                  GELB2180
   39 IC=2                                                              GELB2190
      KST=MA+ML-MC+2                                                    GELB2200
      II=M                                                              GELB2210
      DO 45 I=2,M                                                       GELB2220
      KST=KST-MC                                                        GELB2230
      II=II-1                                                           GELB2240
      J=II-MR                                                           GELB2250
      IF(J)41,41,40                                                     GELB2260
   40 KST=KST+J                                                         GELB2270
   41 DO 43 J=II,NM,M                                                   GELB2280
      TB=R(J)                                                           GELB2290
      MZ=KST+IC-2                                                       GELB2300
      ID=J                                                              GELB2310
      DO 42 JJ=KST,MZ                                                   GELB2320
      ID=ID+1                                                           GELB2330
   42 TB=TB-A(JJ)*R(ID)                                                 GELB2340
   43 R(J)=TB                                                           GELB2350
      IF(IC-MC)44,45,45                                                 GELB2360
   44 IC=IC+1                                                           GELB2370
   45 CONTINUE                                                          GELB2380
   46 RETURN                                                            GELB2390
C                                                                       GELB2400
C                                                                       GELB2410
C     ERROR RETURN                                                      GELB2420
   47 IER=-1                                                            GELB2430
      RETURN                                                            GELB2440
      END                                                               GELB2450