--* From BMT%WATSON.vnet.ibm.com@yktvmh.watson.ibm.com Mon Oct 18 17:41:41 1993 --* Received: from yktvmh.watson.ibm.com by radical.watson.ibm.com (AIX 3.2/UCB 5.64/900524) --* id AA21766; Mon, 18 Oct 1993 17:41:41 -0400 --* Received: from watson.vnet.ibm.com by yktvmh.watson.ibm.com (IBM VM SMTP V2R3) --* with BSMTP id 4987; Mon, 18 Oct 93 17:48:09 EDT --* Received: from YKTVMH by watson.vnet.ibm.com with "VAGENT.V1.0" --* id --* for asbugs@watson; Mon, 18 Oct 93 17:48:09 -0400 --* Received: from YKTVMH by watson.vnet.ibm.com with "VAGENT.V1.0" --* id 2831; Mon, 18 Oct 1993 17:48:07 EDT --* Received: from cyst.watson.ibm.com by yktvmh.watson.ibm.com (IBM VM SMTP V2R3) --* with TCP; Mon, 18 Oct 93 17:48:07 EDT --* Received: from spadserv.watson.ibm.com by cyst.watson.ibm.com (AIX 3.2/UCB 5.64/900528) --* id AA50349; Mon, 18 Oct 1993 17:47:37 -0400 --* Received: by spadserv.watson.ibm.com (AIX 3.2/UCB 5.64/900524) --* id AA23754; Mon, 18 Oct 1993 17:49:35 -0400 --* Date: Mon, 18 Oct 1993 17:49:35 -0400 --* From: bmt@spadserv.watson.ibm.com --* X-External-Networks: yes --* Message-Id: <9310182149.AA23754@spadserv.watson.ibm.com> --* To: asbugs@watson.ibm.com --* Subject: compiler gets segmentation violation [matopfld.as][32.1 (current)] --@ Fixed by: PI Sun Mar 20 20:33:21 EST 1994 --@ Tested by: none --@ Summary: no more a bug #include "aslib.as" #library asdem "asdem" import from asdem Vector ==> Array SI ==> SingleInteger MatrixOpF(R:Field) : MO == Definition where Col ==> Vector R Mat ==> Matrix R MO ==> with inverse : Mat -> Partial Mat ++ inverse(m) fails if m is not invertible / : (Mat,R) -> Mat ++ m/r divide all the entries by r rowEchelon: Mat -> Mat ++ \spad{rowEchelon(m)} returns the row echelon form of the matrix m. rank: Mat -> SingleInteger ++ \spad{rank(m)} returns the rank of the matrix m. nullity: Mat -> SingleInteger ++ \spad{nullity(m)} returns the nullity of the matrix m. This is ++ the dimension of the null space of the matrix m. nullSpace: Mat -> List Col ++ \spad{nullSpace(m)}+ returns a basis for the null space of ++ the matrix m. determinant : Mat -> R ++ determinant of m, error if m is not square Definition ==> add import from SingleInteger import from R import from MatMappings R rowAllZeroes?: (Mat,SingleInteger) -> Boolean colAllZeroes?: (Mat,SingleInteger) -> Boolean rowAllZeroes?(x:Mat,i:SingleInteger) : Boolean == -- determines if the ith row of x consists only of zeroes -- internal function: no check on index i for j in 1..nCols x repeat ~zero? x(i,j) => return false true colAllZeroes?(x:Mat,j:SingleInteger) : Boolean == -- determines if the ith column of x consists only of zeroes -- internal function: no check on index j for i in 1..nRows x repeat zero? x(i,j) => return false true (x:Mat) / (r:R) : Mat == map((s:R):R +-> s/r,x) positivePower(x:Mat,n:SI) : Mat == n=1 => x odd? n => x * positivePower(x,n - 1) y := positivePower(x,n quo 2) y * y (x:Mat) ** (n:SingleInteger) : Mat == not((nr:= nRows x) = nCols x) => error "**: matrix must be square" n=0 => scalarMatrix(nr,1) positivePower(x,n) inverse(x:Mat):Partial Mat == (ndim := nRows x) ~= (nCols x) => error "inverse: matrix must be square" ndim = 2 => ans2 : Mat := zero(ndim, ndim) det := x(1,1)*x(2,2)-x(1,2)*x(2,1) if det = 0 then return failed detinv:=inv det ans2(1,1) := x(2,2)*detinv ans2(1,2) := -x(1,2)*detinv ans2(2,1) := -x(2,1)*detinv ans2(2,2) := x(1,1)*detinv ans2 :: Partial Mat AB : Mat := zero(ndim,ndim + ndim) for i in 1..ndim repeat for j in 1..ndim repeat AB(i,j) := x(i,j) AB(i,i+ndim):= 1 AB := rowEchelon AB zero? AB(ndim,ndim) => failed subMatrix(AB,1,ndim,ndim+1,ndim + ndim) ::Partial Mat rowEchelon(y : Mat) : Mat == -- row echelon form via Gaussian elimination x : Mat := copy y nr:= nRows x nc:= nCols x i : SI := 1 local n : SingleInteger for j in 1..nc repeat i > nr => return x n := 0 -- n = smallest k such that k >= i and x(k,j) ~= 0 for k in i..nr repeat if ~zero? x(k,j) then (n := k; break) zero? n => iterate -- "no non-zeroes" -- put nth row in ith position if i ~= n then swapRows!(x,i,n) -- divide ith row by its first non-zero entry b := inv x(i,j) x(i,j):=1 for k in (j+1)..nc repeat x(i,k) := b * x(i,k) -- perform row operations so that jth column has only one 1 for k in 1..nr repeat if k ~= i and ~zero? x(k,j) then for k1 in (j+1)..nc repeat x(k,k1):= x(k,k1) - x(k,j) * x(i,k1) x(k,j) := 0 -- increment i i := i + 1 x rank(x : Mat) : SingleInteger == y := (rk := nRows x) > (rh := nCols x) => rk := rh transpose x copy x y := rowEchelon y; i := rk while rk > 0 and rowAllZeroes?(y,i) repeat i := i - 1 rk := (rk - 1) rk nullity(x:Mat) : SingleInteger == nCols x - rank x nullSpace(y:Mat) : List Col == x := rowEchelon y nrow := nRows x ncol := nCols x basis : List Col := [] rk := nrow -- compute rank = # rows - # rows of all zeroes while rk > 0 and rowAllZeroes?(x,rk) repeat rk := (rk - 1) -- if maximal rank, return zero vector import from Col ncol <= nrow and rk = ncol => [new(ncol,0)] -- if rank = 0, return standard basis vectors zero? rk => for j in 1..ncol repeat w : Col := new(ncol,0) w.j:= 1 basis := cons(w,basis) basis -- v contains information about initial 1's in the rows of x -- if the ith row has an initial 1 in the jth column, then -- v.j = i; v.j = 0, otherwise v : Vector SingleInteger := new(ncol,0) for i in 1..rk repeat for j in 1..ncol repeat zero? x(i,j) => iterate v.j:=1 break j := ncol while j > 0 repeat w : Col := new(ncol,0) -- if there is no row with an initial 1 in the jth column, -- create a basis vector with a 1 in the jth row if zero? v.j then colAllZeroes?(x,j) => w.j := 1 basis := cons(w,basis) for k in 1..(j-1) repeat if ~zero? v.k then w.k := -x(v.k,j) w.j := 1 basis := cons(w,basis) j := j - 1 basis determinant(y : Mat) : R == (ndim := nRows y) ~= (nCols y) => error "determinant: matrix must be square" -- Gaussian Elimination ndim = 1 => y(1,1) x := copy y ans : R := 1 ndim1:=ndim-1 for i in 1..ndim1 for j in 1..ndim1 repeat if zero? x(i,j) then rown : SI := 0 for k in (i+1)..ndim repeat ~zero? x(k,j) => (rown := k;break) if zero? rown then return 0 swapRows!(x,i,rown); ans := -ans ans := x(i,j) * ans; b := -inv x(i,j) for l in (j+1)..ndim repeat x(i,l) := b * x(i,l) for k in (i+1)..ndim repeat if ~zero? (b := x(k,j)) then for l in (j+1)..ndim repeat x(k,l):= x(k,l) + b * x(i,l) x(ndim,ndim) * ans