(FILECREATED " 7-Oct-85 22:32:11" {QV}<IDL>SOURCES>MATRIX.;9 20678
changes to: (FNS NORM)
previous date: " 5-Sep-85 18:45:25" {QV}<IDL>SOURCES>MATRIX.;8)
(* Copyright (c) 1984, 1985 by Xerox Corporation. All rights reserved.)
(PRETTYCOMPRINT MATRIXCOMS)
(RPAQQ MATRIXCOMS ((* Matrix oriented functions: SWEEP, matrix product etc.)
(FNS INVERT MPROD NORM SWEEP SWEEP.SLTR SWEEPER SWEEP.COMPOSE SWEEP.UNDOPERM)))
(* Matrix oriented functions: SWEEP, matrix product etc.)
(DEFINEQ
(INVERT
[ULAMBDA ((M (EXPECTS MATRIX))
(RETURNS MATRIX))
(* bas: " 8-SEP-80 15:42" posted: "25-NOV-77 21:53")
(* Inverts a matrix via SWEEP adjusting the sign of the
entries on the way out.)
(DPROG ((S (if (AND (NEQ M (UARG M))
(type? SIMARRAY M)
(EQ (fetch AELTTYPE of M)
(QUOTE FLOATING)))
then M
else (COPYIDLARRAY M (QUOTE FLOATING))) MATRIX)
(RETURNS MATRIX))
(if (HASNILS (fetch ELEMENTBLOCK of S))
then (UERROR "Input contains missing data"))
[SETQ S (SWEEPER S 1.0 (SWEEP.SLTR M (QUOTE ALL]
(bind P [GSB ←(SETUP S (if (EQ (fetch FORMAT of S)
(QUOTE SYMMETRIC))
then (QUOTE SYMMETRIC)
else (QUOTE DONTCARE]
until (fetch DONE of GSB)
do (SETQ P (NEXT GSB))
(SETAELT S P (FMINUS (GETAELT S P))) (* Removing minus sign))
(SETTITLE S (MAKETITLE NIL M))
(RETURN S))])
(MPROD
[ULAMBDA ((A (EXPECTS MATRIX (ONEOF MATRIX VECTOR)))
(B (EXPECTS MATRIX (ONEOF MATRIX VECTOR)))
(RETURNS MATRIX))
(* jop: "12-Nov-84 16:20" posted: "14-MAY-78 20:49")
(* Matrix product. If the args are not matrices MPROD
tries to make suitable matrices of them)
(DPROG ((SA (fetch SHAPE of A) ROWINT)
(SB (fetch SHAPE of B) ROWINT)
THEN (C (if (EQ (fetch NELTS of SA)
2)
then (GETRELT SA 2)
elseif (EQ (fetch NELTS of SB)
2)
then (GETRELT SB 1)
else 0) INTEGER)
THEN (AX [if (EQ (fetch NELTS of SA)
2)
then (CONV.SIMARRAY A)
else (create SIMARRAY
using (CONV.SIMARRAY A)
SHAPE ←(if (OR (EQ C 0)
(NEQ (GETRELT SA 1)
C))
then (ROWINTOF (GETRELT SA 1)
1)
else (ROWINTOF 1 (GETRELT SA 1] SIMARRAY)
(BX [if (EQ (fetch NELTS of SB)
2)
then (CONV.SIMARRAY B)
else (create SIMARRAY
using (CONV.SIMARRAY B)
SHAPE ←(if (OR (EQ C 0)
(NEQ (GETRELT SB 1)
C))
then (ROWINTOF 1 (GETRELT SB 1))
else (ROWINTOF (GETRELT SB 1)
1] SIMARRAY)
THEN (SAX (fetch SHAPE of AX) ROWINT)
(SBX (fetch SHAPE of BX) ROWINT)
(MPTYPE (if (OR (EQ (fetch AELTTYPE of AX)
(QUOTE FLOATING))
(EQ (fetch AELTTYPE of BX)
(QUOTE FLOATING)))
then (QUOTE FLOATING)
else (QUOTE INTEGER)) (MEMQ FLOATING INTEGER))
THEN (MP (ALLOC.SARRAY (ROWINTOF (GETRELT SAX 1)
(GETRELT SBX 2))
(QUOTE FULL)
MPTYPE) MATRIX))
(if (NEQ (GETRELT SAX 2)
(GETRELT SBX 1))
then (UERROR "Arguments not conformable"))
(DPROG ((AXROW (FSELECT AX (ROWPTROF 1 (QUOTE ALL))) SELARRAY
(* Used to access the rows)))
[bind (AXROWGSB ←(SETUP AXROW (QUOTE ROWMAJOR)))
(BXGSB ←(SETUP BX (QUOTE COLMAJOR)))
(MPGSB ←(SETUP MP (QUOTE ROWMAJOR)))
(NCOLS ←(GETRELT SBX 2))
declare (AXROWGSB GENSTATEBLOCK (* Setup here to save the GSB))
(BXGSB GENSTATEBLOCK)
(MPGSB GENSTATEBLOCK)
(NCOLS INTEGER (* Number of output columns))
for ROW to (GETRELT SAX 1) declare (ROW IJK)
do (RESETUP BXGSB)
(ADJUST.SELECTION AXROW 1 ROW)
(for COL AELT BELT to NCOLS
declare (COL IJK)
(AELT SCALAR)
(BELT SCALAR)
do (RESETUP AXROWGSB) (* Branch on elementtype to suppress boxes)
(SETAELT MP (NEXT MPGSB)
(SELECTQ MPTYPE
(INTEGER (bind (TOT ← 0) declare (TOT FIXP)
until (fetch DONE of AXROWGSB)
do (if [AND (SETQ AELT (GETAELT AXROW
(NEXT
AXROWGSB)))
(SETQ BELT (GETAELT BX
(NEXT BXGSB]
then (add TOT (ITIMES AELT BELT))
else (SKIP BXGSB (fetch ELTCNT
of AXROWGSB))
(* SKIP and RETURN to avoid unnecessary work.)
(RETURN NIL))
finally (RETURN TOT)))
(FLOATING (bind (TOT ← 0.0) declare (TOT FLOATP)
until (fetch DONE of AXROWGSB)
do (if [AND (SETQ AELT (GETAELT
AXROW
(NEXT AXROWGSB)))
(SETQ BELT (GETAELT
BX
(NEXT BXGSB]
then (fadd TOT (FTIMES AELT BELT))
else (SKIP BXGSB (fetch ELTCNT
of AXROWGSB))
(RETURN NIL))
finally (RETURN TOT)))
(SHOULDNT])
(if LABPROPFLAG
then (AND (EQ (GETRELT SA 1)
(GETRELT SAX 1))
(LAB.COPYDIM A MP 1 1)) (* if SA$1#SAX$1 A's first dimension has been
compressed out)
(if (EQ (fetch NELTS of SB)
1)
then (AND (EQ (GETRELT SB 1)
(GETRELT SBX 2))
(LAB.COPYDIM B MP 1 2))
else (LAB.COPYDIM B MP 2 2)) (* Likewise B's)
(SETTITLE MP (MAKETITLE (QUOTE MATRIX% PRODUCT)
A B)))
(RETURN MP))])
(NORM
[ULAMBDA ((M (EXPECTS MATRIX))
(RETURNS MATRIX))
(* jop: " 5-Sep-85 17:23" posted: " 7-SEP-77 15:33")
(* Norms a matrix by dividing each entry by the sqrt of
its basis {diagonal} elements)
(DPROG ((GENORDER (if (EQ (fetch FORMAT of M)
(QUOTE SYMMETRIC))
then (QUOTE SYMMETRIC)
else (QUOTE ROWMAJOR)) (MEMQ SYMMETRIC ROWMAJOR))
(ASHAPE (fetch SHAPE of M) ROWINT)
THEN (MIN (IMIN (GETRELT ASHAPE 1)
(GETRELT ASHAPE 2)) INTEGER)
THEN (DIAG (create ROWFLOAT
NELTS ← MIN) ROWFLOAT)
(SEL NIL (ONEOF NIL ROWINT) (* If needed, selector for valid rows and columns))
(L 0 INTEGER (* number of diag elements of M gt 0 and non-nil)))
[for I ELT to MIN
do (SETQ ELT (GETAELT M (AELTPTR2 M I I)))
(if (AND ELT (GREATERP ELT 0))
then (add L 1)
(SETRELT DIAG L (SQRT ELT))
(if SEL
then (DPROGN ((SEL ROWINT))
(SETRELT SEL L I)))
elseif (NULL SEL)
then (SETQ SEL (create ROWINT
NELTS ←(DIFFERENCE MIN 1)))
(* The -1 is cause at least the current element is bad)
(for J from 1 to L declare (SEL ROWINT) do (SETRELT SEL J J]
(RETURN
(DPROG ((NORM (ALLOC.SARRAY (ROWINTOF L L)
(fetch FORMAT of M)
(QUOTE FLOATING)) MATRIX)
(VM (if SEL
then (replace NELTS of SEL with L)
(FSELECT M (create ROWPTR
NELTS ← 2
INIT ←(VFROMR SEL)))
else M) ARRAY))
[for I ELT DIAGI (GSBNORM ←(SETUP NORM GENORDER))
(GSBVM ←(SETUP VM GENORDER)) to L declare (DIAGI FLOATING)
do (SETQ DIAGI (GETRELT DIAG I))
(for J to (if (EQ GENORDER (QUOTE SYMMETRIC))
then I
else L)
do (SETAELT NORM (NEXT GSBNORM)
(AND (SETQ ELT (GETAELT VM (NEXT GSBVM)))
(FQUOTIENT ELT (FTIMES DIAGI (GETRELT DIAG J]
[if LABPROPFLAG
then
(LAB.COPYDIMS VM NORM)
(SETTITLE
NORM
(DPROG ((TM (GETTITLE M) TITLE)
(RETURNS TITLE)) (* Picking off the special cases of correlations)
(RETURN (OR (SELECTQ (CAR TM)
(COVARIATIONS (CONS (QUOTE CORRELATIONS)
(CDR TM)))
[SWEEP (if (EQ (CAR (LISTP (CADR TM)))
(QUOTE CROSSPRODUCTS))
then (LIST (QUOTE PARTIAL% CORRELATIONS)
(CADR (CADR TM]
NIL)
(MAKETITLE NIL M))))]
(RETURN NORM))))])
(SWEEP
[ULAMBDA ((M (EXPECTS MATRIX))
(OUTVARS ANY (* M selector for dimension 2 of M))
(INVARS ANY (* Another selector))
(RETURNS MATRIX))
(* bas: " 8-SEP-80 15:56" posted: "23-JUN-77 12:09")
(DPROG ((S (if (AND (NEQ M (UARG M))
(type? SIMARRAY M)
(EQ (fetch AELTTYPE of M)
(QUOTE FLOATING)))
then M
else (COPYIDLARRAY M (QUOTE FLOATING))) MATRIX))
(if (HASNILS (fetch ELEMENTBLOCK of S))
then (UERROR "Input contains missing values"))
[if OUTVARS
then (SETQ S (SWEEPER S 1.0 (SETQ OUTVARS (SWEEP.SLTR M OUTVARS]
[if INVARS
then (SETQ S (SWEEPER S -1.0 (SETQ INVARS (SWEEP.SLTR M INVARS]
[if LABPROPFLAG
then (DPROG ((OV [AND OUTVARS (for I to (fetch NELTS of OUTVARS) declare (OUTVARS ROWINT)
collect (OR (GETLEVLAB M 2 (GETRELT OUTVARS I))
(GETRELT OUTVARS I] LST
(* Vars to go out))
(IV [AND INVARS (for I to (fetch NELTS of INVARS) declare (INVARS ROWINT)
collect (OR (GETLEVLAB M 2 (GETRELT INVARS I))
(GETRELT INVARS I] LST
(* Vars to come in))
(TM (GETTITLE M) TITLE)
THEN (SWT (EQ (CAR (LISTP TM))
(QUOTE SWEEP)) BOOL (* M sweep title))
(COVART (EQ (CAR (LISTP TM))
(QUOTE COVARIATIONS)) BOOL
(* M covar title))
THEN (OL (if SWT
then (CDADR (CADDR TM))
elseif COVART
then (LIST (QUOTE Constant))) LST
(* Currently swept out vars))
(IL (AND SWT (CDADR (CADDDR TM))) LST
(* Curently swept in vars)))
(SETTITLE S (MAKETITLE NIL (if SWT
then (CADR TM)
elseif COVART
then (CONS (QUOTE CROSSPRODUCTS)
(CDR TM))
else TM)
(SWEEP.COMPOSE ", Swept out: " OL IL OV IV)
(SWEEP.COMPOSE ", Swept in: " IL OL IV OV))))]
(RETURN S))])
(SWEEP.SLTR
[DLAMBDA ((M MATRIX)
(SLTR ANY)
(RETURNS ROWINT))
(* bas: " 8-SEP-80 15:39")
(* Processes a sweep selector so that it is confined to
the upper left square section of a non square matrix)
(DPROG ((SV (MAKESLTR M SLTR 2) TTELT)
(MAXV (IMIN (GETRELT (fetch SHAPE of M)
1)
(GETRELT (fetch SHAPE of M)
2)) INTEGER (* Max sweepable variable)))
[RETURN (if (EQ SLTR (QUOTE ALL))
then (GENROW 1 MAXV)
else (DPROG ((V (CONV.ROWINT SV) ROWINT))
(for I to (fetch NELTS of V) when (LESSP MAXV (GETRELT V I))
do (UERROR "Sweep column out of range: " (GETRELT V I)))
(RETURN V))])])
(SWEEPER
[DLAMBDA ((S MATRIX (SATISFIES S:AELTTYPE='FLOATING))
(IO FLOATING (SATISFIES (OR (EQP IO 1.0)
(EQP IO -1.0))) (* 1 sweeps out -1 sweeps in))
(PC ROWINT (* Columns to sweep))
(RETURNS ARRAY))
(* jop: " 5-Sep-85 17:07" posted: "28-JUN-78 14:52")
(* Does the work for SWEEP by sideeffecting S which has been copied from the user argument in the caller.
It behaves slightly differently for symmetric and non-symmetric arrays. For the latter it uses total pivoting to
maximize numerical stability. For the former it just uses sequential diagonal pivots and trusts that the object is
positive definite. The algorithm is based on Dempster 1969, pp 62ff)
(DPROG ((SYM (EQ (fetch FORMAT of S)
(QUOTE SYMMETRIC)) BOOL)
(PFLG NIL BOOL (* Pivot permutation flag))
(PIV NIL FLOATING (* Pivot value))
(SHP (fetch SHAPE of S) ROWINT)
THEN (R (GETRELT SHP 1) INTEGER)
(C (GETRELT SHP 2) INTEGER)
THEN (PR (if SYM
then PC
else (COPYROW PC)) ROWINT (* Rows to sweep))
(BL (fetch NELTS of PC) INTEGER (SATISFIES (IGREATERP BL 0))
(* Number of swept variables)))
(if (AND (NOT SYM)
(IGREATERP (fetch REFCOUNT of PC)
0))
then (SETQ PC (COPYROW PC))) (* If FULL then PC will be changed)
[for K X Y PTR to BL
do [if (OR SYM (EQ K BL))
then (SETQ X (GETRELT PR K))
(SETQ Y (GETRELT PC K))
else (for I PRI (P ← K)
(Q ← K)
(MAXPIV ← 0.0) from K to BL declare (MAXPIV FLOATING)
do (SETQ PRI (GETRELT PR I))
(for J from K to BL
do [SETQ PIV (GETAELT S (AELTPTR2 S PRI (GETRELT PC J]
(OR (FGREATERP PIV 0.0)
(SETQ PIV (FMINUS PIV)))
(* ABS -
in line to avoid the box)
(if (FGREATERP PIV MAXPIV)
then (SETQ P I)
(SETQ Q J)
(SETQ MAXPIV PIV)))
finally (SETQ X (GETRELT PR P))
(if (NEQ K P)
then (SETRELT PR P (GETRELT PR K))
(SETRELT PR K X)
(SETQ PFLG T))
(SETQ Y (GETRELT PC Q))
(if (NEQ K Q)
then (SETRELT PC Q (GETRELT PC K))
(SETRELT PC K Y)
(SETQ PFLG T]
(SETQ PTR (AELTPTR2 S X Y))
[SETQ PIV (FMINUS (SETAELT S PTR (FQUOTIENT -1.0 (DPROG ((V (GETAELT S PTR) FLOATING)
(RETURNS FLOATING))
(if (EQP V 0.0)
then (UERROR
"Linear system has no solution")
else (RETURN V)))]
(* Setting the pivot)
[for I to R unless (EQ I X)
do (* pivot adjustments)
(for J to (if SYM
then I
else C)
unless (EQ J Y)
do (SETQ PTR (AELTPTR2 S I J))
(SETAELT S PTR (FDIFFERENCE (GETAELT S PTR)
(FTIMES (GETAELT S (AELTPTR2 S I Y))
PIV
(GETAELT S (AELTPTR2 S X J]
(* end of pivot adjustments)
(SETQ PIV (FTIMES PIV IO))
[for I to R unless (EQ I X)
do (SETQ PTR (AELTPTR2 S I Y))
(SETAELT S PTR (FTIMES PIV (GETAELT S PTR]
(if (NULL SYM)
then (for I to C unless (EQ I Y)
do (SETQ PTR (AELTPTR2 S X I))
(SETAELT S PTR (FTIMES PIV (GETAELT S PTR]
(RETURN (if PFLG
then (* Pivot choice has scrambled array -
so we unscramble it)
(SWEEP.UNDOPERM S PR PC)
else S)))])
(SWEEP.COMPOSE
[DLAMBDA ((HEAD STRINGP (* Labelling))
(CL LST (* Current var list))
(AL LST (* The other var list))
(IL LST (* Vars to be added))
(OL LST (* Vars to be removed))
(RETURNS LST))
(* bas: "28-JUN-78 13:55" posted: "28-JUN-78 14:16")
(* Composes variable lists for titles of SWEEP results)
[[LAMBDA (V)
(AND V (LIST HEAD (CONS (QUOTE LIST)
V]
(MERGE (SORT (for I in IL collect I unless (FMEMB I AL)))
(for I in CL collect I unless (FMEMB I OL]])
(SWEEP.UNDOPERM
[DLAMBDA ((S ARRAY)
(PR ROWINT (SATISFIES PR:NELTS=PC:NELTS) (* Row permutation))
(PC ROWINT [SATISFIES (for I to PC:NELTS always (for J to PR:NELTS
thereis (EQP PR$J PC$I]
(* Column permutation))
(RETURNS SELARRAY))
(* bas: "30-AUG-78 20:45" posted: "25-NOV-77 19:38")
(* * Selects so as to undo the permutation induced by full pivoting algorithm in SWEEP)
(DPROG ((N (fetch NELTS of PR) INTEGER)
(SHP (fetch SHAPE of S) ROWINT)
THEN (NR (GENROW 1 (GETRELT SHP 1)) ROWINT)
(NC (GENROW 1 (GETRELT SHP 2)) ROWINT))
(for I to N
do (SETRELT NR (GETRELT PC I)
(GETRELT PR I))
(SETRELT NC (GETRELT PR I)
(GETRELT PC I))) (* Inverting each permutation by the other)
[RETURN (FSELECT S (ROWPTROF (VFROMR NR)
(VFROMR NC])])
)
(PUTPROPS MATRIX COPYRIGHT ("Xerox Corporation" 1984 1985))
(DECLARE: DONTCOPY
(FILEMAP (NIL (513 20596 (INVERT 523 . 1774) (MPROD 1776 . 6745) (NORM 6747 . 9989) (SWEEP 9991 .
12807) (SWEEP.SLTR 12809 . 13802) (SWEEPER 13804 . 18464) (SWEEP.COMPOSE 18466 . 19423) (
SWEEP.UNDOPERM 19425 . 20594)))))
STOP