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
Will - Take a look at Roger Koenker's package SparseMatrix, available from CRAN. Look also for some other package from Roger which depends on SparseMatrix, but has a different name. It's a place to look. I don't recall whether it will answer your need or not. - tom blackwell - u michigan medical school - ann arbor - On Wed, 1 Oct 2003, Will Harvey wrote:> 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 >
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
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
_______________________________________________
No banners. No pop-ups. No kidding.
Introducing My Way - http://www.myway.com
I see that the LAPACK routine DGTSV is in the R source, and defined in R_ext/lapack.h, but I don't know how to get to it from R. Andy> -----Original Message----- > From: Roger Koenker [mailto:roger at ysidro.econ.uiuc.edu] > Sent: Wednesday, October 01, 2003 10:33 AM > To: Thomas W Blackwell > Cc: r-help at stat.math.ethz.ch > Subject: Re: [R] Solving a tridiagonal system > > > SparseM is really intended for arbitrary sparse structure, > for banded structural there are much more efficient methods, > some of which are, if I'm not mistaken, now available in lapack. > > > url: www.econ.uiuc.edu/~roger/my.html Roger Koenker > email rkoenker at uiuc.edu Department of Economics > vox: 217-333-4558 University of Illinois > fax: 217-244-6678 > Champaign, IL 61820 > > On Wed, 1 Oct 2003, Thomas W Blackwell wrote: > > > Will - > > > > Take a look at Roger Koenker's package SparseMatrix, > available from > > CRAN. Look also for some other package from Roger which depends on > > SparseMatrix, but has a different name. It's a place to look. I > > don't recall whether it will answer your need or not. > > > > - tom blackwell - u michigan medical school - ann arbor - > > > > On Wed, 1 Oct 2003, Will Harvey wrote: > > > > > 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 > > > > ______________________________________________ > R-help at stat.math.ethz.ch mailing list > https://www.stat.math.ethz.ch/mailman/listinfo/r-help >______________________________________________ R-help at stat.math.ethz.ch mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help