Revision a137ab8c6d6d0a8110a234da6f4de49aa7b548d9 authored by Roger Koenker on 31 August 2015, 23:25:56 UTC, committed by Gabor Csardi on 31 August 2015, 23:25:56 UTC
1 parent a6e1bd2
crq.f
C Comments, bug reports, etc. should be sent to: portnoy@stat.uiuc.edu
C
SUBROUTINE CRQ(M,N,MPLUS,N2,X,Y,C,TZERO,MAXW,STEP,IFT,H,XH,
* WA,WB,WC,WD,WE,WF,IA,NSOL,SOL,LSOL,ICEN,TCEN,LCEN)
C
C M = number of Observations
C N = Number of Parameters
C MPLUS = row dim of X, WA, and WC .GE. M+1
C N2 = N+2
C X is the X matrix (MPLUS BY N) (includes censored obs)
C Y is the y vector (M) (includes censored obs)
C C is the censoring integer(M) vector, 1 if obs is censored,
C 0 otherwise
C TZERO is the initial tau value
C MAXW ( .LE. MPLUS ): Maximal number of calls to weighted RQ
C for possible degeneracies. If exceeded (IFT = 7), dither.
C If MAXW .LE. 0, don't pivot -- use only weighted rq
C STEP: Step size for weighted rq at degeneracy
C IFT exit code:
C 0: ok
C 1: M < 0 or N < 0 OR M < N
C 2: MPLUS .LT. M + 1 or N2 .NE. N+2
C 3: MAXW > MPLUS
C 4: less than N noncensored obs above the tau=0 soln
C 5: possible degeneracy, rq called at tau + step, see IA
C 6: LSOL becomes greater than NSOL
C 7: MAXW exceeded
C 8: weighted rq tries to include infinite basis element
C H is an integer work vector (N) ( = basis indices )
C XH is a double precision work array (N by N)
C WA is a double precision work array (MPLUS by N)
C WB is a double precision work vector (MPLUS)
C WC is a double precision work array (MPLUS by N+2)
C WD is a double precision work vector (MPLUS)
C WE is a double precision work vector (MPLUS)
C WF is a double precision work vector (N)
C IA is an integer work vector (MPLUS)
C NSOL is an (estimated) row dimension of the solution array
C if NSOL < M, solutions returned only at tau = i/(NSOL-1)
C if all solutions are desired, NSOL = 3*M is a good choice
C on output:
C SOL is the coefficient solution array (N+2 by NSOL)
C first and last rows give tau bkpts and quantile
C rows 2:(N+1) give the beta coefficients
C LSOL is the actual dimension of the solution arrays
C LSOL = NSOL if NSOL < M, else LSOL .LE. NSOL
C ICEN (M) indicates status of censored observations
C = 0 if uncensored (or uncrossed censored while pivoting)
C = 1 if crossed censored
C = 2 if deleted as below tau = 0 solution
C = 3 if above last (maximal tau) solution
C TCEN (M) are the corresponding tau value where censored obs is
C firstcrossed: weight = (tau - TCEN(I))/(1 - TCEN(I))
C TCEN = 1 if censored obs is never crossed (ICEN = 3)
C TCEN = 0 if censored obs is deleted (ICEN = 2)
C LCEN = number of censored observations
C WD = list of first MPLUS tau values at which degeneracy
C may have occurred (and weighted RQ was called)
C H(1) = number of calls to weighted RQ (apparent degeneracy)
C
IMPLICIT REAL*8(A-H,O-Z)
INTEGER H(N),ICEN(M),C(M),IA(MPLUS)
LOGICAL APC, DEG
DOUBLE PRECISION ONE
DIMENSION X(M,N),Y(M),SOL(N2,NSOL),WC(MPLUS,N2),XH(N,N)
DIMENSION WA(MPLUS,N),WB(MPLUS),WD(MPLUS),WE(MPLUS)
DIMENSION F(2),WF(N),TCEN(M)
DATA BIG/1.D17/
DATA ZERO/0.00D0/
DATA HALF/0.50D0/
DATA ONE/1.00D0/
DATA TWO/2.00D0/
C
C CHECK DIMENSION PARAMETERS
C
DEG = .FALSE.
IFT=0
N1 = NSOL-1
LCEN = 0
DO 2 I = 1,M
IF(C(I) .EQ. 1) LCEN = LCEN + 1
2 CONTINUE
IF(MPLUS .LT. M+1) IFT = 2
IF(N2 .NE. N+2) IFT = 2
IF(M .LE. ZERO .OR. N. LE. ZERO .OR. M .LT. N) IFT = 1
IF(MAXW .GT. MPLUS) IFT = 3
IF(IFT .NE. 0) RETURN
C
C INITIALIZATION
C
TOLER = 1.D-11
TOL1 = 10.0D0*TOLER
IF(TZERO .LE. ZERO .OR. TZERO .GT. STEP) TZERO = STEP
DO 5 I=1,M
ICEN(I) = 0
TCEN(I) = ONE
WB(I) = Y(I)
DO 5 J=1,N
5 WA(I,J) = X(I,J)
MA = M
C
C GET TAU = O SOLUTION
C DELETE CENSORED OBS BELOW TAU = 0 SOLN
C
10 IF(MA .LE. N) THEN
IFT = 4
RETURN
ENDIF
CALL RQ1(MA,N,MPLUS,N2,WA,WB,TZERO,TOLER,IFT1,
* WF,WE,IA,WC,WD)
K = 0
KL = 0
C
C CHECK FOR CENSORED OBS BELOW TAU = 0 SOLN
C
L = 0
DO 20 I=1,M
IF(ICEN(I) .EQ. 2) GO TO 20
L = L + 1
IF( C(I) .EQ. 1 .AND. WE(L) .LE. TOL1 ) THEN
K = K + 1
ICEN(I) = 2
GO TO 20
ENDIF
IF( DABS(WE(L)) .GE. TOL1) GO TO 20
KL = KL + 1
IF(KL .LE. N) H(KL) = I
20 CONTINUE
IF( K .EQ. 0) GO TO 40
MA = MA - K
L = 0
C
C DELETE CENSORED OBS BELOW TAU = 0 SOLN, AND RETRY
C
DO 30 I=1,M
IF(ICEN(I) .EQ. 2) GO TO 30
L = L + 1
WB(L) = Y(I)
DO 25 J=1,N
25 WA(L,J) = X(I,J)
30 CONTINUE
GO TO 10
40 CONTINUE
C
C SAVE INITIAL SOLUTION, AND INITIALIZE PIVOTING LOOP
C
SOL(1,1) = ZERO
DO 50 J = 1,N
50 SOL(J+1,1) = WF(J)
NWRQ = 0
APC = .FALSE.
LSOL = 2
TAU = TZERO
C
C COMPUTE NEXT TAU
C FIRST CHECK FOR DEGENERACY AT TZERO
C
IF (KL .GT. N) GO TO 500
C
C GET X(H,)
C
DO 70 K = 1,N
I= H(K)
DO 60 J = 1,N
60 XH(K,J) = X(I,J)
70 CONTINUE
C
C GET X(H,) INVERSE
C
CALL DGECO(XH,N,N,IA,V,WC)
IF(V. LT. TOLER) GO TO 500
JOB = 01
CALL DGEDI(XH,N,N,IA,F,WC,JOB)
C
C GET RESIDUALS
C
DO 90 I=1,M
S = Y(I)
DO 80 J = 1,N
80 S = S - WF(J)*X(I,J)
90 WE(I) = S
C
C BEGIN PIVOTING; WD = GRAD UPPER BD, IA = SIGN(GRAD DENOM)
C
200 CONTINUE
CALL GRAD(X,Y,M,N,H,ICEN,TCEN,XH,WE,TOLER,IA,WC,WD)
KL = 0
KM = 1
S = WD(1)
IF(N .EQ. 1) GO TO 230
DO 210 J=2,N
210 S = DMIN1(WD(J),S)
DO 220 J=1,N
IF(WD(J) .GE. S + TOLER) GO TO 220
KL = KL + 1
KM = J
220 CONTINUE
230 CONTINUE
C
C IF ALL POS RESIDS CENSORED, RETURN; ELSE CHECK FOR DEGEN.
C
IF(APC) THEN
SOL(1,LSOL) = DMAX1(S, TAU)
DO 240 J=1,N
240 SOL(J+1,LSOL) = WF(J)
GO TO 600
ENDIF
C
C IF DEGENERACY, CALL WEIGHTED RQ
C
IF (KL .GT. 1) GO TO 500
C
C CHECK FOR INFEASIBILITY CAUSED BY REWEIGHTING
C
IF(S .LT. TAU + TOL1) THEN
LSOL = LSOL - 1
TAU = TAU - .8*STEP
GO TO 500
ENDIF
TAU = S
C
C GET NEW OBSERVATION TO ENTER BASIS
C FIRST DEFINE WD = BASIS INDICATOR = 1 IF I IN H(J)
C
DO 250 I=1,M
250 WD(I) = ZERO
DO 260 J=1,N
260 WD(H(J)) = ONE
KN = 0
D = BIG
KIN = 0
DO 300 I=1,M
IF(ICEN(I) .EQ. 2) GO TO 300
S = ZERO
DO 270 J=1,N
270 S = S + X(I,J)*XH(J,KM)
S = IA(KM)*S
IF(DABS(S).LT.TOL1 .OR. (C(I).EQ.1 .AND. ICEN(I).NE.1)
* .OR. WD(I).EQ.ONE) GO TO 300
S = WE(I)/S
IF(S .LT. TOL1) GO TO 300
IF(S .LE. D - TOL1) THEN
D = S
KIN = I
KN = 0
ELSE
IF(S .LE. D + TOL1) KN = 1
ENDIF
300 CONTINUE
IF(KN .EQ. 1) GO TO 500
H(KM) = KIN
C
C UPDATE SOLUTION
C GET NEW X(H,)
C
DO 310 K = 1,N
I = H(K)
DO 310 J = 1,N
310 XH(K,J) = X(I,J)
C
C GET X(H,) INVERSE
C
315 CALL DGECO(XH,N,N,IA,V,WA)
IF(V. LT. TOLER) GO TO 500
JOB = 01
CALL DGEDI(XH,N,N,IA,F,WA,JOB)
C
C GET BETA-HAT, RESIDS
C
DO 340 K = 1,N
S = ZERO
DO 330 J = 1,N
330 S = S + XH(K,J)*Y(H(J))
340 WF(K) = S
345 DO 360 I=1,M
S = Y(I)
DO 350 J = 1,N
350 S = S - WF(J)*X(I,J)
360 WE(I) = S
C
C SAVE SOLUTION
C
365 DO 380 I = (LSOL-1),N1
IF(NSOL.GE.M .OR. TAU .GT. DFLOAT(I)/DFLOAT(N1) - 10*TOL1) THEN
SOL(1,LSOL) = TAU
DO 370 J = 1,N
370 SOL(J+1,LSOL) = WF(J)
LSOL = LSOL + 1
GO TO 390
ENDIF
380 CONTINUE
390 CONTINUE
C
C CHECK FOR DIM(SOL) EXCEEDED
C
IF(LSOL .GT. NSOL) THEN
IFT = 6
GO TO 600
ENDIF
IF(APC) THEN
LSOL = LSOL - 1
GO TO 600
ENDIF
C
C SET WT FOR CROSSED CENSORED OBSERVATIONS
C APC = .TRUE. IF ALL POS RESID ARE UNCROSSED CENSORED
C
APC = .TRUE.
C if the following are left uncommented: still get zero column
C TAUW = TAU - STEP/TWO
C IF(MAXW.GT.ZERO)THEN
C TAUW=TAU
C ENDIF
DO 400 I=1,M
IF(ICEN(I) .EQ. 2) GO TO 400
IF(WE(I).GE.TOL1 .AND. (C(I).EQ.0 .OR. ICEN(I).EQ.1))
* APC = .FALSE.
IF(WE(I).LE.-TOL1 .AND. C(I).EQ.1 .AND. ICEN(I).EQ.0) THEN
C
C at this point, Y(I) is a crossed censored obs (C_i)
C for the grid method TAUW = TAU - STEP*(B2-Y(I))/(B2-B1)
C where B1 = x_i' beta_j and B2 = x_i' beta_(j+1)
C otherwise (for pivot), TAUW = TAU
C
IF(MAXW.LT.0) THEN
B1 = ZERO
B2 = ZERO
DO 396 J=1,N
B1 = B1 + X(I,J) * SOL(J+1,LSOL - 2)
396 B2 = B2 + X(I,J) * SOL(J+1,LSOL - 1)
A1 = (B2 - Y(I))/(B2-B1)
TAUW = TAU - A1*STEP
ELSE
TAUW = TAU
ENDIF
ICEN(I) = 1
TCEN(I) = TAUW
ENDIF
400 CONTINUE
IF(APC) THEN
IF(DEG) THEN
IFT = 5
LSOL = LSOL - 1
GO TO 600
ENDIF
GO TO 200
ENDIF
C
C RETURN IF TAU > 1
C
IF(TAU .GE. ONE - 10.0D0 * TOL1) GO TO 600
IF(DEG) GO TO 500
IF(MAXW .GT. 0) GO TO 200
C
C POSSIBLE DEGENERACY -- USE WEIGHTED RQ1 TO RESOLVE
C
500 NWRQ = NWRQ + 1
IF(MAXW .GT. 0 .AND. NWRQ .GT. MAXW) THEN
IFT = 7
LSOL = LSOL - 1
GO TO 600
ENDIF
IF(NWRQ .LE. NSOL) SOL(N+2,NWRQ) = TAU
TAU = TAU + STEP
IF(TAU .GE. ONE) TAU = ONE - TOL1
L = 0
L1 = 0
DO 506 J=1,N
506 WF(J) = ZERO
DO 530 I = 1,M
IF(ICEN(I) .EQ. 0) THEN
IF (C(I) .EQ. 0) THEN
L1 = L1 + 1
WB(L1) = Y(I)
DO 510 J=1,N
510 WA(L1,J) = X(I,J)
ELSE
DO 514 J = 1,N
514 WF(J) = WF(J) + X(I,J)
ENDIF
ENDIF
IF(ICEN(I) .EQ. 1) THEN
L1 = L1 + 1
L = L+1
W = (TAU - TCEN(I))/(ONE - TCEN(I))
WB(L1) = W * Y(I)
DO 520 J = 1,N
WA(L1,J) = W * X(I,J)
520 WF(J) = WF(J) + (ONE - W) * X(I,J)
ENDIF
530 CONTINUE
MAL = L1+1
DO 534 J=1,N
534 WA(MAL,J) = WF(J)
WB(MAL) = BIG
CALL RQ1(MAL,N,MPLUS,N2,WA,WB,TAU,TOLER,IFT1,WF,WE,IA,WC,WD)
DEG = .FALSE.
IF(DABS(WE(MAL)) .LE. .0001) THEN
IFT = 8
LSOL = LSOL - 1
GO TO 600
ENDIF
L = 0
K = 0
DO 550 I=1,M
IF(ICEN(I) .EQ. 2) GO TO 550
L = L+1
IF(DABS(WE(L)) .GE. TOL1) GO TO 550
K = K+1
IF(K .LE. N) THEN
H(K) = I
DO 540 J=1,N
540 XH(K,J) = X(I,J)
ENDIF
550 CONTINUE
IF(K .LT. N) THEN
IFT = 8
LSOL = LSOL - 1
GO TO 600
ENDIF
IF(K .GT. N) DEG = .TRUE.
GO TO 345
C
C DEFINE OUTPUT AND RETURN
C
600 H(1) = NWRQ
L = MIN0(NWRQ,MPLUS)
DO 610 I=1,L
610 WD(I) = SOL(N+2,I)
DO 620 I=1,M
IF(ICEN(I) .EQ. 2) GO TO 620
IF(C(I) .EQ.1 .AND. TCEN(I) .EQ. ONE) ICEN(I) = 3
620 CONTINUE
V = DFLOAT(MA)
DO 650 J = 1,N
S = ZERO
DO 640 I = 1,M
IF(ICEN(I) .EQ. 2) GO TO 640
S = S + X(I,J)
640 CONTINUE
650 WF(J) = S/V
DO 670 I = 1,LSOL
S = ZERO
DO 660 J = 1,N
660 S = S + WF(J) * SOL(J+1,I)
670 SOL(N+2,I) = S
RETURN
END
SUBROUTINE GRAD(X,Y,M,N,H,ICEN,TCEN,XH,
* R,TOL,IFLAG,WA,GUP)
C
C X matrix (M BY N)
C Y vector (M)
C M = Number of Observations
C N = Number of Parameters
C H = basis, integer(N) vector
C ICEN (integer(M)) = 2 for deleted obs
C = 1 for crossed non-censored obs
C = 0 otherwise
C TCEN(M) are the corresponding tau values where a censored obs is
C first crossed (ICEN=1): weight = (tau - TCEN(I))/(1 - TCEN(I))
C TCEN = 1 if censored obs is never crossed
C XH = (N by N) X(H,) inverse matrix
C R(M) = residuals
C TOL = zero tolerance (1.D-10 by default)
C IFLAG (M) work vector
C WA (M by N) work array
C returns:
C GUP(N) = upper bounds on tau
C IFLAG(1:N) = sign(grad denom)
C = +1 if bound from lower grad bound
C = -1 if bound from upper grad bound
C
IMPLICIT REAL*8(A-H,O-Z)
INTEGER H(N),ICEN(M),IFLAG(*)
DOUBLE PRECISION ONE
DIMENSION X(M,N),Y(M),GUP(N),XH(N,N)
DIMENSION R(M),TCEN(M),WA(M,N)
DATA ZERO/0.00D0/, ONE/1.0d0/
C
C GET X'(XH-INVERSE)
C
DO 80 I = 1,M
IF(ICEN(I) .EQ. 2) GO TO 80
DO 70 J = 1,N
A = ZERO
DO 60 K = 1,N
60 A = A + X(I,K)*XH(K,J)
70 WA(I,J) = A
80 CONTINUE
C
C GET GRADIENT BOUNDS
C FIRST SET IFLAG = 1 FOR BASIS INDICES (TEMPORARY)
C
T = ZERO
DO 90 I = 1,M
90 IFLAG(I) = 0
DO 95 J = 1,N
95 IFLAG(H(J)) = 1
DO 120 J=1,N
A = ZERO
B = ZERO
C = ZERO
D = ZERO
DO 100 I = 1,M
IF(ICEN(I) .EQ. 2) GO TO 100
IF(ICEN(I) .EQ. 0) THEN
IF(R(I) .GT. TOL) A = A + WA(I,J)
IF(R(I) .LT. - TOL) B = B + WA(I,J)
GO TO 100
ENDIF
IF(IFLAG(I).EQ.1) GO TO 100
C
C IFLAG = 0: NOT BASIS AND ICEN = 1: CROSSED CENSORED
C
IF(R(I) .LT. -TOL) THEN
T = TCEN(I)/(ONE - TCEN(I))
C = C - T*WA(I,J)
GO TO 100
ENDIF
IF(R(I) .GT. TOL) THEN
D = D - WA(I,J)
ENDIF
100 CONTINUE
D = D - C
L = ICEN(H(J))
IF(L .NE. 0) T = TCEN(H(J))/( ONE - TCEN(H(J)) )
S = DFLOAT(L)*(T + ONE) - ONE
E1 = A + B - D - S
E2 = A + B - D + ONE
IF(E1 .GT. ZERO) THEN
GUP(J) = (B + C - S)/E1
IFLAG(J+M) = 1
GO TO 120
ENDIF
IF(E2 .LT. ZERO) THEN
GUP(J) = (B + C)/E2
IFLAG(J+M) = -1
GO TO 120
ENDIF
GUP(J) = - ONE
120 CONTINUE
DO 130 J = 1,N
130 IFLAG(J) = IFLAG(J+M)
RETURN
END
subroutine dgedi(a,lda,n,ipvt,det,work,job)
integer lda,n,ipvt(*),job
double precision a(lda,*),det(2),work(*)
c
c dgedi computes the determinant and inverse of a matrix
c using the factors computed by dgeco or dgefa.
c
c on entry
c
c a double precision(lda, n)
c the output from dgeco or dgefa.
c
c lda integer
c the leading dimension of the array a .
c
c n integer
c the order of the matrix a .
c
c ipvt integer(n)
c the pivot vector from dgeco or dgefa.
c
c work double precision(n)
c work vector. contents destroyed.
c
c job integer
c = 11 both determinant and inverse.
c = 01 inverse only.
c = 10 determinant only.
c
c on return
c
c a inverse of original matrix if requested.
c otherwise unchanged.
c
c det double precision(2)
c determinant of original matrix if requested.
c otherwise not referenced.
c determinant = det(1) * 10.0**det(2)
c with 1.0 .le. dabs(det(1)) .lt. 10.0
c or det(1) .eq. 0.0 .
c
c error condition
c
c a division by zero will occur if the input factor contains
c a zero on the diagonal and the inverse is requested.
c it will not occur if the subroutines are called correctly
c and if dgeco has set rcond .gt. 0.0 or dgefa has set
c info .eq. 0 .
c
c linpack. this version dated 08/14/78 .
c cleve moler, university of new mexico, argonne national lab.
c
c subroutines and functions
c
c blas daxpy,dscal,dswap
c fortran dabs,mod
c
c internal variables
c
double precision t
double precision ten
integer i,j,k,kb,kp1,l,nm1
c
c
c compute determinant
c
if (job/10 .eq. 0) go to 70
det(1) = 1.0d0
det(2) = 0.0d0
ten = 10.0d0
do 50 i = 1, n
if (ipvt(i) .ne. i) det(1) = -det(1)
det(1) = a(i,i)*det(1)
c ...exit
if (det(1) .eq. 0.0d0) go to 60
10 if (dabs(det(1)) .ge. 1.0d0) go to 20
det(1) = ten*det(1)
det(2) = det(2) - 1.0d0
go to 10
20 continue
30 if (dabs(det(1)) .lt. ten) go to 40
det(1) = det(1)/ten
det(2) = det(2) + 1.0d0
go to 30
40 continue
50 continue
60 continue
70 continue
c
c compute inverse(u)
c
if (mod(job,10) .eq. 0) go to 150
do 100 k = 1, n
a(k,k) = 1.0d0/a(k,k)
t = -a(k,k)
call dscal(k-1,t,a(1,k),1)
kp1 = k + 1
if (n .lt. kp1) go to 90
do 80 j = kp1, n
t = a(k,j)
a(k,j) = 0.0d0
call daxpy(k,t,a(1,k),1,a(1,j),1)
80 continue
90 continue
100 continue
c
c form inverse(u)*inverse(l)
c
nm1 = n - 1
if (nm1 .lt. 1) go to 140
do 130 kb = 1, nm1
k = n - kb
kp1 = k + 1
do 110 i = kp1, n
work(i) = a(i,k)
a(i,k) = 0.0d0
110 continue
do 120 j = kp1, n
t = work(j)
call daxpy(n,t,a(1,j),1,a(1,k),1)
120 continue
l = ipvt(k)
if (l .ne. k) call dswap(n,a(1,k),1,a(1,l),1)
130 continue
140 continue
150 continue
return
end
subroutine dgeco(a,lda,n,ipvt,rcond,z)
integer lda,n,ipvt(*)
double precision a(lda,*),z(*)
double precision rcond
c
c dgeco factors a double precision matrix by gaussian elimination
c and estimates the condition of the matrix.
c
c
c if rcond is not needed, dgefa is slightly faster.
c to solve a*x = b , follow dgeco by dgesl.
c to compute inverse(a)*c , follow dgeco by dgesl.
c to compute determinant(a) , follow dgeco by dgedi.
c to compute inverse(a) , follow dgeco by dgedi.
c
c on entry
c
c a double precision(lda, n)
c the matrix to be factored.
c
c lda integer
c the leading dimension of the array a .
c
c n integer
c the order of the matrix a .
c
c on return
c
c a an upper triangular matrix and the multipliers
c which were used to obtain it.
c the factorization can be written a = l*u where
c l is a product of permutation and unit lower
c triangular matrices and u is upper triangular.
c
c ipvt integer(n)
c an integer vector of pivot indices.
c
c rcond double precision
c an estimate of the reciprocal condition of a .
c for the system a*x = b , relative perturbations
c in a and b of size epsilon may cause
c relative perturbations in x of size epsilon/rcond .
c if rcond is so small that the logical expression
c 1.0 + rcond .eq. 1.0
c is true, then a may be singular to working
c precision. in particular, rcond is zero if
c exact singularity is detected or the estimate
c underflows.
c
c z double precision(n)
c a work vector whose contents are usually unimportant.
c if a is close to a singular matrix, then z is
c an approximate null vector in the sense that
c norm(a*z) = rcond*norm(a)*norm(z) .
c
c linpack. this version dated 08/14/78 .
c cleve moler, university of new mexico, argonne national lab.
c
c subroutines and functions
c
c linpack dgefa
c blas daxpy,ddot,dscal,dasum
c fortran dabs,dmax1,dsign
c
c internal variables
c
double precision ddot,ek,t,wk,wkm
double precision anorm,s,dasum,sm,ynorm
integer info,j,k,kb,kp1,l
c
c
c compute 1-norm of a
c
anorm = 0.0d0
do 10 j = 1, n
anorm = dmax1(anorm,dasum(n,a(1,j),1))
10 continue
c
c factor
c
call dgefa(a,lda,n,ipvt,info)
c
c rcond = 1/(norm(a)*(estimate of norm(inverse(a)))) .
c estimate = norm(z)/norm(y) where a*z = y and trans(a)*y = e .
c trans(a) is the transpose of a . the components of e are
c chosen to cause maximum local growth in the elements of w where
c trans(u)*w = e . the vectors are frequently rescaled to avoid
c overflow.
c
c solve trans(u)*w = e
c
ek = 1.0d0
do 20 j = 1, n
z(j) = 0.0d0
20 continue
do 100 k = 1, n
if (z(k) .ne. 0.0d0) ek = dsign(ek,-z(k))
if (dabs(ek-z(k)) .le. dabs(a(k,k))) go to 30
s = dabs(a(k,k))/dabs(ek-z(k))
call dscal(n,s,z,1)
ek = s*ek
30 continue
wk = ek - z(k)
wkm = -ek - z(k)
s = dabs(wk)
sm = dabs(wkm)
if (a(k,k) .eq. 0.0d0) go to 40
wk = wk/a(k,k)
wkm = wkm/a(k,k)
go to 50
40 continue
wk = 1.0d0
wkm = 1.0d0
50 continue
kp1 = k + 1
if (kp1 .gt. n) go to 90
do 60 j = kp1, n
sm = sm + dabs(z(j)+wkm*a(k,j))
z(j) = z(j) + wk*a(k,j)
s = s + dabs(z(j))
60 continue
if (s .ge. sm) go to 80
t = wkm - wk
wk = wkm
do 70 j = kp1, n
z(j) = z(j) + t*a(k,j)
70 continue
80 continue
90 continue
z(k) = wk
100 continue
s = 1.0d0/dasum(n,z,1)
call dscal(n,s,z,1)
c
c solve trans(l)*y = w
c
do 120 kb = 1, n
k = n + 1 - kb
if (k .lt. n) z(k) = z(k) + ddot(n-k,a(k+1,k),1,z(k+1),1)
if (dabs(z(k)) .le. 1.0d0) go to 110
s = 1.0d0/dabs(z(k))
call dscal(n,s,z,1)
110 continue
l = ipvt(k)
t = z(l)
z(l) = z(k)
z(k) = t
120 continue
s = 1.0d0/dasum(n,z,1)
call dscal(n,s,z,1)
c
ynorm = 1.0d0
c
c solve l*v = y
c
do 140 k = 1, n
l = ipvt(k)
t = z(l)
z(l) = z(k)
z(k) = t
if (k .lt. n) call daxpy(n-k,t,a(k+1,k),1,z(k+1),1)
if (dabs(z(k)) .le. 1.0d0) go to 130
s = 1.0d0/dabs(z(k))
call dscal(n,s,z,1)
ynorm = s*ynorm
130 continue
140 continue
s = 1.0d0/dasum(n,z,1)
call dscal(n,s,z,1)
ynorm = s*ynorm
c
c solve u*z = v
c
do 160 kb = 1, n
k = n + 1 - kb
if (dabs(z(k)) .le. dabs(a(k,k))) go to 150
s = dabs(a(k,k))/dabs(z(k))
call dscal(n,s,z,1)
ynorm = s*ynorm
150 continue
if (a(k,k) .ne. 0.0d0) z(k) = z(k)/a(k,k)
if (a(k,k) .eq. 0.0d0) z(k) = 1.0d0
t = -z(k)
call daxpy(k-1,t,a(1,k),1,z(1),1)
160 continue
c make znorm = 1.0
s = 1.0d0/dasum(n,z,1)
call dscal(n,s,z,1)
ynorm = s*ynorm
c
if (anorm .ne. 0.0d0) rcond = ynorm/anorm
if (anorm .eq. 0.0d0) rcond = 0.0d0
return
end
subroutine dgefa(a,lda,n,ipvt,info)
c use numerical_libraries
integer lda,n,ipvt(*),info
double precision a(lda,*)
c
c dgefa factors a double precision matrix by gaussian elimination.
c
c dgefa is usually called by dgeco, but it can be called
c directly with a saving in time if rcond is not needed.
c (time for dgeco) = (1 + 9/n)*(time for dgefa) .
c
c on entry
c
c a double precision(lda, n)
c the matrix to be factored.
c
c lda integer
c the leading dimension of the array a .
c
c n integer
c the order of the matrix a .
c
c on return
c
c a an upper triangular matrix and the multipliers
c which were used to obtain it.
c the factorization can be written a = l*u where
c l is a product of permutation and unit lower
c triangular matrices and u is upper triangular.
c
c ipvt integer(n)
c an integer vector of pivot indices.
c
c info integer
c = 0 normal value.
c = k if u(k,k) .eq. 0.0 . this is not an error
c condition for this subroutine, but it does
c indicate that dgesl or dgedi will divide by zero
c if called. use rcond in dgeco for a reliable
c indication of singularity.
c
c linpack. this version dated 08/14/78 .
c cleve moler, university of new mexico, argonne national lab.
c
c subroutines and functions
c
c blas daxpy,dscal,idamax
c
c internal variables
c
double precision t
c integer idamax,j,k,kp1,l,nm1
integer j,k,kp1,l,nm1
c
c
c gaussian elimination with partial pivoting
c
info = 0
nm1 = n - 1
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,a(k,k),1) + k - 1
ipvt(k) = l
c
c zero pivot implies this column already triangularized
c
if (a(l,k) .eq. 0.0d0) go to 40
c
c interchange if necessary
c
if (l .eq. k) go to 10
t = a(l,k)
a(l,k) = a(k,k)
a(k,k) = t
10 continue
c
c compute multipliers
c
t = -1.0d0/a(k,k)
call dscal(n-k,t,a(k+1,k),1)
c
c row elimination with column indexing
c
do 30 j = kp1, n
t = a(l,j)
if (l .eq. k) go to 20
a(l,j) = a(k,j)
a(k,j) = t
20 continue
call daxpy(n-k,t,a(k+1,k),1,a(k+1,j),1)
30 continue
go to 50
40 continue
info = k
50 continue
60 continue
70 continue
ipvt(n) = n
if (a(n,n) .eq. 0.0d0) info = n
return
end
subroutine daxpy(n,da,dx,incx,dy,incy)
c constant times a vector plus a vector.
c uses unrolled loops for increments equal to one.
c jack dongarra, linpack, 3/11/78.
c modified 12/3/93, array(1) declarations changed to array(*)
c
double precision dx(*),dy(*),da
integer i,incx,incy,ix,iy,m,mp1,n
c
if(n.le.0)return
if (da .eq. 0.0d0) return
if(incx.eq.1.and.incy.eq.1)go to 20
c
c code for unequal increments or equal increments
c not equal to 1
c
ix = 1
iy = 1
if(incx.lt.0)ix = (-n+1)*incx + 1
if(incy.lt.0)iy = (-n+1)*incy + 1
do 10 i = 1,n
dy(iy) = dy(iy) + da*dx(ix)
ix = ix + incx
iy = iy + incy
10 continue
return
c
c code for both increments equal to 1
c
c
c clean-up loop
c
20 m = mod(n,4)
if( m .eq. 0 ) go to 40
do 30 i = 1,m
dy(i) = dy(i) + da*dx(i)
30 continue
if( n .lt. 4 ) return
40 mp1 = m + 1
do 50 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)
50 continue
return
end
double precision function ddot(n,dx,incx,dy,incy)
c
c forms the dot product of two vectors.
c uses unrolled loops for increments equal to one.
c jack dongarra, linpack, 3/11/78.
c modified 12/3/93, array(1) declarations changed to array(*)
c
double precision dx(*),dy(*),dtemp
integer i,incx,incy,ix,iy,m,mp1,n
c
ddot = 0.0d0
dtemp = 0.0d0
if(n.le.0)return
if(incx.eq.1.and.incy.eq.1)go to 20
c
c code for unequal increments or equal increments
c not equal to 1
c
ix = 1
iy = 1
if(incx.lt.0)ix = (-n+1)*incx + 1
if(incy.lt.0)iy = (-n+1)*incy + 1
do 10 i = 1,n
dtemp = dtemp + dx(ix)*dy(iy)
ix = ix + incx
iy = iy + incy
10 continue
ddot = dtemp
return
c
c code for both increments equal to 1
c
c
c clean-up loop
c
20 m = mod(n,5)
if( m .eq. 0 ) go to 40
do 30 i = 1,m
dtemp = dtemp + dx(i)*dy(i)
30 continue
if( n .lt. 5 ) go to 60
40 mp1 = m + 1
do 50 i = mp1,n,5
dtemp = dtemp + dx(i)*dy(i) + dx(i + 1)*dy(i + 1) +
* dx(i + 2)*dy(i + 2) + dx(i + 3)*dy(i + 3) + dx(i + 4)*dy(i + 4)
50 continue
60 ddot = dtemp
return
end
subroutine dscal(n,da,dx,incx)
c scales a vector by a constant.
c uses unrolled loops for increment equal to one.
c jack dongarra, linpack, 3/11/78.
c modified 3/93 to return if incx .le. 0.
c modified 12/3/93, array(1) declarations changed to array(*)
c
double precision da,dx(*)
integer i,incx,m,mp1,n,nincx
c
if( n.le.0 .or. incx.le.0 )return
if(incx.eq.1)go to 20
c
c code for increment not equal to 1
c
nincx = n*incx
do 10 i = 1,nincx,incx
dx(i) = da*dx(i)
10 continue
return
c
c code for increment equal to 1
c
c
c clean-up loop
c
20 m = mod(n,5)
if( m .eq. 0 ) go to 40
do 30 i = 1,m
dx(i) = da*dx(i)
30 continue
if( n .lt. 5 ) return
40 mp1 = m + 1
do 50 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)
50 continue
return
end
integer function idamax(n,dx,incx)
c
c finds the index of element having max. absolute value.
c jack dongarra, linpack, 3/11/78.
c modified 3/93 to return if incx .le. 0.
c modified 12/3/93, array(1) declarations changed to array(*)
c
double precision dx(*),dmax
integer i,incx,ix,n
c
idamax = 0
if( n.lt.1 .or. incx.le.0 ) return
idamax = 1
if(n.eq.1)return
if(incx.eq.1)go to 20
c
c code for increment not equal to 1
c
ix = 1
dmax = dabs(dx(1))
ix = ix + incx
do 10 i = 2,n
if(dabs(dx(ix)).le.dmax) go to 5
idamax = i
dmax = dabs(dx(ix))
5 ix = ix + incx
10 continue
return
c
c code for increment equal to 1
c
20 dmax = dabs(dx(1))
do 30 i = 2,n
if(dabs(dx(i)).le.dmax) go to 30
idamax = i
dmax = dabs(dx(i))
30 continue
return
end
subroutine dswap (n,dx,incx,dy,incy)
c interchanges two vectors.
c uses unrolled loops for increments equal one.
c jack dongarra, linpack, 3/11/78.
c modified 12/3/93, array(1) declarations changed to array(*)
c
double precision dx(*),dy(*),dtemp
integer i,incx,incy,ix,iy,m,mp1,n
c
if(n.le.0)return
if(incx.eq.1.and.incy.eq.1)go to 20
c
c code for unequal increments or equal increments not equal
c to 1
c
ix = 1
iy = 1
if(incx.lt.0)ix = (-n+1)*incx + 1
if(incy.lt.0)iy = (-n+1)*incy + 1
do 10 i = 1,n
dtemp = dx(ix)
dx(ix) = dy(iy)
dy(iy) = dtemp
ix = ix + incx
iy = iy + incy
10 continue
return
c
c code for both increments equal to 1
c
c
c clean-up loop
c
20 m = mod(n,3)
if( m .eq. 0 ) go to 40
do 30 i = 1,m
dtemp = dx(i)
dx(i) = dy(i)
dy(i) = dtemp
30 continue
if( n .lt. 3 ) return
40 mp1 = m + 1
do 50 i = mp1,n,3
dtemp = dx(i)
dx(i) = dy(i)
dy(i) = dtemp
dtemp = dx(i + 1)
dx(i + 1) = dy(i + 1)
dy(i + 1) = dtemp
dtemp = dx(i + 2)
dx(i + 2) = dy(i + 2)
dy(i + 2) = dtemp
50 continue
return
end
double precision function dasum(n,dx,incx)
c
c takes the sum of the absolute values.
c jack dongarra, linpack, 3/11/78.
c modified 3/93 to return if incx .le. 0.
c modified 12/3/93, array(1) declarations changed to array(*)
c
double precision dx(*),dtemp
integer i,incx,m,mp1,n,nincx
c
dasum = 0.0d0
dtemp = 0.0d0
if( n.le.0 .or. incx.le.0 )return
if(incx.eq.1)go to 20
c
c code for increment not equal to 1
c
nincx = n*incx
do 10 i = 1,nincx,incx
dtemp = dtemp + dabs(dx(i))
10 continue
dasum = dtemp
return
c
c code for increment equal to 1
c
c
c clean-up loop
c
20 m = mod(n,6)
if( m .eq. 0 ) go to 40
do 30 i = 1,m
dtemp = dtemp + dabs(dx(i))
30 continue
if( n .lt. 6 ) go to 60
40 mp1 = m + 1
do 50 i = mp1,n,6
dtemp = dtemp + dabs(dx(i)) + dabs(dx(i + 1)) + dabs(dx(i + 2))
* + dabs(dx(i + 3)) + dabs(dx(i + 4)) + dabs(dx(i + 5))
50 continue
60 dasum = dtemp
return
end
Computing file changes ...
