page@swan.ulowell.edu (Bob Page) (11/03/88)
Submitted-by: strovink%galaxy-43@afit-ab.arpa (Mark A. Strovink)
Posting-number: Volume 2, Issue 46
Archive-name: applications/matlab/src.6
# This is a shell archive.
# Remove everything above and including the cut line.
# Then run the rest of the file through sh.
#----cut here-----cut here-----cut here-----cut here----#
#!/bin/sh
# shar: Shell Archiver
# Run the following text with /bin/sh to create:
# src-6
# This archive created: Wed Nov 2 16:20:42 1988
cat << \SHAR_EOF > src-6
SMALL = 1.D0/2.D0**48
C
C DETERMINE WHAT IS TO BE COMPUTED.
C
WANTU = .FALSE.
WANTV = .FALSE.
JOBU = MOD(JOB,100)/10
NCU = N
IF (JOBU .GT. 1) NCU = MIN0(N,P)
IF (JOBU .NE. 0) WANTU = .TRUE.
IF (MOD(JOB,10) .NE. 0) WANTV = .TRUE.
C
C REDUCE X TO BIDIAGONAL FORM, STORING THE DIAGONAL ELEMENTS
C IN S AND THE SUPER-DIAGONAL ELEMENTS IN E.
C
INFO = 0
NCT = MIN0(N-1,P)
NRT = MAX0(0,MIN0(P-2,N))
LU = MAX0(NCT,NRT)
IF (LU .LT. 1) GO TO 190
DO 180 L = 1, LU
LP1 = L + 1
IF (L .GT. NCT) GO TO 30
C
C COMPUTE THE TRANSFORMATION FOR THE L-TH COLUMN AND
C PLACE THE L-TH DIAGONAL IN S(L).
C
SR(L) = WNRM2(N-L+1,XR(L,L),XI(L,L),1)
SI(L) = 0.0D0
IF (CABS1(SR(L),SI(L)) .EQ. 0.0D0) GO TO 20
IF (CABS1(XR(L,L),XI(L,L)) .EQ. 0.0D0) GO TO 10
CALL WSIGN(SR(L),SI(L),XR(L,L),XI(L,L),SR(L),SI(L))
10 CONTINUE
CALL WDIV(1.0D0,0.0D0,SR(L),SI(L),TR,TI)
CALL WSCAL(N-L+1,TR,TI,XR(L,L),XI(L,L),1)
XR(L,L) = FLOP(1.0D0 + XR(L,L))
20 CONTINUE
SR(L) = -SR(L)
SI(L) = -SI(L)
30 CONTINUE
IF (P .LT. LP1) GO TO 60
DO 50 J = LP1, P
IF (L .GT. NCT) GO TO 40
IF (CABS1(SR(L),SI(L)) .EQ. 0.0D0) GO TO 40
C
C APPLY THE TRANSFORMATION.
C
TR = -WDOTCR(N-L+1,XR(L,L),XI(L,L),1,XR(L,J),XI(L,J),1)
TI = -WDOTCI(N-L+1,XR(L,L),XI(L,L),1,XR(L,J),XI(L,J),1)
CALL WDIV(TR,TI,XR(L,L),XI(L,L),TR,TI)
CALL WAXPY(N-L+1,TR,TI,XR(L,L),XI(L,L),1,XR(L,J),
* XI(L,J),1)
40 CONTINUE
C
C PLACE THE L-TH ROW OF X INTO E FOR THE
C SUBSEQUENT CALCULATION OF THE ROW TRANSFORMATION.
C
ER(J) = XR(L,J)
EI(J) = -XI(L,J)
50 CONTINUE
60 CONTINUE
IF (.NOT.WANTU .OR. L .GT. NCT) GO TO 80
C
C PLACE THE TRANSFORMATION IN U FOR SUBSEQUENT BACK
C MULTIPLICATION.
C
DO 70 I = L, N
UR(I,L) = XR(I,L)
UI(I,L) = XI(I,L)
70 CONTINUE
80 CONTINUE
IF (L .GT. NRT) GO TO 170
C
C COMPUTE THE L-TH ROW TRANSFORMATION AND PLACE THE
C L-TH SUPER-DIAGONAL IN E(L).
C
ER(L) = WNRM2(P-L,ER(LP1),EI(LP1),1)
EI(L) = 0.0D0
IF (CABS1(ER(L),EI(L)) .EQ. 0.0D0) GO TO 100
IF (CABS1(ER(LP1),EI(LP1)) .EQ. 0.0D0) GO TO 90
CALL WSIGN(ER(L),EI(L),ER(LP1),EI(LP1),ER(L),EI(L))
90 CONTINUE
CALL WDIV(1.0D0,0.0D0,ER(L),EI(L),TR,TI)
CALL WSCAL(P-L,TR,TI,ER(LP1),EI(LP1),1)
ER(LP1) = FLOP(1.0D0 + ER(LP1))
100 CONTINUE
ER(L) = -ER(L)
EI(L) = +EI(L)
IF (LP1 .GT. N .OR. CABS1(ER(L),EI(L)) .EQ. 0.0D0)
* GO TO 140
C
C APPLY THE TRANSFORMATION.
C
DO 110 I = LP1, N
WORKR(I) = 0.0D0
WORKI(I) = 0.0D0
110 CONTINUE
DO 120 J = LP1, P
CALL WAXPY(N-L,ER(J),EI(J),XR(LP1,J),XI(LP1,J),1,
* WORKR(LP1),WORKI(LP1),1)
120 CONTINUE
DO 130 J = LP1, P
CALL WDIV(-ER(J),-EI(J),ER(LP1),EI(LP1),TR,TI)
CALL WAXPY(N-L,TR,-TI,WORKR(LP1),WORKI(LP1),1,
* XR(LP1,J),XI(LP1,J),1)
130 CONTINUE
140 CONTINUE
IF (.NOT.WANTV) GO TO 160
C
C PLACE THE TRANSFORMATION IN V FOR SUBSEQUENT
C BACK MULTIPLICATION.
C
DO 150 I = LP1, P
VR(I,L) = ER(I)
VI(I,L) = EI(I)
150 CONTINUE
160 CONTINUE
170 CONTINUE
180 CONTINUE
190 CONTINUE
C
C SET UP THE FINAL BIDIAGONAL MATRIX OR ORDER M.
C
M = MIN0(P,N+1)
NCTP1 = NCT + 1
NRTP1 = NRT + 1
IF (NCT .GE. P) GO TO 200
SR(NCTP1) = XR(NCTP1,NCTP1)
SI(NCTP1) = XI(NCTP1,NCTP1)
200 CONTINUE
IF (N .GE. M) GO TO 210
SR(M) = 0.0D0
SI(M) = 0.0D0
210 CONTINUE
IF (NRTP1 .GE. M) GO TO 220
ER(NRTP1) = XR(NRTP1,M)
EI(NRTP1) = XI(NRTP1,M)
220 CONTINUE
ER(M) = 0.0D0
EI(M) = 0.0D0
C
C IF REQUIRED, GENERATE U.
C
IF (.NOT.WANTU) GO TO 350
IF (NCU .LT. NCTP1) GO TO 250
DO 240 J = NCTP1, NCU
DO 230 I = 1, N
UR(I,J) = 0.0D0
UI(I,J) = 0.0D0
230 CONTINUE
UR(J,J) = 1.0D0
UI(J,J) = 0.0D0
240 CONTINUE
250 CONTINUE
IF (NCT .LT. 1) GO TO 340
DO 330 LL = 1, NCT
L = NCT - LL + 1
IF (CABS1(SR(L),SI(L)) .EQ. 0.0D0) GO TO 300
LP1 = L + 1
IF (NCU .LT. LP1) GO TO 270
DO 260 J = LP1, NCU
TR = -WDOTCR(N-L+1,UR(L,L),UI(L,L),1,UR(L,J),
* UI(L,J),1)
TI = -WDOTCI(N-L+1,UR(L,L),UI(L,L),1,UR(L,J),
* UI(L,J),1)
CALL WDIV(TR,TI,UR(L,L),UI(L,L),TR,TI)
CALL WAXPY(N-L+1,TR,TI,UR(L,L),UI(L,L),1,UR(L,J),
* UI(L,J),1)
260 CONTINUE
270 CONTINUE
CALL WRSCAL(N-L+1,-1.0D0,UR(L,L),UI(L,L),1)
UR(L,L) = FLOP(1.0D0 + UR(L,L))
LM1 = L - 1
IF (LM1 .LT. 1) GO TO 290
DO 280 I = 1, LM1
UR(I,L) = 0.0D0
UI(I,L) = 0.0D0
280 CONTINUE
290 CONTINUE
GO TO 320
300 CONTINUE
DO 310 I = 1, N
UR(I,L) = 0.0D0
UI(I,L) = 0.0D0
310 CONTINUE
UR(L,L) = 1.0D0
UI(L,L) = 0.0D0
320 CONTINUE
330 CONTINUE
340 CONTINUE
350 CONTINUE
C
C IF IT IS REQUIRED, GENERATE V.
C
IF (.NOT.WANTV) GO TO 400
DO 390 LL = 1, P
L = P - LL + 1
LP1 = L + 1
IF (L .GT. NRT) GO TO 370
IF (CABS1(ER(L),EI(L)) .EQ. 0.0D0) GO TO 370
DO 360 J = LP1, P
TR = -WDOTCR(P-L,VR(LP1,L),VI(LP1,L),1,VR(LP1,J),
* VI(LP1,J),1)
TI = -WDOTCI(P-L,VR(LP1,L),VI(LP1,L),1,VR(LP1,J),
* VI(LP1,J),1)
CALL WDIV(TR,TI,VR(LP1,L),VI(LP1,L),TR,TI)
CALL WAXPY(P-L,TR,TI,VR(LP1,L),VI(LP1,L),1,VR(LP1,J),
* VI(LP1,J),1)
360 CONTINUE
370 CONTINUE
DO 380 I = 1, P
VR(I,L) = 0.0D0
VI(I,L) = 0.0D0
380 CONTINUE
VR(L,L) = 1.0D0
VI(L,L) = 0.0D0
390 CONTINUE
400 CONTINUE
C
C TRANSFORM S AND E SO THAT THEY ARE REAL.
C
DO 420 I = 1, M
TR = PYTHAG(SR(I),SI(I))
IF (TR .EQ. 0.0D0) GO TO 405
RR = SR(I)/TR
RI = SI(I)/TR
SR(I) = TR
SI(I) = 0.0D0
IF (I .LT. M) CALL WDIV(ER(I),EI(I),RR,RI,ER(I),EI(I))
IF (WANTU) CALL WSCAL(N,RR,RI,UR(1,I),UI(1,I),1)
405 CONTINUE
C ...EXIT
IF (I .EQ. M) GO TO 430
TR = PYTHAG(ER(I),EI(I))
IF (TR .EQ. 0.0D0) GO TO 410
CALL WDIV(TR,0.0D0,ER(I),EI(I),RR,RI)
ER(I) = TR
EI(I) = 0.0D0
CALL WMUL(SR(I+1),SI(I+1),RR,RI,SR(I+1),SI(I+1))
IF (WANTV) CALL WSCAL(P,RR,RI,VR(1,I+1),VI(1,I+1),1)
410 CONTINUE
420 CONTINUE
430 CONTINUE
C
C MAIN ITERATION LOOP FOR THE SINGULAR VALUES.
C
MM = M
ITER = 0
440 CONTINUE
C
C QUIT IF ALL THE SINGULAR VALUES HAVE BEEN FOUND.
C
C ...EXIT
IF (M .EQ. 0) GO TO 700
C
C IF TOO MANY ITERATIONS HAVE BEEN PERFORMED, SET
C FLAG AND RETURN.
C
IF (ITER .LT. MAXIT) GO TO 450
INFO = M
C ......EXIT
GO TO 700
450 CONTINUE
C
C THIS SECTION OF THE PROGRAM INSPECTS FOR
C NEGLIGIBLE ELEMENTS IN THE S AND E ARRAYS. ON
C COMPLETION THE VARIABLE KASE IS SET AS FOLLOWS.
C
C KASE = 1 IF SR(M) AND ER(L-1) ARE NEGLIGIBLE AND L.LT.M
C KASE = 2 IF SR(L) IS NEGLIGIBLE AND L.LT.M
C KASE = 3 IF ER(L-1) IS NEGLIGIBLE, L.LT.M, AND
C SR(L), ..., SR(M) ARE NOT NEGLIGIBLE (QR STEP).
C KASE = 4 IF ER(M-1) IS NEGLIGIBLE (CONVERGENCE).
C
DO 470 LL = 1, M
L = M - LL
C ...EXIT
IF (L .EQ. 0) GO TO 480
TEST = FLOP(DABS(SR(L)) + DABS(SR(L+1)))
ZTEST = FLOP(TEST + DABS(ER(L))/2.0D0)
IF (SMALL*ZTEST .NE. SMALL*TEST) GO TO 460
ER(L) = 0.0D0
C ......EXIT
GO TO 480
460 CONTINUE
470 CONTINUE
480 CONTINUE
IF (L .NE. M - 1) GO TO 490
KASE = 4
GO TO 560
490 CONTINUE
LP1 = L + 1
MP1 = M + 1
DO 510 LLS = LP1, MP1
LS = M - LLS + LP1
C ...EXIT
IF (LS .EQ. L) GO TO 520
TEST = 0.0D0
IF (LS .NE. M) TEST = FLOP(TEST + DABS(ER(LS)))
IF (LS .NE. L + 1) TEST = FLOP(TEST + DABS(ER(LS-1)))
ZTEST = FLOP(TEST + DABS(SR(LS))/2.0D0)
IF (SMALL*ZTEST .NE. SMALL*TEST) GO TO 500
SR(LS) = 0.0D0
C ......EXIT
GO TO 520
500 CONTINUE
510 CONTINUE
520 CONTINUE
IF (LS .NE. L) GO TO 530
KASE = 3
GO TO 550
530 CONTINUE
IF (LS .NE. M) GO TO 540
KASE = 1
GO TO 550
540 CONTINUE
KASE = 2
L = LS
550 CONTINUE
560 CONTINUE
L = L + 1
C
C PERFORM THE TASK INDICATED BY KASE.
C
GO TO (570, 600, 620, 650), KASE
C
C DEFLATE NEGLIGIBLE SR(M).
C
570 CONTINUE
MM1 = M - 1
F = ER(M-1)
ER(M-1) = 0.0D0
DO 590 KK = L, MM1
K = MM1 - KK + L
T1 = SR(K)
CALL RROTG(T1,F,CS,SN)
SR(K) = T1
IF (K .EQ. L) GO TO 580
F = FLOP(-SN*ER(K-1))
ER(K-1) = FLOP(CS*ER(K-1))
580 CONTINUE
IF (WANTV) CALL RROT(P,VR(1,K),1,VR(1,M),1,CS,SN)
IF (WANTV) CALL RROT(P,VI(1,K),1,VI(1,M),1,CS,SN)
590 CONTINUE
GO TO 690
C
C SPLIT AT NEGLIGIBLE SR(L).
C
600 CONTINUE
F = ER(L-1)
ER(L-1) = 0.0D0
DO 610 K = L, M
T1 = SR(K)
CALL RROTG(T1,F,CS,SN)
SR(K) = T1
F = FLOP(-SN*ER(K))
ER(K) = FLOP(CS*ER(K))
IF (WANTU) CALL RROT(N,UR(1,K),1,UR(1,L-1),1,CS,SN)
IF (WANTU) CALL RROT(N,UI(1,K),1,UI(1,L-1),1,CS,SN)
610 CONTINUE
GO TO 690
C
C PERFORM ONE QR STEP.
C
620 CONTINUE
C
C CALCULATE THE SHIFT.
C
SCALE = DMAX1(DABS(SR(M)),DABS(SR(M-1)),DABS(ER(M-1)),
* DABS(SR(L)),DABS(ER(L)))
SM = SR(M)/SCALE
SMM1 = SR(M-1)/SCALE
EMM1 = ER(M-1)/SCALE
SL = SR(L)/SCALE
EL = ER(L)/SCALE
B = FLOP(((SMM1 + SM)*(SMM1 - SM) + EMM1**2)/2.0D0)
C = FLOP((SM*EMM1)**2)
SHIFT = 0.0D0
IF (B .EQ. 0.0D0 .AND. C .EQ. 0.0D0) GO TO 630
SHIFT = FLOP(DSQRT(B**2+C))
IF (B .LT. 0.0D0) SHIFT = -SHIFT
SHIFT = FLOP(C/(B + SHIFT))
630 CONTINUE
F = FLOP((SL + SM)*(SL - SM) - SHIFT)
G = FLOP(SL*EL)
C
C CHASE ZEROS.
C
MM1 = M - 1
DO 640 K = L, MM1
CALL RROTG(F,G,CS,SN)
IF (K .NE. L) ER(K-1) = F
F = FLOP(CS*SR(K) + SN*ER(K))
ER(K) = FLOP(CS*ER(K) - SN*SR(K))
G = FLOP(SN*SR(K+1))
SR(K+1) = FLOP(CS*SR(K+1))
IF (WANTV) CALL RROT(P,VR(1,K),1,VR(1,K+1),1,CS,SN)
IF (WANTV) CALL RROT(P,VI(1,K),1,VI(1,K+1),1,CS,SN)
CALL RROTG(F,G,CS,SN)
SR(K) = F
F = FLOP(CS*ER(K) + SN*SR(K+1))
SR(K+1) = FLOP(-SN*ER(K) + CS*SR(K+1))
G = FLOP(SN*ER(K+1))
ER(K+1) = FLOP(CS*ER(K+1))
IF (WANTU .AND. K .LT. N)
* CALL RROT(N,UR(1,K),1,UR(1,K+1),1,CS,SN)
IF (WANTU .AND. K .LT. N)
* CALL RROT(N,UI(1,K),1,UI(1,K+1),1,CS,SN)
640 CONTINUE
ER(M-1) = F
ITER = ITER + 1
GO TO 690
C
C CONVERGENCE
C
650 CONTINUE
C
C MAKE THE SINGULAR VALUE POSITIVE
C
IF (SR(L) .GE. 0.0D0) GO TO 660
SR(L) = -SR(L)
IF (WANTV) CALL WRSCAL(P,-1.0D0,VR(1,L),VI(1,L),1)
660 CONTINUE
C
C ORDER THE SINGULAR VALUE.
C
670 IF (L .EQ. MM) GO TO 680
C ...EXIT
IF (SR(L) .GE. SR(L+1)) GO TO 680
TR = SR(L)
SR(L) = SR(L+1)
SR(L+1) = TR
IF (WANTV .AND. L .LT. P)
* CALL WSWAP(P,VR(1,L),VI(1,L),1,VR(1,L+1),VI(1,L+1),1)
IF (WANTU .AND. L .LT. N)
* CALL WSWAP(N,UR(1,L),UI(1,L),1,UR(1,L+1),UI(1,L+1),1)
L = L + 1
GO TO 670
680 CONTINUE
ITER = 0
M = M - 1
690 CONTINUE
GO TO 440
700 CONTINUE
RETURN
END
SUBROUTINE WQRDC(XR,XI,LDX,N,P,QRAUXR,QRAUXI,JPVT,WORKR,WORKI,
* JOB)
INTEGER LDX,N,P,JOB
INTEGER JPVT(1)
DOUBLE PRECISION XR(LDX,1),XI(LDX,1),QRAUXR(1),QRAUXI(1),
* WORKR(1),WORKI(1)
C
C WQRDC USES HOUSEHOLDER TRANSFORMATIONS TO COMPUTE THE QR
C FACTORIZATION OF AN N BY P MATRIX X. COLUMN PIVOTING
C BASED ON THE 2-NORMS OF THE REDUCED COLUMNS MAY BE
C PERFORMED AT THE USERS OPTION.
C
C ON ENTRY
C
C X DOUBLE-COMPLEX(LDX,P), WHERE LDX .GE. N.
C X CONTAINS THE MATRIX WHOSE DECOMPOSITION IS TO BE
C COMPUTED.
C
C LDX INTEGER.
C LDX IS THE LEADING DIMENSION OF THE ARRAY X.
C
C N INTEGER.
C N IS THE NUMBER OF ROWS OF THE MATRIX X.
C
C P INTEGER.
C P IS THE NUMBER OF COLUMNS OF THE MATRIX X.
C
C JPVT INTEGER(P).
C JPVT CONTAINS INTEGERS THAT CONTROL THE SELECTION
C OF THE PIVOT COLUMNS. THE K-TH COLUMN X(K) OF X
C IS PLACED IN ONE OF THREE CLASSES ACCORDING TO THE
C VALUE OF JPVT(K).
C
C IF JPVT(K) .GT. 0, THEN X(K) IS AN INITIAL
C COLUMN.
C
C IF JPVT(K) .EQ. 0, THEN X(K) IS A FREE COLUMN.
C
C IF JPVT(K) .LT. 0, THEN X(K) IS A FINAL COLUMN.
C
C BEFORE THE DECOMPOSITION IS COMPUTED, INITIAL COLUMNS
C ARE MOVED TO THE BEGINNING OF THE ARRAY X AND FINAL
C COLUMNS TO THE END. BOTH INITIAL AND FINAL COLUMNS
C ARE FROZEN IN PLACE DURING THE COMPUTATION AND ONLY
C FREE COLUMNS ARE MOVED. AT THE K-TH STAGE OF THE
C REDUCTION, IF X(K) IS OCCUPIED BY A FREE COLUMN
C IT IS INTERCHANGED WITH THE FREE COLUMN OF LARGEST
C REDUCED NORM. JPVT IS NOT REFERENCED IF
C JOB .EQ. 0.
C
C WORK DOUBLE-COMPLEX(P).
C WORK IS A WORK ARRAY. WORK IS NOT REFERENCED IF
C JOB .EQ. 0.
C
C JOB INTEGER.
C JOB IS AN INTEGER THAT INITIATES COLUMN PIVOTING.
C IF JOB .EQ. 0, NO PIVOTING IS DONE.
C IF JOB .NE. 0, PIVOTING IS DONE.
C
C ON RETURN
C
C X X CONTAINS IN ITS UPPER TRIANGLE THE UPPER
C TRIANGULAR MATRIX R OF THE QR FACTORIZATION.
C BELOW ITS DIAGONAL X CONTAINS INFORMATION FROM
C WHICH THE UNITARY PART OF THE DECOMPOSITION
C CAN BE RECOVERED. NOTE THAT IF PIVOTING HAS
C BEEN REQUESTED, THE DECOMPOSITION IS NOT THAT
C OF THE ORIGINAL MATRIX X BUT THAT OF X
C WITH ITS COLUMNS PERMUTED AS DESCRIBED BY JPVT.
C
C QRAUX DOUBLE-COMPLEX(P).
C QRAUX CONTAINS FURTHER INFORMATION REQUIRED TO RECOVER
C THE UNITARY PART OF THE DECOMPOSITION.
C
C JPVT JPVT(K) CONTAINS THE INDEX OF THE COLUMN OF THE
C ORIGINAL MATRIX THAT HAS BEEN INTERCHANGED INTO
C THE K-TH COLUMN, IF PIVOTING WAS REQUESTED.
C
C LINPACK. THIS VERSION DATED 07/03/79 .
C G.W. STEWART, UNIVERSITY OF MARYLAND, ARGONNE NATIONAL LAB.
C
C WQRDC USES THE FOLLOWING FUNCTIONS AND SUBPROGRAMS.
C
C BLAS WAXPY,PYTHAG,WDOTCR,WDOTCI,WSCAL,WSWAP,WNRM2
C FORTRAN DABS,DIMAG,DMAX1,MIN0
C
C INTERNAL VARIABLES
C
INTEGER J,JP,L,LP1,LUP,MAXJ,PL,PU
DOUBLE PRECISION MAXNRM,WNRM2,TT
DOUBLE PRECISION PYTHAG,WDOTCR,WDOTCI,NRMXLR,NRMXLI,TR,TI,FLOP
LOGICAL NEGJ,SWAPJ
C
DOUBLE PRECISION ZDUMR,ZDUMI
DOUBLE PRECISION CABS1
CABS1(ZDUMR,ZDUMI) = DABS(ZDUMR) + DABS(ZDUMI)
C
PL = 1
PU = 0
IF (JOB .EQ. 0) GO TO 60
C
C PIVOTING HAS BEEN REQUESTED. REARRANGE THE COLUMNS
C ACCORDING TO JPVT.
C
DO 20 J = 1, P
SWAPJ = JPVT(J) .GT. 0
NEGJ = JPVT(J) .LT. 0
JPVT(J) = J
IF (NEGJ) JPVT(J) = -J
IF (.NOT.SWAPJ) GO TO 10
IF (J .NE. PL)
* CALL WSWAP(N,XR(1,PL),XI(1,PL),1,XR(1,J),XI(1,J),1)
JPVT(J) = JPVT(PL)
JPVT(PL) = J
PL = PL + 1
10 CONTINUE
20 CONTINUE
PU = P
DO 50 JJ = 1, P
J = P - JJ + 1
IF (JPVT(J) .GE. 0) GO TO 40
JPVT(J) = -JPVT(J)
IF (J .EQ. PU) GO TO 30
CALL WSWAP(N,XR(1,PU),XI(1,PU),1,XR(1,J),XI(1,J),1)
JP = JPVT(PU)
JPVT(PU) = JPVT(J)
JPVT(J) = JP
30 CONTINUE
PU = PU - 1
40 CONTINUE
50 CONTINUE
60 CONTINUE
C
C COMPUTE THE NORMS OF THE FREE COLUMNS.
C
IF (PU .LT. PL) GO TO 80
DO 70 J = PL, PU
QRAUXR(J) = WNRM2(N,XR(1,J),XI(1,J),1)
QRAUXI(J) = 0.0D0
WORKR(J) = QRAUXR(J)
WORKI(J) = QRAUXI(J)
70 CONTINUE
80 CONTINUE
C
C PERFORM THE HOUSEHOLDER REDUCTION OF X.
C
LUP = MIN0(N,P)
DO 210 L = 1, LUP
IF (L .LT. PL .OR. L .GE. PU) GO TO 120
C
C LOCATE THE COLUMN OF LARGEST NORM AND BRING IT
C INTO THE PIVOT POSITION.
C
MAXNRM = 0.0D0
MAXJ = L
DO 100 J = L, PU
IF (QRAUXR(J) .LE. MAXNRM) GO TO 90
MAXNRM = QRAUXR(J)
MAXJ = J
90 CONTINUE
100 CONTINUE
IF (MAXJ .EQ. L) GO TO 110
CALL WSWAP(N,XR(1,L),XI(1,L),1,XR(1,MAXJ),XI(1,MAXJ),1)
QRAUXR(MAXJ) = QRAUXR(L)
QRAUXI(MAXJ) = QRAUXI(L)
WORKR(MAXJ) = WORKR(L)
WORKI(MAXJ) = WORKI(L)
JP = JPVT(MAXJ)
JPVT(MAXJ) = JPVT(L)
JPVT(L) = JP
110 CONTINUE
120 CONTINUE
QRAUXR(L) = 0.0D0
QRAUXI(L) = 0.0D0
IF (L .EQ. N) GO TO 200
C
C COMPUTE THE HOUSEHOLDER TRANSFORMATION FOR COLUMN L.
C
NRMXLR = WNRM2(N-L+1,XR(L,L),XI(L,L),1)
NRMXLI = 0.0D0
IF (CABS1(NRMXLR,NRMXLI) .EQ. 0.0D0) GO TO 190
IF (CABS1(XR(L,L),XI(L,L)) .EQ. 0.0D0) GO TO 130
CALL WSIGN(NRMXLR,NRMXLI,XR(L,L),XI(L,L),NRMXLR,NRMXLI)
130 CONTINUE
CALL WDIV(1.0D0,0.0D0,NRMXLR,NRMXLI,TR,TI)
CALL WSCAL(N-L+1,TR,TI,XR(L,L),XI(L,L),1)
XR(L,L) = FLOP(1.0D0 + XR(L,L))
C
C APPLY THE TRANSFORMATION TO THE REMAINING COLUMNS,
C UPDATING THE NORMS.
C
LP1 = L + 1
IF (P .LT. LP1) GO TO 180
DO 170 J = LP1, P
TR = -WDOTCR(N-L+1,XR(L,L),XI(L,L),1,XR(L,J),
* XI(L,J),1)
TI = -WDOTCI(N-L+1,XR(L,L),XI(L,L),1,XR(L,J),
* XI(L,J),1)
CALL WDIV(TR,TI,XR(L,L),XI(L,L),TR,TI)
CALL WAXPY(N-L+1,TR,TI,XR(L,L),XI(L,L),1,XR(L,J),
* XI(L,J),1)
IF (J .LT. PL .OR. J .GT. PU) GO TO 160
IF (CABS1(QRAUXR(J),QRAUXI(J)) .EQ. 0.0D0)
* GO TO 160
TT = 1.0D0 - (PYTHAG(XR(L,J),XI(L,J))/QRAUXR(J))**2
TT = DMAX1(TT,0.0D0)
TR = FLOP(TT)
TT = FLOP(1.0D0+0.05D0*TT*(QRAUXR(J)/WORKR(J))**2)
IF (TT .EQ. 1.0D0) GO TO 140
QRAUXR(J) = QRAUXR(J)*DSQRT(TR)
QRAUXI(J) = QRAUXI(J)*DSQRT(TR)
GO TO 150
140 CONTINUE
QRAUXR(J) = WNRM2(N-L,XR(L+1,J),XI(L+1,J),1)
QRAUXI(J) = 0.0D0
WORKR(J) = QRAUXR(J)
WORKI(J) = QRAUXI(J)
150 CONTINUE
160 CONTINUE
170 CONTINUE
180 CONTINUE
C
C SAVE THE TRANSFORMATION.
C
QRAUXR(L) = XR(L,L)
QRAUXI(L) = XI(L,L)
XR(L,L) = -NRMXLR
XI(L,L) = -NRMXLI
190 CONTINUE
200 CONTINUE
210 CONTINUE
RETURN
END
SUBROUTINE WQRSL(XR,XI,LDX,N,K,QRAUXR,QRAUXI,YR,YI,QYR,QYI,QTYR,
* QTYI,BR,BI,RSDR,RSDI,XBR,XBI,JOB,INFO)
INTEGER LDX,N,K,JOB,INFO
DOUBLE PRECISION XR(LDX,1),XI(LDX,1),QRAUXR(1),QRAUXI(1),YR(1),
* YI(1),QYR(1),QYI(1),QTYR(1),QTYI(1),BR(1),BI(1),
* RSDR(1),RSDI(1),XBR(1),XBI(1)
C
C WQRSL APPLIES THE OUTPUT OF WQRDC TO COMPUTE COORDINATE
C TRANSFORMATIONS, PROJECTIONS, AND LEAST SQUARES SOLUTIONS.
C FOR K .LE. MIN(N,P), LET XK BE THE MATRIX
C
C XK = (X(JPVT(1)),X(JPVT(2)), ... ,X(JPVT(K)))
C
C FORMED FROM COLUMNNS JPVT(1), ... ,JPVT(K) OF THE ORIGINAL
C N X P MATRIX X THAT WAS INPUT TO WQRDC (IF NO PIVOTING WAS
C DONE, XK CONSISTS OF THE FIRST K COLUMNS OF X IN THEIR
C ORIGINAL ORDER). WQRDC PRODUCES A FACTORED UNITARY MATRIX Q
C AND AN UPPER TRIANGULAR MATRIX R SUCH THAT
C
C XK = Q * (R)
C (0)
C
C THIS INFORMATION IS CONTAINED IN CODED FORM IN THE ARRAYS
C X AND QRAUX.
C
C ON ENTRY
C
C X DOUBLE-COMPLEX(LDX,P).
C X CONTAINS THE OUTPUT OF WQRDC.
C
C LDX INTEGER.
C LDX IS THE LEADING DIMENSION OF THE ARRAY X.
C
C N INTEGER.
C N IS THE NUMBER OF ROWS OF THE MATRIX XK. IT MUST
C HAVE THE SAME VALUE AS N IN WQRDC.
C
C K INTEGER.
C K IS THE NUMBER OF COLUMNS OF THE MATRIX XK. K
C MUST NNOT BE GREATER THAN MIN(N,P), WHERE P IS THE
C SAME AS IN THE CALLING SEQUENCE TO WQRDC.
C
C QRAUX DOUBLE-COMPLEX(P).
C QRAUX CONTAINS THE AUXILIARY OUTPUT FROM WQRDC.
C
C Y DOUBLE-COMPLEX(N)
C Y CONTAINS AN N-VECTOR THAT IS TO BE MANIPULATED
C BY WQRSL.
C
C JOB INTEGER.
C JOB SPECIFIES WHAT IS TO BE COMPUTED. JOB HAS
C THE DECIMAL EXPANSION ABCDE, WITH THE FOLLOWING
C MEANING.
C
C IF A.NE.0, COMPUTE QY.
C IF B,C,D, OR E .NE. 0, COMPUTE QTY.
C IF C.NE.0, COMPUTE B.
C IF D.NE.0, COMPUTE RSD.
C IF E.NE.0, COMPUTE XB.
C
C NOTE THAT A REQUEST TO COMPUTE B, RSD, OR XB
C AUTOMATICALLY TRIGGERS THE COMPUTATION OF QTY, FOR
C WHICH AN ARRAY MUST BE PROVIDED IN THE CALLING
C SEQUENCE.
C
C ON RETURN
C
C QY DOUBLE-COMPLEX(N).
C QY CONNTAINS Q*Y, IF ITS COMPUTATION HAS BEEN
C REQUESTED.
C
C QTY DOUBLE-COMPLEX(N).
C QTY CONTAINS CTRANS(Q)*Y, IF ITS COMPUTATION HAS
C BEEN REQUESTED. HERE CTRANS(Q) IS THE CONJUGATE
C TRANSPOSE OF THE MATRIX Q.
C
C B DOUBLE-COMPLEX(K)
C B CONTAINS THE SOLUTION OF THE LEAST SQUARES PROBLEM
C
C MINIMIZE NORM2(Y - XK*B),
C
C IF ITS COMPUTATION HAS BEEN REQUESTED. (NOTE THAT
C IF PIVOTING WAS REQUESTED IN WQRDC, THE J-TH
C COMPONENT OF B WILL BE ASSOCIATED WITH COLUMN JPVT(J)
C OF THE ORIGINAL MATRIX X THAT WAS INPUT INTO WQRDC.)
C
C RSD DOUBLE-COMPLEX(N).
C RSD CONTAINS THE LEAST SQUARES RESIDUAL Y - XK*B,
C IF ITS COMPUTATION HAS BEEN REQUESTED. RSD IS
C ALSO THE ORTHOGONAL PROJECTION OF Y ONTO THE
C ORTHOGONAL COMPLEMENT OF THE COLUMN SPACE OF XK.
C
C XB DOUBLE-COMPLEX(N).
C XB CONTAINS THE LEAST SQUARES APPROXIMATION XK*B,
C IF ITS COMPUTATION HAS BEEN REQUESTED. XB IS ALSO
C THE ORTHOGONAL PROJECTION OF Y ONTO THE COLUMN SPACE
C OF X.
C
C INFO INTEGER.
C INFO IS ZERO UNLESS THE COMPUTATION OF B HAS
C BEEN REQUESTED AND R IS EXACTLY SINGULAR. IN
C THIS CASE, INFO IS THE INDEX OF THE FIRST ZERO
C DIAGONAL ELEMENT OF R AND B IS LEFT UNALTERED.
C
C THE PARAMETERS QY, QTY, B, RSD, AND XB ARE NOT REFERENCED
C IF THEIR COMPUTATION IS NOT REQUESTED AND IN THIS CASE
C CAN BE REPLACED BY DUMMY VARIABLES IN THE CALLING PROGRAM.
C TO SAVE STORAGE, THE USER MAY IN SOME CASES USE THE SAME
C ARRAY FOR DIFFERENT PARAMETERS IN THE CALLING SEQUENCE. A
C FREQUENTLY OCCURING EXAMPLE IS WHEN ONE WISHES TO COMPUTE
C ANY OF B, RSD, OR XB AND DOES NOT NEED Y OR QTY. IN THIS
C CASE ONE MAY IDENTIFY Y, QTY, AND ONE OF B, RSD, OR XB, WHILE
C PROVIDING SEPARATE ARRAYS FOR ANYTHING ELSE THAT IS TO BE
C COMPUTED. THUS THE CALLING SEQUENCE
C
C CALL WQRSL(X,LDX,N,K,QRAUX,Y,DUM,Y,B,Y,DUM,110,INFO)
C
C WILL RESULT IN THE COMPUTATION OF B AND RSD, WITH RSD
C OVERWRITING Y. MORE GENERALLY, EACH ITEM IN THE FOLLOWING
C LIST CONTAINS GROUPS OF PERMISSIBLE IDENTIFICATIONS FOR
C A SINGLE CALLINNG SEQUENCE.
C
C 1. (Y,QTY,B) (RSD) (XB) (QY)
C
C 2. (Y,QTY,RSD) (B) (XB) (QY)
C
C 3. (Y,QTY,XB) (B) (RSD) (QY)
C
C 4. (Y,QY) (QTY,B) (RSD) (XB)
C
C 5. (Y,QY) (QTY,RSD) (B) (XB)
C
C 6. (Y,QY) (QTY,XB) (B) (RSD)
C
C IN ANY GROUP THE VALUE RETURNED IN THE ARRAY ALLOCATED TO
C THE GROUP CORRESPONDS TO THE LAST MEMBER OF THE GROUP.
C
C LINPACK. THIS VERSION DATED 07/03/79 .
C G.W. STEWART, UNIVERSITY OF MARYLAND, ARGONNE NATIONAL LAB.
C
C WQRSL USES THE FOLLOWING FUNCTIONS AND SUBPROGRAMS.
C
C BLAS WAXPY,WCOPY,WDOTCR,WDOTCI
C FORTRAN DABS,DIMAG,MIN0,MOD
C
C INTERNAL VARIABLES
C
INTEGER I,J,JJ,JU,KP1
DOUBLE PRECISION WDOTCR,WDOTCI,TR,TI,TEMPR,TEMPI
LOGICAL CB,CQY,CQTY,CR,CXB
C
DOUBLE PRECISION ZDUMR,ZDUMI
DOUBLE PRECISION CABS1
CABS1(ZDUMR,ZDUMI) = DABS(ZDUMR) + DABS(ZDUMI)
C
C SET INFO FLAG.
C
INFO = 0
C
C DETERMINE WHAT IS TO BE COMPUTED.
C
CQY = JOB/10000 .NE. 0
CQTY = MOD(JOB,10000) .NE. 0
CB = MOD(JOB,1000)/100 .NE. 0
CR = MOD(JOB,100)/10 .NE. 0
CXB = MOD(JOB,10) .NE. 0
JU = MIN0(K,N-1)
C
C SPECIAL ACTION WHEN N=1.
C
IF (JU .NE. 0) GO TO 80
IF (.NOT.CQY) GO TO 10
QYR(1) = YR(1)
QYI(1) = YI(1)
10 CONTINUE
IF (.NOT.CQTY) GO TO 20
QTYR(1) = YR(1)
QTYI(1) = YI(1)
20 CONTINUE
IF (.NOT.CXB) GO TO 30
XBR(1) = YR(1)
XBI(1) = YI(1)
30 CONTINUE
IF (.NOT.CB) GO TO 60
IF (CABS1(XR(1,1),XI(1,1)) .NE. 0.0D0) GO TO 40
INFO = 1
GO TO 50
40 CONTINUE
CALL WDIV(YR(1),YI(1),XR(1,1),XI(1,1),BR(1),BI(1))
50 CONTINUE
60 CONTINUE
IF (.NOT.CR) GO TO 70
RSDR(1) = 0.0D0
RSDI(1) = 0.0D0
70 CONTINUE
GO TO 290
80 CONTINUE
C
C SET UP TO COMPUTE QY OR QTY.
C
IF (CQY) CALL WCOPY(N,YR,YI,1,QYR,QYI,1)
IF (CQTY) CALL WCOPY(N,YR,YI,1,QTYR,QTYI,1)
IF (.NOT.CQY) GO TO 110
C
C COMPUTE QY.
C
DO 100 JJ = 1, JU
J = JU - JJ + 1
IF (CABS1(QRAUXR(J),QRAUXI(J)) .EQ. 0.0D0)
* GO TO 90
TEMPR = XR(J,J)
TEMPI = XI(J,J)
XR(J,J) = QRAUXR(J)
XI(J,J) = QRAUXI(J)
TR = -WDOTCR(N-J+1,XR(J,J),XI(J,J),1,QYR(J),QYI(J),1)
TI = -WDOTCI(N-J+1,XR(J,J),XI(J,J),1,QYR(J),QYI(J),1)
CALL WDIV(TR,TI,XR(J,J),XI(J,J),TR,TI)
CALL WAXPY(N-J+1,TR,TI,XR(J,J),XI(J,J),1,QYR(J),
* QYI(J),1)
XR(J,J) = TEMPR
XI(J,J) = TEMPI
90 CONTINUE
100 CONTINUE
110 CONTINUE
IF (.NOT.CQTY) GO TO 140
C
C COMPUTE CTRANS(Q)*Y.
C
DO 130 J = 1, JU
IF (CABS1(QRAUXR(J),QRAUXI(J)) .EQ. 0.0D0)
* GO TO 120
TEMPR = XR(J,J)
TEMPI = XI(J,J)
XR(J,J) = QRAUXR(J)
XI(J,J) = QRAUXI(J)
TR = -WDOTCR(N-J+1,XR(J,J),XI(J,J),1,QTYR(J),
* QTYI(J),1)
TI = -WDOTCI(N-J+1,XR(J,J),XI(J,J),1,QTYR(J),
* QTYI(J),1)
CALL WDIV(TR,TI,XR(J,J),XI(J,J),TR,TI)
CALL WAXPY(N-J+1,TR,TI,XR(J,J),XI(J,J),1,QTYR(J),
* QTYI(J),1)
XR(J,J) = TEMPR
XI(J,J) = TEMPI
120 CONTINUE
130 CONTINUE
140 CONTINUE
C
C SET UP TO COMPUTE B, RSD, OR XB.
C
IF (CB) CALL WCOPY(K,QTYR,QTYI,1,BR,BI,1)
KP1 = K + 1
IF (CXB) CALL WCOPY(K,QTYR,QTYI,1,XBR,XBI,1)
IF (CR .AND. K .LT. N)
* CALL WCOPY(N-K,QTYR(KP1),QTYI(KP1),1,RSDR(KP1),RSDI(KP1),1)
IF (.NOT.CXB .OR. KP1 .GT. N) GO TO 160
DO 150 I = KP1, N
XBR(I) = 0.0D0
XBI(I) = 0.0D0
150 CONTINUE
160 CONTINUE
IF (.NOT.CR) GO TO 180
DO 170 I = 1, K
RSDR(I) = 0.0D0
RSDI(I) = 0.0D0
170 CONTINUE
180 CONTINUE
IF (.NOT.CB) GO TO 230
C
C COMPUTE B.
C
DO 210 JJ = 1, K
J = K - JJ + 1
IF (CABS1(XR(J,J),XI(J,J)) .NE. 0.0D0) GO TO 190
INFO = J
C ......EXIT
C ......EXIT
GO TO 220
190 CONTINUE
CALL WDIV(BR(J),BI(J),XR(J,J),XI(J,J),BR(J),BI(J))
IF (J .EQ. 1) GO TO 200
TR = -BR(J)
TI = -BI(J)
CALL WAXPY(J-1,TR,TI,XR(1,J),XI(1,J),1,BR,BI,1)
200 CONTINUE
210 CONTINUE
220 CONTINUE
230 CONTINUE
IF (.NOT.CR .AND. .NOT.CXB) GO TO 280
C
C COMPUTE RSD OR XB AS REQUIRED.
C
DO 270 JJ = 1, JU
J = JU - JJ + 1
IF (CABS1(QRAUXR(J),QRAUXI(J)) .EQ. 0.0D0)
* GO TO 260
TEMPR = XR(J,J)
TEMPI = XI(J,J)
XR(J,J) = QRAUXR(J)
XI(J,J) = QRAUXI(J)
IF (.NOT.CR) GO TO 240
TR = -WDOTCR(N-J+1,XR(J,J),XI(J,J),1,RSDR(J),
* RSDI(J),1)
TI = -WDOTCI(N-J+1,XR(J,J),XI(J,J),1,RSDR(J),
* RSDI(J),1)
CALL WDIV(TR,TI,XR(J,J),XI(J,J),TR,TI)
CALL WAXPY(N-J+1,TR,TI,XR(J,J),XI(J,J),1,RSDR(J),
* RSDI(J),1)
240 CONTINUE
IF (.NOT.CXB) GO TO 250
TR = -WDOTCR(N-J+1,XR(J,J),XI(J,J),1,XBR(J),
* XBI(J),1)
TI = -WDOTCI(N-J+1,XR(J,J),XI(J,J),1,XBR(J),
* XBI(J),1)
CALL WDIV(TR,TI,XR(J,J),XI(J,J),TR,TI)
CALL WAXPY(N-J+1,TR,TI,XR(J,J),XI(J,J),1,XBR(J),
* XBI(J),1)
250 CONTINUE
XR(J,J) = TEMPR
XI(J,J) = TEMPI
260 CONTINUE
270 CONTINUE
280 CONTINUE
290 CONTINUE
RETURN
END
SUBROUTINE MAGIC(A,LDA,N)
C
C ALGORITHMS FOR MAGIC SQUARES TAKEN FROM
C MATHEMATICAL RECREATIONS AND ESSAYS, 12TH ED.,
C BY W. W. ROUSE BALL AND H. S. M. COXETER
C
DOUBLE PRECISION A(LDA,N),T
C
IF (MOD(N,4) .EQ. 0) GO TO 100
IF (MOD(N,2) .EQ. 0) M = N/2
IF (MOD(N,2) .NE. 0) M = N
C
C ODD ORDER OR UPPER CORNER OF EVEN ORDER
C
DO 20 J = 1,M
DO 10 I = 1,M
A(I,J) = 0
10 CONTINUE
20 CONTINUE
I = 1
J = (M+1)/2
MM = M*M
DO 40 K = 1, MM
A(I,J) = K
I1 = I-1
J1 = J+1
IF(I1.LT.1) I1 = M
IF(J1.GT.M) J1 = 1
IF(IDINT(A(I1,J1)).EQ.0) GO TO 30
I1 = I+1
J1 = J
30 I = I1
J = J1
40 CONTINUE
IF (MOD(N,2) .NE. 0) RETURN
C
C REST OF EVEN ORDER
C
T = M*M
DO 60 I = 1, M
DO 50 J = 1, M
IM = I+M
JM = J+M
A(I,JM) = A(I,J) + 2*T
A(IM,J) = A(I,J) + 3*T
A(IM,JM) = A(I,J) + T
50 CONTINUE
60 CONTINUE
M1 = (M-1)/2
IF (M1.EQ.0) RETURN
DO 70 J = 1, M1
CALL RSWAP(M,A(1,J),1,A(M+1,J),1)
70 CONTINUE
M1 = (M+1)/2
M2 = M1 + M
CALL RSWAP(1,A(M1,1),1,A(M2,1),1)
CALL RSWAP(1,A(M1,M1),1,A(M2,M1),1)
M1 = N+1-(M-3)/2
IF(M1.GT.N) RETURN
DO 80 J = M1, N
CALL RSWAP(M,A(1,J),1,A(M+1,J),1)
80 CONTINUE
RETURN
C
C DOUBLE EVEN ORDER
C
100 K = 1
DO 120 I = 1, N
DO 110 J = 1, N
A(I,J) = K
IF (MOD(I,4)/2 .EQ. MOD(J,4)/2) A(I,J) = N*N+1 - K
K = K+1
110 CONTINUE
120 CONTINUE
RETURN
END
SUBROUTINE BASE(X,B,EPS,S,N)
DOUBLE PRECISION X,B,EPS,S(1),T
C
C STORE BASE B REPRESENTATION OF X IN S(1:N)
C
INTEGER PLUS,MINUS,DOT,ZERO,COMMA
DATA PLUS/41/,MINUS/42/,DOT/47/,ZERO/0/,COMMA/48/
L = 1
IF (X .GE. 0.0D0) S(L) = PLUS
IF (X .LT. 0.0D0) S(L) = MINUS
S(L+1) = ZERO
S(L+2) = DOT
X = DABS(X)
IF (X .NE. 0.0D0) K = DLOG(X)/DLOG(B)
IF (X .EQ. 0.0D0) K = 0
IF (X .GT. 1.0D0) K = K + 1
X = X/B**K
IF (B*X .GE. B) K = K + 1
IF (B*X .GE. B) X = X/B
IF (EPS .NE. 0.0D0) M = -DLOG(EPS)/DLOG(B) + 4
IF (EPS .EQ. 0.0D0) M = 54
DO 10 L = 4, M
X = B*X
J = IDINT(X)
S(L) = DFLOAT(J)
X = X - S(L)
10 CONTINUE
S(M+1) = COMMA
IF (K .GE. 0) S(M+2) = PLUS
IF (K .LT. 0) S(M+2) = MINUS
T = DABS(DFLOAT(K))
N = M + 3
IF (T .GE. B) N = N + IDINT(DLOG(T)/DLOG(B))
L = N
20 J = IDINT(DMOD(T,B))
S(L) = DFLOAT(J)
L = L - 1
T = T/B
IF (L .GE. M+3) GO TO 20
RETURN
END
DOUBLE PRECISION FUNCTION URAND(IY)
INTEGER IY
C
C URAND IS A UNIFORM RANDOM NUMBER GENERATOR BASED ON THEORY AND
C SUGGESTIONS GIVEN IN D.E. KNUTH (1969), VOL 2. THE INTEGER IY
C SHOULD BE INITIALIZED TO AN ARBITRARY INTEGER PRIOR TO THE FIRST CALL
C TO URAND. THE CALLING PROGRAM SHOULD NOT ALTER THE VALUE OF IY
C BETWEEN SUBSEQUENT CALLS TO URAND. VALUES OF URAND WILL BE RETURNED
C IN THE INTERVAL (0,1).
C
INTEGER IA,IC,ITWO,M2,M,MIC
DOUBLE PRECISION HALFM,S
DOUBLE PRECISION DATAN,DSQRT
DATA M2/0/,ITWO/2/
IF (M2 .NE. 0) GO TO 20
C
C IF FIRST ENTRY, COMPUTE MACHINE INTEGER WORD LENGTH
C
M = 1
10 M2 = M
M = ITWO*M2
IF (M .GT. M2) GO TO 10
HALFM = M2
C
C COMPUTE MULTIPLIER AND INCREMENT FOR LINEAR CONGRUENTIAL METHOD
C
IA = 8*IDINT(HALFM*DATAN(1.D0)/8.D0) + 5
IC = 2*IDINT(HALFM*(0.5D0-DSQRT(3.D0)/6.D0)) + 1
MIC = (M2 - IC) + M2
C
C S IS THE SCALE FACTOR FOR CONVERTING TO FLOATING POINT
C
S = 0.5D0/HALFM
C
C COMPUTE NEXT RANDOM NUMBER
C
20 IY = IY*IA
C
C THE FOLLOWING STATEMENT IS FOR COMPUTERS WHICH DO NOT ALLOW
C INTEGER OVERFLOW ON ADDITION
C
IF (IY .GT. MIC) IY = (IY - M2) - M2
C
IY = IY + IC
C
C THE FOLLOWING STATEMENT IS FOR COMPUTERS WHERE THE
C WORD LENGTH FOR ADDITION IS GREATER THAN FOR MULTIPLICATION
C
IF (IY/2 .GT. M2) IY = (IY - M2) - M2
C
C THE FOLLOWING STATEMENT IS FOR COMPUTERS WHERE INTEGER
C OVERFLOW AFFECTS THE SIGN BIT
C
IF (IY .LT. 0) IY = (IY + M2) + M2
URAND = DFLOAT(IY)*S
RETURN
END
SUBROUTINE WMUL(AR,AI,BR,BI,CR,CI)
DOUBLE PRECISION AR,AI,BR,BI,CR,CI,T,FLOP
C C = A*B
T = AR*BI + AI*BR
IF (T .NE. 0.0D0) T = FLOP(T)
CR = FLOP(AR*BR - AI*BI)
CI = T
RETURN
END
SUBROUTINE WDIV(AR,AI,BR,BI,CR,CI)
DOUBLE PRECISION AR,AI,BR,BI,CR,CI
C C = A/B
DOUBLE PRECISION S,D,ARS,AIS,BRS,BIS,FLOP
S = DABS(BR) + DABS(BI)
IF (S .EQ. 0.0D0) CALL ERROR(27)
IF (S .EQ. 0.0D0) RETURN
ARS = AR/S
AIS = AI/S
BRS = BR/S
BIS = BI/S
D = BRS**2 + BIS**2
SHAR_EOF
# End of shell archive
exit 0
--
Bob Page, U of Lowell CS Dept. page@swan.ulowell.edu ulowell!page
Have five nice days.