matSweep {lazy.symbolic}

R Documentation

Symbolic Matrix Sweep Operator

Description

Symbolic Matrix Sweep Operator

Usage

matSweep(A, loc = 1:nrow(A), simplify = 2, all = 0, func_level = 0, debug = 0)

Arguments

A	A square matrix
loc	vector of pivot locations
simplify	= 0 not to use Simplify or Expand = 1 to use Simplify once before exit. = 2 to use Simplify at the bottom of each pivot loop.
all	parameter for Simplify
func_level	parameter for Simplify
debug	= 1 to print details

Details

Definition of Sweep
C <- matSweep( A, p ) means:
c[i,j] <- c[i,j] - c[i,p]*c[p,j]/c[p,p]
c[p,] <- c[p,]/c[p,p]
c[,p] <- - c[,p]/c[p,p]

# Inversion of A.
matSweep( A, 1:ncol(A) ) == solve(A)
This is equivalent to Inv(A).
# Sweep is reversible.
matSweep( matSweep( A, loc ), loc ) == A

Use matSwp for numeric matrix.

Value

A matrix of the same size as A.

Examples

# inversion of 2 x 2 A
A <- demomat(2,2,root="a")
Ainv <- matSweep(A)
AAinv <- matTimes(A,Ainv)
Print( Eval(AAinv, A=demomat(2,2)), fuzz=1e-9 )
AinvA <- matTimes(Ainv,A)
Print( Eval(AinvA, A=demomat(2,2)), fuzz=1e-9 )

# inversion of 3 x 3 A
A <- demomat(3,3,root="a")
An <- demomat(3,3); An[1,3] <- 11  # make sure it it fullrank
# An <- demomat(3,3,shape="uppert")
Ainv <- matSweep(A)
AAinv <- matTimes(A,Ainv)
AAinvn <- matReplace(AAinv,A=An)
Print( Eval(AAinvn), fuzz=1e-9 )
AinvA <- matTimes(Ainv,A)
AinvAn <- matReplace(AinvA,A=An)
Print( Eval(AinvAn), fuzz=1e-9 )


# See if it is reversible.
A <- demomat(2,2,root="a")
A1 <- matSweep(A,1)
A11 <- matSweep(A1,1)
Print( Eval(A11,A=demomat(2,2)) )

# See if it is reversible.
A <- demomat(3,3,root="a")
A1 <- matSweep(A,1:2)
A11 <- matSweep(A1,1:2)
Print( Eval(A11,A=demomat(3,3)) )


# Fully Symbolic Partitioned Matrix
Ap=demomat(2,2,root="A", fullsymb=1)
Print(Ap)
invAp=matSweep(Ap)
Print(invAp)

# numeric checking
A11n=matrix(c(1,2,0,2),2,2); A12n=matrix(c(1,0,1,1),2,2)
A21n=matrix(c(2,0,1,1),2,2); A22n=matrix(c(1,1,0,2),2,2)
Apn=rbind(cbind(A11n,A12n),cbind(A21n,A22n))
invApn=solve(Apn)
# change inv to solve
invAps=gsub("inv","solve", invAp, fixed=1)
invApn2=matReplace( invAps, A11=A11n, A12=A12n, A21=A21n, A22=A22n )
Print(invApn-invApn2, fmt="9.6")

---