-
Notifications
You must be signed in to change notification settings - Fork 14
/
Copy pathoptional_err.f90
55 lines (55 loc) · 1.68 KB
/
optional_err.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
module reg_mod
implicit none
contains
subroutine linear_regression(x,y,slope,intcep,ierr)
! regress y(:) on x(:), where size(x) == size(y)
real , intent(in) :: x(:),y(:)
real , intent(out) :: slope,intcep
integer, intent(out), optional :: ierr
integer :: ierr_,n
real :: xmean,ymean
if (present(ierr)) ierr = 0
n = size(x)
ierr_ = findloc([n>1, size(y)==n],.false.,dim=1)
if (ierr_ /= 0) then
if (present(ierr)) then
! if error flag passed, set it and return
ierr = ierr_
return
else
! otherwise stop
print*,"in linear_regression, ierr_ =",ierr_
error stop
end if
end if
xmean = sum(x)/n
ymean = sum(y)/n
slope = sum((x-xmean)*(y-ymean))/sum((x-xmean)**2)
! should check that denominator above is nonzero
intcep = ymean - slope*xmean
end subroutine linear_regression
end module reg_mod
!
program test_reg
use reg_mod, only: linear_regression
implicit none
integer, parameter :: n = 10**3
real :: x(n),y(n),noise(n),slope,intcep
integer :: ierr
call random_seed()
call random_number(x)
call random_number(noise)
y = 5*x + 3 + (noise - 0.5)
call linear_regression(x,y,slope,intcep)
print "('slope, intercept =',2(1x,f0.4))",slope,intcep
! pass error flag ierr
call linear_regression(x(2:),y,slope,intcep,ierr)
if (ierr /= 0) print*,"after linear_regression, ierr =",ierr
! no error flag passed -- stop in procedure upon error
call linear_regression(x(2:),y,slope,intcep)
end program test_reg
! sample output:
! slope, intercept = 4.9824 3.0011
! after linear_regression, ierr = 2
! in linear_regression, ierr_ = 2
! ERROR STOP