LSGRR
Computes the generalized inverse of a real matrix.
Required Arguments
A — NRA by NCA matrix whose generalized inverse is to be computed. (Input)
GINVA — NCA by NRA matrix containing the generalized inverse of A. (Output)
Optional Arguments
NRA — Number of rows in the matrix A. (Input)
Default: NRA = size (A,1).
NCA — Number of columns in the matrix A. (Input)
Default: NCA = size (A,2).
LDA — Leading dimension of A exactly as specified in the dimension statement of the calling program. (Input)
Default: LDA = size (A,1).
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
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
IRANK = 2
GINV
1 2 3
1 0.1000 0.3000 0.0060
2 0.2000 0.6000 -0.0080
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
