[R] Solving a tridiagonal system

[Will Harvey]

Thanks. This works, but still seems too slow for my
application. It looks like my best option is to figure
out how to interface the LAPACK routine DGTSV to R.
The inputs to DGTSV are four vectors: the three
diagonals and the RHS vector.

Will


--- Gabor Grothendieck <ggrothendieck at myway.com>
wrote:
> 
> I just noticed that you defined a to be the lower
> diagonal
> whereas I had it as the upper diagonal so the
> previous email
> should be as follows to correspond to your notation:
> 
> 
> # define test vectors for a, b and c
> va <- -(1:4); vb <- 11:15; vc <- 1:4
> 
> # diag.num is a matrix whose ith diagonal equals i
> diag.num <- -outer(seq(vb),seq(vb),"-")
> 
> x <- diag(vb)
> x[diag.num == -1] <- va
> x[diag.num == 1] <- vc
> 
> 
> From: Gabor Grothendieck <ggrothendieck at myway.com>
> To: <will_harvey03 at yahoo.com>,
> <r-help at stat.math.ethz.ch> 
> Subject: RE: [R] Solving a tridiagonal system 
> 
>  
>  
> 
> The following will create a matrix x with given sub
> diagonal,
> diagonal and super diagonal.
> 
> # define test vectors for a, b and c
> va <- -(1:4); vb <- 11:15; vc <- 1:4
> 
> # diag.num is a matrix whose ith super diagonal
> equals i and 
> # sub diagonal equals -i
> diag.num <- -outer(seq(vb),seq(vb),"-")
> 
> x <- diag(vb)
> x[diag.num == 1] <- va
> x[diag.num == -1] <- vc
> 
> ---
> 
> Date: Wed, 1 Oct 2003 05:56:34 -0700 (PDT) 
> From: Will Harvey <will_harvey03 at yahoo.com>
> To: <r-help at stat.math.ethz.ch> 
> Subject: [R] Solving a tridiagonal system 
> 
> I need to find solutions to a tridiagonal system. By
> this I mean a set of linear equations Ax = d where A
> is a square matrix containing elements A[i,i-1],
> A[i,i] and A[i,i+1] for i in 1:nrow, and zero
> elsewhere. R is probably not the ideal way to do
> this,
> but this is part of a larger problem that requires
> R.
> 
> In my application it is much easier (and much
> faster)
> to generate the diagonal and off-diagonal elements
> of
> A as vectors, i.e. a = A[i,i-1], b = A[i,i] and c =
> A[i,i+1]. So I have three vectors that define A,
> along
> with a solution vector d. The conventional method of
> solving such systems is to use the so-called "Thomas
> algorithm", see e.g. 
>
<http://www.enseeiht.fr/hmf/travaux/CD0001/travaux/optmfn/hi/01pa/hyb74/node24.html>;;.
> This is very easy to code, but much more difficult
> to
> "vectorize". Is anyone aware of a library that
> contains a fast implementation of this algorithm?
> 
> Another alternative is to use backsolve. I can
> easily
> eliminate the lower diagonal a, but I'm still left
> with b and c, whereas backsolve requires a matrix.
> Again, I can write a function to read b and c into a
> matrix, but this requires loops, and is too slow. Is
> there a vectorized way of doing it? Of course, the
> diag command works for b, but what about c? In
> Octave,
> diag allows for an offset, but R apparently does
> not.
> 
> I would appreciate any and all assistance you
> experts
> can offer. Thanks in advance.
> 
> Will Harvey
> 
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
>
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
> 
> 
> 
> 
> message or move to Select a Folder----- Folders
> ------InboxDraftsSentTrashBulk Mail---- My Folders
> ----R
> Back to Inbox | <Prev Next> 
> 
> As AttachmentAs Inline Text 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> Make My Way Your Home Page | Spread the Word
> 
> My Settings: Overview | Search | Email | Chat |
> Portfolio | Calendar | Groups | Profile 
> 
> 
> IMPORTANT: We do not present our users with pop-ups,
> banners or any other non-text advertising. Nor do we
> send email
> to our users. If you see or receive one of these
> items, it is coming from an outside source, either
> as a result of something you
> have previously downloaded or as an "exit" pop-up
> from the site you just visited. It is not coming
> from our site.
> 
> Privacy Policy Terms of Service Partner with Us Our
> Mission Advertise with Us Sign In Sign Out Help
> Center
> 
> © 2002-2003 My Way
> 
> 
> 
> 
> 
> 
> 
> 
> 
> _______________________________________________
> No banners. No pop-ups. No kidding.
> Introducing My Way - http://www.myway.com
> 
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
>
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
> 
>  
>  
>  
>  message or move to Select a Folder----- Folders
> ------InboxDraftsSentTrashBulk Mail---- My Folders
> ----R
>  Back to Inbox | <Prev Next> 
>  
>    As AttachmentAs Inline Text 
>  
>  
> 
> 
> 
>     
>  
>  
> 
> 
> 
> 
> Make My Way Your Home Page  |  Spread the Word
> 
> My Settings: Overview | Search | Email | Chat |
> Portfolio | Calendar | Groups | Profile 
>  
> 
> IMPORTANT: We do not present our users with pop-ups,
> banners or any other non-text advertising. Nor do we
> send email
> to our users. If you see or receive one of these
> items, it is coming from an outside source, either
> as a result of something you
> have previously downloaded or as an "exit" pop-up
> from the site you just visited. It is not coming
> from our site.
> 
> Privacy Policy   Terms of Service   Partner with Us 
>  Our Mission   Advertise with Us   Sign In   Sign
> Out   Help Center
> 
> © 2002-2003 My Way
> 
>  
>  
> 
=== message truncated ===