Merge pull request #30 from reinh-bader/fgsl_devel_1.5.0

changes for FGSL 1.5.0

.travis.yml

@@ -19,10 +19,6 @@ addons:
 
 env:
   matrix:
-    - GCC=6 GSL=2.5
-    - GCC=7 GSL=2.5
-    - GCC=8 GSL=2.5
-    - GCC=9 GSL=2.5
     - GCC=6 GSL=2.6
     - GCC=7 GSL=2.6
     - GCC=8 GSL=2.6
@@ -31,12 +27,12 @@ env:
     - LD_LIBRARY_PATH=/usr/local/lib
     - LD_RUN_PATH=/usr/local/lib
     - PKG_CONFIG_LIBDIR=/usr/local/lib/pkgconfig
-matrix:
-  allow_failures:
-    - env: GCC=6 GSL=2.6
-    - env: GCC=7 GSL=2.6
-    - env: GCC=8 GSL=2.6
-    - env: GCC=9 GSL=2.6
+#matrix:
+#  allow_failures:
+#    - env: GCC=6 GSL=2.6
+#    - env: GCC=7 GSL=2.6
+#    - env: GCC=8 GSL=2.6
+#    - env: GCC=9 GSL=2.6
 before_install:
   - |
     if [ $TRAVIS_OS_NAME == osx ] ; then 

NEWS

@@ -1,3 +1,23 @@
+Changes in 1.5.0 (by R Bader):
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+This release introduces support for gsl 2.6. An overview:
+
+1. added linalg LQ and LDLT routines
+2. augmented the legacy triangular matrix interface by the
+   better designed one currently documented.
+   Note: the old interface will be eventually removed.
+3. added interfaces for banded systems (cholesky, LDLT)
+4. added missing calls from spmatrix module (note: some
+   remain unsupported).
+5. Updated specfunc Hermite polynomial and function interfaces.
+6. Added recursive QR solvers to linalg
+
+Comment:
+
+*  the gsl_bst module is not supported in this release. Note 
+   that the interface is not officially documented.
+
+
 Changes in 1.4.0 (by R Bader):
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 This release introduces support for gsl 2.5. An overview:

README

@@ -1,13 +1,14 @@
 
-This is release 1.4.0 of the Fortran interface 
+This is release 1.5.0 of the Fortran interface 
 to the GNU Scientific Library. 
-It is based on version 2.5 of GSL, and should not be used with 
+It is based on version 2.6 of GSL, and should not be used with 
 earlier GSL releases.
 
 Please use the release 1.0.0 in combination with GSL 1.x
 Please use the release 1.1.0 in combination with GSL 2.1 or 2.2.1
 Please use the release 1.2.0 in combination with GSL 2.3
 Please use the release 1.3.0 in combination with GSL 2.4
+Please use the release 1.4.0 in combination with GSL 2.5
 
 Please consult the NEWS file for a change log.
 
@@ -128,3 +129,4 @@ Releases:
   * 1.2.0: January, 2017
   * 1.3.0: August, 2019
   * 1.4.0: March, 2021
+  * 1.5.0: TBD

README_WORKFLOW

@@ -32,3 +32,16 @@ Some notes on the source management workflow
 * if needed, keep making changes to the code, commit them and check the build results.
 * when done, merge the PR in and delete the local and remote working branch.
 
+(C) Documentation rebuild
+
+* Use Ubuntu-based system
+* Delete contents of html and latex folders in doc via "git rm -r *"
+* Run doxygen Doxyfile in base dir (may want to consider updates to Doxyfile)
+* Go to doc/latex and run make refman.pdf. Copy the file to ..
+* Run git add html in doc subdirectory.
+
+(D) packaging and release
+
+* configure.ac contains the version info and may need an update
+* create a tarball after running configure: make dist
+

api/array.finc

@@ -255,10 +255,42 @@
 !> fgsl_vector_init.
 !> \param type - determine intrinsic type of vector object
 !> \return new object of type fgsl_vector
-  function fgsl_vector_complex_init(type)
+  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
-    fgsl_vector_complex_init%gsl_vector_complex = fgsl_aux_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
@@ -316,6 +348,27 @@
        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)
@@ -366,16 +419,50 @@
 !
 ! matrices (real)
 !
-!> Initialize a GSL matrix object. This is invoked via the generic
+!> 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(type)
+   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
-!> Wrap a rank 2 Fortran array inside a double precision
+!> 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
@@ -418,6 +505,30 @@
     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)
@@ -476,16 +587,50 @@
 !
 ! matrices (complex)
 !
-!> Initialize a GSL matrix object. This is invoked via the generic
+!> 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(type)
+  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
-!> Wrap a rank 2 Fortran array inside a double precision
+!> 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
@@ -529,8 +674,24 @@
     ptr => fp_local(1:n,1:m)
     fgsl_matrix_complex_pointer_align = fgsl_success
   end function fgsl_matrix_complex_pointer_align
-!> 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.
+  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

api/interp.finc

@@ -398,6 +398,18 @@
     fgsl_interp2d_eval_e_extrap = gsl_interp2d_eval_e_extrap(interp%gsl_interp2d, c_loc(xarr), &
     c_loc(yarr), c_loc(zarr), x, y, xa%gsl_interp_accel, ya%gsl_interp_accel, z)
   end function fgsl_interp2d_eval_e_extrap
+  function fgsl_interp2d_eval_extrap_e(interp, xarr, yarr, zarr, x, y, xa, ya, z)
+    type(fgsl_interp2d), intent(in) :: interp
+    real(fgsl_double), dimension(:), intent(in), target, contiguous :: xarr, yarr
+    real(fgsl_double), dimension(:,:), intent(in), target, contiguous :: zarr
+    real(fgsl_double), intent(in) :: x, y
+    type(fgsl_interp_accel), intent(inout) :: xa, ya
+    real(fgsl_double), intent(out) :: z
+    integer(fgsl_int) :: fgsl_interp2d_eval_extrap_e
+    ! I could do more bounds checking here, but seems like overkill
+    fgsl_interp2d_eval_extrap_e = gsl_interp2d_eval_extrap_e(interp%gsl_interp2d, c_loc(xarr), &
+    c_loc(yarr), c_loc(zarr), x, y, xa%gsl_interp_accel, ya%gsl_interp_accel, z)
+  end function fgsl_interp2d_eval_extrap_e
   function fgsl_interp2d_eval_deriv_x(interp, xarr, yarr, zarr, x, y, xa, ya)
     type(fgsl_interp2d), intent(in) :: interp
     real(fgsl_double), dimension(:), intent(in), target, contiguous :: xarr, yarr