Module_LLSQ_Procedures¶
Copyright 2018 IRD
This file is part of statpack.
statpack is free software: you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
statpack is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details.
You can find a copy of the GNU Lesser General Public License in the statpack/doc directory.
MODULE EXPORTING SUBROUTINES AND FUNCTIONS FOR SOLVING LINEAR LEAST SQUARES PROBLEMS.
LATEST REVISION : 25/09/2018
function solve_llsq ( a, b, krank, tol, min_norm )¶
Purpose¶
SOLVE_LLSQ computes a solution to a real linear least squares problem:
Minimize 2-norm || B - A * X ||using an orthogonal factorization with columns pivoting of A. A is an m-by-n matrix which may be rank-deficient. m>=n or n>m is permitted. Here, B is a m right hand side vector and X is a n solution vector.
The function returns the n solution vector X.
A and B are not overwritten by SOLVE_LLSQ.
Arguments¶
- A (INPUT) real(stnd), dimension(:,:)
- On entry, the m-by-n coefficient matrix A.
- B (INPUT) real(stnd), dimension(:)
On entry, the m right hand side vector B.
- The shape of B must verify:
- size( B ) = size( A, 1 ) = m .
- KRANK (INPUT, OPTIONAL) integer(i4b)
- On entry, KRANK=k implies that the first k columns of A are to be forced into the basis, pivoting is performed on the last n-k columns of A. When KRANK >=min(m,n) is used, pivoting is not performed. This is appropriate when A is known to be full rank. If KRANK is absent or is <=0, pivoting is done on all columns of A. This is appropriate if A may not be of full rank (i.e. certain columns of A are linear combinations of other columns).
- TOL (INPUT, OPTIONAL) real(stnd)
On entry, TOL is used to determine the effective rank of A, which is then defined as the order of the largest leading triangular submatrix R11 in the QR factorization (with pivoting) of A whose estimated condition number, in the 1-norm, is less than 1/TOL. TOl must be set to the relative precision of the elements in A and B. If each element is correct to, say, 5 digits then TOL=0.00001 should be used.
TOL must not be greater or equal to 1 or less or equal than 0, otherwise the numerical rank of A is determined and the calculations to determine the condition number are not performed. If TOL is absent, the numerical rank of A is determined.
- MIN_NORM (INPUT, OPTIONAL) logical(lgl)
- On entry, If MIN_NORM=true, the minimun 2-norm solution is computed. If MIN_NORM=false or if MIN_NORM is absent, a solution is computed such that if the j-th column of A is omitted from the basis, X[j] is set to zero.
Further Details¶
The routine first computes a QR factorization with (partial) column pivoting on option (see below):
A * P = Q * R
, here P is n-by-n permutation matrix, R is an upper triangular or trapezoidal (if n>m) matrix and Q is a m-by-m orthogonal matrix.
R can then be partioned by defining R11 as the largest leading submatrix of R whose estimated condition number, in the 1-norm, is less than 1/TOL (or such that abs(R[j,j])>0 if TOL is absent). The order of R11, arank, is the effective rank of A.
This leads to the following partition of R:
[ R11 R12 ]
[ R21 R22 ]
where R21 is zero by construction (since R is an upper triangular or trapezoidal) and R22 is considered to be negligible.
If MIN_NORM is present and has the value true, R12 is annihilated by orthogonal transformations from the right, arriving at the complete orthogonal factorization:
A * P = Q * T * Z
, where Z is a n-by-n orthogonal matrix and T has the form:
[ T11 T12 ]
[ T21 T22 ]
, here T21 (=R21), T12 and T22 (=R22) are zero and T11 is a arank-by-arank upper triangular matrix.
The minimum 2-norm solution is then
X = [ P * Z’ ](:,:arank) * [ inv(T11) * Q1’ * B ]
where inv(T11) is the inverse of T11, Z’ is the transpose of Z and Q1 consists of the first arank columns of Q.
If MIN_NORM is absent or has the value false, a solution is computed as
X = P(:,:arank) * [ inv(R11) * Q1’ * B ]
where inv(R11) is the inverse of R11 and Q1 consists of the first arank columns of Q. In this case, if the j-th column of A is omitted from the basis, X[j] is set to zero.
On input, if KRANK is present and KRANK=k, the first k columns of A are to be forced into the basis. Pivoting is performed on the last n-k columns of A.
When KRANK is present and KRANK>=min(m,n) is used, pivoting is not performed. This is appropriate when A is known to be full rank.
If KRANK is absent or is present with KRANK<=0, pivoting is done on all columns of A.
TOL is an optional argument such that 0<TOL<1. If TOL is not specified, or is outside ]0,1[, the calculations to determine the condition number of A are not performed and crude tests on R(j,j) are performed to determine the numerical rank of A. If TOL is present and is in ]0,1[, the calculations to determine the condition number are performed.
If it is possible that A may not be full rank (i.e., certain columns of A are linear combinations of other columns), then the linearly dependent columns can usually be eliminated by using KRANK=0 and TOL=relative precision of the elements in A and B. If each element is correct to, say, 5 digits then TOL=0.00001 should be used. Also, it may be helpful to scale the columns of A so that all elements are about the same order of magnitude.
On exit, if A or B are empty, the function returns a n-vector filled with nan() value.
function solve_llsq ( a, b, krank, tol, min_norm )¶
Purpose¶
SOLVE_LLSQ computes solutions to real linear least squares problems of the form:
Minimize 2-norm|| B - A * X ||using an orthogonal factorization with columns pivoting of A. A is an m-by-n matrix which may be rank-deficient. m>=n or n>m is permitted.
Several right hand side vectors b can be handled in a single call; they are stored as the columns of the m-by-nb right hand side matrix B.
The function returns the n-by-nb solution matrix X.
A and B are not overwritten by SOLVE_LLSQ.
Arguments¶
- A (INPUT) real(stnd), dimension(:,:)
- On entry, the m-by-n coefficient matrix A.
- B (INPUT) real(stnd), dimension(:,:)
On entry, the m-by-nb right hand side matrix B.
- The shape of B must verify:
- size( B, 1 ) = size( A, 1 ) = m .
- KRANK (INPUT, OPTIONAL) integer(i4b)
- On entry, KRANK=k implies that the first k columns of A are to be forced into the basis, pivoting is performed on the last n-k columns of A. When KRANK >=min(m,n) is used, pivoting is not performed. This is appropriate when A is known to be full rank. If KRANK is absent or is <=0, pivoting is done on all columns of A. This is appropriate if A may not be of full rank (i.e. certain columns of A are linear combinations of other columns).
- TOL (INPUT, OPTIONAL) real(stnd)
On entry, TOL is used to determine the effective rank of A, which is then defined as the order of the largest leading triangular submatrix R11 in the QR factorization (with pivoting) of A whose estimated condition number, in the 1-norm, is less than 1/TOL. TOl must be set to the relative precision of the elements in A and B. If each element is correct to, say, 5 digits then TOL=0.00001 should be used.
TOL must not be greater or equal to 1 or less or equal than 0, otherwise the numerical rank of A is determined and the calculations to determine the condition number are not performed. If TOL is absent, the numerical rank of A is determined.
- MIN_NORM (INPUT, OPTIONAL) logical(lgl)
- On entry, If MIN_NORM=true, minimun 2-norm solutions are computed. If MIN_NORM=false or if MIN_NORM is absent, solutions are computed such that if the j-th column of A is omitted from the basis, X[j,:] is set to zero.
Further Details¶
The routine first computes a QR factorization with (partial) column pivoting on option (see below):
A * P = Q * R
, here P is n-by-n permutation matrix, R is an upper triangular or trapezoidal (if n>m) matrix and Q is a m-by-m orthogonal matrix.
R can then be partioned by defining R11 as the largest leading submatrix of R whose estimated condition number, in the 1-norm, is less than 1/TOL (or such that abs(R[j,j])>0 if TOL is absent). The order of R11, arank, is the effective rank of A.
This leads to the following partition of R:
[ R11 R12 ]
[ R21 R22 ]
where R21 is zero by construction (since R is an upper triangular or trapezoidal) and R22 is considered to be negligible.
If MIN_NORM is present and has the value true, R12 is annihilated by orthogonal transformations from the right, arriving at the complete orthogonal factorization:
A * P = Q * T * Z
, where Z is a n-by-n orthogonal matrix and T has the form:
[ T11 T12 ]
[ T21 T22 ]
, here T21 (=R21), T12 and T22 (=R22) are zero and T11 is a arank-by-arank upper triangular matrix.
The minimum 2-norm solution is then
X = [ P * Z’ ](:,:arank) * [ inv(T11) * Q1’ * B ]
where inv(T11) is the inverse of T11, Z’ is the transpose of Z and Q1 consists of the first arank columns of Q.
If MIN_NORM is absent or has the value false, a solution is computed as
X = P(:,:arank) * [ inv(R11) * Q1’ * B ]
where inv(R11) is the inverse of R11 and Q1 consists of the first arank columns of Q. In this case, if the j-th column of A is omitted from the basis, X[j] is set to zero.
On input, if KRANK is present and KRANK=k, the first k columns of A are to be forced into the basis. Pivoting is performed on the last n-k columns of A.
When KRANK is present and KRANK>=min(m,n) is used, pivoting is not performed. This is appropriate when A is known to be full rank.
If KRANK is absent or is present with KRANK<=0, pivoting is done on all columns of A.
TOL is an optional argument such that 0<TOL<1. If TOL is not specified, or is outside ]0,1[, the calculations to determine the condition number of A are not performed and crude tests on R(j,j) are performed to determine the numerical rank of A. If TOL is present and is in ]0,1[, the calculations to determine the condition number are performed.
If it is possible that A may not be full rank (i.e., certain columns of A are linear combinations of other columns), then the linearly dependent columns can usually be eliminated by using KRANK=0 and TOL=relative precision of the elements in A and B. If each element is correct to, say, 5 digits then TOL=0.00001 should be used. Also, it may be helpful to scale the columns of A so that all elements are about the same order of magnitude.
On exit, if A or B are empty, the function returns a n-by-nb matrix filled with nan() value.
function solve_llsq ( a, b )¶
Purpose¶
SOLVE_LLSQ computes a solution to a real linear least squares problem:
Minimize 2-norm|| B - A * X ||A is an m vector, B is a m right hand side vector and X is a real scalar.
The function returns the solution scalar X.
A and B are not overwritten by SOLVE_LLSQ.
Arguments¶
- A (INPUT) real(stnd), dimension(:)
- On entry, the m coefficient vector A.
- B (INPUT) real(stnd), dimension(:)
On entry, the m right hand side vector B.
- The shape of B must verify:
- size( B ) = size( A ) = m .
Further Details¶
The routine first generates a real elementary reflector H of order m, such that
H * A = D , with H’ * H = I and D’ = ( d 0 )
where d is a scalar. H is represented in the form
H = I + beta * ( v * v’ ) ,
where beta is a real scalar and v is a real m-element vector.
The solution X is then computed as
X = [ H * B ](1) / d
On exit, if A or B are empty, the function returns a nan() value.
function solve_llsq ( a, b )¶
Purpose¶
SOLVE_LLSQ computes solutions to real linear least squares problems of the form:
Minimize 2-norm|| B - A * X ||A is an m vector and several right hand side vectors b can be handled in a single call; they are stored as the columns of the m-by-nb right hand side matrix B.
The function returns the nb solution vector X.
A and B are not overwritten by SOLVE_LLSQ.
Arguments¶
- A (INPUT) real(stnd), dimension(:)
- On entry, the m coefficient vector A.
- B (INPUT) real(stnd), dimension(:,:)
On entry, the m-by-nb right hand side matrix B.
- The shape of B must verify:
- size( B, 1 ) = size( A ) = m .
Further Details¶
The routine first generates a real elementary reflector H of order m, such that
H * A = D , with H’ * H = I and D’ = ( d 0 )
where d is a scalar. H is represented in the form
H = I + beta * ( v * v’ ) ,
where beta is a real scalar and v is a real m-element vector.
The solution vector X is then computed as
X(:) = [ H * B ](1,:) / d
On exit, if A or B are empty, the function returns a nb-vector filled with nan() value.
subroutine llsq_qr_solve ( mat, b, x, rnorm, resid, krank, tol, min_norm )¶
Purpose¶
LLSQ_QR_SOLVE computes a solution to a real linear least squares problem:
Minimize 2-norm|| B - MAT * X ||using an orthogonal factorization with columns pivoting of MAT. MAT is an m-by-n matrix which may be rank-deficient. m>=n or n>m is permitted. Here, B is a m right hand side vector and X is a n solution vector.
MAT and B are not overwritten by LLSQ_QR_SOLVE.
Arguments¶
- MAT (INPUT) real(stnd), dimension(:,:)
- On entry, the m-by-n coefficient matrix MAT.
- B (INPUT) real(stnd), dimension(:)
On entry, the m right hand side vector B.
- The shape of B must verify:
- size( B ) = size( MAT, 1 ) = m .
- X (OUTPUT) real(stnd), dimension(:)
On exit, the n solution vector X.
- The shape of X must verify:
- size( X ) = size( MAT, 2 ) = n .
- RNORM (OUTPUT, OPTIONAL) real(stnd)
- On exit, the 2-norm of the residual vector for the solution vector X.
- RESID (OUTPUT, OPTIONAL) real(stnd), dimension(:)
On exit, the m residual vector for the solution vector X, RESID = B - MAT * X.
- The shape of RESID must verify:
- size( RESID ) = size( B ) = size( MAT, 1 ) = m .
- KRANK (INPUT/OUTPUT, OPTIONAL) integer(i4b)
On entry, KRANK=k implies that the first k columns of MAT are to be forced into the basis, pivoting is performed on the last n-k columns of MAT. When KRANK >=min(m,n) is used, pivoting is not performed. This is appropriate when MAT is known to be full rank. If KRANK is absent or is <=0, pivoting is done on all columns of MAT. This is appropriate if MAT may not be of full rank (i.e. certain columns of MAT are linear combinations of other columns).
On exit, KRANK contains the effective rank of MAT, i.e., the number of independent columns in matrix MAT.
- TOL (INPUT/OUTPUT, OPTIONAL) real(stnd)
On entry, TOL is used to determine the effective rank of MAT, which is then defined as the order of the largest leading triangular submatrix R11 in the QR factorization (with pivoting) of MAT whose estimated condition number, in the 1-norm, is less than 1/TOL. TOl must be set to the relative precision of the elements in MAT and B. If each element is correct to, say, 5 digits then TOL=0.00001 should be used.
TOL must not be greater or equal to 1 or less than 0, otherwise the numerical rank of MAT is determined and the calculations to determine the condition number are not performed. If TOL=0, the numerical rank of MAT is determined, but the condition number is calculated.
On exit, if a condition number is calculated, its reciprocal is returned in TOL. Otherwise, TOL is not changed.
If TOL is absent, the numerical rank of MAT is used and is returned in the optional argument KRANK.
- MIN_NORM (INPUT, OPTIONAL) logical(lgl)
- On entry, If MIN_NORM=true, the minimun 2-norm solution is computed. If MIN_NORM=false or if MIN_NORM is absent, a solution is computed such that if the j-th column of MAT is omitted from the basis, X[j] is set to zero.
Further Details¶
The routine first computes a QR factorization with (partial) column pivoting on option (see below):
MAT * P = Q * R
, here P is n-by-n permutation matrix, R is an upper triangular or trapezoidal (if n>m) matrix and Q is a m-by-m orthogonal matrix.
R can then be partioned by defining R11 as the largest leading submatrix of R whose estimated condition number, in the 1-norm, is less than 1/TOL (or such that abs(R[j,j])>0 if TOL is absent). The order of R11, KRANK, is the effective rank of MAT.
This leads to the following partition of R:
[ R11 R12 ]
[ R21 R22 ]
where R21 is zero by construction (since R is an upper triangular or trapezoidal) and R22 is considered to be negligible.
If MIN_NORM is present and has the value true, R12 is annihilated by orthogonal transformations from the right, arriving at the complete orthogonal factorization:
MAT * P = Q * T * Z
, where Z is a n-by-n orthogonal matrix and T has the form:
[ T11 T12 ]
[ T21 T22 ]
, here T21 (=R21), T12 and T22 (=R22) are zero and T11 is a KRANK-by-KRANK upper triangular matrix.
The minimum 2-norm solution is then
X = [ P * Z’ ](:,:KRANK) * [ inv(T11) * Q1’ * B ]
where inv(T11) is the inverse of T11, Z’ is the transpose of Z and Q1 consists of the first KRANK columns of Q.
If MIN_NORM is absent or has the value false, a solution is computed as
X = P(:,:KRANK) * [ inv(R11) * Q1’ * B ]
where inv(R11) is the inverse of R11 and Q1 consists of the first KRANK columns of Q. In this case, if the j-th column of MAT is omitted from the basis, X[j] is set to zero.
In both cases:
- The 2-norm of the residual vector for the solution X can be obtained through the optional argument RNORM .
- The m residual vector, B - MAT * X, can be obtained through the optional argument RESID.
On input, if KRANK is present and KRANK=k, the first k columns of MAT are to be forced into the basis. Pivoting is performed on the last n-k columns of MAT.
When KRANK is present and KRANK>=min(m,n) is used, pivoting is not performed. This is appropriate when MAT is known to be full rank.
If KRANK is absent or is present with KRANK<=0, pivoting is done on all columns of MAT.
On output, if KRANK is present, it contains the effective rank of MAT, i.e., the order of the submatrix R11.
TOL is an optional argument such that 0<=TOL<1. If TOL is not specified, or is outside [0,1[, the calculations to determine the condition number of MAT are not performed and crude tests on R(j,j) are performed to determine the numerical rank of MAT. If TOL is present and is in [0,1[, the calculations to determine the condition number are performed and its reciprocal is return in TOL.
If it is possible that MAT may not be full rank (i.e., certain columns of MAT are linear combinations of other columns), then the linearly dependent columns can usually be eliminated by using KRANK=0 and TOL=relative precision of the elements in MAT and B. If each element is correct to, say, 5 digits then TOL=0.00001 should be used. Also, it may be helpful to scale the columns of MAT so that all elements are about the same order of magnitude.
subroutine llsq_qr_solve ( mat, b, x, rnorm, resid, krank, tol, min_norm )¶
Purpose¶
LLSQ_QR_SOLVE computes solutions to real linear least squares problems of the form:
Minimize 2-norm|| B - MAT * X ||using an orthogonal factorization with columns pivoting of MAT. MAT is an m-by-n matrix which may be rank-deficient. m>=n or n>m is permitted.
Several right hand side vectors b and solution vectors x can be handled in a single call; they are stored as the columns of the m-by-nb right hand side matrix B and the n-by-nb solution matrix X, respectively.
MAT and B are not overwritten by LLSQ_QR_SOLVE.
Arguments¶
- MAT (INPUT) real(stnd), dimension(:,:)
- On entry, the m-by-n coefficient matrix MAT.
- B (INPUT) real(stnd), dimension(:,:)
On entry, the m-by-nb right hand side matrix B.
The shape of B must verify:
- size( B, 1 ) = size( MAT, 1 ) = m
- size( B, 2 ) = size( X, 2 ) = nb .
- X (OUTPUT) real(stnd), dimension(:,:)
On exit, the n-by-nb solution matrix X.
The shape of X must verify:
- size( X, 1 ) = size( MAT, 2 ) = n
- size( X, 2 ) = size( B, 2 ) = nb .
- RNORM (OUTPUT, OPTIONAL) real(stnd), dimension(:)
On exit, the 2-norm of the residual vectors for the solutions stored columnwise in the matrix X.
- The size of RNORM must verify:
- size( RNORM ) = size( X, 2 ) = size( B, 2 ) = nb .
- RESID (OUTPUT, OPTIONAL) real(stnd), dimension(:,:)
On exit, the residual vectors for the solutions stored columnwise in the matrix X, RESID = B - MAT * X.
The shape of RESID must verify:
- size( RESID, 1 ) = size( B, 1 ) = size( MAT, 1 ) = m
- size( RESID, 2 ) = size( B, 2 ) = size( X, 2 ) = nb .
- KRANK (INPUT/OUTPUT, OPTIONAL) integer(i4b)
On entry, KRANK=k implies that the first k columns of MAT are to be forced into the basis, pivoting is performed on the last n-k columns of MAT. When KRANK >=min(m,n) is used, pivoting is not performed. This is appropriate when MAT is known to be full rank. If KRANK is absent or is <=0, pivoting is done on all columns of MAT. This is appropriate if MAT may not be of full rank (i.e. certain columns of MAT are linear combinations of other columns).
On exit, KRANK contains the effective rank of MAT, i.e., the number of independent columns in matrix MAT.
- TOL (INPUT/OUTPUT, OPTIONAL) real(stnd)
On entry, TOL is used to determine the effective rank of MAT, which is then defined as the order of the largest leading triangular submatrix R11 in the QR factorization (with pivoting) of MAT whose estimated condition number, in the 1-norm, is less than 1/TOL. TOl must be set to the relative precision of the elements in MAT and B. If each element is correct to, say, 5 digits then TOL=0.00001 should be used.
TOL must not be greater or equal to 1 or less than 0, otherwise the numerical rank of MAT is determined and the calculations to determine the condition number are not performed. If TOL=0, the numerical rank of MAT is determined, but the condition number is calculated.
On exit, if a condition number is calculated, its reciprocal is returned in TOL. Otherwise, TOL is not changed.
If TOL is absent, the numerical rank of MAT is used and is returned in the optional argument KRANK.
- MIN_NORM (INPUT, OPTIONAL) logical(lgl)
- On entry, If MIN_NORM=true, minimun 2-norm solutions are computed. If MIN_NORM=false or if MIN_NORM is absent, solutions are computed such that if the j-th column of MAT is omitted from the basis, X[j,:] is set to zero.
Further Details¶
The routine first computes a QR factorization with (partial) column pivoting on option (see below):
MAT * P = Q * R
, here P is n-by-n permutation matrix, R is an upper triangular or trapezoidal (if n>m) matrix and Q is a m-by-m orthogonal matrix.
R can then be partioned by defining R11 as the largest leading submatrix of R whose estimated condition number, in the 1-norm, is less than 1/TOL (or such that abs(R[j,j])>0 if TOL is absent). The order of R11, KRANK, is the effective rank of MAT.
This leads to the following partition of R:
[ R11 R12 ]
[ R21 R22 ]
where R21 is zero by construction (since R is an upper triangular or trapezoidal) and R22 is considered to be negligible.
If MIN_NORM is present and has the value true, R12 is annihilated by orthogonal transformations from the right, arriving at the complete orthogonal factorization:
MAT * P = Q * T * Z
, where Z is a n-by-n orthogonal matrix and T has the form:
[ T11 T12 ]
[ T21 T22 ]
, here T21 (=R21), T12 and T22 (=R22) are zero and T11 is a KRANK-by-KRANK upper triangular matrix.
The minimum 2-norm solution is then
X = [ P * Z’ ](:,:KRANK) * [ inv(T11) * Q1’ * B ]
where inv(T11) is the inverse of T11, Z’ is the transpose of Z and Q1 consists of the first KRANK columns of Q.
If MIN_NORM is absent or has the value false, a solution is computed as
X = P(:,:KRANK) * [ inv(R11) * Q1’ * B ]
where inv(R11) is the inverse of R11 and Q1 consists of the first KRANK columns of Q. In this case, if the j-th column of MAT is omitted from the basis, X[j] is set to zero.
In both cases:
- The 2-norm of the residual vector for the solution in the j-th column of X is given in RNORM[j] if argument RNORM is present.
- The residual matrix, B - MAT * X, can be obtained through the optional argument RESID.
On input, if KRANK is present and KRANK=k, the first k columns of MAT are to be forced into the basis. Pivoting is performed on the last n-k columns of MAT.
When KRANK is present and KRANK>=min(m,n) is used, pivoting is not performed. This is appropriate when MAT is known to be full rank.
If KRANK is absent or is present with KRANK<=0, pivoting is done on all columns of MAT.
On output, if KRANK is present, it contains the effective rank of MAT, i.e., the order of the submatrix R11.
TOL is an optional argument such that 0<=TOL<1. If TOL is not specified, or is outside [0,1[, the calculations to determine the condition number of MAT are not performed and crude tests on R(j,j) are performed to determine the numerical rank of MAT. If TOL is present and is in [0,1[, the calculations to determine the condition number are performed and its reciprocal is return in TOL.
If it is possible that MAT may not be full rank (i.e., certain columns of MAT are linear combinations of other columns), then the linearly dependent columns can usually be eliminated by using KRANK=0 and TOL=relative precision of the elements in MAT and B. If each element is correct to, say, 5 digits then TOL=0.00001 should be used. Also, it may be helpful to scale the columns of MAT so that all elements are about the same order of magnitude.
subroutine llsq_qr_solve ( vec, b, x, rnorm, resid )¶
Purpose¶
LLSQ_QR_SOLVE computes a solution to a real linear least squares problem:
Minimize 2-norm|| B - VEC * X ||VEC is an m vector, B is a m right hand side vector and X is a real scalar.
VEC and B are not overwritten by LLSQ_QR_SOLVE.
Arguments¶
- VEC (INPUT) real(stnd), dimension(:)
- On entry, the m coefficient vector VEC.
- B (INPUT) real(stnd), dimension(:)
On entry, the m right hand side vector B.
- The shape of B must verify:
- size( B ) = size( VEC ) = m .
- X (OUTPUT) real(stnd)
- On exit, the real solution X.
- RNORM (OUTPUT, OPTIONAL) real(stnd)
- On exit, the 2-norm of the residual vector for the solution scalar X.
- RESID (OUTPUT, OPTIONAL) real(stnd), dimension(:)
On exit, the m residual vector for the solution X, RESID = B - MAT * X.
- The shape of RESID must verify:
- size( RESID ) = size( B ) = size( VEC ) = m .
Further Details¶
The routine first generates a real elementary reflector H of order m, such that
H * VEC = D , with H’ * H = I and D’ = ( d 0 )
where d is a scalar. H is represented in the form
H = I + beta * ( v * v’ ) ,
where beta is a real scalar and v is a real m-element vector.
The solution X is then computed as
X = [ H * B ](1) / d
The 2-norm of the residual vector for the solution X can be obtained through the optional argument RNORM as
2-norm|| [ H * B ](2:) ||
The residual vector, B - VEC * X, can be obtained through the optional argument RESID.
subroutine llsq_qr_solve ( vec, b, x, rnorm, resid )¶
Purpose¶
LLSQ_QR_SOLVE computes solutions to real linear least squares problems of the form:
Minimize 2-norm|| B - VEC * X ||VEC is an m vector and several right hand side vectors b and solution scalars x can be handled in a single call; they are stored as the columns of the m-by-nb right hand side matrix B and the nb solution vector X, respectively.
VEC and B are not overwritten by LLSQ_QR_SOLVE.
Arguments¶
- VEC (INPUT) real(stnd), dimension(:)
- On entry, the m coefficient vector VEC.
- B (INPUT) real(stnd), dimension(:,:)
On entry, the m-by-nb right hand side matrix B.
The shape of B must verify:
- size( B, 1 ) = size( VEC ) = m
- size( B, 2 ) = size( X ) = nb .
- X (OUTPUT) real(stnd), dimension(:)
On exit, the nb solution vector X.
- The shape of X must verify:
- size( X ) = size( B, 2 ) = nb .
- RNORM (OUTPUT, OPTIONAL) real(stnd), dimension(:)
On exit, the 2-norms of the residual vectors for the solutions stored in the vector X.
- The size of RNORM must verify:
- size( RNORM ) = size( X ) = size( B, 2 ) = nb .
- RESID (OUTPUT, OPTIONAL) real(stnd), dimension(:,:)
On exit, the residual vectors for the solutions stored in the vector X, RESID = B - VEC * X.
The shape of RESID must verify:
- size( RESID, 1 ) = size( B, 1 ) = size( VEC ) = m
- size( RESID, 2 ) = size( B, 2 ) = size( X ) = nb .
Further Details¶
The routine first generates a real elementary reflector H of order m, such that
H * A = D , with H’ * H = I and D’ = ( d 0 )
where d is a scalar. H is represented in the form
H = I + beta * ( v * v’ ) ,
where beta is a real scalar and v is a real m-element vector.
The solution vector X is then computed as
X(:) = [ H * B ](1,:) / d
The 2-norm of the residual vector for the solution X[j] is given in RNORM[j] if argument RNORM is present.
The residual matrix, B - VEC * X, can be obtained through the optional argument RESID.
subroutine llsq_qr_solve2 ( mat, b, x, rnorm, comp_resid, krank, tol, min_norm, diagr, beta, ip, tau )¶
Purpose¶
LLSQ_QR_SOLVE2 computes a solution to a real linear least squares problem:
Minimize 2-norm|| B - MAT * X ||using a (complete) orthogonal factorization of MAT. MAT is an m-by-n matrix which may be rank-deficient. m>=n or n>m is permitted. Here, B is a m right hand side vector and X is a n solution vector.
MAT and B are overwritten with information generated by LLSQ_QR_SOLVE2. The (complete) orthogonal factorization of MAT is saved in arguments MAT, DIAGR, BETA, IP and TAU on output.
Arguments¶
- MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)
On entry, the m-by-n coefficient matrix MAT.
On exit, MAT has been overwritten by details of its (complete) orthogonal factorization. Other parts of the factorization can be obtained if the optional arguments DIAGR, BETA, IP and TAU are present. See Further Details.
- B (INPUT/OUTPUT) real(stnd), dimension(:)
On entry, the m right hand side vector B.
On exit, if COMP_RESID is present and is equal true, the residual vector B - MAT * X overwrites B on output.
- The shape of B must verify:
- size( B ) = size( MAT, 1 ) = m .
- X (OUTPUT) real(stnd), dimension(:)
On exit, the n solution vector X.
- The shape of X must verify:
- size( X ) = size( MAT, 2 ) = n .
- RNORM (OUTPUT, OPTIONAL) real(stnd)
- On exit, the 2-norm of the residual vector for the solution vector X.
- COMP_RESID (INPUT, OPTIONAL) logical(lgl)
- On entry, if COMP_RESID is present and is equal true, the residual vector B - MAT * X overwrites B on exit.
- KRANK (INPUT/OUTPUT, OPTIONAL) integer(i4b)
On entry, KRANK=k implies that the first k columns of MAT are to be forced into the basis, pivoting is performed on the last n-k columns of MAT. When KRANK >=min(m,n) is used, pivoting is not performed. This is appropriate when MAT is known to be full rank. If KRANK is absent or is <=0, pivoting is done on all columns of MAT. This is appropriate if MAT may not be of full rank (i.e. certain columns of MAT are linear combinations of other columns).
On exit, KRANK contains the effective rank of MAT, i.e., the number of independent columns in matrix MAT.
- TOL (INPUT/OUTPUT, OPTIONAL) real(stnd)
On entry, TOL is used to determine the effective rank of MAT, which is then defined as the order of the largest leading triangular submatrix R11 in the QR factorization (with pivoting) of MAT whose estimated condition number, in the 1-norm, is less than 1/TOL. TOL must be set to the relative precision of the elements in MAT and B. If each element is correct to, say, 5 digits then TOL=0.00001 should be used.
TOL must not be greater or equal to 1 or less than 0, otherwise the numerical rank of MAT is determined and the calculations to determine the condition number are not performed. If TOL=0, the numerical rank of MAT is determined, but the condition number is calculated.
On exit, if a condition number is calculated, its reciprocal is returned in TOL. Otherwise, TOL is not changed.
If TOL is absent, the numerical rank of MAT is used and is returned in the optional argument KRANK.
- MIN_NORM (INPUT, OPTIONAL) logical(lgl)
- On entry: If MIN_NORM=true, a complete orthogonal factorization of MAT and the minimun 2-norm solution are computed. If MIN_NORM=false or if MIN_NORM is absent, a QR factorization with column pivoting of MAT and a solution are computed such that if the j-th column of MAT is omitted from the basis, X[j] is set to zero.
- DIAGR (OUTPUT, OPTIONAL) real(stnd), dimension(:)
On exit, the diagonal elements of the matrix R if MIN_NORM=false or is absent, or the diagonal elements of the matrix T11 if MIN_NORM is present and is equal to true. The diagonal elements of T11 are stored in DIAGR(1:KRANK). See Further Details.
The size of DIAGR must be min( size(MAT,1) , size(MAT,2) ).
- BETA (OUTPUT, OPTIONAL) real(stnd), dimension(:)
On exit, the scalar factors of the elementary reflectors defining Q. See Further Details.
The size of BETA must be min( size(MAT,1) , size(MAT,2) ).
- IP (OUTPUT, OPTIONAL) integer(i4b), dimension(:)
On exit, if IP(j)=k, then the j-th column of MAT*P was the k-th column of MAT. See Further Details.
The size of IP must be size(MAT,2) = n.
- TAU (OUTPUT, OPTIONAL) real(stnd), dimension(:)
On exit, the scalars factors of the elementary reflectors defining Z in the complete orthogonal factorization of MAT if MIN_NORM is present and is equal to true. Otherwise, TAU is set to 0. See Further Details.
The size of TAU must be min( size(MAT,1) , size(MAT,2) ).
Further Details¶
The routine first computes a QR factorization with (partial) column pivoting on option (see below):
MAT * P = Q * R
, here P is n-by-n permutation matrix, R is an upper triangular or trapezoidal (if n>m) matrix and Q is a m-by-m orthogonal matrix.
The matrix Q is represented as a product of elementary reflectors
Q = H(1) * H(2) * … * H(k), where k = min( size(MAT,1) , size(MAT,2) )
Each H(i) has the form
H(i) = I + beta * ( v * v’ ) ,
where beta is a real scalar and v is a real m-element vector with v(1:i-1) = 0. v(i:m) is stored on exit in MAT(i:m,i) and beta in BETA(i).
The matrix P is represented in the array IP as follows: If IP(j) = i then the jth column of P is the ith canonical unit vector.
The elements above the diagonal of the array MAT contain the corresponding elements of the triangular matrix R. The elements on the diagonal of R are stored in the array DIAGR.
R can then be partioned by defining R11 as the largest leading submatrix of R whose estimated condition number, in the 1-norm, is less than 1/TOL (or such that abs(R[j,j])>0 if TOL is absent). The order of R11, KRANK, is the effective rank of MAT.
This leads to the following partition of R:
[ R11 R12 ]
[ R21 R22 ]
where R21 is zero by construction (since R is an upper triangular or trapezoidal) and R22 is considered to be negligible.
If MIN_NORM is present and has the value true, R12 is annihilated by orthogonal transformations from the right, arriving at the complete orthogonal factorization:
MAT * P = Q * T * Z
, where Z is a n-by-n orthogonal matrix and T is a m-by-n matrix has the form:
[ T11 T12 ]
[ T21 T22 ]
, here T21 (=R21), T12 and T22 (=R22) are zero and T11 is a KRANK-by-KRANK upper triangular matrix.
The factorization is obtained by Householder’s method. The kth transformation matrix, Z(k), which is used to introduce zeros into the kth row of R, is given in the form
[ I 0 ]
[ 0 T(k) ]
where
T(k) = I + tau * ( u(k) * u(k)’ ) and u(k)’ = ( 1 0 z(k) )
tau is a scalar and z(k) is an (n-KRANK) element vector. tau and z(k) are chosen to annihilate the elements of the kth row of R12.
The Z n-by-n orthogonal matrix is given by
Z = Z(1) * Z(2) * … * Z(KRANK)
On exit, the scalar tau is returned in the kth element of TAU and the vector u(k) in the kth row of MAT, such that the elements of z(k) are in MAT(k,KRANK+1:n).
On exit, the elements above the diagonal of the array section MAT(1:KRANK,1:KRANK) contain the corresponding elements of the triangular matrix T11. The elements of the diagonal of T11 are stored in the array section DIAGR(1:KRANK). The last part of DIAGR is set to zero. In other words, T11 overwrites R11 and Z overwrites R12 on exit.
The minimum 2-norm solution is then
X = [ P * Z’ ](:,:KRANK) * [ inv(T11) * Q1’ * B ]
where inv(T11) is the inverse of T11, Z’ is the transpose of Z and Q1 consists of the first KRANK columns of Q.
If MIN_NORM is absent or has the value false, a solution is computed as
X = P(:,:KRANK) * [ inv(R11) * Q1’ * B ]
where inv(R11) is the inverse of R11 and Q1 consists of the first KRANK columns of Q. In this case, if the j-th column of MAT is omitted from the basis, X[j] is set to zero and R is not destroyed in MAT.
In both cases:
- The 2-norm of the residual vector for the solution X can be obtained through the optional argument RNORM .
- If COMP_RESID=true, The m residual vector B - MAT * X overwrites B on exit.
On input, if KRANK is present and KRANK=k, the first k columns of MAT are to be forced into the basis. Pivoting is performed on the last n-k columns of MAT.
When KRANK is present and KRANK>=min(m,n) is used, pivoting is not performed. This is appropriate when MAT is known to be full rank.
If KRANK is absent or is present with KRANK<=0, pivoting is done on all columns of MAT.
On output, if KRANK is present, it contains the rank of MAT, i.e., the order of the submatrix R11. This is the same as the order of the submatrix T11 in the complete orthogonal factorization of MAT.
TOL is an optional argument such that 0<=TOL<1. If TOL is not specified, or is outside [0,1[, the calculations to determine the condition number of MAT are not performed and crude tests on R(j,j) are performed to determine the numerical rank of MAT. If TOL is present and is in [0,1[, the calculations to determine the condition number are performed, the effective rank of MAT is determined and the reciprocal of the condition number is returned in TOL.
If it is possible that MAT may not be full rank (i.e., certain columns of MAT are linear combinations of other columns), then the linearly dependent columns can usually be eliminated by using KRANK=0 and TOL=relative precision of the elements in MAT and B. If each element is correct to, say, 5 digits then TOL=0.00001 should be used. Also, it may be helpful to scale the columns of MAT so that all elements are about the same order of magnitude.
subroutine llsq_qr_solve2 ( mat, b, x, rnorm, comp_resid, krank, tol, min_norm, diagr, beta, ip, tau )¶
Purpose¶
LLSQ_QR_SOLVE2 computes solutions to real linear least squares problems of the form:
Minimize 2-norm|| B - MAT * X ||using an orthogonal factorization with columns pivoting of MAT. MAT is an m-by-n matrix which may be rank-deficient. m>=n or n>m is permitted.
Several right hand side vectors b and solution vectors x can be handled in a single call; they are stored as the columns of the m-by-nb right hand side matrix B and the n-by-nb solution matrix X, respectively.
MAT and B are overwritten with information generated by LLSQ_QR_SOLVE2. The (complete) orthogonal factorization of MAT is saved in arguments MAT, DIAGR, BETA, IP and TAU on output.
Arguments¶
- MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)
On entry, the m-by-n coefficient matrix MAT.
On exit, MAT has been overwritten by details of its (complete) orthogonal factorization. Other parts of the factorization can be obtained if the optional arguments DIAGR, BETA, IP and TAU are present. See Further Details.
- B (INPUT/OUTPUT) real(stnd), dimension(:,:)
On entry, the m-by-nb right hand side matrix B.
On exit, if COMP_RESID is present and is equal true, the residual matrix B - MAT * X overwrites B on output.
The shape of B must verify:
- size( B, 1 ) = size( MAT, 1 ) = m
- size( B, 2 ) = size( X, 2 ) = nb .
- X (OUTPUT) real(stnd), dimension(:,:)
On exit, the n-by-nb solution matrix X.
The shape of X must verify:
- size( X, 1 ) = size( MAT, 2 ) = n
- size( X, 2 ) = size( B, 2 ) = nb .
- RNORM (OUTPUT, OPTIONAL) real(stnd), dimension(:)
On exit, the 2-norm of the residual vectors for the solutions stored columnwise in the matrix X.
- The size of RNORM must verify:
- size( RNORM ) = size( X, 2 ) = size( B, 2 ) = nb .
- COMP_RESID (INPUT, OPTIONAL) logical(lgl)
- On entry, if COMP_RESID is present and is equal true, the residual matrix B - MAT * X overwrites B on exit.
- KRANK (INPUT/OUTPUT, OPTIONAL) integer(i4b)
On entry, KRANK=k implies that the first k columns of MAT are to be forced into the basis, pivoting is performed on the last n-k columns of MAT. When KRANK >=min(m,n) is used, pivoting is not performed. This is appropriate when MAT is known to be full rank. If KRANK is absent or is <=0, pivoting is done on all columns of MAT. This is appropriate if MAT may not be of full rank (i.e. certain columns of MAT are linear combinations of other columns).
On exit, KRANK contains the effective rank of MAT, i.e., the number of independent columns in matrix MAT.
- TOL (INPUT/OUTPUT, OPTIONAL) real(stnd)
On entry, TOL is used to determine the effective rank of MAT, which is then defined as the order of the largest leading triangular submatrix R11 in the QR factorization (with pivoting) of MAT whose estimated condition number, in the 1-norm, is less than 1/TOL. TOL must be set to the relative precision of the elements in MAT and B. If each element is correct to, say, 5 digits then TOL=0.00001 should be used.
TOL must not be greater or equal to 1 or less than 0, otherwise the numerical rank of MAT is determined and the calculations to determine the condition number are not performed. If TOL=0, the numerical rank of MAT is determined, but the condition number is calculated.
On exit, if a condition number is calculated, its reciprocal is returned in TOL. Otherwise, TOL is not changed.
If TOL is absent, the numerical rank of MAT is used and is returned in the optional argument KRANK.
- MIN_NORM (INPUT, OPTIONAL) logical(lgl)
- On entry: If MIN_NORM=true, a complete orthogonal factorization of MAT and the minimun 2-norm solution are computed. If MIN_NORM=false or if MIN_NORM is absent, a QR factorization with column pivoting of MAT and a solution are computed such that if the j-th column of MAT is omitted from the basis, X[j] is set to zero.
- DIAGR (OUTPUT, OPTIONAL) real(stnd), dimension(:)
On exit, the diagonal elements of the matrix R if MIN_NORM=false or is absent, or the diagonal elements of the matrix T11 if MIN_NORM is present and is equal to true. The diagonal elements of T11 are stored in DIAGR(1:KRANK). See Further Details.
The size of DIAGR must be min( size(MAT,1) , size(MAT,2) ).
- BETA (OUTPUT, OPTIONAL) real(stnd), dimension(:)
On exit, the scalars factors of the elementary reflectors defining Q. See Further Details.
The size of BETA must be min( size(MAT,1) , size(MAT,2) ).
- IP (OUTPUT, OPTIONAL) integer(i4b), dimension(:)
On exit, if IP(j)=k, then the j-th column of MAT * P was the k-th column of MAT. See Further Details.
The size of IP must be size(MAT,2) = n.
- TAU (OUTPUT, OPTIONAL) real(stnd), dimension(:)
On exit, the scalars factors of the elementary reflectors defining Z in the complete orthogonal factorization of MAT if MIN_NORM is present and is equal to true. Otherwise, TAU is set to 0. See Further Details.
The size of TAU must be min( size(MAT,1) , size(MAT,2) ).
Further Details¶
The routine first computes a QR factorization with (partial) column pivoting on option (see below):
MAT * P = Q * R
, here P is n-by-n permutation matrix, R is an upper triangular or trapezoidal (if n>m) matrix and Q is a m-by-m orthogonal matrix.
The matrix Q is represented as a product of elementary reflectors
Q = H(1) * H(2) * … * H(k), where k = min( size(MAT,1) , size(MAT,2) )
Each H(i) has the form
H(i) = I + beta * ( v * v’ ) ,
where beta is a real scalar and v is a real m-element vector with v(1:i-1) = 0. v(i:m) is stored on exit in MAT(i:m,i) and beta in BETA(i).
The matrix P is represented in the array IP as follows: If IP(j) = i then the jth column of P is the ith canonical unit vector.
The elements above the diagonal of the array MAT contain the corresponding elements of the triangular matrix R. The elements on the diagonal of R are stored in the array DIAGR.
R can then be partioned by defining R11 as the largest leading submatrix of R whose estimated condition number, in the 1-norm, is less than 1/TOL (or such that abs(R[j,j])>0 if TOL is absent). The order of R11, KRANK, is the effective rank of MAT.
This leads to the following partition of R:
[ R11 R12 ]
[ R21 R22 ]
where R21 is zero by construction (since R is an upper triangular or trapezoidal) and R22 is considered to be negligible.
If MIN_NORM is present and has the value true, R12 is annihilated by orthogonal transformations from the right, arriving at the complete orthogonal factorization:
MAT * P = Q * T * Z
, where Z is a n-by-n orthogonal matrix and T has the form:
[ T11 T12 ]
[ T21 T22 ]
, here T21 (=R21), T12 and T22 (=R22) are zero and T11 is a KRANK-by-KRANK upper triangular matrix.
The factorization is obtained by Householder’s method. The kth transformation matrix, Z(k), which is used to introduce zeros into the kth row of R, is given in the form
[ I 0 ]
[ 0 T(k) ]
where
T(k) = I + tau * ( u(k) * u(k)’ ) and u(k)’ = ( 1 0 z(k) )
tau is a scalar and z(k) is an (n-KRANK) element vector. tau and z(k) are chosen to annihilate the elements of the kth row of R12.
The Z n-by-n orthogonal matrix is given by
Z = Z(1) * Z(2) * … * Z(KRANK)
On exit, the scalar tau is returned in the kth element of TAU and the vector u(k) in the kth row of MAT, such that the elements of z(k) are in MAT(k,KRANK+1:n).
On exit, the elements above the diagonal of the array section MAT(1:KRANK,1:KRANK) contain the corresponding elements of the triangular matrix T11. The elements of the diagonal of T11 are stored in the array section DIAGR(1:KRANK). The last part of DIAGR is set to zero. In other words, T11 overwrites R11 and Z overwrites R12 on exit.
The minimum 2-norm solution is then
X = [ P * Z’ ](:,:KRANK) * [ inv(T11) * Q1’ * B ]
where inv(T11) is the inverse of T11, Z’ is the transpose of Z and Q1 consists of the first KRANK columns of Q.
If MIN_NORM is absent or has the value false, a solution is computed as
X = P(:,:KRANK) * [ inv(R11) * Q1’ * B ]
where inv(R11) is the inverse of R11 and Q1 consists of the first KRANK columns of Q. In this case, if the j-th column of MAT is omitted from the basis, X[j] is set to zero and R is not destroyed in MAT.
In both cases:
- The 2-norm of the residual vector for the solution in the j-th column of X is given in RNORM[j] if argument RNORM is present.
- If COMP_RESID=true, The residual matrix B - MAT * X overwrites B on exit.
On input, if KRANK is present and KRANK=k, the first k columns of MAT are to be forced into the basis. Pivoting is performed on the last n-k columns of MAT.
When KRANK is present and KRANK>=min(m,n) is used, pivoting is not performed. This is appropriate when MAT is known to be full rank.
If KRANK is absent or is present with KRANK<=0, pivoting is done on all columns of MAT.
On output, if KRANK is present, it contains the rank of MAT, i.e., the order of the submatrix R11. This is the same as the order of the submatrix T11 in the complete orthogonal factorization of MAT.
TOL is an optional argument such that 0<=TOL<1. If TOL is not specified, or is outside [0,1[, the calculations to determine the condition number of MAT are not performed and crude tests on R(j,j) are performed to determine the numerical rank of MAT. If TOL is present and is in [0,1[, the calculations to determine the condition number are performed, the effective rank of MAT is determined and the reciprocal of the condition number is returned in TOL.
If it is possible that MAT may not be full rank (i.e., certain columns of MAT are linear combinations of other columns), then the linearly dependent columns can usually be eliminated by using KRANK=0 and TOL=relative precision of the elements in MAT and B. If each element is correct to, say, 5 digits then TOL=0.00001 should be used. Also, it may be helpful to scale the columns of MAT so that all elements are about the same order of magnitude.
subroutine llsq_qr_solve2 ( vec, b, x, rnorm, comp_resid, diagr, beta )¶
Purpose¶
LLSQ_QR_SOLVE2 computes a solution to a real linear least squares problem:
Minimize 2-norm|| B - VEC * X ||VEC is an m vector, B is a m right hand side vector and X is a real scalar.
VEC and B are overwritten with information generated by LLSQ_QR_SOLVE2.
Arguments¶
- VEC (INPUT/OUTPUT) real(stnd), dimension(:)
On entry, the m coefficient vector VEC.
On exit, VEC contains the vector v of the Householder reflector H. See Further Details.
- B (INPUT/OUTPUT) real(stnd), dimension(:)
On entry, the m right hand side vector B.
On exit, if COMP_RESID is present and is equal true, the residual vector B - VEC * X overwrites B on output.
- The shape of B must verify:
- size( B ) = size( VEC ) = m .
- X (OUTPUT) real(stnd)
- On exit, the real solution X.
- RNORM (OUTPUT, OPTIONAL) real(stnd)
- On exit, the 2-norm of the residual vector for the solution scalar X.
- COMP_RESID (INPUT, OPTIONAL) logical(lgl)
- On entry, if COMP_RESID is present and is equal true, the residual vector B - VEC * X overwrites B on exit.
- DIAGR (OUTPUT, OPTIONAL) real(stnd)
- On exit, the scalar DIAGR. See Further Details.
- BETA (OUTPUT, OPTIONAL) real(stnd)
- On exit, the scalar factor BETA of the elementary reflector defining H. See Further Details.
Further Details¶
The routine first generates a real elementary reflector H of order m, such that
H * VEC = D , with H’ * H = I and D’ = ( DIAGR 0 )
where DIAGR is scalar. H is represented in the form
H = I + BETA * ( v * v’ ) ,
where BETA is a real scalar and v is a real m-element vector.
The solution X is then computed as
X = [ H * B ](1) / DIAGR
The 2-norm of the residual vector for the solution X can be obtained through the optional argument RNORM as
2-norm|| [ H * B ](2:) ||
If COMP_RESID=true, The residual vector B - VEC * X overwrites B on exit.
subroutine llsq_qr_solve2 ( vec, b, x, rnorm, comp_resid, diagr, beta )¶
Purpose¶
LLSQ_QR_SOLVE2 computes solutions to real linear least squares problems of the form:
Minimize 2-norm|| B - VEC * X ||here VEC is an m vector and several right hand side vectors b and solution scalars x can be handled in a single call; they are stored as the columns of the m-by-nb right hand side matrix B and the nb solution vector X, respectively.
VEC and B are overwritten with information generated by LLSQ_QR_SOLVE2.
Arguments¶
- VEC (INPUT/OUTPUT) real(stnd), dimension(:)
On entry, the m coefficient vector VEC.
On exit, VEC contains the vector v of the Householder reflector H. See Further Details.
- B (INPUT/OUTPUT) real(stnd), dimension(:,:)
On entry, the m-by-nb right hand side matrix B.
On exit, if COMP_RESID is present and is equal true, the residual matrix B - VEC * X overwrites B on output.
The shape of B must verify:
- size( B, 1 ) = size( VEC ) = m
- size( B, 2 ) = size( X ) = nb .
- X (OUTPUT) real(stnd), dimension(:)
On exit, the nb solution vector X.
- The shape of X must verify:
- size( X ) = size( B, 2 ) = nb .
- RNORM (OUTPUT, OPTIONAL) real(stnd), dimension(:)
On exit, the 2-norms of the residual vectors for the solutions stored in the vector X.
- The size of RNORM must verify:
- size( RNORM ) = size( X ) = size( B, 2 ) = nb .
- COMP_RESID (INPUT, OPTIONAL) logical(lgl)
- On entry, if COMP_RESID is present and is equal true, the residual matrix B - VEC * X overwrites B on exit.
- DIAGR (OUTPUT, OPTIONAL) real(stnd)
- On exit, the scalar DIAGR. See Further Details.
- BETA (OUTPUT, OPTIONAL) real(stnd)
- On exit, the scalar factor BETA of the elementary reflector defining H. See Further Details.
Further Details¶
The routine first generates a real elementary reflector H of order m, such that
H * VEC = D , with H’ * H = I and D’ = ( DIAGR 0 )
where DIAGR is scalar. H is represented in the form
H = I + BETA * ( v * v’ ) ,
where BETA is a real scalar and v is a real m-element vector.
The solution vector X is then computed as
X(:) = [ H * B ](1,:) / DIAGR
The 2-norm of the residual vector for the solution X[j] is given in RNORM[j] if argument RNORM is present.
If COMP_RESID=true, The residual matrix B - VEC * X overwrites B on exit.
subroutine qr_solve ( mat, diagr, beta, b, x, rnorm, comp_resid )¶
Purpose¶
QR_SOLVE solves overdetermined or underdetermined real linear systems
MAT * X = Bwith an m-by-n matrix MAT, using a QR factorization of MAT as computed by QR_CMP. m>=n or n>m is permitted, but it is assumed that MAT has full rank. B is a m right hand side vector and X is a n solution vector.
It is assumed that QR_CMP has been used to compute the QR factorization of MAT before QR_SOLVE.
Arguments¶
- MAT (INPUT) real(stnd), dimension(:,:)
- On entry, the QR factorization of the real coefficient matrix MAT as returned by QR_CMP. The elements above the diagonal of the array contain the corresponding elements of R. The elements on and below the diagonal, with the array BETA, represent the orthogonal matrix Q in the QR decomposition of MAT, as a product of elementary reflectors, as returned by QR_CMP.
- DIAGR (INPUT) real(stnd), dimension(:)
On entry, the diagonal elements of the matrix R in the QR decomposition of MAT.
The size of DIAGR must be min( size(MAT,1) , size(MAT,2) ).
- BETA (INPUT) real(stnd), dimension(:)
On entry, the scalars factors of the elementary reflectors defining Q, as returned by QR_CMP.
The size of BETA must be min( size(MAT,1) , size(MAT,2) ),
- B (INPUT/OUTPUT) real(stnd), dimension(:)
On entry, the m right hand side vector B.
On exit, if COMP_RESID is present and is equal true, the residual vector B - MAT * X overwrites B.
- The shape of B must verify:
- size( B ) = size( MAT, 1 ) = m .
- X (OUTPUT) real(stnd), dimension(:)
On exit, the n solution vector X.
- The shape of X must verify:
- size( X ) = size( MAT, 2 ) = n .
- RNORM (OUTPUT, OPTIONAL) real(stnd)
- On exit, the 2-norm of the residual vector for the solution vector X.
- COMP_RESID (INPUT, OPTIONAL) logical(lgl)
- On entry, if COMP_RESID is present and is equal true, the residual vector B - MAT * X overwrites B on exit.
Further Details¶
It is assumed that QR_CMP has been used to compute the QR factorization of MAT before calling QR_SOLVE.
If m>=n: the subroutine find the least squares solution of an overdetermined system, i.e., solve the least squares problem
Minimize 2-norm|| B - MAT * X ||
If m<n: the subroutine find a solution of an underdetermined system
MAT * X = B
The 2-norm of the residual vector for the solution X can be obtained through the optional argument RNORM .
If COMP_RESID=true, The m residual vector B - MAT * X overwrites B on exit.
MAT, DIAGR, BETA are not modified by this routine and can be left in place for successive calls with different right-hand side vectors B.
subroutine qr_solve ( mat, diagr, beta, b, x, rnorm, comp_resid )¶
Purpose¶
QR_SOLVE solves overdetermined or underdetermined real linear systems of the form:
MAT * X = Bwith an m-by-n matrix MAT, using a QR factorization of MAT as computed by QR_CMP. m>=n or n>m is permitted, but it is assumed that MAT has full rank.
Several right hand side vectors b and solution vectors x can be handled in a single call; they are stored as the columns of the m-by-nb right hand side matrix B and the n-by-nb solution matrix X, respectively.
It is assumed that QR_CMP has been used to compute the QR factorization of MAT before QR_SOLVE.
Arguments¶
- MAT (INPUT) real(stnd), dimension(:,:)
- On entry, the QR factorization of the real coefficient matrix MAT as returned by QR_CMP. The elements above the diagonal of the array contain the corresponding elements of R. The elements on and below the diagonal, with the array BETA, represent the orthogonal matrix Q in the QR decomposition of MAT, as a product of elementary reflectors, as returned by QR_CMP.
- DIAGR (INPUT) real(stnd), dimension(:)
On entry, the diagonal elements of the matrix R in the QR decomposition of MAT.
The size of DIAGR must be min( size(MAT,1) , size(MAT,2) ).
- BETA (INPUT) real(stnd), dimension(:)
On entry, the scalars factors of the elementary reflectors defining Q, as returned by QR_CMP.
The size of BETA must be min( size(MAT,1) , size(MAT,2) ),
- B (INPUT/OUTPUT) real(stnd), dimension(:,:)
On entry, the m-by-nb right hand side matrix B.
On exit, if COMP_RESID is present and is equal true, the residual matrix B - MAT * X overwrites B on output.
The shape of B must verify:
- size( B, 1 ) = size( MAT, 1 ) = m
- size( B, 2 ) = size( X, 2 ) = nb .
- X (OUTPUT) real(stnd), dimension(:,:)
On exit, the n-by-nb solution matrix X.
The shape of X must verify:
- size( X, 1 ) = size( MAT, 2 ) = n
- size( X, 2 ) = size( B, 2 ) = nb .
- RNORM (OUTPUT, OPTIONAL) real(stnd), dimension(:)
On exit, the 2-norm of the residual vectors for the solutions stored columnwise in the matrix X.
- The size of RNORM must verify:
- size( RNORM ) = size( X, 2 ) = size( B, 2 ) = nb .
- COMP_RESID (INPUT, OPTIONAL) logical(lgl)
- On entry, if COMP_RESID is present and is equal true, the residual matrix B - MAT * X overwrites B on exit.
Further Details¶
It is assumed that QR_CMP has been used to compute the QR factorization of MAT before calling QR_SOLVE.
If m>=n: the subroutine find the least squares solutions of overdetermined systems, i.e., solve least squares problems of the form
Minimize 2-norm|| B - MAT * X ||
If m<n: the subroutine find solutions of underdetermined systems of the form
MAT * X = B
In both cases, several right hand side vectors b and solution vectors x can be handled in a single call; they are stored as the columns of the m-by-nb right hand side matrix B and the n-by-nb solution matrix X, respectively.
The 2-norm of the residual vector for the solution in the j-th column of X is given in RNORM[j] if argument RNORM is present.
If COMP_RESID=true, The residual matrix B - MAT * X overwrites B on exit.
MAT, DIAGR, BETA are not modified by this routine and can be left in place for successive calls with different right-hand side matrices B.
subroutine qr_solve2 ( mat, diagr, beta, ip, krank, b, x, rnorm, comp_resid, tau )¶
Purpose¶
QR_SOLVE2 solves overdetermined or underdetermined real linear systems
MAT * X = Bwith an m-by-n matrix MAT, using a (complete) orthogonal factorization of MAT as computed by QR_CMP2. m>=n or n>m is permitted and MAT may be rank-deficient. B is a m right hand side vector and X is a n solution vector.
It is assumed that QR_CMP2 has been used to compute the (complete) orthogonal factorization of MAT before QR_SOLVE2.
Arguments¶
- MAT (INPUT) real(stnd), dimension(:,:)
- On entry, details of the (complete) orthogonal factorization of the real coefficient matrix MAT as returned by QR_CMP2. See description of QR_CMP2 subroutine.
- DIAGR (INPUT) real(stnd), dimension(:)
On entry, the diagonal elements of the matrix R in the QR decomposition with column pivoting of MAT if TAU is absent or the diagonal elements of the matrix T11 in the complete orthogonal factorization of MAT if TAU is present, as computed by QR_CMP2. The diagonal elements of T11 are stored in DIAGR(1:KRANK). See description of QR_CMP2 subroutine.
The size of DIAGR must be min( size(MAT,1) , size(MAT,2) ).
- BETA (INPUT) real(stnd), dimension(:)
On entry, the scalars factors of the elementary reflectors defining Q, as returned by QR_CMP2. See description of QR_CMP2 subroutine.
The size of BETA must be min( size(MAT,1) , size(MAT,2) ),
- IP (INPUT) integer(i4b), dimension(:)
On entry, the permutation P in the (complete) orthogonal factorization of MAT, as returned by QR_CMP2. See description of QR_CMP2 subroutine.
The size of IP must be size(MAT,2) = n.
- KRANK (INPUT) integer(i4b)
- On entry, KRANK contains the effective rank of MAT, as returned by QR_CMP2. See description of QR_CMP2 subroutine.
- B (INPUT/OUTPUT) real(stnd), dimension(:)
On entry, the m right hand side vector B.
On exit, if COMP_RESID is present and is equal true, the residual vector B - MAT * X overwrites B.
- The shape of B must verify:
- size( B ) = size( MAT, 1 ) = m .
- X (OUTPUT) real(stnd), dimension(:)
On exit, the n solution vector X.
- The shape of X must verify:
- size( X ) = size( MAT, 2 ) = n .
- RNORM (OUTPUT, OPTIONAL) real(stnd)
- On exit, the 2-norm of the residual vector for the solution vector X.
- COMP_RESID (INPUT, OPTIONAL) logical(lgl)
- On entry, if COMP_RESID is present and is equal true, the residual vector B - MAT * X overwrites B on exit.
- TAU (INPUT, OPTIONAL) real(stnd), dimension(:)
On entry, if TAU is present, a complete orthogonal factorization of MAT has been computed by QR_CMP2 and TAU contains the scalars factors of the elementary reflectors defining Z in this decomposition. Otherwise, only a QR factorization with column pivoting of MAT has been computed by QR_CMP2. See description of QR_CMP2 subroutine.
The size of TAU must be min( size(MAT,1) , size(MAT,2) ).
Further Details¶
It is assumed that QR_CMP2 has been used to compute the (complete) orthogonal factorization (TAU present) or the QR factorization with column pivoting (TAU absent) of MAT before calling QR_SOLVE2.
If m>=n: the subroutine find the least squares solution of an overdetermined system, i.e., solve the least squares problem
Minimize 2-norm|| B - MAT * X ||
If m<n: the subroutine find a solution of an underdetermined system
MAT * X = B
In both cases, the minimun 2-norm solution is computed if TAU is present. Otherwise, a solution is computed such that if the j-th column of MAT is omitted from the basis, X[j] is set to zero.
The 2-norm of the residual vector for the solution X can be obtained through the optional argument RNORM .
If COMP_RESID=true, The m residual vector B - MAT * X overwrites B on exit.
MAT, DIAGR, BETA, IP, KRANK and TAU are not modified by this routine and can be left in place for successive calls with different right-hand side vectors B.
subroutine qr_solve2 ( mat, diagr, beta, ip, krank, b, x, rnorm, comp_resid, tau )¶
Purpose¶
QR_SOLVE2 solves overdetermined or underdetermined real linear systems of the form:
MAT * X = Bwith an m-by-n matrix MAT, using a (complete) orthogonal factorization of MAT as computed by QR_CMP2. m>=n or n>m is permitted and MAT may be rank-deficient.
Several right hand side vectors b and solution vectors x can be handled in a single call; they are stored as the columns of the m-by-nb right hand side matrix B and the n-by-nb solution matrix X, respectively.
It is assumed that QR_CMP2 has been used to compute the (complete) orthogonal factorization of MAT before QR_SOLVE2.
Arguments¶
- MAT (INPUT) real(stnd), dimension(:,:)
- On entry, details of the (complete) orthogonal factorization of the real coefficient matrix MAT as returned by QR_CMP2. See description of QR_CMP2 subroutine.
- DIAGR (INPUT) real(stnd), dimension(:)
On entry, the diagonal elements of the matrix R in the QR decomposition with column pivoting of MAT if TAU is absent or the diagonal elements of the matrix T11 in the complete orthogonal factorization of MAT if TAU is present, as computed by QR_CMP2. The diagonal elements of T11 are stored in DIAGR(1:KRANK). See description of QR_CMP2 subroutine.
The size of DIAGR must be min( size(MAT,1) , size(MAT,2) ).
- BETA (INPUT) real(stnd), dimension(:)
On entry, the scalars factors of the elementary reflectors defining Q, as returned by QR_CMP2. See description of QR_CMP2 subroutine.
The size of BETA must be min( size(MAT,1) , size(MAT,2) ),
- IP (INPUT) integer(i4b), dimension(:)
On entry, the permutation P in the (complete) orthogonal factorization of MAT, as returned by QR_CMP2. See description of QR_CMP2 subroutine.
The size of IP must be size(MAT,2) = n.
- KRANK (INPUT) integer(i4b)
- On entry, KRANK contains the effective rank of MAT, as returned by QR_CMP2. See description of QR_CMP2 subroutine.
- B (INPUT/OUTPUT) real(stnd), dimension(:,:)
On entry, the m-by-nb right hand side matrix B.
On exit, if COMP_RESID is present and is equal true, the residual matrix B - MAT * X overwrites B on output.
The shape of B must verify:
- size( B, 1 ) = size( MAT, 1 ) = m
- size( B, 2 ) = size( X, 2 ) = nb .
- X (OUTPUT) real(stnd), dimension(:,:)
On exit, the n-by-nb solution matrix X.
The shape of X must verify:
- size( X, 1 ) = size( MAT, 2 ) = n
- size( X, 2 ) = size( B, 2 ) = nb .
- RNORM (OUTPUT, OPTIONAL) real(stnd), dimension(:)
On exit, the 2-norm of the residual vectors for the solutions stored columnwise in the matrix X.
- The size of RNORM must verify:
- size( RNORM ) = size( X, 2 ) = size( B, 2 ) = nb .
- COMP_RESID (INPUT, OPTIONAL) logical(lgl)
- On entry, if COMP_RESID is present and is equal true, the residual matrix B - MAT * X overwrites B on exit.
- TAU (INPUT, OPTIONAL) real(stnd), dimension(:)
On entry, if TAU is present, a complete orthogonal factorization of MAT has been computed by QR_CMP2 and TAU contains the scalars factors of the elementary reflectors defining Z in this decomposition. Otherwise, only a QR factorization with column pivoting of MAT has been computed by QR_CMP2. See description of QR_CMP2 subroutine.
The size of TAU must be min( size(MAT,1) , size(MAT,2) ).
Further Details¶
It is assumed that QR_CMP2 has been used to compute the (complete) orthogonal factorization (TAU present) or the QR factorization with column pivoting (TAU absent) of MAT before calling QR_SOLVE2.
If m>=n: the subroutine find the least squares solutions of overdetermined systems, i.e., solve least squares problems of the form
Minimize 2-norm|| B - MAT * X ||
If m<n: the subroutine find solutions of underdetermined systems of the form
MAT * X = B
In both cases, several right hand side vectors b and solution vectors x can be handled in a single call; they are stored as the columns of the m-by-nb right hand side matrix B and the n-by-nb solution matrix X, respectively.
In both cases, the minimun 2-norm solutions are computed if TAU is present. Otherwise, solutions are computed such that if the j-th column of MAT is omitted from the basis, X[j,:] is set to zero.
The 2-norm of the residual vector for the solution in the j-th column of X is given in RNORM[j] if argument RNORM is present.
If COMP_RESID=true, The residual matrix B - MAT * X overwrites B on exit.
MAT, DIAGR, BETA, IP, KRANK and TAU are not modified by this routine and can be left in place for successive calls with different right-hand side matrices B.
subroutine llsq_svd_solve ( mat, b, failure, x, singvalues, krank, rnorm, tol, mul_size, maxiter, max_francis_steps, perfect_shift )¶
Purpose¶
LLSQ_SVD_SOLVE computes the minimum norm solution to a real linear least squares problem:
Minimize 2-norm|| B - MAT * X ||using the singular value decomposition (SVD) of MAT. MAT is an m-by-n matrix which may be rank-deficient.
Several right hand side vectors b and solution vectors x can be handled in a single call; they are stored as the columns of the m-by-nrhs right hand side matrix B and the n-by-nrhs solution matrix X, respectively.
The effective rank of MAT, KRANK, is determined by treating as zero those singular values which are less than TOL times the largest singular value.
Arguments¶
- MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)
On entry, the m-by-n matrix MAT.
On exit, MAT is destroyed. If m>=n, MAT(:n,:n) is overwritten with the right singular vectors of MAT, stored columnwise.
- B (INPUT/OUTPUT) real(stnd), dimension(:,:)
On entry, the m-by-nrhs right hand side matrix B.
On exit, B is destroyed. If m>KRANK, the residual sum-of-squares for the solution in the i-th column is given by the sum of squares of elements KRANK+1:m in that column.
The shape of B must verify:
- size( B, 1 ) = size( MAT, 1 ) = m
- size( B, 2 ) = size( X, 2 ) = nrhs .
- FAILURE (OUTPUT) logical(lgl)
If:
- FAILURE= false : indicates successful exit
- FAILURE= true : indicates that the SVD algorithm did not converge and that full accuracy was not attained in the bidiagonal SVD of an intermediate bidiagonal form BD of MAT.
- X (OUTPUT) real(stnd), dimension(:,:)
On exit, the n-by-nrhs solution matrix X.
The shape of X must verify:
- size( X, 1 ) = size( MAT, 2 ) = n
- size( X, 2 ) = size( B, 2 ) = nrhs .
- SINGVALUES (OUTPUT, OPTIONAL) real(stnd), dimension(:)
The singular values of MAT in decreasing order. The condition number of MAT in the 2-norm is
SINGVALUES(1)/SINGVALUES(min(m,n)).The size of SINGVALUES must verify: size( SINGVALUES ) = min(m,n) .
- KRANK (OUTPUT, OPTIONAL) integer(i4b)
- On exit, the effective rank of MAT, i.e., the number of singular values which are greater than TOL * SINGVALUES(1).
- RNORM (OUTPUT, OPTIONAL) real(stnd), dimension(:)
On exit, the 2-norms of the residual vectors for the solutions stored columnwise in the matrix X.
- The size of RNORM must verify:
- size( RNORM ) = size( X, 2 ) = size( B, 2 ) = nrhs .
- TOL (INPUT, OPTIONAL) real(stnd)
- On entry, TOL is used to determine the effective rank of MAT. Singular values SINGVALUES(i) <= TOL * SINGVALUES(1) are treated as zero. If TOL is absent, machine precision is used instead.
- MUL_SIZE (INPUT, OPTIONAL) integer(i4b)
Internal parameter. MUL_SIZE must verify 1 <= MUL_SIZE <= max(m,n), otherwise a default value is used. MUL_SIZE can be increased or decreased to improve the performance of the algorithm used in LLSQ_SVD_SOLVE. Maximum performance will be obtained when a real matrix of size MUL_SIZE * max(m,n) and kind stnd fits in the cache of the processors.
The default value is 32.
- MAXITER (INPUT, OPTIONAL) integer(i4b)
MAXITER controls the maximum number of QR sweeps in the bidiagonal SVD phase of the SVD algorithm. The bidiagonal SVD algorithm of an intermediate bidiagonal form B of MAT fails to converge if the number of QR sweeps exceeds MAXITER * min(m,n). Convergence usually occurs in about 2 * min(m,n) QR sweeps.
The default is 10.
- MAX_FRANCIS_STEPS (INPUT,OPTIONAL) integer(i4b)
MAX_FRANCIS_STEPS controls the maximum number of Francis sets (eg QR sweeps) of Givens rotations which must be saved before applying them with a wavefront algorithm to accumulate the singular vectors in the bidiagonal SVD algorithm. MAX_FRANCIS_STEPS is a strictly positive integer, otherwise the default value is used.
The default is 10.
- PERFECT_SHIFT (INPUT,OPTIONAL) logical(lgl)
PERFECT_SHIFT determines if a perfect shift strategy is used in the implicit QR algorithm in order to minimize the number of QR sweeps in the bidiagonal SVD algorithm.
The default is true.
Further Details¶
This subroutine is adapted from the routine DGELSS in LAPACK. If OPENMP is used, the algorithm is parallelized.
For further details, on using the SVD for solving a least square problem, see the references (1) or (2).
- Golub, G.H., and Van Loan, C.F., 1996: Matrix Computations. The Johns
- Hopkins University Press, Baltimore, Maryland.
- Lawson, C.L., and Hanson, R.J., 1974: Solving least square problems.
- Prentice-Hall.
subroutine llsq_svd_solve ( mat, b, failure, x, singvalues, krank, rnorm, tol, mul_size, maxiter, max_francis_steps, perfect_shift )¶
Purpose¶
LLSQ_SVD_SOLVE computes the minimum norm solution to a real linear least squares problem:
Minimize 2-norm|| B - MAT * X ||using the singular value decomposition (SVD) of MAT. MAT is an m-by-n matrix which may be rank-deficient, B is a m right hand side vector and X is a n solution vector.
The effective rank of MAT, KRANK, is determined by treating as zero those singular values which are less than TOL times the largest singular value.
Arguments¶
- MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)
On entry, the m-by-n matrix MAT.
On exit, MAT is destroyed. If m>=n, MAT(:n,:n) is overwritten with the right singular vectors of MAT, stored columnwise.
- B (INPUT/OUTPUT) real(stnd), dimension(:)
On entry, the m right hand side vector B.
On exit, B is destroyed. If m>KRANK, the residual sum-of-squares for the solution X is given by the sum of squares of elements KRANK+1:m of B .
- The size of B must verify:
- size( B ) = size( MAT, 1 ) = m .
- FAILURE (OUTPUT) logical(lgl)
If:
- FAILURE = false : indicates successful exit
- FAILURE = true : indicates that the SVD algorithm did not converge and that full accuracy was not attained in the bidiagonal SVD of an intermediate bidiagonal form BD of MAT.
- X (OUTPUT) real(stnd), dimension(:)
On exit, the n solution vector X.
- The size of X must verify:
- size( X ) = size( MAT, 2 ) = n .
- SINGVALUES (OUTPUT, OPTIONAL) real(stnd), dimension(:)
The singular values of MAT in decreasing order. The condition number of MAT in the 2-norm is
SINGVALUES(1)/SINGVALUES(min(m,n)).The size of SINGVALUES must verify : size( SINGVALUES ) = min(m,n) .
- KRANK (OUTPUT, OPTIONAL) integer(i4b)
- On exit, the effective rank of MAT, i.e., the number of singular values which are greater than TOL * SINGVALUES(1).
- RNORM (OUTPUT, OPTIONAL) real(stnd)
- On exit, the 2-norm of the residual vector for the solution vector X.
- TOL (INPUT, OPTIONAL) real(stnd)
- On entry, TOL is used to determine the effective rank of MAT. Singular values SINGVALUES(i) <= TOL * SINGVALUES(1) are treated as zero. If TOL is absent, machine precision is used instead.
- MUL_SIZE (INPUT, OPTIONAL) integer(i4b)
Internal parameter. MUL_SIZE must verify 1 <= MUL_SIZE <= max(m,n), otherwise a default value is used. MUL_SIZE can be increased or decreased to improve the performance of the algorithm used in LLSQ_SVD_SOLVE. Maximum performance will be obtained when a real matrix of size MUL_SIZE * max(m,n) and kind stnd fits in the cache of the processors.
The default value is 32.
- MAXITER (INPUT, OPTIONAL) integer(i4b)
MAXITER controls the maximum number of QR sweeps in the bidiagonal SVD phase of the SVD algorithm. The bidiagonal SVD algorithm of an intermediate bidiagonal form B of MAT fails to converge if the number of QR sweeps exceeds MAXITER * min(m,n). Convergence usually occurs in about 2 * min(m,n) QR sweeps.
The default is 10.
- MAX_FRANCIS_STEPS (INPUT,OPTIONAL) integer(i4b)
MAX_FRANCIS_STEPS controls the maximum number of Francis sets (eg QR sweeps) of Givens rotations which must be saved before applying them with a wavefront algorithm to accumulate the singular vectors in the bidiagonal SVD algorithm. MAX_FRANCIS_STEPS is a strictly positive integer, otherwise the default value is used.
The default is 10.
- PERFECT_SHIFT (INPUT,OPTIONAL) logical(lgl)
PERFECT_SHIFT determines if a perfect shift strategy is used in the implicit QR algorithm in order to minimize the number of QR sweeps in the bidiagonal SVD algorithm.
The default is true.
Further Details¶
This subroutine is adapted from the routine DGELSS in LAPACK. If OPENMP is used, the algorithm is parallelized.
For further details on using the SVD for solving a least square problem, see the references (1) or (2).
- Golub, G.H., and Van Loan, C.F., 1996: Matrix Computations. The Johns
- Hopkins University Press, Baltimore, Maryland.
- Lawson, C.L., and Hanson, R.J., 1974: Solving least square problems.
- Prentice-Hall.