ken@turtleva.UUCP (Ken Turkowski) (08/19/83)
echo x - matrix_solve
mkdir matrix_solve
echo x - matrix_solve/gausselim.3
sed 's/^x//' >matrix_solve/gausselim.3 <<'!Funky!Stuff!'
x.TH GAUSSELIM 3
x.SH NAME
xgausselim -
xsolution of simultaneous linear equations by Gaussian elimination
x.SH SYNOPSIS
x.B gausselim(A, x, b, n)
x.br
xint
x.B n;
x.br
xdouble
x.B *A, *x, *b;
xint n;
x.SH DESCRIPTION
x.I Gausselim
xsolves simultaneous linear equations
xusing Gaussian elimination
xwith the robust techniques
xof
x.I partial pivoting
xand
x.I row equilibration.
xSuccessful solution is indicated by a return value of 1,
xwhereas ill-conditioning,
xsuch as singularity,
xis indicated with a return value of 0.
x.LP
xThe equation to be solved is
x.RS
xAx = b
x.RE
xwhere
x.B A
xis the (n x n) coefficient matrix,
x.B b
xis the driving column vector,
xand
x.B x
xis the solution vector.
x.LP
xBy C language convention,
xthe ordering of the
x.B A
xmatrix
xis row-sequential,
xi.e.
xthe matrix entries are ordered thusly:
xA[0][0],
xA[0][1],
xA[0][2], ... ,
xA[1][0],
xA[1][1], ... ,
xA[n-1][n-1].
xNote that this is the opposite convention
xassumed by Fortran.
x.LP
xThe driving vector
x.B b
xand the solution vector
x.B x
xare considered to be (1 x n) column vectors.
xIf row vectors are used,
xthe
x.B A
xmatrix must first be transposed.
xEquivalently,
xthe line containing
x.RS
x# define Ael(i,j) A[(i*n)+j] /* column vectors */
x.RE
xcould be changed to
x.RS
x# define Ael(i,j) A[(j*n)+i] /* row vectors */
x.RE
x.SH EXAMPLE
x# define N 4
x.br
xdouble A[N][N], x[N], b[n];
x.br
xgetmatrix(A, N);
x.br
xgetvector(x, N);
x.br
xgausselim(A, x, b, N);
x.br
xputsolution(x, N);
x.SH AUTHOR
xKen Turkowski,
xCADLINC, Inc.
x.SH SEE ALSO
xlu (1)
x.SH BUGS
xAlthough the algorithm is very robust
x(in terms of being able to handle "difficult" matrices),
xfull pivoting would result in a more robust solution.
x.LP
xThere is no computation of the magnitude of the residue
x(i.e. error)
xof the result,
xalthough this is relatively straightforward to implement.
!Funky!Stuff!
echo x - matrix_solve/lu.3
sed 's/^x//' >matrix_solve/lu.3 <<'!Funky!Stuff!'
x.TH LU 3
x.SH NAME
xlu_decompose, lu_solve -
xsolution of simultaneous linear equations by LU decomposition
x.SH SYNOPSIS
x.B lu_decompose(a, n)
x.br
xdouble
x.B *a;
x.br
xint
x.B n;
x.LP
x.B lu_solve(x, b, n)
x.br
xdouble
x.B *x, *b;
xint
x.B n;
x.SH DESCRIPTION
xThis module contains two routines
xwhich work together
xto solve simultaneous linear equations.
xIt significantly reduces the computation involved
xwhen multiple equations are to be solved
xusing the same coefficient matrix
xbut different driving vectors.
xIt is particularly well suited for matrix inversion
xby sequentially solving for a sequence of driving vectors
xwith all 0's except one 1.
x.LP
x.I LU decomposition
xof the
x.I (n x n)
xmatrix
x.B A
xis performed by the routine
x.I lu_decompose()
xusing the robust techniques
xof
x.I partial pivoting
xand
x.I row equilibration.
xSuccessful solution is indicated by a return value of 1,
xwhereas ill-conditioning,
xsuch as singularity,
xis indicated with a return value of 0.
x.LP
xOnce the
x.B A
xmatrix has been decomposed,
xthe routine
x.I lu_solve()
xmay be called several times with different driving vectors.
x.LP
xThe equation to be solved is
x.RS
x.B Ax
x=
x.B b
x.RE
xwhere
x.B A
xis the (n x n) coefficient matrix,
x.B b
xis the driving column-vector,
xand
x.B x
xis the solution column-vector.
x.LP
xBy C language convention,
xthe ordering of the
x.B A
xmatrix
xis row-sequential,
xi.e.
xthe matrix entries are ordered thusly:
xA[0][0],
xA[0][1],
xA[0][2], ... ,
xA[1][0],
xA[1][1], ... ,
xA[n-1][n-1].
xNote that this is the opposite convention
xassumed by Fortran.
x.LP
xThe driving vector
x.B b
xand the solution vector
x.B x
xare considered to be (1 x n) column vectors.
xIf row vectors are used,
xthe
x.B A
xmatrix must first be transposed.
xEquivalently,
xthe line containing
x.RS
x# define ael(i,j) a[i*n+j] /* column vectors */
x.RE
xcould be changed to
x.RS
x# define ael(i,j) a[j*n+i] /* row vectors */
x.RE
x.SH EXAMPLE
x# define N 4
x.br
xdouble A[N][N], x[N], b[n];
x.br
xlu_decompose(A, N);
xgetmatrix(A, N);
x.br
xwhile (getvector(x, N)) {
x.RS
xlu_solve(x, b, N);
x.br
xputsolution(x, N);
x.RE
x}
x.SH AUTHOR
xKen Turkowski,
xCADLINC, Inc.
x.SH SEE ALSO
xlu (1)
x.SH BUGS
xAlthough the algorithm is very robust
x(in terms of being able to handle "difficult" matrices),
xfull pivoting would result in a more robust solution.
x.LP
xThere is no computation of the magnitude of the residue
x(i.e. error)
xof the result,
xalthough this is relatively straightforward to implement.
!Funky!Stuff!
echo x - matrix_solve/matrixinv.3
sed 's/^x//' >matrix_solve/matrixinv.3 <<'!Funky!Stuff!'
x.TH MATRIXINV 3
x.SH NAME
xmatrixinv -
xinvert a matrix
x.SH SYNOPSIS
x.B matrixinv(A, Ainv, n)
x.br
xdouble
x.B *A, *Ainv;
x.br
xint
x.B n;
x.SH DESCRIPTION
xThis routine
x.I matrixinv(),
xinverts a matrix
xusing the LU decomposition routines
x.I lu_decompose()
xand
x.I lu_solve()
x( see lu(3) ).
xIt returns 1 if the inversion was successful
xand 0 if the matrix is singular.
x.LP
x.B A
xis the (n x n) matrix to be inverted,
x.B Ainv
xis the (n x n) inverted matrix,
xand
x.B n
xis the order of the matrices.
x.SH AUTHOR
xKen Turkowski,
xCADLINC, Inc.
x.SH SEE ALSO
xlu(3), gausselim(3)
x.SH BUGS
xThere should be a calculation of the condition number
xto indicate how well- or ill- conditioned the matrix is.
!Funky!Stuff!
echo x - matrix_solve/matrixinv.c
cat >matrix_solve/matrixinv.c <<'!Funky!Stuff!'
/* Matrixinv() inverts the matrix A using LU decomposition. Arguments:
* A - the (n x n) matrix to be inverted
* Ainv - the (n x n) inverted matrix
* n - the order of the matrices A and Ainv
*/
# define MAXN 16
# define Ael(i,j) A[i*n+j]
matrixinv(A, Ainv, n)
double *A;
register double *Ainv;
register int n;
{
register int i, j;
double b[MANX], temp;
/* Decompose matrix into L and U triangular matrices */
if (lu_decompose(A, n) == 0)
return(0); /* Singular */
/* Invert matrix by solving n simultaneous equations n times */
for (i = 0; i < n; i++) {
for(j = 0; j < n; j++)
b[j] = 0.0;
b[i] = 1.0;
lu_solve(Ainv[i], b, n); /* Into a row of Ainv: fix later */
}
/* Transpose matrix */
for (i = 0; i < n; i++) {
for (j = 0; j < i; j++) {
temp = Ainv[i][j];
Ainv[i][j] = Ainv[j][i];
Ainv[j][i] = temp;
}
}
return(1);
}
!Funky!Stuff!