lib_ext_math_m Module

Mathematical functions downloaded from John Burkardt page fSolve and brent

  • https://people.sc.fsu.edu/~jburkardt/f_src/fsolve/fsolve.html
  • https://people.math.sc.edu/Burkardt/f_src/brent/brent.html

Modified:

  • 06 April 2010
  • 01 February 2024 - introducing object with virtual fcn

Author:

  • Original FORTRAN77 version by Jorge More, Burton Garbow, Kenneth Hillstrom.
  • FORTRAN90 version by John Burkardt.
  • Fortran03 version by Jacek Kosek and Damien Furfaro

Reference:

  • Jorge More, Burton Garbow, Kenneth Hillstrom,
  • User Guide for MINPACK-1,
  • Technical Report ANL-80-74,
  • Argonne National Laboratory, 1980.

Uses

  • module~~lib_ext_math_m~~UsesGraph module~lib_ext_math_m lib_ext_math_m module~krn_global_tools_m krn_global_tools_m module~lib_ext_math_m->module~krn_global_tools_m ieee_arithmetic ieee_arithmetic module~krn_global_tools_m->ieee_arithmetic iso_c_binding iso_c_binding module~krn_global_tools_m->iso_c_binding iso_fortran_env iso_fortran_env module~krn_global_tools_m->iso_fortran_env

Used by

  • module~~lib_ext_math_m~~UsedByGraph module~lib_ext_math_m lib_ext_math_m module~cmp_boundary_calc_m cmp_boundary_calc_m module~cmp_boundary_calc_m->module~lib_ext_math_m module~lib_he_thermo_m lib_He_thermo_m module~cmp_boundary_calc_m->module~lib_he_thermo_m module~cmp_boundary_init_m cmp_boundary_init_m module~cmp_boundary_calc_m->module~cmp_boundary_init_m module~cmp_circulator_init_m cmp_circulator_init_m module~cmp_circulator_init_m->module~lib_ext_math_m module~cmp_circulator_init_m->module~lib_he_thermo_m module~cmp_junction_calc_m cmp_junction_calc_m module~cmp_junction_calc_m->module~lib_ext_math_m module~cmp_junction_init_m cmp_junction_init_m module~cmp_junction_calc_m->module~cmp_junction_init_m module~cmp_junction_calc_m->module~lib_he_thermo_m module~cmp_junction_init_m->module~lib_ext_math_m module~cmp_junction_init_m->module~lib_he_thermo_m module~lib_he_thermo_m->module~lib_ext_math_m module~cmp_boundary_init_m->module~lib_he_thermo_m module~cmp_channel_flux_m cmp_channel_flux_m module~cmp_channel_flux_m->module~lib_he_thermo_m module~cmp_channel_init_m cmp_channel_init_m module~cmp_channel_flux_m->module~cmp_channel_init_m module~cmp_channel_init_m->module~lib_he_thermo_m module~cmp_channel_source_terms_m cmp_channel_source_terms_m module~cmp_channel_source_terms_m->module~lib_he_thermo_m module~cmp_channel_source_terms_m->module~cmp_channel_init_m module~cmp_circulator_calc_m cmp_circulator_calc_m module~cmp_circulator_calc_m->module~cmp_circulator_init_m module~cmp_ffsrclink_calc_m cmp_FFsrcLink_calc_m module~cmp_ffsrclink_calc_m->module~lib_he_thermo_m program~reims_p reims_p program~reims_p->module~cmp_boundary_calc_m program~reims_p->module~cmp_junction_calc_m program~reims_p->module~cmp_channel_flux_m program~reims_p->module~cmp_circulator_calc_m program~reims_p->module~cmp_ffsrclink_calc_m module~cmp_channel_calc_m cmp_channel_calc_m program~reims_p->module~cmp_channel_calc_m module~cmp_channel_calc_m->module~cmp_channel_init_m module~cmp_channel_calc_m->module~cmp_channel_source_terms_m

Abstract Interfaces

abstract interface

  • private subroutine fSolve_function(me, n, x, fVec)

    User-supplied system of equations

    Arguments

    Type IntentOptional Attributes Name
    class(fSolve_t), intent(in) :: me

    User-defined object transmitted to the solver

    integer, intent(in) :: n

    number of functions and variables
    real(kind=dp), intent(inout) :: x(n)

    result or initial point where the jacobian is evaluated
    real(kind=dp), intent(inout) :: fVec(n)

    the functions evaluated at 'x'

abstract interface

  • private function brent_function(me, x)

    user-supplied function, of the form which evaluates the function whose zero is being sought.

    Arguments

    Type IntentOptional Attributes Name
    class(brent_t), intent(in) :: me

    User-defined object transmitted to the solver

    real(kind=dp), intent(in) :: x

    'x' parameter of the function f(x)

    Return Value real(kind=dp)

    returned value


Derived Types

type, public, abstract ::  fSolve_t

Procedure 'fcn' of this class is the parameter to fSolve subroutine.

Read more…

Type-Bound Procedures

procedure(fSolve_function), public, deferred :: fcn

type, public, abstract ::  brent_t

class with the function f(x)

Type-Bound Procedures

procedure(brent_function), public, deferred :: f

Functions

private function eNorm(n, x)

Computes the Euclidean norm of a vector.

Read more…

Arguments

Type IntentOptional Attributes Name
integer :: n

Length of the vector
real(kind=dp) :: x(n)

Vector whose norm is desired

Return Value real(kind=dp)

Euclidean norm of the vector

public function zero(obj, a, b, machep, t)

Seeks the root of a function F(X) in an interval [A,B].

Read more…

Arguments

Type IntentOptional Attributes Name
class(brent_t), intent(in) :: obj

external object with the user-supplied subroutine which calculates the functions.
real(kind=dp), intent(in) :: a

the endpoints of the change of sign interval
real(kind=dp), intent(in) :: b

the endpoints of the change of sign interval
real(kind=dp), intent(in) :: machep

an estimate for the relative machine precision.
real(kind=dp), intent(in) :: t

a positive error tolerance

Return Value real(kind=dp)

the estimated value of a zero of the function F


Subroutines

private subroutine dogleg(n, r, lr, diag, qtb, delta, x)

Finds the minimizing combination of Gauss-Newton and gradient steps.

Read more…

Arguments

Type IntentOptional Attributes Name
integer :: n

In - order of the matrix R
real(kind=dp) :: r(lr)

In - upper triangular matrix R stored by rows.
integer :: lr

In - size of the R array, which must be no less than (N*(N+1))/2.
real(kind=dp) :: diag(n)

In - diagonal elements of the matrix D.
real(kind=dp) :: qtb(n)

In - first N elements of the vector Q'* B.
real(kind=dp) :: delta

In - positive upper bound on the euclidean norm of D*X(1:N).
real(kind=dp) :: x(n)

Out - desired convex combination of the Gauss-Newton direction and the scaled gradient direction.

private subroutine fdJac1(obj, n, x, fVec, fjac, ldfjac, ml, mu, epsFcn)

Estimates a jacobian matrix using forward differences.

Read more…

Arguments

Type IntentOptional Attributes Name
class(fSolve_t), intent(in) :: obj

External object with the user-supplied subroutine which calculates the functions.
integer, intent(in) :: n

Number of functions and variables
real(kind=dp), intent(inout) :: x(n)

Point where the jacobian is evaluated
real(kind=dp), intent(in) :: fVec(n)

Functions evaluated at 'x'
real(kind=dp), intent(inout) :: fjac(ldfjac,n)

The N by N approximate
integer, intent(in) :: ldfjac

Leading dimension of FJAC, which must not be less than N.
integer, intent(in) :: ml

Number of sub-diagonals within the band of the jacobian matrix If the jacobian is not banded, set ML and MU to N-1.
integer, intent(in) :: mu

Number of super-diagonals within the band of the jacobian matrix If the jacobian is not banded, set ML and MU to N-1.
real(kind=dp), intent(in) :: epsFcn

Is used in determining a suitable step length for the forward-difference approximation.

Read more…

public subroutine fsolve(obj, n, x, fVec, tol, info)

Seeks a zero of N non-linear equations in N variables

Read more…

Arguments

Type IntentOptional Attributes Name
class(fSolve_t), intent(in) :: obj

External object with the user-supplied subroutine which calculates the functions.
integer, intent(in) :: n

Number of functions and variables
real(kind=dp), intent(inout) :: x(n)

Initial estimate of the solution vector / the estimate of the solution vector
real(kind=dp), intent(out) :: fVec(n)

The functions evaluated at the output X
real(kind=dp), intent(in) :: tol

Satisfactory termination occurs when the algorithm estimates that the relative error between X and the solution is at most TOL. TOL should be non-negative.
integer, intent(out) :: info

Status flag

Read more…

private subroutine hybrd(obj, n, x, fVec, xTol, maxFev, ml, mu, epsFcn, diag, mode, factor, info, nFev, fjac, ldfjac, r, lr, qtf)

Seeks a zero of N non-linear equations in N variables.

Read more…

Arguments

Type IntentOptional Attributes Name
class(fSolve_t), intent(in) :: obj

External object with the user-supplied subroutine which calculates the functions.
integer, intent(in) :: n

The number of functions and variables
real(kind=dp), intent(inout) :: x(n)

Input initial estimate of the solution vector / output X final estimate of the solution vector
real(kind=dp), intent(out) :: fVec(n)

The functions evaluated at the output X
real(kind=dp), intent(in) :: xTol

Termination occurs when the relative error between two consecutive iterates is at most XTOL. XTOL should be non-negative.
integer, intent(in) :: maxFev

Termination occurs when the number of calls to FCN is at least MAXFEV by the end of an iteration
integer, intent(in) :: ml

Specify the number of sub-diagonals within the band of the jacobian matrix. If the jacobian is not banded, set ML to at least n - 1.
integer, intent(in) :: mu

Specify the number of super-diagonals within the band of the jacobian matrix. If the jacobian is not banded, set MU to at least n - 1.
real(kind=dp), intent(in) :: epsFcn

is used in determining a suitable step length for the forward-difference approximation. This approximation assumes that the relative errors in the functions are of the order of EPSFCN. If EPSFCN is less than the machine precision, it is assumed that the relative errors in the functions are of the order of the machine precision.
real(kind=dp), intent(inout) :: diag(n)

If MODE = 1, then DIAG is set internally. If MODE = 2, then DIAG must contain positive entries that serve as multiplicative scale factors for the variables.
integer, intent(in) :: mode

variables will be scaled internally. 2, scaling is specified by the input DIAG vector.
real(kind=dp), intent(in) :: factor

determines the initial step bound. This bound is set to the product of FACTOR and the euclidean norm of DIAG*X if non-zero, or else to FACTOR itself. In most cases, FACTOR should lie in the interval (0.1, 100) with 100 the recommended value.
integer, intent(out) :: info

Error flag

Read more…
integer, intent(out) :: nFev

The number of calls to FCN. produced by the QR factorization of the final approximate jacobian.
real(kind=dp), intent(out) :: fjac(ldfjac,n)

Out - N by N array which contains the orthogonal matrix Q
integer, intent(in) :: ldfjac

The leading dimension of FJAC. LDFJAC must be at least N.
real(kind=dp), intent(out) :: r(lr)

Upper triangular matrix produced by the QR factorization of the final approximate jacobian, stored row-wise.
integer, intent(in) :: lr

Size of the R array, which must be no less than (N*(N+1))/2.
real(kind=dp), intent(out) :: qtf(n)

Contains the vector Q'*FVEC.

private subroutine qform(m, n, q, ldq)

Produces the explicit QR factorization of a matrix.

Read more…

Arguments

Type IntentOptional Attributes Name
integer :: m

In - number of rows of A and the order of Q
integer :: n

In - number of columns of A
real(kind=dp) :: q(ldq,m)

In/Out - lower trapezoid in the first min(M,N) / accumulated into a square matrix

Read more…
integer :: ldq

In - not less than M which specifies the leading dimension of the array Q

private subroutine qrFac(m, n, a, lda, pivot, ipvt, lipvt, rdiag, acnorm)

computes a QR factorization using Householder transformations.

Read more…

Arguments

Type IntentOptional Attributes Name
integer :: m

In - number of rows of A.
integer :: n

In - number of columns of A.
real(kind=dp) :: a(lda,n)

In/Out - before/after QR factosiation size should be the M by N array.

Read more…
integer :: lda

In - leading dimension of A, which must be no less than M.
logical :: pivot

Out - is TRUE if column pivoting is to be carried out
integer :: ipvt(lipvt)

Out - defines the permutation matrix P such that AP = QR.

Read more…
integer :: lipvt

In - dimension of IPVT, which should be N if pivoting is used.
real(kind=dp) :: rdiag(n)

Out - contains the diagonal elements of R.
real(kind=dp) :: acnorm(n)

private subroutine r1mpyq(m, n, a, lda, v, w)

computes A*Q, where Q is the product of Householder transformations

Read more…

Arguments

Type IntentOptional Attributes Name
integer :: m

In - number of rows of A
integer :: n

In - number of columns of A
real(kind=dp) :: a(lda,n)

In/Out - the M by N array. On input, the matrix A to be postmultiplied by the orthogonal matrix Q. On output, the value of A*Q.
integer :: lda

In - leading dimension of A, which must not be less than M.
real(kind=dp) :: v(n)

In - contain the information necessary to recover the Givens rotations GV
real(kind=dp) :: w(n)

In - contain the information necessary to recover the Givens rotations GW

private subroutine r1updt(m, n, s, ls, u, v, w, sing)

re-triangularizes a matrix after a rank one update

Read more…

Arguments

Type IntentOptional Attributes Name
integer :: m

In - number of rows of S.
integer :: n

In - number of columns of S. N must not exceed M.

real(kind=dp) :: s(ls)

In/Out - On input, the lower trapezoidal matrix S stored by columns. On output S contains the lower trapezoidal matrix produced as described above.
integer :: ls

In - length of the S array. LS must be at least (N(2M-N+1))/2.
real(kind=dp) :: u(m)

In - the U vector.
real(kind=dp) :: v(n)

In/Out - On input, V must contain the vector V. On output V contains the information necessary to recover the Givens rotations GV described above.
real(kind=dp) :: w(m)

Out - contains information necessary to recover the Givens rotations GW described above.
logical :: sing

Out - is set to TRUE if any of the diagonal elements of the output S are zero. Otherwise SING is set FALSE.