Web pdp-10.trailing-edge.com

Trailing-Edge - PDP-10 Archives - decuslib20-02 - decus/20-0026/dgels.ssp
There are 2 other files named dgels.ssp in the archive. Click here to see a list.
```C                                                                       DELS  10
C     ..................................................................DELS  20
C                                                                       DELS  30
C        SUBROUTINE DGELS                                               DELS  40
C                                                                       DELS  50
C        PURPOSE                                                        DELS  60
C           TO SOLVE A SYSTEM OF SIMULTANEOUS LINEAR EQUATIONS WITH     DELS  70
C           SYMMETRIC COEFFICIENT MATRIX UPPER TRIANGULAR PART OF WHICH DELS  80
C           IS ASSUMED TO BE STORED COLUMNWISE.                         DELS  90
C                                                                       DELS 100
C        USAGE                                                          DELS 110
C           CALL DGELS(R,A,M,N,EPS,IER,AUX)                             DELS 120
C                                                                       DELS 130
C        DESCRIPTION OF PARAMETERS                                      DELS 140
C           R      - DOUBLE PRECISION M BY N RIGHT HAND SIDE MATRIX     DELS 150
C                    (DESTROYED). ON RETURN R CONTAINS THE SOLUTION OF  DELS 160
C                    THE EQUATIONS.                                     DELS 170
C           A      - UPPER TRIANGULAR PART OF THE SYMMETRIC DOUBLE      DELS 180
C                    PRECISION M BY M COEFFICIENT MATRIX.  (DESTROYED)  DELS 190
C           M      - THE NUMBER OF EQUATIONS IN THE SYSTEM.             DELS 200
C           N      - THE NUMBER OF RIGHT HAND SIDE VECTORS.             DELS 210
C           EPS    - SINGLE PRECISION INPUT CONSTANT WHICH IS USED AS   DELS 220
C                    RELATIVE TOLERANCE FOR TEST ON LOSS OF             DELS 230
C                    SIGNIFICANCE.                                      DELS 240
C           IER    - RESULTING ERROR PARAMETER CODED AS FOLLOWS         DELS 250
C                    IER=0  - NO ERROR,                                 DELS 260
C                    IER=-1 - NO RESULT BECAUSE OF M LESS THAN 1 OR     DELS 270
C                             PIVOT ELEMENT AT ANY ELIMINATION STEP     DELS 280
C                             EQUAL TO 0,                               DELS 290
C                    IER=K  - WARNING DUE TO POSSIBLE LOSS OF SIGNIFI-  DELS 300
C                             CANCE INDICATED AT ELIMINATION STEP K+1,  DELS 310
C                             WHERE PIVOT ELEMENT WAS LESS THAN OR      DELS 320
C                             EQUAL TO THE INTERNAL TOLERANCE EPS TIMES DELS 330
C                             ABSOLUTELY GREATEST MAIN DIAGONAL         DELS 340
C                             ELEMENT OF MATRIX A.                      DELS 350
C           AUX    - DOUBLE PRECISION AUXILIARY STORAGE ARRAY           DELS 360
C                    WITH DIMENSION M-1.                                DELS 370
C                                                                       DELS 380
C        REMARKS                                                        DELS 390
C           UPPER TRIANGULAR PART OF MATRIX A IS ASSUMED TO BE STORED   DELS 400
C           COLUMNWISE IN M*(M+1)/2 SUCCESSIVE STORAGE LOCATIONS, RIGHT DELS 410
C           HAND SIDE MATRIX R COLUMNWISE IN N*M SUCCESSIVE STORAGE     DELS 420
C           LOCATIONS. ON RETURN SOLUTION MATRIX R IS STORED COLUMNWISE DELS 430
C           TOO.                                                        DELS 440
C           THE PROCEDURE GIVES RESULTS IF THE NUMBER OF EQUATIONS M IS DELS 450
C           GREATER THAN 0 AND PIVOT ELEMENTS AT ALL ELIMINATION STEPS  DELS 460
C           ARE DIFFERENT FROM 0. HOWEVER WARNING IER=K - IF GIVEN -    DELS 470
C           INDICATES POSSIBLE LOSS OF SIGNIFICANCE. IN CASE OF A WELL  DELS 480
C           SCALED MATRIX A AND APPROPRIATE TOLERANCE EPS, IER=K MAY BE DELS 490
C           INTERPRETED THAT MATRIX A HAS THE RANK K. NO WARNING IS     DELS 500
C           GIVEN IN CASE M=1.                                          DELS 510
C           ERROR PARAMETER IER=-1 DOES NOT NECESSARILY MEAN THAT       DELS 520
C           MATRIX A IS SINGULAR, AS ONLY MAIN DIAGONAL ELEMENTS        DELS 530
C           ARE USED AS PIVOT ELEMENTS. POSSIBLY SUBROUTINE DGELG (WHICHDELS 540
C           WORKS WITH TOTAL PIVOTING) WOULD BE ABLE TO FIND A SOLUTION.DELS 550
C                                                                       DELS 560
C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED                  DELS 570
C           NONE                                                        DELS 580
C                                                                       DELS 590
C        METHOD                                                         DELS 600
C           SOLUTION IS DONE BY MEANS OF GAUSS-ELIMINATION WITH         DELS 610
C           PIVOTING IN MAIN DIAGONAL, IN ORDER TO PRESERVE             DELS 620
C           SYMMETRY IN REMAINING COEFFICIENT MATRICES.                 DELS 630
C                                                                       DELS 640
C     ..................................................................DELS 650
C                                                                       DELS 660
SUBROUTINE DGELS(R,A,M,N,EPS,IER,AUX)                             DELS 670
C                                                                       DELS 680
C                                                                       DELS 690
DIMENSION A(1),R(1),AUX(1)                                        DELS 700
DOUBLE PRECISION R,A,AUX,PIV,TB,TOL,PIVI                          DELS 710
IF(M)24,24,1                                                      DELS 720
C                                                                       DELS 730
C     SEARCH FOR GREATEST MAIN DIAGONAL ELEMENT                         DELS 740
1 IER=0                                                             DELS 750
PIV=0.D0                                                          DELS 760
L=0                                                               DELS 770
DO 3 K=1,M                                                        DELS 780
L=L+K                                                             DELS 790
TB=DABS(A(L))                                                     DELS 800
IF(TB-PIV)3,3,2                                                   DELS 810
2 PIV=TB                                                            DELS 820
I=L                                                               DELS 830
J=K                                                               DELS 840
3 CONTINUE                                                          DELS 850
TOL=EPS*PIV                                                       DELS 860
C     MAIN DIAGONAL ELEMENT A(I)=A(J,J) IS FIRST PIVOT ELEMENT.         DELS 870
C     PIV CONTAINS THE ABSOLUTE VALUE OF A(I).                          DELS 880
C                                                                       DELS 890
C                                                                       DELS 900
C     START ELIMINATION LOOP                                            DELS 910
LST=0                                                             DELS 920
NM=N*M                                                            DELS 930
LEND=M-1                                                          DELS 940
DO 18 K=1,M                                                       DELS 950
C                                                                       DELS 960
C     TEST ON USEFULNESS OF SYMMETRIC ALGORITHM                         DELS 970
IF(PIV)24,24,4                                                    DELS 980
4 IF(IER)7,5,7                                                      DELS 990
5 IF(PIV-TOL)6,6,7                                                  DELS1000
6 IER=K-1                                                           DELS1010
7 LT=J-K                                                            DELS1020
LST=LST+K                                                         DELS1030
C                                                                       DELS1040
C     PIVOT ROW REDUCTION AND ROW INTERCHANGE IN RIGHT HAND SIDE R      DELS1050
PIVI=1.D0/A(I)                                                    DELS1060
DO 8 L=K,NM,M                                                     DELS1070
LL=L+LT                                                           DELS1080
TB=PIVI*R(LL)                                                     DELS1090
R(LL)=R(L)                                                        DELS1100
8 R(L)=TB                                                           DELS1110
C                                                                       DELS1120
C     IS ELIMINATION TERMINATED                                         DELS1130
IF(K-M)9,19,19                                                    DELS1140
C                                                                       DELS1150
C     ROW AND COLUMN INTERCHANGE AND PIVOT ROW REDUCTION IN MATRIX A.   DELS1160
C     ELEMENTS OF PIVOT COLUMN ARE SAVED IN AUXILIARY VECTOR AUX.       DELS1170
9 LR=LST+(LT*(K+J-1))/2                                             DELS1180
LL=LR                                                             DELS1190
L=LST                                                             DELS1200
DO 14 II=K,LEND                                                   DELS1210
L=L+II                                                            DELS1220
LL=LL+1                                                           DELS1230
IF(L-LR)12,10,11                                                  DELS1240
10 A(LL)=A(LST)                                                      DELS1250
TB=A(L)                                                           DELS1260
GO TO 13                                                          DELS1270
11 LL=L+LT                                                           DELS1280
12 TB=A(LL)                                                          DELS1290
A(LL)=A(L)                                                        DELS1300
13 AUX(II)=TB                                                        DELS1310
14 A(L)=PIVI*TB                                                      DELS1320
C                                                                       DELS1330
C     SAVE COLUMN INTERCHANGE INFORMATION                               DELS1340
A(LST)=LT                                                         DELS1350
C                                                                       DELS1360
C     ELEMENT REDUCTION AND SEARCH FOR NEXT PIVOT                       DELS1370
PIV=0.D0                                                          DELS1380
LLST=LST                                                          DELS1390
LT=0                                                              DELS1400
DO 18 II=K,LEND                                                   DELS1410
PIVI=-AUX(II)                                                     DELS1420
LL=LLST                                                           DELS1430
LT=LT+1                                                           DELS1440
DO 15 LLD=II,LEND                                                 DELS1450
LL=LL+LLD                                                         DELS1460
L=LL+LT                                                           DELS1470
15 A(L)=A(L)+PIVI*A(LL)                                              DELS1480
LLST=LLST+II                                                      DELS1490
LR=LLST+LT                                                        DELS1500
TB=DABS(A(LR))                                                    DELS1510
IF(TB-PIV)17,17,16                                                DELS1520
16 PIV=TB                                                            DELS1530
I=LR                                                              DELS1540
J=II+1                                                            DELS1550
17 DO 18 LR=K,NM,M                                                   DELS1560
LL=LR+LT                                                          DELS1570
18 R(LL)=R(LL)+PIVI*R(LR)                                            DELS1580
C     END OF ELIMINATION LOOP                                           DELS1590
C                                                                       DELS1600
C                                                                       DELS1610
C     BACK SUBSTITUTION AND BACK INTERCHANGE                            DELS1620
19 IF(LEND)24,23,20                                                  DELS1630
20 II=M                                                              DELS1640
DO 22 I=2,M                                                       DELS1650
LST=LST-II                                                        DELS1660
II=II-1                                                           DELS1670
L=A(LST)+.5D0                                                     DELS1680
DO 22 J=II,NM,M                                                   DELS1690
TB=R(J)                                                           DELS1700
LL=J                                                              DELS1710
K=LST                                                             DELS1720
DO 21 LT=II,LEND                                                  DELS1730
LL=LL+1                                                           DELS1740
K=K+LT                                                            DELS1750
21 TB=TB-A(K)*R(LL)                                                  DELS1760
K=J+L                                                             DELS1770
R(J)=R(K)                                                         DELS1780
22 R(K)=TB                                                           DELS1790
23 RETURN                                                            DELS1800
C                                                                       DELS1810
C                                                                       DELS1820
C     ERROR RETURN                                                      DELS1830
24 IER=-1                                                            DELS1840
RETURN                                                            DELS1850
END                                                               DELS1860

```