LSGRR

TOL — Scalar containing the tolerance used to determine when a singular value (from the singular value decomposition of A) is negligible. (Input)
If TOL is positive, then a singular value σi considered negligible if σi ≤ TOL . If TOL is negative, then a singular value σi considered negligible if σi ≤ ∣TOL∣ * ∥A∥∞. In this case, |TOL| generally contains an estimate of the level of the relative error in the data.
Default: TOL = 1.0e-5 for single precision and 1.0d-10 for double precision.

IRANK — Scalar containing an estimate of the rank of A. (Output)

LDGINV — Leading dimension of GINVA exactly as specified in the dimension statement of the calling program. (Input)
Default: LDGINV = size (GINV,1).

FORTRAN 90 Interface

Generic: CALL LSGRR (A, GINVA [, …])

Specific: The specific interface names are S_LSGRR and D_LSGRR.

FORTRAN 77 Interface

Single: CALL LSGRR (NRA, NCA, A, LDA, TOL, IRANK, GINVA, LDGINV)

Double: The double precision name is DLSGRR.

ScaLAPACK Interface

Generic: CALL LSGRR (A0, GINVA0 [, …])

Specific: The specific interface names are S_LSGRR and D_LSGRR.

See the ScaLAPACK Usage Notes below for a description of the arguments for distributed computing.

Description

Let k = IRANK, the rank of A; let n = NRA, the number of rows in A; let p = NCA, the number of columns in A; and let

be the generalized inverse of A.

To compute the Moore-Penrose generalized inverse, the routine LSVRR is first used to compute the singular value decomposition of A. A singular value decomposition of A consists of an n × n orthogonal matrix U, a p × p orthogonal matrix V and a diagonal matrix ∑ = diag(σ1,…, σm), m = min(n, p), such that UT AV = [∑, 0] if n ≤ p and UT AV = [∑, 0]T if n ≥ p. Only the first p columns of U are computed. The rank k is estimated by counting the number of nonnegligible σi.

The matrices U and V can be partitioned as U = (U1, U2) and V = (V1, V2) where both U1 and V1 are k × k matrices. Let∑1 = diag(σ1, …, σk). The Moore-Penrose generalized inverse of A is

A† = V1∑1-1U1T

The underlying code of routine LSGRR is based on either LINPACK, LAPACK, or ScaLAPACK code depending upon which supporting libraries are used during linking. For a detailed explanation see Using ScaLAPACK, LAPACK, LINPACK, and EISPACK in the Introduction section of this manual.

Comments

1. Workspace may be explicitly provided, if desired, by use of L2GRR/DL2GRR. The reference is:

CALL L2GRR (NRA, NCA, A, LDA, TOL, IRANK, GINVA, LDGINV, WKA, WK)

The additional arguments are as follows:

WKA — Work vector of length NRA * NCA used as workspace for the matrix A. If A is not needed, WKA and A can share the same storage locations.

WK — Work vector of length LWK where LWK is equal to NRA2 + NCA2 + min(NRA + 1, NCA) + NRA + NCA + max(NRA, NCA) − 2.

2. Informational error

Type	Code	Description
4	1	Convergence cannot be achieved for all the singular values and their corresponding singular vectors.

ScaLAPACK Usage Notes

The arguments which differ from the standard version of this routine are:

A0 — MXLDA by MXCOL local matrix containing the local portions of the distributed matrix A. A contains the matrix for which the generalized inverse is to be computed. (Input)

GINVA0 — MXLDG by MXCOLG local matrix containing the local portions of the distributed matrix GINVA. GINVA contains the generalized inverse of matrix A. (Output)

All other arguments are global and are the same as described for the standard version of the routine. In the argument descriptions above, MXLDA, MXCOL, MXLDG, and MXCOLG can be obtained through a call to SCALAPACK_GETDIM (see Utilities) after a call to SCALAPACK_SETUP (see Chapter 11, ”Utilities”) has been made. See the ScaLAPACK Example below.

Examples

Example

This example computes the generalized inverse of a 3 × 2 matrix A. The rank k = IRANK and the inverse

are printed.

USE IMSL_LIBRARIES

! Declare variables

PARAMETER (NRA=3, NCA=2, LDA=NRA, LDGINV=NCA)

REAL A(LDA,NCA), GINV(LDGINV,NRA)

! Set values for A

! A = ( 1 0 )

! ( 1 1 )

! ( 100 -50 )

DATA A/1., 1., 100., 0., 1., -50./

! Compute generalized inverse

TOL = AMACH(4)

TOL = 10.*TOL

CALL LSGRR (A, GINV,TOL=TOL, IRANK=IRANK)

! Print results

CALL UMACH (2, NOUT)

WRITE (NOUT, *) ’IRANK = ’, IRANK

CALL WRRRN (’GINV’, GINV)

END

Output

IRANK = 2

GINV

1 2 3

1 0.1000 0.3000 0.0060

2 0.2000 0.6000 -0.0080

ScaLAPACK Example

This example computes the generalized inverse of a 6 × 4 matrix A as a distributed example. The rank k = IRANK and the inverse

are printed.

USE MPI_SETUP_INT

USE IMSL_LIBRARIES

USE SCALAPACK_SUPPORT

IMPLICIT NONE

INCLUDE ‘mpif.h’

! Declare variables

INTEGER IRANK, LDA, NCA, NRA, DESCA(9), DESCG(9), &

LDGINV, MXLDG, MXCOLG, NOUT

INTEGER INFO, MXCOL, MXLDA

REAL TOL, AMACH

REAL, ALLOCATABLE :: A(:,:),GINVA(:,:)

REAL, ALLOCATABLE :: A0(:,:), GINVA0(:,:)

PARAMETER (NRA=6, NCA=4, LDA=NRA, LDGINV=NCA)

! Set up for MPI

MP_NPROCS = MP_SETUP()

IF(MP_RANK .EQ. 0) THEN

ALLOCATE (A(LDA,NCA), GINVA(NCA,NRA))

! Set values for A

A(1,:) = (/ 1.0, 2.0, 1.0, 4.0/)

A(2,:) = (/ 3.0, 2.0, 1.0, 3.0/)

A(3,:) = (/ 4.0, 3.0, 1.0, 4.0/)

A(4,:) = (/ 2.0, 1.0, 3.0, 1.0/)

A(5,:) = (/ 1.0, 5.0, 2.0, 2.0/)

A(6,:) = (/ 1.0, 2.0, 2.0, 3.0/)

ENDIF

! Set up a 1D processor grid and define

! its context ID, MP_ICTXT

CALL SCALAPACK_SETUP(NRA, NCA, .TRUE., .TRUE.)

! Get the array descriptor entities MXLDA,

! MXCOL, MXLDG, and MXCOLG

CALL SCALAPACK_GETDIM(NRA, NCA, MP_MB, MP_NB, MXLDA, MXCOL)

CALL SCALAPACK_GETDIM(NCA, NRA, MP_NB, MP_MB, MXLDG, MXCOLG)

! Set up the array descriptors

CALL DESCINIT(DESCA, NRA, NCA, MP_MB, MP_NB, 0, 0, MP_ICTXT, MXLDA, &

INFO)

CALL DESCINIT(DESCG, NCA, NRA, MP_NB, MP_MB, 0, 0, MP_ICTXT, MXLDG, &

INFO)

! Allocate space for the local arrays

ALLOCATE (A0(MXLDA,MXCOL), GINVA0(MXLDG,MXCOLG))

! Map input array to the processor grid

CALL SCALAPACK_MAP(A, DESCA, A0)

! Compute the generalized inverse

TOL = AMACH(4)

TOL = 10. * TOL

CALL LSGRR (A0, GINVA0, TOL=TOL, IRANK=IRANK)

! Unmap the results from the distributed

! array back to a non-distributed array.

! After the unmap, only Rank=0 has the full

! array.

CALL SCALAPACK_UNMAP(GINVA0, DESCG, GINVA)

! Print results.

! Only Rank=0 has the solution, GINVA

IF(MP_RANK .EQ. 0) THEN

CALL UMACH (2, NOUT)

WRITE (NOUT, *) ‘IRANK = ‘,IRANK

CALL WRRRN (‘GINVA’, GINVA)

ENDIF

! Exit ScaLAPACK usage

CALL SCALAPACK_EXIT(MP_ICTXT)

! Shut down MPI

MP_NPROCS = MP_SETUP(‘FINAL’)

END