LSGRR

Computes the generalized inverse of a real matrix.

Required Arguments

ANRA by NCA matrix whose generalized inverse is to be computed.   (Input)

GINVANCA 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(np), 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. Let1 = diag(σ1, , σk). The Moore-Penrose generalized inverse of A is

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

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  AA 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  GINVAGINVA  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.

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

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 


Visual Numerics, Inc.
Visual Numerics - Developers of IMSL and PV-WAVE
http://www.vni.com/
PHONE: 713.784.3131
FAX:713.781.9260