home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l292 / 1.ddi / COMPLEX.FOR < prev    next >
Encoding:
Text File  |  1990-02-20  |  7.8 KB  |  284 lines
  1.       FUNCTION CAngle (CNum1)
  2.       INCLUDE 'STDHDR.FOR'
  3.       COMPLEX CNum1
  4.       REAL NearZero
  5.       PARAMETER (NearZero = 1E-20)
  6.  
  7.  
  8.       IF (ABS(REAL(CNum1)) .LT. NearZero) THEN
  9.         IF (ABS(IMAG(CNum1)) .LT. NearZero)  CAngle = 0
  10.         IF (IMAG(CNum1) .LT. -(NearZero))  CAngle = 3.0 * pi / 2.0
  11.         IF (IMAG(CNum1) .GT. NearZero) CAngle = pi / 2.0
  12.       ELSE
  13.         IF (ABS(IMAG(CNum1)) .LT. NearZero) THEN
  14.           IF (REAL(CNum1) .LT. -(NearZero))  CAngle = pi
  15.           IF (REAL(CNum1) .GT. NearZero)  CAngle = 0
  16.         ELSE
  17.           IF (REAL(CNum1) .GT. NearZero .AND.
  18.      +        IMAG(CNum1) .GT. NearZero)
  19.      +           CAngle = ATAN(IMAG(CNum1) / REAL(CNum1))
  20.           IF (REAL(CNum1) .LT. -NearZero .AND.
  21.      +        IMAG(CNum1) .GT. NearZero)
  22.      +           CAngle = pi - ATAN(-IMAG(CNum1) / REAL(CNum1))
  23.           IF (REAL(CNum1) .LT. -NearZero .AND.
  24.      +        IMAG(CNum1) .LT. -NearZero)
  25.      +           CAngle = pi + ATAN(IMAG(CNum1) / REAL(CNum1))
  26.           IF (REAL(CNum1) .GT. NearZero .AND.
  27.      +        IMAG(CNum1) .LT. -NearZero)
  28.      +         CAngle = 2.0 * pi - ATAN(-IMAG(CNum1) / REAL(CNum1))
  29.         END IF
  30.       END IF
  31.       END !FUNCTION CAngle
  32.  
  33.       SUBROUTINE CSwap (n1, n2)
  34.       REAL n1, n2 ,temp
  35.         temp = n1
  36.         n1 =  n2
  37.         n2 = temp
  38.       END !SUB CSwap
  39.  
  40.  
  41.       FUNCTION CMag (CNum1)
  42.       COMPLEX CNum1
  43.         CMag = SQRT((REAL(CNum1) ** 2) + (IMAG(CNum1) ** 2))
  44.       END !FUNCTION CMag
  45.  
  46.  
  47.  
  48.       SUBROUTINE CMatAdd (m1, m2, numrow, numcol, m3)
  49.       INCLUDE 'STDHDR.FOR'
  50.       COMPLEX m1(0:maxr,0:maxc), m2(0:maxr, 0:maxc), m3(0:maxr, 0:maxc)
  51.       INTEGER numrow,numcol
  52.       INTEGER i, j
  53.  
  54.       DO i = 0, numrow - 1
  55.         DO j = 0, numcol - 1
  56.           CALL ComplexMath(m1(i, j), 0, m2(i, j), m3(i, j))
  57.         END DO
  58.       END DO
  59.       END !SUB CMatAdd
  60.  
  61.  
  62.  
  63.       SUBROUTINE ComplexMatInvert (matdata, numcol, det, invary)
  64.       INCLUDE 'STDHDR.FOR'
  65.       REAL matdata(0:maxr,0:maxc),  invary(0:maxr, 0:maxc)
  66.       INTEGER numcol
  67.       REAL det
  68.  
  69.       INTEGER i, j, k, l, lerow, lecol, l1, iErr
  70.       REAL piv, t, leval
  71.       INTEGER pivlst[ALLOCATABLE](:,:)
  72.       LOGICAL pivchk[ALLOCATABLE](:)
  73.       ALLOCATE(pivlst(0:maxc,0:1),STAT=iErr)
  74.       CALL CheckMem(iErr)
  75.       ALLOCATE(pivchk(0:maxc),STAT=iErr)
  76.       CALL CheckMem(iErr)
  77.  
  78.       det = 1.0
  79.  
  80.       DO i = 0, numcol - 1
  81.         pivchk(i) = .false.
  82.         DO j = 0, numcol - 1
  83.           invary(i, j) = matdata(i, j)
  84.         END DO
  85.       END DO
  86.  
  87.       DO i = 0, numcol - 1
  88.         leval = 0.0
  89.         DO j = 0, numcol - 1
  90.            IF (.NOT. pivchk(j)) THEN
  91.              DO k = 0, numcol - 1
  92.                IF (.NOT. pivchk(k)) THEN
  93.                  IF (ABS(invary(j, k)) .GT. leval) THEN
  94.                    lerow = j
  95.                    lecol = k
  96.                    leval = ABS(invary(j, k))
  97.                  END IF
  98.                END IF
  99.              END DO
  100.            END IF
  101.          END DO
  102.  
  103.          pivchk(lecol) = .TRUE.
  104.          pivlst(i, 0) = lerow
  105.          pivlst(i, 1) = lecol
  106.  
  107.         IF (lerow .NE. lecol) THEN
  108.            det = -det
  109.            DO l = 0, numcol - 1
  110.              CALL CSwap (invary(lerow, l), invary(lecol, l))
  111.            END DO
  112.          END IF
  113.  
  114.          piv = invary(lecol, lecol)
  115.          det = det * piv
  116.          IF (det .GT. 1E30) THEN
  117.            det = 1E30       ! {prevents numeric overflow}
  118.          END IF
  119.          invary(lecol, lecol) = 1.0
  120.          DO l = 0, numcol - 1
  121.            invary(lecol, l) = invary(lecol, l) / piv
  122.          END DO
  123.          DO l1 = 0, numcol - 1
  124.            IF (l1 .NE. lecol)  THEN
  125.              t = invary(l1, lecol)
  126.              invary(l1, lecol) = 0
  127.              DO l = 0, numcol - 1
  128.                invary(l1, l) = invary(l1, l) - invary(lecol, l) * t
  129.              END DO
  130.            END IF
  131.          END DO
  132.        END DO
  133.  
  134.        DO i = 0, numcol - 1
  135.          l = numcol - i - 1
  136.          IF (pivlst(l, 0) .NE. pivlst(l, 1)) THEN
  137.            lerow = pivlst(l, 0)
  138.            lecol = pivlst(l, 1)
  139.            DO k = 0, numcol - 1
  140.              CALL CSwap( invary(k, lerow), invary(k, lecol))
  141.            END DO
  142.          END IF
  143.        END DO
  144.        DEALLOCATE(pivlst,STAT=iErr)
  145.        CALL CheckDealloc(iErr)
  146.        DEALLOCATE(pivchk,STAT=iErr)
  147.        CALL CheckDealloc(iErr)
  148.       END !SUB
  149.  
  150.  
  151.  
  152.  
  153.  
  154.       SUBROUTINE CMatInvert (matdata, numcol, det, invary)
  155.       INCLUDE 'STDHDR.FOR'
  156.       COMPLEX matdata(0:maxr,0:maxc), invary(0:maxr, 0:maxc)
  157.       REAL det
  158.       REAL coef[ALLOCATABLE](:,:), inv[ALLOCATABLE](:,:)
  159.       INTEGER numcol, i,j,nn,ii,jj,iErr
  160.       ALLOCATE (coef(0:maxr, 0:maxc),inv(0:maxr, 0:maxc), STAT=iErr)
  161.  
  162.       DO i = 0, numcol - 1
  163.         ii = 2*i
  164.         DO j = 0, numcol - 1
  165.           jj = 2*j
  166.           coef(ii, jj) = REAL(matdata(i, j))
  167.           coef(ii, jj+ 1) = -IMAG(matdata(i, j))
  168.           coef(ii + 1, jj ) = IMAG(matdata(i, j))
  169.           coef(ii + 1, jj + 1) = REAL(matdata(i, j))
  170.          END DO
  171.       END DO
  172.       nn = 2 * numcol
  173.  
  174.       CALL ComplexMatInvert(coef, nn, det, inv)
  175.  
  176.       DO i = 0, numcol - 1
  177.         ii = 2*i
  178.         DO j = 0, numcol - 1
  179.           jj = 2*j
  180.           invary(i, j) = CMPLX(inv(ii, jj),inv(ii + 1, jj))
  181.         END DO
  182.       END DO
  183.       DEALLOCATE (coef,inv, STAT=iErr)
  184.       END !SUB CMatInvert
  185.  
  186.  
  187.       SUBROUTINE CMatPrint (m1, m, n, matname)
  188.       INCLUDE 'STDHDR.FOR'
  189.       COMPLEX m1(0:maxr,0:maxc)
  190.       INTEGER m,n
  191.       CHARACTER * 80 matname
  192.       INTEGER i, j
  193.  
  194.       WRITE (*,*) matname
  195.       WRITE (*,*)
  196.       DO i = 0, m - 1
  197.         DO j = 0, n - 1
  198.           WRITE (*,432) '[',REAL(m1(i,j)), IMAG(m1(i,j)),'j]'
  199.         END DO
  200.         WRITE (*,*)
  201.       END DO
  202.       WRITE (*,*)
  203. 432    FORMAT (1X, A1, F7.2, 2X, F7.2, A2 \)
  204.       END !SUBROUTINE
  205.  
  206.  
  207.  
  208.       SUBROUTINE CMatProd (m1, m2, l, m, n, m3)
  209.       INCLUDE 'STDHDR.FOR'
  210.       COMPLEX m1(0:maxr,0:maxc), m2(0:maxr, 0:maxc), m3(0:maxr, 0:maxc)
  211.       INTEGER l,m,n
  212.       COMPLEX temp,temp2
  213.       INTEGER i,j,k
  214.  
  215.       DO i = 0, l - 1
  216.         DO j = 0, n - 1
  217.           m3(i,j) = CMPLX(0,0)
  218.           DO k = 0, m - 1
  219.             CALL ComplexMath(m1(i, k), 2, m2(k, j), temp)
  220.             CALL ComplexMath(m3(i, j), 0, temp, temp2)
  221.             m3(i,j) = temp2
  222.           END DO
  223.         END DO
  224.       END DO
  225.       END !SUBROUTINE CMatProd
  226.  
  227.  
  228.  
  229.       SUBROUTINE CMatScalarProd (m1, sval, numrow, numcol, m2)
  230.       INCLUDE 'STDHDR.FOR'
  231.       COMPLEX m1(0:maxr,0:maxc), sval, m2(0:maxr, 0:maxc)
  232.       INTEGER i, j
  233.  
  234.       DO i = 0, numrow - 1
  235.         DO j = 0, numcol - 1
  236.           CALL ComplexMath(sval, 2, m1(i, j), m2(i, j))
  237.         END DO
  238.       END DO
  239.       END !SUB CMatScalarProd
  240.  
  241.  
  242.  
  243.       SUBROUTINE CMatTranspose (m1, numrow, numcol, m2)
  244.       INCLUDE 'STDHDR.FOR'
  245.       COMPLEX m1(0:maxr,0:maxc),  m2(0:maxr, 0:maxc)
  246.       INTEGER   i, j
  247.  
  248.       DO i = 0, numrow - 1
  249.         DO j = 0, numcol - 1
  250.           m2(j, i) = m1(i, j)
  251.         END DO
  252.       END DO
  253.       END !SUBROUTINE CMatTranspose
  254.  
  255.  
  256.  
  257.  
  258.  
  259.       SUBROUTINE ComplexMath (c1, op, c2, result)
  260.       COMPLEX c1,c2,result
  261.       INTEGER op
  262.       REAL tempr, tempi
  263.  
  264.       SELECT CASE (op)
  265.          CASE (0)
  266.            result = c1 + c2
  267.          CASE (1)
  268.            result = c1 - c2
  269.          CASE (2)
  270.            tempr = REAL(c1) * REAL(c2) - IMAG(c1) * IMAG(c2)
  271.            tempi = REAL(c1) * IMAG(c2) + IMAG(c1) * REAL(c2)
  272.            result = CMPLX( tempr, tempi)
  273.          CASE (3)
  274.            tempr = (REAL(c1)*REAL(c2)+IMAG(c1)*IMAG(c2))/
  275.      +                    ((REAL(c2) ** 2) + (IMAG(c2) ** 2))
  276.            tempi = (IMAG(c1) * REAL(c2) - REAL(c1) * IMAG(c2))/
  277.      +                    ((REAL(c2) ** 2) + (IMAG(c2) ** 2))
  278.            result = CMPLX( tempr, tempi)
  279.         END SELECT
  280.  
  281.       END !SUBROUTINE ComplexMath
  282.  
  283.  
  284.