DECOMP Subroutine

public subroutine DECOMP(N, K, Eta, Tol, Irank, Insing)

NAME

decomp(3f) - [M_datapac:STATISTICS] decomposes a weighted data matrix
(utility routine used by other routines)

SYNOPSIS

   SUBROUTINE DECOMP(N,K,Eta,Tol,Irank,Insing)

    INTEGER :: N
    INTEGER :: K
    REAL(kind=wp) :: Eta
    REAL(kind=wp) :: Tol
    INTEGER :: Irank
    INTEGER :: Insing

DESCRIPTION

DECOMP(3f) decomposes the weighted data matrix q which originally =
the N by K data matrix x times the square root of the weights (in w).

The original q is decomposed into a new q times the inverse of a
diagonal matrix d times the diagonal matrix d times an upper triangular
matrix r.

The new N by K q has orthogonal columns.

A second output from DECOMP(3f) is the rank and status (non-singular
or singular) of the data matrix X.

A third output from DECOMP(3f) is the numerically optimal pivot points
for the decomposition.

OPTIONS

 X   description of parameter
 Y   description of parameter

EXAMPLES

Sample program:

program demo_decomp
use M_datapac, only : decomp
implicit none
! call decomp(x,y)
end program demo_decomp

Results:

AUTHOR

The original DATAPAC library was written by James Filliben of the Statistical
Engineering Division, National Institute of Standards and Technology.

MAINTAINER

John Urban, 2022.05.31

LICENSE

CC0-1.0

Arguments

Type IntentOptional Attributes Name
integer :: N
integer :: K
real(kind=wp) :: Eta
real(kind=wp) :: Tol
integer :: Irank
integer :: Insing

Common Blocks

CAUPLT (subroutine)
CHSCDF (subroutine)
CODE (subroutine)
DEMOD (subroutine)
DEXPLT (subroutine)
EV1PLT (subroutine)
EV2PLT (subroutine)
EXPPLT (subroutine)
EXTREM (subroutine)
EXTREM (subroutine)
FREQ (subroutine)
GAMPLT (subroutine)
GEOPLT (subroutine)
HFNPLT (subroutine)
INVXWX (subroutine)
LAMPLT (subroutine)
LGNPLT (subroutine)
LOGPLT (subroutine)
MEDIAN (subroutine)
norout (subroutine)
NORPLT (subroutine)
PARPLT (subroutine)
PLOTU (subroutine)
POIPLT (subroutine)
RUNS (subroutine)
SAMPP (subroutine)
SPCORR (subroutine)
TAIL (subroutine)
TPLT (subroutine)
TRIM (subroutine)
UNIPLT (subroutine)
WEIB (subroutine)
WEIPLT (subroutine)
WIND (subroutine)
"> common /BLOCK2_real32/

Type Attributes Name Initial
real :: WS(15000)

INVXWX (subroutine)
"> common /BLOCK3_real32/

Type Attributes Name Initial
real :: DUM1(3000)
real :: DUM2(3000)

Source Code

SUBROUTINE DECOMP(N,K,Eta,Tol,Irank,Insing)
INTEGER :: N
INTEGER :: K
REAL(kind=wp) :: Eta
REAL(kind=wp) :: Tol
INTEGER :: Irank
INTEGER :: Insing

REAL(kind=wp) :: D, dis, dn, DUM1, DUM2, hold, Q, R, risj, tol2, WS
INTEGER i, ip, IPIvot, iqarg, iqarg1, iqarg2, irarg, irarg1, irarg2, is, ism1, isp1, j, l, m
LOGICAL fsum
DIMENSION Q(10000) , R(2500) , D(50) , IPIvot(50)
COMMON /BLOCK2_real32/ WS(15000)
COMMON /BLOCK3_real32/ DUM1(3000) , DUM2(3000)
EQUIVALENCE (Q(1),WS(1))          !     Q--USED AND CHANGED
EQUIVALENCE (R(1),WS(10001))      !     R--DEFINED
EQUIVALENCE (D(1),WS(12501))      !     D--PERMANENTLY DEFINED
EQUIVALENCE (IPIvot(1),WS(12551)) !     IPIVOT--PERMANENTLY DEFINED
!
!-----START POINT-----------------------------------------------------
!
!     ZERO-OUT SOME VARIABLES, VECTORS, AND ARRAYS
!
      Insing = 0
      Irank = 0
      DO j = 1 , K
         D(j) = 0.0_wp
         DO i = 1 , K
            irarg = (i-1)*K + j
            R(irarg) = 0.0_wp
         ENDDO
      ENDDO
!
      tol2 = Tol*Tol
      DO j = 1 , K
         IPIvot(j) = j
      ENDDO
      DO is = 1 , K
!
!     BEGIN STEP NUMBER      IS      IN THE DECOMPOSITION
!
         IF ( is==1 ) fsum = .TRUE.
 50      dis = 0.0_wp
         ip = is
!
!     BEGIN THE PIVOT SEARCH
!
         DO j = is , K
            m = IPIvot(j)
            IF ( fsum ) THEN
               DO l = 1 , N
                  iqarg = (l-1)*K + m
                  DUM1(l) = Q(iqarg)
                  DUM2(l) = Q(iqarg)
               ENDDO
!
               CALL DOT(DUM1,DUM2,1,N,0.0_wp,D(j))
            ENDIF
!
            IF ( dis<D(j) ) THEN
               dis = D(j)
               ip = j
            ENDIF
         ENDDO
!
!     END THE PIVOT SEARCH
!
         m = IPIvot(ip)
         IF ( fsum ) dn = dis
         IF ( dis<Eta*dn ) THEN
            fsum = .TRUE.
         ELSE
            fsum = .FALSE.
         ENDIF
         IF ( fsum ) GOTO 50
         IF ( ip/=is ) THEN
!
!     BEGIN COLUMN INTERCHANGES
!
            D(ip) = D(is)
            IPIvot(ip) = IPIvot(is)
            IPIvot(is) = m
            IF ( is/=1 ) THEN
               ism1 = is - 1
               DO i = 1 , ism1
                  irarg1 = (i-1)*K + ip
                  irarg2 = (i-1)*K + is
                  hold = R(irarg1)
                  R(irarg1) = R(irarg2)
                  R(irarg2) = hold
               ENDDO
            ENDIF
         ENDIF
!
!     END COLUMN INTERCHANGES
!
         DO l = 1 , N
            iqarg = (l-1)*K + m
            DUM1(l) = Q(iqarg)
            DUM2(l) = Q(iqarg)
         ENDDO
!
         CALL DOT(DUM1,DUM2,1,N,0.0_wp,D(is))
!
         dis = D(is)
         IF ( dis<=tol2*D(1) ) RETURN
         IF ( dis/=0.0_wp ) THEN
            isp1 = is + 1
            IF ( isp1<=K ) THEN
!
!     BEGIN ORTHOGONALIZATION
!
               DO j = isp1 , K
                  ip = IPIvot(j)
                  DO l = 1 , N
                     iqarg1 = (l-1)*K + m
                     iqarg2 = (l-1)*K + ip
                     DUM1(l) = Q(iqarg1)
                     DUM2(l) = Q(iqarg2)
                  ENDDO
!
                  irarg = (is-1)*K + j
                  CALL DOT(DUM1,DUM2,1,N,0.0_wp,R(irarg))
                  R(irarg) = R(irarg)/dis
!
                  risj = R(irarg)
                  DO i = 1 , N
                     iqarg1 = (i-1)*K + ip
                     iqarg2 = (i-1)*K + m
                     Q(iqarg1) = Q(iqarg1) - risj*Q(iqarg2)
                  ENDDO
                  D(j) = D(j) - dis*risj*risj
               ENDDO
            ENDIF
!
!     END ORTHOGONALIZATION
!
            Irank = is
         ELSE
            Insing = 0
            RETURN
         ENDIF
      ENDDO
!
!     END STEP NUMBER     IS     INTHE DECOMPOSITION
!
      Insing = 1
END SUBROUTINE DECOMP