subroutine aptsolv (a, s, nrows, w, nrowm, tol, q, b, qm, qr,
& sm, sr, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTSOLV
c
c call aptsolv (a, s, nrows, w, nrowm, tol, q, b, qm, qr,
c & sm, sr, nerr)
c
c Version: aptsolv Updated 1990 December 3 14:20.
c aptsolv Originated 1990 October 18 15:20.
c
c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c Purpose: To solve a set of nrows equations a * q = s for the nrows
c unknowns "q", given the nrows by nrows coefficient matrix "a"
c and the nrows right-side constants "s". The method used is
c Gauss-Jordan elimination with partial pivoting. In addition
c to "q", the inverse matrix b = a**-1, the recalculated unknowns
c qm = b * s, the residuals of the unknowns qr = b * s - q,
c the recalculated right-side constants sm = a * q, and the
c residuals of the right-side constants sr = a * q - s, are
c calculated and returned. Failure to obtain a good solution
c is indicated by a non-zero value of nerr, but results may be
c usable with nerr = 2 or 3, which indicate numerical truncation
c error.
c
c Input: a, s, nrows, w, nrowm, tol.
c
c Output: q, b, qm, qr, sm, sr, nerr.
c
c Glossary:
c
c a Input Matrix of coefficients of equations to be solved.
c Shape in memory is nrowm by nrows(+).
c Data is in first nrows rows and first nrows columns.
c
c b Output The inverse of matrix "a".
c Shape in memory is nrowm by nrows(+).
c Data is in first nrows rows and first nrows columns.
c
c nerr Output Error indicator. 0 if a good solution was obtained.
c 1 if the matrix of coefficients "a" is singular.
c The determinant of the matrix "a" is zero.
c 2 if any residual "qr" exceeds 1.e-6 times the
c maximum q value. Results may be usable.
c 3 if any residual "sr" exceeds 1.e-6 times the
c maximum r value. Results may be usable.
c
c nrowm Input First dimension of array "a". Increment in array "a"
c between successive column values.
c
c nrows Input Number of equations and number of unknowns.
c Size of square matrix "a".
c Size of arrays q, qm, qr, s, sm, sr.
c
c q Output Unknowns, obtained by Gauss-Jordan elimination, with
c partial pivoting. Size nrows.
c
c qm Output Unknowns calculated from b * s. Size nrows.
c
c qr Output The residuals of the unknowns b * s - q.
c Small if the solution is good. Size nrows.
c
c s Input Solution vector. Right-hand side of equations.
c Size nrows.
c
c sm Output Solution vector calculated from a * q. Size nrows.
c
c sr Output The residuals of the right-side constants a * q - s.
c Small if the solution is good. Size nrows.
c
c w Input Working memory space for matrix calculations.
c Shape in memory is nrowm by nrows + 1(+).
c Data is in first nrows rows and in first nrows + 1
c columns.
c
c tol Input Numerical tolerance limit.
c On Cray computers, recommend 1.e-5 to 1.e-11.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.
c.... Dimensioned arguments.
dimension a (nrowm,nrows) ! Coefficient matrix.
dimension b (nrowm,nrows) ! Inverse matrix by elimination.
dimension q (nrowm) ! Unknowns by elimination.
dimension qm (nrowm) ! Unknowns from b * s.
dimension qr (nrowm) ! Residuals qr = qm - q.
dimension s (nrowm) ! Solution vector s = a * q.
dimension sm (nrowm) ! Estimate sm = a * q.
dimension sr (nrowm) ! Residuals sr = sm - s.
dimension w (nrowm,nrows+1) ! Working coefficient matrix.
c.... Local variables.
common /laptsolv/ adiv ! Matrix element at pivot row, column.
common /laptsolv/ amax ! Maximum matrix element in a column.
common /laptsolv/ amaxabs ! Abs value of max column element.
common /laptsolv/ amult ! Initial pivot column element of row.
common /laptsolv/ atemp ! Temporary matrix element.
common /laptsolv/ npivot ! Index of pivot point row and column.
common /laptsolv/ nq ! Index of q = row index.
common /laptsolv/ nr ! Row index.
common /laptsolv/ nrmax ! Index of row containing amax.
common /laptsolv/ qmax ! Maximum absolute value of q.
common /laptsolv/ smax ! Maximum absolute value of s.
cbugc***DEBUG begins.
cbug common /laptsolv/ alabq (12)
cbug data (alabq(n), n = 1, 12) / 'q1', 'q2', 'q3', 'q4', 'q5', 'q6',
cbug & 'q7', 'q8', 'q9', 'q10', 'q11',
cbug & 'q12' /
cbug common /laptsolv/ ns ! Row index.
cbug common /laptsolv/ u (12,12) ! Check matrix b * a. Identity?
cbug 9901 format (/ 'aptsolv solving equation set. tol=',1pe13.5)
cbug 9902 format (/ 'nrows=',i3,'. matrix a and solution vector s:' /)
cbug 9903 format (1x,1p12e13.5)
cbug write ( 3, 9901) tol
cbug write ( 3, 9902) nrows
cbug do 108 nr = 1, nrows ! Loop over rows.
cbug write ( 3, 9903) (a(nr,nc), nc = 1, nrows), s(nr)
cbug 108 continue ! End of loop over rows.
cbugc***DEBUG ends.
c.... Initialize.
nerr = 0
ncols = nrows + 1
c.... Store coefficients "a" and solution vector "s" in working matrix "w".
do 102 nr = 1, nrows ! Loop over rows.
do 100 nc = 1, nrows ! Loop over first nrows columns.
w(nr,nc) = a(nr,nc)
100 continue ! End of loop over columns.
w(nr,ncols) = s(nr)
102 continue ! End of loop over rows.
c.... Form initial inverse matrix "b".
do 106 nr = 1, nrows ! Loop over rows.
do 104 nc = 1, nrows ! Loop over first nrows columns.
b(nr,nc) = 0.0
104 continue ! End of loop over columns.
b(nr,nr) = 1.0
106 continue ! End of loop over rows.
cbugcbugc***DEBUG begins.
cbugcbug 9904 format (/ 'nrows=',i3,'. inverse matrix b:' /)
cbugcbug write ( 3, 9904) nrows
cbugcbug do 110 nr = 1, nrows ! Loop over rows.
cbugcbug write ( 3, 9903) (b(nr,nc), nc = 1, nrows)
cbugcbug 110 continue ! End of loop over rows.
cbugcbugc***DEBUG ends.
c.... Do Gauss-Jordan elimination with partial pivoting.
do 140 npivot = 1, nrows ! Loop over pivot rows.
c.... Find max element in column npivot. Save row index.
amax = w(npivot,npivot)
amaxabs = abs (amax)
nrmax = npivot
if (npivot .lt. nrows) then
do 112 nr = npivot + 1, nrows ! Loop over pivot rows.
if (abs (w(nr,npivot)) .gt. amaxabs) then
amax = w(nr,npivot)
amaxabs = abs (amax)
nrmax = nr
endif
112 continue ! End of loop over pivot rows.
endif
cbugcbugc***DEBUG begins.
cbugcbug 9905 format (/ 10("--------") //
cbugcbug & "largest lower triangular element in column",i3,
cbugcbug & " is",1pe13.5," in row",i3,".")
cbugcbug write ( 3, 9905) npivot, amax, nrmax
cbugcbugc***DEBUG ends.
if (amax .eq. 0.0) then
nerr = 1
cbugcbugc***DEBUG begins.
cbugcbug 9906 format ('aptsolv results: coefficient matrix a is singular.',
cbugcbug & ' no solution.' /
cbugcbug & 2x,5("*****")," failed ",5("*****"))
cbugcbug write ( 3, '(/ 10("--------"))')
cbugcbug write ( 6, 9906)
cbugcbug write ( 3, '(/)')
cbugcbug write ( 3, 9906)
cbugcbugc***DEBUG ends.
go to 210
endif
c.... Swap maximum element to row npivot.
if (nrmax .ne. npivot) then ! Swap w row.
do 114 nc = npivot, ncols ! Loop over pivit columns.
atemp = w(npivot,nc)
w(npivot,nc) = w(nrmax, nc)
w(nrmax, nc) = atemp
114 continue ! End of loop over pivot columns.
do 116 nc = 1, nrows ! Loop over first nrows columns.
atemp = b(npivot,nc)
b(npivot,nc) = b(nrmax, nc)
b(nrmax, nc) = atemp
116 continue ! End of loop over columns.
endif ! Tested nrmax vs. npivot.
cbugcbugc***DEBUG begins.
cbugcbug 9907 format (/ "swapped rows",i3," and",i3," of matrix.")
cbugcbug if (nrmax .ne. npivot) then
cbugcbug write ( 3, 9907) nrmax, npivot
cbugcbug write ( 3, 9902) nrows
cbugcbug do 118 nr = 1, nrows ! Loop over rows.
cbugcbug write ( 3, 9903) (w(nr,nc), nc = 1, ncols)
cbugcbug 118 continue ! End of loop over rows.
cbugcbug write ( 3, 9904) nrows
cbugcbug do 120 nr = 1, nrows ! Loop over rows.
cbugcbug write ( 3, 9903) (b(nr,nc), nc = 1, nrows)
cbugcbug 120 continue ! End of loop over rows.
cbugcbug endif
cbugcbugc***DEBUG ends.
c.... Make element w(npivot,npivot) = 1.0.
adiv = w(npivot,npivot)
if (adiv .ne. 1.0) then ! Divide column by adiv.
w(npivot,npivot) = 1.0
do 122 nc = npivot + 1, ncols ! Loop over pivot columns.
w(npivot,nc) = w(npivot,nc) / adiv
122 continue ! End of loop over pivot columns.
do 124 nc = 1, nrows ! Loop over first nrows columns.
b(npivot,nc) = b(npivot,nc) / adiv
124 continue ! End of loop over columns.
endif ! Tested adiv.
cbugcbugc***DEBUG begins.
cbugcbug 9908 format (/ "divided row",i3," by",1pe13.5,".")
cbugcbug if (abs (adiv - 1.0) .gt. tol) then
cbugcbug write ( 3, 9908) npivot, adiv
cbugcbug write ( 3, 9902) nrows
cbugcbug do 126 nr = 1, nrows ! Loop over rows
cbugcbug write ( 3, 9903) (w(nr,nc), nc = 1, ncols)
cbugcbug 126 continue ! End of loop over rows.
cbugcbug write ( 3, 9904) nrows
cbugcbug do 128 nr = 1, nrows ! Loop over rows.
cbugcbug write ( 3, 9903) (b(nr,nc), nc = 1, nrows)
cbugcbug 128 continue ! End of loop over rows.
cbugcbug endif
cbugcbugc***DEBUG ends.
c.... Make other elements in pivot column zero.
cbugcbugc***DEBUG begins.
cbugcbug 9909 format (" row",i3,", column",i3," reduced from",1pe13.5," to",
cbugcbug & 1pe13.5,", truncated to zero.")
cbugcbug 9910 format ("subtracted",1pe13.5," * row",i3," from row",i3,".")
cbugcbug write ( 3, '(/)')
cbugcbugc***DEBUG ends.
amult = 0.0
do 134 nr = 1, nrows ! Loop over rows.
if (nr .ne. npivot) then ! Zero out row element.
amult = w(nr,npivot)
if (amult .ne. 0.0) then ! Zero out row element.
cbugcbugc***DEBUG begins.
cbugcbug write ( 3, 9910) amult, npivot, nr
cbugcbugc***DEBUG ends.
w(nr,npivot) = 0.0
do 130 nc = npivot + 1, ncols ! Loop over pivot columns.
atemp = abs (w(nr,nc))
w(nr,nc) = w(nr,nc) - amult * w(npivot,nc)
if (abs(w(nr,nc)) .lt. tol * atemp) then
cbugcbugc***DEBUG begins.
cbugcbug if (w(nr,nc) .ne. 0.0) then
cbugcbug write ( 3, 9909) nr, nc, atemp, w(nr,nc)
cbugcbug endif
cbugcbugc***DEBUG ends.
w(nr,nc) = 0.0
endif
130 continue ! End of loop over pivot columns.
do 132 nc = 1, nrows ! Loop over first nrows columns.
atemp = abs (b(nr,nc))
b(nr,nc) = b(nr,nc) - amult * b(npivot,nc)
if (abs(b(nr,nc)) .lt. tol * atemp) then
cbugcbugc***DEBUG begins.
cbugcbug if (b(nr,nc) .ne. 0.0) then
cbugcbug write ( 3, 9909) nr, nc, atemp, b(nr,nc)
cbugcbug endif
cbugcbugc***DEBUG ends.
b(nr,nc) = 0.0
endif
132 continue ! End of loop over columns.
endif ! Tested amult.
endif ! Tested nr vs. npivot.
134 continue ! End of loop over rows.
cbugcbugc***DEBUG begins.
cbugcbug write ( 3, 9902) nrows
cbugcbug do 136 nr = 1, nrows ! Loop over rows.
cbugcbug write ( 3, 9903) (w(nr,nc), nc = 1, ncols)
cbugcbug 136 continue ! End of loop over rows.
cbugcbug write ( 3, 9904) nrows
cbugcbug do 138 nr = 1, nrows ! Loop over rows.
cbugcbug write ( 3, 9903) (b(nr,nc), nc = 1, nrows)
cbugcbug 138 continue ! End of loop over rows.
cbugcbugc***DEBUG ends.
140 continue ! End of loop over pivot rows.
c.... Get the residuals.
qmax = 0.0
smax = 0.0
do 146 nq = 1, nrows ! Loop over rows.
q(nq) = w(nq,ncols)
qmax = amax1 (abs (q(nq)), qmax)
smax = amax1 (abs (s(nq)), smax)
146 continue ! End of loop over rows.
do 152 nr = 1, nrows ! Loop over rows.
qm(nr) = 0.0
sm(nr) = 0.0
do 150 nc = 1, nrows ! Loop over first nrows columns.
qm(nr) = qm(nr) + b(nr,nc) * s(nc)
sm(nr) = sm(nr) + a(nr,nc) * q(nc)
cbugcbugc***DEBUG begins.
cbugcbug u(nr,nc) = 0.0
cbugcbug do 148 nu = 1,nrows ! Loop over rows.
cbugcbug u(nr,nc) = u(nr,nc) + b(nu,nc) * a(nr,nu)
cbugcbug 148 continue ! End of loop over rows.
cbugcbugc***DEBUG ends.
150 continue ! End of loop over columns.
qr(nr) = qm(nr) - q(nr)
if (abs (qr(nr)) .lt. tol * abs (q(nr))) then
qr(nr) = 0.0
endif
if (abs (qr(nr)) .gt. 1.e-6 * qmax) nerr = 2
sr(nr) = sm(nr) - s(nr)
if (abs (sr(nr)) .lt. tol * abs (s(nr))) then
sr(nr) = 0.0
endif
if (abs (sr(nr)) .gt. 1.e-6 * smax) nerr = 3
152 continue ! End of loop over rows.
210 continue
cbugc***DEBUG begins.
cbug 9911 format ('nerr =',i2,'. 1 = singular matrix.',
cbug & ' 2 = large q residuals. 3 = large s residuals.')
cbug 9912 format (/ 'unknowns q:' //
cbug & ' nq',8x,'q (fixed)',7x,'q (floating)' /
cbug & (1x,a2,0pf20.12,1pe20.12))
cbug 9913 format (/ 'unknowns q and residuals:' //
cbug & ' nq',4x,'q',19x,'b * s',15x,'b * s - q' /
cbug & (1x,a2,1p3e20.12))
cbug 9914 format (/ 'solution vector s and residuals:' //
cbug & ' ns',4x,'s',19x,'a * q',15x,'a * q - s' /
cbug & (i3,1p3e20.12))
cbug 9915 format (/ 'matrix product u = b * a:' /)
cbug 9916 format (1x,1p12e13.5)
cbug 9920 format (/ 'aptsolv results. nerr = ',i2)
cbug
cbug write ( 3, '(/ 10("--------"))')
cbug write ( 3, 9920) nerr
cbug
cbug if (nerr .ne. 0) then
cbug write ( 6 9911) nerr
cbug write ( 3, '(/)')
cbug write ( 3, 9911) nerr
cbug if (nerr .eq. 1) go to 220
cbug endif
cbug
cbug write ( 6, 9912) (alabq(nq), q(nq), q(nq), nq = 1, nrows)
cbug write ( 3, 9912) (alabq(nq), q(nq), q(nq), nq = 1, nrows)
cbug write ( 3, 9913) (alabq(nq), q(nq), qm(nq), qr(nq), nq = 1, nrows)
cbug write ( 3, 9914) (ns, s(ns), sm(ns), sr(ns), ns = 1, nrows)
cbug write ( 3, 9915)
cbug
cbug do 160 nr = 1, nrows ! Loop over rows.
cbug do 158 nc = 1, nrows ! Loop over first nrows columns.
cbug if (abs (u(nr,nc)) .lt. tol) then
cbug u(nr,nc) = 0.0
cbug endif
cbug 158 continue ! End of loop over columns.
cbug write ( 3, 9916) (u(nr,nc), nc = 1, nrows)
cbug 160 continue ! End of loop over rows.
cbugc***DEBUG ends.
220 return
c.... End of subroutine aptsolv. (+1 line.)
end
UCRL-WEB-209832