-
Notifications
You must be signed in to change notification settings - Fork 9
/
tridiagmatrix.R
68 lines (64 loc) · 2.39 KB
/
tridiagmatrix.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
## Copyright (c) 2016, James P. Howard, II <jh@jameshoward.us>
##
## Redistribution and use in source and binary forms, with or without
## modification, are permitted provided that the following conditions are
## met:
##
## Redistributions of source code must retain the above copyright
## notice, this list of conditions and the following disclaimer.
##
## Redistributions in binary form must reproduce the above copyright
## notice, this list of conditions and the following disclaimer in
## the documentation and/or other materials provided with the
## distribution.
##
## THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
## "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
## LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
## A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
## HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
## SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
## LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
## DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
## THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
## (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
## OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#' @title Solve a tridiagonal matrix
#'
#' @description
#' use the tridiagonal matrix algorithm to solve a tridiagonal matrix
#'
#' @param D vector of entries on the main diagonal
#' @param L vector of entries below the main diagonal
#' @param U vector of entries above the main diagonal
#' @param b vector of the right-hand side of the linear system
#'
#' @details
#' \code{tridiagmatrix} uses the tridiagonal matrix algorithm to solve a
#' tridiagonal matrix.
#'
#' @return the solution vector
#'
#' @family linear
#'
#' @export
tridiagmatrix <- function(L, D, U, b) {
n <- length(D)
L <- c(NA, L)
## The forward sweep
U[1] <- U[1] / D[1]
b[1] <- b[1] / D[1]
for(i in 2:(n - 1)) {
U[i] <- U[i] / (D[i] - L[i] * U[i - 1])
b[i] <- (b[i] - L[i] * b[i - 1]) /
(D[i] - L[i] * U[i - 1])
}
b[n] <- (b[n] - L[n] * b[n - 1]) /
(D[n] - L[n] * U[n - 1])
## The backward sweep
x <- rep.int(0, n)
x[n] <- b[n]
for(i in (n - 1):1)
x[i] <- b[i] - U[i] * x[i + 1]
return(x)
}