Solves a rectangular least-squares system of linear equations Ax≅b using singular value decomposition
With optional arguments, any of several related computations can be performed. These extra tasks include computing the rank of A, the orthogonal m×m and n×n matrices U and V, and the m×n diagonal matrix of singular values, S.
Required Arguments
A — Array of size m×n containing the matrix. (Input [/Output]) If the packaged option lin_sol_svd_overwrite_input is used, this array is not saved on output.
B — Array of size m×nb containing the right-hand side matrix. (Input [/Output] If the packaged option lin_sol_svd_overwrite_input is used, this array is not saved on output.
X— Array of size n×nb containing the solution matrix. (Output)
Optional Arguments
MROWS= m (Input) Uses array A(1:m, 1:n) for the input matrix. Default: m = size (A, 1)
NCOLS= n (Input) Uses array A(1:m, 1:n) for the input matrix. Default: n = size(A, 2)
NRHS= nb (Input) Uses the array b(1:, 1:nb) for the input right-hand side matrix. Default: nb = size(b, 2) Note that b must be a rank-2 array.
RANK= k (Output) Number of singular values that are at least as large as the value Small. It will satisfy k <= min(m, n).
u= u(:,:) (Output) Array of the same type and kind as A(1:m, 1:n). It contains the m×m orthogonal matrix U of the singular value decomposition.
s= s(:) (Output) Array of the same precision as A(1:m, 1:n). This array is real even when the matrix data is complex. It contains the m×n diagonal matrix S in a rank-1 array. The singular values are nonnegative and ordered non-increasing.
v= v(:,:) (Output) Array of the same type and kind as A(1:m, 1:n). It contains the n×n orthogonal matrix V.
iopt= iopt(:) (Input) Derived type array with the same precision as the input matrix. Used for passing optional data to the routine. The options are as follows:
Packaged Options for lin_sol_svd
Option Prefix = ?
Option Name
Option Value
s_, d_, c_, z_
lin_sol_svd_set_small
1
s_, d_, c_, z_
lin_sol_svd_overwrite_input
2
s_, d_, c_, z_
lin_sol_svd_safe_reciprocal
3
s_, d_, c_, z_
lin_sol_svd_scan_for_NaN
4
iopt(IO)= ?_options(?_lin_sol_svd_set_small, Small) Replaces with zero a diagonal term of the matrix S if it is smaller in magnitude than the value Small. This determines the approximate rank of the matrix, which is returned as the “rank=” optional argument. A solution is approximated based on this replacement. Default: the smallest number that can be safely reciprocated
iopt(IO)= ?_options(?_lin_sol_svd_overwrite_input, ?_dummy) Does not save the input arrays A(:,:) and b(:,:).
iopt(IO)= ?_options(?_lin_sol_svd_safe_reciprocal, safe) Replaces a denominator term with safe if it is smaller in magnitude than the value safe. Default: the smallest number that can be safely reciprocated
iopt(IO)= ?_options(?_lin_sol_svd_scan_for_NaN, ?_dummy) Examines each input array entry to find the first value such that
isNaN(a(i,j)) .or. isNan(b(i,j)) ==.true.
See the isNaN() function, Chapter 10. Default: Does not scan for NaNs
FORTRAN 90 Interface
Generic: CALLLIN_SOL_SVD (A, B, X[, …])
Specific: The specific interface names are S_LIN_SOL_SVD, D_LIN_SOL_SVD, C_LIN_SOL_SVD, and Z_LIN_SOL_SVD.
Description
Routine LIN_SOL_SVD solves a rectangular system of linear algebraic equations in a least-squares sense. It computes the factorization of A known as the singular value decomposition. This decomposition has the following form:
A = USVT
The matrices U and V are orthogonal. The matrix S is diagonal with the diagonal terms non-increasing. See Golub and Van Loan (1989, Chapters 5.4 and 5.5) for further details.
Fatal, Terminal, and Warning Error Messages
See the messages.gls file for error messages for LIN_SOL_SVD. These error messages are numbered 401–412; 421–432; 441–452; 461–472.
Examples
Example 1: Least-squares solution of a Rectangular System
The least-squares solution of a rectangular m×n system Ax≅b is obtained. The use of lin_sol_lsq is more efficient in this case since the matrix is of full rank. This example anticipates a problem where the matrix A is poorly conditioned or not of full rank; thus, lin_sol_svd is the appropriate routine. Also, see operator_ex13, in Chapter 10.
write (*,*) 'Example 1 for LIN_SOL_SVD is correct.'
end if
end
Output
Example 1 for LIN_SOL_SVD is correct.
Example 2: Polar Decomposition of a Square Matrix
A polar decomposition of an n×n random matrix is obtained. This decomposition satisfies A = PQ, where P is orthogonal and Q is self-adjoint and positive definite.
Given the singular value decomposition
the polar decomposition follows from the matrix products
This example uses the optional arguments “u=”, “s=”, and “v=”, then array intrinsic functions to calculate P and Q. Also, see operator_ex14, in Chapter 10.
if (sum(abs(matmul(p,transpose(p)) - ident))/sum(abs(p)) &
<= sqrt(epsilon(one))) then
if (sum(abs(a - matmul(p,q)))/sum(abs(a)) &
<= sqrt(epsilon(one))) then
write (*,*) 'Example 2 for LIN_SOL_SVD is correct.'
end if
end if
end
Output
Example 2 for LIN_SOL_SVD is correct.
Example 3: Reduction of an Array of Black and White
An n×n array A contains entries that are either 0 or 1. The entry is chosen so that as a two-dimensional object with origin at the point (1, 1), the array appears as a black circle of radius n/4 centered at the point (n/2, n/2).
A singular value decomposition
is computed, where S is of low rank. Approximations using fewer of these nonzero singular values and vectors suffice to reconstruct A. Also, see operator_ex15, supplied with the product examples.
if ((i-n/2)**2 + (j-n/2)**2 <= (n/4)**2) a(i,j) = one
end do
end do
! Compute the singular value decomposition.
call lin_sol_svd(a, b, x, nrhs=0,&
s=s, u=u, v=v)
! How many terms, to the nearest integer, exactly
! match the circle?
c = zero; k = count(s > half)
do i=1, k
c = c + spread(u(1:n,i),2,n)*spread(v(1:n,i),1,n)*s(i)
if (count(int(c-a) /= 0) == 0) exit
end do
if (i < k) then
write (*,*) 'Example 3 for LIN_SOL_SVD is correct.'
end if
end
Output
Example 3 for LIN_SOL_SVD is correct.
Example 4: Laplace Transform Solution
This example illustrates the solution of a linear least-squares system where the matrix is poorly conditioned. The problem comes from solving the integral equation:
The unknown function f(t) = 1 is computed. This problem is equivalent to the numerical inversion of the Laplace Transform of the function g(s) using real values of t and s, solving for a function that is nonzero only on the unit interval. The evaluation of the integral uses the following approximate integration rule:
The points are chosen equally spaced by using the following:
The points are computed so that the range of g(s) is uniformly sampled. This requires the solution of m equations
for j = 1, …, n and i = 1, …, m. Fortran 90 array operations are used to solve for the collocation points as a single series of steps. Newton's method,
is applied to the array function
where the following is true:
Note the coefficient matrix for the solution values
whose entry at the intersection of row i and column j is equal to the value
is explicitly integrated and evaluated as an array operation. The solution analysis of the resulting linear least-squares system
is obtained by computing the singular value decomposition
An approximate solution is computed with the transformed right-hand side
followed by using as few of the largest singular values as possible to minimize the following squared error residual:
This determines an optimal value k to use in the approximate solution
Also, see operator_ex16, supplied with the product examples.
