COND
Computes the condition number of a matrix.
Function Return Value
Computes condition number of matrix A. This is a scalar for the case where A is rank-2 or a sparse matrix. It is a rank-1 array when A is a dense rank-3 array. (Output)
Required Argument
A — Matrix for which the condition number is to be computed. The matrix may be real, double, complex, double-complex, or one of the computational sparse matrix derived types, ?_hbc_sparse. For an array of type real, double, complex, or double-complex the array may be of rank-2 or rank-3.
For a dense rank-3 array, each rank-2 array section (for fixed third subscript) is a separate problem. In this case, the output is a rank-1 array of condition numbers for each problem. (Input)
Optional Arguments, Packaged Options
NORM_CHOICE — Integer indicating the type of norm to be used in computing the condition number.
NORM_CHOICE | CONDITION Number | Square Matrix | Rectangular Matrix | ||
|---|---|---|---|---|---|
Dense | Sparse | Dense | Sparse | ||
1 | L1 | Yes | Yes | No | No |
2 (Default) | L2 | Yes | Yes | Yes | No |
huge(1) | L∞ | Yes | Yes | No | No |
This function uses LIN_SOL_SVD (see Chapter 1, “Linear Systems”).
The option and derived type names are given in the following tables:
Option Names for COND | Option Value |
|---|---|
?_cond_set_small | 1 |
?_cond_for_lin_sol_svd | 2 |
Name of Unallocated Option Array to Use for Setting Options | Use | Derived Type |
|---|---|---|
?_cond_options(:) | Use when setting options for calls hereafter. | ?_options |
?_cond_options_once(:) | Use when setting options for next call only. | ?_options |
For a description on how to use these options, see Matrix Optional Data Changes. See LIN_SOL_SVD in Chapter 1, “Linear Systems” for the specific options for these routines.
FORTRAN 90 Interface
COND (A [, …])
Description
The mathematical definitions of the condition numbers which this routine estimates are:
COND can be used with either dense or sparse matrices as follows:
Square Matrix | Rectangular Matrix | |||
|---|---|---|---|---|
Dense | Sparse | Dense | Sparse | |
L1 | Yes | Yes | No | No |
L2 | Yes | Yes | Yes | No |
L∞ | Yes | Yes | No | No |
The generic function COND can be used with either dense or sparse square matrices. This function uses LIN_SOL_SVD for dense square and rectangular matrices in computing κ2(A) = s1/sn. The function uses LIN_SOL_GEN for dense square matrices in computing κ1(A) andκ∞(A). For sparse square matrices, the values returned for κ1(A) andκ∞(A) are provided by the SuperLU linear equation solver. The condition number κ2(A) = s1/sn is computed by an algorithm that first approximates s1 by computing the singular values of the κ × κ bidiagonal matrix obtained using the Lanczos method found in Golub and Van Loan, Ed. 3, p. 495. Here κ is set using the value A%Options%Cond_Iteration_Max, which has the default value of 30. The value sn-2 is obtained using the power method, Golub and Van Loan, p. 330, iterating with the inverse matrix. For complex matricesis replaced by. The dominant eigenvalue of this inverse matrix is sn-2. The number of iterations is limited by the parameter value κ or relative accuracy equal to the cube root of machine epsilon. Some timing tests indicate that computing κ2(A) for sparse matrices by this algorithm typically requires about twice the time as for a single linear solve using the defined operator A .ix. b.
For computation of κ2(A) with rectangular sparse matrices one can use a dense matrix representation for the matrix. This is not recommended except for small problem sizes. For overdetermined systems of sparse least-squares equations Ax ≅ b a related square system is given by
One can form C, which has more than twice the number of non-zeros as A. But C is still sparse. One can use the condition number of C as an estimate of the accuracy for the solution vector x and the residual vector r. Note that this version of the condition number is not the same as the l2 condition number of A but is relevant to determining the accuracy of the least-squares system.
Examples
Dense Matrix Example (operator_ex02.f90)
use wrrrn_int
use linear_operators
integer, parameter :: N=3
real (kind(1.e0)) A(N,N)
real (kind(1.e0)) C1, C2, CINF
DATA A/2.0, 2.0, -4.0, 0.0, -1.0, 2.0, 0.0, 0.0, 5.0/
CINF = COND (A, norm_choice=huge(1))
C1 = COND (A, norm_choice=1)
C2 = COND (A)
call wrrrn ( 'A', A)
write (*,*) 'L1 condition number= ', C1
write (*,*) 'L2 condition number= ', C2
write (*,*) 'L infinity condition number= ', CINF
end
Output
A
1 2 3
1 2.000 0.000 0.000
2 2.000 -1.000 0.000
3 -4.000 2.000 5.000
L1 condition number= 12.0
L2 condition number= 10.405088
L infinity condition number= 22.0
Sparse Matrix Example
use wrrrn_int
use linear_operators
type (s_sparse) S
type (s_hbc_sparse) H
integer, parameter :: N=3
real (kind(1.e0)) X(N,N)
real (kind(1.e0)) C1, C2, CINF
S = s_entry (1, 1, 2.0)
S = s_entry (2, 1, 2.0)
S = s_entry (3, 1, -4.0)
S = s_entry (3, 2, 2.0)
S = s_entry (2, 2, -1.0)
S = s_entry (3, 3, 5.0)
H = S ! sparse
X = H ! dense equivalent of H
CINF = COND (H, norm_choice=huge(1))
C1 = COND (H, norm_choice=1)
C2 = COND (H)
call wrrrn ( 'H', X)
write (*,*) 'L1 condition number= ', C1
write (*,*) 'L2 condition number= ', C2
write (*,*) 'L infinity condition number= ', CINF
end
Output
H
1 2 3
1 2.000 0.000 0.000
2 2.000 -1.000 0.000
3 -4.000 2.000 5.000
L1 condition number= 12.0
L2 condition number= 10.405088
L infinity condition number= 22.0
Parallel Example (parallel_ex02.f90)
use linear_operators
use mpi_setup_int
implicit none
! This is the equivalent of Parallel Example 2 for .i. and det() with box
! data types, operators and functions.
integer, parameter :: n=32, nr=4
integer J
real(kind(1e0)) :: one=1e0
real(kind(1e0)), dimension(nr) :: err, det_A, det_i
real(kind(1e0)), dimension(n,n,nr) :: A, inv, R, S
! Setup for MPI.
MP_NPROCS=MP_SETUP()
! Generate a random matrix.
A = rand(A)
! Compute the matrix inverse and its determinant.
inv = .i.A; det_A = det(A)
! Compute the determinant for the inverse matrix.
det_i = det(inv)
! Check the quality of both left and right inverses.
DO J=1,nr; R(:,:,J)=EYE(N); END DO
S=R; R=R-(A .x. inv); S=S-(inv .x. A)
err = (norm(R)+norm(S))/cond(A)
if (ALL(err <= sqrt(epsilon(one)) .and. &
abs(det_A*det_i - one) <= sqrt(epsilon(one)))&
.and. MP_RANK == 0) &
write (*,*) 'Parallel Example 2 is correct.'
! See to any error messages and quit MPI.
MP_NPROCS=MP_SETUP('Final')
end
