I suggest taking a look at the FEAPpv finite element code written by Robert
Taylor to accompany this book:
http://www.amazon.com/The-Finite-Element-Method-Fundamentals/dp/1856176339
You can download the fortran source code from this page:
http://www.ce.berkeley.edu/projects/feap/feappv/
The code contains two subroutines, datri
and dasol
that perform the factorization and solution, respectively, of a matrix stored in skyline format; you first call datri
to the factor the matrix and then call dasol
to perform the forward and backward substitutions.
There are sufficiently detailed comments in both subroutines to be able to use them. Here is a simple example with five equations and a symmetric matrix (entries below the diagonal are not included):
program test
C
C Solve A*U = b
C [2 -2 0 0 -1]
C [0 3 -2 0 0]
C A = [0 0 5 -3 0]
C [0 0 0 10 4]
C [0 0 0 0 10]
C
C The solution, U, is: 636.,619.,292.,74.,34.
C
parameter (neqs = 5)
real*8 au(7), ad(neqs), b(neqs), energy
integer jp(neqs)
data au/-2.,-2.,-3.,-1.,0.,0.,4./
data ad/2.,3.,5.,10.,10./
data b/0.,1.,0.,0.,0./
data jp/0,1,2,3,7/
call datri(0, au, ad, jp, neqs, neqs)
call dasol(0, au, ad, b, jp, neqs, neqs, energy)
write(6,*) 'solution', (b(i), i=1,neqs)
end
edit
link under the tags. $\endgroup$