Google
 

Trailing-Edge - PDP-10 Archives - decuslib20-02 - decus/20-0026/dllsq.ssp
There are 2 other files named dllsq.ssp in the archive. Click here to see a list.
C                                                                       DLLS  10
C     ..................................................................DLLS  20
C                                                                       DLLS  30
C        SUBROUTINE DLLSQ                                               DLLS  40
C                                                                       DLLS  50
C        PURPOSE                                                        DLLS  60
C           TO SOLVE LINEAR LEAST SQUARES PROBLEMS, I.E. TO MINIMIZE    DLLS  70
C           THE EUCLIDEAN NORM OF B-A*X, WHERE A IS A M BY N MATRIX     DLLS  80
C           WITH M NOT LESS THAN N. IN THE SPECIAL CASE M=N SYSTEMS OF  DLLS  90
C           LINEAR EQUATIONS MAY BE SOLVED.                             DLLS 100
C                                                                       DLLS 110
C        USAGE                                                          DLLS 120
C           CALL DLLSQ (A,B,M,N,L,X,IPIV,EPS,IER,AUX)                   DLLS 130
C                                                                       DLLS 140
C        DESCRIPTION OF PARAMETERS                                      DLLS 150
C           A      - DOUBLE PRECISION M BY N COEFFICIENT MATRIX         DLLS 160
C                    (DESTROYED).                                       DLLS 170
C           B      - DOUBLE PRECISION M BY L RIGHT HAND SIDE MATRIX     DLLS 180
C                    (DESTROYED).                                       DLLS 190
C           M      - ROW NUMBER OF MATRICES A AND B.                    DLLS 200
C           N      - COLUMN NUMBER OF MATRIX A, ROW NUMBER OF MATRIX X. DLLS 210
C           L      - COLUMN NUMBER OF MATRICES B AND X.                 DLLS 220
C           X      - DOUBLE PRECISION N BY L SOLUTION MATRIX.           DLLS 230
C           IPIV   - INTEGER OUTPUT VECTOR OF DIMENSION N WHICH         DLLS 240
C                    CONTAINS INFORMATIONS ON COLUMN INTERCHANGES       DLLS 250
C                    IN MATRIX A. (SEE REMARK NO.3).                    DLLS 260
C           EPS    - SINGLE PRECISION INPUT PARAMETER WHICH SPECIFIES   DLLS 270
C                    A RELATIVE TOLERANCE FOR DETERMINATION OF RANK OF  DLLS 280
C                    MATRIX A.                                          DLLS 290
C           IER    - A RESULTING ERROR PARAMETER.                       DLLS 300
C           AUX    - A DOUBLE PRECISION AUXILIARY STORAGE ARRAY OF      DLLS 310
C                    DIMENSION MAX(2*N,L). ON RETURN FIRST L LOCATIONS  DLLS 320
C                    OF AUX CONTAIN THE RESULTING LEAST SQUARES.        DLLS 330
C                                                                       DLLS 340
C        REMARKS                                                        DLLS 350
C           (1) NO ACTION BESIDES ERROR MESSAGE IER=-2 IN CASE          DLLS 360
C               M LESS THAN N.                                          DLLS 370
C           (2) NO ACTION BESIDES ERROR MESSAGE IER=-1 IN CASE          DLLS 380
C               OF A ZERO-MATRIX A.                                     DLLS 390
C           (3) IF RANK K OF MATRIX A IS FOUND TO BE LESS THAN N BUT    DLLS 400
C               GREATER THAN 0, THE PROCEDURE RETURNS WITH ERROR CODE   DLLS 410
C               IER=K INTO CALLING PROGRAM. THE LAST N-K ELEMENTS OF    DLLS 420
C               VECTOR IPIV DENOTE THE USELESS COLUMNS IN MATRIX A.     DLLS 430
C               THE REMAINING USEFUL COLUMNS FORM A BASE OF MATRIX A.   DLLS 440
C           (4) IF THE PROCEDURE WAS SUCCESSFUL, ERROR PARAMETER IER    DLLS 450
C               IS SET TO 0.                                            DLLS 460
C                                                                       DLLS 470
C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED                  DLLS 480
C           NONE                                                        DLLS 490
C                                                                       DLLS 500
C        METHOD                                                         DLLS 510
C           HOUSEHOLDER TRANSFORMATIONS ARE USED TO TRANSFORM MATRIX A  DLLS 520
C           TO UPPER TRIANGULAR FORM. AFTER HAVING APPLIED THE SAME     DLLS 530
C           TRANSFORMATION TO THE RIGHT HAND SIDE MATRIX B, AN          DLLS 540
C           APPROXIMATE SOLUTION OF THE PROBLEM IS COMPUTED BY          DLLS 550
C           BACK SUBSTITUTION. FOR REFERENCE, SEE                       DLLS 560
C           G. GOLUB, NUMERICAL METHODS FOR SOLVING LINEAR LEAST        DLLS 570
C           SQUARES PROBLEMS, NUMERISCHE MATHEMATIK, VOL.7,             DLLS 580
C           ISS.3 (1965), PP.206-216.                                   DLLS 590
C                                                                       DLLS 600
C     ..................................................................DLLS 610
C                                                                       DLLS 620
      SUBROUTINE DLLSQ(A,B,M,N,L,X,IPIV,EPS,IER,AUX)                    DLLS 630
C                                                                       DLLS 640
      DIMENSION A(1),B(1),X(1),IPIV(1),AUX(1)                           DLLS 650
      DOUBLE PRECISION A,B,X,AUX,PIV,H,SIG,BETA,TOL                     DLLS 660
C                                                                       DLLS 670
C     ERROR TEST                                                        DLLS 680
      IF(M-N)30,1,1                                                     DLLS 690
C                                                                       DLLS 700
C     GENERATION OF INITIAL VECTOR S(K) (K=1,2,...,N) IN STORAGE        DLLS 710
C     LOCATIONS AUX(K) (K=1,2,...,N)                                    DLLS 720
    1 PIV=0.D0                                                          DLLS 730
      IEND=0                                                            DLLS 740
      DO 4 K=1,N                                                        DLLS 750
      IPIV(K)=K                                                         DLLS 760
      H=0.D0                                                            DLLS 770
      IST=IEND+1                                                        DLLS 780
      IEND=IEND+M                                                       DLLS 790
      DO 2 I=IST,IEND                                                   DLLS 800
    2 H=H+A(I)*A(I)                                                     DLLS 810
      AUX(K)=H                                                          DLLS 820
      IF(H-PIV)4,4,3                                                    DLLS 830
    3 PIV=H                                                             DLLS 840
      KPIV=K                                                            DLLS 850
    4 CONTINUE                                                          DLLS 860
C                                                                       DLLS 870
C     ERROR TEST                                                        DLLS 880
      IF(PIV)31,31,5                                                    DLLS 890
C                                                                       DLLS 900
C     DEFINE TOLERANCE FOR CHECKING RANK OF A                           DLLS 910
    5 SIG=DSQRT(PIV)                                                    DLLS 920
      TOL=SIG*ABS(EPS)                                                  DLLS 930
C                                                                       DLLS 940
C                                                                       DLLS 950
C     DECOMPOSITION LOOP                                                DLLS 960
      LM=L*M                                                            DLLS 970
      IST=-M                                                            DLLS 980
      DO 21 K=1,N                                                       DLLS 990
      IST=IST+M+1                                                       DLLS1000
      IEND=IST+M-K                                                      DLLS1010
      I=KPIV-K                                                          DLLS1020
      IF(I)8,8,6                                                        DLLS1030
C                                                                       DLLS1040
C     INTERCHANGE K-TH COLUMN OF A WITH KPIV-TH IN CASE KPIV.GT.K       DLLS1050
    6 H=AUX(K)                                                          DLLS1060
      AUX(K)=AUX(KPIV)                                                  DLLS1070
      AUX(KPIV)=H                                                       DLLS1080
      ID=I*M                                                            DLLS1090
      DO 7 I=IST,IEND                                                   DLLS1100
      J=I+ID                                                            DLLS1110
      H=A(I)                                                            DLLS1120
      A(I)=A(J)                                                         DLLS1130
    7 A(J)=H                                                            DLLS1140
C                                                                       DLLS1150
C     COMPUTATION OF PARAMETER SIG                                      DLLS1160
    8 IF(K-1)11,11,9                                                    DLLS1170
    9 SIG=0.D0                                                          DLLS1180
      DO 10 I=IST,IEND                                                  DLLS1190
   10 SIG=SIG+A(I)*A(I)                                                 DLLS1200
      SIG=DSQRT(SIG)                                                    DLLS1210
C                                                                       DLLS1220
C     TEST ON SINGULARITY                                               DLLS1230
      IF(SIG-TOL)32,32,11                                               DLLS1240
C                                                                       DLLS1250
C     GENERATE CORRECT SIGN OF PARAMETER SIG                            DLLS1260
   11 H=A(IST)                                                          DLLS1270
      IF(H)12,13,13                                                     DLLS1280
   12 SIG=-SIG                                                          DLLS1290
C                                                                       DLLS1300
C     SAVE INTERCHANGE INFORMATION                                      DLLS1310
   13 IPIV(KPIV)=IPIV(K)                                                DLLS1320
      IPIV(K)=KPIV                                                      DLLS1330
C                                                                       DLLS1340
C     GENERATION OF VECTOR UK IN K-TH COLUMN OF MATRIX A AND OF         DLLS1350
C     PARAMETER BETA                                                    DLLS1360
      BETA=H+SIG                                                        DLLS1370
      A(IST)=BETA                                                       DLLS1380
      BETA=1.D0/(SIG*BETA)                                              DLLS1390
      J=N+K                                                             DLLS1400
      AUX(J)=-SIG                                                       DLLS1410
      IF(K-N)14,19,19                                                   DLLS1420
C                                                                       DLLS1430
C     TRANSFORMATION OF MATRIX A                                        DLLS1440
   14 PIV=0.D0                                                          DLLS1450
      ID=0                                                              DLLS1460
      JST=K+1                                                           DLLS1470
      KPIV=JST                                                          DLLS1480
      DO 18 J=JST,N                                                     DLLS1490
      ID=ID+M                                                           DLLS1500
      H=0.D0                                                            DLLS1510
      DO 15 I=IST,IEND                                                  DLLS1520
      II=I+ID                                                           DLLS1530
   15 H=H+A(I)*A(II)                                                    DLLS1540
      H=BETA*H                                                          DLLS1550
      DO 16 I=IST,IEND                                                  DLLS1560
      II=I+ID                                                           DLLS1570
   16 A(II)=A(II)-A(I)*H                                                DLLS1580
C                                                                       DLLS1590
C     UPDATING OF ELEMENT S(J) STORED IN LOCATION AUX(J)                DLLS1600
      II=IST+ID                                                         DLLS1610
      H=AUX(J)-A(II)*A(II)                                              DLLS1620
      AUX(J)=H                                                          DLLS1630
      IF(H-PIV)18,18,17                                                 DLLS1640
   17 PIV=H                                                             DLLS1650
      KPIV=J                                                            DLLS1660
   18 CONTINUE                                                          DLLS1670
C                                                                       DLLS1680
C     TRANSFORMATION OF RIGHT HAND SIDE MATRIX B                        DLLS1690
   19 DO 21 J=K,LM,M                                                    DLLS1700
      H=0.D0                                                            DLLS1710
      IEND=J+M-K                                                        DLLS1720
      II=IST                                                            DLLS1730
      DO 20 I=J,IEND                                                    DLLS1740
      H=H+A(II)*B(I)                                                    DLLS1750
   20 II=II+1                                                           DLLS1760
      H=BETA*H                                                          DLLS1770
      II=IST                                                            DLLS1780
      DO 21 I=J,IEND                                                    DLLS1790
      B(I)=B(I)-A(II)*H                                                 DLLS1800
   21 II=II+1                                                           DLLS1810
C     END OF DECOMPOSITION LOOP                                         DLLS1820
C                                                                       DLLS1830
C                                                                       DLLS1840
C     BACK SUBSTITUTION AND BACK INTERCHANGE                            DLLS1850
      IER=0                                                             DLLS1860
      I=N                                                               DLLS1870
      LN=L*N                                                            DLLS1880
      PIV=1.D0/AUX(2*N)                                                 DLLS1890
      DO 22 K=N,LN,N                                                    DLLS1900
      X(K)=PIV*B(I)                                                     DLLS1910
   22 I=I+M                                                             DLLS1920
      IF(N-1)26,26,23                                                   DLLS1930
   23 JST=(N-1)*M+N                                                     DLLS1940
      DO 25 J=2,N                                                       DLLS1950
      JST=JST-M-1                                                       DLLS1960
      K=N+N+1-J                                                         DLLS1970
      PIV=1.D0/AUX(K)                                                   DLLS1980
      KST=K-N                                                           DLLS1990
      ID=IPIV(KST)-KST                                                  DLLS2000
      IST=2-J                                                           DLLS2010
      DO 25 K=1,L                                                       DLLS2020
      H=B(KST)                                                          DLLS2030
      IST=IST+N                                                         DLLS2040
      IEND=IST+J-2                                                      DLLS2050
      II=JST                                                            DLLS2060
      DO 24 I=IST,IEND                                                  DLLS2070
      II=II+M                                                           DLLS2080
   24 H=H-A(II)*X(I)                                                    DLLS2090
      I=IST-1                                                           DLLS2100
      II=I+ID                                                           DLLS2110
      X(I)=X(II)                                                        DLLS2120
      X(II)=PIV*H                                                       DLLS2130
   25 KST=KST+M                                                         DLLS2140
C                                                                       DLLS2150
C                                                                       DLLS2160
C     COMPUTATION OF LEAST SQUARES                                      DLLS2170
   26 IST=N+1                                                           DLLS2180
      IEND=0                                                            DLLS2190
      DO 29 J=1,L                                                       DLLS2200
      IEND=IEND+M                                                       DLLS2210
      H=0.D0                                                            DLLS2220
      IF(M-N)29,29,27                                                   DLLS2230
   27 DO 28 I=IST,IEND                                                  DLLS2240
   28 H=H+B(I)*B(I)                                                     DLLS2250
      IST=IST+M                                                         DLLS2260
   29 AUX(J)=H                                                          DLLS2270
      RETURN                                                            DLLS2280
C                                                                       DLLS2290
C     ERROR RETURN IN CASE M LESS THAN N                                DLLS2300
   30 IER=-2                                                            DLLS2310
      RETURN                                                            DLLS2320
C                                                                       DLLS2330
C     ERROR RETURN IN CASE OF ZERO-MATRIX A                             DLLS2340
   31 IER=-1                                                            DLLS2350
      RETURN                                                            DLLS2360
C                                                                       DLLS2370
C     ERROR RETURN IN CASE OF RANK OF MATRIX A LESS THAN N              DLLS2380
   32 IER=K-1                                                           DLLS2390
      RETURN                                                            DLLS2400
      END                                                               DLLS2410