pjs269@tijc02.UUCP (Paul Schmidt ) (04/21/86)
This message is the code I wrote for my Numerical
Computing class at East Tennessee State University.
Instead of doing the "old rm" now that the course is over,
I am posting it on the net hoping that I may save someone
some time.
The bulk of the code exists as routines that I put into
the library "libMAT.a". These routines have a header
block telling what parameters are used and a small description
of what it does. These routines have not been fully debugged
but were run with various test data with no problems.
In all the library routines I define a matrix as:
double **matrix;
You reference an element in the matrix by matrix[x][y].
The mat_new function allocates the space for the matrix and
mat_free frees the space used by the matrix.
Some of the algorithms are from the book Numerical Analysis by
Burden, Faires, and Reynolds. This book, or other numerical
computing books can give a complete explanation of each routine.
I have also supplied the following programs that use the library
routines:
add - add two matrices.
mult - multiply two matrices.
gauss - Gaussian elimination routine with backward
substitution to solve a system of equations.
SOR - Successive Over-relaxation algorithm used to
solve a system of equations.
least_sq - Find the least square equation of degree n for
a set of data points.
w_least - Use weighted least sqaure algorithm to find a
polynomial of degree n for a set of data points.
If anyone has any problems with this code or has any comments
please send me mail.
Good luck,
Paul J. Schmidt
Texas Instruments
P.O. Drawer 1255 M/S 3520
Johnson City, TN 37601
E-Mail:
decvax!mcnc!rti-sel!tijc02!pjs269
#--------------------------- CUT HERE -----------------------------------
#!/bin/sh
# Cut above the preceeding line, or cut here if you must.
# This is a shar archive. Extract with sh, not csh.
# The rest of this file will extract:
# Makefile
# README
# SOR.c
# add.c
# gauss.c
# least_sq.c
# mat_SOR.c
# mat_add.c
# mat_backward.c
# mat_copy.c
# mat_free.c
# mat_gauss.c
# mat_least.c
# mat_mult.c
# mat_new.c
# mat_print.c
# mat_remove.c
# mat_scale.c
# mat_scan.c
# mat_trans.c
# mat_w_least.c
# matrix.h
# mult.c
# norm_2.c
# norm_inf.c
# row_add.c
# row_scale.c
# row_swap.c
# w_least.c
sed 's/^X//' > Makefile << '/*EOF'
X#
X# Makefile for the matrix routines
X#
X
XCFLAGS= -O
XLIB_NAME=MAT
XLIB=lib$(LIB_NAME).a
X
XFILES=$(LIB)(mat_least.o) $(LIB)(mat_w_least.o)\
X $(LIB)(mat_SOR.o) $(LIB)(mat_gauss.o)\
X $(LIB)(mat_add.o) $(LIB)(mat_backward.o)\
X $(LIB)(mat_copy.o) $(LIB)(mat_free.o)\
X $(LIB)(mat_mult.o) $(LIB)(mat_new.o)\
X $(LIB)(mat_print.o) $(LIB)(mat_remove.o) $(LIB)(mat_scale.o)\
X $(LIB)(mat_scan.o) $(LIB)(mat_trans.o)\
X $(LIB)(row_add.o) $(LIB)(row_scale.o) $(LIB)(row_swap.o)\
X $(LIB)(norm_2.o) $(LIB)(norm_inf.o)
X
Xall: lib SOR add gauss least_sq mult w_least
X#
X# make the library.
X#
X
Xlib: $(FILES)
X
XSOR: SOR.c $(LIB) matrix.h
X cc $(CFLAGS) -o SOR SOR.c $(LIB) -lm
X
Xgauss: gauss.c $(LIB) matrix.h
X cc $(CFLAGS) -o gauss gauss.c $(LIB)
X
Xleast_sq: least_sq.c $(LIB) matrix.h
X cc $(CFLAGS) -o least_sq least_sq.c $(LIB)
X
Xw_least: w_least.c $(LIB) matrix.h
X cc $(CFLAGS) -o w_least w_least.c $(LIB)
X
Xadd: add.c $(LIB) matrix.h
X cc $(CFLAGS) -o add add.c $(LIB)
X
Xmult: mult.c $(LIB) matrix.h
X cc $(CFLAGS) -o mult mult.c $(LIB)
/*EOF
ls -l Makefile
sed 's/^X//' > README << '/*EOF'
XThis file describes the code I wrote for my Numerical
XComputing class at East Tennessee State University.
XInstead of doing the "old rm" now that the course is over,
XI am posting it on the net hoping that I may save someone
Xsome time.
X
XThe bulk of the code exists as routines that I put into
Xthe library "libMAT.a". These routines have a header
Xblock telling what parameters are used and a small description
Xof what it does. These routines have not been fully debugged
Xbut were run with various test data with no problems.
X
XIn all the library routines I define a matrix as:
X double **matrix;
XYou reference an element in the matrix by matrix[x][y].
XThe mat_new function allocates the space for the matrix and
Xmat_free frees the space used by the matrix.
X
XSome of the algorithms are from the book Numerical Analysis by
XBurden, Faires, and Reynolds. This book, or other numerical
Xcomputing books can give a complete explanation of each routine.
X
XI have also supplied the following programs that use the library
Xroutines:
X add - add two matrices.
X mult - multiply two matrices.
X gauss - Gaussian elimination routine with backward
X substitution to solve a system of equations.
X SOR - Successive Over-relaxation algorithm used to
X solve a system of equations.
X least_sq - Find the least square equation of degree n for
X a set of data points.
X w_least - Use weighted least sqaure algorithm to find a
X polynomial of degree n for a set of data points.
X
XIf anyone has any problems with this code or has any comments
Xplease send me mail.
X
XGood luck,
X Paul J. Schmidt
X Texas Instruments
X P.O. Drawer 1255 M/S 3520
X Johnson City, TN 37601
X
XE-Mail:
X decvax!mcnc!rti-sel!tijc02!pjs269
X
/*EOF
ls -l README
sed 's/^X//' > SOR.c << '/*EOF'
X/*
X** SOR
X**
X** Program to use Successive Over-Relaxation to solve a system
X** of equations.
X**
X** Author: Paul Schmidt
X** Date: 19 March 1986
X** Class: Numerical & Statistical Methods
X** Number: 222-4267-201
X** Instructor: Jerry Sayers
X** Reference: Numerical Analysis;
X** Burden, Faires, and Reynolds;
X** page 391
X*/
X
X#include <stdio.h>
X#include "matrix.h"
X
Xmain()
X{
X extern double **mat_new();
X
X double **matrix, **x, error, omega;
X int m, n;
X int max_iterations;
X
X fprintf(stderr, "Input size: ");
X scanf("%d %d", &m, &n);
X matrix = mat_new(m, n);
X fprintf(stderr, "Input matrix: \n");
X mat_scan("%g", matrix, m, n);
X printf("System of equations:\n");
X mat_print("%10.6g ", matrix, m, n);
X x = mat_new(m, 1);
X fprintf(stderr, "Initial guess: \n");
X mat_scan("%g", x, m, 1);
X printf("Input initial approximation:\n");
X mat_print("%10.6g ", x, m, 1);
X fprintf(stderr, "What is the tolerance? ");
X scanf("%lg", &error);
X fprintf(stderr, "What is the omega? ");
X scanf("%lg", &omega);
X fprintf(stderr, "What is the maximum number of iterations? ");
X scanf("%d", &max_iterations);
X printf("Tolerance = %lg, omega = %lg, max iterations = %d\n", error, omega, max_iterations);
X if (mat_SOR(matrix, m, x, error, omega, max_iterations) != OK)
X printf("Maximum nuber of iterations exceeded.\n");
X printf("Best approximation:\n");
X mat_print("%10.6g ", x, m, 1);
X mat_free(matrix, m, n);
X exit(0);
X}
/*EOF
ls -l SOR.c
sed 's/^X//' > add.c << '/*EOF'
X/*
X** Program to add two matrices.
X**
X** Author: Paul Schmidt
X** Date: February 1986
X*/
X
X#include <stdio.h>
X
Xmain()
X{
X extern double **mat_new();
X
X double **a, **b, **sum;
X int m, n, m1, n1;
X
X fprintf(stderr,"Input size: ");
X scanf("%d %d", &m, &n);
X a = mat_new(m, n);
X fprintf(stderr,"Input matrix: \n");
X mat_scan("%g",a , m, n);
X fprintf(stderr,"Input size: ");
X scanf("%d %d", &m1, &n1);
X if ((m != m1) || (n != n1))
X {
X mat_free(a, m, n);
X fprintf(stderr,"add error: Matrices not the same size\n");
X exit(-1);
X }
X b = mat_new(m, n);
X fprintf(stderr,"Input matrix: \n");
X mat_scan("%g", b, m, n);
X sum = mat_new(m, n);
X mat_add(sum, a, b, m, n, (double)1.0);
X mat_free(a, m, n);
X mat_free(b, m, n);
X fprintf(stderr,"Matrix sum:\n");
X printf("%d %d\n", m, n);
X mat_print("%10.6g ", sum, m, n);
X mat_free(sum, m, n);
X exit(0);
X}
/*EOF
ls -l add.c
sed 's/^X//' > gauss.c << '/*EOF'
X/*
X** Gauss
X**
X** Program to use gaussian elimination to solve matrices.
X**
X** Author: Paul Schmidt
X** Date: 5 February 1986
X** Class: Numerical & Statistical Methods
X** Number: 222-4267-201
X** Instructor: Jerry Sayers
X** Reference: Numerical Analysis
X** Burden, Faires, and Reynolds
X** pages 268-269
X**
X*/
X
X#include <stdio.h>
X#include "matrix.h"
X
Xmain()
X{
X extern double **mat_new();
X
X double **matrix;
X int m, n;
X
X fprintf(stderr,"Input size: ");
X scanf("%d %d", &m, &n);
X matrix = mat_new(m, n);
X fprintf(stderr,"Input matrix: \n");
X mat_scan("%g", matrix, m, n);
X if (mat_gauss(matrix, m) == NON_UNIQUE)
X printf("Matrix results in non_unique solutions.\n");
X printf("Reduced matrix:\n");
X mat_print("%-10.6g ", matrix, m, n);
X mat_free(matrix, m, n);
X exit(0);
X}
/*EOF
ls -l gauss.c
sed 's/^X//' > least_sq.c << '/*EOF'
X/*
X** least_sq
X**
X** Program to approximate a polynomial by using least squares.
X**
X** Author: Paul Schmidt
X** Date: 5 February 1986
X** Class: Numerical & Statistical Methods
X** Number: 222-4267-201
X** Instructor: Jerry Sayers
X*/
X
X#include <stdio.h>
X#include "matrix.h"
X
Xmain()
X{
X double **x, **y, **coef, **resid;
X int number, degree, i;
X
X extern double **mat_new();
X
X fprintf(stderr, "Program to calculate a polynomial using least squares.\n");
X fprintf(stderr, "What is the degree of the polynomial? ");
X scanf("%d", °ree);
X fprintf(stderr, "How many data points? ");
X scanf("%d", &number);
X x = mat_new(number, 1);
X y = mat_new(number, 1);
X fprintf(stderr, "Input data points (x y):\n");
X for (i=0; i<number; i++)
X {
X scanf("%g", &(x[i][0]));
X scanf("%g", &(y[i][0]));
X }
X
X /* Echo the data points */
X printf("%-10s %-10s\n", " X ", " Y ");
X for (i=0; i<number; i++)
X printf("%- 010.6g %- 010.6g\n", x[i][0], y[i][0]);
X
X /* Allocate space for the coefficients and residuals and call the
X least squares routine. */
X resid = mat_new(number, 1);
X coef = mat_new(degree+1, 1);
X if (mat_least(x, y, number, degree, coef, resid) != OK)
X {
X fprintf(stderr, "Least sqaure calculation failed\n");
X exit(-1);
X }
X
X printf("\nThe approximating polynomial:\n");
X for (i=0; i<=degree; i++)
X {
X if (i==0)
X printf("Y = %g ", coef[i][0]);
X if (i==1)
X printf("%gX ", coef[i][0]);
X if (i>1)
X printf("%gX^%d ", coef[i][0], i);
X if (i != degree)
X printf("+ ");
X }
X printf("\n");
X
X printf("\nThe residuals for the data points are:\n");
X mat_print("%- 010.6g", resid, number, 1);
X
X /* Free the matrices */
X mat_free(x, number, 1);
X mat_free(y, number, 1);
X mat_free(coef, degree+1, 1);
X mat_free(resid, number,1);
X
X /* Exit the program successfully */
X exit(0);
X}
/*EOF
ls -l least_sq.c
sed 's/^X//' > mat_SOR.c << '/*EOF'
X/*
X** mat_SOR
X**
X** Iterative algorithm using Succesive Over-Relaxation to solve
X** a linear system.
X**
X** Input:
X** a - The system of equations.
X** m - The number of unknowns and equations.
X** xo - A nx1 matrix for the initial approximation.
X** error - Error bound used for stopping criterium.
X** omega - Relaxation factor.
X** max_it - Maximum number of iterations.
X**
X** Output:
X** xo - The nxn matrix to solve the system.
X**
X** Author: Paul Schmidt
X** Date: 19 March 1986
X** Class: Numerical & Statistical Methods
X** Number: 222-4267-201
X** Instructor: Jerry Sayers
X** Reference: Numerical Analysis
X** Burden, Faires, and Reynolds
X** page 391
X*/
X
X#include <math.h>
X#include "matrix.h"
X
Xmat_SOR(a, m, xo, error, omega, max_it)
Xdouble **a;
Xint m;
Xdouble **xo;
Xdouble error;
Xdouble omega;
Xint max_it;
X{
X extern double **mat_new();
X extern double norm_2();
X int i, j, k;
X double sum1, sum2;
X double **x, **diff;
X
X /* Allocate space for the local matrices */
X x = mat_new(m, 1);
X diff = mat_new(m, 1);
X
X /* Iterate up to max_it times */
X for (k = 0; k < max_it; k++)
X {
X for (i = 0; i < m; i++)
X {
X /* Calculate the new xi */
X sum1 = 0;
X for (j = 0; j < i; j++)
X sum1 += a[i][j]*x[j][0];
X sum2 = 0;
X for (j = i+1; j < m; j++)
X sum2 += a[i][j]*xo[j][0];
X x[i][0] = (1.0 - omega)*xo[i][0] + omega*(-sum1-sum2+a[i][m])/a[i][i];
X
X }
X
X#ifdef DEBUG
X printf("Next approximation:\n");
X mat_print("%10.6g ", x, m, 1);
X#endif
X
X /* Check for stopping criteria */
X mat_add(diff, x, xo, m, 1, (double)-1.0);
X if (norm_2(diff, m, 1) < error)
X {
X /* Copy the new x into the xo matrix */
X for (i = 0; i < m; i++)
X xo[i][0] = x[i][0];
X /* Values approximated successfully */
X mat_free(x, m, 1);
X mat_free(diff, m, 1);
X return(OK);
X }
X
X /* Copy the new x into the xo matrix */
X for (i = 0; i < m; i++)
X xo[i][0] = x[i][0];
X }
X
X /* Maximum iterations exceeded */
X mat_free(x, m, 1);
X mat_free(diff, m, 1);
X return(MAX_ITERATIONS);
X}
/*EOF
ls -l mat_SOR.c
sed 's/^X//' > mat_add.c << '/*EOF'
X/*
X** mat_add
X**
X** Add two matrices.
X**
X** Input:
X** a, b - the two matrices to add.
X** m, n - the size of the matrices.
X** scale - scale b before add.
X**
X** Output:
X** sum - the sum of two matrices.
X*/
X
Xmat_add(sum, a, b, m, n, scale)
Xdouble **sum, **a, **b;
Xint m, n;
Xdouble scale;
X{
X int i, j;
X
X for (i = 0; i < m; i++)
X for (j = 0; j < n; j++)
X sum[i][j] = a[i][j] + scale * b[i][j];
X}
/*EOF
ls -l mat_add.c
sed 's/^X//' > mat_backward.c << '/*EOF'
X/*
X** mat_backward
X**
X** Do backward substitution on a upper triangular matrix.
X**
X** Input:
X** matrix - matrix to do backward substitiution.
X** m - number of unknowns.
X**
X** Output:
X** matrix - The reduced matrix.
X**
X** Author: Paul Schmidt
X** Date: 5 February 1986
X** Class: Numerical & Statistical Methods
X** Number: 222-4267-201
X** Instructor: Jerry Sayers
X**
X*/
X
X#include "matrix.h"
X
Xmat_backward(matrix, m)
Xdouble **matrix;
Xint m;
X{
X int i, j;
X
X for (i = m-1; i > 0; i--)
X {
X
X /* Scale the row so the first element is 1 */
X /* This row should have been reduced to all zeros except two */
X#ifdef DEBUG
X printf("(matrix[%d]) <-- (%g*matrix[%d])\n", i, 1.0/matrix[i][i], i);
X#endif
X matrix[i][m] /= matrix[i][i];
X matrix[i][i] = (double)1.0;
X#ifdef DEBUG
X mat_print("%10.6g ",matrix, m, m+1);
X#endif
X
X for (j = 0; j < i; j++)
X {
X
X /* subtract the row from every one above */
X#ifdef DEBUG
X printf("(matrix[%d] - %g*matrix[%d]) --> (matrix[%d])\n", j, matrix[j][i], i, j);
X#endif
X matrix[j][m] += -matrix[j][i] * matrix[i][m];
X matrix[j][i] = (double)0.0;
X
X }
X
X#ifdef DEBUG
X mat_print("%10.6g ",matrix, m, m+1);
X#endif
X }
X
X /* Scale the top row */
X#ifdef DEBUG
X printf("(matrix[%d]) <-- (%g*matrix[%d])\n", 0, 1.0/matrix[0][0], 0);
X#endif
X matrix[0][m] /= matrix[0][0];
X matrix[0][0] = (double)1.0;
X#ifdef DEBUG
X mat_print("%10.6g ",matrix, m, m+1);
X#endif
X
X return(OK);
X}
/*EOF
ls -l mat_backward.c
sed 's/^X//' > mat_copy.c << '/*EOF'
X/*
X** mat_copy
X**
X** copy a smaller matrix into a matrix (usually augmented)
X**
X** Input:
X** old - small matrix to copy.
X** m, n - size of Old (mxn).
X** row, col - row and column to start copying to.
X**
X** Output:
X** new - the matrix being copied into.
X*/
X
Xmat_copy(new, old, m, n, row, col)
Xdouble **new, **old;
Xint m, n, row, col;
X{
X int i, j;
X
X for (i = 0; i < m; i++)
X for (j = 0; j < n; j++)
X new[row+i][col+j] = old[i][j];
X}
/*EOF
ls -l mat_copy.c
sed 's/^X//' > mat_free.c << '/*EOF'
X/*
X** mat_free
X**
X** Free the space used by a matrix.
X**
X** Input:
X** matrix - The matrix.
X** m, n - The size of the matrix.
X**
X** Output:
X**
X** Author: Paul Schmidt
X** Date: 5 February 1986
X** Class: Numerical & Statistical Methods
X** Number: 222-4267-201
X** Instructor: Jerry Sayers
X**
X*/
X
Xmat_free(matrix, m, n)
Xdouble **matrix;
Xint m, n;
X{
X int i;
X
X for (i=0; i<m; i++)
X free(matrix[i]);
X free(matrix);
X}
/*EOF
ls -l mat_free.c
sed 's/^X//' > mat_gauss.c << '/*EOF'
X/*
X** mat_gauss
X**
X** gaussian elimination of a matrix using partial pivoting and
X** backward substitiution.
X**
X** Input:
X** matrix - The matrix to reduce.
X** m - number of unknowns.
X**
X** Output:
X** matrix - The reduced matrix.
X**
X** Author: Paul Schmidt
X** Date: 5 February 1986
X** Class: Numerical & Statistical Methods
X** Number: 222-4267-201
X** Instructor: Jerry Sayers
X** Reference: Numerical Analysis
X** Burden, Faires, and Reynolds
X** pages 268-269
X**
X*/
X
X#include <math.h>
X#include "matrix.h"
X
Xmat_gauss(matrix, m)
Xdouble **matrix;
Xint m;
X{
X int i, p, best_p;
X double scale;
X
X#ifdef DEBUG
X printf("Solving:\n");
X mat_print("%10.6g ", matrix, m, m+1);
X#endif
X for (i = 0; i < (m-1); i++)
X {
X
X /* Find the maximum first row element and swap rows */
X best_p = i;
X for (p = (i+1); p < m; p++)
X if (_ABS(matrix[p][i]) > _ABS(matrix[best_p][i]))
X best_p = p;
X if (matrix[best_p][i] == (double)0.0)
X return(NON_UNIQUE);
X if (best_p != i)
X {
X row_swap(matrix, i, best_p);
X#ifdef DEBUG
X printf("(matrix[%d]) <--> (matrix[%d])\n", i, best_p);
X mat_print("%10.6g ", matrix, m, m+1);
X#endif
X }
X
X /* Do Gaussian elmination to eliminate the lower triangle */
X for (p = (i+1); p < m; p++)
X {
X scale = matrix[p][i]/matrix[i][i];
X row_add(matrix, p, i, m+1, -scale);
X /* eliminate roundoff error. */
X matrix[p][i] = 0;
X#ifdef DEBUG
X printf("(matrix[%d] - %g*matrix[%d]) --> (matrix[%d])\n", p, scale, i, p);
X#endif
X }
X#ifdef DEBUG
X mat_print("%10.6g ", matrix, m, m+1);
X#endif
X }
X if (matrix[m-1][m-1] == 0)
X return(NON_UNIQUE);
X
X /* Do backward substitiution */
X mat_backward(matrix, m);
X
X return(OK);
X}
/*EOF
ls -l mat_gauss.c
sed 's/^X//' > mat_least.c << '/*EOF'
X/*
X** mat_least
X**
X** Least squares appoximating polynomial of degree n.
X**
X** Input:
X** x - nx1 matrix of x data points.
X** y - nx1 matrix of y data points.
X** n - number of data points.
X** degree - degree of the approximating polynomial.
X**
X** Output:
X** coef - degreex1 matrix of polynomial coefficients.
X** resid - degreex1 matrix of residuals.
X**
X** Author: Paul Schmidt
X** Date: 19 February 1986
X** Class: Numerical & Statistical Methods
X** Number: 222-4267-201
X** Instructor: Jerry Sayers
X**
X*/
X
X#include <math.h>
X#include "matrix.h"
X
Xmat_least(x, y, n, degree, coef, resid)
Xdouble **x, **y;
Xint n, degree;
Xdouble **coef, **resid;
X{
X double **A; /* matrix of the x^i for i=0..degree, j=1..n */
X double **At; /* Transpose of A */
X double **AtA; /* Transpose of A premultiplies A */
X double **Aty; /* Transpose of A premultiplies y */
X double **Aug; /* Augmneted matrix [AtA|Aty] */
X double **Acoef; /* A premultiplies coef */
X
X double xtemp;
X int i,j;
X int code;
X
X extern double **mat_new();
X
X /* Build matrix A which is:
X [ 1 x0 x0^2 ... x0^n ]
X ...
X [ 1 xi xi^2 ... xi^n ]
X */
X degree++;
X A = mat_new(n, degree);
X for(j=0; j<n; j++)
X {
X xtemp = 1.0;
X for (i=0; i<degree; i++)
X {
X A[j][i] = xtemp;
X xtemp *= x[j][0];
X }
X }
X
X#ifdef DEBUG
X printf("\nA matrix:\n");
X mat_print("%- 010.6g ", A, n, degree);
X#endif
X
X /* Calculate the transpose of A (refered to as At) */
X At = mat_new(degree, n);
X mat_trans(At, A, n, degree);
X
X /* Calculate At premultiply A (AtA) */
X AtA = mat_new(degree, degree);
X mat_mult(AtA, At, A, degree, n, degree);
X
X#ifdef DEBUG
X printf("\nAtA matrix:\n");
X mat_print("%- 010.6g ", AtA, degree, degree);
X#endif
X
X /* Calculate At premultiply y (Aty) */
X Aty = mat_new(degree, 1);
X mat_mult(Aty, At, y, degree, n, 1);
X
X#ifdef DEBUG
X printf("\nAty matrix:\n");
X mat_print("%- 010.6g ", Aty, degree, 1);
X#endif
X
X /* Make an augmented matrix to solve AtA(coef) = Aty */
X Aug = mat_new(degree, degree+1);
X mat_copy(Aug, AtA, degree, degree, 0, 0);
X mat_copy(Aug, Aty, degree, 1, 0, degree);
X
X /* Call the routine to solve the augmented matrix */
X code = mat_gauss(Aug, degree);
X
X if(code == OK)
X {
X /* Remove the solution from the matrix */
X mat_remove(coef, degree, 1, Aug, degree, degree+1, 0, degree);
X
X /* Calculate A premultiply coef (Acoef) */
X Acoef = mat_new(n, 1);
X mat_mult(Acoef, A, coef, n, degree, 1);
X
X /* Calculate the residuals using: y - coef = resid */
X mat_add(resid, y, Acoef, n, 1, (double)-1.0);
X }
X
X /* Free up the temporary matrices used in this function */
X mat_free(A, n, degree);
X mat_free(At, degree, n);
X mat_free(AtA, degree, degree);
X mat_free(Aty, degree, 1);
X mat_free(Aug, degree, degree+1);
X mat_free(Acoef, n, 1);
X
X return(code);
X}
/*EOF
ls -l mat_least.c
sed 's/^X//' > mat_mult.c << '/*EOF'
X/*
X** mat_mult
X**
X** Multiply two matrices.
X**
X** Input:
X** a, b - the two matrices to multiply.
X** m, n, p - the size of the matrices. (mxn and nxp)
X**
X** Output:
X** prod - the product of two matrices.
X*/
X
Xmat_mult(prod, a, b, m, n, p)
Xdouble **prod, **a, **b;
Xint m, n, p;
X{
X int i, j, k;
X
X for (i = 0; i < m; i++)
X for (j = 0; j < p; j++)
X {
X prod[i][j] = 0;
X for (k = 0; k < n; k++)
X prod[i][j] += a[i][k] * b[k][j];
X }
X}
/*EOF
ls -l mat_mult.c
sed 's/^X//' > mat_new.c << '/*EOF'
X/*
X** mat_new
X**
X** Create a new mxn matrix.
X**
X** Input:
X** m, n - The size of the matrix.
X**
X** Output:
X** matrix - The matrix.
X**
X** Author: Paul Schmidt
X** Date: 5 February 1986
X** Class: Numerical & Statistical Methods
X** Number: 222-4267-201
X** Instructor: Jerry Sayers
X**
X*/
X
Xdouble **mat_new(m, n)
Xint m, n;
X{
X double **matrix;
X int i;
X
X matrix = (double **)malloc(sizeof(double *) * m);
X for (i=0; i<m; i++)
X matrix[i] = (double *)calloc(n, sizeof(double));
X return(matrix);
X}
/*EOF
ls -l mat_new.c
sed 's/^X//' > mat_print.c << '/*EOF'
X/*
X** mat_print
X**
X** print the matrix
X**
X** Input:
X** format - format of the numbers.
X** matrix - matrix with size of mxn.
X** m, n - size of matrix.
X**
X*/
X
Xmat_print(format, matrix, m, n)
Xchar *format;
Xdouble **matrix;
Xint m, n;
X{
X int i, j;
X
X for(i = 0; i < m; i++)
X {
X for (j = 0; j < n; j++)
X printf(format, matrix[i][j]);
X printf("\n");
X }
X}
/*EOF
ls -l mat_print.c
sed 's/^X//' > mat_remove.c << '/*EOF'
X/*
X** mat_remove
X**
X** Remove a part of a matrix and put it into new.
X**
X** Input:
X** m_new, n_new - size of new (mxn).
X** old - small matrix to copy.
X** m_old, n_old - size of Old (mxn).
X** row, col - row and column to start copying from.
X**
X** Output:
X** new - the matrix being copied into.
X*/
X
Xmat_remove(new, m_new, n_new, old, m_old, n_old, row, col)
Xdouble **new, **old;
Xint m_new, n_new, m_old, n_old, row, col;
X{
X int i, j;
X
X for (i = 0; i < m_new; i++)
X for (j = 0; j < n_new; j++)
X new[i][j] = old[row+i][col+j];
X}
/*EOF
ls -l mat_remove.c
sed 's/^X//' > mat_scale.c << '/*EOF'
X/*
X** mat_scale
X**
X** Scale a matrix.
X**
X** Input:
X** matrix - matrix to scale.
X** m, n - size of matrix (mxn).
X** scale - scaling factor.
X**
X** Output:
X** new - new matrix that has been scaled.
X**
X*/
X
Xmat_scale(new, matrix, m, n, scale)
Xdouble **new, **matrix;
Xint m, n;
Xdouble scale;
X{
X int i, j;
X
X for (i = 0; i < m; i++)
X for (j = 0; j < n; j++)
X new[i][j] = scale * matrix[i][j];
X}
/*EOF
ls -l mat_scale.c
sed 's/^X//' > mat_scan.c << '/*EOF'
X/*
X** mat_scan
X**
X** Read in an mxn matrix.
X**
X** Input:
X** Format - format to read the values in.
X** matrix - an mxn matrix.
X** m, n - the size of the matrix.
X**
X** Output:
X** matrix - the matrix with all the values set.
X**
X** Author: Paul Schmidt
X** Date: 5 February 1986
X** Class: Numerical & Statistical Methods
X** Number: 222-4267-201
X** Instructor: Jerry Sayers
X**
X*/
X
Xmat_scan(format, matrix, m, n)
Xchar *format;
Xdouble **matrix;
Xint m, n;
X{
X int i, j;
X
X for(i = 0; i < m; i++)
X for (j = 0; j < n; j++)
X scanf(format, &(matrix[i][j]));
X}
/*EOF
ls -l mat_scan.c
sed 's/^X//' > mat_trans.c << '/*EOF'
X/*
X** mat_trans
X**
X** Transpose a matrix.
X**
X** Input:
X** a - the matrix to transpose
X** m, n - the size of the matrix. (mxn)
X**
X** Output:
X** at - the transpose of a.
X*/
X
Xmat_trans(at, a, m, n)
Xdouble **at, **a;
Xint m, n;
X{
X int i, j;
X
X for (i = 0; i < m; i++)
X for (j = 0; j < n; j++)
X at[j][i] = a[i][j];
X}
/*EOF
ls -l mat_trans.c
sed 's/^X//' > mat_w_least.c << '/*EOF'
X/*
X** mat_w_least
X**
X** Weighted least squares appoximating polynomial of degree n.
X**
X** Input:
X** x - nx1 matrix of x data points.
X** y - nx1 matrix of y data points.
X** w - nxn matrix of w weights.
X** n - number of data points.
X** degree - degree of the approximating polynomial.
X**
X** Output:
X** coef - degreex1 matrix of polynomial coefficients.
X** resid - degreex1 matrix of residuals.
X**
X** Author: Paul Schmidt
X** Date: 19 February 1986
X** Class: Numerical & Statistical Methods
X** Number: 222-4267-201
X** Instructor: Jerry Sayers
X**
X*/
X
X#include <math.h>
X#include "matrix.h"
X
Xmat_w_least(x, y, w, n, degree, coef, resid)
Xdouble **x, **y, **w;
Xint n, degree;
Xdouble **coef, **resid;
X{
X double **A; /* matrix of the x^i for i=0..degree, j=1..n */
X double **At; /* Transpose of A */
X double **Atw; /* Transpose of A premultiplies w */
X double **AtwA; /* Atw premultiplies A */
X double **Atwy; /* Atw premultiplies y */
X double **Aug; /* Augmneted matrix [AtwA|Atwy] */
X double **Acoef; /* A premultiplies coef */
X
X double xtemp;
X int i,j;
X int code;
X
X extern double **mat_new();
X
X /* Build matrix A which is:
X [ 1 x0 x0^2 ... x0^n ]
X ...
X [ 1 xi xi^2 ... xi^n ]
X */
X degree++;
X A = mat_new(n, degree);
X for(j=0; j<n; j++)
X {
X xtemp = 1.0;
X for (i=0; i<degree; i++)
X {
X A[j][i] = xtemp;
X xtemp *= x[j][0];
X }
X }
X
X#ifdef DEBUG
X printf("\nA matrix:\n");
X mat_print("%- 010.6g ", A, n, degree);
X#endif
X
X /* Calculate the transpose of A (refered to as At) */
X At = mat_new(degree, n);
X mat_trans(At, A, n, degree);
X
X#ifdef DEBUG
X printf("\nAt matrix:\n");
X mat_print("%- 010.6g ", At, degree, n);
X#endif
X
X /* Calculate At premultiply w (Atw) */
X Atw = mat_new(degree, n);
X mat_mult(Atw, At, w, degree, n, n);
X
X#ifdef DEBUG
X printf("\nAtw matrix:\n");
X mat_print("%- 010.6g ", Atw, degree, n);
X#endif
X
X /* Calculate AtwA */
X AtwA = mat_new(degree, degree);
X mat_mult(AtwA, Atw, A, degree, n, degree);
X
X#ifdef DEBUG
X printf("\nAtwA matrix:\n");
X mat_print("%- 010.6g ", AtwA, degree, degree);
X#endif
X
X /* Calculate At premultiply y (Atwy) */
X Atwy = mat_new(degree, 1);
X mat_mult(Atwy, Atw, y, degree, n, 1);
X
X#ifdef DEBUG
X printf("\nAtwy matrix:\n");
X mat_print("%- 010.6g ", Atwy, degree, 1);
X#endif
X
X /* Make an augmented matrix to solve AtwA(coef) = Atwy */
X Aug = mat_new(degree, degree+1);
X mat_copy(Aug, AtwA, degree, degree, 0, 0);
X mat_copy(Aug, Atwy, degree, 1, 0, degree);
X
X /* Call the routine to solve the augmented matrix */
X code = mat_gauss(Aug, degree);
X
X if(code == OK)
X {
X /* Remove the solution from the matrix */
X mat_remove(coef, degree, 1, Aug, degree, degree+1, 0, degree);
X
X /* Calculate A premultiply coef (Acoef) */
X Acoef = mat_new(n, 1);
X mat_mult(Acoef, A, coef, n, degree, 1);
X
X /* Calculate the residuals using: y - coef = resid */
X mat_add(resid, y, Acoef, n, 1, (double)-1.0);
X }
X
X /* Free up the temporary matrices used in this function */
X mat_free(A, n, degree);
X mat_free(At, degree, n);
X mat_free(Atw, degree, n);
X mat_free(AtwA, degree, degree);
X mat_free(Atwy, degree, 1);
X mat_free(Aug, degree, degree+1);
X mat_free(Acoef, n, 1);
X
X return(code);
X}
/*EOF
ls -l mat_w_least.c
sed 's/^X//' > matrix.h << '/*EOF'
X#define OK 0 /* No error return value */
X#define NON_UNIQUE 1 /* Non-unique answer to a reduced matrix */
X#define MAX_ITERATIONS 2 /* Maximum nuber of iterations excedded */
X#define FALSE 0
X#define TRUE 1
/*EOF
ls -l matrix.h
sed 's/^X//' > mult.c << '/*EOF'
X/*
X** mult
X**
X** Program to multiply two matrices.
X**
X** Author: Paul Schmidt
X** Date: February 1986
X**
X*/
X
X#include <stdio.h>
X
Xmain()
X{
X extern double **mat_new();
X
X double **a, **b, **prod;
X int m, n, m1, n1;
X
X fprintf(stderr,"Input size: ");
X scanf("%d %d", &m, &n);
X a = mat_new(m, n);
X fprintf(stderr,"Input matrix: \n");
X mat_scan("%g",a , m, n);
X fprintf(stderr,"Input size: ");
X scanf("%d %d", &m1, &n1);
X if (n != m1)
X {
X mat_free(a, m, n);
X fprintf(stderr,"mult error: Matrices not the correct size\n");
X exit(-1);
X }
X b = mat_new(m1, n1);
X fprintf(stderr,"Input matrix: \n");
X mat_scan("%g", b, m1, n1);
X prod = mat_new(m, n1);
X mat_mult(prod, a, b, m, n, n1);
X mat_free(a, m, n);
X mat_free(b, m1, n1);
X fprintf(stderr,"Matrix product:\n");
X printf("%d %d\n", m, n1);
X mat_print("%10.6g ", prod, m, n1);
X mat_free(prod, m, n1);
X exit(0);
X}
/*EOF
ls -l mult.c
sed 's/^X//' > norm_2.c << '/*EOF'
X/*
X** norm_2
X**
X** Find the l2 norm of a matrix.
X**
X** Input:
X** A - nxm matrix to calculate norm
X** n, m - size of the matrix.
X**
X** Output:
X** value of the norm.
X**
X** Author: Paul Schmidt
X** Date: 19 March 1986
X** Class: Numerical & Statistical Methods
X** Number: 222-4267-201
X** Instructor: Jerry Sayers
X*/
X
X#include <math.h>
X#include "matrix.h"
X
Xdouble norm_2(a, m, n)
Xdouble **a;
Xint m, n;
X{
X int i, j;
X double sum;
X
X if (n != 1)
X {
X fprintf("norm_2 called with n > 1. Not implemented for more\n");
X return((double)0.0);
X }
X sum = 0.0;
X for (i = 0; i < m; i++)
X {
X sum += a[i][0]*a[i][0];
X }
X return(sqrt(sum));
X}
/*EOF
ls -l norm_2.c
sed 's/^X//' > norm_inf.c << '/*EOF'
X/*
X** norm_inf
X**
X** Find the infinity norm of a matrix.
X**
X** Input:
X** A - nxm matrix to calculate norm
X** n, m - size of the matrix.
X**
X** Output:
X** value of the norm.
X**
X** Author: Paul Schmidt
X** Date: 19 March 1986
X** Class: Numerical & Statistical Methods
X** Number: 222-4267-201
X** Instructor: Jerry Sayers
X*/
X
X#include <math.h>
X#include "matrix.h"
X
Xdouble norm_inf(a, m, n)
Xdouble **a;
Xint m, n;
X{
X int i, j;
X double sum, max;
X
X max = 0.0;
X for (i = 0; i < m; i++)
X {
X sum = 0.0;
X for (j = 0; j < n; j++)
X sum += _ABS(a[i][j]);
X if (sum > max)
X max = sum;
X }
X#ifdef DEBUG
X printf("The infinity norm of:\n");
X mat_print("%10.6g ", a, m, n);
X printf("is: %g\n", max);
X#endif DEBUG
X return(max);
X}
/*EOF
ls -l norm_inf.c
sed 's/^X//' > row_add.c << '/*EOF'
X/*
X** row_add
X**
X** Add a scaled row to a row an replace that row.
X**
X** Input:
X** matrix - matrix to operate on.
X** i - row i to replace.
X** j - row j to scale and add to row i.
X** cols - number of columns.
X** scale - scale of row j.
X**
X** Output:
X**
X** Author: Paul Schmidt
X** Date: 5 February 1986
X** Class: Numerical & Statistical Methods
X** Number: 222-4267-201
X** Instructor: Jerry Sayers
X**
X*/
X
Xrow_add(matrix, i, j, cols, scale)
Xdouble **matrix;
Xint i, j, cols;
Xdouble scale;
X{
X int index;
X
X for (index = 0; index < cols; index++)
X matrix[i][index] += scale * matrix[j][index];
X}
/*EOF
ls -l row_add.c
sed 's/^X//' > row_scale.c << '/*EOF'
X/*
X** row_scale
X**
X** Routine to scale a row in a matrix.
X**
X** Input:
X** matrix - The matrix to operate on.
X** row - row to scale.
X** cols - Number of columns.
X** scale - Scale factor.
X**
X** Output:
X** matrix The altered matrix.
X**
X*/
X
X
Xrow_scale(matrix, row, cols, scale)
Xdouble **matrix;
Xint row, cols;
Xdouble scale;
X{
X int i;
X
X for (i = 0; i < cols; i++)
X matrix[row][i] *= scale;
X}
/*EOF
ls -l row_scale.c
sed 's/^X//' > row_swap.c << '/*EOF'
X/*
X** row_swap
X**
X** Routine to swap two rows in a matrix.
X**
X** Input:
X** matrix The matrix to operate on.
X** i, j Rows to swap.
X**
X** Output:
X** matrix The altered matrix.
X**
X** Author: Paul Schmidt
X** Date: 5 February 1986
X** Class: Numerical & Statistical Methods
X** Number: 222-4267-201
X** Instructor: Jerry Sayers
X**
X*/
X
X
Xrow_swap(matrix, i, j)
Xdouble **matrix;
Xint i, j;
X{
X double *temp;
X
X temp = matrix[i];
X matrix[i] = matrix[j];
X matrix[j] = temp;
X}
/*EOF
ls -l row_swap.c
sed 's/^X//' > w_least.c << '/*EOF'
X/*
X** w_least
X**
X** Program to approximate a polynomial by using weighted least squares.
X**
X** Author: Paul Schmidt
X** Date: 19 March 1986
X** Class: Numerical & Statistical Methods
X** Number: 222-4267-201
X** Instructor: Jerry Sayers
X*/
X
X#include <stdio.h>
X#include "matrix.h"
X
Xmain()
X{
X double **x, **y, **w, **coef, **resid;
X int number, degree, i;
X
X extern double **mat_new();
X
X fprintf(stderr, "Program to calculate a polynomial using weighted least squares.\n");
X fprintf(stderr, "What is the degree of the polynomial? ");
X scanf("%d", °ree);
X fprintf(stderr, "How many data points? ");
X scanf("%d", &number);
X x = mat_new(number, 1);
X y = mat_new(number, 1);
X w = mat_new(number, number);
X fprintf(stderr, "Input data points (x y):\n");
X for (i=0; i<number; i++)
X {
X scanf("%g", &(x[i][0]));
X scanf("%g", &(y[i][0]));
X }
X fprintf(stderr, "Input weights:\n");
X for (i=0; i<number; i++)
X {
X scanf("%g", &(w[i][i]));
X }
X
X /* Echo the data points */
X printf("%-10s %-10s %-10s\n", " X ", " Y ", " WEIGHT ");
X for (i=0; i<number; i++)
X printf("%- 010.6g %- 010.6g %- 010.6g\n", x[i][0], y[i][0], w[i][i]);
X
X /* Allocate space for the coefficients and residuals and call the
X least squares routine. */
X resid = mat_new(number, 1);
X coef = mat_new(degree+1, 1);
X if (mat_w_least(x, y, w, number, degree, coef, resid) != OK)
X {
X fprintf(stderr, "Least sqaure calculation failed\n");
X exit(-1);
X }
X
X printf("\nThe approximating polynomial:\n");
X for (i=0; i<=degree; i++)
X {
X if (i==0)
X printf("Y = %g ", coef[i][0]);
X if (i==1)
X printf("%gX ", coef[i][0]);
X if (i>1)
X printf("%gX^%d ", coef[i][0], i);
X if (i != degree)
X printf("+ ");
X }
X printf("\n");
X
X printf("\nThe residuals for the data points are:\n");
X mat_print("%- 010.6g", resid, number, 1);
X
X /* Free the matrices */
X mat_free(x, number, 1);
X mat_free(y, number, 1);
X mat_free(w, number, number);
X mat_free(coef, degree+1, 1);
X mat_free(resid, number,1);
X
X /* Exit the program successfully */
X exit(0);
X}
/*EOF
ls -l w_least.c
exit