program linsolve real a(6,6),at(6,6),b(6),bt(6),bc(6) integer ipiv(6) c c Driver program to test linear equation solvers c Sets up a 6x6 system of linear equations c This is more specific than might be desired. A better form c would declare these arrays to be allocatable. c c c a - the array (2-dimensional)of linear equation coefficients c the first index gives equation number c the second index gives the unknown number c c at - array to keep a copy of a c c b - array (1-dimensional) containing the right hand sides c of the linear equations c c bt - array to keep a copy of bt, in the end it contains the solutions c to the equations. c c bc - an array that should match the original right hand side array c generated by multiplying the solutions times the appropriate c coefficients c c ipiv- special "pivoting" array used to provide an improved solution c in the Gauss Elimination c c c Set up the A matrix and Right Hand Side Generation c c We start by loading a desired solution in the bt array c and creating a matrix "a" with random element values c do 20 i=1,6 bt(i)=float(i) do 20 j=1,6 c a(i,j)=float(i+j-(i*j)/7*4+1) a(i,j)=rand() 20 continue c c Generate the RHS c do 30 i=1,6 b(i)=0. do 30 j=1,6 30 b(i)=b(i)+a(i,j)*bt(j) write(6,*) 'Right Hand Side of the linear equation Ax=b' write(6,2000)(b(i),i=1,6) write(6,*) ' Matrix A' write(6,2000)((a(i,j),j=1,6),i=1,6) 2000 format(1x,6f10.4) c c LINPACK or LAPACK ROUTINES overwrite the A matrix and b array c If you need them again make a copy c do 40 i=1,6 bt(i)=b(i) do 40 j=1,6 40 at(i,j)=a(i,j) c c Insert calls to LINPACK or LAPACK routines here c call sgefa(at,6,6,ipiv,info) c c Check status by writing the info variable from the LU routine c write(6,*) ' INFO = ',info call sgesl(at,6,6,ipiv,bt,0) c c Write out the solution (this assumes you used bt in the last LINPACK c call. c write(6,*) ' Solution of the Linear System:' write(6,2000) (bt(i),i=1,6) c c Check Solution by Multipling Ax c do 50 i=1,6 bc(i)=0. do 50 j=1,6 bc(i)=bc(i)+a(i,j)*bt(j) 50 continue write(6,*)' Ax = ' write(6,2000) (bc(i),i=1,6) stop end