         TITLE    'MAT-B00,08/22/73,DWG702985'
         PAGE
*
*
*  E X T E R N A L    C O M M U N I C A T I O N
*
*
*  DEFINITIONS
*
         DEF      DMATDIV           DYADIC MATRIX DIVIDE OP ROUTINE
         DEF      MMATINV           MONADIC MATRIX INVERT OP ROUTINE
         DEF      MAT@              START OF PROCEDURE
*
*  REFERENCES
*
         REF      ALOCHNW           ALLOCATE HEADER AND N WORDS
         REF      ALOCRS            ALLOCATE RESULT DATA BLOCK
         REF      DREF              DE-REF
         REF      DTYPEF            DYADIC TYPE CHECK: FLOT RESULT
         REF      DXRETURN          DYADIC EXECUTION DRIVER RETURN
         REF      ERLENGTH          LENGTH ERROR
         REF      ERRANK            RANK ERROR
         REF      ERSING            SINGULAR MATRIX ERROR EXIT
         REF      EXECUTE           EXECUTE XSEG
         REF      FLOT0             FLOATING POINT 0.0
         REF      FLOT1             FLOATING POINT 1.0
         REF      FSQRT             F SQUARE ROOT EVAL
         REF      GENLOAD           GEN LOAD BY RSTYPE
         REF      GIVEBACK          SHRINK DATA BLOCK
         REF      GXSEGINI          GEN XSEG INITIALIZATION
         REF      LFARG             LEFT ARG PNTR
         REF      LOOPLOC           LOOP LOCATION
         REF      MTYPEF            MONADIC TYPE CHECK: FLOT RESULT
         REF      MXRETURN          MONADIC EXECUTION DRIVER RETURN
         REF      OPBREAK           OP BREAK HANDLER
         REF      RESULT            RESULT PNTR
         REF      RETURN            RETURN ADR CELL
         REF      RSADR             RESULT ADDRESS
         REF      RSRANK            RESULT RANK
         REF      RSSIZE            RESULT SIZE
         REF      RTARG             RIGHT ARG PNTR
         REF      SETADR            SET UP ARG ADR CELL
         REF      XSEGBASE          XSEG BASE
         REF      XSEGBRK           XSEG BREAK FLAG
         REF      XSETUP            SET ADR FOR INDEXED ACCESS
         PAGE
*
*
*  A S S E M B L Y    P A R A M E T E R S
*
*
         SYSTEM   SIG5F
PROGSECT CSECT    1
*
*  REGISTERS
*
N        EQU      1                 XSEG EXECUTION-TIME INDEX
N1       EQU      4                 XSEG EXECUTION REG
X        EQU      1                 GENERAL INDEX REG
K        EQU      2                 SUBSCRIPT REG
T        EQU      2                 TYPE
XL       EQU      3                 XSEG LOC
A        EQU      4                 ARG ADR/INDEX
LX       EQU      5                 INDEX LINK REG
LX7      EQU      7                 INDEX LINK REG
GIJ      EQU      1                 INDEX TO ACCESS G(I,J)
GIK      EQU      2                 INDEX TO ACCESS G(I,K)
GKI      EQU      3                 INDEX TO ACCESS G(K,I)
VI       EQU      2                 INDEX TO ACCESS V(I)
RI       EQU      2                 INDEX TO ACCESS R(I)
FJ       EQU      3                 INDEX TO ACCESS F(J)
FI       EQU      3                 INDEX TO ACCESS F(I)
X0       EQU      2                 EVEN/ODD INDEX PAIR
X1       EQU      3                   *
AIJ      EQU      4                 INDEX TO ACCESS A(I,K)
AF       EQU      6                 ACCUM FOR FLOT VALUES
AF1      EQU      7                     *
BF       EQU      8                 2ND ACCUM FOR FLOT VALUES
BF1      EQU      9                     *
BUF      EQU      7                 BUFFER FOR MOVING DATA/CODE GROUPS
R        EQU      8                 GENERAL WORK REG
S        EQU      11                SIZE
NI       EQU      12                I-LOOP COUNT
NJ       EQU      13                J-LOOP COUNT
L3       EQU      12                LINK REG
L2       EQU      13                LINK REG
L1       EQU      14                LINK REG
         PAGE
*
*
*  P R O C S
*
*
*              (THIS MODULE RESIDES IN AN OVERLAY, SO THE USUAL
*              TEMP ALLOCATION PROCEDURE WILL NOT WORK).
*
TLOC     SET      0
MATTEMPS EQU      XSEGBASE+50
*
TEMP     CNAME    1
DTEMP    CNAME    2
         PROC
         DO1      NAME=2
TLOC     SET      TLOC+(TLOC&1)
         DISP     TLOC+50
LF       EQU      MATTEMPS+TLOC
TLOC     SET      TLOC+NAME
         PEND
*
*
EVEN     CNAME    0
ODD      CNAME    1
         PROC
LF       EQU      %
         ERROR,1,(CF(2)+NAME)&1   'REGISTER HAS WRONG PARITY'
         PEND
*
*
EQUAL    CNAME
         PROC
LF       EQU      %
         ERROR,1,1-(CF(2)=CF(3))  'REGISTERS MUST BE EQUAL'
         PEND
*
*
NB       CNAME    X'680'
NBGE     CNAME    X'681'
NBLE     CNAME    X'682'
NBE      CNAME    X'683'
NBL      CNAME    X'691'
NBG      CNAME    X'692'
NBNE     CNAME    X'693'
         PROC
         ERROR,1,(AF>=0)+(NUM(AF)>1)  'AF MUST BE NEG CONST ADR'
LF       GEN,12,20  NAME-1,AF
         PEND
         PAGE
*
*
*  XSEG GEN PROCS
*
*
         OPEN     GEN
GEN      CNAME
         OPEN     M,N,MN,I
         PROC
LF       EQU      %
         ERROR,1,1-(NUM(CF)=3)  'WRONG NUMBER OF CF ARGS'
M        SET      CF(2)
N        SET      CF(3)
MN       SET      M+N
         ERROR,1,1-(NUM(AF)=(M>0)+(N>0)) 'WRONG NUMBER OF AF ARGS'
         DO       M>0
I          DO       N
             LW,AF(1)+M+I-1  AF(2)+I-1
           FIN
I          DO       MN*(MN<3)
             STW,AF(1)+I-1   I-1,XL
           ELSE
             LCI         MN
             STM,AF(1)   0,XL
           FIN
         ELSE
I          DO       MN*(MN<3)
             LW,BUF      AF(1)+I-1
             STW,BUF     I-1,XL
           ELSE
             LCI         N
             LM,BUF      AF(1)
             STM,BUF     0,XL
           FIN
         FIN
         AI,XL    MN
         PEND
         CLOSE    M,N,MN,I
         PAGE
*
*
*  M A T R I X    I N V E R T / D I V I D E    O P    R O U T I N E S
*
*
         USECT    PROGSECT
MAT@     RES      0                 START OF PROCEDURE
*
*
MMATINV  EQU      %                 MONADIC MATRIX INVERT
         BAL,L1   MTYPEF            ALLOW NUMERIC ARG; SET FLOT RESULT
         LW,A     RTARG             GET ARG'S 1ST DIMEN (M)
         LW,R     2,A               CONTINUE LIKE MATRIX DIVIDE,
         B        MATOP1              PRETENDING THERE WAS AN M-BY-M
*                                     LEFT ARG.
*
*
DMATDIV  EQU      %                 DYADIC MATRIX INVERT
         BAL,L1   DTYPEF            ALLOW NUMERIC ARGS; SET FLOT RESULT
         LI,X     1
         LB,R    *LFARG,X           GET RANK OF LEFT ARG
         LW,A     LFARG
         CLM,R    1OR2              MUST BE 1 OR 2 (VECTOR OR MATRIX)
         BCS,9    ERRANK
         BCR,3    1Z1               IF RANK =1, SET P=1
         LW,R     3,A               RANK =2, P = 2'ND DIMEN
         BEZ      ERLENGTH
MATOP1   EQU      %
1Z1      STW,R    PVAL
         LW,R     2,A               M = 1'ST DIMEN (LFARG IS M-BY-P)
         STW,R    MVAL
         LW,A     RTARG
         LI,X     1
         LB,X    *RTARG,X           GET RIGHT ARG RANK
         CI,X     2
         BNE      ERRANK            MUST BE 2 (MATRIX)
         STW,X    RSRANK            SET RESULT TO MATRIX
         CW,R     2,A               RIGHT ARG MUST BE M-BY-N
         BNE      ERLENGTH
         LW,S     3,A                 WITH   0<N<=M.
         BEZ      ERLENGTH
         CW,S     MVAL
         BG       ERLENGTH
         STW,S    NVAL
         AW,S     PVAL              COMPUTE ROW SIZE OF G
         STW,S    ROWSIZE             (M BY N+P).
         ODD,S
         MW,S     MVAL              COMPUTE SIZE OF G
*                                   NOTE: THIS SIZE IS >0 SINCE
*                                     M*(N+P) >  M*N >= N*N > 0.
         STW,S    RSSIZE            WE'LL ALOCATE THE G DATA BLOCK
         BAL,L1   ALOCRS              IN RESULT POSITION. AT THE MOMENT,
*                                     IT'S BIGGER THAN THE FINAL RESULT;
*                                     WE'LL PARE IT DOWN LATER.
         LW,R     NVAL
         STW,R    2,A               SET RESULT DIMENS:
         LW,R     PVAL                N-BY-P.
         STW,R    3,A
         BDR,R    1Z5               IF P=1,
         LI,X     1                   REDUCE RESULT RANK TO 1
         STB,X   *RESULT,X            (2ND DIMEN BECOMES GAP WORD).
1Z5      LW,A     LFARG
         BEZ      1Z2               IF LFARG EXISTS (IE., MATRIX DIV),
         LI,A     0                   SET ITS ADR FOR INDEXED LOAD.
         BAL,LX   XSETUP
1Z2      LI,A     1                 SET RTADR FOR INDEXED LOAD
         BAL,LX   XSETUP
         LI,X     -1                SET RSADR FOR SEQUENCIAL STORE
         BAL,LX   SETADR
         BAL,LX   GXSEGINI          INIT XSEG CODE (RESULT ISN'T NULL)
         GEN,0,3  CODE1             GEN XSEG TO BUILD G;
         LI,A     1                   EITHER G = B,A  (CATENATED
         BAL,L1   GENLOAD                OR  G = B,I   COLUMN-WISE)
         GEN,0,6  CODE2               WHERE  A = LFARG,
         LW,R     RSADR                      B = RTARG,
         AWM,R    -6,XL                      I = IDENTITY MATRIX
         STW,XL   LOOPLOC
         LW,A     LFARG             TEST DIVIDE/INVERT
         BNEZ     1Z3
         GEN,0,8  CODE3             INVERT: GEN CODE TO BUILD IDENTITY
         AWM,XL   -5,XL               MATRIX, ROW BY ROW.
         AWM,XL   -7,XL
         LW,R     PVAL
         STW,R    TESTVAL
         B        1Z4
1Z3      LI,A     0                 DIVIDE: GEN CODE TO APPEND ROWS
         BAL,L1   GENLOAD             OF LFARG.
         GEN,0,3  CODE5
1Z4      LW,R     RSADR
         AWM,R    -3,XL
         B        EXECUTE           EXECUTE XSEG
*
*
*              XSEG CODE:
*
*        LCW,N    RSSIZE         0  INIT STORE INDEX / RESULT COUNT
CODE1    LI,K     0              1  INIT LEFT INDEX
         STW,K    KTEMP          2  INIT RIGHT INDEX
         LW,N1    NVAL           3  INIT LEFT LOOP COUNT
*        LOAD/CNV RTARG(K)       4  LEFT LOOP: MOVE A ROW OF 'B'
CODE2    STD,AF   0  RESULT(N)        TO RESULT.
         AI,N     1                 BUMP STORE AND
         AI,K     1                   LEFT INDECES.
         BDR,N1   XSEGBASE+4        COUNT LEFT LOOP
         XW,K     KTEMP             SAVE LEFT INDEX, GET RIGHT INDEX
         LW,N1    PVAL              INIT RIGHT LOOP COUNT
*
*           FOR INVERT:
*
*LOOPLOC EQU      %
CODE3    CW,N1    TESTVAL           RIGHT LOOP: COMPUTE A ROW OF 'I'
         NBNE     -4  (%+3)
         LD,AF    FLOT1             ONE ELEMENT IS 1.0,
         NB       -3  (%+2)
         LD,AF    FLOT0               ALL OTHERS ARE 0.0.
         STD,AF   0  RESULT(N)      MOVE TO RESULT (G)
         BIR,N    CODE4             BUMP RESULT INDEX
         B        XSEGDONE          EXIT WHEN G FILLED
CODE4    BDR,N1  *LOOPLOC        *  COUNT RIGHT LOOP
         MTW,-1   TESTVAL        *  END OF ROW, ADVANCE POSITION OF 1.0
         XW,K     KTEMP          *  SAVE RIGHT INDEX (UNUSED), GET LEFT
         B        XSEGBASE+3     *  GO DO LEFT LOOP
*
*           FOR DIVIDE:
*
*LOOPLOC LOAD/CNV LFARG(K)          RIGHT LOOP: MOVE A ROW OF 'A'
CODE5    STD,AF   0  RESULT(N)        TO RESULT.
         BIR,N    CODE6             BUMP STORE INDEX
         B        XSEGDONE          EXIT WHEN G FILLED
CODE6    AI,K     1              *  BUMP  RIGHT INDEX
         BDR,N1  *LOOPLOC        *  COUNT RIGHT LOOP
         XW,K     KTEMP          *  SAVE RIGHT INDEX, GET LEFT
         B        XSEGBASE+3     *  GO DO LEFT LOOP
         PAGE
*
*
XSEGDONE EQU      %                 WE NOW HAVE OUR INITIAL G MATRIX
         LI,A     0                 DON'T ALLOW BREAKS DURING ALOC'S
         STW,A    XSEGBRK             AND DREF'S. (FLAG WAS SET
*                                     BY 'EXECUTE').
         XW,A     LFARG
         BEZ      2Z1               IF LFARG EXITS,
         BAL,LX7  DREF                GET RID OF IT.
         LI,A     0
2Z1      XW,A     RTARG             GET RID OF RTARG
         BAL,LX7  DREF
         LW,S     MVAL              ALOC DB FOR ROW-MAX VECTOR
         SLS,S    1                   R(1),...,R(M).
         BAL,LX7  ALOCHNW
         STW,A    LFARG
         LI,R     X'0100'           SET IT TO LEGIT TYPE
         STH,R   *LFARG
         LW,S     NVAL              ALOC DB FOR SCALING FACTOR VECTOR,
         SLS,S    1                   F(1),...,F(N).
         BAL,LX7  ALOCHNW
         STW,A    RTARG
         LI,R     X'0100'           SET IT TO LEGIT TYPE
         STH,R   *RTARG
         LI,R     OPBREAK           NOW THAT ALL MEMORY FIDDLING IS
         STW,R    XSEGBRK             DONE, IT'S SAFE TO ALLOW BREAKS.
         AI,A     2                 GET PAST RTARG HEADER
         SLS,A    -1
         STW,A    F1ADR             SET DBLWORD ADR OF F(1)
         LW,A     LFARG
         AI,A     2                 GET PAST LFARG HEADER
         STW,A    V1ADR             SET WORD ADR OF V(1) (INTG.VECTOR)
         SLS,A    -1
         STW,A    R1ADR             SET DBLWORD ADR OF R(1)
         LW,A     RESULT
         AI,A     4                 GET PAST RESULT HEADER AND DIMENS
         SLS,A    -1
         STW,A    G11ADR            SET DBLWORD ADR OF G(1,1)
         PAGE
*
*
*  DO SCALING
*
*           ROW SCALING: FOR I=1 TO M,
*                        SET R(I) = MAX OF ABS(G(I,J)) FOR J=1 TO N
*
         LW,RI    R1ADR
         LW,GIJ   G11ADR
         LW,NI    MVAL              I LOOP COUNT = M
         LD,AF    FLOT0             INIT BIGGEST ELMT = 0.0
3Z1      STD,AF   BIGGEST           REMEMBER BIGGEST ELEMENT
         LD,AF    FLOT0             I LOOP: INIT R(I)=0.0
         STD,AF   0,RI
         LW,NJ    NVAL              J LOOP COUNT = N
3Z2      LAD,AF   0,GIJ             J LOOP
         CD,AF    0,RI              SET R(I)=MAX(R(I),ABS(G(I,J)))
         BLE      3Z3
         STD,AF   0,RI
3Z3      AI,GIJ   1                 MOVE GIJ ALONG ROW
         BDR,NJ   3Z2               COUNT J LOOP
         LD,AF    BIGGEST           IS THIS STILL THE BIGGEST ?
         CD,AF    0,RI              IF NOT,
         BGE      3Z31                UPDATE IT.
         LD,AF    0,RI
3Z31     AW,GIJ   PVAL              MOVE ACCROSS LAST P COLS, TO NXT ROW
         AI,RI    1                 MOVE RI TO NEXT ELMT
         BDR,NI   3Z1               COUNT I LOOP
         FML,AF   EPSILON           SET PIVOT TEST VALUE
         STD,AF   PIVOTMIN            = BIGGEST ELEMENT * EPSILON.
*
*           COLUMN SCALING: FOR J=1 TO N,
*                           SET F(J) = MAX OF ABS(G(I,J))/R(I)
*                                      FOR I = 1 TO M.
*
         LW,FJ    F1ADR
         LW,GIJ   G11ADR
         LW,NJ    NVAL              J LOOP COUNT = N
3Z4      STW,GIJ  G1JADR            J LOOP: REMEMBER ADR OF G(1,J)
         LD,AF    FLOT0
         STD,AF   0,FJ              INIT F(J)=0
         LW,RI    R1ADR             INIT R TO R(1)
         LW,NI    MVAL              I LOOP COUNT = M
3Z5      LAD,AF   0,GIJ             I LOOP
         FDL,AF   0,RI              A = ABS(G(I,J))/R(I)
         CD,AF    0,FJ              SET F(J)= MAX(F(J),A)
         BLE      3Z6
         STD,AF   0,FJ
3Z6      AW,GIJ   ROWSIZE           BUMP G DOWN COLUMN
         AI,RI    1                 BUMP R
         BDR,NI   3Z5               COUNT I LOOP
         AI,FJ    1                 BUMP TO NEXT F(J)
         LW,GIJ   G1JADR            RECALL ADR OF G(1,J)
         AI,GIJ   1                 BUMP TO NEXT COL
         BDR,NJ   3Z4               COUNT J LOOP
*
*           APPLY SCALE FACTORS TO MATRIX COLUMNS:
*                   FOR I=1 TO M, AND J=1 TO N,
*                   SET G(I,J) = G(I,J)/F(J)
*
         LW,GIJ   G11ADR
         LW,NI    MVAL              I LOOP COUNT = M
3Z7      LW,FJ    F1ADR
         LW,NJ    NVAL              J LOOP COUNT = N
3Z8      LD,AF    0,GIJ
         FDL,AF   0,FJ              SET G(I,J) = G(I,J)/F(J)
         STD,AF   0,GIJ
         AI,GIJ   1                 BUMP G ALONG ROW
         AI,FJ    1                 BUMP TO NEXT F
         BDR,NJ   3Z8               COUNT J LOOP
         AW,GIJ   PVAL              BUMP G PAST LAST P COLS, TO NXT ROW
         BDR,NI   3Z7               COUNT I LOOP
*
*           INITIALIZE COLUMN PERMUTATION VECTOR:
*              SET V(I) = I-1 FOR I = 1 TO N.
*
         LW,X     NVAL              I LOOP COUNT = N
3Z9      STW,X   *V1ADR,X           V(I+1)=I
         BDR,X    3Z9               COUNT I LOOP
         STW,X   *V1ADR             V(1)=0
         PAGE
*
*
*  PERFORM HOUSEHOLDER TRANSFORMATIONS,
*  LEADING TO UPPER TRIANGULAR SYSTEM.
*
*           INITIALIZE K LOOP PARAMETERS
*                K RUNS FROM 1 TO N.
*
         LW,GIJ   G11ADR            SET INIT (K=1) VALUES FOR:
         STW,GIJ  GK1ADR                ADR(G(K,1)),
         STW,GIJ  GKKADR                ADR(G(K,K)),
         LW,R     MVAL
         STW,R    MK1VAL                M-K+1,
         LW,R     NVAL
         STW,R    NK1VAL                N-K+1,
         LI,R     0
         STW,R    K1VAL                 K-1.
*
*           SELECT PIVOT ELEMENT: FIND S, T SUCH THAT
*                    ABS(G(S,T)) = MAX(ABS(G(I,J)))
*                                  FOR I = K TO M, J = K TO N.
*
KLOOP    LW,GIJ   GK1ADR            START G AT G(K,1)
         LD,AF    FLOT0             INIT ABS(G(S,T)) = 0
         STD,AF   ABSGST
         LW,NI    MK1VAL            I LOOP COUNT = M-K+1
4Z2      AW,GIJ   K1VAL             I LOOP: SET G TO G(I,K)
         LW,NJ    NK1VAL            J LOOP COUNT = N-K+1
4Z3      LAD,AF   0,GIJ             J LOOP:
         CD,AF    ABSGST            IF ABS(G(I,J)) > ABSGST,
         BLE      4Z4
         STD,AF   ABSGST              SET ABSGST = ABS(G(I,J)),
         STW,GIJ  GSTADR              AND REMEMBER S AND T.
4Z4      AI,GIJ   1                 BUMP G ALONG ROW
         BDR,NJ   4Z3               COUNT J LOOP
         AW,GIJ   PVAL              BUMP G PAST LAST P COLS TO NEXT ROW
         BDR,NI   4Z2               COUNT I LOOP
         LD,AF    ABSGST
         CD,AF    PIVOTMIN          IF PIVOT ELMT DANGEROUSLY LOW,
         BLE      ERSING              PROCLAIM THE MATRIX SINGULAR.
*
*           PERMUTE COLUMNS K/T: FOR I = 1 TO M,
*                    EXCHANGE G(I,K) AND G(I,T),
*                    EXCHANGE V(K) AND V(T).
*
         LW,X1    GSTADR            GET ADR(G(S,T))
         LI,X0    0
         EVEN,X0                    X0/X1 ARE AN EVEN/ODD PAIR
         EQUAL,X0+1,X1
         SW,X1    G11ADR            SEPARATE S,T INDEX INTO DISTINCT
         DW,X0    ROWSIZE             ROW & COL INDECES:
         STW,X0   T1VAL               X0 = T-1,
*                                     X1 = S-1.
         CW,X0    K1VAL             IF T=K, NO EXCHANGE NEEDED
         BE       4Z6
         LW,GIJ   K1VAL             GET K-1
         LW,R    *V1ADR,X0          EXCHANGE V(T) AND V(K)
         XW,R    *V1ADR,GIJ
         STW,R   *V1ADR,X0
         AW,X0    G11ADR            COMPUTE ADR(G(1,T)),
         AW,GIJ   G11ADR                AND ADR(G(1,K)).
         LW,NI    MVAL              I LOOP COUNT = M
4Z5      LD,AF    0,GIJ             I LOOP
         LD,BF    0,X0              EXCHANGE G(I,K) AND G(I,T)
         STD,AF   0,X0
         STD,BF   0,GIJ
         AW,GIJ   ROWSIZE           BUMP G DOWN COLUMN
         AW,X0    ROWSIZE               *
         BDR,NI   4Z5               COUNT I LOOP
4Z6      EQU      %
*
*           PERMUTE ROWS S/K: FOR J = 1 TO N+P,
*                    EXCHANGE G(S,J) AND G(K,J).
*
         CW,X1    K1VAL             IF S=K, NO EXCHANGE NECESSARY
         BE       4Z8
         LW,X1    GSTADR
         SW,X1    T1VAL             GET ADR(G(S,1)),
         LW,GIJ   GK1ADR            AND ADR(G(K,1)).
         LW,NJ    ROWSIZE           J LOOP COUNT = N+P
4Z7      LD,AF    0,GIJ             J LOOP
         LD,BF    0,X1              EXCHANGE G(S,J) AND G(K,J)
         STD,AF   0,X1
         STD,BF   0,GIJ
         AI,GIJ   1                 BUMP G ALONG ROW
         AI,X1    1                     *
         BDR,NJ   4Z7               COUNT J LOOP
4Z8      EQU      %
*
*
*           APPLY K'TH TRANSFORMATION:
*
*                    S = SQRT(G(K,K)**2 + G(K+1,K)**2 +...+ G(M,K)**2)
*                    U = G(K,K) + S*SGN(G(K,K))
*                    B = 1/(S*S + S*ABS(G(K,K)))
*                    FOR J = K+1 TO N+P,
*                        Y = B*U*G(K,J)+B*(G(K+1,J)*G(K+1,K) +...
*                                     ...+ G(M,J)*G(M,K))
*                        G(I,J):=G(I,J)-Y*G(I,K) FOR I=K+1 TO M
*                        G(K,J):=G(K,J)-U*Y
*                    G(K,K):=-S*SGN(G(K,K)) = G(K,K)-U
*
         LD,AF    FLOT0             INIT SUM (FOR S**2) TO 0.0
         LW,NI    MK1VAL            I LOOP COUNT = M-K+1
         LW,GIJ   GKKADR            START AT G(K,K)
5Z0      STD,AF   S2VAL             I LOOP
         LD,AF    0,GIJ             ADD IN G(I,K)**2
         FML,AF   0,GIJ
         FAL,AF   S2VAL
         AW,GIJ   ROWSIZE           BUMP G DOWN COLUMN
         BDR,NI   5Z0               COUNT I LOOP
         BAL,LX   FSQRT             COMPUTE S=SQRT(S**2)
         STD,AF   SVAL              SAVE S
         FAL,AF   ABSGKK            S+ABS(G(K,K))
         FML,AF   SVAL              S**2 + S*ABS(G(K,K)) = 1/B
         LD,BF    FLOT1
         FDL,BF   AF                B
         STD,BF   BVAL
         LW,GIJ   GKKADR            GET ADR(G(K,K))
         LD,AF    0,GIJ
         BGEZ     5Z1               IF G(K,K)<0,
         FSL,AF   SVAL                U = G(K,K)-S,
         LD,BF    SVAL                G(K,K):=S.
         B        5Z2               IF G(K,K)>0,
5Z1      FAL,AF   SVAL                U = G(K,K)+S,
         LCD,BF   SVAL                G(K,K):=-S.
5Z2      STD,AF   UVAL              U = G(K,K)+S*SGN(G(K,K))
         STD,BF   0,GIJ             G(K,K) := -S*SGN(G(K,K))
         LW,NJ    NK1VAL
         AW,NJ    PVAL              J LOOP COUNT = N+P-K
         B        5Z8
5Z3      STW,GIJ  GKJADR            J LOOP: SET ADR(G(K,J)),
         LW,GIK   GKKADR                    AND ADR(G(K,K)).
         LD,AF    0,GIJ             START COMPUTATION OF Y:
         FML,AF   UVAL                AF = U*G(K,J).
         LW,NI    MK1VAL            I LOOP COUNT = M-K (POSSIBLY 0)
         B        5Z5               START I LOOP WITH I=K+1
5Z4      STD,AF   YVAL              I LOOP:
         LD,AF    0,GIK               AF=AF+G(I,K)*G(I,J).
         FML,AF   0,GIJ
         FAL,AF   YVAL
5Z5      AW,GIK   ROWSIZE           BUMP G DOWN COLUMN
         AW,GIJ   ROWSIZE
         BDR,NI   5Z4               COUNT I LOOP
         FML,AF   BVAL              Y = B*AF
         STD,AF   YVAL
         LW,GIK   GKKADR            RESTORE INIT G ADR'S
         LW,GIJ   GKJADR
         LCD,AF   UVAL              G(K,J):=G(K,J)-Y*U
         LW,NI    MK1VAL            I LOOP COUNT = M-K (POSSIBLY 0)
         B        5Z7               START I LOOP WITH I=K+1
5Z6      LCD,AF   0,GIK
5Z7      FML,AF   YVAL              G(I,J):=G(I,J)-Y*G(I,K)
         FAL,AF   0,GIJ
         STD,AF   0,GIJ
         AW,GIK   ROWSIZE           BUMP G DOWN COLUMN
         AW,GIJ   ROWSIZE
         BDR,NI   5Z6               COUNT I LOOP
         LW,GIJ   GKJADR            UPDATE ADR(G(K,J)) BY
5Z8      AI,GIJ   1                   BUMPING IT ALONG ROW K.
         BDR,NJ   5Z3               COUNT J LOOP
*
*
*           UPDATE K LOOP PARAMETERS
*
         MTW,-1   NK1VAL            UPDATE N-K+1 VALUE
*                                   IF N-K+1=0 (K>N),
         BEZ      SOLVE               WE'RE DONE.
         MTW,-1   MK1VAL            K<=N: UPDATE M-K+1,
         MTW,1    K1VAL                          K-1,
         LW,S     ROWSIZE
         AWM,S    GK1ADR                         ADR(G(K,1)),
         AI,S     1
         AWM,S    GKKADR                     AND ADR(G(K,K)).
         B        KLOOP             CONTINUE K LOOP
         PAGE
*
*
*  SOLVE UPPER TRIANGULAR SYSTEM
*
*
*           FOR K = N TO 1, AND J = N+1 TO N+P, SET
*                    G(K,J) := (G(K,J) - G(K,K+1)*G(K+1,J) - ...
*                                  ... - G(K,N)*G(N,J))/G(K,K)
*
SOLVE    MTW,1    NK1VAL            RESTORE N-K+1 VALUE
         LW,GIJ   GKKADR            START G(K,N) AT G(N,N)
6Z1      STW,GIJ  GKNADR            REMEMBER ADR(G(K,N))
         AI,GIJ   1                 START G(K,J) AT G(K,N+1)
         LW,NJ    PVAL              J LOOP COUNT = P
6Z2      STW,GIJ  GKJADR            REMEMBER ADR(G(K,J))
         LW,GKI   GKKADR            START G(K,I) AT G(K,K)
         LD,AF    0,GIJ             START SUM WITH G(K,J)
         LW,NI    NK1VAL            I LOOP COUNT = N-K (POSSIBLY 0)
         B        6Z4               START WITH I=K+1
6Z3      STD,AF   SUM
         LCD,AF   0,GKI             SUM := SUM-G(K,I)*G(I,J)
         FML,AF   0,GIJ
         FAL,AF   SUM
6Z4      AI,GKI   1                 BUMP G(K,I) ALONG ROW
         AW,GIJ   ROWSIZE           BUMP G(I,J) DOWN COLUMN
         BDR,NI   6Z3               COUNT I LOOP
         LW,GIJ   GKKADR            COMPUTE NEW G(K,J)
         FDL,AF   0,GIJ               = SUM/G(K,K).
         LW,GIJ   GKJADR
         STD,AF   0,GIJ             STORE NEW G(K,J)
         AI,GIJ   1                 UPDATE ADR(G(K,J))
         BDR,NJ   6Z2               COUNT J LOOP
         MTW,-1   K1VAL             COUNT K LOOP (UPDATE K-1)
         BLZ      NORMALYZ          IF K=0, WE'RE DONE
         MTW,1    NK1VAL            K>0: UPDATE N-K+1,
         LCW,R    ROWSIZE
         AI,R     -1
         AWM,R    GKKADR                        ADR(G(K,K)),
         LW,GIJ   GKNADR                    AND ADR(G(K,N)).
         SW,GIJ   ROWSIZE
         B        6Z1
         PAGE
*
*
*  NORMALIZE THE SOLUTION TO CORRECT FOR
*  COLUMN SCALING AND INTERCHANGES.
*
*           FOR I = 1 TO N, AND J = 1 TO P,
*              EXCHANGE G(I,N+J) WITH G(L,N+J),
*              AND COPY V(I) TO V(L), WHERE
*              V(L)=I.
*
NORMALYZ LW,AIJ   G11ADR            START 'A' AT A(1,1)
         LW,NI    NVAL              I LOOP COUNT = N
         LI,VI    0                 INIT TO V(1)
7Z1      CW,VI   *V1ADR,VI          IF V(I)=I, NO EXCHANGE NEEDED
         BNE      7Z10
         AW,AIJ   ROWSIZE             (IN WHICH CASE, BUMP TO NEXT ROW)
         B        7Z3
7Z10     LW,GIJ   VI                L:=I
7Z11     LW,GIJ  *V1ADR,GIJ         L:=V(L)
         CW,VI   *V1ADR,GIJ         SEARCH V FOR I: IF I.NE.V(L),
         BNE      7Z11                SET L:=V(L) AND TRY AGAIN.
         LW,R    *V1ADR,VI          COPY V(I)
         STW,R   *V1ADR,GIJ           TO V(L).
         ODD,GIJ
         MW,GIJ   ROWSIZE           COMPUTE ADR(G(L,N+1))
         AW,GIJ   NVAL
         AW,GIJ   G11ADR
         AW,AIJ   NVAL              MOVE 'A' TO A(I,N+1)
         LW,NJ    PVAL              J LOOP COUNT = P
7Z2      LD,AF    0,GIJ             EXCHANGE G(L,N+J)
         LD,BF    0,AIJ               AND A(I,N+J).
         STD,AF   0,AIJ
         STD,BF   0,GIJ
         AI,GIJ   1                 BUMP G ALONG ROW
         AI,AIJ   1                 BUMP A ALONG ROW
         BDR,NJ   7Z2               COUNT J LOOP
7Z3      AI,VI    1                 BUMP TO NEXT V(I)
         BDR,NI   7Z1               COUNT I LOOP
*
*
*           FOR I = 1 TO N, AND J = 1 TO P,
*              SET A(I,J) = G(I,N+J)/F(I).
*           'A' IS THE N-BY-P CORNER OF G IN
*           WHICH THE FINAL RESULT APPEARS.
*
         LW,AIJ   G11ADR            START A AT A(1,1)
         LW,GIJ   G11ADR            START G AT G(1,1)
         LW,FI    F1ADR             START F AT F(1)
         LW,NI    NVAL              I LOOP COUNT = N
7Z4      AW,GIJ   NVAL              START G(I,J) AT G(I,N+1)
         LW,NJ    PVAL              J LOOP COUNT = P
7Z5      LD,AF    0,GIJ             GET G(I,N+J)
         FDL,AF   0,FI              DIVIDE BY F(I)
         STD,AF   0,AIJ             STORE AS A(I,J)
         AI,GIJ   1                 BUMP G AND A
         AI,AIJ   1                   ALONG ROW.
         BDR,NJ   7Z5               COUNT J LOOP
         AI,FI    1                 BUMP TO NEXT F(I)
         BDR,NI   7Z4               COUNT I LOOP
*
*
*           NOW GIVE BACK THE M*(N+P)-N*P DOUBLEWORDS OF 'G' THAT
*           ARE NO LONGER NEEDED.
*
         LI,A     0                 MUSTN'T ALLOW BREAKS DURING MEMORY
         STW,A    XSEGBRK             MANAGEMENT OPERATIONS TO FOLLOW.
         LW,S     NVAL              FINAL RESULT SIZE = N*P
         MW,S     PVAL
         XW,S     RSSIZE            COMPUTE NR OF DBLWORDS TO GIVE BACK
         SW,S     RSSIZE              = OLD G SIZE - RESULT SIZE.
         SLS,S    1                 CHANGE TO NR OF WORDS
         LW,A     RESULT
         BAL,LX7  GIVEBACK          GIVE THEM BACK
         LI,A     DXRETURN          IF DYADIC ENTRY USED, SIMPLY EXIT
         CW,A     RETURN              DYADICALLY; BOTH LFARG (V(I))
         BE       DXRETURN            AND RTARG (F(J)) WILL BE DEREFFED.
         LI,A     0                 OTHERWISE, WE HAVE TO DEREF
         XW,A     LFARG               LFARG OURSELVES.
         BAL,LX7  DREF
         B        MXRETURN
         PAGE
*
*
*  DATA / TEMPS
*
*
         BOUND    8
1OR2     DATA     1,2               RANGE 1 TO 2, FOR ARG RANK TEST
EPSILON  DATA     X'34100000',0     16**(-13) FOR SINGULARITY TEST
*
*
V1ADR    TEMP                       WORD ADR OF V(1)
R1ADR    TEMP                       DBLWORD ADR OF R(1)
F1ADR    TEMP                       DBLWORD ADR OF F(1)
G11ADR   TEMP                       DBLWORD ADR OF G(1,1)
GK1ADR   TEMP                       DBLWORD ADR OF G(K,1)
GKKADR   TEMP                       DBLWORD ADR OF G(K,K)
GKJADR   TEMP                       DBLWORD ADR OF G(K,J)
GKNADR   TEMP                       DBLWORD ADR OF G(K,N)
G1JADR   TEMP                       DBLWORD ADR OF G(1,J)
GSTADR   TEMP                       DBLWORD ADR OF G(S,T)
*
MVAL     TEMP                       VALUE OF M
NVAL     TEMP                       VALUE OF N
PVAL     TEMP                       VALUE OF P
K1VAL    TEMP                       VALUE OF K-1
T1VAL    TEMP                       VALUE OF T-1
MK1VAL   TEMP                       VALUE OF M-K+1
NK1VAL   TEMP                       VALUE OF N-K+1
ROWSIZE  TEMP                       VALUE OF N+P
TESTVAL  TEMP                       INDEX OF NEXT DIAGONAL POS. (XSEG)
KTEMP    TEMP                       LEFT/RIGHT ARG INDEX        (XSEG)
*
ABSGST   DTEMP                      VALUE OF ABS(G(S,T))
ABSGKK   EQU      ABSGST                   = ABS(G(K,K)) AFTER EXCHANGES
S2VAL    DTEMP                      SUM(I=K,M) OF G(I,K)**2
SVAL     EQU      S2VAL             SQRT OF ABOVE
BVAL     DTEMP                      1/(S**2 + S*ABS(G(K,K)))
UVAL     DTEMP                      G(K,K)  + S*SGN(G(K,K))
YVAL     DTEMP                      B*U*G(K,J)+B*SUM(G(I,K)*G(I,J))
SUM      DTEMP                      G(K,J)-SUM(G(K,I)*G(I,J))
PIVOTMIN DTEMP                      PIVOT TEST VALUE
BIGGEST  DTEMP                      BIGGEST ELEMENT OF MATRIX
*
*
7Z       END

