home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l292 / 1.ddi / MATMATH.FOR < prev    next >
Encoding:
Text File  |  1990-02-20  |  9.7 KB  |  384 lines
  1.  
  2.       SUBROUTINE SwapMatVals(x,y)
  3.       REAL x,y,temp
  4.         temp = x
  5.         x = y
  6.         y = temp
  7.       END !SUB swap
  8.  
  9.       SUBROUTINE CopyMatrixVectors (xin, xout, n, r)
  10.       INCLUDE 'STDHDR.FOR'
  11.  
  12.       REAL xin(0:maxr,0:maxc),xout(0:maxc)
  13.       INTEGER n,r,i
  14.  
  15.       DO i = 0, n - 1
  16.          xout(i) = xin(r, i)
  17.       END DO
  18.  
  19.       END !SUB
  20.  
  21.  
  22.  
  23.  
  24.       SUBROUTINE Copyto2DVector (xin, xout, n, r)
  25.       INCLUDE 'STDHDR.FOR'
  26.       REAL xout(0:maxr,0:maxc),xin(0:maxc)
  27.       INTEGER n,r,i
  28.  
  29.       DO i = 0, n - 1
  30.          xout(r, i) = xin(i)
  31.       END DO
  32.  
  33.       END !SUB
  34.  
  35.  
  36.  
  37.       SUBROUTINE CyclicJacobi (a, n, eigenvalues, v, count, success)
  38.       INCLUDE 'STDHDR.FOR'
  39.       REAL a(0:maxr,0:maxc), eigenvalues(0:maxc)
  40.       REAL v(0:maxr,0:maxc), tolerance, nearzero
  41.       REAL atemp1(0:maxc), atemp2(0: maxc), sumsq
  42.       INTEGER Maxcount, p, q
  43.       INTEGER n,count
  44.       LOGICAL success
  45.  
  46.       PARAMETER (nearzero = 1E-08)   ! {change to 1.0e-15 FOR 8087 system}
  47.       PARAMETER (tolerance = 1E-08)
  48.  
  49.       Maxcount = count
  50.       count = 0
  51.       DO p = 0, n - 1
  52.          DO q = 0, n - 1
  53.             v(p, q) = 0.0
  54.          END DO
  55.          v(p, p) = 1.0
  56.       END DO
  57.  
  58.       DO WHILE (.NOT. success .OR. count .LE. Maxcount)
  59.         count = count + 1
  60.         sumsq = 0.0
  61.         DO p = 0, n - 1
  62.            sumsq = sumsq + a(p, p) ** 2
  63.         END DO
  64.         success = .TRUE.
  65.         DO p = 0, n - 2
  66.            DO q = p + 1, n - 1
  67.               IF (ABS(a(p, q) / sumsq) .GT. tolerance ) THEN
  68.                 success = .FALSE.
  69.                 CALL CopyMatrixVectors(a, atemp1, n, p)
  70.                 CALL CopyMatrixVectors(a, atemp2, n, q)
  71.                 CALL JacobiRotation(v, a, p, q, atemp1, atemp2, n)
  72.                 CALL Copyto2DVector(atemp1, a, n, p)
  73.                 CALL Copyto2DVector(atemp2, a, n, q)
  74.               END IF
  75.           END DO
  76.         END DO
  77.       END DO
  78.       IF (success) THEN
  79.         DO p = 0, n - 1
  80.         eigenvalues(p) = a(p, p)
  81.         END DO
  82.       END IF
  83.       END !SUB cycjacob
  84.  
  85.  
  86.  
  87.       SUBROUTINE JacobiRotation (v, m, p, q, mp, mq, n)
  88.       INCLUDE 'STDHDR.FOR'
  89.       REAL v(0:maxr,0:maxc),m(0:maxr,0:maxc)
  90.       REAL mp(0:maxc), mq(0:maxc), tolerance, nearzero
  91.       REAL g, h, mpp, mqq, mpq, mpj, mqj, tau
  92.       REAL tempr, angle, c, s, cs, cSquared, sSquared
  93.       REAL PiOver4, PiOver2
  94.       INTEGER n,p,q, j
  95.  
  96.       PARAMETER (nearzero = 1E-08)   ! {change to 1.0e-15 FOR 8087 system}
  97.       PARAMETER (tolerance = 1E-08)
  98.  
  99.       PiOver4 = pi / 4.0
  100.       PiOver2 = pi / 2.0
  101.  
  102.       mpp = mp(p)
  103.       mpq = mp(q)
  104.       mqq = mq(q)
  105.       IF (ABS(mpp - mqq) .GT. nearzero) THEN
  106.          angle = ATAN(2.0 * mpq / (mpp - mqq)) / 2.0
  107.         IF (angle .GT. PiOver4) angle = angle - PiOver2
  108.       ELSE
  109.        IF (mpq .GT. 0 ) THEN
  110.            angle = PiOver4
  111.         ELSE
  112.           angle = -PiOver4
  113.         END IF
  114.       END IF
  115.       c = COS(angle)
  116.       s = SIN(angle)
  117.       cSquared = c ** 2
  118.       sSquared = s ** 2
  119.       cs = c * s
  120.       tempr = 2.0 * mpq * cs
  121.       mp(p) = mpp * cSquared + tempr + mqq * sSquared
  122.       mq(q) = mpp * sSquared - tempr + mqq * cSquared
  123.       mp(q) = 0.0
  124.       mq(p) = 0.0
  125.       tau = s / (1.0 + c)
  126.       DO j = 0, n - 1
  127.          IF ( j .NE. p .AND. j .NE. q) THEN
  128.              mpj = mp(j) * c + mq(j) * s
  129.              mqj = mq(j) * c - mp(j) * s
  130.              mp(j) = mpj
  131.              m(j, p) = mpj
  132.              mq(j) = mqj
  133.              m(j, q) = mqj
  134.          END IF
  135.       END DO
  136.       DO j = 0, n - 1
  137.          g = v(j, p)
  138.          h = v(j, q)
  139.          v(j, p) = c * g + s * h
  140.          v(j, q) = -s * g + c * h
  141.       END DO
  142.  
  143.       END !SUB JacobRot
  144.  
  145.  
  146.  
  147.       SUBROUTINE MatAdd (m1, m2, numrow, numcol, m3)
  148.       INCLUDE 'STDHDR.FOR'
  149.       REAL m1(0:maxr,0:maxc), m2(0:maxr,0:maxc), m3(0:maxr,0:maxc)
  150.       INTEGER numrow,numcol,i,j
  151.  
  152.  
  153.       DO i = 0, numrow - 1
  154.          DO j = 0, numcol - 1
  155.             m3(i, j) = m1(i, j) + m2(i, j)
  156.          END DO
  157.       END DO
  158.       END !SUB MatAdd
  159.  
  160.  
  161.  
  162.       SUBROUTINE MatDeter (m1, numrow, det)
  163.       INCLUDE 'STDHDR.FOR'
  164.       REAL m1(0:maxr,0:maxc)
  165.       INTEGER i,j,k,f, lx, nxt,iErr
  166.  
  167.       REAL bignum, eval, temp1, temp2, rac
  168.       REAL tempmat[ALLOCATABLE](:,:)
  169.       ALLOCATE(tempmat(0:maxc,0:maxc),STAT=iErr)
  170.       CALL CheckMem(iErr)
  171.  
  172.       DO i = 0, numrow - 1
  173.          DO j = 0, numrow - 1
  174.             tempmat(i, j) = m1(i, j)
  175.          END DO
  176.       END DO
  177.  
  178.       f = 1
  179.       DO i = 0, numrow - 1
  180.         bignum = 0
  181.         DO k = i, numrow - 1
  182.           eval = ABS(tempmat(k, i))
  183.           IF (eval - bignum .GT. 0) THEN
  184.             bignum = eval
  185.             lx = k
  186.           END IF
  187.         END DO
  188.         IF (i - lx .NE. 0 ) f = -f
  189.         DO j = 0, numrow - 1
  190.           temp1 = tempmat(i, j)
  191.           tempmat(i, j) = tempmat(lx, j)
  192.           tempmat(lx, j) = temp1
  193.         END DO
  194.         pivot = tempmat(i, i)
  195.         nxt = i + 1
  196.         DO j = nxt, numrow - 1
  197.           rac = tempmat(j, i) / pivot
  198.           DO k = i, numrow - 1
  199.             tempmat(j, k) = tempmat(j, k) - rac * tempmat(i, k)
  200.           END DO
  201.         END DO
  202.       END DO
  203.       temp2 = 1
  204.       DO i = 0, numrow - 1
  205.          temp2 = temp2 * tempmat(i, i)
  206.       END DO
  207.       det = temp2 * f
  208.       DEALLOCATE(tempmat,STAT=iErr)
  209.       CALL CheckDealloc(iErr)
  210.       END !SUB       matdet
  211.  
  212.  
  213.  
  214.       SUBROUTINE MatInvert (matdata, numcol, det, invary)
  215.       INCLUDE 'STDHDR.FOR'
  216.       REAL matdata(0:maxr,0:maxc), invary(0:maxr,0:maxc)
  217.       REAL det,piv, t, leval
  218.       INTEGER numcol, i, j, k, l, lerow, lecol, l1, iErr
  219.       INTEGER pivlst[ALLOCATABLE](:,:)
  220.       LOGICAL pivchk[ALLOCATABLE](:)
  221.  
  222.       ALLOCATE(pivlst(0:maxc,0:1),STAT=iErr)
  223.       CALL CheckMem(iErr)
  224.       ALLOCATE(pivchk(0:maxc),STAT=iErr)
  225.       CALL CheckMem(iErr)
  226.  
  227.       det = 1
  228.       DO i = 0, numcol - 1
  229.          pivchk(i) = .FALSE.
  230.          DO j = 0, numcol - 1
  231.            invary(i, j) = matdata(i, j)
  232.          END DO
  233.       END DO
  234.  
  235.       DO i = 0, numcol - 1
  236.         leval = 0.0
  237.         DO j = 0, numcol - 1
  238.            IF (.NOT. pivchk(j)) THEN
  239.              DO k = 0, numcol - 1
  240.                IF (.NOT. pivchk(k)) THEN
  241.                  IF ( ABS(invary(j, k)) .GT. leval ) THEN
  242.                    lerow = j
  243.                    lecol = k
  244.                    leval = ABS(invary(j, k))
  245.                  END IF
  246.                END IF
  247.              END DO
  248.            END IF
  249.         END DO
  250.         pivchk(lecol) = .TRUE.
  251.         pivlst(i, 0) = lerow
  252.         pivlst(i, 1) = lecol
  253.         IF (lerow .NE. lecol) THEN
  254.           det = -det
  255.           DO l = 0, numcol - 1
  256.             call SwapMatVals( invary(lerow, l), invary(lecol, l))
  257.           END DO
  258.         END IF
  259.  
  260.         piv = invary(lecol, lecol)
  261.         det = det * piv
  262.         IF (det .GT. 1E30 )THEN
  263.           det = 1      ! {prevents numeric overflow}
  264.         END IF
  265.         invary(lecol, lecol) = 1.0
  266.         DO l = 0, numcol - 1
  267.           invary(lecol, l) = invary(lecol, l) / piv
  268.         END DO
  269.         DO l1 = 0, numcol - 1
  270.           IF (l1 .NE. lecol) THEN
  271.             t = invary(l1, lecol)
  272.             invary(l1, lecol) = 0
  273.             DO l = 0, numcol - 1
  274.               invary(l1, l) = invary(l1, l) - invary(lecol, l) * t
  275.             END DO
  276.           END IF
  277.         END DO
  278.       END DO
  279.  
  280.       DO i = 0, numcol - 1
  281.         l = numcol - i - 1
  282.         IF (pivlst(l, 0) .NE. pivlst(l, 1)) THEN
  283.           lerow = pivlst(l, 0)
  284.           lecol = pivlst(l, 1)
  285.           DO k = 0, numcol-1
  286.             call SWAPMatVals( invary(k, lerow), invary(k, lecol))
  287.           END DO
  288.        END IF
  289.       END DO
  290.       DEALLOCATE(pivlst,STAT=iErr)
  291.       CALL CheckDealloc(iErr)
  292.       DEALLOCATE(pivchk,STAT=iErr)
  293.       CALL CheckDealloc(iErr)
  294.       END !SUB MatInvert
  295.  
  296.  
  297.  
  298.       SUBROUTINE MatPrint (m1, m, n, matname)
  299.       INCLUDE 'STDHDR.FOR'
  300.       REAL m1(0:maxr,0:maxc)
  301.       INTEGER m,n
  302.       INTEGER i,j
  303.       CHARACTER* (*) matname
  304.  
  305.       WRITE (*, 10)  matname
  306. 10    FORMAT (1X, A80)
  307.       DO i = 0, m - 1
  308.          DO j = 0, n - 1
  309.            WRITE (*,20) m1(i, j)
  310. 20         FORMAT ( 5x,  F9.3 \)
  311.          END DO
  312.          WRITE (*,*)
  313.       END DO
  314.       END !SUB matPrint
  315.  
  316.  
  317.  
  318.       SUBROUTINE MatProd (m1, m2, l, m, n, m3)
  319.       INCLUDE 'STDHDR.FOR'
  320.       REAL m1(0:maxr,0:maxc), m2(0:maxr,0:maxc), m3(0:maxr,0:maxc)
  321.       INTEGER l,m,n,i,j,k
  322.  
  323.       DO i = 0, l - 1
  324.          DO j = 0, n - 1
  325.             m3(i, j) = 0
  326.             DO k = 0, m - 1
  327.                m3(i, j) = m3(i, j) + m1(i, k) * m2(k, j)
  328.             END DO
  329.          END DO
  330.       END DO
  331.       END !SUB MatProd
  332.  
  333.  
  334.  
  335.       SUBROUTINE MatScalarProd (m1, sval, numrow, numcol, m2)
  336.       INCLUDE 'STDHDR.FOR'
  337.       REAL m1,m2,sval
  338.       DIMENSION m1(0:maxr,0:maxc), m2(0:maxr,0:maxc)
  339.       INTEGER numrow,numcol
  340.       INTEGER i,j
  341.  
  342.       DO i = 0, numrow - 1
  343.          DO j = 0, numcol - 1
  344.             m2(i, j) = sval * m1(i, j)
  345.          END DO
  346.       END DO
  347.       END !SUB  MatScalarProduct
  348.  
  349.  
  350.  
  351.       SUBROUTINE MatTranspose (m1, numrow, numcol, m2)
  352.       INCLUDE 'STDHDR.FOR'
  353.       REAL m1,m2
  354.       DIMENSION m1(0:maxr,0:maxc), m2(0:maxr,0:maxc)
  355.       INTEGER numrow,numcol
  356.       INTEGER i,j
  357.  
  358.       DO i = 0, numrow - 1
  359.          DO j = 0, numcol - 1
  360.             m2(j, i) = m1(i, j)
  361.          END DO
  362.       END DO
  363.       END !SUB MatTranspose
  364.  
  365.  
  366.  
  367.       SUBROUTINE NormalizeEigenVectors (a, n)
  368.       INCLUDE 'STDHDR.FOR'
  369.       REAL a, max
  370.       DIMENSION a(0:maxr,0:maxc)
  371.       INTEGER n, i, j
  372.  
  373.       DO j = 0, n - 1
  374.         max = a(0, j)
  375.         DO i = 1, n - 1
  376.            IF (ABS(a(i, j)) .GT. ABS(max)) max = a(i, j)
  377.         END DO
  378.         DO i = 0, n - 1
  379.           a(i, j) = a(i, j) / max
  380.         END DO
  381.       END DO
  382.       END !SUB NormEigVec
  383.  
  384.