!-*-f90-*-
!
!  API: Array support
!> \page "Comments on vectors and matrices"
!> Please go to api/array.finc for the API documentation.
!> Since array processing is one of the strengths of Fortran, FGSL focuses on
!> leveraging Fortran-style array processing for those GSL routines which
!> require arguments of type <CODE>fgsl_vector*</CODE> or <CODE>fgsl_matrix*</CODE>.
!
! vectors (extended precision real)
!
!> Initialize a GSL vector object. This is invoked via the generic
!> fgsl_vector_init.
!> \param[in] array. The result variable's block is aliased to this
!! contiguous array or a section of it. The actual argument must be
!! a CONTIGUOUS array with the TARGET attribute. It can be of type
!! integer(fgsl_int) or real(fgsl_double).
!> \param[in] stride. If present, the stride between subsequent array
!! elements of the function result. Otherwise, the value one is assumed.
!> \param[inout] status. If present, the exit status.
  function fgsl_vector_init(array, stride, stat)
    real(fgsl_double), target, contiguous, intent(in) :: array(:)
    integer(fgsl_size_t), intent(in), optional :: stride
    integer(fgsl_int), intent(inout), optional :: stat
    type(fgsl_vector) :: fgsl_vector_init
    integer(fgsl_size_t) :: stride_local, array_size, section_size
    integer(fgsl_int) :: stat_local

    try : block

      if (present(stride)) then
         stride_local = stride
      else
         stride_local = 1_fgsl_size_t
      end if
      if (stride_local <= 0) then
         stat_local = fgsl_einval
         exit try
      end if

      array_size = size(array,dim=1,kind=fgsl_size_t)
      section_size = (array_size - 1) / stride_local + 1
      
      fgsl_vector_init%gsl_vector = fgsl_aux_vector_double_init()
      stat_local = fgsl_aux_vector_double_align(c_loc(array), &
           array_size, fgsl_vector_init%gsl_vector, &
           section_size, 0_fgsl_size_t, stride_local)
    end block try
    if ( present(stat) ) stat = stat_local
    if ( .not. present(stat) .and. stat_local /= fgsl_success ) &
         call fgsl_error("aligning failed", __FILE__, __LINE__, stat_local)
  end function fgsl_vector_init
  function fgsl_vector_int_init(array, stride, stat)
    integer(fgsl_int), target, contiguous, intent(in) :: array(:)
    integer(fgsl_size_t), intent(in), optional :: stride
    integer(fgsl_int), intent(inout), optional :: stat
    type(fgsl_vector_int) :: fgsl_vector_int_init
    integer(fgsl_size_t) :: stride_local, array_size, section_size
    integer(fgsl_int) :: stat_local

    try : block

      if (present(stride)) then
         stride_local = stride
      else
         stride_local = 1_fgsl_size_t
      end if
      if (stride_local <= 0) then
         stat_local = fgsl_einval
         exit try
      end if

      array_size = size(array,dim=1,kind=fgsl_size_t)
      section_size = (array_size - 1) / stride_local + 1
      
      fgsl_vector_int_init%gsl_vector_int = fgsl_aux_vector_int_init()
      stat_local = fgsl_aux_vector_int_align(c_loc(array), &
           array_size, fgsl_vector_int_init%gsl_vector_int, &
           section_size, 0_fgsl_size_t, stride_local)
    end block try
    if ( present(stat) ) stat = stat_local
    if ( .not. present(stat) .and. stat_local /= fgsl_success ) &
         call fgsl_error("aligning failed", __FILE__, __LINE__, stat_local)
  end function fgsl_vector_int_init
!> Legacy specific fgsl_vector_init of  for GSL vector initialization 
!> \param type - determine intrinsic type of vector object
!> \return new object of type fgsl_vector
  function fgsl_vector_init_legacy(type)
    real(fgsl_double), intent(in) :: type
    type(fgsl_vector) :: fgsl_vector_init_legacy
    fgsl_vector_init_legacy%gsl_vector = fgsl_aux_vector_double_init()
  end function fgsl_vector_init_legacy
!> Legacy function to wrap a rank 1 Fortran array slice inside a double
!!  precision real GSL vector object. This is invoked via the generic
!! fgsl_vector_align.
!! It is recommended to update codes using this to use the
!! new fgsl_vector_init specific instead
!> \param array - requires the actual argument to have the
!!  TARGET attribute. Otherwise being passed by reference is
!!  not guaranteed by the Fortran standard.
!> \param len - number of elements of the rank 1 array
!> \param fvec - previously initialized GSL vector object
!> \param size - number of elements from array wrapped inside fvec
!> \param offset - index of first element of array to be mapped to fvec
!> \param stride - stride in array for successive elements of fvec
!> \return Status
  function fgsl_vector_align(array, len, fvec, size, offset, stride)
    integer(fgsl_size_t), intent(in) :: len, size, offset, stride
    real(fgsl_double), dimension(len), target, intent(in) :: array
    type(fgsl_vector), intent(inout) :: fvec
    integer(fgsl_int) :: fgsl_vector_align
!
    fgsl_vector_align = fgsl_aux_vector_double_align(c_loc(array), len, &
         fvec%gsl_vector, size, offset, stride)
  end function fgsl_vector_align
!> Function to associate a Fortran pointer with a GSL vector object.
!> \param[in] fvec. double precision real GSL vector
!> The function result is a null pointer if the object is invalid,
!! otherwise it points to the data described by the fvec object  
  function fgsl_vector_to_fptr(fvec)
    type(fgsl_vector), intent(in) :: fvec
    real(fgsl_double), pointer :: fgsl_vector_to_fptr(:)
    real(fgsl_double), pointer, contiguous :: fptr_local(:)
    integer(fgsl_size_t) :: size, stride
    type(c_ptr) :: cp

    if ( fgsl_vector_status(fvec) ) then
       size = fgsl_aux_vector_double_size(fvec%gsl_vector)
       stride = fgsl_aux_vector_double_stride(fvec%gsl_vector)
       if (stride == 0) then
          fgsl_vector_to_fptr => null()
       else
          cp = gsl_vector_ptr(fvec%gsl_vector,0_fgsl_size_t)
          call c_f_pointer(cp, fptr_local, [ size*stride ])
          fgsl_vector_to_fptr => fptr_local(1:size*stride:stride)
       end if
    else
       fgsl_vector_to_fptr => null()
    end if
  end function fgsl_vector_to_fptr
  function fgsl_vector_int_to_fptr(fvec)
    type(fgsl_vector_int), intent(in) :: fvec
    integer(fgsl_int), pointer :: fgsl_vector_int_to_fptr(:)
    integer(fgsl_int), pointer, contiguous :: fptr_local(:)
    integer(fgsl_size_t) :: size, stride
    type(c_ptr) :: cp

    if ( fgsl_vector_int_status(fvec) ) then
       size = fgsl_aux_vector_int_size(fvec%gsl_vector_int)
       stride = fgsl_aux_vector_int_stride(fvec%gsl_vector_int)
       if (stride == 0) then
          fgsl_vector_int_to_fptr => null()
       else
          cp = gsl_vector_int_ptr(fvec%gsl_vector_int,0_fgsl_size_t)
          call c_f_pointer(cp, fptr_local, [ size*stride ])
          fgsl_vector_int_to_fptr => fptr_local(1:size*stride:stride)
       end if
    else
       fgsl_vector_int_to_fptr => null()
    end if
  end function fgsl_vector_int_to_fptr
!> Legacy function to associate a Fortran pointer with the data stored inside
!> a GSL vector object. Codes should be updated to use fgsl_vector_ptr.
!> This is invoked via the generic fgsl_vector_align. Objects of type
!> <CODE>gsl_vector</CODE> which are returned by GSL routines often are
!> persistent subobjects of other GSL objects. A Fortran pointer aligned with
!> a subobject hence will remain up-to-date throughout the lifetime of the
!> object; it may become undefined once the object ceases to exist.
!> \param ptr - rank 1 Fortran pointer
!> \param fvec - double precision real GSL vector
!> \return Status
  function fgsl_vector_pointer_align(ptr, fvec)
    real(fgsl_double), pointer, intent(out) :: ptr(:)
    type(fgsl_vector), intent(in) :: fvec
    integer(fgsl_int) :: fgsl_vector_pointer_align
!
    real(fgsl_double), pointer :: fp_local(:)
    type(c_ptr) :: cp
    integer(fgsl_size_t) :: size, stride
! tests
!    real(fgsl_double) :: cc(3)
    size = fgsl_aux_vector_double_size(fvec%gsl_vector)
    stride = fgsl_aux_vector_double_stride(fvec%gsl_vector)
    if (stride == 0) then
       fgsl_vector_pointer_align = fgsl_einval
    else
       cp = gsl_vector_ptr(fvec%gsl_vector,0_fgsl_size_t)
!       cc(1) = gsl_vector_get(fvec%gsl_vector,0_c_size_t)
!       cc(2) = gsl_vector_get(fvec%gsl_vector,1_c_size_t)
!       cc(3) = gsl_vector_get(fvec%gsl_vector,2_c_size_t)
       call c_f_pointer(cp, fp_local, (/ size*stride /))
!       write(6, *) 'size, stride, fp_local: ',size,stride,fp_local(1:3),cc(1:3)
       ptr => fp_local(1:size*stride:stride)
       fgsl_vector_pointer_align = fgsl_success
    end if
  end function fgsl_vector_pointer_align
!> The assignment operator (see interface/generics.finc) is overloaded to enable
!> copying of the content of a GSL vector into a Fortran array.
  subroutine fgsl_vector_to_array(result, source)
    real(fgsl_double), intent(inout) :: result(:)
    type(fgsl_vector), intent(in) :: source
!
    integer(fgsl_size_t) :: i, n, k
    k = size(result)
    n = min(k,fgsl_aux_vector_double_size(source%gsl_vector))
!    write(6,*) 'result length: ',size(result)
!    write(6,*) 'vector length: ', &
!    fgsl_aux_vector_double_size(source%gsl_vector)
    do i=1,n
       result(i) = gsl_vector_get(source%gsl_vector,i-1)
    end do
    do i=n+1,size(result)
       result(i) = 0.0_fgsl_double
    end do
  end subroutine fgsl_vector_to_array
!> Free the resources inside a GSL vector object previously established
!> by a call to fgsl_vector_init(). This is invoked via the generic
!> fgsl_vector_free.
  subroutine fgsl_vector_free(fvec)
    type(fgsl_vector), intent(inout) :: fvec
!    call gsl_vector_free(fvec%gsl_vector)
    call fgsl_aux_vector_double_free(fvec%gsl_vector)
  end subroutine fgsl_vector_free
  subroutine fgsl_vector_int_free(fvec)
    type(fgsl_vector_int), intent(inout) :: fvec
    call fgsl_aux_vector_int_free(fvec%gsl_vector_int)
  end subroutine fgsl_vector_int_free
  subroutine fgsl_vector_c_ptr(res, src)
    type(c_ptr), intent(in) :: src
    type(fgsl_vector), intent(out) :: res
    res%gsl_vector = src
  end subroutine fgsl_vector_c_ptr
  function fgsl_vector_status(vector)
    type(fgsl_vector), intent(in) :: vector
    logical :: fgsl_vector_status
    fgsl_vector_status = .true.
    if (.not. c_associated(vector%gsl_vector)) fgsl_vector_status = .false.
  end function fgsl_vector_status
!> Inquire the size of a double precision real GSL vector object.
  function fgsl_vector_int_status(vector)
    type(fgsl_vector_int), intent(in) :: vector
    logical :: fgsl_vector_int_status
    fgsl_vector_int_status = .true.
    if (.not. c_associated(vector%gsl_vector_int)) fgsl_vector_int_status = .false.
  end function fgsl_vector_int_status
  function fgsl_sizeof_vector(w)
    type(fgsl_vector), intent(in) :: w
    integer(fgsl_size_t) :: fgsl_sizeof_vector
    fgsl_sizeof_vector = gsl_aux_sizeof_vector()
  end function fgsl_sizeof_vector
!
! vectors (complex)
!
!> Initialize a complex GSL vector object. This is invoked via the generic
!> fgsl_vector_init.
!> \param type - determine intrinsic type of vector object
!> \return new object of type fgsl_vector
  function fgsl_vector_complex_init_legacy(type)
    complex(fgsl_double_complex), intent(in) :: type
    type(fgsl_vector_complex) :: fgsl_vector_complex_init_legacy
    fgsl_vector_complex_init_legacy%gsl_vector_complex = fgsl_aux_vector_complex_init()
  end function fgsl_vector_complex_init_legacy
  function fgsl_vector_complex_init(array, stride, stat)
    complex(fgsl_double), target, contiguous, intent(in) :: array(:)
    integer(fgsl_size_t), intent(in), optional :: stride
    integer(fgsl_int), intent(inout), optional :: stat
    type(fgsl_vector_complex) :: fgsl_vector_complex_init
    integer(fgsl_size_t) :: stride_local, array_size, section_size
    integer(fgsl_int) :: stat_local

    try : block

      if (present(stride)) then
         stride_local = stride
      else
         stride_local = 1_fgsl_size_t
      end if
      if (stride_local <= 0) then
         stat_local = fgsl_einval
         exit try
      end if

      array_size = size(array,dim=1,kind=fgsl_size_t)
      section_size = (array_size - 1) / stride_local + 1
      
      fgsl_vector_complex_init%gsl_vector_complex = fgsl_aux_vector_complex_init()
      stat_local = fgsl_aux_vector_complex_align(c_loc(array), &
           array_size, fgsl_vector_complex_init%gsl_vector_complex, &
           section_size, 0_fgsl_size_t, stride_local)
    end block try
    if ( present(stat) ) stat = stat_local
    if ( .not. present(stat) .and. stat_local /= fgsl_success ) &
         call fgsl_error("aligning failed", __FILE__, __LINE__, stat_local)
  end function fgsl_vector_complex_init
!> Wrap a rank 1 Fortran array slice inside a double precision complex
!> real GSL vector object. This is invoked via the generic
!> fgsl_vector_align.
!> \param array - requires the actual argument to have the
!>  TARGET attribute. Otherwise being passed by reference is
!>  not guaranteed by the Fortran standard.
!> \param len - number of elements of the rank 1 array
!> \param fvec - previously initialized complex GSL vector object
!> \param size - number of elements from array wrapped inside fvec
!> \param offset - index of first element of array to be mapped to fvec
!> \param stride - stride in array for successive elements of fvec
!> \return Status
  function fgsl_vector_complex_align(array, len, fvec, size, offset, stride)
    integer(fgsl_size_t), intent(in) :: len, size, offset, stride
    complex(fgsl_double_complex), dimension(len), target, intent(in) :: array
    type(fgsl_vector_complex), intent(inout) :: fvec
    integer(fgsl_int) :: fgsl_vector_complex_align
!
    fgsl_vector_complex_align = &
         fgsl_aux_vector_complex_align(c_loc(array), len, &
         fvec%gsl_vector_complex, size, offset, stride)
  end function fgsl_vector_complex_align
!> Associate a Fortran pointer with the data stored inside a GSL vector object.
!> This is invoked via the generic fgsl_vector_align. Objects of type
!> <CODE>gsl_vector_complex</CODE> which are returned by GSL routines often are
!> persistent subobjects of other GSL objects. A Fortran pointer aligned with
!> a subobject hence will remain up-to-date throughout the lifetime of the
!> object; it may become undefined once the object ceases to exist.
!> \param ptr - rank 1 Fortran pointer
!> \param fvec - double precision complex GSL vector
!> \return Status
  function fgsl_vector_complex_pointer_align(ptr, fvec)
    complex(fgsl_double_complex), pointer, intent(out) :: ptr(:)
    type(fgsl_vector_complex), intent(in) :: fvec
    integer(fgsl_int) :: fgsl_vector_complex_pointer_align
!
    complex(fgsl_double_complex), pointer :: fp_local(:)
    type(c_ptr) :: cp
    integer(fgsl_size_t) :: size, stride
! tests
!    real(fgsl_double) :: cc(3)
    size = fgsl_aux_vector_complex_size(fvec%gsl_vector_complex)
    stride = fgsl_aux_vector_complex_stride(fvec%gsl_vector_complex)
    if (stride == 0) then
       fgsl_vector_complex_pointer_align = fgsl_einval
    else
       cp = gsl_vector_complex_ptr(fvec%gsl_vector_complex,0_fgsl_size_t)
!       cc(1) = gsl_vector_complex_get(fvec%gsl_vector_complex,0_c_size_t)
!       cc(2) = gsl_vector_complex_get(fvec%gsl_vector_complex,1_c_size_t)
!       cc(3) = gsl_vector_complex_get(fvec%gsl_vector_complex,2_c_size_t)
       call c_f_pointer(cp, fp_local, (/ size*stride /))
!       write(6, *) 'size, stride, fp_local: ',size,stride,fp_local(1:3),cc(1:3)
       ptr => fp_local(1:size*stride:stride)
       fgsl_vector_complex_pointer_align = fgsl_success
    end if
  end function fgsl_vector_complex_pointer_align
  function fgsl_vector_complex_to_fptr(fvec)
    type(fgsl_vector_complex), intent(in) :: fvec
    complex(fgsl_double), pointer :: fgsl_vector_complex_to_fptr(:)
    complex(fgsl_double), pointer, contiguous :: fptr_local(:)
    integer(fgsl_size_t) :: size, stride
    type(c_ptr) :: cp

    if ( fgsl_vector_complex_status(fvec) ) then
       size = fgsl_aux_vector_complex_size(fvec%gsl_vector_complex)
       stride = fgsl_aux_vector_complex_stride(fvec%gsl_vector_complex)
       if (stride == 0) then
          fgsl_vector_complex_to_fptr => null()
       else
          cp = gsl_vector_complex_ptr(fvec%gsl_vector_complex,0_fgsl_size_t)
          call c_f_pointer(cp, fptr_local, [ size*stride ])
          fgsl_vector_complex_to_fptr => fptr_local(1:size*stride:stride)
       end if
    else
       fgsl_vector_complex_to_fptr => null()
    end if
  end function fgsl_vector_complex_to_fptr
!> The assignment operator (see interface/generics.finc) is overloaded to enable
!> copying of the content of a complex GSL vector into a Fortran array.
  subroutine fgsl_vector_complex_to_array(result, source)
    complex(fgsl_double_complex), intent(inout) :: result(:)
    type(fgsl_vector_complex), intent(in) :: source
!    type(gsl_complex) :: aux
!
    integer(fgsl_size_t) :: i, n, k
    k = size(result)
    n = min(k,fgsl_aux_vector_complex_size(source%gsl_vector_complex))
!    write(6,*) 'result length: ',size(result)
!    write(6,*) 'vector_complex length: ', &
!         fgsl_aux_vector_complex_size(source%gsl_vector_complex)
    do i=1,n
       result(i) = gsl_vector_complex_get(source%gsl_vector_complex,i-1)
!       aux = gsl_vector_complex_get(source%gsl_vector_complex,i-1)
!       result(i) = aux
!       write(6, *) 'i=',i,' res = ',result(i)
    end do
    do i=n+1,size(result)
       result(i) = 0.0_fgsl_double
    end do
  end subroutine fgsl_vector_complex_to_array
!> Free the resources inside a complex GSL vector object previously established
!> by a call to fgsl_vector_complex_init(). This is invoked via the generic
!> fgsl_vector_free.
  subroutine fgsl_vector_complex_free(fvec)
    type(fgsl_vector_complex), intent(inout) :: fvec
    call fgsl_aux_vector_complex_free(fvec%gsl_vector_complex)
  end subroutine fgsl_vector_complex_free
  subroutine fgsl_vector_complex_c_ptr(res, src)
    type(c_ptr), intent(in) :: src
    type(fgsl_vector_complex), intent(out) :: res
    res%gsl_vector_complex = src
  end subroutine fgsl_vector_complex_c_ptr
  function fgsl_vector_complex_status(vector_complex)
    type(fgsl_vector_complex), intent(in) :: vector_complex
    logical :: fgsl_vector_complex_status
    fgsl_vector_complex_status = .true.
    if (.not. c_associated(vector_complex%gsl_vector_complex)) fgsl_vector_complex_status = .false.
  end function fgsl_vector_complex_status
!> Inquire the size of a double precision complex GSL vector object.
  function fgsl_sizeof_vector_complex(w)
    type(fgsl_vector_complex), intent(in) :: w
    integer(fgsl_size_t) :: fgsl_sizeof_vector_complex
    fgsl_sizeof_vector_complex = gsl_aux_sizeof_vector_complex()
  end function fgsl_sizeof_vector_complex
!
! matrices (real)
!
!> Legacy function to initialize a GSL matrix object. This is invoked via the generic
!> fgsl_matrix_init.
!> \param type - determine intrinsic type of vector object
!> \return new object of type fgsl_matrix.
   function fgsl_matrix_init_legacy(type)
    real(fgsl_double), intent(in) :: type
    type(fgsl_matrix) :: fgsl_matrix_init_legacy
    fgsl_matrix_init_legacy%gsl_matrix = fgsl_aux_matrix_double_init()
  end function fgsl_matrix_init_legacy
!> Initialize a rank 2 Fortran array to become associated with a double precision
!> GSL matrix object. This is invoked via the generic fgsl_matrix_init.
!> \param array - requires the actual argument to have the
!>  TARGET and CONTIGUOUS attributes. 
!> \param n - number of rows in array
!> \param m - number of columns in array
!> \param fmat - double precision GSL matrix object, which is allocated
!> \return Status
  function fgsl_matrix_init(array, n, m, stat)
    integer(fgsl_size_t), intent(in), optional :: n, m
    real(fgsl_double), dimension(:,:), target, contiguous, intent(in) :: array
    type(fgsl_matrix) :: fgsl_matrix_init
    integer(fgsl_int), optional :: stat
    integer :: stat_local
    integer(fgsl_size_t) :: mloc, nloc
!
    fgsl_matrix_init%gsl_matrix = fgsl_aux_matrix_double_init()
    if ( present(n) ) then
       nloc = n
    else
       nloc =  size(array,1,KIND=c_size_t)
    end if
    if ( present(m) ) then
       mloc = m
    else
       mloc =  size(array,2,KIND=c_size_t)
    end if
    stat_local = &
         fgsl_aux_matrix_double_align(c_loc(array), size(array,1,KIND=c_size_t), &
         nloc, mloc, fgsl_matrix_init%gsl_matrix)
   if ( present(stat) ) stat = stat_local
   if ( .not. present(stat) .and. stat_local /= fgsl_success ) &
        call fgsl_error("aligning failed", __FILE__, __LINE__, stat_local)    
  end function fgsl_matrix_init
!> Legacy specific to wrap a rank 2 Fortran array inside a double precision
!> real GSL matrix object. This is invoked via the generic
!> fgsl_matrix_align.
!> \param array - requires the actual argument to have the
!>  TARGET attribute. Otherwise being passed by reference is
!>  not guaranteed by the Fortran standard.
!> \param lda - leading dimension of the rank 2 array
!> \param n - number of rows in array
!> \param m - number of columns in array
!> \param fmat - previously initialized double precision GSL matrix object
!> \return Status
   function fgsl_matrix_align(array, lda, n, m, fmat)
    integer(fgsl_size_t), intent(in) :: lda, n, m
    real(fgsl_double), dimension(lda, m), target, intent(in) :: array
    type(fgsl_matrix), intent(inout) :: fmat
    integer(fgsl_int) :: fgsl_matrix_align
!
    fgsl_matrix_align = fgsl_aux_matrix_double_align(c_loc(array), lda, &
         n, m, fmat%gsl_matrix)
  end function fgsl_matrix_align
!> Associate a Fortran pointer with the data stored inside a GSL matrix object.
!> This is invoked via the generic fgsl_matrix_align. Objects of type
!> <CODE>gsl_matrix</CODE> which are returned by GSL routines often are
!> persistent subobjects of other GSL objects. A Fortran pointer aligned with
!> a subobject hence will remain up-to-date throughout the lifetime of the
!> object; it may become undefined once the object ceases to exist.
!> \param ptr - rank 2 Fortran pointer
!> \param fmat - double precision real GSL matrix
!> \return Status
  function fgsl_matrix_pointer_align(ptr, fmat)
    real(fgsl_double), pointer, intent(out) :: ptr(:,:)
    type(fgsl_matrix), intent(in) :: fmat
    integer(fgsl_int) :: fgsl_matrix_pointer_align
!
    real(fgsl_double), pointer :: fp_local(:,:)
    type(c_ptr) :: cp
    integer(fgsl_size_t) :: m, n, lda
    call fgsl_aux_matrix_double_size(fmat%gsl_matrix, lda, n, m)
    cp = gsl_matrix_ptr(fmat%gsl_matrix,0_fgsl_size_t,0_fgsl_size_t)
    call c_f_pointer(cp, fp_local, (/ lda , m /))
    ptr => fp_local(1:n,1:m)
    fgsl_matrix_pointer_align = fgsl_success
  end function fgsl_matrix_pointer_align
!> Associate a Fortran pointer with the data stored inside a GSL matrix object.
!> This is invoked via the generic fgsl_matrix_to_fptr. Objects of type
!> <CODE>gsl_matrix</CODE> which are returned by GSL routines often are
!> persistent subobjects of other GSL objects. A Fortran pointer aligned with
!> a subobject hence will remain up-to-date throughout the lifetime of the
!> object; it may become undefined once the object ceases to exist.
!> \param fmat - GSL matrix
!> \return rank 2 Fortran pointer
  function fgsl_matrix_to_fptr(fmat)
    real(fgsl_double), pointer :: fgsl_matrix_to_fptr(:,:)
    type(fgsl_matrix), intent(in) :: fmat
!
    real(fgsl_double), pointer :: fp_local(:,:)
    type(c_ptr) :: cp
    integer(fgsl_size_t) :: m, n, lda
    if ( fgsl_matrix_status(fmat) ) then
       call fgsl_aux_matrix_double_size(fmat%gsl_matrix, lda, n, m)
       cp = gsl_matrix_ptr(fmat%gsl_matrix,0_fgsl_size_t,0_fgsl_size_t)
       call c_f_pointer(cp, fp_local, (/ lda , m /))
       fgsl_matrix_to_fptr => fp_local(1:n,1:m)
    else
       fgsl_matrix_to_fptr => null()
    end if
  end function fgsl_matrix_to_fptr
!> The assignment operator (see interface/generics.finc) is overloaded to enable
!> copying of the content of a GSL matrix into a rank 2 Fortran array.
  subroutine fgsl_matrix_to_array(result, source)
    real(fgsl_double), intent(inout) :: result(:,:)
    type(fgsl_matrix), intent(in) :: source
!
    integer(fgsl_size_t) :: i, j, kl, m, n, ml, nl, lda
    call fgsl_aux_matrix_double_size(source%gsl_matrix, lda, n, m)

    kl = size(result,1)
    nl = min(kl,n)
    kl = size(result,2)
    ml = min(kl,m)
!    write(6, *) 'Number of rows: ', nl, n
!    write(6, *) 'Number of cols: ', ml, m
    do j=1,ml
       do i=1,nl
          result(i,j) = gsl_matrix_get(source%gsl_matrix,j-1,i-1)
       end do
    end do
    do j=1,ml
       do i=nl+1,size(result,1)
          result(i,j) = 0.0_fgsl_double
       end do
    end do
    do j=ml+1,size(result,2)
       do i=1,size(result,1)
          result(i,j) = 0.0_fgsl_double
       end do
    end do
  end subroutine fgsl_matrix_to_array
!> Free the resources inside a GSL matrix object previously established
!> by a call to fgsl_matrix_init(). This is invoked via the generic
!> fgsl_matrix_free.
  subroutine fgsl_matrix_free(fvec)
    type(fgsl_matrix), intent(inout) :: fvec
    call fgsl_aux_matrix_double_free(fvec%gsl_matrix)
  end subroutine fgsl_matrix_free
  subroutine fgsl_matrix_c_ptr(res, src)
    type(c_ptr), intent(in) :: src
    type(fgsl_matrix), intent(out) :: res
    res%gsl_matrix = src
  end subroutine fgsl_matrix_c_ptr
  function fgsl_matrix_status(matrix)
    type(fgsl_matrix), intent(in) :: matrix
    logical :: fgsl_matrix_status
    fgsl_matrix_status = .true.
    if (.not. c_associated(matrix%gsl_matrix)) fgsl_matrix_status = .false.
  end function fgsl_matrix_status
!> Inquire the number of elements in a double precision real GSL matrix object.
  function fgsl_sizeof_matrix(w)
    type(fgsl_matrix), intent(in) :: w
    integer(fgsl_size_t) :: fgsl_sizeof_matrix
    fgsl_sizeof_matrix = gsl_aux_sizeof_matrix()
  end function fgsl_sizeof_matrix
!
! matrices (complex)
!
!> Legacy specifit to initialize a GSL matrix object. This is invoked via the generic
!> fgsl_matrix_init.
!> \param type - determine intrinsic type of vector object
!> \return new object of type fgsl_matrix.
  function fgsl_matrix_complex_init_legacy(type)
    complex(fgsl_double_complex), intent(in) :: type
    type(fgsl_matrix_complex) :: fgsl_matrix_complex_init_legacy
    fgsl_matrix_complex_init_legacy%gsl_matrix_complex = fgsl_aux_matrix_complex_init()
  end function fgsl_matrix_complex_init_legacy
!> Initialize a rank 2 Fortran array to become associated with a double precision
!> complex GSL matrix object. This is invoked via the generic
!> fgsl_matrix_init.
!> \param array - requires the actual argument to have the
!>  TARGET and CONTIGUOUS attributes. 
!> \param n - number of rows in array
!> \param m - number of columns in array
!> \param fmat - double precision complex GSL matrix object, which is allocated
!> \return Status
  function fgsl_matrix_complex_init(array, n, m, stat)
    integer(fgsl_size_t), intent(in), optional :: n, m
    complex(fgsl_double_complex), dimension(:,:), target, contiguous, intent(in) :: array
    type(fgsl_matrix_complex) :: fgsl_matrix_complex_init
    integer(fgsl_int), optional :: stat
    integer(fgsl_int) :: stat_local
    integer(fgsl_size_t) :: mloc, nloc
!
    fgsl_matrix_complex_init%gsl_matrix_complex = fgsl_aux_matrix_complex_init()
    if ( present(n) ) then
       nloc = n
    else
       nloc =  size(array,1,KIND=c_size_t)
    end if
    if ( present(m) ) then
       mloc = m
    else
       mloc =  size(array,2,KIND=c_size_t)
    end if
    stat_local = fgsl_aux_matrix_complex_align(c_loc(array), size(array,1,KIND=c_size_t), &
         nloc, mloc, fgsl_matrix_complex_init%gsl_matrix_complex)
    if ( present(stat) ) stat = stat_local
    if ( .not. present(stat) .and. stat_local /= fgsl_success ) &
         call fgsl_error("aligning failed", __FILE__, __LINE__, stat_local)
  end function fgsl_matrix_complex_init
!> Legacy function to wrap a rank 2 Fortran array inside a double precision
!> complex GSL matrix object. This is invoked via the generic
!> fgsl_matrix_align.
!> \param array - requires the actual argument to have the
!>  TARGET attribute. Otherwise being passed by reference is
!>  not guaranteed by the Fortran standard.
!> \param lda - leading dimension of the rank 2 array
!> \param n - number of rows in array
!> \param m - number of columns in array
!> \param fmat - previously initialized double precision complex GSL matrix object
!> \return Status
  function fgsl_matrix_complex_align(array, lda, n, m, fmat)
    integer(fgsl_size_t), intent(in) :: lda, n, m
    complex(fgsl_double_complex), dimension(lda, m), target, intent(in) :: array
    type(fgsl_matrix_complex), intent(inout) :: fmat
    integer(fgsl_int) :: fgsl_matrix_complex_align
!
    fgsl_matrix_complex_align = &
         fgsl_aux_matrix_complex_align(c_loc(array), lda, &
         n, m, fmat%gsl_matrix_complex)
  end function fgsl_matrix_complex_align
!> Associate a Fortran pointer with the data stored inside a complex GSL matrix object.
!> This is invoked via the generic fgsl_matrix_align. Objects of type
!> <CODE>gsl_matrix_complex</CODE> which are returned by GSL routines often are
!> persistent subobjects of other GSL objects. A Fortran pointer aligned with
!> a subobject hence will remain up-to-date throughout the lifetime of the
!> object; it may become undefined once the object ceases to exist.
!> \param ptr - rank 2 Fortran pointer
!> \param fmat - double precision complex GSL matrix
!> \return Status
  function fgsl_matrix_complex_pointer_align(ptr, fmat)
    complex(fgsl_double_complex), pointer, intent(out) :: ptr(:,:)
    type(fgsl_matrix_complex), intent(in) :: fmat
    integer(fgsl_int) :: fgsl_matrix_complex_pointer_align
!
    complex(fgsl_double_complex), pointer :: fp_local(:,:)
    type(c_ptr) :: cp
    integer(fgsl_size_t) :: m, n, lda
    call fgsl_aux_matrix_complex_size(fmat%gsl_matrix_complex, lda, n, m)
    cp = gsl_matrix_complex_ptr(fmat%gsl_matrix_complex,0_fgsl_size_t,0_fgsl_size_t)
    call c_f_pointer(cp, fp_local, (/ lda , m /))
    ptr => fp_local(1:n,1:m)
    fgsl_matrix_complex_pointer_align = fgsl_success
  end function fgsl_matrix_complex_pointer_align
  function fgsl_matrix_complex_to_fptr(fmat)
    complex(fgsl_double), pointer :: fgsl_matrix_complex_to_fptr(:,:)
    type(fgsl_matrix_complex), intent(in) :: fmat
!
    complex(fgsl_double), pointer :: fp_local(:,:)
    type(c_ptr) :: cp
    integer(fgsl_size_t) :: m, n, lda
    if ( fgsl_matrix_complex_status(fmat) ) then
       call fgsl_aux_matrix_complex_size(fmat%gsl_matrix_complex, lda, n, m)
       cp = gsl_matrix_complex_ptr(fmat%gsl_matrix_complex,0_fgsl_size_t,0_fgsl_size_t)
       call c_f_pointer(cp, fp_local, (/ lda , m /))
       fgsl_matrix_complex_to_fptr => fp_local(1:n,1:m)
    else
       fgsl_matrix_complex_to_fptr => null()
    end if
  end function fgsl_matrix_complex_to_fptr
  !> The assignment operator (see interface/generics.finc) is overloaded to enable
  !> copying of the content of a complex GSL matrix into a rank 2 Fortran array.
  subroutine fgsl_matrix_complex_to_array(result, source)
    complex(fgsl_double_complex), intent(inout) :: result(:,:)
    type(fgsl_matrix_complex), intent(in) :: source
!
    integer(fgsl_size_t) :: i, j, kl, m, n, ml, nl, lda
    call fgsl_aux_matrix_complex_size(source%gsl_matrix_complex, lda, n, m)

    kl = size(result,1)
    nl = min(kl,n)
    kl = size(result,2)
    ml = min(kl,m)
!    write(6, *) 'Number of rows: ', nl, n
!    write(6, *) 'Number of cols: ', ml, m
    do j=1,ml
       do i=1,nl
          result(i,j) = gsl_matrix_complex_get(source%gsl_matrix_complex,j-1,i-1)
       end do
    end do
    do j=1,ml
       do i=nl+1,size(result,1)
          result(i,j) = 0.0_fgsl_double_complex
       end do
    end do
    do j=ml+1,size(result,2)
       do i=1,size(result,1)
          result(i,j) = 0.0_fgsl_double_complex
       end do
    end do
  end subroutine fgsl_matrix_complex_to_array
!> Free the resources inside a complex GSL matrix object previously established
!> by a call to fgsl_matrix_complex_init(). This is invoked via the generic
!> fgsl_matrix_free.
 subroutine fgsl_matrix_complex_free(fvec)
    type(fgsl_matrix_complex), intent(inout) :: fvec
    call fgsl_aux_matrix_complex_free(fvec%gsl_matrix_complex)
  end subroutine fgsl_matrix_complex_free
  subroutine fgsl_matrix_complex_c_ptr(res, src)
    type(c_ptr), intent(in) :: src
    type(fgsl_matrix_complex), intent(out) :: res
    res%gsl_matrix_complex = src
  end subroutine fgsl_matrix_complex_c_ptr
  function fgsl_matrix_complex_status(matrix_complex)
    type(fgsl_matrix_complex), intent(in) :: matrix_complex
    logical :: fgsl_matrix_complex_status
    fgsl_matrix_complex_status = .true.
    if (.not. c_associated(matrix_complex%gsl_matrix_complex)) fgsl_matrix_complex_status = .false.
  end function fgsl_matrix_complex_status
!> Inquire the number of elements in a double precision complex GSL matrix object.
  function fgsl_sizeof_matrix_complex(w)
    type(fgsl_matrix_complex), intent(in) :: w
    integer(fgsl_size_t) :: fgsl_sizeof_matrix_complex
    fgsl_sizeof_matrix_complex = gsl_aux_sizeof_matrix_complex()
  end function fgsl_sizeof_matrix_complex
  function fgsl_vector_get_size(vec)
    type(fgsl_vector), intent(in) :: vec
    integer(fgsl_size_t) :: fgsl_vector_get_size
    fgsl_vector_get_size = fgsl_aux_vector_double_size(vec%gsl_vector)
  end function fgsl_vector_get_size
  function fgsl_vector_get_stride(vec)
    type(fgsl_vector), intent(in) :: vec
    integer(fgsl_size_t) :: fgsl_vector_get_stride
    fgsl_vector_get_stride = fgsl_aux_vector_double_stride(vec%gsl_vector)
  end function fgsl_vector_get_stride
  function fgsl_matrix_get_size1(matr)
    type(fgsl_matrix), intent(in) :: matr
    integer(fgsl_size_t) :: tda, size1, size2
    integer(fgsl_size_t) :: fgsl_matrix_get_size1
    call fgsl_aux_matrix_double_size(matr%gsl_matrix, tda, size2,&
    size1)
    fgsl_matrix_get_size1 = size1
  end function fgsl_matrix_get_size1
  function fgsl_matrix_get_size2(matr)
    type(fgsl_matrix), intent(in) :: matr
    integer(fgsl_size_t) :: tda, size1, size2
    integer(fgsl_size_t) :: fgsl_matrix_get_size2
    call fgsl_aux_matrix_double_size(matr%gsl_matrix, tda, size2,&
    size1)
    fgsl_matrix_get_size2 = size2
  end function fgsl_matrix_get_size2
  function fgsl_matrix_get_tda(matr)
    type(fgsl_matrix), intent(in) :: matr
    integer(fgsl_size_t) :: tda, size1, size2
    integer(fgsl_size_t) :: fgsl_matrix_get_tda
    call fgsl_aux_matrix_double_size(matr%gsl_matrix, tda, size2,&
    size1)
    fgsl_matrix_get_tda = tda
  end function fgsl_matrix_get_tda