subroutine aptmprd (nrows, smat1, smat2, tol, smat, nerr) ccbeg. cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c SUBROUTINE APTMPRD c c call aptmprd (nrows, smat1, smat2, tol, smat, nerr) c c Version: aptmprd Updated 1989 December 29 11:40. c aptmprd Originated 1989 November 2 14:10. c c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123. c c c Purpose: To find the matrix product smat of the two square matices c smat1 and smat2. Each must be sized nrows by nrows. c Considered as operators, smat is equivalent to the sequential c application of smat2, then smat1. Components of smat within c tol of 0, 1, or -1 will be adjusted to those values. c If nrows is not positive, nerr = 1 will be returned. c c Input: smat1, smat2. c c Output: smat, nerr. c c Glossary: c c nerr Output Indicates an input error, if not 0. c 1 if nrows is non-positive. c c nrows Input Number of rows = number of columns in smat1, smat2, c and smat. c c smat Output A square matrix, with size smat(nrows,nrows). c matrix product of smat1, smat2. c c smat1 Input A square matrix, with size smat1(nrows,nrows). c c smat2 Input A square matrix, with size smat2(nrows,nrows). c c tol Input Numerical tolerance limit. Used to test and adjust c matrix elements. If the row and column vectors of c smat1 and smat2 are unit vectors, then c On Cray computers, recommend 1.e-5 to 1.e-11. c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ccend. c.... Dimensioned arguments. c---- Matrix product of smat1, smat2. dimension smat (nrows,nrows) c---- A square matrix. dimension smat1 (nrows,nrows) c---- A square matrix. dimension smat2 (nrows,nrows) c.... Local variables. c---- Index of matrix row. common /laptmprd/ i c---- Index of matrix column. common /laptmprd/ j c---- One index of matrix element. Comment out to allow use in index register. c____ common /laptmprd/ n cbugc***DEBUG begins. cbug 9901 format (/ 'aptmprd smat1=',3(/ 1p3e22.14)) cbug 9902 format (/ ' smat2=',3(/ 1p3e22.14)) cbug write ( 3, 9901) ((smat1(i,j), j = 1, 3), i = 1, 3) cbug write ( 3, 9902) ((smat2(i,j), j = 1, 3), i = 1, 3) cbugc***DEBUG ends. c.... Test for input errors. nerr = 0 if (nrows .lt. 1) then nerr = 1 cbugc***DEBUG begins. cbug write ( 3, '(/ "aptmprd fatal. bad nrows=",i3)') nrows cbugc***DEBUG ends. go to 210 endif c.... Find the matrix product of smat1 and smat2. The ij'th element of smat c.... is the dot product of the i'th row vector of smat1 and the j'th column c.... vector of smat2. Adjust any components within tol of 0, 1, or -1. c---- Loop over matrix rows. do 130 i = 1, nrows c---- Loop over matrix columns. do 120 j = 1, nrows smat(i,j) = 0.0 c---- Loop over terms of dot product. do 110 n = 1, nrows smat(i,j) = smat(i,j) + smat1(i,n) * smat2(n,j) c---- End of loop over terms of dot product. 110 continue c---- End of loop over matrix columns. 120 continue c---- End of loop over matrix rows. 130 continue c.... See if elements should be tested for adjustment. c---- Test elements for adjustment. if (tol .gt. 0.0) then c---- Loop over matrix rows. do 150 i = 1, nrows c---- Loop over matrix columns. do 140 j = 1, nrows if (abs (smat(i,j)) .lt. tol) then smat(i,j) = 0.0 elseif ((abs (abs (smat(i,j)) - 1.0)) .lt. tol) then smat(i,j) = sign (1.0, smat(i,j)) endif c---- End of loop over matrix columns. 140 continue c---- End of loop over matrix rows. 150 continue c---- Tested tol. endif cbugc***DEBUG begins. cbug 9911 format (/ ' smat=',3(/ 1p3e22.14)) cbug write ( 3, 9911) ((smat(i,j), j = 1, 3), i = 1, 3) cbugc***DEBUG ends. 210 return c.... End of subroutine aptmprd. (+1 line.) end UCRL-WEB-209832