c c Numerical Analysis: c The Mathematics of Scientific Computing c D.R. Kincaid & E.W. Cheney c Brooks/Cole Publ., 1990 c c Section 4.5 c c Example of Gaussian elimination followed by iterative improvement c c c file: ex2s45.f c parameter (n=4,ia=4) dimension a(ia,n),b(n),r(n),l(n),s(n),x(n) dimension aa(ia,n),e(n) data (a(1,j),j=1,n) /420.,210.,140.,105./ data (a(2,j),j=1,n) /210.,140.,105.,84./ data (a(3,j),j=1,n) /140.,105.,84.,70./ data (a(4,j),j=1,n) /105.,84.,70.,60./ data (b(i),i=1,n) /875.,539.,399.,319./ c print * print *,' Gaussian elimination followed by iterative improvement' print *,' Section 3.5, Kincaid-Cheney' print * c c save a c do 2 i=1,n do 2 j=1,n aa(i,j) = a(i,j) 2 continue c call gauss(a,ia,n,l,s) call solve(a,ia,n,l,b,x) print *,' Initial solution is:' print 5,x print * c do 4 k=1,5 call residual(aa,ia,n,b,x,r) print *,' At iteration',k print *,' Residual vector is:' print 5,r call solve(a,ia,n,l,r,e) print *,' Error vector is:' print 5,e c do 3 i=1,n x(i) = x(i) + e(i) 3 continue c print *,' Solution is:' print 5,x print * 4 continue c 5 format (3x,4(e13.6,2x)) stop end c subroutine residual(a,ia,n,b,x,r) c c computes residual vector r = b-Ax c dimension a(ia,n),b(n),r(n),x(n) double precision sum do 3 i=1,n sum = 0.0 do 2 j=1,n sum = sum + a(i,j)*x(j) 2 continue r(i) = b(i) - sum 3 continue c return end c subroutine gauss(a,ia,n,l,s) c c gaussian elimination with scaled row pivoting c dimension a(ia,n),l(n),s(n) real max do 2 i=1,n l(i) = i s(i) = 0.0 do 2 j=1,n s(i) = max(s(i),abs(a(i,j))) 2 continue do 4 k=1,n-1 rmax = 0.0 do 3 i=k,n r = abs(a(l(i),k))/s(l(i)) if (r .gt. rmax) then j = i rmax = r endif 3 continue c lk = l(j) l(j) = l(k) l(k) = lk do 4 i = k+1,n xmult = a(l(i),k)/a(lk,k) a(l(i),k) = xmult do 4 j = k+1,n a(l(i),j) = a(l(i),j) - xmult*a(lk,j) 4 continue c return end c subroutine solve(a,ia,n,l,b,x) c c back substitution LUx = b c dimension a(ia,n),l(n),b(n),x(n),z(4) c do 3 i=1,n sum = b(l(i)) do 2 j=1,i-1 sum = sum - a(l(i),j)*z(j) 2 continue z(i) = sum 3 continue do 5 i=n,1,-1 sum = z(i) do 4 j=i+1,n sum = sum - a(l(i),j)*x(j) 4 continue x(i) = sum/a(l(i),i) 5 continue c return end