home *** CD-ROM | disk | FTP | other *** search
/ The Fred Fish Collection 1.5 / ffcollection-1-5-1992-11.iso / ff_progs / libs / matlab.lzh / MATLAB / MATLAB.LZH / Source / MatLab / matfn4.for < prev    next >
Encoding:
Text File  |  1991-04-13  |  5.3 KB  |  195 lines
  1.       SUBROUTINE MATFN4
  2.       IMPLICIT NONE
  3. C
  4. C EVALUATE FUNCTIONS INVOLVING QR DECOMPOSITION (LEAST SQUARES)
  5. C
  6.       INCLUDE MATLAB$KOM:SIZEPARMS.INC
  7.       INCLUDE MATLAB$KOM:VSTK.KOM
  8.       INCLUDE MATLAB$KOM:IOP.KOM
  9.       INCLUDE MATLAB$KOM:COM.KOM
  10. C
  11.       DOUBLE PRECISION T, TOL, EPS
  12.       INTEGER QUOTE, J, JB, JOB, K, INFO, M, M2, MM, MN, N, N2
  13.       INTEGER L, L2, L3, L4, LE, LL, LS
  14. C
  15.       DOUBLE PRECISION FLOP, DFLOAT
  16. C
  17.       DATA QUOTE / 49 /
  18. C
  19. C
  20.       IF (DDT.EQ.1) WRITE (WTE, 100) FIN
  21. 100   FORMAT (' MATFN4', I4)
  22. C
  23.       L = LSTK(TOP)
  24.       M = MSTK(TOP)
  25.       N = NSTK(TOP)
  26.       IF (FIN.EQ.-1) GO TO 10
  27.       IF (FIN.EQ.-2) GO TO 20
  28.       GO TO 40
  29. C
  30. C ***      RECTANGULAR MATRIX RIGHT DIVISION, A/A2
  31. 10    CONTINUE
  32.       L2 = LSTK(TOP+1)
  33.       M2 = MSTK(TOP+1)
  34.       N2 = NSTK(TOP+1)
  35.       TOP = TOP+1
  36.       IF (N.GT.1 .AND. N.NE.N2) CALL ERROR (11)
  37.       IF (ERR.GT.0) RETURN
  38.       CALL STACK1 (QUOTE)
  39.       IF (ERR.GT.0) RETURN
  40.       LL = L2+M2*N2
  41.       CALL WCOPY (M*N, STKR(L), STKI(L), 1, STKR(LL), STKI(LL), 1)
  42.       CALL WCOPY (M*N+M2*N2, STKR(L2), STKI(L2), 1, STKR(L), STKI(L), 1)
  43.       LSTK(TOP) = L+M2*N2
  44.       MSTK(TOP) = M
  45.       NSTK(TOP) = N
  46.       CALL STACK1 (QUOTE)
  47.       IF (ERR.GT.0) RETURN
  48.       TOP = TOP-1
  49.       M = N2
  50.       N = M2
  51.       GO TO 20
  52. C
  53. C ***      RECTANGULAR MATRIX LEFT DIVISION A BACKSLASH A2
  54. 20    CONTINUE
  55.       L2 = LSTK(TOP+1)
  56.       M2 = MSTK(TOP+1)
  57.       N2 = NSTK(TOP+1)
  58.       IF (M2*N2.GT.1) GO TO 21
  59.       M2 = M
  60.       N2 = M
  61.       ERR = L2+M*M-LSTK(BOT)
  62.       IF (ERR.GT.0) CALL ERROR (17)
  63.       IF (ERR.GT.0) RETURN
  64.       CALL WSET (M*M-1, 0.0D0, 0.0D0, STKR(L2+1), STKI(L2+1), 1)
  65.       CALL WCOPY (M, STKR(L2), STKI(L2), 0, STKR(L2), STKI(L2), M+1)
  66. 21    CONTINUE
  67.       IF (M2.NE.M) CALL ERROR (12)
  68.       IF (ERR.GT.0) RETURN
  69.       L3 = L2+MAX0 (M, N)*N2
  70.       L4 = L3+N
  71.       ERR = L4+N-LSTK(BOT)
  72.       IF (ERR.GT.0) CALL ERROR (17)
  73.       IF (ERR.GT.0) RETURN
  74.       IF (M.GT.N) GO TO 23
  75.       DO 22 JB = 1, N2
  76.         J = N+1-JB
  77.         LS = L2+(J-1)*M
  78.         LL = L2+(J-1)*N
  79.         CALL WCOPY (M, STKR(LS), STKI(LS), -1, STKR(LL), STKI(LL), -1)
  80. 22    CONTINUE
  81. 23    CONTINUE
  82.       DO 24 J = 1, N
  83.         BUF(J) = 0
  84. 24    CONTINUE
  85.       CALL WQRDC (STKR(L), STKI(L), M, M, N, STKR(L4), STKI(L4),
  86.      .            BUF, STKR(L3), STKI(L3), 1)
  87.       K = 0
  88.       EPS = STKR(VSIZE-4)
  89.       T = DABS (STKR(L))+DABS (STKI(L))
  90.       TOL = FLOP (DFLOAT (MAX0 (M, N))*EPS*T)
  91.       MN = MIN0 (M, N)
  92.       DO 27 J = 1, MN
  93.         LS = L+J-1+(J-1)*M
  94.         T = DABS (STKR(LS))+DABS (STKI(LS))
  95.         IF (T.GT.TOL) K = J
  96. 27    CONTINUE
  97.       IF (K.LT.MN) WRITE (WTE, 28) K, TOL
  98.       IF (K.LT.MN .AND. WIO.NE.0) WRITE (WIO, 28) K, TOL
  99. 28    FORMAT (' RANK DEFICIENT,  RANK = ', I4, ',  TOL = ', 1PD13.4)
  100.       MN = MAX0 (M, N)
  101.       DO 29 J = 1, N2
  102.         LS = L2+(J-1)*MN
  103.         CALL WQRSL (STKR(L), STKI(L), M, M, K, STKR(L4), STKI(L4),
  104.      .              STKR(LS), STKI(LS), T, T, STKR(LS), STKI(LS),
  105.      .              STKR(LS), STKI(LS), T, T, T, T, 100, INFO)
  106.         LL = LS+K
  107.         CALL WSET (N-K, 0.0D0, 0.0D0, STKR(LL), STKI(LL), 1)
  108. 29    CONTINUE
  109.       DO 31 J = 1, N
  110.         BUF(J) = -BUF(J)
  111. 31    CONTINUE
  112.       DO 35 J = 1, N
  113.         IF (BUF(J).GT.0) GO TO 35
  114.         K = -BUF(J)
  115.         BUF(J) = K
  116.    33   CONTINUE
  117.         IF (K.EQ.J) GO TO 34
  118.         LS = L2+J-1
  119.         LL = L2+K-1
  120.         CALL WSWAP (N2, STKR(LS), STKI(LS), MN, STKR(LL), STKI(LL), MN)
  121.         BUF(K) = -BUF(K)
  122.         K = BUF(K)
  123.         GO TO 33
  124. C
  125. 34      CONTINUE
  126. 35    CONTINUE
  127.       DO 36 J = 1, N2
  128.         LS = L2+(J-1)*MN
  129.         LL = L+(J-1)*N
  130.         CALL WCOPY (N, STKR(LS), STKI(LS), 1, STKR(LL), STKI(LL), 1)
  131. 36    CONTINUE
  132.       MSTK(TOP) = N
  133.       NSTK(TOP) = N2
  134.       IF (FIN.EQ.-1) CALL STACK1 (QUOTE)
  135.       IF (ERR.GT.0) RETURN
  136.       GO TO 99
  137. C
  138. C ***      QR
  139. 40    CONTINUE
  140.       MM = MAX0 (M, N)
  141.       LS = L+MM*MM
  142.       IF (LHS.EQ.1 .AND. FIN.EQ.1) LS = L
  143.       LE = LS+M*N
  144.       L4 = LE+MM
  145.       ERR = L4+MM-LSTK(BOT)
  146.       IF (ERR.GT.0) CALL ERROR (17)
  147.       IF (ERR.GT.0) RETURN
  148.       IF (LS.NE.L) CALL WCOPY (M*N, STKR(L), STKI(L),
  149.      .                         1, STKR(LS), STKI(LS), 1)
  150.       JOB = 1
  151.       IF (LHS.LT.3) JOB = 0
  152.       DO 42 J = 1, N
  153.         BUF(J) = 0
  154. 42    CONTINUE
  155.       CALL WQRDC (STKR(LS), STKI(LS), M, M, N, STKR(L4),
  156.      .            STKI(L4), BUF, STKR(LE), STKI(LE), JOB)
  157.       IF (LHS.EQ.1 .AND. FIN.EQ.1) GO TO 99
  158.       CALL WSET (M*M, 0.0D0, 0.0D0, STKR(L), STKI(L), 1)
  159.       CALL WSET (M, 1.0D0, 0.0D0, STKR(L), STKI(L), M+1)
  160.       DO 43 J = 1, M
  161.         LL = L+(J-1)*M
  162.         CALL WQRSL (STKR(LS), STKI(LS), M, M, N, STKR(L4), STKI(L4),
  163.      .              STKR(LL), STKI(LL), STKR(LL), STKI(LL), T, T,
  164.      .              T, T, T, T, T, T, 10000, INFO)
  165. 43    CONTINUE
  166.       IF (FIN.EQ.2) GO TO 99
  167.       NSTK(TOP) = M
  168.       DO 45 J = 1, N
  169.         LL = LS+J+(J-1)*M
  170.         CALL WSET (M-J, 0.0D0, 0.0D0, STKR(LL), STKI(LL), 1)
  171. 45    CONTINUE
  172.       IF (TOP+1.GE.BOT) CALL ERROR (18)
  173.       IF (ERR.GT.0) RETURN
  174.       TOP = TOP+1
  175.       LSTK(TOP) = LS
  176.       MSTK(TOP) = M
  177.       NSTK(TOP) = N
  178.       IF (LHS.EQ.2) GO TO 99
  179.       CALL WSET (N*N, 0.0D0, 0.0D0, STKR(LE), STKI(LE), 1)
  180.       DO 47 J = 1, N
  181.         LL = LE+BUF(J)-1+(J-1)*N
  182.         STKR(LL) = 1.0D0
  183. 47    CONTINUE
  184.       IF (TOP+1.GE.BOT) CALL ERROR (18)
  185.       IF (ERR.GT.0) RETURN
  186.       TOP = TOP+1
  187.       LSTK(TOP) = LE
  188.       MSTK(TOP) = N
  189.       NSTK(TOP) = N
  190.       GO TO 99
  191. C
  192. 99    CONTINUE
  193.       RETURN
  194.       END
  195.