# $LU$ factorization

Our task is to implement a factorization routine that given A in a suitably efficient data structure returns the factors $L$ and $U$ where $L$ is unit lower triangular and $U$ is upper triangular. In order to guarantee that pivoting is not required you may restrict $A$ to be a diagonally dominant matrix. You must also write routines that solve systems $Lv = f$ and $Uv = f$ for arbitary vectors $v$ and $f$ and $L$ and $U$ matrices that result from your factorization routine. Your code must be as efficient as possible in terms of the number of comptuations and the storage required.

Let $A\in\mathbb{R}^{n\times n}$ be a nonsymmetric and nonsingular matrix with zero/nonzero element structure that has nonzero elements on the main diagonal, i.e., $\alpha_{i,j}\neq 0$, the first superdiagonal, i.e., $\alpha_{i,i+1}\neq 0$, the first subdiagonal, i.e., $\alpha_{i-1,i}\neq 0$, the fourth superdiagonal, i.e., $\alpha_{i,i+4}\neq 0$, and the fourth subdiagonal, i.e., $\alpha_{i-4,i}\neq 0$. For $n = 15$ this has the form $$A = \begin{pmatrix} * & * & 0 & 0 & * & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ * & * & * & 0 & 0 & * & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & * & * & * & 0 & 0 & * & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & * & * & * & 0 & 0 & * & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ * & 0 & 0 & * & * & * & 0 & 0 & * & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & * & 0 & 0 & * & * & * & 0 & 0 & * & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & * & 0 & 0 & * & * & * & 0 & 0 & * & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & * & 0 & 0 & * & * & * & 0 & 0 & * & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & * & 0 & 0 & * & * & * & 0 & 0 & * & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & * & 0 & 0 & * & * & * & 0 & 0 & * & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & * & 0 & 0 & * & * & * & 0 & 0 & * \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & * & 0 & 0 & * & * & * & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & * & 0 & 0 & * & * & * & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & * & 0 & 0 & * & * & * \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & * & 0 & 0 & * & * \\ \end{pmatrix}$$

Again, I have to perform the $LU$ factorization in as efficiently as possible, so far what I have done is store this matrix $A$ in a storage efficient manner. But I am not sure how to get $L$ and $U$ in the same manner, any suggestions is greatly appreciated. Here is my code so far:

            /*
* File:   main.cpp
* Author: Morgan Weiss
*
* Created on February 15, 2016, 8:21 PM
*/

#include <iostream>
#include <cmath>

double** banded_matrix(int n);
double get(double** A, int n, int i, int j);
void set(double** A, int n, int i, int j, double val);
void create(double** A, int n, int flag);

int main()
{
int N = 4;

//Initialize A
double **A = banded_matrix(N);
create(A,N,0);

//Initialize L
double L[N][N];
for (int i = 0; i < N; i++) {
for (int j = 0; j < N; j++){
L[i][j] = 0;
}
}

//Initialize U
double U[N][N];
for (int i = 0; i < N; i++) {
for (int j = 0; j < N; j++){
U[i][j] = 0;
}
}

//Set L[i][i] = 1
for (int i = 0; i < N; i++) {
L[i][i] = 1;
}

//Solve for L and U
for (int j = 0; j < N; j++) //loop over columns
{
//Compute U
for (int i = 0; i <= j; i++)
{
double sum = 0;
if (i > 0)
{
for (int k = 0; k < i; k++){
sum += L[i][k]*U[k][j];
}
}
U[i][j] = get(A,N,i,j) - sum;
}

//Compute L
for (int i = j+1; i < N; i++)
{
double sum = 0;
for (int k = 0; k < j; k++){
sum += L[i][k]*U[k][j];
}

L[i][j] = (get(A,N,i,j) - sum)/U[j][j];

}
}

//Initialize result matrix R
double R[N][N];
for (int i = 0; i < N; i++) {
for (int j = 0; j < N; j++){
R[i][j] = 0;
for (int k = 0; k < N; k++){
R[i][j] += L[i][k]*U[k][j];
}
}
}

//Print R
std::cout << "R = " << std::endl;
for (int i = 0; i < N; i++) {
for (int j = 0; j < N; j++) {
std::cout << R[i][j] << "   ";
}
std::cout << std::endl;
}
std::cout << std::endl;

//Print A
std::cout << "A = " << std::endl;
for (int i = 0; i < N; i++) {
for (int j = 0; j < N; j++) {
std::cout << get(A,N,i,j) << "   ";
}
std::cout << std::endl;
}
std::cout << std::endl;

//Print U
std::cout << "U = " << std::endl;
for (int i = 0; i < N; i++) {
for (int j = 0; j < N; j++) {
std::cout << U[i][j] << "   ";
}
std::cout << std::endl;
}
std::cout << std::endl;

//Print L
std::cout << "L = " << std::endl;
for (int i = 0; i < N; i++) {
for (int j = 0; j < N; j++) {
std::cout << L[i][j] << "   ";
}
std::cout << std::endl;
}
std::cout << std::endl;

return 0;
}

void create(double** A, int n, int flag) {
if (flag == 0) {
for (int i=0; i<n-1; i++) {
for (int j=0; j<n-1; j++) {
if (i==j) {
set(A,n,i,j,2.);
} else if (abs(i-j)==1 || abs(i-j)==4) {
set(A,n,i,j,1.);
}
}
}
} else if (flag==1) {}
}

double** banded_matrix(int n) {
double **data = new double *[n];
data = new double;
for(int i = 1; i <= 3; i++){
data[i] = new double;
}
for(int i = 4; i <= n-3; i++){
data[i] = new double;
}
for(int i = n-4; i <= n-2; i++){
data[i] = new double;
}
data[n-1] = new double;
return data;
}

double get(double** A, int n, int i, int j) {
double result;
// main tridiagonal
if(abs(i-j) < 2){
if(i == 0){
result = A[i][j];
} else if(i <= 3){
result = A[i][j-(i-2)];
} else if(i >= 3){
result = A[i][j-(i-1)];
}
// 4th superdiagonal
} else if (j - i == 4){
if(i == 0 || i == n-1) {
result = A[i];
} else if((i >= 1 && i <= 3) || (i >= n - 4 && i <= n - 2)) {
result = A[i];
} else {
result = A[i];
}
// 4th subdiagonal
} else if (i - j == 4){
result = A[i];
// the rest of the matrix
} else {
result = 0.;
}
return result;
}

void set(double** A, int n, int i, int j, double val) {
// main tridiagonal
if(abs(i-j) < 2){
if(i == 0){
A[i][j] = val;
} else if(i <= 3){
A[i][j-(i-2)] = val;
} else if(i >= 3){
A[i][j-(i-1)] = val;
}
// 4th superdiagonal
} else if (j - i == 4){
if(i == 0 || i == n-1) {
A[i] = val;
} else if((i >= 1 && i <= 3) || (i >= n - 4 && i <= n - 2)) {
A[i] = val;
} else {
A[i] = val;
}
// 4th subdiagonalf
} else if (i - j == 4){
A[i] = val;
// the rest of the matrix
} else {
std::cout << "cannot set element (" << i <<"," << j <<") in matrix" << std::endl;
}
}


This is my best attempt, if anyone has any suggestions on the way I implemented this, I would appreciate it.

• Just a brief comment - the matrices L and U can not have analogous format to A, e.g. L has lower triangular form of A, instead you should consider a band format, at least for L and U, say the non zero values can occur for L at L_{i,i}, L_{i,i-1}, L_{i,i-2}, L_{i,i-3}, L_{i,i-4}. – Peter Frolkovič Feb 17 '16 at 8:38

I am all for creating your own algorithm in order to do whatever. But I have to ask, why not use a linear algebra package written in C++ that can do LU factorization? I can see having to write your own if there are specific design requirements that the packages do not meet. If there aren't any special requirements, then why create more work for yourself?

Usually, these packages have gone through the fire and flames of testing to make them efficient. I will link a few to these in the answer:

http://math.nist.gov/tnt/jama_doxygen/class_JAMA__LU.html

http://www.gnu.org/software/gsl/

For some low level stuff - http://www.netlib.org/blas/

Here is a list of some more Linear Algebra packages:

https://en.wikipedia.org/wiki/Comparison_of_linear_algebra_libraries

• It's homework for a class (see the tag), so using one of these libraries would be cheating... – Christian Clason Feb 18 '16 at 7:18
• Oh, I see, I apologize. Best of luck to you. I am fairly new to the site. Even then, you can take a look at how the open source programs implemented the same thing to get a rough idea of how you will program yours. I would imagine that the two will be similar. There are only so many ways to can do the LU factorization. – philm Feb 18 '16 at 13:38
• Thank you for the suggestion but yes I need to write my own algorithm for it – Wolfy Feb 18 '16 at 18:20

Have you tried running your code? I think it will segfault. The routine "banded_matrix" twice sets the same values of data, and then in the 3rd loop sets a couple of non-existing elements.

• Yes, I did run the code. I get correct results, for when n = 4, 8 , and 15. I don't think it is necessary to make n very large – Wolfy Feb 22 '16 at 19:10