FNLMath : Linear Systems : LIN_SOL_SVD
LIN_SOL_SVD
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: CALL LIN_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.
 
use lin_sol_svd_int
use rand_gen_int
implicit none
! This is Example 1 for LIN_SOL_SVD.
integer, parameter :: m=128, n=32
real(kind(1d0)), parameter :: one=1d0
real(kind(1d0)) A(m,n), b(m,1), x(n,1), y(m*n), err
! Generate a random matrix and right-hand side.
call rand_gen(y)
A = reshape(y,(/m,n/))
call rand_gen(b(1:m,1))
! Compute the least-squares solution matrix of Ax=b.
call lin_sol_svd(A, b, x)
! Check that the residuals are orthogonal to the
! column vectors of A.
err = sum(abs(matmul(transpose(A),b-matmul(A,x))))/sum(abs(A))
if (err <= sqrt(epsilon(one))) then
 
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.
 
use lin_sol_svd_int
use rand_gen_int
implicit none
! This is Example 2 for LIN_SOL_SVD.
integer i
integer, parameter :: n=32
real(kind(1d0)), parameter :: one=1.0d0, zero=0.0d0
real(kind(1d0)) a(n,n), b(n,0), ident(n,n), p(n,n), q(n,n), &
s_d(n), u_d(n,n), v_d(n,n), x(n,0), y(n*n)
! Generate a random matrix.
call rand_gen(y)
a = reshape(y,(/n,n/))
! Compute the singular value decomposition.
call lin_sol_svd(a, b, x, nrhs=0, s=s_d, &
u=u_d, v=v_d)
! Compute the (left) orthogonal factor.
p = matmul(u_d,transpose(v_d))
! Compute the (right) self-adjoint factor.
q = matmul(v_d*spread(s_d,1,n),transpose(v_d))
ident=zero
do i=1, n
ident(i,i) = one
end do
! Check the results.
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.
use lin_sol_svd_int
use rand_gen_int
use error_option_packet
implicit none
! This is Example 3 for LIN_SOL_SVD.
integer i, j, k
integer, parameter :: n=32
real(kind(1e0)), parameter :: half=0.5e0, one=1e0, zero=0e0
real(kind(1e0)) a(n,n), b(n,0), x(n,0), s(n), u(n,n), &
v(n,n), c(n,n)
! Fill in value one for points inside the circle.
a = zero
do i=1, n
do j=1, n
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.
use lin_sol_svd_int
use rand_gen_int
use error_option_packet
implicit none
! This is Example 4 for LIN_SOL_SVD.
integer i, j, k
integer, parameter :: m=64, n=16
real(kind(1e0)), parameter :: one=1e0, zero=0.0e0
real(kind(1e0)) :: g(m), s(m), t(n+1), a(m,n), b(m,1), &
f(n,1), U_S(m,m), V_S(n,n), S_S(n), &
rms, oldrms
real(kind(1e0)) :: delta_g, delta_t
delta_g = one/real(m+1,kind(one))
! Compute which collocation equations to solve.
do i=1,m
g(i)=i*delta_g
end do
! Compute equally spaced quadrature points.
delta_t =one/real(n,kind(one))
do j=1,n+1
t(j)=(j-1)*delta_t
end do
! Compute collocation points.
s=m
solve_equations: do
s=s-(exp(-s)-(one-s*g))/(g-exp(-s))
if (sum(abs((one-exp(-s))/s - g)) <= &
epsilon(one)*sum(g)) &
exit solve_equations
end do solve_equations
! Evaluate the integrals over the quadrature points.
a = (exp(-spread(t(1:n),1,m)*spread(s,2,n)) &
- exp(-spread(t(2:n+1),1,m)*spread(s,2,n))) / &
spread(s,2,n)
b(1:,1)=g
! Compute the singular value decomposition.
call lin_sol_svd(a, b, f, nrhs=0, &
rank=k, u=U_S, v=V_S, s=S_S)
! Singular values that are larger than epsilon determine
! the rank=k.
k = count(S_S > epsilon(one))
oldrms = huge(one)
g = matmul(transpose(U_S), b(1:m,1))
! Find the minimum number of singular values that gives a good
! approximation to f(t) = 1.
do i=1,k
f(1:n,1) = matmul(V_S(1:,1:i), g(1:i)/S_S(1:i))
f = f - one
rms = sum(f**2)/n
if (rms > oldrms) exit
oldrms = rms
end do
write (*,"( ' Using this number of singular values, ', &
&i4 / ' the approximate R.M.S. error is ', 1pe12.4)") &
i-1, oldrms
if (sqrt(oldrms) <= delta_t**2) then
write (*,*) 'Example 4 for LIN_SOL_SVD is correct.'
end if
end
Output
 
Example 4 for LIN_SOL_SVD is correct.