|
| 1 | +! Matrix pseudo-inversion example: function, subroutine, and operator interfaces |
| 2 | +program example_pseudoinverse |
| 3 | + use stdlib_linalg, only: pinv, pseudoinvert, operator(.pinv.), linalg_state_type |
| 4 | + implicit none(type,external) |
| 5 | + |
| 6 | + real :: A(15,5), Am1(5,15) |
| 7 | + type(linalg_state_type) :: state |
| 8 | + integer :: i, j |
| 9 | + real, parameter :: tol = sqrt(epsilon(0.0)) |
| 10 | + |
| 11 | + ! Generate random matrix A (15x15) |
| 12 | + call random_number(A) |
| 13 | + |
| 14 | + ! Pseudo-inverse: Function interfcae |
| 15 | + Am1 = pinv(A, err=state) |
| 16 | + print *, 'Max error (function) : ', maxval(abs(A-matmul(A, matmul(Am1,A)))) |
| 17 | + |
| 18 | + ! User threshold |
| 19 | + Am1 = pinv(A, rtol=0.001, err=state) |
| 20 | + print *, 'Max error (rtol=0.001): ', maxval(abs(A-matmul(A, matmul(Am1,A)))) |
| 21 | + |
| 22 | + ! Pseudo-inverse: Subroutine interface |
| 23 | + call pseudoinvert(A, Am1, err=state) |
| 24 | + |
| 25 | + print *, 'Max error (subroutine): ', maxval(abs(A-matmul(A, matmul(Am1,A)))) |
| 26 | + |
| 27 | + ! Operator interface |
| 28 | + Am1 = .pinv.A |
| 29 | + |
| 30 | + print *, 'Max error (operator) : ', maxval(abs(A-matmul(A, matmul(Am1,A)))) |
| 31 | + |
| 32 | +end program example_pseudoinverse |
0 commit comments