Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Homework Help: 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
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook

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