Merge branch 'sparse_sellc' of https://github.com/jalvesz/stdlib into sparse_sellc

Author: jalvesz

doc/specs/stdlib_linalg.md

@@ -239,37 +239,34 @@ Pure function.
 
 ### Description
 
-Construct the identity matrix.
+Constructs the identity matrix.
 
 ### Syntax
 
-`I = ` [[stdlib_linalg(module):eye(function)]] `(dim1 [, dim2])`
+`I = ` [[stdlib_linalg(module):eye(function)]] `(dim1 [, dim2] [, mold])`
 
 ### Arguments
 
-`dim1`: Shall be a scalar of default type `integer`.
-This is an `intent(in)` argument. 
-
-`dim2`: Shall be a scalar of default type `integer`.
-This is an `intent(in)` and `optional` argument. 
+- `dim1`: A scalar of type `integer`. This is an `intent(in)` argument and specifies the number of rows.
+- `dim2`: A scalar of type `integer`. This is an optional `intent(in)` argument specifying the number of columns. If not provided, the matrix is square (`dim1 = dim2`).
+- `mold`: A scalar of any supported `integer`, `real`, or `complex` type. This is an optional `intent(in)` argument. If provided, the returned identity matrix will have the same type and kind as `mold`. If not provided, the matrix will be of type `real(real64)` by default.
 
 ### Return value
 
-Return the identity matrix, i.e. a matrix with ones on the main diagonal and zeros elsewhere. The return value is of type `integer(int8)`.
-The use of `int8` was suggested to save storage.
+Returns the identity matrix, with ones on the main diagonal and zeros elsewhere. 
+
+- By default, the return value is of type `real(real64)`, which is recommended for arithmetic safety.
+- If the `mold` argument is provided, the return value will match the type and kind of `mold`, allowing for arbitrary `integer`, `real`, or `complex` return types.
 
-#### Warning
+### Example
 
-Since the result of `eye` is of `integer(int8)` type, one should be careful about using it in arithmetic expressions. For example:
 ```fortran
-!> Be careful
-A = eye(2,2)/2     !! A == 0.0
-!> Recommend
-A = eye(2,2)/2.0   !! A == diag([0.5, 0.5])
+!> Return default type (real64)
+A = eye(2,2)/2             !! A == diag([0.5_dp, 0.5_dp])
+!> Return 32-bit complex
+A = eye(2,2, mold=(0.0,0.0))/2   !! A == diag([(0.5,0.5), (0.5,0.5)])
 ```
 
-### Example
-
 ```fortran
 {!example/linalg/example_eye1.f90!}
 ```
@@ -509,6 +506,37 @@ Returns a `logical` scalar that is `.true.` if the input matrix is skew-symmetri
 {!example/linalg/example_is_skew_symmetric.f90!}
 ```
 
+## `hermitian` - Compute the Hermitian version of a rank-2 matrix
+
+### Status
+
+Experimental
+
+### Description
+
+Compute the Hermitian version of a rank-2 matrix. 
+For `complex` matrices, the function returns the conjugate transpose (`conjg(transpose(a))`). 
+For `real` or `integer` matrices, the function returns the transpose (`transpose(a)`).
+
+### Syntax
+
+`h = ` [[stdlib_linalg(module):hermitian(interface)]] `(a)`
+
+### Arguments
+
+`a`: Shall be a rank-2 array of type `integer`, `real`, or `complex`. The input matrix `a` is not modified.
+
+### Return value
+
+Returns a rank-2 array of the same shape and type as `a`. If `a` is of type `complex`, the Hermitian matrix is computed as `conjg(transpose(a))`. 
+For `real` or `integer` types, it is equivalent to the intrinsic `transpose(a)`.
+
+### Example
+
+```fortran
+{!example/linalg/example_hermitian.f90!}
+```
+
 ## `is_hermitian` - Checks if a matrix is Hermitian
 
 ### Status
@@ -1459,6 +1487,142 @@ If `err` is not present, exceptions trigger an `error stop`.
 {!example/linalg/example_inverse_function.f90!}
 ```
 
+## `pinv` - Moore-Penrose pseudo-inverse of a matrix
+
+### Status
+
+Experimental
+
+### Description
+
+This function computes the Moore-Penrose pseudo-inverse of a `real` or `complex` matrix.  
+The pseudo-inverse, \( A^{+} \), generalizes the matrix inverse and satisfies the conditions:
+- \( A \cdot A^{+} \cdot A = A \)
+- \( A^{+} \cdot A \cdot A^{+} = A^{+} \)
+- \( (A \cdot A^{+})^T = A \cdot A^{+} \)
+- \( (A^{+} \cdot A)^T = A^{+} \cdot A \)
+
+The computation is based on singular value decomposition (SVD). Singular values below a relative 
+tolerance threshold \( \text{rtol} \cdot \sigma_{\max} \), where \( \sigma_{\max} \) is the largest 
+singular value, are treated as zero.
+
+### Syntax
+
+`b =` [[stdlib_linalg(module):pinv(interface)]] `(a, [, rtol, err])`
+
+### Arguments
+
+`a`: Shall be a rank-2, `real` or `complex` array of shape `[m, n]` containing the coefficient matrix. 
+It is an `intent(in)` argument.
+
+`rtol` (optional): Shall be a scalar `real` value specifying the relative tolerance for singular value cutoff.  
+If `rtol` is not provided, the default relative tolerance is \( \text{rtol} = \text{max}(m, n) \cdot \epsilon \),  
+where \( \epsilon \) is the machine precision for the element type of `a`. It is an `intent(in)` argument.
+
+`err` (optional): Shall be a `type(linalg_state_type)` value. It is an `intent(out)` argument.
+
+### Return value
+
+Returns an array value of the same type, kind, and rank as `a` with shape `[n, m]`, that contains the pseudo-inverse matrix \( A^{+} \).
+
+Raises `LINALG_ERROR` if the underlying SVD did not converge.
+Raises `LINALG_VALUE_ERROR` if `a` has invalid size.
+If `err` is not present, exceptions trigger an `error stop`.
+
+### Example
+
+```fortran
+{!example/linalg/example_pseudoinverse.f90!}
+```
+
+## `pseudoinvert` - Moore-Penrose pseudo-inverse of a matrix
+
+### Status
+
+Experimental
+
+### Description
+
+This subroutine computes the Moore-Penrose pseudo-inverse of a `real` or `complex` matrix.
+The pseudo-inverse \( A^{+} \) is a generalization of the matrix inverse and satisfies the following properties:
+- \( A \cdot A^{+} \cdot A = A \)
+- \( A^{+} \cdot A \cdot A^{+} = A^{+} \)
+- \( (A \cdot A^{+})^T = A \cdot A^{+} \)
+- \( (A^{+} \cdot A)^T = A^{+} \cdot A \)
+
+The computation is based on singular value decomposition (SVD). Singular values below a relative 
+tolerance threshold \( \text{rtol} \cdot \sigma_{\max} \), where \( \sigma_{\max} \) is the largest 
+singular value, are treated as zero.
+
+On return, matrix `pinva` `[n, m]` will store the pseudo-inverse of `a` `[m, n]`.
+
+### Syntax
+
+`call ` [[stdlib_linalg(module):pseudoinvert(interface)]] `(a, pinva [, rtol] [, err])`
+
+### Arguments
+
+`a`: Shall be a rank-2, `real` or `complex` array containing the coefficient matrix.  
+It is an `intent(in)` argument.  
+
+`pinva`: Shall be a rank-2 array of the same kind as `a`, and size equal to that of `transpose(a)`.  
+On output, it contains the Moore-Penrose pseudo-inverse of `a`.
+
+`rtol` (optional): Shall be a scalar `real` value specifying the relative tolerance for singular value cutoff.  
+If not provided, the default threshold is \( \text{max}(m, n) \cdot \epsilon \), where \( \epsilon \) is the 
+machine precision for the element type of `a`.
+
+`err` (optional): Shall be a `type(linalg_state_type)` value. It is an `intent(out)` argument.
+
+### Return value
+
+Computes the Moore-Penrose pseudo-inverse of the matrix \( A \), \( A^{+} \), and returns it in matrix `pinva`.
+
+Raises `LINALG_ERROR` if the underlying SVD did not converge.
+Raises `LINALG_VALUE_ERROR` if `pinva` and `a` have degenerate or incompatible sizes.
+If `err` is not present, exceptions trigger an `error stop`.
+
+### Example
+
+```fortran
+{!example/linalg/example_pseudoinverse.f90!}
+```
+
+## `.pinv.` - Moore-Penrose Pseudo-Inverse operator
+
+### Status
+
+Experimental
+
+### Description
+
+This operator returns the Moore-Penrose pseudo-inverse of a `real` or `complex` matrix \( A \).
+The pseudo-inverse \( A^{+} \) is computed using Singular Value Decomposition (SVD), and singular values 
+below a given threshold are treated as zero.
+
+This interface is equivalent to the function [[stdlib_linalg(module):pinv(interface)]].
+
+### Syntax
+
+`b = ` [[stdlib_linalg(module):operator(.pinv.)(interface)]] `a`
+
+### Arguments
+
+`a`: Shall be a rank-2 array of any `real` or `complex` kinds, with arbitrary dimensions \( m \times n \). It is an `intent(in)` argument.
+
+### Return value
+
+Returns a rank-2 array with the same type, kind, and rank as `a`, that contains the Moore-Penrose pseudo-inverse of `a`.
+
+If an exception occurs, or if the input matrix is degenerate (e.g., rank-deficient), the returned matrix will contain `NaN`s.
+For more detailed error handling, it is recommended to use the `subroutine` or `function` interfaces.
+
+### Example
+
+```fortran
+{!example/linalg/example_pseudoinverse.f90!}
+```
+
 ## `get_norm` - Computes the vector norm of a generic-rank array.
 
 ### Status

example/linalg/CMakeLists.txt

@@ -6,6 +6,7 @@ ADD_EXAMPLE(diag5)
 ADD_EXAMPLE(eye1)
 ADD_EXAMPLE(eye2)
 ADD_EXAMPLE(is_diagonal)
+ADD_EXAMPLE(hermitian)
 ADD_EXAMPLE(is_hermitian)
 ADD_EXAMPLE(is_hessenberg)
 ADD_EXAMPLE(is_skew_symmetric)
@@ -16,6 +17,7 @@ ADD_EXAMPLE(inverse_operator)
 ADD_EXAMPLE(inverse_function)
 ADD_EXAMPLE(inverse_inplace)
 ADD_EXAMPLE(inverse_subroutine)
+ADD_EXAMPLE(pseudoinverse)
 ADD_EXAMPLE(outer_product)
 ADD_EXAMPLE(eig)
 ADD_EXAMPLE(eigh)

example/linalg/example_eye1.f90

@@ -5,10 +5,10 @@ program example_eye1
   real :: a(3, 3)
   real :: b(2, 3)  !! Matrix is non-square.
   complex :: c(2, 2)
-  I = eye(2)              !! [1,0; 0,1]
-  A = eye(3)              !! [1.0,0.0,0.0; 0.0,1.0,0.0; 0.0,0.0,1.0]
-  A = eye(3, 3)            !! [1.0,0.0,0.0; 0.0,1.0,0.0; 0.0,0.0,1.0]
-  B = eye(2, 3)            !! [1.0,0.0,0.0; 0.0,1.0,0.0]
-  C = eye(2, 2)            !! [(1.0,0.0),(0.0,0.0); (0.0,0.0),(1.0,0.0)]
+  I = eye(2)                !! [1,0; 0,1]
+  A = eye(3)                !! [1.0,0.0,0.0; 0.0,1.0,0.0; 0.0,0.0,1.0]
+  A = eye(3, 3)             !! [1.0,0.0,0.0; 0.0,1.0,0.0; 0.0,0.0,1.0]
+  B = eye(2, 3)             !! [1.0,0.0,0.0; 0.0,1.0,0.0]
+  C = eye(2, 2)             !! [(1.0,0.0),(0.0,0.0); (0.0,0.0),(1.0,0.0)]
   C = (1.0, 1.0)*eye(2, 2)  !! [(1.0,1.0),(0.0,0.0); (0.0,0.0),(1.0,1.0)]
 end program example_eye1

example/linalg/example_hermitian.f90

@@ -0,0 +1,19 @@
+! Example program demonstrating the usage of hermitian
+program example_hermitian
+  use stdlib_linalg, only: hermitian
+  implicit none
+
+  complex, dimension(2, 2) :: A, AT
+
+  ! Define input matrices
+  A = reshape([(1,2),(3,4),(5,6),(7,8)],[2,2])
+
+  ! Compute Hermitian matrices
+  AT = hermitian(A)
+
+  print *, "Original Complex Matrix:"
+  print *, A
+  print *, "Hermitian Complex Matrix:"
+  print *, AT
+  
+end program example_hermitian

example/linalg/example_pseudoinverse.f90

@@ -0,0 +1,32 @@
+! Matrix pseudo-inversion example: function, subroutine, and operator interfaces
+program example_pseudoinverse
+  use stdlib_linalg, only: pinv, pseudoinvert, operator(.pinv.), linalg_state_type
+  implicit none(type,external)
+
+  real :: A(15,5), Am1(5,15)
+  type(linalg_state_type) :: state
+  integer :: i, j
+  real, parameter :: tol = sqrt(epsilon(0.0))
+  
+  ! Generate random matrix A (15x15)
+  call random_number(A)
+
+  ! Pseudo-inverse: Function interfcae
+  Am1 = pinv(A, err=state)
+  print *, 'Max error (function)  : ', maxval(abs(A-matmul(A, matmul(Am1,A))))
+  
+  ! User threshold
+  Am1 = pinv(A, rtol=0.001, err=state)
+  print *, 'Max error (rtol=0.001): ', maxval(abs(A-matmul(A, matmul(Am1,A))))  
+
+  ! Pseudo-inverse: Subroutine interface
+  call pseudoinvert(A, Am1, err=state)
+
+  print *, 'Max error (subroutine): ', maxval(abs(A-matmul(A, matmul(Am1,A))))
+
+  ! Operator interface
+  Am1 = .pinv.A
+  
+  print *, 'Max error (operator)  : ', maxval(abs(A-matmul(A, matmul(Am1,A))))
+
+end program example_pseudoinverse

src/CMakeLists.txt

@@ -32,6 +32,7 @@ set(fppFiles
     stdlib_linalg_determinant.fypp
     stdlib_linalg_qr.fypp
     stdlib_linalg_inverse.fypp
+    stdlib_linalg_pinv.fypp
     stdlib_linalg_norms.fypp
     stdlib_linalg_state.fypp
     stdlib_linalg_svd.fypp 

src/stdlib_io.fypp

@@ -167,7 +167,7 @@ contains
           read (s,*,iostat=ios,iomsg=iomsg) d(i, :)
 
           if (ios/=0) then 
-             write(msgout,1) trim(iomsg),i,trim(filename) 
+             write(msgout,1) trim(iomsg),size(d,2),i,trim(filename) 
              call error_stop(msg=trim(msgout))
           end if          
 
@@ -178,7 +178,7 @@ contains
           read (s,fmt_,iostat=ios,iomsg=iomsg) d(i, :)
 
           if (ios/=0) then 
-             write(msgout,1) trim(iomsg),i,trim(filename) 
+             write(msgout,1) trim(iomsg),size(d,2),i,trim(filename) 
              call error_stop(msg=trim(msgout))
           end if             
 
@@ -187,7 +187,7 @@ contains
 
       close(s)
 
-      1 format('loadtxt: error <',a,'> reading line ',i0,' of ',a,'.')
+      1 format('loadtxt: error <',a,'> reading ',i0,' values from line ',i0,' of ',a,'.')
 
     end subroutine loadtxt_${t1[0]}$${k1}$
   #:endfor
@@ -230,14 +230,14 @@ contains
                 iostat=ios,iomsg=iomsg) d(i, :)
 
         if (ios/=0) then 
-           write(msgout,1) trim(iomsg),i,trim(filename) 
+           write(msgout,1) trim(iomsg),size(d,2),i,trim(filename) 
            call error_stop(msg=trim(msgout))
         end if           
 
       end do
       close(s)
 
-      1 format('savetxt: error <',a,'> writing line ',i0,' of ',a,'.')
+      1 format('savetxt: error <',a,'> writing ',i0,' values to line ',i0,' of ',a,'.')
 
     end subroutine savetxt_${t1[0]}$${k1}$
   #:endfor