1. Limited time only! Sign up for a free 30min personal tutor trial with Chegg Tutors
    Dismiss Notice
Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Using fortran to Solve linear equations using gauss elimination and back substitution

  1. Feb 16, 2012 #1
    1. The problem statement, all variables and given/known data

    What will be the value of the variable ipvt and AMD when the input value of lud to the following subroutine (solver) is zero ?

    subroutine Solver (A,B,N,lud,Z)
    integer lda,N,ipvt(1000),info,lud,IDAMAX
    &j,k,kp1,l,nm1,kb

    double precision A(1000,1000),B(1000),Z(1000),t,AMD(1000,1000)
    common/ludcmp/ipvt,AMD

    nm1=N-1
    do 5 i=1,N
    Z(i)=B(i)
    5 continue

    if (lud.eq.0) goto 99

    do 6 i=1,N
    do 6 j=1,N
    AMD(i,j)=A(i,j)
    6 continue

    info = 0
    IF (nm1 .LT. 1) GO TO 70
    DO 60 k = 1, nm1
    kp1 = k + 1
    C
    C FIND l = PIVOT INDEX
    C
    l = IDAMAX(N-k+1,AMD(k,k),1) + k - 1
    ipvt(k) = l
    C
    C ZERO PIVOT IMPLIES THIS COLUMN ALREADY TRIANGULARIZED
    C
    IF (AMD(l,k) .EQ. 0.0D0) GO TO 40
    C
    C INTERCHANGE IF NECESSARY
    C
    IF (l .EQ.k) GO TO 10
    t = AMD(l,k)
    AMD(l,k) = AMD(k,k)
    AMD(k,k) = t
    10 CONTINUE
    C
    C COMPUTE MULTIPLIERS
    C
    t = -1.0D0/AMD(k,k)
    CALL DSCAL(N-k,t,AMD(k+1,k),1)
    C
    C ROW ELIMINATION WITH COLUMN INDEXING
    C
    DO 30 j= kp1, N
    t = AMD(l,j)
    IF (l .EQ.k) GO TO 20
    AMD(l,j) = AMD(k,j)
    AMD(k,j) = t
    20 CONTINUE
    CALL DAXPY(N-k,t,AMD(k+1,k),1,AMD(k+1,j),1)
    30 CONTINUE
    GO TO 50
    40 CONTINUE
    info= k
    50 CONTINUE
    60 CONTINUE
    70 CONTINUE
    Ipvt(N) = N
    IF (AMD(N,N) .EQ. 0.0D0) info = N

    if (nm1 .lt. 1) go to 130
    do 120 k = 1, nm1
    l = ipvt(k)
    t = Z(l)
    if (l .eq. k) go to 110
    Z(l) = Z(k)
    Z(k) = t
    110 continue
    call DAXPY (N-k,t,AMD(k+1,k),1,Z(k+1),1)
    120 continue
    130 continue
    c
    c now solve u*x = y
    c
    do 140 kb = 1, N
    k = N + 1 - kb
    Z(k) = Z(k)/AMD(k,k)
    t = -Z(k)
    call DAXPY (k-1,t,AMD(1,k),1,Z(1),1)
    140 continue
    return
    end

    SUBROUTINE DAXPY(N,DA,DX,INCX,DY,INCY)
    * .. Scalar Arguments ..
    DOUBLE PRECISION DA
    INTEGER INCX,INCY,N
    * ..
    * .. Array Arguments ..
    DOUBLE PRECISION DX(*),DY(*)
    * ..
    *
    * Purpose
    * =======
    *
    * DAXPY constant times a vector plus a vector.
    * uses unrolled loops for increments equal to one.
    *
    * Further Details
    * ===============
    *
    * jack dongarra, linpack, 3/11/78.
    * modified 12/3/93, array(1) declarations changed to array(*)
    *
    * =====================================================================
    *
    * .. Local Scalars ..
    INTEGER I,IX,IY,M,MP1
    * ..
    * .. Intrinsic Functions ..
    INTRINSIC MOD
    * ..
    IF (N.LE.0) RETURN
    IF (DA.EQ.0.0d0) RETURN
    IF (INCX.EQ.1 .AND. INCY.EQ.1) THEN
    *
    * code for both increments equal to 1
    *
    *
    * clean-up loop
    *
    M = MOD(N,4)
    IF (M.NE.0) THEN
    DO I = 1,M
    DY(I) = DY(I) + DA*DX(I)
    END DO
    END IF
    IF (N.LT.4) RETURN
    MP1 = M + 1
    DO I = MP1,N,4
    DY(I) = DY(I) + DA*DX(I)
    DY(I+1) = DY(I+1) + DA*DX(I+1)
    DY(I+2) = DY(I+2) + DA*DX(I+2)
    DY(I+3) = DY(I+3) + DA*DX(I+3)
    END DO
    ELSE
    *
    * code for unequal increments or equal increments
    * not equal to 1
    *
    IX = 1
    IY = 1
    IF (INCX.LT.0) IX = (-N+1)*INCX + 1
    IF (INCY.LT.0) IY = (-N+1)*INCY + 1
    DO I = 1,N
    DY(IY) = DY(IY) + DA*DX(IX)
    IX = IX + INCX
    IY = IY + INCY
    END DO
    END IF
    RETURN
    END

    SUBROUTINE DSCAL(N,DA,DX,INCX)
    * .. Scalar Arguments ..
    DOUBLE PRECISION DA
    INTEGER INCX,N
    * ..
    * .. Array Arguments ..
    DOUBLE PRECISION DX(*)
    * ..
    *
    * Purpose
    * =======
    *
    * DSCAL scales a vector by a constant.
    * uses unrolled loops for increment equal to one.
    *
    * Further Details
    * ===============
    *
    * jack dongarra, linpack, 3/11/78.
    * modified 3/93 to return if incx .le. 0.
    * modified 12/3/93, array(1) declarations changed to array(*)
    *
    * =====================================================================
    *
    * .. Local Scalars ..
    INTEGER I,M,MP1,NINCX
    * ..
    * .. Intrinsic Functions ..
    INTRINSIC MOD
    * ..
    IF (N.LE.0 .OR. INCX.LE.0) RETURN
    IF (INCX.EQ.1) THEN
    *
    * code for increment equal to 1
    *
    *
    * clean-up loop
    *
    M = MOD(N,5)
    IF (M.NE.0) THEN
    DO I = 1,M
    DX(I) = DA*DX(I)
    END DO
    IF (N.LT.5) RETURN
    END IF
    MP1 = M + 1
    DO I = MP1,N,5
    DX(I) = DA*DX(I)
    DX(I+1) = DA*DX(I+1)
    DX(I+2) = DA*DX(I+2)
    DX(I+3) = DA*DX(I+3)
    DX(I+4) = DA*DX(I+4)
    END DO
    ELSE
    *
    * code for increment not equal to 1
    *
    NINCX = N*INCX
    DO I = 1,NINCX,INCX
    DX(I) = DA*DX(I)
    END DO
    END IF
    RETURN
    END

    INTEGER FUNCTION IDAMAX(N,DX,INCX)
    * .. Scalar Arguments ..
    INTEGER INCX,N
    * ..
    * .. Array Arguments ..
    DOUBLE PRECISION DX(*)
    * ..
    *
    * Purpose
    * =======
    *
    * IDAMAX finds the index of element having max. absolute value.
    *
    * Further Details
    * ===============
    *
    * jack dongarra, linpack, 3/11/78.
    * modified 3/93 to return if incx .le. 0.
    * modified 12/3/93, array(1) declarations changed to array(*)
    *
    * =====================================================================
    *
    * .. Local Scalars ..
    DOUBLE PRECISION DMAX
    INTEGER I,IX
    * ..
    * .. Intrinsic Functions ..
    INTRINSIC DABS
    * ..
    IDAMAX = 0
    IF (N.LT.1 .OR. INCX.LE.0) RETURN
    IDAMAX = 1
    IF (N.EQ.1) RETURN
    IF (INCX.EQ.1) THEN
    *
    * code for increment equal to 1
    *
    DMAX = DABS(DX(1))
    DO I = 2,N
    IF (DABS(DX(I)).GT.DMAX) THEN
    IDAMAX = I
    DMAX = DABS(DX(I))
    END IF
    END DO
    ELSE
    *
    * code for increment not equal to 1
    *
    IX = 1
    DMAX = DABS(DX(1))
    IX = IX + INCX
    DO I = 2,N
    IF (DABS(DX(IX)).GT.DMAX) THEN
    IDAMAX = I
    DMAX = DABS(DX(IX))
    END IF
    IX = IX + INCX
    END DO
    END IF
    RETURN
    END




    2. Relevant equations
    Solving linear equations using
    1- Gauss elimination
    2- Backsubstitution


    3. The attempt at a solution

    What will be the value of the variable ipvt and AMD when the input value of lud to the subroutine (solver) is zero ?
     
  2. jcsd
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Can you offer guidance or do you also need help?
Draft saved Draft deleted



Similar Discussions: Using fortran to Solve linear equations using gauss elimination and back substitution
Loading...