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 / matfn3.for < prev    next >
Encoding:
Text File  |  1991-04-13  |  7.0 KB  |  279 lines
  1.       SUBROUTINE MATFN3
  2.       IMPLICIT NONE
  3. C
  4. C EVALUATE FUNCTIONS INVOLVING SINGULAR VALUE DECOMPOSITION
  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.       INTEGER I, J, JB, JOB, K, M, MN, N
  12.       INTEGER L, L1, L2, LD, LS, LI, LJ, LU, LV, LL
  13.       LOGICAL FRO, INF
  14.       DOUBLE PRECISION P, S, T, TOL, EPS
  15. C
  16.       INTEGER IWAMAX
  17.       DOUBLE PRECISION WDOTCR, WDOTCI, PYTHAG, WNRM2, WASUM
  18.       DOUBLE PRECISION FLOP, DFLOAT
  19. C
  20. C
  21.       IF (DDT.EQ.1) WRITE (WTE, 100) FIN
  22. 100   FORMAT (' MATFN3', I4)
  23. C
  24.       IF (FIN.EQ.1 .AND. RHS.EQ.2) TOP = TOP-1
  25.       L = LSTK(TOP)
  26.       M = MSTK(TOP)
  27.       N = NSTK(TOP)
  28.       MN = M*N
  29.       GO TO (50, 70, 10, 30, 70), FIN
  30. C
  31. C ***      COND
  32. 10    CONTINUE
  33.       LD = L+M*N
  34.       L1 = LD+MIN0 (M+1, N)
  35.       L2 = L1+N
  36.       ERR = L2+MIN0 (M, N)-LSTK(BOT)
  37.       IF (ERR.GT.0) CALL ERROR (17)
  38.       IF (ERR.GT.0) RETURN
  39.       CALL WSVDC (STKR(L), STKI(L), M, M, N, STKR(LD), STKI(LD),
  40.      .            STKR(L1), STKI(L1), T, T, 1, T, T, 1, STKR(L2),
  41.      .            STKI(L2), 0, ERR)
  42.       IF (ERR.NE.0) CALL ERROR (24)
  43.       IF (ERR.GT.0) RETURN
  44.       S = STKR(LD)
  45.       LD = LD+MIN0 (M, N)-1
  46.       T = STKR(LD)
  47.       IF (T.EQ.0.0D0) GO TO 13
  48.       STKR(L) = FLOP (S/T)
  49.       STKI(L) = 0.0D0
  50.       MSTK(TOP) = 1
  51.       NSTK(TOP) = 1
  52.       GO TO 99
  53. C
  54. 13    CONTINUE
  55.       WRITE (WTE, 14)
  56.       IF (WIO.NE.0) WRITE (WIO, 14)
  57. 14    FORMAT (' CONDITION IS INFINITE')
  58.       MSTK(TOP) = 0
  59.       GO TO 99
  60. C
  61. C ***      NORM
  62. 30    CONTINUE
  63.       P = 2.0D0
  64.       INF = .FALSE.
  65.       IF (RHS.NE.2) GO TO 31
  66.       FRO = IDINT (STKR(L)).EQ.15 .AND. MN.GT.1
  67.       INF = IDINT (STKR(L)).EQ.18 .AND. MN.GT.1
  68.       IF (.NOT.FRO) P = STKR(L)
  69.       TOP = TOP-1
  70.       L = LSTK(TOP)
  71.       M = MSTK(TOP)
  72.       N = NSTK(TOP)
  73.       MN = M*N
  74.       IF (FRO) M = MN
  75.       IF (FRO) N = 1
  76. 31    CONTINUE
  77.       IF (M.GT.1 .AND. N.GT.1) GO TO 40
  78.       IF (P.EQ.1.0D0) GO TO 36
  79.       IF (P.EQ.2.0D0) GO TO 38
  80.       I = IWAMAX (MN, STKR(L), STKI(L), 1)+L-1
  81.       S = DABS (STKR(I))+DABS (STKI(I))
  82.       IF (INF .OR. S.EQ.0.0D0) GO TO 49
  83.       T = 0.0D0
  84.       DO 33 I = 1, MN
  85.         LS = L+I-1
  86.         T = FLOP (T+(PYTHAG (STKR(LS), STKI(LS))/S)**P)
  87. 33    CONTINUE
  88.       IF (P.NE.0.0D0) P = 1.0D0/P
  89.       S = FLOP (S*T**P)
  90.       GO TO 49
  91. C
  92. 36    CONTINUE
  93.       S = WASUM (MN, STKR(L), STKI(L), 1)
  94.       GO TO 49
  95. C
  96. 38    CONTINUE
  97.       S = WNRM2 (MN, STKR(L), STKI(L), 1)
  98.       GO TO 49
  99. C
  100. C ***      MATRIX NORM
  101. 40    CONTINUE
  102.       IF (INF) GO TO 43
  103.       IF (P.EQ.1.0D0) GO TO 46
  104.       IF (P.NE.2.0D0) CALL ERROR (23)
  105.       IF (ERR.GT.0) RETURN
  106.       LD = L+M*N
  107.       L1 = LD+MIN0 (M+1, N)
  108.       L2 = L1+N
  109.       ERR = L2+MIN0 (M, N)-LSTK(BOT)
  110.       IF (ERR.GT.0) CALL ERROR (17)
  111.       IF (ERR.GT.0) RETURN
  112.       CALL WSVDC (STKR(L), STKI(L), M, M, N, STKR(LD), STKI(LD),
  113.      .            STKR(L1), STKI(L1), T, T, 1, T, T, 1, STKR(L2),
  114.      .            STKI(L2), 0, ERR)
  115.       IF (ERR.NE.0) CALL ERROR (24)
  116.       IF (ERR.GT.0) RETURN
  117.       S = STKR(LD)
  118.       GO TO 49
  119. C
  120. 43    CONTINUE
  121.       S = 0.0D0
  122.       DO 45 I = 1, M
  123.         LI = L+I-1
  124.         T = WASUM (N, STKR(LI), STKI(LI), M)
  125.         S = DMAX1 (S, T)
  126. 45    CONTINUE
  127.       GO TO 49
  128. C
  129. 46    CONTINUE
  130.       S = 0.0D0
  131.       DO 48 J = 1, N
  132.         LJ = L+(J-1)*M
  133.         T = WASUM (M, STKR(LJ), STKI(LJ), 1)
  134.         S = DMAX1 (S, T)
  135. 48    CONTINUE
  136.       GO TO 49
  137. C
  138. 49    CONTINUE
  139.       STKR(L) = S
  140.       STKI(L) = 0.0D0
  141.       MSTK(TOP) = 1
  142.       NSTK(TOP) = 1
  143.       GO TO 99
  144. C
  145. C ***      SVD
  146. 50    CONTINUE
  147.       IF (LHS.NE.3) GO TO 52
  148.       K = M
  149.       IF (RHS.EQ.2) K = MIN0 (M, N)
  150.       LU = L+M*N
  151.       LD = LU+M*K
  152.       LV = LD+K*N
  153.       L1 = LV+N*N
  154.       L2 = L1+N
  155.       ERR = L2+MIN0 (M, N)-LSTK(BOT)
  156.       IF (ERR.GT.0) CALL ERROR (17)
  157.       IF (ERR.GT.0) RETURN
  158.       JOB = 11
  159.       IF (RHS.EQ.2) JOB = 21
  160.       CALL WSVDC (STKR(L), STKI(L), M, M, N, STKR(LD), STKI(LD),
  161.      .            STKR(L1), STKI(L1), STKR(LU), STKI(LU), M,
  162.      .            STKR(LV), STKI(LV), N, STKR(L2), STKI(L2), JOB, ERR)
  163.       DO 53 JB = 1, N
  164.         DO 51 I = 1, K
  165.           J = N+1-JB
  166.           LL = LD+I-1+(J-1)*K
  167.           IF (I.NE.J) STKR(LL) = 0.0D0
  168.           STKI(LL) = 0.0D0
  169.           LS = LD+I-1
  170.           IF (I.EQ.J) STKR(LL) = STKR(LS)
  171.           LS = L1+I-1
  172.           IF (ERR.NE.0 .AND. I.EQ.J-1) STKR(LL) = STKR(LS)
  173. 51      CONTINUE
  174. 53    CONTINUE
  175.       IF (ERR.NE.0) CALL ERROR (24)
  176.       ERR = 0
  177.       CALL WCOPY (M*K+K*N+N*N, STKR(LU), STKI(LU), 1,
  178.      .            STKR(L), STKI(L), 1)
  179.       MSTK(TOP) = M
  180.       NSTK(TOP) = K
  181.       IF (TOP+1.GE.BOT) CALL ERROR (18)
  182.       IF (ERR.GT.0) RETURN
  183.       TOP = TOP+1
  184.       LSTK(TOP) = L+M*K
  185.       MSTK(TOP) = K
  186.       NSTK(TOP) = N
  187.       IF (TOP+1.GE.BOT) CALL ERROR (18)
  188.       IF (ERR.GT.0) RETURN
  189.       TOP = TOP+1
  190.       LSTK(TOP) = L+M*K+K*N
  191.       MSTK(TOP) = N
  192.       NSTK(TOP) = N
  193.       GO TO 99
  194. C
  195. 52    CONTINUE
  196.       LD = L+M*N
  197.       L1 = LD+MIN0 (M+1, N)
  198.       L2 = L1+N
  199.       ERR = L2+MIN0 (M, N)-LSTK(BOT)
  200.       IF (ERR.GT.0) CALL ERROR (17)
  201.       IF (ERR.GT.0) RETURN
  202.       CALL WSVDC (STKR(L), STKI(L), M, M, N, STKR(LD), STKI(LD),
  203.      .            STKR(L1), STKI(L1), T, T, 1, T, T, 1, STKR(L2),
  204.      .            STKI(L2), 0, ERR)
  205.       IF (ERR.NE.0) CALL ERROR (24)
  206.       IF (ERR.GT.0) RETURN
  207.       K = MIN0 (M, N)
  208.       CALL WCOPY (K, STKR(LD), STKI(LD), 1, STKR(L), STKI(L), 1)
  209.       MSTK(TOP) = K
  210.       NSTK(TOP) = 1
  211.       GO TO 99
  212. C
  213. C ***      PINV AND RANK
  214. 70    CONTINUE
  215.       TOL = -1.0D0
  216.       IF (RHS.NE.2) GO TO 71
  217.       TOL = STKR(L)
  218.       TOP = TOP-1
  219.       L = LSTK(TOP)
  220.       M = MSTK(TOP)
  221.       N = NSTK(TOP)
  222. 71    CONTINUE
  223.       LU = L+M*N
  224.       LD = LU+M*M
  225.       IF (FIN.EQ.5) LD = L+M*N
  226.       LV = LD+M*N
  227.       L1 = LV+N*N
  228.       IF (FIN.EQ.5) L1 = LD+N
  229.       L2 = L1+N
  230.       ERR = L2+MIN0 (M, N)-LSTK(BOT)
  231.       IF (ERR.GT.0) CALL ERROR (17)
  232.       IF (ERR.GT.0) RETURN
  233.       IF (FIN.EQ.2) JOB = 11
  234.       IF (FIN.EQ.5) JOB = 0
  235.       CALL WSVDC (STKR(L), STKI(L), M, M, N, STKR(LD), STKI(LD),
  236.      .            STKR(L1), STKI(L1), STKR(LU), STKI(LU), M,
  237.      .            STKR(LV), STKI(LV), N, STKR(L2), STKI(L2), JOB, ERR)
  238.       IF (ERR.NE.0) CALL ERROR (24)
  239.       IF (ERR.GT.0) RETURN
  240.       EPS = STKR(VSIZE-4)
  241.       IF (TOL.LT.0.0D0) TOL = FLOP (DFLOAT (MAX0 (M, N))*EPS*STKR(LD))
  242.       MN = MIN0 (M, N)
  243.       K = 0
  244.       DO 72 J = 1, MN
  245.         LS = LD+J-1
  246.         S = STKR(LS)
  247.         IF (S.LE.TOL) GO TO 73
  248.         K = J
  249.         LL = LV+(J-1)*N
  250.         IF (FIN.EQ.2) CALL WRSCAL (N, 1.0D0/S, STKR(LL), STKI(LL), 1)
  251. 72    CONTINUE
  252. 73    CONTINUE
  253.       IF (FIN.EQ.5) GO TO 78
  254.       DO 77 J = 1, M
  255.         DO 76 I = 1, N
  256.           LL = L+I-1+(J-1)*N
  257.           L1 = LV+I-1
  258.           L2 = LU+J-1
  259.           STKR(LL) = WDOTCR (K, STKR(L2), STKI(L2),
  260.      .                       M, STKR(L1), STKI(L1), N)
  261.           STKI(LL) = WDOTCI (K, STKR(L2), STKI(L2),
  262.      .                       M, STKR(L1), STKI(L1), N)
  263. 76      CONTINUE
  264. 77    CONTINUE
  265.       MSTK(TOP) = N
  266.       NSTK(TOP) = M
  267.       GO TO 99
  268. C
  269. 78    CONTINUE
  270.       STKR(L) = DFLOAT (K)
  271.       STKI(L) = 0.0D0
  272.       MSTK(TOP) = 1
  273.       NSTK(TOP) = 1
  274.       GO TO 99
  275. C
  276. 99    CONTINUE
  277.       RETURN
  278.       END
  279.