Download: matrix.zip.
See: ludcmp.cpp lusolve.cpp (test program).
LU decomposition is a method of solving linear matrix equations,
This program uses the Matrix class.
ludcmp.cpp/*\file ludcmp.cpp
*
*! Finds solution to set of linear equations A x = b by LU decomposition.
*
* Chapter 2, Programs 3-5, Fig. 2.8-2.10
* Gerald/Wheatley, APPLIED NUMERICAL ANALYSIS (fourth edition)
* Addison-Wesley, 1989
*
*/
#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
#include <math.h>
#include "matrix.h"
int ludcmp(Matrix &a, int order[]);
void solvlu(const Matrix &a, const Vec &b, Vec &x, const int order[]);
static int pivot(Matrix &a, int order[], int jcol);
#define TINY 1e-20
//! finds LU decomposition of Matrix
/*!
* The function ludcmp computes the lower L and upper U triangular
* matrices equivalent to the A Matrix, such that L U = A. These
* matrices are returned in the space of A, in compact form.
* The U Matrix has ones on its diagonal. Partial pivoting is used
* to give maximum valued elements on the diagonal of L. The order of
* the rows after pivoting is returned in the integer vector "order".
* This should be used to reorder the right-hand-side vectors before
* solving the system A x = b.
*
*
* \param a - n by n Matrix of coefficients
* \param order - integer vector holding row order after pivoting.
*
*/
int ludcmp(Matrix &a, int order[])
{
int i, j, k, n, nm1;
int flag = 1; /* changes sign with each row interchange */
double sum, diag;
n = a.rows();
assert(a.cols()==n);
/* establish initial ordering in order vector */
for (i=0; i<n; i++) order[i] = i;
/* do pivoting for first column and check for singularity */
if (pivot(a,order,0)) flag = -flag;
diag = 1.0/a[0][0];
for (i=1; i<n; i++) a[0][i] *= diag;
/*
* Now complete the computing of L and U elements.
* The general plan is to compute a column of L's, then
* call pivot to interchange rows, and then compute
* a row of U's.
*/
nm1 = n - 1;
for (j=1; j<nm1; j++) {
/* column of L's */
for (i=j; i<n; i++) {
sum = 0.0;
for (k=0; k<j; k++) sum += a[i][k]*a[k][j];
a[i][j] -= sum;
}
/* pivot, and check for singularity */
if (pivot(a,order,j)) flag = -flag;
/* row of U's */
diag = 1.0/a[j][j];
for (k=j+1; k<n; k++) {
sum = 0.0;
for (i=0; i<j; i++) sum += a[j][i]*a[i][k];
a[j][k] = (a[j][k]-sum)*diag;
}
}
/* still need to get last element in L Matrix */
sum = 0.0;
for (k=0; k<nm1; k++) sum += a[nm1][k]*a[k][nm1];
a[nm1][nm1] -= sum;
return flag;
}
//! Find pivot element
/*!
* The function pivot finds the largest element for a pivot in "jcol"
* of Matrix "a", performs interchanges of the appropriate
* rows in "a", and also interchanges the corresponding elements in
* the order vector.
*
*
* \param a - n by n Matrix of coefficients
* \param order - integer vector to hold row ordering
* \param jcol - column of "a" being searched for pivot element
*
*/
int pivot(Matrix &a, int order[], int jcol)
{
int i, ipvt,n;
double big, anext;
n = a.rows();
/*
* Find biggest element on or below diagonal.
* This will be the pivot row.
*/
ipvt = jcol;
big = fabs(a[ipvt][ipvt]);
for (i = ipvt+1; i<n; i++) {
anext = fabs(a[i][jcol]);
if (anext>big) {
big = anext;
ipvt = i;
}
}
assert(fabs(big)>TINY); // otherwise Matrix is singular
/*
* Interchange pivot row (ipvt) with current row (jcol).
*/
if (ipvt==jcol) return 0;
a.swaprows(jcol,ipvt);
i = order[jcol];
order[jcol] = order[ipvt];
order[ipvt] = i;
return 1;
}
//! This function is used to find the solution to a system of equations,
/*! A x = b, after LU decomposition of A has been found.
* Within this routine, the elements of b are rearranged in the same way
* that the rows of a were interchanged, using the order vector.
* The solution is returned in x.
*
*
* \param a - the LU decomposition of the original coefficient Matrix.
* \param b - the vector of right-hand sides
* \param x - the solution vector
* \param order - integer array of row order as arranged during pivoting
*
*/
void solvlu(const Matrix &a, const Vec &b, Vec &x, const int order[])
{
int i,j,n;
double sum;
n = a.rows();
/* rearrange the elements of the b vector. x is used to hold them. */
for (i=0; i<n; i++) {
j = order[i];
x[i] = b[j];
}
/* do forward substitution, replacing x vector. */
x[0] /= a[0][0];
for (i=1; i<n; i++) {
sum = 0.0;
for (j=0; j<i; j++) sum += a[i][j]*x[j];
x[i] = (x[i]-sum)/a[i][i];
}
/* now get the solution vector, x[n-1] is already done */
for (i=n-2; i>=0; i--) {
sum = 0.0;
for (j=i+1; j<n; j++) sum += a[i][j] * x[j];
x[i] -= sum;
}
}
lusolve.cpp/*\file lusolve.cpp
*
* Solve linear equations by LU decomposition
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "matrix.h"
#include <string>
using namespace std;
char text[128];
void print_vector(double *v, int n);
void lusolve(Matrix &mx, FILE *in);
int ludcmp(Matrix &a, int order[]);
void solvlu(const Matrix &a, const Vec &b, Vec &x, const int order[]);
int main(int argc, char *argv[])
{
Matrix mx;
FILE *in;
char *filename;
if (argc>1) {
filename = argv[1];
printf("opening: %s\n",filename);
}
else {
printf("This program requires a data file.\n");
printf("Enter filename: ");
gets_s(text,72);
if (*text==0 || *text=='\n') return -1;
filename = text;
}
errno_t nerr = fopen_s(&in, filename,"rt");
if (nerr) {
perror(filename);
return -1;
}
string title;
while (read_matrix(mx,title,in)==0) {
printf("\n%s\n",title.c_str());
printf("\nMatrix A\n");
mx.print();
lusolve(mx, in);
}
return 0;
}
void lusolve(Matrix &mx, FILE *in)
{
int i, j, k, n, flag, index, nrhs;
int *ipvt;
double sum;
Vec b, x, save;
double det;
/*
* Create matrix a as a copy of input matrix.
*/
assert(mx.rows()==mx.cols()); // input must be square matrix
n = mx.rows();
Matrix a = mx;
ipvt = ivector(n);
flag = ludcmp(a,ipvt);
/* Calculate determinant */
det = flag;
for (i=0; i<n; i++) det *= a[i][i];
printf("\ndeterminant: %g\n",det);
printf("pivot vector: ");
for (i=0; i<n; i++) printf("%3d",ipvt[i]);
printf("\n");
/* print LU Matrix */
Matrix lower = a;
Matrix upper = a;
Matrix product(n, n);
for (i=0; i<n; i++) {
for (j=i+1; j<n; j++) lower[i][j] = 0.0;
for (j=0; j<i; j++) upper[i][j] = 0.0;
upper[i][i] = 1.0;
}
printf("\nLower triangular Matrix\n");
lower.print();
printf("\nUpper triangular Matrix\n");
upper.print();
/* Multiply the lower and upper matrices */
for (i=0; i<n; i++) {
for (j=0; j<n; j++) {
sum = 0.0;
for (k=0; k<n; k++) sum += lower[i][k] * upper[k][j];
product[i][j] = sum;
}
}
printf("\n\nProduct of L U\n\n");
product.print();
/*
* Read in the right-hand side vectors.
*/
fscanf_s(in," %d",&nrhs);
fgets(text,72,in);
x.create(n);
b.create(n);
index = 0;
while (nrhs--) {
for (i=0; i<n; i++) fscanf_s(in," %lf", &b[i]);
if (feof(in)) break;
printf("\nRight-hand-side number %d\n\n",++index);
b.print();
save = b;
if (fabs(det)<1e-10) {
printf("\nCoefficient Matrix is singular.\n");
continue;
}
solvlu(a,b,x,ipvt);
printf("\nSolution vector\n\n");
x.print();
k = 0;
for (i=0; i<n; i++) {
sum = 0.0;
for (j=0; j<n; j++) sum += mx[i][j]*x[j];
if (fabs(sum-save[i])>1e-8) k++;
}
printf("\nsolution %s.\n",(k? "does not check": "checks"));
}
free(ipvt);
}
void print_vector(double *v, int n)
{
while (n--) printf("%10.4f", *v++);
printf("\n");
}
The following output was produced from the command:
lusolve eqn5.dat > output.txt
Here are the results:
opening: eqn5.dat
Example set of equations used in text, beginning on p. 113
Matrix A
4 -2 1
-3 -1 4
1 -1 3
determinant: -18
pivot vector: 0 1 2
Lower triangular matrix
4 0 0
-3 -2.5 0
1 -0.5 1.8
Upper triangular matrix
1 -0.5 0.25
0 1 -1.9
0 0 1
Product of L U
4 -2 1
-3 -1 4
1 -1 3
Right-hand-side number 1
15.0000 8.0000 13.0000
Solution vector
2.0000 -2.0000 3.0000
solution checks.
Example set of equations used in text, beginning on p. 120
Matrix A
0 2 0 1
2 2 3 2
4 -3 0 1
6 1 -6 -5
determinant: -234
pivot vector: 3 2 1 0
Lower triangular matrix
6 0 0 0
4 -3.667 0 0
2 1.667 6.818 0
0 2 2.182 1.56
Upper triangular matrix
1 0.1667 -1 -0.8333
0 1 -1.091 -1.182
0 0 1 0.8267
0 0 0 1
Product of L U
6 1 -6 -5
4 -3 0 1
2 2 3 2
0 2 0 1
Right-hand-side number 1
0.0000 -2.0000 -7.0000 6.0000
Solution vector
-0.5000 1.0000 0.3333 -2.0000
solution checks.
Test data used in Chapter 2, Program 2.1 See page 184 (Fig 2.4).
Matrix A
6 1 -6 -5
2 2 3 2
4 -3 0 1
0 2 0 1
determinant: 234
pivot vector: 0 2 1 3
Lower triangular matrix
6 0 0 0
4 -3.667 0 0
2 1.667 6.818 0
0 2 2.182 1.56
Upper triangular matrix
1 0.1667 -1 -0.8333
0 1 -1.091 -1.182
0 0 1 0.8267
0 0 0 1
Product of L U
6 1 -6 -5
4 -3 0 1
2 2 3 2
0 2 0 1
Right-hand-side number 1
6.0000 -2.0000 -7.0000 0.0000
Solution vector
-0.5000 1.0000 0.3333 -2.0000
solution checks.
Set of equations for problem 2-5 (page 201)
Matrix A
2 1 1 -2
4 0 2 1
3 2 2 0
1 3 2 0
determinant: 11
pivot vector: 1 3 2 0
Lower triangular matrix
4 0 0 0
1 3 0 0
3 2 -0.5 0
2 1 -0.5 -1.833
Upper triangular matrix
1 0 0.5 0.25
0 1 0.5 -0.08333
0 0 1 1.167
0 0 0 1
Product of L U
4 0 2 1
1 3 2 0
3 2 2 0
2 1 1 -2
Right-hand-side number 1
0.0000 8.0000 7.0000 3.0000
Solution vector
3.1818 2.3636 -3.6364 2.5455
solution checks.
Set of equations given for problem 2-11 (page 201)
Matrix A
3 2 -1 -4
1 -1 3 -1
2 1 -3 0
0 -1 8 -5
determinant: -1.33227e-014
pivot vector: 0 1 3 2
Lower triangular matrix
3 0 0 0
1 -1.667 0 0
0 -1 6 0
2 -0.3333 -3-4.441e-016
Upper triangular matrix
1 0.6667 -0.3333 -1.333
0 1 -2 -0.2
0 0 1 -0.8667
0 0 0 1
Product of L U
3 2 -1 -4
1 -1 3 -1
0 -1 8 -5
2 1 -3 0
Right-hand-side number 1
10.0000 -4.0000 16.0000 3.0000
Coefficient matrix is singular.
Maintained by John Loomis, updated Wed Feb 07 23:01:57 2007