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
NORM_CHOICE	CONDITION Number	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 matrices is 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

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

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