| The following routine headers describe these routines:
LIB$ABS
LIB$ACOPY
LIB$MAT_ADD
LIB$MAT_SUB
LIB$MAT_ADR
LIB$MAT_IDN
LIB$MAT_INI
LIB$MAT_INV
LIB$MAT_MUL
LIB$MAT_SCA
LIB$MAT_TRN
LIB$SQRT
.sbttl LIB$MAT_ADD/SUB/MUL Matrix Add/Subtract/Multiply
;---
;
; Calling sequence:
;
; CALL LIB$MAT_ADD( src1.rx.da, src2.rx.da, dest.wx.da )
; CALL LIB$MAT_SUB( src1.rx.da, src2.rx.da, dest.wx.da )
; CALL LIB$MAT_MUL( src1.rx.da, src2.rx.da, dest.wx.da )
;
; Inputs:
;
; src1 Address of descriptor for first matrix
; src2 Address of descriptor for second matrix
; dest Address of descriptor for sum matrix
;
; Description:
;
; These routines add/subtract/multiply two matrices, giving a third.
;
; dest = src1 + src2 (lib$mat_add)
; dest = src1 - src2 (lib$mat_sub)
; dest = src1 * src2 (lib$mat_mul)
;
; Any two of the three descriptors may be identically equal.
; The result of other overlap is undefined.
;
; The following data types are supported:
; W, L, F, D, G, H
;
; Returned status:
;
; ss$_normal Normal succesful completion
; lib$_invarg Unsupported datatype
; lib$_invarg Unsupported descriptor class
; lib$_invarg The array must have two dimensions
; lib$_invarg Error converting descriptor from class A to NCA
; lib$_invarg Differing datatypes
;
; Signalled status:
;
; ss$_fltovf_f Floating overflow fault
; ss$_opcdec Reserved opcode (only for data types G and H)
;
;---
.sbttl LIB$MAT_SCA Matrix Scalar Multiplication
;---
;
; Calling sequence:
;
; CALL LIB$MAT_SCA( src1.rx.rv, src2.rx.da, dest.wx.da )
;
; Inputs:
;
; src1 Address of scalar source item
; src2 Address of descriptor for source matrix
; dest Address of descriptor for destination matrix
;
; Description:
;
; This routine multiplies a matrix by a scalar, giving another matrix.
;
; The two matrices must have equivalent datatypes.
; The matrices may identically overlap.
;
; Returned status:
;
; ss$_normal Normal succesful completion
; Any error that can be returned by lib$$cvt_dsc.
; lib$_invarg Unsupported datatype
; lib$_invarg Unsupported descriptor class
; lib$_invarg The array must have two dimensions
;
;---
.sbttl LIB$MAT_TRN Matrix Transposition
;---
;
; Calling sequence:
;
; CALL LIB$MAT_TRN( src.rx.da, dest.wx.da )
;
; Inputs:
;
; src Address of descriptor for source matrix
; dest Address of descriptor for destination matrix
;
; Description:
;
; This routine transposes one matrix, giving a second.
;
; The matrices must have equivalent datatypes.
;
; Implementation:
;
; The array descriptor is converted to NCA descriptors, if needed.
;
; Returned status:
;
; ss$_normal Normal succesful completion
; Any error that can be returned by lib$acopy or lib$$cvt_dsc.
; lib$_invarg Unsupported datatype
; lib$_invarg Unsupported descriptor class
; lib$_invarg The array must have two dimensions
;
;---
.sbttl LIB$ACOPY Copy arrays
;---
;
; LIB$ACOPY Copy Multi-Dimensional Arrays
;
; Calling sequence:
;
; CALL LIB$ACOPY( srcdesc.dx, dstdesc.dx )
;
; Inputs:
;
; srcdesc(ap) Address of source array descriptor (A or NCA)
; dstdesc(ap) Address of destination array descriptor (A or NCA)
;
; Description:
;
; This routine copies the source array to the destination array.
; The arrays are passed by descriptor.
; If dsc$v_fl_bounds is 0, indices range from 0 to the multiplier-1.
; If dsc$v_fl_bounds is 1, indices range from the lower to upper bounds.
; Only items with indices common to both arrays are copied.
;
; Arrays with different dimensions are in error.
; The result of overlapping source and destination is undefined.
;
; The source and destination datatypes, lengths and scales must match.
;
; Implementation:
;
; The array descriptors are converted to NCA descriptors, if needed.
; Stack space is allocated for strides, maximum indices, etc.
; Lower bounds are maximized, and upper bounds are minimized.
; Maximum indices are computed as the difference in the bounds.
; Revert amounts are computed as stride x maximum index.
; All indices are varied in a loop (the first index most rapidly), and
; the array elements moved (by a MOVC3 instruction).
;
; Using the strides and the revert amounts, we are able to easily advance
; from one element to the next (with one subscript changed by 1), or to
; back up to the beginning of a row.
;
; addr( SRC(x1,...,xi+1,...,xdim) ) =
; addr( SRC(x1,...,xi,...,xdim) ) + src_stride[i]
; addr( SRC(x1,...,max_index[i],...,xdim) ) =
; addr( SRC(x1,...,0,...,xdim) ) + src_stride[i] * max_index[i] =
; addr( SRC(x1,...,0,...,xdim) ) + src_revert[i]
;
; And similarly for the destination array.
;
; Register usage:
;
; R10 = Address of source array descriptor
; R11 = Address of destination array descriptor
;
; Futures:
;
; The dsc$w_length field as assumed to be the length in bytes of the
; data items; thus, no check of the datatype is done.
; Note that we'd really like to vary the last index most rapidly.
; We may eventually special case lengths of 1,2,4 or 8.
; If so, and the smallest stride equals the length, autoincrement
; could be used in the loop.
; To get fancy, we could delete subscripts that don't really vary.
; Also, some adjacent subscripts could be combined.
; A user-callback could be used to handle data conversions.
;---
.sbttl LIB$MAT_INI Matrix Initialization
;---
;
; Calling sequence:
;
; CALL LIB$MAT_INI( [ src.rx.dr ], dest.wx.da )
;
; Inputs:
;
;
src = 4 ; Address of source item
dest = 8 ; Address of descriptor for destination matrix
;
; Description:
;
; This routine initializes all elements of a matrix with a given value.
; If the value is not specified, the elements are set to zero.
;
; The following data types are supported:
; W, L, Q, O, F, D, G, H
;
; Implementation:
;
; The array descriptor is converted to NCA descriptors, if needed.
;
; Returned status:
;
; ss$_normal Normal succesful completion
; lib$_invarg Unsupported datatype
; lib$_invarg Unsupported descriptor class
; lib$_invarg The array must have two dimensions
; lib$_invarg Error converting descriptor from class A to NCA
;
; Signalled status:
;
; ss$_opcdec Reserved opcode (only for data types O and H)
;
;---
.sbttl LIB$MAT_IDN Set Matrix to the Identity Matrix
;---
;
; Calling sequence:
;
; CALL LIB$MAT_IDN( dest.wx.da )
;
; Inputs:
;
dest = 4 ; Address of descriptor for destination matrix
;
; Description:
;
; This routine a square matrix to the identity matrix.
;
; The following data types are supported:
; W, L, Q, O, F, D, G, H
;
; Implementation:
;
; The array descriptor is converted to NCA descriptors, if needed.
;
; Returned status:
;
; ss$_normal Normal succesful completion
; lib$_invarg Unsupported datatype
; lib$_invarg Unsupported descriptor class
; lib$_invarg The array must have two dimensions
; lib$_invarg Error converting descriptor from class A to NCA
; lib$_matmussqu The matrix is not square
;
; Signalled status:
;
; ss$_opcdec Reserved opcode (only for data types O and H)
;
;---
.sbttl LIB$MAT_ADR Matrix Address
;---
;
; Calling sequence:
;
; CALL LIB$MAT_ADR( dsc.rx.da, idx.rl.v, ... )
;
; Inputs:
;
; dsc Address of array descriptor.
; idx One or more optional indices, one for each dimension.
;
; Outputs:
;
; R0 Address of the element in the array
;
; Description:
;
; This routine computes the address of an element in an array.
; If any error occurs, a zero is returned, unless it was due to
; an index out of bounds, in which case a subscript range trap occurs.
;
;---
.sbttl LIB$ABS Absolute value
;---
;
; Calling sequence:
;
; CALL LIB$ABS( arg.mx.dx )
;
; Inputs:
;
arg = 4 ; Argument, passed by descriptor
;
; Description:
;
; This routine replaces a value with its absolute value.
; The following data types are supported:
; W, L, F, D, G, H
;
; Returned status:
;
; ss$_normal Normal succesful completion
; lib$_invarg Unsupported datatype
;
; Signalled status:
;
; ss$_fltovf_f Floating overflow fault
; ss$_opcdec Reserved opcode (only for data types G and H)
;
;---
.sbttl LIB$SQRT Square root
;---
;
; Calling sequence:
;
; CALL LIB$SQRT( arg.mx.dx )
;
; Inputs:
;
arg = 4 ; Argument, passed by descriptor
;
; Description:
;
; This routine replaces a value with its square root.
; The following data types are supported:
; W, L, F, D, G, H
;
; Returned status:
;
; ss$_normal Normal succesful completion
; lib$_invarg Unsupported datatype
;
; Signalled status:
;
; ss$_fltovf_f Floating overflow fault
; ss$_opcdec Reserved opcode (only for data types G and H)
;
;---
|
| The header for LIB$MAT_INV (invert matrix in place) happened to be missing.
Here it is.
.sbttl LIB$MAT_INV Invert Matrix
;---
;
; Calling sequence:
;
; CALL LIB$MAT_INV( matrix.dx [, determinant.ax ] )
;
; Inputs:
;
; matrix(ap) Address of array descriptor (A or NCA)
; determinant(ap) Address of where to store determinant (optional)
; The determinant must be the same datatype as the array.
;
; Description:
;
; This routine inverts a matrix in place.
;
; The following data types are supported:
; W, L, F, D, G, H, FC, DC, GC, HC
; The following array classes are supported:
; A, NCA
;
; Returned status:
;
; ss$_normal Normal successful completion
; lib$_matmussqu The matrix is not square
; lib$_invarg Unsupported datatype
; lib$_matissing Matrix is singular
; lib$_invarg Unsupported descriptor class
; lib$_invarg The array must have two dimensions
; lib$_invarg Error converting descriptor from class A to NCA
;
; Signalled status:
;
; ss$_fltovf_f Floating overflow fault
; ss$_opcdec Reserved opcode (only for data types G and H)
;
; Note that it is possible for the determinant to overflow, even if
; the inverse matrix is correctly computed. However, the determinant
; is computed only if it is requested.
;
; Implementation Notes:
;
; The inversion algorithm is from COLLECTED ALGORITHMS FROM CACM,
; algorithms 230 and 231. It inverts a matrix in place using the
; Gauss-Jordan method with complete matrix pivoting.
;
; After some data checks and initialization, a copy of the descriptor
; is put on the stack (converted to NCA format). This is modified for
; ease of use.
;
; Because of the multiplicity of data types, macroes is used to expand
; into code for each supported data type.
; Since VAX-11 MACRO allows no support for symbolic registers, these
; macroes are also used to define names for symbolic registers (a few
; of which are actually stack locations).
;
; A complete description of the algorithm is given in comments in a
; Pascal-like language. Loops that do not depend on the order in which
; the loop index assumes values are written "for x in {lo..hi}".
;
; In the expansion of the heart1 macro (for each datatype), an "accuracy"
; parameter is included. This can be changed at compile-time.
; If this is zero, divides are done by an inverse and multiplies.
; If this is non-zero, divides are done by divides.
; Also, this parameter has an important effect on integer data types,
; described below.
;
; For integer data types, a naive implementation would result in a
; matrix consisting exclusively of ones and zeros, due to division
; being used for the multiplicative inverse. Instead, the inverse
; of an integer is computed; when a number is multiplied by it's
; inverse, the result equals one (mod 2^32). This integer inverse
; exists for all odd integers. This provides a (fairly) reasonable
; result, as the product of the matrix and its inverse (if it exists)
; equals the unit matrix. The inverse matrix exists if and only if
; the determinant of the original matrix is odd. Permutation matrices
; are correctly inverted.
; Note that this can be changed by use of the "accuracy" parameter.
;---
|