|
| 1 | +! Demonstrate expert subroutine interface with pre-allocated arrays |
| 2 | +program example_lstsq2 |
| 3 | + use stdlib_linalg_constants, only: dp,ilp |
| 4 | + use stdlib_linalg, only: solve_lstsq, lstsq_space, linalg_state_type |
| 5 | + implicit none |
| 6 | + |
| 7 | + integer, allocatable :: x(:),y(:) |
| 8 | + real(dp), allocatable :: A(:,:),b(:),coef(:),real_space(:),singvals(:) |
| 9 | + integer(ilp), allocatable :: int_space(:) |
| 10 | + integer(ilp) :: lrwork,liwork,arank |
| 11 | + |
| 12 | + ! Data set |
| 13 | + x = [1, 2, 2] |
| 14 | + y = [5, 13, 25] |
| 15 | + |
| 16 | + ! Fit three points using a parabola, least squares method |
| 17 | + ! A = [1 x x**2] |
| 18 | + A = reshape([[1,1,1],x,x**2],[3,3]) |
| 19 | + b = y |
| 20 | + |
| 21 | + ! Get storage sizes for the arrays and pre-allocate data |
| 22 | + call lstsq_space(A,b,lrwork,liwork) |
| 23 | + allocate(coef(size(x)),real_space(lrwork),int_space(liwork),singvals(minval(shape(A)))) |
| 24 | + |
| 25 | + ! Solve coefficients of y = coef(1) + x*coef(2) + x^2*coef(3) |
| 26 | + ! with no internal allocations |
| 27 | + call solve_lstsq(A,b,x=coef, & |
| 28 | + real_storage=real_space, & |
| 29 | + int_storage=int_space, & |
| 30 | + singvals=singvals, & |
| 31 | + overwrite_a=.true., & |
| 32 | + rank=arank) |
| 33 | + |
| 34 | + print *, 'parabola: ',coef |
| 35 | + ! parabola: -0.42857142857141695 1.1428571428571503 4.2857142857142811 |
| 36 | + print *, 'rank: ',arank |
| 37 | + ! rank: 2 |
| 38 | + |
| 39 | + |
| 40 | +end program example_lstsq2 |
0 commit comments