[Fortran] 纯文本查看 复制代码
include 'link_fnl_shared.h'
use lin_sol_gen_int
use rand_gen_int
use error_option_packet
implicit none
! This is Example 1 for LIN_SOL_GEN.
integer, parameter :: n=32
real(kind(1e0)), parameter :: one=1e0
real(kind(1e0)) err
real(kind(1e0)) A(n,n), b(n,n), x(n,n), res(n,n), y(n**2)
! Generate a random matrix.
call rand_gen(y)
A = reshape(y,(/n,n/))
! Generate random right-hand sides.
call rand_gen(y)
b = reshape(y,(/n,n/))
! Compute the solution matrix of Ax=b.
call lin_sol_gen(A, b, x)
! Check the results for small residuals.
res = b - matmul(A,x)
err = maxval(abs(res))/sum(abs(A)+abs(b))
if (err <= sqrt(epsilon(one))) then
write (*,*) 'Example 1 for LIN_SOL_GEN is correct.'
end if
end