[R] writing spdiags function for R
moreno
mcoco at staffmail.ed.ac.uk
Fri Aug 30 16:46:10 CEST 2013
I am posting here the brilliant solutions, gently provided by Prof JC Nash
<nashjc at uottawa.ca>, to me; so that people struggling in the future with
the same issue can find a way through.
FYI, compared to the original Matlab implementation: 1) "it does not handle
the case with more than one input, and
2) (m > n) matrices give the B matrix columns in a different order, but the
d vector of indices will also be changed
accordingly, so the set of columns is OK, just ordered differently" cif. JN.
Copied below the .R version of the spdiags code, a Fortran implementation of
it, and an R wrapper to run the .f
thanks again JN, your help was really invaluable :)
#########################################################################################
R version
#########################################################################################
spdiagsj <- function(A) { # A is a matrix
m <- dim(A)[[1]]
n <- dim(A)[[2]]
k <- min(m, n) # length of diagonals
Bdata<-NULL # start with nothing in B matrix (as vector)
jb<-0 # column index of "last" column saved for B
d<-NULL # index vector of diagonals from A
# d contains 0 for the principal diagonal, -i for i'th lower
# diagonal (prefaced with zeros), +j for j'th upper diagonal
# (suffixed by zeros)
q<-(m-1)+n # There are m-1 subdiagonals and n-1 superdiagonals + main
diagonal
if (m > n) { # tall matrix
Adata <- as.vector(t(A)) # convert to vector BY ROWS
} else { # fat or square matrix (m <= n)
Adata <- as.vector(A) # convert to vector BY COLUMNS
}
# Augment the data with columns of zeros fore and aft
Adata<-c(rep(0,(k-1)*k), Adata, rep(0,(k-1)*k))
cat("Augmented Adata with ",length(Adata)," elements:\n")
print(Adata)
for (i in 1:q) {
tv <- c(Adata[[k*(i-1)+1]], rep(0,k-1)) # top element of augmented
column
# plus enough zeros to pad it out (some zeros may be replaced
below)
# i.e., start at 1, then m+1, 2*m+1 etc.
qx<-min((q-i), (k-1)) #
cat(" qx=",qx,"\n")
if (qx > 0) { # qx will be 0 when we are at last
superdiagonal,i.e., i == q
for (j in 1:qx) { # get the rest of the diagonal elements
tv[[j+1]] <- Adata[[k*(i+j-1)+j+1]]
}
}
if (any(tv != 0)){ # check for non-zeros, if there are, then save
jb<-jb+1 # next column of B
d<-c(d,(i-k)) # record the index
Bdata<-c(Bdata, tv) # save the diagonal as column of B in
vector form
}
}
if (m > n) d <- -d # reset index
cat("Bdata:"); print(Bdata)
B <- matrix(Bdata, nrow=k, byrow=FALSE) # convert to matrix form
result<-list(B=B, d=d)
}
cat("Matlab example 1\n")
dta <- c(0, 5, 0, 10, 0, 0, 0, 0, 6, 0, 11, 0, 3, 0, 0, 7, 0, 12, 1, 4, 0,
0, 8, 0, 0, 2, 5, 0, 0, 9)
A1 <- matrix(dta, nrow=5, ncol=6, byrow=TRUE)
print(A1)
res1<-spdiagsj(A1)
print(res1)
tmp<-readline("Next")
cat("Matlab example 2\n")
n<-10 # choose 10 for an example
A2<-matrix(rep(0, n*n), nrow=n, ncol=n)
for (i in 1:n) {
for (j in 1:n) {
if (i == j) A2[i, j] <- -2
if ( (i == (j-1)) || (i == (j+1))) A2[i,j] <- 1
}
}
print(A2)
res2<-spdiagsj(A2)
print(res2)
tmp<-readline("Next")
cat("Matlab example 3\n")
dta3 <- c(11, 0, 13, 0, 0, 22, 0, 24, 0, 0, 33, 0, 41, 0, 0, 44,
0, 52, 0, 0, 0, 0, 63, 0, 0, 0, 0, 74)
A3 <- matrix(dta3, nrow=7, ncol=4, byrow=TRUE)
print(A3)
res3<-spdiagsj(A3)
print(res3)
tmp<-readline("Next")
cat("try transpose\n")
A3T<-t(A3)
print(A3T)
res3T<-spdiagsj(A3T)
print(res3T)
tmp<-readline("Next")
cat("Example 5B \n")
dta5b1<-c(6, 0, 13, 0, 0, 0, 7, 0, 14, 0, 1, 0, 8, 0, 15, 0, 2, 0, 9, 0, 0,
0, 3, 0, 10)
A5b1<-matrix(dta5b1, nrow=5, ncol=5, byrow=TRUE)
print(A5b1)
res5b1<-spdiagsj(A5b1)
print(res5b1)
#########################################################################################
Fortran version
#########################################################################################
subroutine jspd(m, n, k, Adata, jb, Bdata, d, tv, na, nb, nd)
C Central part of spdiags for R
C m and n are row and column sizes of A (underlying matrix)
C jb will be number of returned diagonals
C returns jb, Bdata, d
integer m, n, na, nb, nd, jb, d(nd)
integer i, j, k, kend, q, mn, kk1, js, je, qx
double precision Adata(na), Bdata(nb), tv(k)
LOGICAL not0
C k = min(m, n)
C ?? check if k=1
C Bdata<-NULL # start with nothing in B matrix (as vector)
jb = 0
kend=0
C column index of "last" column saved for B
C d = NULL # index vector of diagonals from A
C # d contains 0 for the principal diagonal, -i for i'th lower
C # diagonal (prefaced with zeros), +j for j'th upper diagonal
C # (suffixed by zeros)
q = (m-1)+n
C There are m-1 subdiagonals and n-1 superdiagonals + main diagonal
C assume we have already built Adata for tall or fat matrix
C # Augment the data with columns of zeros fore and aft
mn = m*n
C print *,"Original Adata"
C print 1000, (Adata(i), i=1,mn)
do 10 i=1,mn
Bdata(i) = Adata(i)
10 continue
kk1 = (k-1)*k
do 15 i=1,kk1
Adata(i)=0.0
15 continue
js = kk1+1
je = kk1+mn
do 20 i=js,je
Adata(i) = Bdata(i-js+1)
20 continue
js = je+1
je = je+kk1
do 25 i=js,je
Adata(i) = 0.0
25 continue
C print *, "Augmented Adata with ", je," elements"
C print 1000,(Adata(i), i=1,je)
1000 FORMAT(1H ,25f4.0)
do 100 i=1,q
qx = k*(i-1)+1
C print *,"Top element is no ",qx,"=",Adata(qx)
tv(1) = Adata(qx)
do 30 j=1,(k-1)
tv(j+1) = 0.0
30 continue
qx = min(q-i, k-1)
C print *,"qx=",qx
if (qx .gt. 0) then
C qx will be 0 when we are at last superdiagonal,i.e., i == q
do 35 j=1,qx
tv(j+1) = Adata(k*(i+j-1)+j+1)
35 continue
endif
not0 = .FALSE.
do 40 j=1,k
if (tv(j) .ne. 0.0) not0 = .TRUE.
40 continue
C if (.NOT. not0) print *," Zero column for i=",i
if (not0) then
C print *,"Nonzero column for i=",i
jb = jb+1
d(jb) = (i-k)
C print *, " New jb=",jb," d element =",d(jb)
C print 1000,(tv(j), j=1,k)
do 45 j=1,k
Bdata(kend+j) = tv(j)
45 continue
kend = kend + k
C save the diagonal as column of B in vector form
endif
100 continue
C print *,"Bdata ", jb, " columns"
C qx = k*jb
C print 1000,(Bdata(j), j=1,qx)
return
end
#########################################################################################
R Wrapper
#########################################################################################
spDIAGS <- function(A) { # A is a matrix
m <- dim(A)[[1]]
n <- dim(A)[[2]]
k <- min(m, n) # length of diagonals
if (m > n) { # tall matrix
Adata <- as.vector(t(A)) # convert to vector BY ROWS
} else { # fat or square matrix (m <= n)
Adata <- as.vector(A) # convert to vector BY COLUMNS
}
# print(Adata)
la<-length(Adata)
na<-m*n+2*k*(k-1)
Adata <- c(Adata, rep(0,(na-la)))
nb<-(m+n)*k
nd<-m+n-1
jb<-0
Bdata<-rep(0,nb)
d<-rep(0,nd)
tv<-rep(0,k)
tres<-.Fortran("jspd",m=as.integer(m), n=as.integer(n), k=as.integer(k),
Adata=as.double(Adata), jb=as.integer(jb),
Bdata=as.double(Bdata), d=as.integer(d), tv=as.double(tv),
na=as.integer(na), nb=as.integer(nb), nd=as.integer(nd) )
# print(str(tres))
jb<-tres$jb
d<-tres$d[1:jb]
Bdata<-tres$Bdata[1:(jb*k)]
# print(Bdata)
if (m > n) d <- -d # reset index
B <- matrix(Bdata, nrow=k, byrow=FALSE) # convert to matrix form
result<-list(B=B, d=d)
}
#########################################################################################
--
View this message in context: http://r.789695.n4.nabble.com/writing-spdiags-function-for-R-tp4552626p4675010.html
Sent from the R help mailing list archive at Nabble.com.
More information about the R-help
mailing list