NAG Library Routine Document
f08kcf (dgelsd)
1
Purpose
f08kcf (dgelsd) computes the minimum norm solution to a real linear least squares problem
2
Specification
Fortran Interface
Subroutine f08kcf ( |
m, n, nrhs, a, lda, b, ldb, s, rcond, rank, work, lwork, iwork, info) |
Integer, Intent (In) | :: | m, n, nrhs, lda, ldb, lwork | Integer, Intent (Inout) | :: | iwork(*) | Integer, Intent (Out) | :: | rank, info | Real (Kind=nag_wp), Intent (In) | :: | rcond | Real (Kind=nag_wp), Intent (Inout) | :: | a(lda,*), b(ldb,*), s(*) | Real (Kind=nag_wp), Intent (Out) | :: | work(max(1,lwork)) |
|
C Header Interface
#include <nagmk26.h>
void |
f08kcf_ (const Integer *m, const Integer *n, const Integer *nrhs, double a[], const Integer *lda, double b[], const Integer *ldb, double s[], const double *rcond, Integer *rank, double work[], const Integer *lwork, Integer iwork[], Integer *info) |
|
The routine may be called by its
LAPACK
name dgelsd.
3
Description
f08kcf (dgelsd) uses the singular value decomposition (SVD) of , where is a real by matrix which may be rank-deficient.
Several right-hand side vectors and solution vectors can be handled in a single call; they are stored as the columns of the by right-hand side matrix and the by solution matrix .
The problem is solved in three steps:
1. |
reduce the coefficient matrix to bidiagonal form with Householder transformations, reducing the original problem into a ‘bidiagonal least squares problem’ (BLS); |
2. |
solve the BLS using a divide-and-conquer approach; |
3. |
apply back all the Householder transformations to solve the original least squares problem. |
The effective rank of
is determined by treating as zero those singular values which are less than
rcond times the largest singular value.
4
References
Anderson E, Bai Z, Bischof C, Blackford S, Demmel J, Dongarra J J, Du Croz J J, Greenbaum A, Hammarling S, McKenney A and Sorensen D (1999)
LAPACK Users' Guide (3rd Edition) SIAM, Philadelphia
http://www.netlib.org/lapack/lug
Golub G H and Van Loan C F (1996) Matrix Computations (3rd Edition) Johns Hopkins University Press, Baltimore
5
Arguments
- 1: – IntegerInput
-
On entry: , the number of rows of the matrix .
Constraint:
.
- 2: – IntegerInput
-
On entry: , the number of columns of the matrix .
Constraint:
.
- 3: – IntegerInput
-
On entry: , the number of right-hand sides, i.e., the number of columns of the matrices and .
Constraint:
.
- 4: – Real (Kind=nag_wp) arrayInput/Output
-
Note: the second dimension of the array
a
must be at least
.
On entry: the by coefficient matrix .
On exit: the contents of
a are destroyed.
- 5: – IntegerInput
-
On entry: the first dimension of the array
a as declared in the (sub)program from which
f08kcf (dgelsd) is called.
Constraint:
.
- 6: – Real (Kind=nag_wp) arrayInput/Output
-
Note: the second dimension of the array
b
must be at least
.
On entry: the by right-hand side matrix .
On exit:
b is overwritten by the
by
solution matrix
. If
and
, the residual sum of squares for the solution in the
th column is given by the sum of squares of elements
in that column.
- 7: – IntegerInput
-
On entry: the first dimension of the array
b as declared in the (sub)program from which
f08kcf (dgelsd) is called.
Constraint:
.
- 8: – Real (Kind=nag_wp) arrayOutput
-
Note: the dimension of the array
s
must be at least
.
On exit: the singular values of in decreasing order.
- 9: – Real (Kind=nag_wp)Input
-
On entry: used to determine the effective rank of . Singular values are treated as zero. If , machine precision is used instead.
- 10: – IntegerOutput
-
On exit: the effective rank of , i.e., the number of singular values which are greater than .
- 11: – Real (Kind=nag_wp) arrayWorkspace
-
On exit: if
,
contains the minimum value of
lwork required for optimal performance.
- 12: – IntegerInput
-
On entry: the dimension of the array
work as declared in the (sub)program from which
f08kcf (dgelsd) is called.
The exact minimum amount of workspace needed depends on
m,
n and
nrhs. As long as
lwork is at least
where
is equal to the maximum size of the subproblems at the bottom of the computation tree (usually about
),
and
, the code will execute correctly.
If
, a workspace query is assumed; the routine only calculates the optimal size of the
work array and the minimum size of the
iwork array, and returns these values as the first entries of the
work and
iwork arrays, and no error message related to
lwork is issued.
Suggested value:
for optimal performance,
lwork should generally be larger than the minimum required as set out above. Consider increasing
lwork by at least
, where
is the optimal
block size.
Constraint:
.
- 13: – Integer arrayWorkspace
-
Note: the dimension of the array
iwork
must be at least
, where
is at least
.
On exit: if , returns the minimum .
- 14: – IntegerOutput
On exit:
unless the routine detects an error (see
Section 6).
6
Error Indicators and Warnings
If , argument had an illegal value. An explanatory message is output, and execution of the program is terminated.
-
The algorithm for computing the SVD failed to converge; off-diagonal elements of an intermediate bidiagonal form did not converge to zero.
7
Accuracy
See Section 4.5 of
Anderson et al. (1999) for details.
8
Parallelism and Performance
f08kcf (dgelsd) is threaded by NAG for parallel execution in multithreaded implementations of the NAG Library.
f08kcf (dgelsd) makes calls to BLAS and/or LAPACK routines, which may be threaded within the vendor library used by this implementation. Consult the documentation for the vendor library for further information.
Please consult the
X06 Chapter Introduction for information on how to control and interrogate the OpenMP environment used within this routine. Please also consult the
Users' Note for your implementation for any additional implementation-specific information.
The complex analogue of this routine is
f08kqf (zgelsd).
10
Example
This example solves the linear least squares problem
for the solution,
, of minimum norm, where
A tolerance of is used to determine the effective rank of .
Note that the block size (NB) of assumed in this example is not realistic for such a small problem, but should be suitable for large problems.
10.1
Program Text
Program Text (f08kcfe.f90)
10.2
Program Data
Program Data (f08kcfe.d)
10.3
Program Results
Program Results (f08kcfe.r)