pinvSVD_rel Function

private pure function pinvSVD_rel(A, method, tol) result(Apinv)

Calculates the pseudoinverse of a matrix A using the SVD.

Arguments

Type IntentOptional Attributes Name
real(kind=rk), intent(in), dimension(:, :), contiguous :: A
character(len=*), intent(in), optional :: method
real(kind=rk), intent(in), optional :: tol

Return Value real(kind=rk), dimension(size(A,2), size(A,1))


Calls

proc~~pinvsvd_rel~~CallsGraph proc~pinvsvd_rel forinv::pinvSVD_rel svd svd proc~pinvsvd_rel->svd

Called by

proc~~pinvsvd_rel~~CalledByGraph proc~pinvsvd_rel forinv::pinvSVD_rel proc~pinv_rel forinv::pinv_rel proc~pinv_rel->proc~pinvsvd_rel interface~inv forinv::inv interface~inv->proc~pinv_rel program~test1 test1 program~test1->interface~inv program~test2 test2 program~test2->interface~inv program~test3 test3 program~test3->interface~inv

Source Code

   pure function pinvSVD_rel(A, method, tol) result(Apinv)
#elif defined (IMPURE)
   impure function pinvSVD_rel(A, method, tol) result(Apinv)
#endif

      ! Inputs:
      real(rk), dimension(:, :), contiguous, intent(in)  :: A     ! Input matrix A
      real(rk),     intent(in), optional                 :: tol
      character(*), intent(in), optional                 :: method

      ! Outputs:
      real(rk), dimension(size(A,2), size(A,1))          :: Apinv ! Pseudoinverse of A

      ! Local variables
      real(rk), dimension(size(A,1), size(A,1))          :: U    ! Left singular vectors
      real(rk), dimension(size(A,2), size(A,2))          :: VT   ! Right singular vectors
      real(rk), dimension(min(size(A,1), size(A,2)))     :: S    ! Singular values
      integer                                            :: m, n, i, j, irank, rank

      m = size(A, 1)
      n = size(A, 2)

      call svd(A, U,S,VT, method)

      if (.not. present(tol)) then
         rank = min(m,n)
      else
         rank = count(S > tol)
      end if

      Apinv = 0.0_rk

      do irank = 1, rank
         do j = 1, m
            do i = 1, n
               Apinv(i, j) = Apinv(i, j) + VT(irank, i) * U(j, irank) / S(irank)
            end do
         end do
      end do

   end function pinvSVD_rel