symsolvels()

subroutine, private ompi_kronlininv::symsolvels ( real(dp), dimension(:,:), intent(in)  A, real(dp), dimension(:,:), intent(in)  B, real(dp), dimension(:,:), intent(out)  sol  )

private

Solves a linear system AX = B for symmetric A, real numbers

Parameters

[in]	A	coefficients matrix (linear operator)
[in]	B	matrix
[out]	sol	solution to the linear system (matrix)

Date

5/8/2016 - Initial Version

Definition at line 1245 of file ompi_kronlininv.f08.

    real(dp),intent(in) :: a(:,:),b(:,:)
    real(dp),intent(out) :: sol(:,:)
    integer,parameter :: wmax=3000 ! max size allowed for work
    integer,allocatable :: ipiv(:)
    integer :: n,nrhs,lda,ldb,lwork
    character :: uplo
    real(dp),allocatable :: work(:),a2(:,:)
    integer :: info
    ! External procedures defined in LAPACK
    external dsysv
    uplo = 'L'
    n = size(a,1)
    nrhs = size(b,2)
    allocate(a2(n,n))
    a2 = a
    lda = size(a,1)
    allocate(ipiv(n))
    ldb = size(b,1)
    ! copy B to sol to avoid it to be overwritten
    sol = b
    ! apparently work must be allocated for the query to dsysv to work...
    allocate(work(1)) 
    lwork = -1
    call dsysv(uplo,n,nrhs,a2,lda,ipiv,sol,ldb,work,lwork,info)
    lwork = min(wmax,int(work(1)))
    !! now lwork is assigned, so actually calculate eigvec/val
    deallocate(work)
    allocate(work(lwork))
    call dsysv(uplo,n,nrhs,a2,lda,ipiv,sol,ldb,work,lwork,info)
    if (info/=0) then
       write(*,*) "linear system solver failed...'"
       print*,'info: ',info
       stop
    end if

Here is the caller graph for this function: