Example manual page

Note that this page is composed of extracts from several different sections of the actual manual. Data fitting, linear systems and general matrix/vector routines are actually described in different sections.

General

(General information, not included in this example.)

Routines

Least squares polynomial fit to 1D data by solving overdetermined system of equations

Polynomial is of form c1+c2.x+c3.x^2+...

      subroutine bm_polyfitg(x,w,npts,n,eps,yc,fail,res,a)

Data in

      real x(npts) : x ordinate array
      real w(npts) : weighting array
      integer npts : no of data pts (x,y pairs)
      integer n : order+1 of polynomial to fit i.e. n=number of coefficients
      real eps : numerical bandwidth used in matrix inversion

Data in/out

      real yc(npts) : y(x) data array. Coefficients returned in first n elements

Data out

      logical fail : true if system of equations was singular
      real res : residue from least squares fit

Working space

      real a(npts,n) : overdetermined system of equations

Solution of matrix equation Ax=y

A:nr by nc, x:nc, y:nr; nr >= nc (i.e. can be overdetermined) uses Givens rotations to make upper triangular, then back substitution

      subroutine bm_matsolg(a,nr,nc,nrm,eps,yx,sing,res,det)

Data in

      real a(nrm,nc) : matrix A (values corrupted in finding solution!)
      integer nr : no of rows
      integer nc : no of cols (less than or equal to nr)
      integer nrm : dimensioned no of rows in A
      real eps : numerical bandwidth used to declare a number insignificant

Upper triangular matrix solver

N.B. Matrix values corrupted in solution.

      subroutine bm_uptrisol(a,nac,narmax,eps,yx,sing)

Data in

      real a(narmax,nac) : upper triangular matrix
      integer nac : number of columns used
      integer narmax : dimensioned rows in A
      real eps : numerical bandwidth used to declare a number insignificant

Data in/out

      real yx(*) : rhs vector (nar values) in; lhs vector (nac values) out
      logical sing : true if matrix too close to singular so cannot solve

Simple display of matrix

Prints to standard output

      subroutine bm_matshow(a,nr,nc,nrm)

Data in

      real a(nrm,*) : matrix A
      integer nr,nc : rows and columns of A
      integer nrm : dimensioned rows of A

Simple display of vector components

      subroutine bm_vecshow(n,a)

Example of use

      parameter (nmax=25,eps=1.0e-6)
      dimension a(nmax,nmax),b(nmax,nmax)
      dimension x(nmax),y(nmax)
      dimension wrk1(nmax),wrk2(nmax),wrk3(nmax),wrk4(nmax)
      logical sing,stuck,hitmax
      external yeqAx,yeqATx ! sparse multiplication routines

c inverse (rotations)
      read *,n ! rows and columns
      read *,((a(i,j),j=1,n),i=1,n)   ! elements
      print *,'Inverse of '
      call bm_matshow(a,n,n,nmax)
      print *,'='
      call bm_matinvg(a,n,nmax,eps,sing,det,b,nmax)
      if (sing) print *,'Matrix singular!'
      call bm_matshow(b,n,n,nmax)
      print *,'Original matrix had determinant=',det

c solution of `sparse' matrix equation (conjugate gradients)
      read *,n   ! rows and columns
      read *,((a(i,j),j=1,n),i=1,n)   ! elements
      read *,(y(i),i=1,n)   ! rhs vector
      read *,(x(i),i=1,n)   ! guess at solution vector
      read *,tol,maxits   ! accuracy,max iterations
      print *,'Solution of '
      call bm_matshow(a,n,n,nmax)
      print *,'*x='
      call bm_vecshow(n,y)
      print *,':'
      call bm_sparsematsolcg(a,n,y,tol,maxits,eps,yeqAx,yeqATx,
     &   x,sing,stuck,hitmax,res,nits,wrk1,wrk2,wrk3,wrk4)
      if (sing) print *,'Matrix singular!'
      if (stuck) print *,'Finished but inaccurate!'
      if (hitmax) print *,'Reached max iterations!'
      print *,'Solution='
      call bm_vecshow(n,x)
      print *,'with residue=',res,' after ',nits,' iterations'
      end

c `Sparse' matrix multiplication routines, setup for non-sparse matrix
      subroutine yeqAx(n,a,x,y)
      parameter (nmax=25)
      dimension a(nmax,*),x(*),y(*)
      do i=1,n
         cum=0.0
         do j=1,n
            cum=cum+a(i,j)*x(j)
         enddo
         y(i)=cum
      enddo
      return
      end

      subroutine yeqATx(n,a,x,y)
      parameter (nmax=25)
      dimension a(nmax,*),x(*),y(*)
      do i=1,n
         cum=0.0
         do j=1,n
            cum=cum+a(j,i)*x(j)
         enddo
         y(i)=cum
      enddo
      return
      end