1.
program solver
! solve f(x)=0 by direct iteration
implicit none
real, external :: f
real x

  x = solve (f, -1.5)
  write(*,*) x, f(x)

contains

real function solve (f, x0)
  interface
    real function f(x)
    real, intent(in) :: x
    end function
  end interface
  real, intent(in) :: x0
  real:: x, y
  integer :: i=0, maxiter = 10

  x=x0
  do
    y = f(x) 
    write(*,*) x, y
    if (abs(y) < 0.0001) then ! accuracy achieved
      solve=x
      exit
    end if
    x=y
    i=i+1
    if (i > maxiter) exit
  end do
  solve = x
end function solve

end program solver

real function f(x)
real, intent(in) :: x
   if (x < 1) then
      f = -(1-x)**(1/3.0)
   else
      f = (x-1)**(1/3.0)
   end if
end function f


2.

program lintest
use linear
implicit none
  integer, parameter :: n=4
  real, dimension(n,n) :: a
  real, dimension(n) :: b, c, x
  integer, dimension(n) :: order

  a(1,:) = (/1, 2, 3, 4 /)
  a(2,:) = (/1, 3, 1, 2 /)  
  a(3,:) = (/2, 1, 1, 1 /)
  a(4,:) = (/3, 1, 0, 1 /)

  b = (/20, 11, 6, 4/)

  call lu(a, order, n)
  call solve (a, b, order, n, x)

  write(*,*) x
end program

  5.7500010       7.9999990       19.750002      -19.250002  

diagonal dominance added:

 1.6951844      0.72678351      0.18614596     -0.16475782   

3.

program jacobi
! solve a set of linear equations
! using Jacobi iteration
  integer, parameter :: n=4
  real, parameter :: epsilon = 0.00001
  integer i, j, iter
  real d, s, x(n), y(n), A(n,n), M(n), b(n)

  A(1,:) = (/ 11, 2, 3, 4 /)
  A(2,:) = (/ 1, 13, 1, 2 /)
  A(3,:) = (/ 2, 1, 11, 1 /)
  A(4,:) = (/ 3, 1, 0, 11 /)
  b = (/20, 11, 6, 4/)

  ! matrix M, only diagonal elements needed
  do i=1,n 
     M(i)=A(i,i)
  end do

  ! replace A by A-M
  do i=1,n 
    A(i,i)=0
  end do

  ! initial solution
  x=0.0

  ! iterate until the successive solutions
  ! x and y differ less than required accuracy
  d=1.0
  iter=0
    do while (d > epsilon)

    ! right hand side matrix
    do i=1,n
       s=0.0
       do j=1,n 
          s=s-A(i,j)*x(j)
       end do
      y(i) = (s+b(i))/M(i)
    end do
    ! difference of successive iterates
    d=0.0
    do i=1,n 
       d=d+(x(i)-y(i))**2
    end do
    ! update the solution
    x=y

    write (6,*) x, d

   iter = iter+1
   if (iter > 100) exit

 end do
 end


> a.out
 20.0000000 3.66666675 6.00000000 4.00000000 465.444458
 -21.3333340 -7.66666651 -41.6666679 -59.6666679 8162.44482
 399.000000 64.4444427 116.000000 75.6666718 225054.016
 -759.555542 -218.444443 -932.111084 -1257.44446 4297999.50
 8283.00000 1405.85181 3001.00000 2501.11108 114002240.
 -21799.1484 -5425.07422 -20466.9629 -26250.8516 2.32901786E+09
 177274.438 31592.9375 75280.2188 70826.5156 5.95921633E+10
 -572312.625 -131398.891 -456962.313 -563412.250 1.27398799E+12
 3887353.75 718703.500 1839442.38 1848340.75 3.17013256E+13
 -14349077.0 -3141155.50 -10341746.0 -12380761.0 6.98314572E+14
 86830616.0 16484119.0 44220076.0 46188392.0 1.70298191E+16
 -350382016. -74475824.0 -236333744. -276975968. 3.82574246E+17
 1.96585677E+09 380222560. 1.05221581E+09 1.12562189E+09 9.19935321E+18
 -8.41957990E+09 -1.75643878E+09 -5.43755776E+09 -6.27779277E+09 2.09350321E+20
 4.49367204E+10 8.80424141E+09 2.48733921E+10 2.70151782E+10 4.98559850E+21
 -2.00289370E+11 -4.12801556E+10 -1.25692862E+11 -1.43614394E+11 1.14428928E+23
 1.03409648E+12 2.04403671E+11 5.85473262E+11 6.42148270E+11 2.70724906E+24
 -4.73381994E+12 -9.67955448E+11 -2.91474480E+12 -3.30669307E+12 6.24881565E+25
 2.39069183E+13 4.75398354E+12 1.37422889E+13 1.51694155E+13 1.47185582E+27
 -1.11412492E+14 -2.26626802E+13 -6.77372337E+13 -7.64747364E+13 3.41005791E+28
 5.54436020E+14 1.10699737E+14 3.21962392E+14 3.56900172E+14 8.00819362E+29
 -2.61488722E+15 -5.30066241E+14 -1.57647191E+15 -1.77400779E+15 1.86000128E+31
 1.28855794E+16 2.57979170E+15 7.53384830E+15 8.37472774E+15 4.35930449E+32
 -6.12600377E+16 -1.23896277E+16 -3.67256769E+16 -4.12365308E+16 1.01418380E+34
 2.99902414E+17 6.01529239E+16 1.76146230E+17 1.96169745E+17 2.37376944E+35
 -1.43342355E+18 -2.89462723E+17 -8.56127513E+17 -9.59860182E+17 5.52864386E+36
 6.98674878E+18 1.40309050E+18 4.11616994E+18 4.58973344E+18 1.29285768E+38
A floating overflow exception was detected.
 Error occurs at or near line 45 of MAIN__    

> a.out
 1.81818187 0.846153855 0.545454562 0.363636374 4.45151377
 1.38334394 0.608391583 0.104895093 -0.209154502 0.767796993
 1.75501359 0.763851523 0.257643193 -6.89475834E-02 0.205296084
 1.63410521 0.701941431 0.163188085 -0.184444755 4.07130606E-02
 1.71312106 0.736276627 0.201299354 -0.145841554 1.03650866E-02
 1.68244684 0.721327901 0.180302069 -0.170512706 2.21392396E-03
 1.69986260 0.729098201 0.189481005 -0.160788044 5.42508264E-04
 1.69241023 0.725556374 0.184724063 -0.166244179 1.20480239E-04
 1.69633567 0.727334917 0.186897025 -0.163889736 2.88374831E-05
 1.69456351 0.726503611 0.185807586 -0.165121987 6.53695997E-06



4.
mopdify the function simplify to handle also negative values

function simplify (q)
  ! reduce the rational number q to its simplest form
  ! by applying the Euclidean algorithm
  ! the value of the function is the reduced rational number
  type (rational) simplify
  type (rational), intent(in):: q
  integer  a, b, sig
  a=q%nominator
  b=q%denominator
  sig = 1
  if (a*b < 0) then
     sig = -1
     a = abs(a); b=abs(b)
  end if
  do
    if (a < b) then
       c=a; a=b; b=c; 
    end if
    do while (a >= b) 
      a = a-b 
    end do
    if (a == 0) exit
  end do
  simplify%nominator = sig * q%nominator/b
  simplify%denominator = q%denominator/b
end function


program rattest
  use ratpack
  implicit none
  type (rational) :: a, b, c, d, e, f
  a=rat(1,2)
  b=rat(5,3)
  c=rat(3,4)
  d = a+b
  e = 1+c
  f=  b+1
  write(*,*) 'a=', nomin(a), '/', denom(a)
  write(*,*) 'd=', nomin(d), '/', denom(d)
  write(*,*) 'e=', nomin(e), '/', denom(e)
  write(*,*) 'f=', nomin(f), '/', denom(f)

end program

 a=           1 /           2
 d=          13 /           6
 e=           7 /           4
 f=           8 /           3