home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l292 / 1.ddi / GJ.FOR < prev    next >
Encoding:
Text File  |  1990-02-20  |  2.8 KB  |  118 lines
  1.       SUBROUTINE PrintVector(v,n)
  2.       INCLUDE 'stdhdr.for'
  3.       REAL n, v(0:maxv)
  4.       INTEGER i
  5.       DO i=0,n-1
  6.         PRINT *,v(i)
  7.       END DO
  8.       END
  9.  
  10.  
  11.       SUBROUTINE GaussSwap (Swapx ,Swapy )
  12.       REAL Swapx, Swapy,temp
  13.  
  14.        temp = Swapx
  15.        Swapx = Swapy
  16.        Swapy = temp
  17.       RETURN
  18.       END
  19.  
  20.  
  21.  
  22.  
  23.       SUBROUTINE GaussJordan( coefary,  constary,numcol,
  24.      +        solcoef, invary, det)
  25.  
  26.        INCLUDE 'stdhdr.for'
  27.        REAL coefary(0:maxc, 0:maxc), constary(0:maxc )
  28.        REAL det, solcoef( 0:maxc), invary(0:maxc,0:maxc)
  29.        REAL leval, piv, t
  30.        INTEGER  numcol, i, j, lerow, lecol, k, l, iErr
  31.        INTEGER pivlst[ALLOCATABLE](:,:)
  32.        LOGICAL pivchk[ALLOCATABLE](:)
  33.  
  34.        ALLOCATE(pivlst(0:maxc,0:1),STAT=iErr)
  35.        CALL CheckMem(iErr)
  36.        ALLOCATE(pivchk(0:maxc),STAT=iErr)
  37.        CALL CheckMem(iErr)
  38.  
  39.        det = 1.0
  40.        DO i  = 0, numcol  - 1
  41.     pivchk(i) = .FALSE.
  42.      DO j  = 0, numcol  - 1
  43.        invary(i, j ) = coefary(i, j )
  44.      END DO
  45.        END DO
  46.  
  47.        DO i  = 0, numcol  - 1
  48.       leval = 0.0
  49.       DO j    = 0, numcol  - 1
  50.         IF ( .not.(pivchk(j) )) THEN
  51.          DO k  = 0, numcol    - 1
  52.            IF (.not.(pivchk(k))) THEN
  53.            IF (ABS(invary(j, k )) .GT. leval) THEN
  54.               lerow  = j
  55.               lecol  = k
  56.               leval = ABS(invary(j, k ))
  57.             END IF
  58.            END IF
  59.          END DO
  60.         END IF
  61.     END DO
  62.  
  63.      pivchk(lecol ) = .TRUE.
  64.      pivlst(i, 0) = lerow
  65.      pivlst(i, 1) = lecol
  66.  
  67.      IF (lerow  .NE. lecol ) THEN
  68.        det = -det
  69.        DO l  = 0, numcol  - 1
  70.          call GaussSwap( invary(lerow, l ), invary(lecol, l ))
  71.        END DO
  72.        call GaussSwap (constary(lerow ), constary(lecol ))
  73.      END IF
  74.  
  75.      piv = invary(lecol, lecol )
  76.      det = det * piv
  77.      IF (det .GT. 1.0E30) det = 1.0E30 !prevents numeric overflow
  78.      invary(lecol, lecol ) = 1.0
  79.      DO l  = 0, numcol  - 1
  80.        invary(lecol, l ) = invary(lecol, l ) / piv
  81.      END DO
  82.      constary(lecol ) = constary(lecol ) / piv
  83.  
  84.      DO l1    = 0, numcol  - 1
  85.        IF (l1  .NE. lecol ) THEN
  86.          t = invary(l1, lecol )
  87.          invary(l1, lecol ) = 0.0
  88.          DO l  = 0, numcol    - 1
  89.            invary(l1, l ) = invary(l1, l ) - invary(lecol, l ) * t
  90.          END DO
  91.          constary(l1 ) = constary(l1 ) - constary(lecol ) * t
  92.        END IF
  93.      END DO
  94.        END DO
  95.  
  96.       DO i  = 0, numcol  - 1
  97.      l  = numcol  - i  - 1
  98.      IF (pivlst(l, 0) .NE. pivlst(l, 1)) THEN
  99.        lerow  = pivlst(l, 0)
  100.        lecol  = pivlst(l, 1)
  101.        DO k  = 0, numcol  - 1
  102.          CALL GaussSwap (invary(k, lerow ), invary(k, lecol ))
  103.        END DO
  104.       END IF
  105.       END DO
  106.  
  107.       DO i  = 0, numcol  - 1
  108.       solcoef(i ) = constary(i )
  109.       END DO
  110.       DEALLOCATE(pivlst,STAT=iErr)
  111.       CALL CheckMem(iErr)
  112.       DEALLOCATE(pivchk,STAT=iErr)
  113.       CALL CheckMem(iErr)
  114.  
  115.       END  !GaussJordan
  116.  
  117.  
  118.