relibc_openlibm/slatec/dlssud.f

*DECK DLSSUD
      SUBROUTINE DLSSUD (A, X, B, N, M, NRDA, U, NRDU, IFLAG, MLSO,
     +   IRANK, ISCALE, Q, DIAG, KPIVOT, S, DIV, TD, ISFLG, SCALES)
C***BEGIN PROLOGUE  DLSSUD
C***SUBSIDIARY
C***PURPOSE  Subsidiary to DBVSUP and DSUDS
C***LIBRARY   SLATEC
C***TYPE      DOUBLE PRECISION (LSSUDS-S, DLSSUD-D)
C***AUTHOR  Watts, H. A., (SNLA)
C***DESCRIPTION
C
C    DLSSUD solves the underdetermined system of equations  A Z = B,
C    where A is N by M and N .LE. M.  In particular, if rank A equals
C    IRA, a vector X and a matrix U are determined such that X is the
C    UNIQUE solution of smallest length, satisfying A X = B, and the
C    columns of U form an orthonormal basis for the null space of A,
C    satisfying A U = 0 .  Then all solutions Z are given by
C              Z = X + C(1)*U(1) + ..... + C(M-IRA)*U(M-IRA)
C    where U(J) represents the J-th column of U and the C(J) are
C    arbitrary constants.
C    If the system of equations are not compatible, only the least
C    squares solution of minimal length is computed.
C
C *********************************************************************
C   INPUT
C *********************************************************************
C
C     A -- Contains the matrix of N equations in M unknowns, A remains
C          unchanged, must be dimensioned NRDA by M.
C     X -- Solution array of length at least M.
C     B -- Given constant vector of length N, B remains unchanged.
C     N -- Number of equations, N greater or equal to 1.
C     M -- Number of unknowns, M greater or equal to N.
C     NRDA -- Row dimension of A, NRDA greater or equal to N.
C     U -- Matrix used for solution, must be dimensioned NRDU by
C          (M - rank of A).
C          (storage for U may be ignored when only the minimal length
C           solution X is desired)
C     NRDU -- Row dimension of U, NRDU greater or equal to M.
C             (if only the minimal length solution is wanted,
C              NRDU=0 is acceptable)
C     IFLAG -- Status indicator
C           =0  for the first call (and for each new problem defined by
C               a new matrix A) when the matrix data is treated as exact
C           =-K for the first call (and for each new problem defined by
C               a new matrix A) when the matrix data is assumed to be
C               accurate to about K digits.
C           =1  for subsequent calls whenever the matrix A has already
C               been decomposed (problems with new vectors B but
C               same matrix A can be handled efficiently).
C     MLSO -- =0 if only the minimal length solution is wanted.
C             =1 if the complete solution is wanted, includes the
C                linear space defined by the matrix U.
C     IRANK -- Variable used for the rank of A, set by the code.
C     ISCALE -- Scaling indicator
C               =-1 if the matrix A is to be pre-scaled by
C               columns when appropriate.
C               If the scaling indicator is not equal to -1
C               no scaling will be attempted.
C            For most problems scaling will probably not be necessary.
C     Q -- Matrix used for the transformation, must be dimensioned
C            NRDA by M.
C     DIAG,KPIVOT,S, -- Arrays of length at least N used for internal
C      DIV,TD,SCALES    storage (except for SCALES which is M).
C     ISFLG -- Storage for an internal variable.
C
C *********************************************************************
C   OUTPUT
C *********************************************************************
C
C     IFLAG -- Status indicator
C            =1 if solution was obtained.
C            =2 if improper input is detected.
C            =3 if rank of matrix is less than N.
C               To continue, simply reset IFLAG=1 and call DLSSUD again.
C            =4 if the system of equations appears to be inconsistent.
C               However, the least squares solution of minimal length
C               was obtained.
C     X -- Minimal length least squares solution of A Z = B
C     IRANK -- Numerically determined rank of A, must not be altered
C              on succeeding calls with input values of IFLAG=1.
C     U -- Matrix whose M-IRANK columns are mutually orthogonal unit
C          vectors which span the null space of A. This is to be ignored
C          when MLSO was set to zero or IFLAG=4 on output.
C     Q -- Contains the strictly upper triangular part of the reduced
C           matrix and transformation information.
C     DIAG -- Contains the diagonal elements of the triangular reduced
C             matrix.
C     KPIVOT -- Contains the pivotal information.  The row interchanges
C               performed on the original matrix are recorded here.
C     S -- Contains the solution of the lower triangular system.
C     DIV,TD -- Contains transformation information for rank
C               deficient problems.
C     SCALES -- Contains the column scaling parameters.
C
C *********************************************************************
C
C***SEE ALSO  DBVSUP, DSUDS
C***REFERENCES  H. A. Watts, Solving linear least squares problems
C                 using SODS/SUDS/CODS, Sandia Report SAND77-0683,
C                 Sandia Laboratories, 1977.
C***ROUTINES CALLED  D1MACH, DDOT, DOHTRL, DORTHR, J4SAVE, XERMAX,
C                    XERMSG, XGETF, XSETF
C***REVISION HISTORY  (YYMMDD)
C   750601  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
C   900328  Added TYPE section.  (WRB)
C   910408  Updated the AUTHOR and REFERENCES sections.  (WRB)
C   920501  Reformatted the REFERENCES section.  (WRB)
C***END PROLOGUE  DLSSUD
      INTEGER J4SAVE
      DOUBLE PRECISION DDOT, D1MACH
      INTEGER I, IFLAG, IRANK, IRP, ISCALE, ISFLG, J, JR, K, KP,
     1     KPIVOT(*), L, M, MAXMES, MJ, MLSO, N, NFAT, NFATAL, NMIR,
     2     NRDA, NRDU, NU
      DOUBLE PRECISION A(NRDA,*), B(*), DIAG(*), DIV(*), GAM, GAMMA,
     1     Q(NRDA,*), RES, S(*), SCALES(*), SS, TD(*), U(NRDU,*), URO,
     2     X(*)
C
C     ******************************************************************
C
C          MACHINE PRECISION (COMPUTER UNIT ROUNDOFF VALUE) IS DEFINED
C          BY THE FUNCTION D1MACH.
C
C     ******************************************************************
C
C     BEGIN BLOCK PERMITTING ...EXITS TO 310
C        BEGIN BLOCK PERMITTING ...EXITS TO 80
C***FIRST EXECUTABLE STATEMENT  DLSSUD
            URO = D1MACH(4)
C
            IF (N .LT. 1 .OR. M .LT. N .OR. NRDA .LT. N) GO TO 70
            IF (NRDU .NE. 0 .AND. NRDU .LT. M) GO TO 70
               IF (IFLAG .GT. 0) GO TO 60
C
                  CALL XGETF(NFATAL)
                  MAXMES = J4SAVE(4,0,.FALSE.)
                  ISFLG = -15
                  IF (IFLAG .EQ. 0) GO TO 10
                     ISFLG = IFLAG
                     NFAT = -1
                     IF (NFATAL .EQ. 0) NFAT = 0
                     CALL XSETF(NFAT)
                     CALL XERMAX(1)
   10             CONTINUE
C
C                 COPY MATRIX A INTO MATRIX Q
C
                  DO 30 K = 1, M
                     DO 20 J = 1, N
                        Q(J,K) = A(J,K)
   20                CONTINUE
   30             CONTINUE
C
C                 USE ORTHOGONAL TRANSFORMATIONS TO REDUCE Q TO LOWER
C                 TRIANGULAR FORM
C
                  CALL DORTHR(Q,N,M,NRDA,IFLAG,IRANK,ISCALE,DIAG,KPIVOT,
     1                        SCALES,DIV,TD)
C
                  CALL XSETF(NFATAL)
                  CALL XERMAX(MAXMES)
                  IF (IRANK .EQ. N) GO TO 40
C
C                    FOR RANK DEFICIENT PROBLEMS USE ADDITIONAL
C                    ORTHOGONAL TRANSFORMATIONS TO FURTHER REDUCE Q
C
                     IF (IRANK .NE. 0)
     1                  CALL DOHTRL(Q,N,NRDA,DIAG,IRANK,DIV,TD)
C     ...............EXIT
                     GO TO 310
   40             CONTINUE
C
C                 STORE DIVISORS FOR THE TRIANGULAR SOLUTION
C
                  DO 50 K = 1, N
                     DIV(K) = DIAG(K)
   50             CONTINUE
C        .........EXIT
                  GO TO 80
   60          CONTINUE
C        ......EXIT
               IF (IFLAG .EQ. 1) GO TO 80
   70       CONTINUE
C
C           INVALID INPUT FOR DLSSUD
            IFLAG = 2
            CALL XERMSG ('SLATEC', 'DLSSUD',
     +         'INVALID IMPUT PARAMETERS.', 2, 1)
C     ......EXIT
            GO TO 310
   80    CONTINUE
C
C
         IF (IRANK .GT. 0) GO TO 130
C
C           SPECIAL CASE FOR THE NULL MATRIX
            DO 110 K = 1, M
               X(K) = 0.0D0
               IF (MLSO .EQ. 0) GO TO 100
                  U(K,K) = 1.0D0
                  DO 90 J = 1, M
                     IF (J .NE. K) U(J,K) = 0.0D0
   90             CONTINUE
  100          CONTINUE
  110       CONTINUE
            DO 120 K = 1, N
               IF (B(K) .GT. 0.0D0) IFLAG = 4
  120       CONTINUE
         GO TO 300
  130    CONTINUE
C           BEGIN BLOCK PERMITTING ...EXITS TO 180
C
C              COPY CONSTANT VECTOR INTO S AFTER FIRST INTERCHANGING
C              THE ELEMENTS ACCORDING TO THE PIVOTAL SEQUENCE
C
               DO 140 K = 1, N
                  KP = KPIVOT(K)
                  X(K) = B(KP)
  140          CONTINUE
               DO 150 K = 1, N
                  S(K) = X(K)
  150          CONTINUE
C
               IRP = IRANK + 1
               NU = 1
               IF (MLSO .EQ. 0) NU = 0
C           ...EXIT
               IF (IRANK .EQ. N) GO TO 180
C
C              FOR RANK DEFICIENT PROBLEMS WE MUST APPLY THE
C              ORTHOGONAL TRANSFORMATION TO S
C              WE ALSO CHECK TO SEE IF THE SYSTEM APPEARS TO BE
C              INCONSISTENT
C
               NMIR = N - IRANK
               SS = DDOT(N,S(1),1,S(1),1)
               DO 170 L = 1, IRANK
                  K = IRP - L
                  GAM = ((TD(K)*S(K)) + DDOT(NMIR,Q(IRP,K),1,S(IRP),1))
     1                  /(TD(K)*DIV(K))
                  S(K) = S(K) + GAM*TD(K)
                  DO 160 J = IRP, N
                     S(J) = S(J) + GAM*Q(J,K)
  160             CONTINUE
  170          CONTINUE
               RES = DDOT(NMIR,S(IRP),1,S(IRP),1)
C           ...EXIT
               IF (RES
     1             .LE. SS*(10.0D0*MAX(10.0D0**ISFLG,10.0D0*URO))**2)
     2            GO TO 180
C
C              INCONSISTENT SYSTEM
               IFLAG = 4
               NU = 0
  180       CONTINUE
C
C           APPLY FORWARD SUBSTITUTION TO SOLVE LOWER TRIANGULAR SYSTEM
C
            S(1) = S(1)/DIV(1)
            IF (IRANK .LT. 2) GO TO 200
            DO 190 K = 2, IRANK
               S(K) = (S(K) - DDOT(K-1,Q(K,1),NRDA,S(1),1))/DIV(K)
  190       CONTINUE
  200       CONTINUE
C
C           INITIALIZE X VECTOR AND THEN APPLY ORTHOGONAL TRANSFORMATION
C
            DO 210 K = 1, M
               X(K) = 0.0D0
               IF (K .LE. IRANK) X(K) = S(K)
  210       CONTINUE
C
            DO 230 JR = 1, IRANK
               J = IRP - JR
               MJ = M - J + 1
               GAMMA = DDOT(MJ,Q(J,J),NRDA,X(J),1)/(DIAG(J)*Q(J,J))
               DO 220 K = J, M
                  X(K) = X(K) + GAMMA*Q(J,K)
  220          CONTINUE
  230       CONTINUE
C
C           RESCALE ANSWERS AS DICTATED
C
            DO 240 K = 1, M
               X(K) = X(K)*SCALES(K)
  240       CONTINUE
C
            IF (NU .EQ. 0 .OR. M .EQ. IRANK) GO TO 290
C
C              INITIALIZE U MATRIX AND THEN APPLY ORTHOGONAL
C              TRANSFORMATION
C
               L = M - IRANK
               DO 280 K = 1, L
                  DO 250 I = 1, M
                     U(I,K) = 0.0D0
                     IF (I .EQ. IRANK + K) U(I,K) = 1.0D0
  250             CONTINUE
C
                  DO 270 JR = 1, IRANK
                     J = IRP - JR
                     MJ = M - J + 1
                     GAMMA = DDOT(MJ,Q(J,J),NRDA,U(J,K),1)
     1                       /(DIAG(J)*Q(J,J))
                     DO 260 I = J, M
                        U(I,K) = U(I,K) + GAMMA*Q(J,I)
  260                CONTINUE
  270             CONTINUE
  280          CONTINUE
  290       CONTINUE
  300    CONTINUE
  310 CONTINUE
C
      RETURN
      END