- #1
mmokhtar
- 1
- 0
Homework Statement
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
Homework Equations
Solving linear equations using
1- Gauss elimination
2- Backsubstitution
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 ?