diff --git a/doc/specs/stdlib_specialmatrices.md b/doc/specs/stdlib_specialmatrices.md index 6c86f3bd1..d50bce46a 100644 --- a/doc/specs/stdlib_specialmatrices.md +++ b/doc/specs/stdlib_specialmatrices.md @@ -12,7 +12,8 @@ The `stdlib_specialmatrices` module provides derived types and specialized drive These include: - Tridiagonal matrices -- Symmetric Tridiagonal matrices (not yet supported) +- Symmetric Tridiagonal matrices +- Hermitian Tridiagonal matrices - Circulant matrices (not yet supported) - Toeplitz matrices (not yet supported) - Hankel matrices (not yet supported) @@ -76,6 +77,107 @@ Tridiagonal matrices are available with all supported data types as `tridiagonal {!example/specialmatrices/example_tridiagonal_dp_type.f90!} ``` +### Symmetric tridiagonal matrices {#SymTridiagonal} + +#### Status + +Experimental + +#### Description + +Symmetric tridiagonal matrices are ubiquituous in scientific computing and often appear when discretizing 1D differential operators. +A generic symmetric tridiagonal matrix has the following structure: +$$ + A + = + \begin{bmatrix} + a_1 & b_1 \\ + b_1 & a_2 & b_2 \\ + & \ddots & \ddots & \ddots \\ + & & b_{n-2} & a_{n-1} & b_{n-1} \\ + & & & b_{n-1} & a_n + \end{bmatrix}. +$$ +Hence, only one vector of size `n` and two of size `n-1` need to be stored to fully represent the matrix. +This particular structure also lends itself to specialized implementations for many linear algebra tasks. +Interfaces to the most common ones will soon be provided by `stdlib_specialmatrices`. +Symmetric tridiagonal matrices are available with all supported data types as `symtridiagonal__type`, for example: + +- `symtridiagonal_sp_type` : Symmetric tridiagonal matrix of size `n` with `real`/`single precision` data. +- `symtridiagonal_dp_type` : Symmetric tridiagonal matrix of size `n` with `real`/`double precision` data. +- `symtridiagonal_xdp_type` : Symmetric tridiagonal matrix of size `n` with `real`/`extended precision` data. +- `symtridiagonal_qp_type` : Symmetric tridiagonal matrix of size `n` with `real`/`quadruple precision` data. +- `symtridiagonal_csp_type` : Symmetric tridiagonal matrix of size `n` with `complex`/`single precision` data. +- `symtridiagonal_cdp_type` : Symmetric tridiagonal matrix of size `n` with `complex`/`double precision` data. +- `symtridiagonal_cxdp_type` : Symmetric tridiagonal matrix of size `n` with `complex`/`extended precision` data. +- `symtridiagonal_cqp_type` : Symmetric tridiagonal matrix of size `n` with `complex`/`quadruple precision` data. + +#### Syntax + +- To construct a symmetric tridiagonal matrix from already allocated arrays `dv` (main diagonal, size `n`, only its real part is being referenced) and `ev` (upper diagonal, size `n-1`): + +`A = ` [[stdlib_specialmatrices(module):symtridiagonal(interface)]] `(dv, ev)` + +- To construct a symmetric tridiagonal matrix of size `n x n` with constant diagonal elements `dv`, and `ev`: + +`A = ` [[stdlib_specialmatrices(module):symtridiagonal(interface)]] `(dv, ev, n)` + +#### Example + +```fortran +{!example/specialmatrices/example_symtridiagonal_dp_type.f90!} +``` + +### Hermitian tridiagonal matrices {#HermTridiagonal} + +#### Status + +Experimental + +#### Description + +Hermitian tridiagonal matrices are ubiquituous in scientific computing. +A generic hermitian tridiagonal matrix has the following structure: +$$ + A + = + \begin{bmatrix} + a_1 & b_1 \\ + \bar{b}_1 & a_2 & b_2 \\ + & \ddots & \ddots & \ddots \\ + & & \bar{b}_{n-2} & a_{n-1} & b_{n-1} \\ + & & & \bar{b}_{n-1} & a_n + \end{bmatrix}, +$$ +where \( a_i \in \mathbb{R} \), \( b_i \in \mathbb{C} \) and the overbar denotes the complex conjugate operation. +Hence, only one vector of size `n` and two of size `n-1` need to be stored to fully represent the matrix. +This particular structure also lends itself to specialized implementations for many linear algebra tasks. +Interfaces to the most common ones will soon be provided by `stdlib_specialmatrices`. +Hermitian tridiagonal matrices are available with all supported data types as `hermtridiagonal__type`, for example: + +- `hermtridiagonal_csp_type` : Hermitian tridiagonal matrix of size `n` with `complex`/`single precision` data. +- `hermtridiagonal_cdp_type` : Hermitian tridiagonal matrix of size `n` with `complex`/`double precision` data. +- `hermtridiagonal_cxdp_type` : Hermitian tridiagonal matrix of size `n` with `complex`/`extended precision` data. +- `hermtridiagonal_cqp_type` : Hermitian tridiagonal matrix of size `n` with `complex`/`quadruple precision` data. + +#### Syntax + +- To construct a hermitian tridiagonal matrix from already allocated arrays `dv` (main diagonal, size `n`) and `ev` (upper diagonal, size `n-1`): + +`A = ` [[stdlib_specialmatrices(module):hermtridiagonal(interface)]] `(dv, ev)` + +- To construct a hermitian tridiagonal matrix of size `n x n` with constant diagonal elements `dv`, and `ev`: + +`A = ` [[stdlib_specialmatrices(module):hermtridiagonal(interface)]] `(dv, ev, n)` + +Note that only the real parts of the diagonal elements `dv` are being used to construct the corresponding Hermitian matrix. + +#### Example + +```fortran +{!example/specialmatrices/example_hermtridiagonal_cdp_type.f90!} +``` + ## Specialized drivers for linear algebra tasks Below is a list of all the specialized drivers for linear algebra tasks currently provided by the `stdlib_specialmatrices` module. @@ -111,7 +213,7 @@ With the exception of `extended precision` and `quadruple precision`, all the ty - `op` (optional) : In-place operator identifier. Shall be a character(1) argument. It can have any of the following values: `N`: no transpose, `T`: transpose, `H`: hermitian or complex transpose. @warning -Due to limitations of the underlying `lapack` driver, currently `alpha` and `beta` can only take one of the values `[-1, 0, 1]` for `tridiagonal` and `symtridiagonal` matrices. See `lagtm` for more details. +Due to limitations of the underlying `lapack` driver, currently `alpha` and `beta` can only take one of the values `[-1, 0, 1]` for `tridiagonal`, `symtridiagonal` and `hermtridiagonal` matrices. See `lagtm` for more details. @endwarning #### Examples @@ -192,17 +294,17 @@ Experimental The definition of all standard artihmetic operators have been overloaded to be applicable for the matrix types defined by `stdlib_specialmatrices`: -- Overloading the `+` operator for adding two matrices of the same type and kind. -- Overloading the `-` operator for subtracting two matrices of the same type and kind. +- Overloading the `+` operator for adding two matrices of the same class and kind. +- Overloading the `-` operator for subtracting two matrices of the same class and kind. - Overloading the `*` for scalar-matrix multiplication. #### Syntax -- Adding two matrices of the same type: +- Adding two matrices of the same class: `C = A` [[stdlib_specialmatrices(module):operator(+)(interface)]] `B` -- Subtracting two matrices of the same type: +- Subtracting two matrices of the same class: `C = A` [[stdlib_specialmatrices(module):operator(-)(interface)]] `B` @@ -211,5 +313,5 @@ The definition of all standard artihmetic operators have been overloaded to be a `B = alpha` [[stdlib_specialmatrices(module):operator(*)(interface)]] `A` @note -For addition (`+`) and subtraction (`-`), matrices `A`, `B` and `C` all need to be of the same type and kind. For scalar multiplication (`*`), `A` and `B` need to be of the same type and kind, while `alpha` is either `real` or `complex` (with the same kind again) depending on the type being used. +For addition (`+`) and subtraction (`-`), matrices `A`, and `B` need to be of the same class and kind. For scalar multiplication (`*`), `A` and `B` need to be of the same class and kind, while `alpha` is either `real` or `complex` (with the same kind again) depending on the type being used. @endnote diff --git a/doc/specs/stdlib_system.md b/doc/specs/stdlib_system.md index 96eebb2e8..56c6f70e6 100644 --- a/doc/specs/stdlib_system.md +++ b/doc/specs/stdlib_system.md @@ -174,7 +174,7 @@ The result is a real value representing the elapsed time in seconds, measured fr ### Syntax -`delta_t = ` [[stdlib_system(module):elapsed(subroutine)]] `(process)` +`delta_t = ` [[stdlib_system(module):elapsed(interface)]] `(process)` ### Arguments @@ -212,7 +212,7 @@ in case of process hang or delay. ### Syntax -`call ` [[stdlib_system(module):wait(subroutine)]] `(process [, max_wait_time])` +`call ` [[stdlib_system(module):wait(interface)]] `(process [, max_wait_time])` ### Arguments @@ -243,7 +243,7 @@ This is especially useful for monitoring asynchronous processes and retrieving t ### Syntax -`call ` [[stdlib_system(module):update(subroutine)]] `(process)` +`call ` [[stdlib_system(module):update(interface)]] `(process)` ### Arguments @@ -269,7 +269,7 @@ This interface is useful when a process needs to be forcefully stopped, for exam ### Syntax -`call ` [[stdlib_system(module):kill(subroutine)]] `(process, success)` +`call ` [[stdlib_system(module):kill(interface)]] `(process, success)` ### Arguments @@ -298,7 +298,7 @@ It ensures that the requested sleep duration is honored on both Windows and Unix ### Syntax -`call ` [[stdlib_system(module):sleep(subroutine)]] `(millisec)` +`call ` [[stdlib_system(module):sleep(interface)]] `(millisec)` ### Arguments @@ -324,7 +324,7 @@ This function is highly efficient and works during the compilation phase, avoidi ### Syntax -`result = ` [[stdlib_system(module):is_windows(function)]] `()` +`result = ` [[stdlib_system(module):is_windows(interface)]] `()` ### Return Value @@ -359,7 +359,7 @@ If the OS cannot be identified, the function returns `OS_UNKNOWN`. ### Syntax -`os = [[stdlib_system(module):get_runtime_os(function)]]()` +`os = ` [[stdlib_system(module):get_runtime_os(function)]] `()` ### Class @@ -396,7 +396,7 @@ This caching mechanism ensures negligible overhead for repeated calls, unlike `g ### Syntax -`os = [[stdlib_system(module):OS_TYPE(function)]]()` +`os = ` [[stdlib_system(module):OS_TYPE(function)]]`()` ### Class @@ -418,6 +418,85 @@ Returns one of the `integer` `OS_*` parameters representing the OS type, from th --- +## `FS_ERROR` - Helper function for error handling + +### Status + +Experimental + +### Description + +A helper function for returning the `type(state_type)` with the flag `STDLIB_FS_ERROR` set. + +### Syntax + +`err = FS_ERROR([a1,a2,a3,a4...... a20])` + +### Class +Pure Function + +### Arguments + +`a1,a2,a3.....a20`(optional): They are of type `class(*), dimension(..), optional, intent(in)`. +An arbitrary list of `integer`, `real`, `complex`, `character` or `string_type` variables. Numeric variables may be provided as either scalars or rank-1 (array) inputs. + +### Behavior + +Formats all the arguments into a nice error message, utilizing the constructor of [[stdlib_system(module):state_type(type)]] + +### Return values + +`type(state_type)` + +### Example + +```fortran +{!example/system/example_fs_error.f90!} +``` + +--- + +## `FS_ERROR_CODE` - Helper function for error handling (with error code) + +### Status + +Experimental + +### Description + +A helper function for returning the `type(state_type)` with the flag `STDLIB_FS_ERROR` set. +It also formats and prefixes the `code` passed to it as the first argument. + +### Syntax + +`err = FS_ERROR_CODE(code [, a1,a2,a3,a4...... a19])` + +### Class +Pure Function + +### Arguments + +`code`: An `integer` code. + +`a1,a2,a3.....a19`(optional): They are of type `class(*), dimension(..), optional, intent(in)`. +An arbitrary list of `integer`, `real`, `complex`, `character` or `string_type` variables. Numeric variables may be provided as either scalars or rank-1 (array) inputs. + +### Behavior + +Formats all the arguments into a nice error message, utilizing the constructor of [[stdlib_system(module):state_type(type)]] + +### Return values + +`type(state_type)` + +### Example + +```fortran +{!example/system/example_fs_error.f90!} +``` + +--- + ## `is_directory` - Test if a path is a directory ### Status @@ -431,7 +510,7 @@ It is designed to work across multiple platforms. On Windows, paths with both fo ### Syntax -`result = [[stdlib_system(module):is_directory(function)]] (path)` +`result = ` [[stdlib_system(module):is_directory(function)]]`(path)` ### Class @@ -471,7 +550,7 @@ It reads as an empty file. The null device's path varies by operating system: ### Syntax -`path = [[stdlib_system(module):null_device(function)]]()` +`path = ` [[stdlib_system(module):null_device(function)]]`()` ### Class @@ -506,7 +585,7 @@ The function provides an optional error-handling mechanism via the `state_type` ### Syntax -`call [[stdlib_system(module):delete_file(subroutine)]] (path [, err])` +`call ` [[stdlib_system(module):delete_file(subroutine)]]` (path [, err])` ### Class Subroutine @@ -532,3 +611,175 @@ The file is removed from the filesystem if the operation is successful. If the o ```fortran {!example/system/example_delete_file.f90!} ``` + +## `join_path` - Joins the provided paths according to the OS + +### Status + +Experimental + +### Description + +This interface joins the paths provided to it according to the platform specific path-separator. +i.e `\` for windows and `/` for others + +### Syntax + +`res = ` [[stdlib_system(module):join_path(interface)]] ` (p1, p2)` + +`res = ` [[stdlib_system(module):join_path(interface)]] ` (p)` + +### Class +Pure function + +### Arguments + +`p1, p2`: Shall be a character string or `type(string_type)`. It is an `intent(in)` argument. + or +`p`: Shall be a list of character strings or list of `type(string_type)`. It is an `intent(in)` argument. + +### Return values + +The resultant path, either a character string or `type(string_type)`. + +## `operator(/)` + +Alternative syntax to`join_path` using an overloaded operator. Join two paths according to the platform specific path-separator. + +### Status + +Experimental + +### Syntax + +`p = lval / rval` + +### Class + +Pure function. + +### Arguments + +`lval`: A character string or `type(string_type)`. It is an `intent(in)` argument. + +`rval`: A character string or `type(string_type)`. It is an `intent(in)` argument. + +### Result value + +The result is an `allocatable` character string or `type(string_type)` + +#### Example + +```fortran +{!example/system/example_path_join.f90!} +``` + +## `split_path` - splits a path immediately following the last separator + +### Status + +Experimental + +### Description + +This subroutine splits a path immediately following the last separator after removing the trailing separators +splitting it into most of the times a directory and a file name. + +### Syntax + +`call `[[stdlib_system(module):split_path(interface)]]`(p, head, tail)` + +### Class +Subroutine + +### Arguments + +`p`: A character string or `type(string_type)` containing the path to be split. It is an `intent(in)` argument. +`head`: The first part of the path. Either a character string or `type(string_type)`. It is an `intent(out)` argument. +`tail`: The rest part of the path. Either a character string or `type(string_type)`. It is an `intent(out)` argument. + +### Behavior + +- If `p` is empty, `head` is set to `.` and `tail` is left empty. +- If `p` consists entirely of path-separators, `head` is set to the path-separator and `tail` is left empty. +- `head` ends with a path-separator if and only if `p` appears to be a root directory or child of one. + +### Return values + +The splitted path. `head` and `tail`. + +### Example + +```fortran +{!example/system/example_path_split_path.f90!} +``` + +## `base_name` - The last part of a path + +### Status + +Experimental + +### Description + +This function returns the last part of a path after removing trailing path separators. + +### Syntax + +`res = ` [[stdlib_system(module):base_name(interface)]]`(p)` + +### Class +Function + +### Arguments + +`p`: the path, a character string or `type(string_type)`. It is an `intent(in)` argument. + +### Behavior + +- The `tail` of `[[stdlib_system(module):split_path(interface)]]` is exactly what is returned. Same Behavior. + +### Return values + +A character string or `type(string_type)`. + +### Example + +```fortran +{!example/system/example_path_base_name.f90!} +``` + +## `dir_name` - Everything except the last part of the path + +### Status + +Experimental + +### Description + +This function returns everything except the last part of a path. + +### Syntax + +`res = ` [[stdlib_system(module):dir_name(interface)]]`(p)` + +### Class +Function + +### Arguments + +`p`: the path, a character string or `type(string_type)`. It is an `intent(in)` argument. + +### Behavior + +- The `head` of `[[stdlib_system(module):split_path(interface)]]` is exactly what is returned. Same Behavior. + +### Return values + +A character string or `type(string_type)`. + +### Example + +```fortran +{!example/system/example_path_dir_name.f90!} +``` diff --git a/example/specialmatrices/CMakeLists.txt b/example/specialmatrices/CMakeLists.txt index 1f691bbde..c8b141518 100644 --- a/example/specialmatrices/CMakeLists.txt +++ b/example/specialmatrices/CMakeLists.txt @@ -1,2 +1,4 @@ ADD_EXAMPLE(specialmatrices_dp_spmv) ADD_EXAMPLE(tridiagonal_dp_type) +ADD_EXAMPLE(symtridiagonal_dp_type) +ADD_EXAMPLE(hermtridiagonal_cdp_type) diff --git a/example/specialmatrices/example_hermtridiagonal_cdp_type.f90 b/example/specialmatrices/example_hermtridiagonal_cdp_type.f90 new file mode 100644 index 000000000..e857992b7 --- /dev/null +++ b/example/specialmatrices/example_hermtridiagonal_cdp_type.f90 @@ -0,0 +1,18 @@ +program example_hermtridiagonal_matrix + use stdlib_linalg_constants, only: dp + use stdlib_specialmatrices + implicit none + + integer, parameter :: n = 5 + type(hermtridiagonal_cdp_type) :: A + complex(dp) :: dv(n), ev(n - 1) + real(dp) :: data(n, 2) + + ! Generate random tridiagonal elements. + call random_number(data); dv = cmplx(data(:, 1), 0*data(:, 2), kind=dp) + call random_number(data); ev = cmplx(data(:n - 1, 1), data(:n - 1, 2), kind=dp) + + ! Create the corresponding Hermitian tridiagonal matrix. + A = hermtridiagonal(dv, ev) + +end program example_hermtridiagonal_matrix diff --git a/example/specialmatrices/example_symtridiagonal_dp_type.f90 b/example/specialmatrices/example_symtridiagonal_dp_type.f90 new file mode 100644 index 000000000..38a2f4988 --- /dev/null +++ b/example/specialmatrices/example_symtridiagonal_dp_type.f90 @@ -0,0 +1,17 @@ +program example_symtridiagonal_matrix + use stdlib_linalg_constants, only: dp + use stdlib_specialmatrices + implicit none + + integer, parameter :: n = 5 + type(symtridiagonal_dp_type) :: A + real(dp) :: dv(n), ev(n - 1) + + ! Generate random tridiagonal elements. + call random_number(dv) + call random_number(ev) + + ! Create the corresponding Symmetric tridiagonal matrix. + A = symtridiagonal(dv, ev) + +end program example_symtridiagonal_matrix diff --git a/example/system/CMakeLists.txt b/example/system/CMakeLists.txt index a2a7525c9..57ec0c737 100644 --- a/example/system/CMakeLists.txt +++ b/example/system/CMakeLists.txt @@ -11,3 +11,9 @@ ADD_EXAMPLE(process_5) ADD_EXAMPLE(process_6) ADD_EXAMPLE(process_7) ADD_EXAMPLE(sleep) +ADD_EXAMPLE(fs_error) +ADD_EXAMPLE(path_join) +ADD_EXAMPLE(path_split_path) +ADD_EXAMPLE(path_base_name) +ADD_EXAMPLE(path_dir_name) + diff --git a/example/system/example_fs_error.f90 b/example/system/example_fs_error.f90 new file mode 100644 index 000000000..29ad3e213 --- /dev/null +++ b/example/system/example_fs_error.f90 @@ -0,0 +1,23 @@ +! Demonstrate usage of `FS_ERROR`, `FS_ERROR_CODE` +program example_fs_error + use stdlib_system, only: FS_ERROR, FS_ERROR_CODE + use stdlib_error, only: state_type, STDLIB_FS_ERROR + implicit none + + type(state_type) :: err0, err1 + + err0 = FS_ERROR("Could not create directory", "`temp.dir`", "- Already exists") + + if (err0%state == STDLIB_FS_ERROR) then + ! Error encountered: Filesystem Error: Could not create directory `temp.dir` - Already exists + print *, err0%print() + end if + + err1 = FS_ERROR_CODE(1, "Could not create directory", "`temp.dir`", "- Already exists") + + if (err1%state == STDLIB_FS_ERROR) then + ! Error encountered: Filesystem Error: code - 1, Could not create directory `temp.dir` - Already exists + print *, err1%print() + end if + +end program example_fs_error diff --git a/example/system/example_path_base_name.f90 b/example/system/example_path_base_name.f90 new file mode 100644 index 000000000..0e5f46cb1 --- /dev/null +++ b/example/system/example_path_base_name.f90 @@ -0,0 +1,16 @@ +! Usage of base_name +program example_path_base_name + use stdlib_system, only: base_name, OS_TYPE, OS_WINDOWS + character(len=:), allocatable :: p1 + + if(OS_TYPE() == OS_WINDOWS) then + p1 = 'C:\Users' + else + p1 = '/home' + endif + + print *, 'base name of '// p1 // ' -> ' // base_name(p1) + ! base name of C:\Users -> Users + ! OR + ! base name of /home -> home +end program example_path_base_name diff --git a/example/system/example_path_dir_name.f90 b/example/system/example_path_dir_name.f90 new file mode 100644 index 000000000..aff61077b --- /dev/null +++ b/example/system/example_path_dir_name.f90 @@ -0,0 +1,16 @@ +! Usage of dir_name +program example_path_dir_name + use stdlib_system, only: dir_name, OS_TYPE, OS_WINDOWS + character(len=:), allocatable :: p1, head, tail + + if(OS_TYPE() == OS_WINDOWS) then + p1 = 'C:\Users' ! C:\Users + else + p1 = '/home' ! /home + endif + + print *, 'dir_name of '// p1 // ' -> ' // dir_name(p1) + ! dir_name of C:\Users -> C:\ + ! OR + ! dir_name of /home -> / +end program example_path_dir_name diff --git a/example/system/example_path_join.f90 b/example/system/example_path_join.f90 new file mode 100644 index 000000000..cdeb24a90 --- /dev/null +++ b/example/system/example_path_join.f90 @@ -0,0 +1,23 @@ +! Usage of join_path, operator(/) +program example_path_join + use stdlib_system, only: join_path, operator(/), OS_TYPE, OS_WINDOWS + character(len=:), allocatable :: p1, p2, p3 + character(len=20) :: parr(4) + + if(OS_TYPE() == OS_WINDOWS) then + p1 = 'C:'/'Users'/'User1'/'Desktop' + p2 = join_path('C:\Users\User1', 'Desktop') + parr = [character(len=20) :: 'C:', 'Users', 'User1', 'Desktop'] + p3 = join_path(parr) + else + p1 = ''/'home'/'User1'/'Desktop' + p2 = join_path('/home/User1', 'Desktop') + parr = [character(len=20) :: '', 'home', 'User1', 'Desktop'] + p3 = join_path(parr) + end if + + ! (p1 == p2 == p3) = '/home/User1/Desktop' OR 'C:\Users\User1\Desktop' + print *, p1 ! /home/User1/Desktop OR 'C:\Users\User1\Desktop' + print *, "p1 == p2: ", p1 == p2 ! T + print *, "p2 == p3: ", p2 == p3 ! T +end program example_path_join diff --git a/example/system/example_path_split_path.f90 b/example/system/example_path_split_path.f90 new file mode 100644 index 000000000..00054d786 --- /dev/null +++ b/example/system/example_path_split_path.f90 @@ -0,0 +1,25 @@ +! Usage of split_path +program example_path_split_path + use stdlib_system, only: join_path, split_path, OS_TYPE, OS_WINDOWS + character(len=:), allocatable :: p1, head, tail + + if(OS_TYPE() == OS_WINDOWS) then + p1 = join_path('C:\Users\User1', 'Desktop') ! C:\Users\User1\Desktop + else + p1 = join_path('/home/User1', 'Desktop') ! /home/User1/Desktop + endif + + call split_path(p1, head, tail) + ! head = /home/User1 OR C:\Users\User1, tail = Desktop + print *, p1 // " -> " // head // " + " // tail + ! C:\Users\User1\Desktop -> C:\Users\User1 + Desktop + ! OR + ! /home/User1/Desktop -> /home/User1 + Desktop + + call split_path(head, p1, tail) + ! p1 = /home OR C:\Users, tail = User1 + print *, head // " -> " // p1 // " + " // tail + ! C:\Users\User1 -> C:\Users + User1 + ! OR + ! /home/User1 -> /home + User1 +end program example_path_split_path diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index e22a37e81..7d52e7fcd 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -56,7 +56,7 @@ set(fppFiles stdlib_specialfunctions_gamma.fypp stdlib_specialfunctions.fypp stdlib_specialmatrices.fypp - stdlib_specialmatrices_tridiagonal.fypp + stdlib_specialmatrices_tridiagonal.fypp stdlib_stats.fypp stdlib_stats_corr.fypp stdlib_stats_cov.fypp @@ -93,7 +93,6 @@ set(fppFiles # Preprocessed files to contain preprocessor directives -> .F90 set(cppFiles stdlib_linalg_constants.fypp - stdlib_linalg_blas.fypp stdlib_linalg_lapack.fypp ) @@ -118,6 +117,7 @@ set(SRC stdlib_sorting_radix_sort.f90 stdlib_system_subprocess.c stdlib_system_subprocess.F90 + stdlib_system_path.f90 stdlib_system.F90 stdlib_sparse.f90 stdlib_specialfunctions_legendre.f90 diff --git a/src/stdlib_specialmatrices.fypp b/src/stdlib_specialmatrices.fypp index 4db497fe1..3776b1967 100644 --- a/src/stdlib_specialmatrices.fypp +++ b/src/stdlib_specialmatrices.fypp @@ -14,7 +14,7 @@ module stdlib_specialmatrices LINALG_INTERNAL_ERROR, LINALG_VALUE_ERROR implicit none private - public :: tridiagonal + public :: tridiagonal, symtridiagonal, hermtridiagonal public :: spmv public :: dense, transpose, hermitian public :: operator(*), operator(+), operator(-) @@ -35,6 +35,26 @@ module stdlib_specialmatrices end type #:endfor + !--> Symmetric Tridiagonal matrices + #:for k1, t1, s1 in (KINDS_TYPES) + type, extends(tridiagonal_${s1}$_type), public :: symtridiagonal_${s1}$_type + !! Base type to define a `symtridiagonal` matrix. + private + #:if t1.startswith('real') + logical(lk) :: is_posdef + #:endif + end type + #:endfor + + !--> Hermitian Tridiagonal matrices + #:for k1, t1, s1 in (C_KINDS_TYPES) + type, extends(tridiagonal_${s1}$_type), public :: hermtridiagonal_${s1}$_type + !! Base type to de fine a `hermtridiagonal` matrix. + private + logical(lk) :: is_posdef + end type + #:endfor + !-------------------------------- !----- ----- !----- CONSTRUCTORS ----- @@ -126,6 +146,176 @@ module stdlib_specialmatrices #:endfor end interface + interface symtridiagonal + !! ([Specifications](../page/specs/stdlib_specialmatrices.html#SymTridiagonal)) This + !! interface provides different methods to construct a `symtridiagonal` + !! matrix. Only the non-zero elements of \( A \) are stored, i.e. + !! + !! \[ + !! A + !! = + !! \begin{bmatrix} + !! a_1 & b_1 \\ + !! b_1 & a_2 & b_2 \\ + !! & \ddots & \ddots & \ddots \\ + !! & & b_{n-2} & a_{n-1} & b_{n-1} \\ + !! & & & b_{n-1} & a_n + !! \end{bmatrix}. + !! \] + !! + !! #### Syntax + !! + !! - Construct a real `symtridiagonal` matrix from rank-1 arrays: + !! + !! ```fortran + !! integer, parameter :: n + !! real(dp), allocatable :: dv(:), ev(:) + !! type(symtridiagonal_rdp_type) :: A + !! integer :: i + !! + !! ev = [(i, i=1, n-1)]; dv = [(2*i, i=1, n)] + !! A = SymTridiagonal(dv, ev) + !! ``` + !! + !! - Construct a real `symtridiagonal` matrix with constant diagonals: + !! + !! ```fortran + !! integer, parameter :: n + !! real(dp), parameter :: a = 1.0_dp, b = 1.0_dp + !! type(symtridiagonal_rdp_type) :: A + !! + !! A = SymTridiagonal(a, b, n) + !! ``` + #:for k1, t1, s1 in (KINDS_TYPES) + pure module function initialize_symtridiagonal_pure_${s1}$(dv, ev) result(A) + !! Construct a `symtridiagonal` matrix from the rank-1 arrays + !! `dl`, `dv` and `du`. + ${t1}$, intent(in) :: dv(:), ev(:) + !! SymTridiagonal matrix elements. + type(symtridiagonal_${s1}$_type) :: A + !! Corresponding SymTridiagonal matrix. + end function + + pure module function initialize_constant_symtridiagonal_pure_${s1}$(dv, ev, n) result(A) + !! Construct a `symtridiagonal` matrix with constant elements. + ${t1}$, intent(in) :: dv, ev + !! SymTridiagonal matrix elements. + integer(ilp), intent(in) :: n + !! Matrix dimension. + type(symtridiagonal_${s1}$_type) :: A + !! Corresponding SymTridiagonal matrix. + end function + + module function initialize_symtridiagonal_impure_${s1}$(dv, ev, err) result(A) + !! Construct a `symtridiagonal` matrix from the rank-1 arrays + !! `dl`, `dv` and `du`. + ${t1}$, intent(in) :: dv(:), ev(:) + !! Tridiagonal matrix elements. + type(linalg_state_type), intent(out) :: err + !! Error handling. + type(symtridiagonal_${s1}$_type) :: A + !! Corresponding SymTridiagonal matrix. + end function + + module function initialize_constant_symtridiagonal_impure_${s1}$(dv, ev, n, err) result(A) + !! Construct a `symtridiagonal` matrix with constant elements. + ${t1}$, intent(in) :: dv, ev + !! Tridiagonal matrix elements. + integer(ilp), intent(in) :: n + !! Matrix dimension. + type(linalg_state_type), intent(out) :: err + !! Error handling. + type(symtridiagonal_${s1}$_type) :: A + !! Corresponding SymTridiagonal matrix. + end function + #:endfor + end interface + + interface hermtridiagonal + !! ([Specifications](../page/specs/stdlib_specialmatrices.html#HermTridiagonal)) This + !! interface provides different methods to construct a `hermtridiagonal` + !! matrix. Only the non-zero elements of \( A \) are stored, i.e. + !! + !! \[ + !! A + !! = + !! \begin{bmatrix} + !! a_1 & b_1 \\ + !! \bar{b}_1 & a_2 & b_2 \\ + !! & \ddots & \ddots & \ddots \\ + !! & & \bar{b}_{n-2} & a_{n-1} & b_{n-1} \\ + !! & & & \bar{b}_{n-1} & a_n + !! \end{bmatrix}. + !! \] + !! + !! #### Syntax + !! + !! - Construct a complex `hermtridiagonal` matrix from rank-1 arrays: + !! + !! ```fortran + !! integer, parameter :: n + !! complex(dp), allocatable :: dv(:), ev(:) + !! type(hermtridiagonal_cdp_type) :: A + !! integer :: i + !! + !! ev = [(i, i=1, n-1)]; dv = [(2*i, i=1, n)] + !! A = HermTridiagonal(dv, ev) + !! ``` + !! + !! - Construct a complex `hermtridiagonal` matrix with constant diagonals: + !! + !! ```fortran + !! integer, parameter :: n + !! complex(dp), parameter :: a = 1.0_dp, b = 1.0_dp + !! type(hermtridiagonal_rdp_type) :: A + !! + !! A = HermTridiagonal(a, b, n) + !! ``` + #:for k1, t1, s1 in (C_KINDS_TYPES) + pure module function initialize_hermtridiagonal_pure_${s1}$(dv, ev) result(A) + !! Construct a `hermtridiagonal` matrix from the rank-1 arrays + !! `dl`, `dv` and `du`. + ${t1}$, intent(in) :: dv(:), ev(:) + !! HermTridiagonal matrix elements. + type(hermtridiagonal_${s1}$_type) :: A + !! Corresponding HermTridiagonal matrix. + end function + + pure module function initialize_constant_hermtridiagonal_pure_${s1}$(dv, ev, n) result(A) + !! Construct a `hermtridiagonal` matrix with constant elements. + ${t1}$, intent(in) :: dv, ev + !! HermTridiagonal matrix elements. + integer(ilp), intent(in) :: n + !! Matrix dimension. + type(hermtridiagonal_${s1}$_type) :: A + !! Corresponding HermTridiagonal matrix. + end function + + module function initialize_hermtridiagonal_impure_${s1}$(dv, ev, err) result(A) + !! Construct a `hermtridiagonal` matrix from the rank-1 arrays + !! `dl`, `dv` and `du`. + ${t1}$, intent(in) :: dv(:), ev(:) + !! Tridiagonal matrix elements. + type(linalg_state_type), intent(out) :: err + !! Error handling. + type(hermtridiagonal_${s1}$_type) :: A + !! Corresponding HermTridiagonal matrix. + end function + + module function initialize_constant_hermtridiagonal_impure_${s1}$(dv, ev, n, err) result(A) + !! Construct a `hermtridiagonal` matrix with constant elements. + ${t1}$, intent(in) :: dv, ev + !! Tridiagonal matrix elements. + integer(ilp), intent(in) :: n + !! Matrix dimension. + type(linalg_state_type), intent(out) :: err + !! Error handling. + type(Hermtridiagonal_${s1}$_type) :: A + !! Corresponding HermTridiagonal matrix. + end function + #:endfor + end interface + !---------------------------------- !----- ----- !----- LINEAR ALGEBRA ----- @@ -142,12 +332,22 @@ module stdlib_specialmatrices #:for k1, t1, s1 in (KINDS_TYPES) #:for rank in RANKS module subroutine spmv_tridiag_${rank}$d_${s1}$(A, x, y, alpha, beta, op) - type(tridiagonal_${s1}$_type), intent(in) :: A + !! Matrix-vector kernel for matrices extended from the `tridiagonal` class. + class(tridiagonal_${s1}$_type), intent(in) :: A + !! Input matrix. ${t1}$, intent(in), contiguous, target :: x${ranksuffix(rank)}$ + !! Vector(s) to be multiplied by the matrix `A`. ${t1}$, intent(inout), contiguous, target :: y${ranksuffix(rank)}$ + !! Resulting vector. real(${k1}$), intent(in), optional :: alpha + !! Real scaling parameter for `Ax` real(${k1}$), intent(in), optional :: beta + !! Real scaling parameter for `y` (right-hand side) character(1), intent(in), optional :: op + !! Which operation to perform: + !! - `op = "N"` : `y = alpha * A @ x + beta * y` + !! - `op = "T"` : `y = alpha * A.T @ x + beta * y` + !! - `op = "H"` : `y = alpha * A.H @ x + beta * y` (for complex-valued matrices only). end subroutine #:endfor #:endfor @@ -166,7 +366,7 @@ module stdlib_specialmatrices #:for k1, t1, s1 in (KINDS_TYPES) pure module function tridiagonal_to_dense_${s1}$(A) result(B) !! Convert a `tridiagonal` matrix to its dense representation. - type(tridiagonal_${s1}$_type), intent(in) :: A + class(tridiagonal_${s1}$_type), intent(in) :: A !! Input Tridiagonal matrix. ${t1}$, allocatable :: B(:, :) !! Corresponding dense matrix. @@ -184,6 +384,18 @@ module stdlib_specialmatrices !! Input matrix. type(tridiagonal_${s1}$_type) :: B end function + pure module function transpose_symtridiagonal_${s1}$(A) result(B) + type(symtridiagonal_${s1}$_type), intent(in) :: A + !! Input matrix. + type(symtridiagonal_${s1}$_type) :: B + end function + #:if t1.startswith('complex') + pure module function transpose_hermtridiagonal_${s1}$(A) result(B) + type(hermtridiagonal_${s1}$_type), intent(in) :: A + !! Input matrix. + type(hermtridiagonal_${s1}$_type) :: B + end function + #:endif #:endfor end interface @@ -198,6 +410,18 @@ module stdlib_specialmatrices !! Input matrix. type(tridiagonal_${s1}$_type) :: B end function + pure module function hermitian_symtridiagonal_${s1}$(A) result(B) + type(symtridiagonal_${s1}$_type), intent(in) :: A + !! Input matrix. + type(symtridiagonal_${s1}$_type) :: B + end function + #:if t1.startswith('complex') + pure module function hermitian_hermtridiagonal_${s1}$(A) result(B) + type(hermtridiagonal_${s1}$_type), intent(in) :: A + !! Input matrix. + type(hermtridiagonal_${s1}$_type) :: B + end function + #:endif #:endfor end interface @@ -222,32 +446,146 @@ module stdlib_specialmatrices ${t1}$, intent(in) :: alpha type(tridiagonal_${s1}$_type) :: B end function + + pure module function scalar_multiplication_symtridiagonal_${s1}$(alpha, A) result(B) + ${t1}$, intent(in) :: alpha + type(symtridiagonal_${s1}$_type), intent(in) :: A + type(symtridiagonal_${s1}$_type) :: B + end function + pure module function scalar_multiplication_bis_symtridiagonal_${s1}$(A, alpha) result(B) + type(symtridiagonal_${s1}$_type), intent(in) :: A + ${t1}$, intent(in) :: alpha + type(symtridiagonal_${s1}$_type) :: B + end function + + #:if t1.startswith("complex") + pure module function scalar_multiplication_hermtridiagonal_${s1}$(alpha, A) result(B) + ${t1}$, intent(in) :: alpha + type(hermtridiagonal_${s1}$_type), intent(in) :: A + type(tridiagonal_${s1}$_type) :: B + end function + pure module function scalar_multiplication_bis_hermtridiagonal_${s1}$(A, alpha) result(B) + type(hermtridiagonal_${s1}$_type), intent(in) :: A + ${t1}$, intent(in) :: alpha + type(tridiagonal_${s1}$_type) :: B + end function + pure module function real_scalar_multiplication_hermtridiagonal_${s1}$(alpha, A) result(B) + real(${k1}$), intent(in) :: alpha + type(hermtridiagonal_${s1}$_type), intent(in) :: A + type(hermtridiagonal_${s1}$_type) :: B + end function + pure module function real_scalar_multiplication_bis_hermtridiagonal_${s1}$(A, alpha) result(B) + type(hermtridiagonal_${s1}$_type), intent(in) :: A + real(${k1}$), intent(in) :: alpha + type(hermtridiagonal_${s1}$_type) :: B + end function + #:endif #:endfor end interface interface operator(+) !! Overload the `+` operator for matrix-matrix addition. The two matrices need to - !! be of the same type and kind. + !! be of the class type and kind. !! [Specifications](../page/specs/stdlib_specialmatrices.html#operators) #:for k1, t1, s1 in (KINDS_TYPES) - pure module function matrix_add_tridiagonal_${s1}$(A, B) result(C) + pure module function matrix_add_tridiag_tridiag_${s1}$(A, B) result(C) + type(tridiagonal_${s1}$_type), intent(in) :: A, B + type(tridiagonal_${s1}$_type) :: C + end function + pure module function matrix_add_tridiag_symtridiag_${s1}$(A, B) result(C) + type(tridiagonal_${s1}$_type), intent(in) :: A + type(symtridiagonal_${s1}$_type), intent(in) :: B + type(tridiagonal_${s1}$_type) :: C + end function + pure module function matrix_add_symtridiag_tridiag_${s1}$(A, B) result(C) + type(symtridiagonal_${s1}$_type), intent(in) :: A + type(tridiagonal_${s1}$_type), intent(in) :: B + type(tridiagonal_${s1}$_type) :: C + end function + pure module function matrix_add_symtridiag_symtridiag_${s1}$(A, B) result(C) + type(symtridiagonal_${s1}$_type), intent(in) :: A, B + type(symtridiagonal_${s1}$_type) :: C + end function + + #:if t1.startswith("complex") + pure module function matrix_add_tridiag_hermtridiag_${s1}$(A, B) result(C) type(tridiagonal_${s1}$_type), intent(in) :: A + type(hermtridiagonal_${s1}$_type), intent(in) :: B + type(tridiagonal_${s1}$_type) :: C + end function + pure module function matrix_add_hermtridiag_tridiag_${s1}$(A, B) result(C) + type(hermtridiagonal_${s1}$_type), intent(in) :: A type(tridiagonal_${s1}$_type), intent(in) :: B type(tridiagonal_${s1}$_type) :: C end function + pure module function matrix_add_symtridiag_hermtridiag_${s1}$(A, B) result(C) + type(symtridiagonal_${s1}$_type), intent(in) :: A + type(hermtridiagonal_${s1}$_type), intent(in) :: B + type(tridiagonal_${s1}$_type) :: C + end function + pure module function matrix_add_hermtridiag_symtridiag_${s1}$(A, B) result(C) + type(hermtridiagonal_${s1}$_type), intent(in) :: A + type(symtridiagonal_${s1}$_type), intent(in) :: B + type(tridiagonal_${s1}$_type) :: C + end function + pure module function matrix_add_hermtridiag_hermtridiag_${s1}$(A, B) result(C) + type(hermtridiagonal_${s1}$_type), intent(in) :: A, B + type(hermtridiagonal_${s1}$_type) :: C + end function + #:endif #:endfor end interface interface operator(-) !! Overload the `-` operator for matrix-matrix subtraction. The two matrices need to - !! be of the same type and kind. + !! be of the same class and kind. !! [Specifications](../page/specs/stdlib_specialmatrices.html#operators) #:for k1, t1, s1 in (KINDS_TYPES) - pure module function matrix_sub_tridiagonal_${s1}$(A, B) result(C) + pure module function matrix_sub_tridiag_tridiag_${s1}$(A, B) result(C) + type(tridiagonal_${s1}$_type), intent(in) :: A, B + type(tridiagonal_${s1}$_type) :: C + end function + pure module function matrix_sub_tridiag_symtridiag_${s1}$(A, B) result(C) type(tridiagonal_${s1}$_type), intent(in) :: A + type(symtridiagonal_${s1}$_type), intent(in) :: B + type(tridiagonal_${s1}$_type) :: C + end function + pure module function matrix_sub_symtridiag_tridiag_${s1}$(A, B) result(C) + type(symtridiagonal_${s1}$_type), intent(in) :: A type(tridiagonal_${s1}$_type), intent(in) :: B type(tridiagonal_${s1}$_type) :: C end function + pure module function matrix_sub_symtridiag_symtridiag_${s1}$(A, B) result(C) + type(symtridiagonal_${s1}$_type), intent(in) :: A, B + type(symtridiagonal_${s1}$_type) :: C + end function + + #:if t1.startswith("complex") + pure module function matrix_sub_tridiag_hermtridiag_${s1}$(A, B) result(C) + type(tridiagonal_${s1}$_type), intent(in) :: A + type(hermtridiagonal_${s1}$_type), intent(in) :: B + type(tridiagonal_${s1}$_type) :: C + end function + pure module function matrix_sub_hermtridiag_tridiag_${s1}$(A, B) result(C) + type(hermtridiagonal_${s1}$_type), intent(in) :: A + type(tridiagonal_${s1}$_type), intent(in) :: B + type(tridiagonal_${s1}$_type) :: C + end function + pure module function matrix_sub_symtridiag_hermtridiag_${s1}$(A, B) result(C) + type(symtridiagonal_${s1}$_type), intent(in) :: A + type(hermtridiagonal_${s1}$_type), intent(in) :: B + type(tridiagonal_${s1}$_type) :: C + end function + pure module function matrix_sub_hermtridiag_symtridiag_${s1}$(A, B) result(C) + type(hermtridiagonal_${s1}$_type), intent(in) :: A + type(symtridiagonal_${s1}$_type), intent(in) :: B + type(tridiagonal_${s1}$_type) :: C + end function + pure module function matrix_sub_hermtridiag_hermtridiag_${s1}$(A, B) result(C) + type(hermtridiagonal_${s1}$_type), intent(in) :: A, B + type(hermtridiagonal_${s1}$_type) :: C + end function + #:endif #:endfor end interface diff --git a/src/stdlib_specialmatrices_tridiagonal.fypp b/src/stdlib_specialmatrices_tridiagonal.fypp index 8fabaa5df..9bf4b7271 100644 --- a/src/stdlib_specialmatrices_tridiagonal.fypp +++ b/src/stdlib_specialmatrices_tridiagonal.fypp @@ -7,6 +7,7 @@ submodule (stdlib_specialmatrices) tridiagonal_matrices use stdlib_linalg_lapack, only: lagtm character(len=*), parameter :: this = "tridiagonal matrices" + contains !-------------------------------- @@ -15,6 +16,8 @@ submodule (stdlib_specialmatrices) tridiagonal_matrices !----- ----- !-------------------------------- + ! ----- Tridiagonal matrices ----- + #:for k1, t1, s1 in (KINDS_TYPES) pure module function initialize_tridiagonal_pure_${s1}$(dl, dv, du) result(A) !! Construct a `tridiagonal` matrix from the rank-1 arrays @@ -137,6 +140,238 @@ submodule (stdlib_specialmatrices) tridiagonal_matrices end function #:endfor + !----- Symmetric Tridiagonal matrices ----- + + #:for k1, t1, s1 in (KINDS_TYPES) + pure module function initialize_symtridiagonal_pure_${s1}$(dv, ev) result(A) + !! Construct a `symtridiagonal` matrix from the rank-1 arrays + !! `dv` and `ev`. + ${t1}$, intent(in) :: dv(:), ev(:) + !! symtridiagonal matrix elements. + type(symtridiagonal_${s1}$_type) :: A + !! Corresponding symtridiagonal matrix. + + ! Internal variables. + integer(ilp) :: n + type(linalg_state_type) :: err0 + + ! Sanity check. + n = size(dv, kind=ilp) + if (n <= 0) then + err0 = linalg_state_type(this, LINALG_VALUE_ERROR, "Matrix size needs to be positive, n = ", n, ".") + call linalg_error_handling(err0) + endif + if (size(ev, kind=ilp) /= n-1) then + err0 = linalg_state_type(this, LINALG_VALUE_ERROR, "Vector ev does not have the correct length.") + call linalg_error_handling(err0) + endif + + ! Description of the matrix. + A%n = n + ! Matrix elements. + A%dl = ev ; A%dv = dv ; A%du = ev + end function + + pure module function initialize_constant_symtridiagonal_pure_${s1}$(dv, ev, n) result(A) + !! Construct a `symtridiagonal` matrix with constant elements. + ${t1}$, intent(in) :: dv, ev + !! symtridiagonal matrix elements. + integer(ilp), intent(in) :: n + !! Matrix dimension. + type(symtridiagonal_${s1}$_type) :: A + !! Corresponding symtridiagonal matrix. + + ! Internal variables. + integer(ilp) :: i + type(linalg_state_type) :: err0 + + ! Description of the matrix. + A%n = n + if (n <= 0) then + err0 = linalg_state_type(this, LINALG_VALUE_ERROR, "Matrix size needs to be positive, n = ", n, ".") + call linalg_error_handling(err0) + endif + ! Matrix elements. + A%dl = [(ev, i = 1, n-1)] + A%dv = [(dv, i = 1, n-1)] + A%du = [(ev, i = 1, n-1)] + end function + + module function initialize_symtridiagonal_impure_${s1}$(dv, ev, err) result(A) + !! Construct a `symtridiagonal` matrix from the rank-1 arrays + !! `dl` and `ev`. + ${t1}$, intent(in) :: dv(:), ev(:) + !! symtridiagonal matrix elements. + type(linalg_state_type), intent(out) :: err + !! Error handling. + type(symtridiagonal_${s1}$_type) :: A + !! Corresponding symtridiagonal matrix. + + ! Internal variables. + integer(ilp) :: n + type(linalg_state_type) :: err0 + + ! Sanity check. + n = size(dv, kind=ilp) + if (n <= 0) then + err0 = linalg_state_type(this, LINALG_VALUE_ERROR, "Matrix size needs to be positive, n = ", n, ".") + call linalg_error_handling(err0, err) + endif + if (size(ev, kind=ilp) /= n-1) then + err0 = linalg_state_type(this, LINALG_VALUE_ERROR, "Vector ev does not have the correct length.") + call linalg_error_handling(err0, err) + endif + + ! Description of the matrix. + A%n = n + ! Matrix elements. + A%dl = ev ; A%dv = dv ; A%du = ev + end function + + module function initialize_constant_symtridiagonal_impure_${s1}$(dv, ev, n, err) result(A) + !! Construct a `symtridiagonal` matrix with constant elements. + ${t1}$, intent(in) :: dv, ev + !! symtridiagonal matrix elements. + integer(ilp), intent(in) :: n + !! Matrix dimension. + type(linalg_state_type), intent(out) :: err + !! Error handling + type(symtridiagonal_${s1}$_type) :: A + !! Corresponding symtridiagonal matrix. + + ! Internal variables. + integer(ilp) :: i + type(linalg_state_type) :: err0 + + ! Description of the matrix. + A%n = n + if (n <= 0) then + err0 = linalg_state_type(this, LINALG_VALUE_ERROR, "Matrix size needs to be positive, n = ", n, ".") + call linalg_error_handling(err0, err) + endif + ! Matrix elements. + A%dl = [(ev, i = 1, n)] + A%dv = [(dv, i = 1, n-1)] + A%du = [(ev, i = 1, n)] + end function + #:endfor + + !----- Hermitian Tridiagonal matrices ----- + + #:for k1, t1, s1 in (C_KINDS_TYPES) + pure module function initialize_hermtridiagonal_pure_${s1}$(dv, ev) result(A) + !! Construct a `hermtridiagonal` matrix from the rank-1 arrays + !! `dl` and `ev`. + ${t1}$, intent(in) :: dv(:), ev(:) + !! hermtridiagonal matrix elements. + type(hermtridiagonal_${s1}$_type) :: A + !! Corresponding hermtridiagonal matrix. + + ! Internal variables. + integer(ilp) :: n + type(linalg_state_type) :: err0 + + ! Sanity check. + n = size(dv, kind=ilp) + if (n <= 0) then + err0 = linalg_state_type(this, LINALG_VALUE_ERROR, "Matrix size needs to be positive, n = ", n, ".") + call linalg_error_handling(err0) + endif + if (size(ev, kind=ilp) /= n-1) then + err0 = linalg_state_type(this, LINALG_VALUE_ERROR, "Vector ev does not have the correct length.") + call linalg_error_handling(err0) + endif + + ! Description of the matrix. + A%n = n + ! Matrix elements. + A%dl = conjg(ev) ; A%dv = dv%re ; A%du = ev + end function + + pure module function initialize_constant_hermtridiagonal_pure_${s1}$(dv, ev, n) result(A) + !! Construct a `hermtridiagonal` matrix with constant elements. + ${t1}$, intent(in) :: dv, ev + !! hermtridiagonal matrix elements. + integer(ilp), intent(in) :: n + !! Matrix dimension. + type(hermtridiagonal_${s1}$_type) :: A + !! Corresponding hermtridiagonal matrix. + + ! Internal variables. + integer(ilp) :: i + type(linalg_state_type) :: err0 + + ! Description of the matrix. + A%n = n + if (n <= 0) then + err0 = linalg_state_type(this, LINALG_VALUE_ERROR, "Matrix size needs to be positive, n = ", n, ".") + call linalg_error_handling(err0) + endif + ! Matrix elements. + A%dl = [(conjg(ev), i = 1, n-1)] + A%dv = [(dv%re, i = 1, n-1)] + A%du = [(ev, i = 1, n-1)] + end function + + module function initialize_hermtridiagonal_impure_${s1}$(dv, ev, err) result(A) + !! Construct a `hermtridiagonal` matrix from the rank-1 arrays + !! `dl`, `dv` and `du`. + ${t1}$, intent(in) :: dv(:), ev(:) + !! symtridiagonal matrix elements. + type(linalg_state_type), intent(out) :: err + !! Error handling. + type(hermtridiagonal_${s1}$_type) :: A + !! Corresponding hermtridiagonal matrix. + + ! Internal variables. + integer(ilp) :: n + type(linalg_state_type) :: err0 + + ! Sanity check. + n = size(dv, kind=ilp) + if (n <= 0) then + err0 = linalg_state_type(this, LINALG_VALUE_ERROR, "Matrix size needs to be positive, n = ", n, ".") + call linalg_error_handling(err0, err) + endif + if (size(ev, kind=ilp) /= n-1) then + err0 = linalg_state_type(this, LINALG_VALUE_ERROR, "Vector ev does not have the correct length.") + call linalg_error_handling(err0, err) + endif + + ! Description of the matrix. + A%n = n + ! Matrix elements. + A%dl = conjg(ev) ; A%dv = dv%re ; A%du = ev + end function + + module function initialize_constant_hermtridiagonal_impure_${s1}$(dv, ev, n, err) result(A) + !! Construct a `hermtridiagonal` matrix with constant elements. + ${t1}$, intent(in) :: dv, ev + !! symtridiagonal matrix elements. + integer(ilp), intent(in) :: n + !! Matrix dimension. + type(linalg_state_type), intent(out) :: err + !! Error handling + type(hermtridiagonal_${s1}$_type) :: A + !! Corresponding hermtridiagonal matrix. + + ! Internal variables. + integer(ilp) :: i + type(linalg_state_type) :: err0 + + ! Description of the matrix. + A%n = n + if (n <= 0) then + err0 = linalg_state_type(this, LINALG_VALUE_ERROR, "Matrix size needs to be positive, n = ", n, ".") + call linalg_error_handling(err0, err) + endif + ! Matrix elements. + A%dl = [(conjg(ev), i = 1, n)] + A%dv = [(dv%re, i = 1, n-1)] + A%du = [(ev, i = 1, n)] + end function + #:endfor + !----------------------------------------- !----- ----- !----- MATRIX-VECTOR PRODUCT ----- @@ -147,7 +382,7 @@ submodule (stdlib_specialmatrices) tridiagonal_matrices #:for k1, t1, s1 in (KINDS_TYPES) #:for rank in RANKS module subroutine spmv_tridiag_${rank}$d_${s1}$(A, x, y, alpha, beta, op) - type(tridiagonal_${s1}$_type), intent(in) :: A + class(tridiagonal_${s1}$_type), intent(in) :: A ${t1}$, intent(in), contiguous, target :: x${ranksuffix(rank)}$ ${t1}$, intent(inout), contiguous, target :: y${ranksuffix(rank)}$ real(${k1}$), intent(in), optional :: alpha @@ -166,6 +401,7 @@ submodule (stdlib_specialmatrices) tridiagonal_matrices alpha_ = 1.0_${k1}$ ; if (present(alpha)) alpha_ = alpha beta_ = 0.0_${k1}$ ; if (present(beta)) beta_ = beta op_ = "N" ; if (present(op)) op_ = op + if (op_ == "H") op_ = "C" ! Prepare Lapack arguments. n = A%n ; ldx = n ; ldy = n ; y = 0.0_${k1}$ @@ -175,6 +411,7 @@ submodule (stdlib_specialmatrices) tridiagonal_matrices ! Pointer trick. xmat(1:n, 1:nrhs) => x ; ymat(1:n, 1:nrhs) => y call lagtm(op_, n, nrhs, alpha_, A%dl, A%dv, A%du, xmat, ldx, beta_, ymat, ldy) + nullify(xmat, ymat) #:else call lagtm(op_, n, nrhs, alpha_, A%dl, A%dv, A%du, x, ldx, beta_, y, ldy) #:endif @@ -188,10 +425,12 @@ submodule (stdlib_specialmatrices) tridiagonal_matrices !----- ----- !------------------------------------- + !----- To dense matrix ----- + #:for k1, t1, s1 in (KINDS_TYPES) pure module function tridiagonal_to_dense_${s1}$(A) result(B) !! Convert a `tridiagonal` matrix to its dense representation. - type(tridiagonal_${s1}$_type), intent(in) :: A + class(tridiagonal_${s1}$_type), intent(in) :: A !! Input tridiagonal matrix. ${t1}$, allocatable :: B(:, :) !! Corresponding dense matrix. @@ -216,19 +455,35 @@ submodule (stdlib_specialmatrices) tridiagonal_matrices end function #:endfor + !----- Matrix transposition ----- + #:for k1, t1, s1 in (KINDS_TYPES) pure module function transpose_tridiagonal_${s1}$(A) result(B) type(tridiagonal_${s1}$_type), intent(in) :: A - !! Input matrix. type(tridiagonal_${s1}$_type) :: B B = tridiagonal(A%du, A%dv, A%dl) end function + + pure module function transpose_symtridiagonal_${s1}$(A) result(B) + type(symtridiagonal_${s1}$_type), intent(in) :: A + type(symtridiagonal_${s1}$_type) :: B + B = symtridiagonal(A%dv, A%du) + end function + + #:if t1.startswith('complex') + pure module function transpose_hermtridiagonal_${s1}$(A) result(B) + type(hermtridiagonal_${s1}$_type), intent(in) :: A + type(hermtridiagonal_${s1}$_type) :: B + B = hermtridiagonal(A%dv, A%dl) + end function + #:endif #:endfor + !----- Hermitian operator ----- + #:for k1, t1, s1 in (KINDS_TYPES) pure module function hermitian_tridiagonal_${s1}$(A) result(B) type(tridiagonal_${s1}$_type), intent(in) :: A - !! Input matrix. type(tridiagonal_${s1}$_type) :: B #:if t1.startswith("complex") B = tridiagonal(conjg(A%du), conjg(A%dv), conjg(A%dl)) @@ -236,8 +491,28 @@ submodule (stdlib_specialmatrices) tridiagonal_matrices B = tridiagonal(A%du, A%dv, A%dl) #:endif end function + + pure module function hermitian_symtridiagonal_${s1}$(A) result(B) + type(symtridiagonal_${s1}$_type), intent(in) :: A + type(symtridiagonal_${s1}$_type) :: B + #:if t1.startswith("complex") + B = symtridiagonal(conjg(A%dv), conjg(A%du)) + #:else + B = A + #:endif + end function + + #:if t1.startswith("complex") + pure module function hermitian_hermtridiagonal_${s1}$(A) result(B) + type(hermtridiagonal_${s1}$_type), intent(in) :: A + type(hermtridiagonal_${s1}$_type) :: B + B = A + end function + #:endif #:endfor + !----- Scalar multiplication ----- + #:for k1, t1, s1 in (KINDS_TYPES) pure module function scalar_multiplication_tridiagonal_${s1}$(alpha, A) result(B) ${t1}$, intent(in) :: alpha @@ -251,27 +526,226 @@ submodule (stdlib_specialmatrices) tridiagonal_matrices type(tridiagonal_${s1}$_type), intent(in) :: A ${t1}$, intent(in) :: alpha type(tridiagonal_${s1}$_type) :: B + B = scalar_multiplication_tridiagonal_${s1}$(alpha, A) + end function + + pure module function scalar_multiplication_symtridiagonal_${s1}$(alpha, A) result(B) + ${t1}$, intent(in) :: alpha + type(symtridiagonal_${s1}$_type), intent(in) :: A + type(symtridiagonal_${s1}$_type) :: B + B = symtridiagonal(A%dv, A%du) + B%dl = alpha*B%dl; B%dv = alpha*B%dv; B%du = B%dl + end function + + pure module function scalar_multiplication_bis_symtridiagonal_${s1}$(A, alpha) result(B) + type(symtridiagonal_${s1}$_type), intent(in) :: A + ${t1}$, intent(in) :: alpha + type(symtridiagonal_${s1}$_type) :: B + B = symtridiagonal(A%dv, A%du) + B%dl = alpha*B%dl; B%dv = alpha*B%dv; B%du = B%dl + end function + + #:if t1.startswith("complex") + pure module function scalar_multiplication_hermtridiagonal_${s1}$(alpha, A) result(B) + ${t1}$, intent(in) :: alpha + type(hermtridiagonal_${s1}$_type), intent(in) :: A + type(tridiagonal_${s1}$_type) :: B B = tridiagonal(A%dl, A%dv, A%du) B%dl = alpha*B%dl; B%dv = alpha*B%dv; B%du = alpha*B%du end function + + pure module function scalar_multiplication_bis_hermtridiagonal_${s1}$(A, alpha) result(B) + type(hermtridiagonal_${s1}$_type), intent(in) :: A + ${t1}$, intent(in) :: alpha + type(tridiagonal_${s1}$_type) :: B + B = tridiagonal(A%dl, A%dv, A%du) + B%dl = alpha*B%dl; B%dv = alpha*B%dv; B%du = alpha*B%du + end function + + pure module function real_scalar_multiplication_hermtridiagonal_${s1}$(alpha, A) result(B) + real(${k1}$), intent(in) :: alpha + type(hermtridiagonal_${s1}$_type), intent(in) :: A + type(hermtridiagonal_${s1}$_type) :: B + B = hermtridiagonal(A%dv, A%du) + B%dl = alpha*B%dl; B%dv = alpha*B%dv; B%du = conjg(B%dl) + end function + + pure module function real_scalar_multiplication_bis_hermtridiagonal_${s1}$(A, alpha) result(B) + type(hermtridiagonal_${s1}$_type), intent(in) :: A + real(${k1}$), intent(in) :: alpha + type(hermtridiagonal_${s1}$_type) :: B + B = hermtridiagonal(A%dv, A%du) + B%dl = alpha*B%dl; B%dv = alpha*B%dv; B%du = conjg(B%dl) + end function + #:endif + #:endfor + + !----- Matrix addition ----- + + #:for k1, t1, s1 in (KINDS_TYPES) + ! Tridiag + Tridiag = Tridiag + pure module function matrix_add_tridiag_tridiag_${s1}$(A, B) result(C) + type(tridiagonal_${s1}$_type), intent(in) :: A, B + type(tridiagonal_${s1}$_type) :: C + C = tridiagonal(A%dl, A%dv, A%du) + C%dl = C%dl + B%dl ; C%dv = C%dv + B%dv ; C%du = C%du + B%du + end function + + ! Tridiag + SymTridiag = Tridiag + pure module function matrix_add_tridiag_symtridiag_${s1}$(A, B) result(C) + type(tridiagonal_${s1}$_type), intent(in) :: A + type(symtridiagonal_${s1}$_type), intent(in) :: B + type(tridiagonal_${s1}$_type) :: C + C = tridiagonal(A%dl, A%dv, A%du) + C%dl = C%dl + B%dl ; C%dv = C%dv + B%dv ; C%du = C%du + B%du + end function + + ! SymTridiag + Tridiag = Tridiag + pure module function matrix_add_symtridiag_tridiag_${s1}$(A, B) result(C) + type(symtridiagonal_${s1}$_type), intent(in) :: A + type(tridiagonal_${s1}$_type), intent(in) :: B + type(tridiagonal_${s1}$_type) :: C + C = tridiagonal(A%dl, A%dv, A%du) + C%dl = C%dl + B%dl ; C%dv = C%dv + B%dv ; C%du = C%du + B%du + end function + + ! SymTridiag + SymTridiag = SymTridiag + pure module function matrix_add_symtridiag_symtridiag_${s1}$(A, B) result(C) + type(symtridiagonal_${s1}$_type), intent(in) :: A, B + type(symtridiagonal_${s1}$_type) :: C + C = symtridiagonal(A%dv, A%du) + C%dl = C%dl + B%dl ; C%dv = C%dv + B%dv ; C%du = C%dl + end function + + #:if t1.startswith("complex") + ! Tridiag + HermTridiag = Tridiag + pure module function matrix_add_tridiag_hermtridiag_${s1}$(A, B) result(C) + type(tridiagonal_${s1}$_type), intent(in) :: A + type(hermtridiagonal_${s1}$_type), intent(in) :: B + type(tridiagonal_${s1}$_type) :: C + C = tridiagonal(A%dl, A%dv, A%du) + C%dl = C%dl + B%dl ; C%dv = C%dv + B%dv ; C%du = C%du + B%du + end function + + ! HermTridiag + Tridiag = Tridiag + pure module function matrix_add_hermtridiag_tridiag_${s1}$(A, B) result(C) + type(hermtridiagonal_${s1}$_type), intent(in) :: A + type(tridiagonal_${s1}$_type), intent(in) :: B + type(tridiagonal_${s1}$_type) :: C + C = tridiagonal(A%dl, A%dv, A%du) + C%dl = C%dl + B%dl ; C%dv = C%dv + B%dv ; C%du = C%du + B%du + end function + + ! SymTridiag + HermTridiag = Tridiag + pure module function matrix_add_symtridiag_hermtridiag_${s1}$(A, B) result(C) + type(symtridiagonal_${s1}$_type), intent(in) :: A + type(hermtridiagonal_${s1}$_type), intent(in) :: B + type(tridiagonal_${s1}$_type) :: C + C = tridiagonal(A%dl, A%dv, A%du) + C%dl = C%dl + B%dl ; C%dv = C%dv + B%dv ; C%du = C%du + B%du + end function + + ! HermTridiag + SymTridiag = Tridiag + pure module function matrix_add_hermtridiag_symtridiag_${s1}$(A, B) result(C) + type(hermtridiagonal_${s1}$_type), intent(in) :: A + type(symtridiagonal_${s1}$_type), intent(in) :: B + type(tridiagonal_${s1}$_type) :: C + C = tridiagonal(A%dl, A%dv, A%du) + C%dl = C%dl + B%dl ; C%dv = C%dv + B%dv ; C%du = C%du + B%du + end function + + ! HermTridiag + HermTridiag = HermTridiag + pure module function matrix_add_hermtridiag_hermtridiag_${s1}$(A, B) result(C) + type(hermtridiagonal_${s1}$_type), intent(in) :: A, B + type(hermtridiagonal_${s1}$_type) :: C + C = hermtridiagonal(A%dv, A%du) + C%dl = C%dl + B%dl ; C%dv = C%dv + B%dv ; C%du = conjg(C%dl) + end function + #:endif #:endfor + !----- Matrix subtraction ----- + #:for k1, t1, s1 in (KINDS_TYPES) - pure module function matrix_add_tridiagonal_${s1}$(A, B) result(C) + ! Tridiag - Tridiag = Tridiag + pure module function matrix_sub_tridiag_tridiag_${s1}$(A, B) result(C) + type(tridiagonal_${s1}$_type), intent(in) :: A, B + type(tridiagonal_${s1}$_type) :: C + C = tridiagonal(A%dl, A%dv, A%du) + C%dl = C%dl - B%dl ; C%dv = C%dv - B%dv ; C%du = C%du - B%du + end function + + ! Tridiag - SymTridiag = Tridiag + pure module function matrix_sub_tridiag_symtridiag_${s1}$(A, B) result(C) type(tridiagonal_${s1}$_type), intent(in) :: A + type(symtridiagonal_${s1}$_type), intent(in) :: B + type(tridiagonal_${s1}$_type) :: C + C = tridiagonal(A%dl, A%dv, A%du) + C%dl = C%dl - B%dl ; C%dv = C%dv - B%dv ; C%du = C%du - B%du + end function + + ! SymTridiag - Tridiag = Tridiag + pure module function matrix_sub_symtridiag_tridiag_${s1}$(A, B) result(C) + type(symtridiagonal_${s1}$_type), intent(in) :: A type(tridiagonal_${s1}$_type), intent(in) :: B type(tridiagonal_${s1}$_type) :: C C = tridiagonal(A%dl, A%dv, A%du) - C%dl = C%dl + B%dl; C%dv = C%dv + B%dv; C%du = C%du + B%du + C%dl = C%dl - B%dl ; C%dv = C%dv - B%dv ; C%du = C%du - B%du end function - pure module function matrix_sub_tridiagonal_${s1}$(A, B) result(C) + ! SymTridiag - SymTridiag = SymTridiag + pure module function matrix_sub_symtridiag_symtridiag_${s1}$(A, B) result(C) + type(symtridiagonal_${s1}$_type), intent(in) :: A, B + type(symtridiagonal_${s1}$_type) :: C + C = symtridiagonal(A%dv, A%du) + C%dl = C%dl - B%dl ; C%dv = C%dv - B%dv ; C%du = C%dl + end function + + #:if t1.startswith("complex") + ! Tridiag - HermTridiag = Tridiag + pure module function matrix_sub_tridiag_hermtridiag_${s1}$(A, B) result(C) type(tridiagonal_${s1}$_type), intent(in) :: A + type(hermtridiagonal_${s1}$_type), intent(in) :: B + type(tridiagonal_${s1}$_type) :: C + C = tridiagonal(A%dl, A%dv, A%du) + C%dl = C%dl - B%dl ; C%dv = C%dv - B%dv ; C%du = C%du - B%du + end function + + ! HermTridiag - Tridiag = Tridiag + pure module function matrix_sub_hermtridiag_tridiag_${s1}$(A, B) result(C) + type(hermtridiagonal_${s1}$_type), intent(in) :: A type(tridiagonal_${s1}$_type), intent(in) :: B type(tridiagonal_${s1}$_type) :: C C = tridiagonal(A%dl, A%dv, A%du) - C%dl = C%dl - B%dl; C%dv = C%dv - B%dv; C%du = C%du - B%du + C%dl = C%dl - B%dl ; C%dv = C%dv - B%dv ; C%du = C%du - B%du + end function + + ! SymTridiag - HermTridiag = Tridiag + pure module function matrix_sub_symtridiag_hermtridiag_${s1}$(A, B) result(C) + type(symtridiagonal_${s1}$_type), intent(in) :: A + type(hermtridiagonal_${s1}$_type), intent(in) :: B + type(tridiagonal_${s1}$_type) :: C + C = tridiagonal(A%dl, A%dv, A%du) + C%dl = C%dl - B%dl ; C%dv = C%dv - B%dv ; C%du = C%du - B%du + end function + + ! HermTridiag - SymTridiag = Tridiag + pure module function matrix_sub_hermtridiag_symtridiag_${s1}$(A, B) result(C) + type(hermtridiagonal_${s1}$_type), intent(in) :: A + type(symtridiagonal_${s1}$_type), intent(in) :: B + type(tridiagonal_${s1}$_type) :: C + C = tridiagonal(A%dl, A%dv, A%du) + C%dl = C%dl - B%dl ; C%dv = C%dv - B%dv ; C%du = C%du - B%du + end function + + ! HermTridiag - HermTridiag = HermTridiag + pure module function matrix_sub_hermtridiag_hermtridiag_${s1}$(A, B) result(C) + type(hermtridiagonal_${s1}$_type), intent(in) :: A, B + type(hermtridiagonal_${s1}$_type) :: C + C = hermtridiagonal(A%dv, A%du) + C%dl = C%dl - B%dl ; C%dv = C%dv - B%dv ; C%du = conjg(C%dl) end function + #:endif #:endfor end submodule diff --git a/src/stdlib_system.F90 b/src/stdlib_system.F90 index a9c3e4d55..bf8c9f0c7 100644 --- a/src/stdlib_system.F90 +++ b/src/stdlib_system.F90 @@ -2,7 +2,8 @@ module stdlib_system use, intrinsic :: iso_c_binding, only : c_int, c_long, c_ptr, c_null_ptr, c_int64_t, c_size_t, & c_f_pointer use stdlib_kinds, only: int64, dp, c_bool, c_char -use stdlib_strings, only: to_c_char +use stdlib_strings, only: to_c_char, to_string +use stdlib_string_type, only: string_type use stdlib_error, only: state_type, STDLIB_SUCCESS, STDLIB_FS_ERROR implicit none private @@ -83,7 +84,15 @@ module stdlib_system public :: kill public :: elapsed public :: is_windows - + +!! Public path related functions and interfaces +public :: path_sep +public :: join_path +public :: operator(/) +public :: split_path +public :: base_name +public :: dir_name + !! version: experimental !! !! Tests if a given path matches an existing directory. @@ -133,6 +142,21 @@ module stdlib_system !! On Windows, this is `NUL`. On UNIX-like systems, this is `/dev/null`. !! public :: null_device + +!! version: experimental +!! +!! A helper function for returning the `type(state_type)` with the flag `STDLIB_FS_ERROR` set. +!! ([Specification](../page/specs/stdlib_system.html#FS_ERROR)) +!! +public :: FS_ERROR + +!! version: experimental +!! +!! A helper function for returning the `type(state_type)` with the flag `STDLIB_FS_ERROR` set. +!! It also formats and prefixes the `code` passed to it as the first argument +!! ([Specification](../page/specs/stdlib_system.html#FS_ERROR_CODE)) +!! +public :: FS_ERROR_CODE ! CPU clock ticks storage integer, parameter, private :: TICKS = int64 @@ -550,6 +574,141 @@ end function process_get_ID end interface +interface join_path + !! version: experimental + !! + !!### Summary + !! join the paths provided according to the OS-specific path-separator + !! ([Specification](../page/specs/stdlib_system.html#join_path)) + !! + module function join2_char_char(p1, p2) result(path) + character(:), allocatable :: path + character(*), intent(in) :: p1, p2 + end function join2_char_char + + module function join2_char_string(p1, p2) result(path) + character(:), allocatable :: path + character(*), intent(in) :: p1 + type(string_type), intent(in) :: p2 + end function join2_char_string + + module function join2_string_char(p1, p2) result(path) + type(string_type) :: path + type(string_type), intent(in) :: p1 + character(*), intent(in) :: p2 + end function join2_string_char + + module function join2_string_string(p1, p2) result(path) + type(string_type) :: path + type(string_type), intent(in) :: p1, p2 + end function join2_string_string + + module function joinarr_char(p) result(path) + character(:), allocatable :: path + character(*), intent(in) :: p(:) + end function joinarr_char + + module function joinarr_string(p) result(path) + type(string_type) :: path + type(string_type), intent(in) :: p(:) + end function joinarr_string +end interface join_path + +interface operator(/) + !! version: experimental + !! + !!### Summary + !! A binary operator to join the paths provided according to the OS-specific path-separator + !! ([Specification](../page/specs/stdlib_system.html#operator(/))) + !! + module function join_op_char_char(p1, p2) result(path) + character(:), allocatable :: path + character(*), intent(in) :: p1, p2 + end function join_op_char_char + + module function join_op_char_string(p1, p2) result(path) + character(:), allocatable :: path + character(*), intent(in) :: p1 + type(string_type), intent(in) :: p2 + end function join_op_char_string + + module function join_op_string_char(p1, p2) result(path) + type(string_type) :: path + type(string_type), intent(in) :: p1 + character(*), intent(in) :: p2 + end function join_op_string_char + + module function join_op_string_string(p1, p2) result(path) + type(string_type) :: path + type(string_type), intent(in) :: p1, p2 + end function join_op_string_string +end interface operator(/) + +interface split_path + !! version: experimental + !! + !!### Summary + !! splits the path immediately following the final path-separator + !! separating into typically a directory and a file name. + !! ([Specification](../page/specs/stdlib_system.html#split_path)) + !! + !!### Description + !! If the path is empty `head`='.' and tail='' + !! If the path only consists of separators, `head` is set to the separator and tail is empty + !! If the path is a root directory, `head` is set to that directory and tail is empty + !! `head` ends with a path-separator iff the path appears to be a root directory or a child of the root directory + module subroutine split_path_char(p, head, tail) + character(*), intent(in) :: p + character(:), allocatable, intent(out) :: head, tail + end subroutine split_path_char + + module subroutine split_path_string(p, head, tail) + type(string_type), intent(in) :: p + type(string_type), intent(out) :: head, tail + end subroutine split_path_string +end interface split_path + +interface base_name + !! version: experimental + !! + !!### Summary + !! returns the base name (last component) of the provided path + !! ([Specification](../page/specs/stdlib_system.html#base_name)) + !! + !!### Description + !! The value returned is the `tail` of the interface `split_path` + module function base_name_char(p) result(base) + character(:), allocatable :: base + character(*), intent(in) :: p + end function base_name_char + + module function base_name_string(p) result(base) + type(string_type) :: base + type(string_type), intent(in) :: p + end function base_name_string +end interface base_name + +interface dir_name + !! version: experimental + !! + !!### Summary + !! returns everything but the last component of the provided path + !! ([Specification](../page/specs/stdlib_system.html#dir_name)) + !! + !!### Description + !! The value returned is the `head` of the interface `split_path` + module function dir_name_char(p) result(dir) + character(:), allocatable :: dir + character(*), intent(in) :: p + end function dir_name_char + + module function dir_name_string(p) result(dir) + type(string_type) :: dir + type(string_type), intent(in) :: p + end function dir_name_string +end interface dir_name + + contains integer function get_runtime_os() result(os) @@ -770,4 +929,42 @@ subroutine delete_file(path, err) end if end subroutine delete_file +pure function FS_ERROR_CODE(code,a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,& + a11,a12,a13,a14,a15,a16,a17,a18,a19) result(state) + + type(state_type) :: state + !> Platform specific error code + integer, intent(in) :: code + !> Optional rank-agnostic arguments + class(*), intent(in), optional, dimension(..) :: a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,& + a11,a12,a13,a14,a15,a16,a17,a18,a19 + + character(32) :: code_msg + + write(code_msg, "('code - ', i0, ',')") code + + state = state_type(STDLIB_FS_ERROR, code_msg,a1,a2,a3,a4,a5,a6,a7,a8,& + a9,a10,a11,a12,a13,a14,a15,a16,a17,a18,a19) +end function FS_ERROR_CODE + +pure function FS_ERROR(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,& + a12,a13,a14,a15,a16,a17,a18,a19,a20) result(state) + + type(state_type) :: state + !> Optional rank-agnostic arguments + class(*), intent(in), optional, dimension(..) :: a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,& + a11,a12,a13,a14,a15,a16,a17,a18,a19,a20 + + state = state_type(STDLIB_FS_ERROR, a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,& + a13,a14,a15,a16,a17,a18,a19,a20) +end function FS_ERROR + +character function path_sep() + if (OS_TYPE() == OS_WINDOWS) then + path_sep = '\' + else + path_sep = '/' + end if +end function path_sep + end module stdlib_system diff --git a/src/stdlib_system_path.f90 b/src/stdlib_system_path.f90 new file mode 100644 index 000000000..c2ad1aec8 --- /dev/null +++ b/src/stdlib_system_path.f90 @@ -0,0 +1,170 @@ +submodule(stdlib_system) stdlib_system_path + use stdlib_ascii, only: reverse + use stdlib_strings, only: chomp, find, join + use stdlib_string_type, only: string_type, char, move +contains + module function join2_char_char(p1, p2) result(path) + character(:), allocatable :: path + character(*), intent(in) :: p1, p2 + + path = trim(p1) // path_sep() // trim(p2) + end function join2_char_char + + module function join2_char_string(p1, p2) result(path) + character(:), allocatable :: path + character(*), intent(in) :: p1 + type(string_type), intent(in) :: p2 + + path = join_path(p1, char(p2)) + end function join2_char_string + + module function join2_string_char(p1, p2) result(path) + type(string_type) :: path + type(string_type), intent(in) :: p1 + character(*), intent(in) :: p2 + character(:), allocatable :: join_char + + join_char = join_path(char(p1), p2) + + call move(join_char, path) + end function join2_string_char + + module function join2_string_string(p1, p2) result(path) + type(string_type) :: path + type(string_type), intent(in) :: p1, p2 + character(:), allocatable :: join_char + + join_char = join_path(char(p1), char(p2)) + + call move(join_char, path) + end function join2_string_string + + module function joinarr_char(p) result(path) + character(:), allocatable :: path + character(*), intent(in) :: p(:) + + path = join(p, path_sep()) + end function joinarr_char + + module function joinarr_string(p) result(path) + type(string_type) :: path + type(string_type), intent(in) :: p(:) + + path = join(p, path_sep()) + end function joinarr_string + + module function join_op_char_char(p1, p2) result(path) + character(:), allocatable :: path + character(*), intent(in) :: p1, p2 + + path = join_path(p1, p2) + end function join_op_char_char + + module function join_op_char_string(p1, p2) result(path) + character(:), allocatable :: path + character(*), intent(in) :: p1 + type(string_type), intent(in) :: p2 + + path = join_path(p1, p2) + end function join_op_char_string + + module function join_op_string_char(p1, p2) result(path) + type(string_type) :: path + type(string_type), intent(in) :: p1 + character(*), intent(in) :: p2 + + path = join_path(p1, p2) + end function join_op_string_char + + module function join_op_string_string(p1, p2) result(path) + type(string_type) :: path + type(string_type), intent(in) :: p1, p2 + + path = join_path(p1, p2) + end function join_op_string_string + + module subroutine split_path_char(p, head, tail) + character(*), intent(in) :: p + character(:), allocatable, intent(out) :: head, tail + character(:), allocatable :: temp + integer :: i + character(len=1) :: sep + sep = path_sep() + + ! Empty string, return (.,'') + if (trim(p) == '') then + head = '.' + tail = '' + return + end if + + ! Remove trailing path separators + temp = trim(chomp(trim(p), sep)) + + if (temp == '') then + head = sep + tail = '' + return + end if + + i = find(reverse(temp), sep) + + ! if no `pathsep`, then it probably was a root dir like `C:\` + if (i == 0) then + head = temp // sep + tail = '' + return + end if + + head = temp(:len(temp)-i) + + ! child of a root directory + if (find(head, sep) == 0) then + head = head // sep + end if + + tail = temp(len(temp)-i+2:) + end subroutine split_path_char + + module subroutine split_path_string(p, head, tail) + type(string_type), intent(in) :: p + type(string_type), intent(out) :: head, tail + + character(:), allocatable :: head_char, tail_char + + call split_path(char(p), head_char, tail_char) + + call move(head_char, head) + call move(tail_char, tail) + end subroutine split_path_string + + module function base_name_char(p) result(base) + character(:), allocatable :: base, temp + character(*), intent(in) :: p + + call split_path(p, temp, base) + end function base_name_char + + module function base_name_string(p) result(base) + type(string_type) :: base + type(string_type), intent(in) :: p + type(string_type) :: temp + + call split_path(p, temp, base) + end function base_name_string + + module function dir_name_char(p) result(dir) + character(:), allocatable :: dir, temp + character(*), intent(in) :: p + + call split_path(p, dir, temp) + end function dir_name_char + + module function dir_name_string(p) result(dir) + type(string_type) :: dir + type(string_type), intent(in) :: p + type(string_type) :: temp + + call split_path(p, dir, temp) + end function dir_name_string +end submodule stdlib_system_path diff --git a/test/linalg/test_linalg_specialmatrices.fypp b/test/linalg/test_linalg_specialmatrices.fypp index c445ec0da..5a09fb894 100644 --- a/test/linalg/test_linalg_specialmatrices.fypp +++ b/test/linalg/test_linalg_specialmatrices.fypp @@ -1,9 +1,11 @@ #:include "common.fypp" #:set R_KINDS_TYPES = list(zip(REAL_KINDS, REAL_TYPES, REAL_SUFFIX)) -#:set KINDS_TYPES = R_KINDS_TYPES +#:set C_KINDS_TYPES = list(zip(CMPLX_KINDS, CMPLX_TYPES, CMPLX_SUFFIX)) +#:set KINDS_TYPES = R_KINDS_TYPES+C_KINDS_TYPES module test_specialmatrices use testdrive, only : new_unittest, unittest_type, error_type, check, skip_test use stdlib_kinds + use stdlib_constants use stdlib_linalg, only: hermitian use stdlib_linalg_state, only: linalg_state_type use stdlib_math, only: all_close @@ -12,19 +14,141 @@ module test_specialmatrices contains - !> Collect all exported unit tests subroutine collect_suite(testsuite) !> Collection of tests type(unittest_type), allocatable, intent(out) :: testsuite(:) testsuite = [ & - new_unittest('tridiagonal', test_tridiagonal), & - new_unittest('tridiagonal error handling', test_tridiagonal_error_handling) & + new_unittest('tridiagonal arithmetic', test_tridiagonal_arithmetic), & + new_unittest('tridiagonal spmv kernel', test_tridiagonal_spmv), & + new_unittest('tridiagonal error handling', test_tridiagonal_error_handling), & + new_unittest('symtridiagonal arithmetic', test_symtridiagonal_arithmetic), & + new_unittest('symtridiagonal spmv kernel', test_symtridiagonal_spmv), & + new_unittest('hermtridiagonal arithmetic', test_hermtridiagonal_arithmetic), & + new_unittest('hermtridiagonal spmv kernel', test_hermtridiagonal_spmv) & ] end subroutine - subroutine test_tridiagonal(error) + !---------------------------------------- + !----- ----- + !----- TRIDIAGONAL MATRICES ----- + !----- ----- + !---------------------------------------- + + subroutine test_tridiagonal_arithmetic(error) + !> Error handling + type(error_type), allocatable, intent(out) :: error + #:for k1, t1, s1 in (KINDS_TYPES) + block + integer, parameter :: wp = ${k1}$ + integer, parameter :: n = 5 + type(tridiagonal_${s1}$_type) :: A, B, C + type(symtridiagonal_${s1}$_type) :: D + ${t1}$, allocatable :: Amat(:, :), Bmat(:, :), Cmat(:, :), Dmat(:, :) + ${t1}$, allocatable :: dl(:), dv(:), du(:), ev(:) + ${t1}$, parameter :: alpha = #{if t1.startswith('real')}# 2.0_wp #{else}# cmplx(2.0_wp, 2.0_wp, kind=wp) #{endif}# + integer :: i, j + #:if t1.startswith('complex') + type(hermtridiagonal_${s1}$_type) :: E + ${t1}$, allocatable :: Emat(:, :) + real(wp), allocatable :: data(:, :) + #:endif + + ! Initialize A and B matrices. + allocate(dl(n-1), dv(n), du(n-1), ev(n-1)) + #:if t1.startswith('real') + call random_number(dl) ; call random_number(dv) ; call random_number(du) + #:else + allocate(data(n, 2)) + call random_number(data) ; dl = cmplx(data(:n-1, 1), data(:n-1, 2), kind=wp) + call random_number(data) ; dv = cmplx(data(:, 1), data(:, 2), kind=wp) + call random_number(data) ; du = cmplx(data(:n-1, 1), data(:n-1, 2), kind=wp) + #:endif + A = tridiagonal(dl, dv, du) ; Amat = dense(A) + + #:if t1.startswith('real') + call random_number(dl) ; call random_number(dv) ; call random_number(du) + #:else + call random_number(data) ; dl = cmplx(data(:n-1, 1), data(:n-1, 2), kind=wp) + call random_number(data) ; dv = cmplx(data(:, 1), data(:, 2), kind=wp) + call random_number(data) ; du = cmplx(data(:n-1, 1), data(:n-1, 2), kind=wp) + #:endif + B = tridiagonal(dl, dv, du) ; Bmat = dense(B) + + #:if t1.startswith('real') + call random_number(dv) ; call random_number(ev) + #:else + call random_number(data) ; dv = cmplx(data(:, 1), data(:, 2), kind=wp) + call random_number(data) ; ev = cmplx(data(:n-1, 1), data(:n-1, 2), kind=wp) + #:endif + D = symtridiagonal(dv, ev) ; Dmat = dense(D) + + #:if t1.startswith('complex') + call random_number(data) ; dv = cmplx(data(:, 1), data(:, 2), kind=wp) + call random_number(data) ; ev = cmplx(data(:n-1, 1), data(:n-1, 2), kind=wp) + E = hermtridiagonal(dv, ev) ; Emat = dense(E) + #:endif + + ! Matrix addition. + C = A + B ; Cmat = dense(C) ! Tridiag + Tridiag = Tridiag + call check(error, all_close(Cmat, Amat+Bmat), .true.) + if (allocated(error)) return + + C = A + D ; Cmat = dense(C) ! Tridiag + SymTridiag = Tridiag + call check(error, all_close(Cmat, Amat+Dmat), .true.) + if (allocated(error)) return + + C = D + A ; Cmat = dense(C) ! SymTridiag + Tridiag = Tridiag + call check(error, all_close(Cmat, Amat+Dmat), .true.) + if (allocated(error)) return + + #:if t1.startswith("complex") + C = A + E ; Cmat = dense(C) ! Tridiag + HermTridiag = Tridiag + call check(error, all_close(Cmat, Amat+Emat), .true.) + if (allocated(error)) return + + C = E + A ; Cmat = dense(C) ! HermTridiag + Tridiag = Tridiag + call check(error, all_close(Cmat, Amat+Emat), .true.) + if (allocated(error)) return + #:endif + + ! Matrix subtraction. + C = A - B ; Cmat = dense(C) ! Tridiag - Tridiag = Tridiag + call check(error, all_close(Cmat, Amat-Bmat), .true.) + if (allocated(error)) return + + C = A - D ; Cmat = dense(C) ! Tridiag - SymTridiag = Tridiag + call check(error, all_close(Cmat, Amat-Dmat), .true.) + if (allocated(error)) return + + C = D - A ; Cmat = dense(C) ! SymTridiag - Tridiag = Tridiag + call check(error, all_close(Cmat, Dmat-Amat), .true.) + if (allocated(error)) return + + #:if t1.startswith("complex") + C = A - E ; Cmat = dense(C) ! Tridiag - HermTridiag = Tridiag + call check(error, all_close(Cmat, Amat-Emat), .true.) + if (allocated(error)) return + + C = E - A ; Cmat = dense(C) ! HermTridiag - Tridiag = Tridiag + call check(error, all_close(Cmat, Emat-Amat), .true.) + if (allocated(error)) return + #:endif + + ! Matrix scalar multiplication + C = alpha * A ; Cmat = dense(C) + call check(error, all_close(Cmat, alpha * Amat), .true.) + if (allocated(error)) return + + C = A * alpha ; Cmat = dense(C) + call check(error, all_close(Cmat, alpha * Amat), .true.) + if (allocated(error)) return + end block + #:endfor + end subroutine + + subroutine test_tridiagonal_spmv(error) !> Error handling type(error_type), allocatable, intent(out) :: error #:for k1, t1, s1 in (KINDS_TYPES) @@ -33,17 +157,33 @@ contains integer, parameter :: n = 5 type(tridiagonal_${s1}$_type) :: A ${t1}$, allocatable :: Amat(:,:), dl(:), dv(:), du(:) + #: if t1.startswith('complex') + real(wp), allocatable :: data(:, :) + #:endif ${t1}$, allocatable :: x(:) ${t1}$, allocatable :: y1(:), y2(:) ! Initialize matrix. allocate(dl(n-1), dv(n), du(n-1)) + #:if t1.startswith('real') call random_number(dl) ; call random_number(dv) ; call random_number(du) + #:else + allocate(data(n, 2), source=0.0_wp) + call random_number(data) ; dl = cmplx(data(:n-1, 1), data(:n-1, 2), kind=wp) + call random_number(data) ; dv = cmplx(data(:, 1), data(:, 2), kind=wp) + call random_number(data) ; du = cmplx(data(:n-1, 1), data(:n-1, 2), kind=wp) + #:endif A = tridiagonal(dl, dv, du) ; Amat = dense(A) ! Random vectors. + #:if t1.startswith('real') allocate(x(n), source = 0.0_wp) ; call random_number(x) allocate(y1(n), source = 0.0_wp) ; allocate(y2(n), source=0.0_wp) + #:else + allocate(x(n), source=zero_c${k1}$) + call random_number(data) ; x = cmplx(data(:, 1), data(:, 2), kind=wp) + allocate(y1(n), source = zero_c${k1}$) ; allocate(y2(n), source=zero_c${k1}$) + #:endif ! Test y = A @ x y1 = matmul(Amat, x) ; call spmv(A, x, y2) @@ -94,6 +234,266 @@ contains #:endfor end subroutine + !-------------------------------------------------- + !----- ----- + !----- SYMMETRIC TRIDIAGONAL MATRICES ----- + !----- ----- + !-------------------------------------------------- + + subroutine test_symtridiagonal_arithmetic(error) + !> Error handling + type(error_type), allocatable, intent(out) :: error + #:for k1, t1, s1 in (KINDS_TYPES) + block + integer, parameter :: wp = ${k1}$ + integer, parameter :: n = 5 + type(symtridiagonal_${s1}$_type) :: A, B, C + type(tridiagonal_${s1}$_type) :: D + ${t1}$, allocatable :: Amat(:, :), Bmat(:, :), Cmat(:, :) + ${t1}$, allocatable :: dv(:), ev(:) + ${t1}$, parameter :: alpha = #{if t1.startswith("real")}# 2.0_wp #{else}# cmplx(2.0_wp, 2.0_wp, kind=wp) #{endif}# + integer :: i, j + #:if t1.startswith('complex') + type(hermtridiagonal_${s1}$_type) :: E + ${t1}$, allocatable :: Emat(:, :) + real(wp), allocatable :: data(:, :) + #:endif + + ! Initialize A and B matrices. + allocate(ev(n-1), dv(n)) + #:if t1.startswith('real') + call random_number(ev) ; call random_number(dv) + #:else + allocate(data(n, 2)) + call random_number(data) ; dv = cmplx(data(:, 1), data(:, 2), kind=wp) + call random_number(data) ; ev = cmplx(data(:n-1, 1), data(:n-1, 2), kind=wp) + #:endif + A = symtridiagonal(dv, ev) ; Amat = dense(A) + + #:if t1.startswith('real') + call random_number(ev) ; call random_number(dv) + #:else + call random_number(data) ; dv = cmplx(data(:, 1), data(:, 2), kind=wp) + call random_number(data) ; ev = cmplx(data(:n-1, 1), data(:n-1, 2), kind=wp) + #:endif + B = symtridiagonal(dv, ev) ; Bmat = dense(B) + + #:if t1.startswith('complex') + call random_number(data) ; dv = cmplx(data(:, 1), data(:, 2), kind=wp) + call random_number(data) ; ev = cmplx(data(:n-1, 1), data(:n-1, 2), kind=wp) + E = hermtridiagonal(dv, ev) ; Emat = dense(E) + #:endif + + ! Matrix addition. + C = A + B ; Cmat = dense(C) ! SymTridiag + SymTridiag = SymTridiag + call check(error, all_close(Cmat, Amat+Bmat), .true.) + if (allocated(error)) return + + #:if t1.startswith("complex") + D = A + E ; Cmat = dense(D) ! SymTridiag + HermTridiag = Tridiag + call check(error, all_close(Cmat, Amat+Emat), .true.) + if (allocated(error)) return + + D = E + A ; Cmat = dense(D) ! HermTridiag + SymTridiag = Tridiag + call check(error, all_close(Cmat, Amat+Emat), .true.) + if (allocated(error)) return + #:endif + + ! Matrix subtraction. + C = A - B ; Cmat = dense(C) ! SymTridiag - SymTridiag = SymTridiag + call check(error, all_close(Cmat, Amat-Bmat), .true.) + if (allocated(error)) return + + #:if t1.startswith("complex") + D = A - E ; Cmat = dense(D) ! SymTridiag - HermTridiag = Tridiag + call check(error, all_close(Cmat, Amat-Emat), .true.) + if (allocated(error)) return + + D = E - A ; Cmat = dense(D) ! HermTridiag - SymTridiag = Tridiag + call check(error, all_close(Cmat, Emat-Amat), .true.) + if (allocated(error)) return + #:endif + + ! Matrix scalar multiplication + C = alpha * A ; Cmat = dense(C) + call check(error, all_close(Cmat, alpha * Amat), .true.) + if (allocated(error)) return + + C = A * alpha ; Cmat = dense(C) + call check(error, all_close(Cmat, alpha * Amat), .true.) + if (allocated(error)) return + end block + #:endfor + end subroutine + + subroutine test_symtridiagonal_spmv(error) + !> Error handling + type(error_type), allocatable, intent(out) :: error + #:for k1, t1, s1 in (KINDS_TYPES) + block + integer, parameter :: wp = ${k1}$ + integer, parameter :: n = 5 + type(symtridiagonal_${s1}$_type) :: A + ${t1}$, allocatable :: Amat(:,:), dv(:), ev(:) + #:if t1.startswith('complex') + real(wp), allocatable :: data(:, :) + #:endif + ${t1}$, allocatable :: x(:) + ${t1}$, allocatable :: y1(:), y2(:) + + ! Initialize matrix. + allocate(ev(n-1), dv(n)) + #:if t1.startswith('real') + call random_number(dv) ; call random_number(ev) + #:else + allocate(data(n, 2), source=0.0_wp) + call random_number(data) ; dv = cmplx(data(:, 1), data(:, 2), kind=wp) + call random_number(data) ; ev = cmplx(data(:n-1, 1), data(:n-1, 2), kind=wp) + #:endif + A = symtridiagonal(dv, ev) ; Amat = dense(A) + + ! Random vectors. + #:if t1.startswith('real') + allocate(x(n), source = 0.0_wp) ; call random_number(x) + allocate(y1(n), source = 0.0_wp) ; allocate(y2(n), source=0.0_wp) + #:else + allocate(x(n), source=zero_c${k1}$) + call random_number(data) ; x = cmplx(data(:, 1), data(:, 2), kind=wp) + allocate(y1(n), source = zero_c${k1}$) ; allocate(y2(n), source=zero_c${k1}$) + #:endif + + ! Test y = A @ x + y1 = matmul(Amat, x) ; call spmv(A, x, y2) + call check(error, all_close(y1, y2), .true.) + if (allocated(error)) return + + ! Test y = A.T @ x + y1 = 0.0_wp ; y2 = 0.0_wp + y1 = matmul(transpose(Amat), x) ; call spmv(A, x, y2, op="T") + call check(error, all_close(y1, y2), .true.) + if (allocated(error)) return + + #:if t1.startswith('complex') + ! Test y = A.H @ x + y1 = 0.0_wp ; y2 = 0.0_wp + y1 = matmul(hermitian(Amat), x) ; call spmv(A, x, y2, op="H") + call check(error, all_close(y1, y2), .true.) + if (allocated(error)) return + #:endif + end block + #:endfor + end subroutine + + !-------------------------------------- + !----- ----- + !----- HERMITIAN MATRICES ----- + !----- ----- + !-------------------------------------- + + subroutine test_hermtridiagonal_arithmetic(error) + !> Error handling + type(error_type), allocatable, intent(out) :: error + #:for k1, t1, s1 in (C_KINDS_TYPES) + block + integer, parameter :: wp = ${k1}$ + integer, parameter :: n = 5 + type(hermtridiagonal_${s1}$_type) :: A, B, C + type(tridiagonal_${s1}$_type) :: D + ${t1}$, allocatable :: Amat(:, :), Bmat(:, :), Cmat(:, :), Dmat(:, :) + ${t1}$, allocatable :: dv(:), ev(:) + ${t1}$, parameter :: alpha = cmplx(2.0_wp, 2.0_wp, kind=wp) + real(${k1}$), parameter :: beta = 2.0_wp + integer :: i, j + real(wp), allocatable :: data(:, :) + + ! Initialize A and B matrices. + allocate(ev(n-1), dv(n)) + allocate(data(n, 2)) + call random_number(data) ; dv = cmplx(data(:, 1), data(:, 2), kind=wp) + call random_number(data) ; ev = cmplx(data(:n-1, 1), data(:n-1, 2), kind=wp) + A = hermtridiagonal(dv, ev) ; Amat = dense(A) + + call random_number(data) ; dv = cmplx(data(:, 1), data(:, 2), kind=wp) + call random_number(data) ; ev = cmplx(data(:n-1, 1), data(:n-1, 2), kind=wp) + B = hermtridiagonal(dv, ev) ; Bmat = dense(B) + + ! Matrix addition. + C = A + B ; Cmat = dense(C) + call check(error, all_close(Cmat, Amat+Bmat), .true.) + if (allocated(error)) return + + ! Matrix subtraction. + C = A - B ; Cmat = dense(C) + call check(error, all_close(Cmat, Amat-Bmat), .true.) + if (allocated(error)) return + + ! Matrix scalar multiplication + D = alpha * A ; Dmat = dense(D) + call check(error, all_close(Dmat, alpha * Amat), .true.) + if (allocated(error)) return + + D = A * alpha ; Dmat = dense(D) + call check(error, all_close(Dmat, alpha * Amat), .true.) + if (allocated(error)) return + + C = beta * A ; Cmat = dense(C) + call check(error, all_close(Cmat, beta * Amat), .true.) + if (allocated(error)) return + + C = A * beta ; Cmat = dense(C) + call check(error, all_close(Cmat, beta * Amat), .true.) + if (allocated(error)) return + end block + #:endfor + end subroutine + + subroutine test_hermtridiagonal_spmv(error) + !> Error handling + type(error_type), allocatable, intent(out) :: error + #:for k1, t1, s1 in (C_KINDS_TYPES) + block + integer, parameter :: wp = ${k1}$ + integer, parameter :: n = 5 + type(hermtridiagonal_${s1}$_type) :: A + ${t1}$, allocatable :: Amat(:,:), dv(:), ev(:) + real(wp), allocatable :: data(:, :) + ${t1}$, allocatable :: x(:) + ${t1}$, allocatable :: y1(:), y2(:) + + ! Initialize matrix. + allocate(ev(n-1), dv(n)) + allocate(data(n, 2), source=0.0_wp) + call random_number(data) ; dv = cmplx(data(:, 1), data(:, 2), kind=wp) + call random_number(data) ; ev = cmplx(data(:n-1, 1), data(:n-1, 2), kind=wp) + A = hermtridiagonal(dv, ev) ; Amat = dense(A) + + ! Random vectors. + allocate(x(n), source=zero_c${k1}$) + call random_number(data) ; x = cmplx(data(:, 1), data(:, 2), kind=wp) + allocate(y1(n), source = zero_c${k1}$) ; allocate(y2(n), source=zero_c${k1}$) + + ! Test y = A @ x + y1 = matmul(Amat, x) ; call spmv(A, x, y2) + call check(error, all_close(y1, y2), .true.) + if (allocated(error)) return + + ! Test y = A.T @ x + y1 = 0.0_wp ; y2 = 0.0_wp + y1 = matmul(transpose(Amat), x) ; call spmv(A, x, y2, op="T") + call check(error, all_close(y1, y2), .true.) + if (allocated(error)) return + + #:if t1.startswith('complex') + ! Test y = A.H @ x + y1 = 0.0_wp ; y2 = 0.0_wp + y1 = matmul(hermitian(Amat), x) ; call spmv(A, x, y2, op="H") + call check(error, all_close(y1, y2), .true.) + if (allocated(error)) return + #:endif + end block + #:endfor + end subroutine + end module @@ -109,7 +509,7 @@ program tester stat = 0 testsuites = [ & - new_testsuite("sparse", collect_suite) & + new_testsuite("Tridiagonal matrices", collect_suite) & ] do is = 1, size(testsuites) diff --git a/test/system/CMakeLists.txt b/test/system/CMakeLists.txt index b7623ea83..0ab568a18 100644 --- a/test/system/CMakeLists.txt +++ b/test/system/CMakeLists.txt @@ -2,3 +2,4 @@ ADDTEST(filesystem) ADDTEST(os) ADDTEST(sleep) ADDTEST(subprocess) +ADDTEST(path) diff --git a/test/system/test_filesystem.f90 b/test/system/test_filesystem.f90 index 4cf1690e4..838ced263 100644 --- a/test/system/test_filesystem.f90 +++ b/test/system/test_filesystem.f90 @@ -1,7 +1,7 @@ module test_filesystem use testdrive, only : new_unittest, unittest_type, error_type, check, skip_test - use stdlib_system, only: is_directory, delete_file - use stdlib_error, only: state_type + use stdlib_system, only: is_directory, delete_file, FS_ERROR, FS_ERROR_CODE + use stdlib_error, only: state_type, STDLIB_FS_ERROR implicit none @@ -13,6 +13,7 @@ subroutine collect_suite(testsuite) type(unittest_type), allocatable, intent(out) :: testsuite(:) testsuite = [ & + new_unittest("fs_error", test_fs_error), & new_unittest("fs_is_directory_dir", test_is_directory_dir), & new_unittest("fs_is_directory_file", test_is_directory_file), & new_unittest("fs_delete_non_existent", test_delete_file_non_existent), & @@ -21,6 +22,26 @@ subroutine collect_suite(testsuite) ] end subroutine collect_suite + subroutine test_fs_error(error) + type(error_type), allocatable, intent(out) :: error + type(state_type) :: s1, s2 + character(:), allocatable :: msg + + msg = "code - 10, Cannot create File temp.txt - File already exists" + s1 = FS_ERROR_CODE(10, "Cannot create File temp.txt -", "File already exists") + + call check(error, s1%state == STDLIB_FS_ERROR .and. s1%message == msg, & + "FS_ERROR_CODE: Could not construct the state with code correctly") + if (allocated(error)) return + + msg = "Cannot create File temp.txt - File already exists" + s2 = FS_ERROR("Cannot create File temp.txt -", "File already exists") + + call check(error, s2%state == STDLIB_FS_ERROR .and. s2%message == msg, & + "FS_ERROR: Could not construct state without code correctly") + if (allocated(error)) return + end subroutine test_fs_error + ! Test `is_directory` for a directory subroutine test_is_directory_dir(error) type(error_type), allocatable, intent(out) :: error diff --git a/test/system/test_path.f90 b/test/system/test_path.f90 new file mode 100644 index 000000000..8d892b928 --- /dev/null +++ b/test/system/test_path.f90 @@ -0,0 +1,149 @@ +module test_path + use testdrive, only : new_unittest, unittest_type, error_type, check, skip_test + use stdlib_system, only: join_path, operator(/), split_path, OS_TYPE, OS_WINDOWS + implicit none +contains + !> Collect all exported unit tests + subroutine collect_suite(testsuite) + !> Collection of tests + type(unittest_type), allocatable, intent(out) :: testsuite(:) + + testsuite = [ & + new_unittest('test_join_path', test_join_path), & + new_unittest('test_join_path_operator', test_join_path_op), & + new_unittest('test_split_path', test_split_path) & + ] + end subroutine collect_suite + + subroutine checkpath(error, funcname, expected, got) + type(error_type), allocatable, intent(out) :: error + character(len=*), intent(in) :: funcname + character(len=*), intent(in) :: expected + character(len=:), allocatable :: got + character(len=:), allocatable :: message + + message = "'"//funcname//"'"//" error: Expected '"// expected // "' but got '" // got // "'" + call check(error, expected == got, message) + + end subroutine checkpath + + subroutine test_join_path(error) + type(error_type), allocatable, intent(out) :: error + character(len=:), allocatable :: path + character(len=20) :: paths(5) + + if (OS_TYPE() == OS_WINDOWS) then + path = join_path('C:\Users', 'Alice') + call checkpath(error, 'join_path', 'C:\Users\Alice', path) + if (allocated(error)) return + + paths = [character(20) :: 'C:','Users','Bob','Pictures','2025'] + path = join_path(paths) + + call checkpath(error, 'join_path', 'C:\Users\Bob\Pictures\2025', path) + if (allocated(error)) return + + path = join_path('"C:\Users\John Doe"', 'Pictures\2025') ! path with spaces + call checkpath(error, 'join_path', '"C:\Users\John Doe"\Pictures\2025', path) + if (allocated(error)) return + else + path = join_path('/home', 'Alice') + call checkpath(error, 'join_path', '/home/Alice', path) + if (allocated(error)) return + + paths = [character(20) :: '','home','Bob','Pictures','2025'] + path = join_path(paths) + + call checkpath(error, 'join_path', '/home/Bob/Pictures/2025', path) + if (allocated(error)) return + end if + end subroutine test_join_path + + !> Test the operator + subroutine test_join_path_op(error) + type(error_type), allocatable, intent(out) :: error + character(len=:), allocatable :: path + + if (OS_TYPE() == OS_WINDOWS) then + path = 'C:'/'Users'/'Alice'/'Desktop' + call checkpath(error, 'join_path operator', 'C:\Users\Alice\Desktop', path) + if (allocated(error)) return + else + path = ''/'home'/'Alice'/'.config' + call checkpath(error, 'join_path operator', '/home/Alice/.config', path) + if (allocated(error)) return + end if + end subroutine test_join_path_op + + subroutine test_split_path(error) + type(error_type), allocatable, intent(out) :: error + character(len=:), allocatable :: head, tail + + call split_path('', head, tail) + call checkpath(error, 'split_path-head', '.', head) + if (allocated(error)) return + call checkpath(error, 'split_path-tail', '', tail) + if (allocated(error)) return + + if (OS_TYPE() == OS_WINDOWS) then + call split_path('\\\\', head, tail) + call checkpath(error, 'split_path-head', '\', head) + if (allocated(error)) return + call checkpath(error, 'split_path-tail', '', tail) + if (allocated(error)) return + + call split_path('C:\', head, tail) + call checkpath(error, 'split_path-head', 'C:\', head) + if (allocated(error)) return + call checkpath(error, 'split_path-tail', '', tail) + if (allocated(error)) return + + call split_path('C:\Users\Alice\\\\\', head, tail) + call checkpath(error, 'split_path-head', 'C:\Users', head) + if (allocated(error)) return + call checkpath(error, 'split_path-tail', 'Alice', tail) + if (allocated(error)) return + else + call split_path('/////', head, tail) + call checkpath(error, 'split_path-head', '/', head) + if (allocated(error)) return + call checkpath(error, 'split_path-tail', '', tail) + if (allocated(error)) return + + call split_path('/home/Alice/foo/bar.f90///', head, tail) + call checkpath(error, 'split_path-head', '/home/Alice/foo', head) + if (allocated(error)) return + call checkpath(error, 'split_path-tail', 'bar.f90', tail) + if (allocated(error)) return + end if + end subroutine test_split_path + +end module test_path + +program tester + use, intrinsic :: iso_fortran_env, only : error_unit + use testdrive, only : run_testsuite, new_testsuite, testsuite_type + use test_path, only : collect_suite + + implicit none + + integer :: stat, is + type(testsuite_type), allocatable :: testsuites(:) + character(len=*), parameter :: fmt = '("#", *(1x, a))' + + stat = 0 + + testsuites = [ & + new_testsuite("path", collect_suite) & + ] + + do is = 1, size(testsuites) + write(error_unit, fmt) "Testing:", testsuites(is)%name + call run_testsuite(testsuites(is)%collect, error_unit, stat) + end do + + if (stat > 0) then + write(error_unit, '(i0, 1x, a)') stat, "test(s) failed!" + error stop + end if +end program tester