Google
 

Trailing-Edge - PDP-10 Archives - decuslib10-02 - 43,50145/llsq.ssp
There are 2 other files named llsq.ssp in the archive. Click here to see a list.
C                                                                       LLSQ  10
C     ..................................................................LLSQ  20
C                                                                       LLSQ  30
C        SUBROUTINE LLSQ                                                LLSQ  40
C                                                                       LLSQ  50
C        PURPOSE                                                        LLSQ  60
C           TO SOLVE LINEAR LEAST SQUARES PROBLEMS, I.E. TO MINIMIZE    LLSQ  70
C           THE EUCLIDEAN NORM OF B-A*X, WHERE A IS A M BY N MATRIX     LLSQ  80
C           WITH M NOT LESS THAN N. IN THE SPECIAL CASE M=N SYSTEMS OF  LLSQ  90
C           LINEAR EQUATIONS MAY BE SOLVED.                             LLSQ 100
C                                                                       LLSQ 110
C        USAGE                                                          LLSQ 120
C           CALL LLSQ (A,B,M,N,L,X,IPIV,EPS,IER,AUX)                    LLSQ 130
C                                                                       LLSQ 140
C        DESCRIPTION OF PARAMETERS                                      LLSQ 150
C           A      - M BY N COEFFICIENT MATRIX (DESTROYED).             LLSQ 160
C           B      - M BY L RIGHT HAND SIDE MATRIX (DESTROYED).         LLSQ 170
C           M      - ROW NUMBER OF MATRICES A AND B.                    LLSQ 180
C           N      - COLUMN NUMBER OF MATRIX A, ROW NUMBER OF MATRIX X. LLSQ 190
C           L      - COLUMN NUMBER OF MATRICES B AND X.                 LLSQ 200
C           X      - N BY L SOLUTION MATRIX.                            LLSQ 210
C           IPIV   - INTEGER OUTPUT VECTOR OF DIMENSION N WHICH         LLSQ 220
C                    CONTAINS INFORMATIONS ON COLUMN INTERCHANGES       LLSQ 230
C                    IN MATRIX A. (SEE REMARK NO.3).                    LLSQ 240
C           EPS    - INPUT PARAMETER WHICH SPECIFIES A RELATIVE         LLSQ 250
C                    TOLERANCE FOR DETERMINATION OF RANK OF MATRIX A.   LLSQ 260
C           IER    - A RESULTING ERROR PARAMETER.                       LLSQ 270
C           AUX    - AUXILIARY STORAGE ARRAY OF DIMENSION MAX(2*N,L).   LLSQ 280
C                    ON RETURN FIRST L LOCATIONS OF AUX CONTAIN THE     LLSQ 290
C                    RESULTING LEAST SQUARES.                           LLSQ 300
C                                                                       LLSQ 310
C        REMARKS                                                        LLSQ 320
C           (1) NO ACTION BESIDES ERROR MESSAGE IER=-2 IN CASE          LLSQ 330
C               M LESS THAN N.                                          LLSQ 340
C           (2) NO ACTION BESIDES ERROR MESSAGE IER=-1 IN CASE          LLSQ 350
C               OF A ZERO-MATRIX A.                                     LLSQ 360
C           (3) IF RANK K OF MATRIX A IS FOUND TO BE LESS THAN N BUT    LLSQ 370
C               GREATER THAN 0, THE PROCEDURE RETURNS WITH ERROR CODE   LLSQ 380
C               IER=K INTO CALLING PROGRAM. THE LAST N-K ELEMENTS OF    LLSQ 390
C               VECTOR IPIV DENOTE THE USELESS COLUMNS IN MATRIX A.     LLSQ 400
C               THE REMAINING USEFUL COLUMNS FORM A BASE OF MATRIX A.   LLSQ 410
C           (4) IF THE PROCEDURE WAS SUCCESSFUL, ERROR PARAMETER IER    LLSQ 420
C               IS SET TO 0.                                            LLSQ 430
C                                                                       LLSQ 440
C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED                  LLSQ 450
C           NONE                                                        LLSQ 460
C                                                                       LLSQ 470
C        METHOD                                                         LLSQ 480
C           HOUSEHOLDER TRANSFORMATIONS ARE USED TO TRANSFORM MATRIX A  LLSQ 490
C           TO UPPER TRIANGULAR FORM. AFTER HAVING APPLIED THE SAME     LLSQ 500
C           TRANSFORMATION TO THE RIGHT HAND SIDE MATRIX B, AN          LLSQ 510
C           APPROXIMATE SOLUTION OF THE PROBLEM IS COMPUTED BY          LLSQ 520
C           BACK SUBSTITUTION. FOR REFERENCE, SEE                       LLSQ 530
C           G. GOLUB, NUMERICAL METHODS FOR SOLVING LINEAR LEAST        LLSQ 540
C           SQUARES PROBLEMS, NUMERISCHE MATHEMATIK, VOL.7,             LLSQ 550
C           ISS.3 (1965), PP.206-216.                                   LLSQ 560
C                                                                       LLSQ 570
C     ..................................................................LLSQ 580
C                                                                       LLSQ 590
      SUBROUTINE LLSQ(A,B,M,N,L,X,IPIV,EPS,IER,AUX)                     LLSQ 600
C                                                                       LLSQ 610
      DIMENSION A(1),B(1),X(1),IPIV(1),AUX(1)                           LLSQ 620
C                                                                       LLSQ 630
C     ERROR TEST                                                        LLSQ 640
      IF(M-N)30,1,1                                                     LLSQ 650
C                                                                       LLSQ 660
C     GENERATION OF INITIAL VECTOR S(K) (K=1,2,...,N) IN STORAGE        LLSQ 670
C     LOCATIONS AUX(K) (K=1,2,...,N)                                    LLSQ 680
    1 PIV=0.                                                            LLSQ 690
      IEND=0                                                            LLSQ 700
      DO 4 K=1,N                                                        LLSQ 710
      IPIV(K)=K                                                         LLSQ 720
      H=0.                                                              LLSQ 730
      IST=IEND+1                                                        LLSQ 740
      IEND=IEND+M                                                       LLSQ 750
      DO 2 I=IST,IEND                                                   LLSQ 760
    2 H=H+A(I)*A(I)                                                     LLSQ 770
      AUX(K)=H                                                          LLSQ 780
      IF(H-PIV)4,4,3                                                    LLSQ 790
    3 PIV=H                                                             LLSQ 800
      KPIV=K                                                            LLSQ 810
    4 CONTINUE                                                          LLSQ 820
C                                                                       LLSQ 830
C     ERROR TEST                                                        LLSQ 840
      IF(PIV)31,31,5                                                    LLSQ 850
C                                                                       LLSQ 860
C     DEFINE TOLERANCE FOR CHECKING RANK OF A                           LLSQ 870
    5 SIG=SQRT(PIV)                                                     LLSQ 880
      TOL=SIG*ABS(EPS)                                                  LLSQ 890
C                                                                       LLSQ 900
C                                                                       LLSQ 910
C     DECOMPOSITION LOOP                                                LLSQ 920
      LM=L*M                                                            LLSQ 930
      IST=-M                                                            LLSQ 940
      DO 21 K=1,N                                                       LLSQ 950
      IST=IST+M+1                                                       LLSQ 960
      IEND=IST+M-K                                                      LLSQ 970
      I=KPIV-K                                                          LLSQ 980
      IF(I)8,8,6                                                        LLSQ 990
C                                                                       LLSQ1000
C     INTERCHANGE K-TH COLUMN OF A WITH KPIV-TH IN CASE KPIV.GT.K       LLSQ1010
    6 H=AUX(K)                                                          LLSQ1020
      AUX(K)=AUX(KPIV)                                                  LLSQ1030
      AUX(KPIV)=H                                                       LLSQ1040
      ID=I*M                                                            LLSQ1050
      DO 7 I=IST,IEND                                                   LLSQ1060
      J=I+ID                                                            LLSQ1070
      H=A(I)                                                            LLSQ1080
      A(I)=A(J)                                                         LLSQ1090
    7 A(J)=H                                                            LLSQ1100
C                                                                       LLSQ1110
C     COMPUTATION OF PARAMETER SIG                                      LLSQ1120
    8 IF(K-1)11,11,9                                                    LLSQ1130
    9 SIG=0.                                                            LLSQ1140
      DO 10 I=IST,IEND                                                  LLSQ1150
   10 SIG=SIG+A(I)*A(I)                                                 LLSQ1160
      SIG=SQRT(SIG)                                                     LLSQ1170
C                                                                       LLSQ1180
C     TEST ON SINGULARITY                                               LLSQ1190
      IF(SIG-TOL)32,32,11                                               LLSQ1200
C                                                                       LLSQ1210
C     GENERATE CORRECT SIGN OF PARAMETER SIG                            LLSQ1220
   11 H=A(IST)                                                          LLSQ1230
      IF(H)12,13,13                                                     LLSQ1240
   12 SIG=-SIG                                                          LLSQ1250
C                                                                       LLSQ1260
C     SAVE INTERCHANGE INFORMATION                                      LLSQ1270
   13 IPIV(KPIV)=IPIV(K)                                                LLSQ1280
      IPIV(K)=KPIV                                                      LLSQ1290
C                                                                       LLSQ1300
C     GENERATION OF VECTOR UK IN K-TH COLUMN OF MATRIX A AND OF         LLSQ1310
C     PARAMETER BETA                                                    LLSQ1320
      BETA=H+SIG                                                        LLSQ1330
      A(IST)=BETA                                                       LLSQ1340
      BETA=1./(SIG*BETA)                                                LLSQ1350
      J=N+K                                                             LLSQ1360
      AUX(J)=-SIG                                                       LLSQ1370
      IF(K-N)14,19,19                                                   LLSQ1380
C                                                                       LLSQ1390
C     TRANSFORMATION OF MATRIX A                                        LLSQ1400
   14 PIV=0.                                                            LLSQ1410
      ID=0                                                              LLSQ1420
      JST=K+1                                                           LLSQ1430
      KPIV=JST                                                          LLSQ1440
      DO 18 J=JST,N                                                     LLSQ1450
      ID=ID+M                                                           LLSQ1460
      H=0.                                                              LLSQ1470
      DO 15 I=IST,IEND                                                  LLSQ1480
      II=I+ID                                                           LLSQ1490
   15 H=H+A(I)*A(II)                                                    LLSQ1500
      H=BETA*H                                                          LLSQ1510
      DO 16 I=IST,IEND                                                  LLSQ1520
      II=I+ID                                                           LLSQ1530
   16 A(II)=A(II)-A(I)*H                                                LLSQ1540
C                                                                       LLSQ1550
C     UPDATING OF ELEMENT S(J) STORED IN LOCATION AUX(J)                LLSQ1560
      II=IST+ID                                                         LLSQ1570
      H=AUX(J)-A(II)*A(II)                                              LLSQ1580
      AUX(J)=H                                                          LLSQ1590
      IF(H-PIV)18,18,17                                                 LLSQ1600
   17 PIV=H                                                             LLSQ1610
      KPIV=J                                                            LLSQ1620
   18 CONTINUE                                                          LLSQ1630
C                                                                       LLSQ1640
C     TRANSFORMATION OF RIGHT HAND SIDE MATRIX B                        LLSQ1650
   19 DO 21 J=K,LM,M                                                    LLSQ1660
      H=0.                                                              LLSQ1670
      IEND=J+M-K                                                        LLSQ1680
      II=IST                                                            LLSQ1690
      DO 20 I=J,IEND                                                    LLSQ1700
      H=H+A(II)*B(I)                                                    LLSQ1710
   20 II=II+1                                                           LLSQ1720
      H=BETA*H                                                          LLSQ1730
      II=IST                                                            LLSQ1740
      DO 21 I=J,IEND                                                    LLSQ1750
      B(I)=B(I)-A(II)*H                                                 LLSQ1760
   21 II=II+1                                                           LLSQ1770
C     END OF DECOMPOSITION LOOP                                         LLSQ1780
C                                                                       LLSQ1790
C                                                                       LLSQ1800
C     BACK SUBSTITUTION AND BACK INTERCHANGE                            LLSQ1810
      IER=0                                                             LLSQ1820
      I=N                                                               LLSQ1830
      LN=L*N                                                            LLSQ1840
      PIV=1./AUX(2*N)                                                   LLSQ1850
      DO 22 K=N,LN,N                                                    LLSQ1860
      X(K)=PIV*B(I)                                                     LLSQ1870
   22 I=I+M                                                             LLSQ1880
      IF(N-1)26,26,23                                                   LLSQ1890
   23 JST=(N-1)*M+N                                                     LLSQ1900
      DO 25 J=2,N                                                       LLSQ1910
      JST=JST-M-1                                                       LLSQ1920
      K=N+N+1-J                                                         LLSQ1930
      PIV=1./AUX(K)                                                     LLSQ1940
      KST=K-N                                                           LLSQ1950
      ID=IPIV(KST)-KST                                                  LLSQ1960
      IST=2-J                                                           LLSQ1970
      DO 25 K=1,L                                                       LLSQ1980
      H=B(KST)                                                          LLSQ1990
      IST=IST+N                                                         LLSQ2000
      IEND=IST+J-2                                                      LLSQ2010
      II=JST                                                            LLSQ2020
      DO 24 I=IST,IEND                                                  LLSQ2030
      II=II+M                                                           LLSQ2040
   24 H=H-A(II)*X(I)                                                    LLSQ2050
      I=IST-1                                                           LLSQ2060
      II=I+ID                                                           LLSQ2070
      X(I)=X(II)                                                        LLSQ2080
      X(II)=PIV*H                                                       LLSQ2090
   25 KST=KST+M                                                         LLSQ2100
C                                                                       LLSQ2110
C                                                                       LLSQ2120
C     COMPUTATION OF LEAST SQUARES                                      LLSQ2130
   26 IST=N+1                                                           LLSQ2140
      IEND=0                                                            LLSQ2150
      DO 29 J=1,L                                                       LLSQ2160
      IEND=IEND+M                                                       LLSQ2170
      H=0.                                                              LLSQ2180
      IF(M-N)29,29,27                                                   LLSQ2190
   27 DO 28 I=IST,IEND                                                  LLSQ2200
   28 H=H+B(I)*B(I)                                                     LLSQ2210
      IST=IST+M                                                         LLSQ2220
   29 AUX(J)=H                                                          LLSQ2230
      RETURN                                                            LLSQ2240
C                                                                       LLSQ2250
C     ERROR RETURN IN CASE M LESS THAN N                                LLSQ2260
   30 IER=-2                                                            LLSQ2270
      RETURN                                                            LLSQ2280
C                                                                       LLSQ2290
C     ERROR RETURN IN CASE OF ZERO-MATRIX A                             LLSQ2300
   31 IER=-1                                                            LLSQ2310
      RETURN                                                            LLSQ2320
C                                                                       LLSQ2330
C     ERROR RETURN IN CASE OF RANK OF MATRIX A LESS THAN N              LLSQ2340
   32 IER=K-1                                                           LLSQ2350
      RETURN                                                            LLSQ2360
      END                                                               LLSQ2370