*DASUM DOUBLE PRECISION FUNCTION DASUM(N,DX,INCX) C***BEGIN PROLOGUE DASUM C***DATE WRITTEN 791001 (YYMMDD) C***REVISION DATE 820801 (YYMMDD) C***CATEGORY NO. D1A3A C***KEYWORDS ADD,BLAS,DOUBLE PRECISION,LINEAR ALGEBRA,MAGNITUDE,SUM, C VECTOR C***AUTHOR LAWSON, C. L., (JPL) C HANSON, R. J., (SNLA) C KINCAID, D. R., (U. OF TEXAS) C KROGH, F. T., (JPL) C***PURPOSE SUM OF MAGNITUDES OF D.P. VECTOR COMPONENTS C***DESCRIPTION C B L A S SUBPROGRAM C DESCRIPTION OF PARAMETERS C --INPUT-- C N NUMBER OF ELEMENTS IN INPUT VECTOR(S) C DX DOUBLE PRECISION VECTOR WITH N ELEMENTS C INCX STORAGE SPACING BETWEEN ELEMENTS OF DX C --OUTPUT-- C DASUM DOUBLE PRECISION RESULT (ZERO IF N .LE. 0) C RETURNS SUM OF MAGNITUDES OF DOUBLE PRECISION DX. C DASUM = SUM FROM 0 TO N-1 OF DABS(DX(1+I*INCX)) C***REFERENCES LAWSON C.L., HANSON R.J., KINCAID D.R., KROGH F.T., C *BASIC LINEAR ALGEBRA SUBPROGRAMS FOR FORTRAN USAGE*, C ALGORITHM NO. 539, TRANSACTIONS ON MATHEMATICAL C SOFTWARE, VOLUME 5, NUMBER 3, SEPTEMBER 1979, 308-323 C***ROUTINES CALLED (NONE) C***END PROLOGUE DASUM C...SCALAR ARGUMENTS INTEGER + INCX,N C...ARRAY ARGUMENTS DOUBLE PRECISION + DX(*) C...LOCAL SCALARS INTEGER + I,M,MP1,NS C...INTRINSIC FUNCTIONS INTRINSIC + DABS,MOD C***FIRST EXECUTABLE STATEMENT DASUM DASUM = 0.D0 IF(N.LE.0)RETURN IF(INCX.EQ.1)GOTO 20 C CODE FOR INCREMENTS NOT EQUAL TO 1. NS = N*INCX DO 10 I=1,NS,INCX DASUM = DASUM + DABS(DX(I)) 10 CONTINUE RETURN C CODE FOR INCREMENTS EQUAL TO 1. C CLEAN-UP LOOP SO REMAINING VECTOR LENGTH IS A MULTIPLE OF 6. 20 M = MOD(N,6) IF( M .EQ. 0 ) GO TO 40 DO 30 I = 1,M DASUM = DASUM + DABS(DX(I)) 30 CONTINUE IF( N .LT. 6 ) RETURN 40 MP1 = M + 1 DO 50 I = MP1,N,6 DASUM = DASUM + DABS(DX(I)) + DABS(DX(I+1)) + DABS(DX(I+2)) 1 + DABS(DX(I+3)) + DABS(DX(I+4)) + DABS(DX(I+5)) 50 CONTINUE RETURN END *DAXPY SUBROUTINE DAXPY(N,DA,DX,INCX,DY,INCY) C***BEGIN PROLOGUE DAXPY C***DATE WRITTEN 791001 (YYMMDD) C***REVISION DATE 820801 (YYMMDD) C***CATEGORY NO. D1A7 C***KEYWORDS BLAS,DOUBLE PRECISION,LINEAR ALGEBRA,TRIAD,VECTOR C***AUTHOR LAWSON, C. L., (JPL) C HANSON, R. J., (SNLA) C KINCAID, D. R., (U. OF TEXAS) C KROGH, F. T., (JPL) C***PURPOSE D.P COMPUTATION Y = A*X + Y C***DESCRIPTION C B L A S SUBPROGRAM C DESCRIPTION OF PARAMETERS C --INPUT-- C N NUMBER OF ELEMENTS IN INPUT VECTOR(S) C DA DOUBLE PRECISION SCALAR MULTIPLIER C DX DOUBLE PRECISION VECTOR WITH N ELEMENTS C INCX STORAGE SPACING BETWEEN ELEMENTS OF DX C DY DOUBLE PRECISION VECTOR WITH N ELEMENTS C INCY STORAGE SPACING BETWEEN ELEMENTS OF DY C --OUTPUT-- C DY DOUBLE PRECISION RESULT (UNCHANGED IF N .LE. 0) C OVERWRITE DOUBLE PRECISION DY WITH DOUBLE PRECISION DA*DX + DY. C FOR I = 0 TO N-1, REPLACE DY(LY+I*INCY) WITH DA*DX(LX+I*INCX) + C DY(LY+I*INCY), WHERE LX = 1 IF INCX .GE. 0, ELSE LX = (-INCX)*N C AND LY IS DEFINED IN A SIMILAR WAY USING INCY. C***REFERENCES LAWSON C.L., HANSON R.J., KINCAID D.R., KROGH F.T., C *BASIC LINEAR ALGEBRA SUBPROGRAMS FOR FORTRAN USAGE*, C ALGORITHM NO. 539, TRANSACTIONS ON MATHEMATICAL C SOFTWARE, VOLUME 5, NUMBER 3, SEPTEMBER 1979, 308-323 C***ROUTINES CALLED (NONE) C***END PROLOGUE DAXPY C...SCALAR ARGUMENTS DOUBLE PRECISION + DA INTEGER + INCX,INCY,N C...ARRAY ARGUMENTS DOUBLE PRECISION + DX(*),DY(*) C...LOCAL SCALARS INTEGER + I,IX,IY,M,MP1,NS C...INTRINSIC FUNCTIONS INTRINSIC + MOD C***FIRST EXECUTABLE STATEMENT DAXPY IF(N.LE.0.OR.DA.EQ.0.D0) RETURN IF(INCX.EQ.INCY) IF(INCX-1) 5,20,60 5 CONTINUE C CODE FOR NONEQUAL OR NONPOSITIVE INCREMENTS. 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 CODE FOR BOTH INCREMENTS EQUAL TO 1 C CLEAN-UP LOOP SO REMAINING VECTOR LENGTH IS A MULTIPLE OF 4. 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 C CODE FOR EQUAL, POSITIVE, NONUNIT INCREMENTS. 60 CONTINUE NS = N*INCX DO 70 I=1,NS,INCX DY(I) = DA*DX(I) + DY(I) 70 CONTINUE RETURN END *DCHEX SUBROUTINE DCHEX(R,LDR,P,K,L,Z,LDZ,NZ,C,S,JOB) C***BEGIN PROLOGUE DCHEX C***DATE WRITTEN 780814 (YYMMDD) C***REVISION DATE 820801 (YYMMDD) C***CATEGORY NO. D7B C***KEYWORDS CHOLESKY DECOMPOSITION,DOUBLE PRECISION,EXCHANGE, C LINEAR ALGEBRA,LINPACK,MATRIX,POSITIVE DEFINITE C***AUTHOR STEWART, G. W., (U. OF MARYLAND) C***PURPOSE UPDATES THE CHOLESKY FACTORIZATION A=TRANS(R)*R OF A C POSITIVE DEFINITE MATRIX A OF ORDER P UNDER DIAGONAL C PERMUTATIONS OF THE FORM TRANS(E)*A*E WHERE E IS A C PERMUTATION MATRIX. C***DESCRIPTION C DCHEX UPDATES THE CHOLESKY FACTORIZATION C A = TRANS(R)*R C OF A POSITIVE DEFINITE MATRIX A OF ORDER P UNDER DIAGONAL C PERMUTATIONS OF THE FORM C TRANS(E)*A*E C WHERE E IS A PERMUTATION MATRIX. SPECIFICALLY, GIVEN C AN UPPER TRIANGULAR MATRIX R AND A PERMUTATION MATRIX C E (WHICH IS SPECIFIED BY K, L, AND JOB), DCHEX DETERMINES C AN ORTHOGONAL MATRIX U SUCH THAT C U*R*E = RR, C WHERE RR IS UPPER TRIANGULAR. AT THE USERS OPTION, THE C TRANSFORMATION U WILL BE MULTIPLIED INTO THE ARRAY Z. C IF A = TRANS(X)*X, SO THAT R IS THE TRIANGULAR PART OF THE C QR FACTORIZATION OF X, THEN RR IS THE TRIANGULAR PART OF THE C QR FACTORIZATION OF X*E, I.E. X WITH ITS COLUMNS PERMUTED. C FOR A LESS TERSE DESCRIPTION OF WHAT DCHEX DOES AND HOW C IT MAY BE APPLIED, SEE THE LINPACK GUIDE. C THE MATRIX Q IS DETERMINED AS THE PRODUCT U(L-K)*...*U(1) C OF PLANE ROTATIONS OF THE FORM C ( C(I) S(I) ) C ( ) , C ( -S(I) C(I) ) C WHERE C(I) IS DOUBLE PRECISION. THE ROWS THESE ROTATIONS OPERATE C ON ARE DESCRIBED BELOW. C THERE ARE TWO TYPES OF PERMUTATIONS, WHICH ARE DETERMINED C BY THE VALUE OF JOB. C 1. RIGHT CIRCULAR SHIFT (JOB = 1). C THE COLUMNS ARE REARRANGED IN THE FOLLOWING ORDER. C 1,...,K-1,L,K,K+1,...,L-1,L+1,...,P. C U IS THE PRODUCT OF L-K ROTATIONS U(I), WHERE U(I) C ACTS IN THE (L-I,L-I+1)-PLANE. C 2. LEFT CIRCULAR SHIFT (JOB = 2). C THE COLUMNS ARE REARRANGED IN THE FOLLOWING ORDER C 1,...,K-1,K+1,K+2,...,L,K,L+1,...,P. C U IS THE PRODUCT OF L-K ROTATIONS U(I), WHERE U(I) C ACTS IN THE (K+I-1,K+I)-PLANE. C ON ENTRY C R DOUBLE PRECISION(LDR,P), WHERE LDR .GE. P. C R CONTAINS THE UPPER TRIANGULAR FACTOR C THAT IS TO BE UPDATED. ELEMENTS OF R C BELOW THE DIAGONAL ARE NOT REFERENCED. C LDR INTEGER. C LDR IS THE LEADING DIMENSION OF THE ARRAY R. C P INTEGER. C P IS THE ORDER OF THE MATRIX R. C K INTEGER. C K IS THE FIRST COLUMN TO BE PERMUTED. C L INTEGER. C L IS THE LAST COLUMN TO BE PERMUTED. C L MUST BE STRICTLY GREATER THAN K. C Z DOUBLE PRECISION(LDZ,N)Z), WHERE LDZ .GE. P. C Z IS AN ARRAY OF NZ P-VECTORS INTO WHICH THE C TRANSFORMATION U IS MULTIPLIED. Z IS C NOT REFERENCED IF NZ = 0. C LDZ INTEGER. C LDZ IS THE LEADING DIMENSION OF THE ARRAY Z. C NZ INTEGER. C NZ IS THE NUMBER OF COLUMNS OF THE MATRIX Z. C JOB INTEGER. C JOB DETERMINES THE TYPE OF PERMUTATION. C JOB = 1 RIGHT CIRCULAR SHIFT. C JOB = 2 LEFT CIRCULAR SHIFT. C ON RETURN C R CONTAINS THE UPDATED FACTOR. C Z CONTAINS THE UPDATED MATRIX Z. C C DOUBLE PRECISION(P). C C CONTAINS THE COSINES OF THE TRANSFORMING ROTATIONS. C S DOUBLE PRECISION(P). C S CONTAINS THE SINES OF THE TRANSFORMING ROTATIONS. C LINPACK. THIS VERSION DATED 08/14/78 . C G. W. STEWART, UNIVERSITY OF MARYLAND, ARGONNE NATIONAL LAB. C***REFERENCES DONGARRA J.J., BUNCH J.R., MOLER C.B., STEWART G.W., C *LINPACK USERS GUIDE*, SIAM, 1979. C***ROUTINES CALLED DROTG C***END PROLOGUE DCHEX C...SCALAR ARGUMENTS INTEGER + JOB,K,L,LDR,LDZ,NZ,P C...ARRAY ARGUMENTS DOUBLE PRECISION + C(*),R(LDR,*),S(*),Z(LDZ,*) C...LOCAL SCALARS DOUBLE PRECISION + T,T1 INTEGER + I,II,IL,IU,J,JJ,KM1,KP1,LM1,LMK C...EXTERNAL SUBROUTINES EXTERNAL + DROTG C...INTRINSIC FUNCTIONS INTRINSIC + MAX0,MIN0 C***FIRST EXECUTABLE STATEMENT DCHEX KM1 = K - 1 KP1 = K + 1 LMK = L - K LM1 = L - 1 C PERFORM THE APPROPRIATE TASK. GO TO (10,130), JOB C RIGHT CIRCULAR SHIFT. 10 CONTINUE C REORDER THE COLUMNS. DO 20 I = 1, L II = L - I + 1 S(I) = R(II,L) 20 CONTINUE DO 40 JJ = K, LM1 J = LM1 - JJ + K DO 30 I = 1, J R(I,J+1) = R(I,J) 30 CONTINUE R(J+1,J+1) = 0.0D0 40 CONTINUE IF (K .EQ. 1) GO TO 60 DO 50 I = 1, KM1 II = L - I + 1 R(I,K) = S(II) 50 CONTINUE 60 CONTINUE C CALCULATE THE ROTATIONS. T = S(1) DO 70 I = 1, LMK T1 = S(I) CALL DROTG(S(I+1),T,C(I),T1) S(I) = T1 T = S(I+1) 70 CONTINUE R(K,K) = T DO 90 J = KP1, P IL = MAX0(1,L-J+1) DO 80 II = IL, LMK I = L - II T = C(II)*R(I,J) + S(II)*R(I+1,J) R(I+1,J) = C(II)*R(I+1,J) - S(II)*R(I,J) R(I,J) = T 80 CONTINUE 90 CONTINUE C IF REQUIRED, APPLY THE TRANSFORMATIONS TO Z. IF (NZ .LT. 1) GO TO 120 DO 110 J = 1, NZ DO 100 II = 1, LMK I = L - II T = C(II)*Z(I,J) + S(II)*Z(I+1,J) Z(I+1,J) = C(II)*Z(I+1,J) - S(II)*Z(I,J) Z(I,J) = T 100 CONTINUE 110 CONTINUE 120 CONTINUE GO TO 260 C LEFT CIRCULAR SHIFT 130 CONTINUE C REORDER THE COLUMNS DO 140 I = 1, K II = LMK + I S(II) = R(I,K) 140 CONTINUE DO 160 J = K, LM1 DO 150 I = 1, J R(I,J) = R(I,J+1) 150 CONTINUE JJ = J - KM1 S(JJ) = R(J+1,J+1) 160 CONTINUE DO 170 I = 1, K II = LMK + I R(I,L) = S(II) 170 CONTINUE DO 180 I = KP1, L R(I,L) = 0.0D0 180 CONTINUE C REDUCTION LOOP. DO 220 J = K, P IF (J .EQ. K) GO TO 200 C APPLY THE ROTATIONS. IU = MIN0(J-1,L-1) DO 190 I = K, IU II = I - K + 1 T = C(II)*R(I,J) + S(II)*R(I+1,J) R(I+1,J) = C(II)*R(I+1,J) - S(II)*R(I,J) R(I,J) = T 190 CONTINUE 200 CONTINUE IF (J .GE. L) GO TO 210 JJ = J - K + 1 T = S(JJ) CALL DROTG(R(J,J),T,C(JJ),S(JJ)) 210 CONTINUE 220 CONTINUE C APPLY THE ROTATIONS TO Z. IF (NZ .LT. 1) GO TO 250 DO 240 J = 1, NZ DO 230 I = K, LM1 II = I - KM1 T = C(II)*Z(I,J) + S(II)*Z(I+1,J) Z(I+1,J) = C(II)*Z(I+1,J) - S(II)*Z(I,J) Z(I,J) = T 230 CONTINUE 240 CONTINUE 250 CONTINUE 260 CONTINUE RETURN END *DCOPY SUBROUTINE DCOPY(N,DX,INCX,DY,INCY) C***BEGIN PROLOGUE DCOPY C***DATE WRITTEN 791001 (YYMMDD) C***REVISION DATE 820801 (YYMMDD) C***CATEGORY NO. D1A5 C***KEYWORDS BLAS,COPY,DOUBLE PRECISION,LINEAR ALGEBRA,VECTOR C***AUTHOR LAWSON, C. L., (JPL) C HANSON, R. J., (SNLA) C KINCAID, D. R., (U. OF TEXAS) C KROGH, F. T., (JPL) C***PURPOSE D.P. VECTOR COPY Y = X C***DESCRIPTION C B L A S SUBPROGRAM C DESCRIPTION OF PARAMETERS C --INPUT-- C N NUMBER OF ELEMENTS IN INPUT VECTOR(S) C DX DOUBLE PRECISION VECTOR WITH N ELEMENTS C INCX STORAGE SPACING BETWEEN ELEMENTS OF DX C DY DOUBLE PRECISION VECTOR WITH N ELEMENTS C INCY STORAGE SPACING BETWEEN ELEMENTS OF DY C --OUTPUT-- C DY COPY OF VECTOR DX (UNCHANGED IF N .LE. 0) C COPY DOUBLE PRECISION DX TO DOUBLE PRECISION DY. C FOR I = 0 TO N-1, COPY DX(LX+I*INCX) TO DY(LY+I*INCY), C WHERE LX = 1 IF INCX .GE. 0, ELSE LX = (-INCX)*N, AND LY IS C DEFINED IN A SIMILAR WAY USING INCY. C***REFERENCES LAWSON C.L., HANSON R.J., KINCAID D.R., KROGH F.T., C *BASIC LINEAR ALGEBRA SUBPROGRAMS FOR FORTRAN USAGE*, C ALGORITHM NO. 539, TRANSACTIONS ON MATHEMATICAL C SOFTWARE, VOLUME 5, NUMBER 3, SEPTEMBER 1979, 308-323 C***ROUTINES CALLED (NONE) C***END PROLOGUE DCOPY C...SCALAR ARGUMENTS INTEGER + INCX,INCY,N C...ARRAY ARGUMENTS DOUBLE PRECISION + DX(*),DY(*) C...LOCAL SCALARS INTEGER + I,IX,IY,M,MP1,NS C...INTRINSIC FUNCTIONS INTRINSIC + MOD C***FIRST EXECUTABLE STATEMENT DCOPY IF(N.LE.0)RETURN IF(INCX.EQ.INCY) IF(INCX-1) 5,20,60 5 CONTINUE C CODE FOR UNEQUAL OR NONPOSITIVE INCREMENTS. 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) = DX(IX) IX = IX + INCX IY = IY + INCY 10 CONTINUE RETURN C CODE FOR BOTH INCREMENTS EQUAL TO 1 C CLEAN-UP LOOP SO REMAINING VECTOR LENGTH IS A MULTIPLE OF 7. 20 M = MOD(N,7) IF( M .EQ. 0 ) GO TO 40 DO 30 I = 1,M DY(I) = DX(I) 30 CONTINUE IF( N .LT. 7 ) RETURN 40 MP1 = M + 1 DO 50 I = MP1,N,7 DY(I) = DX(I) DY(I + 1) = DX(I + 1) DY(I + 2) = DX(I + 2) DY(I + 3) = DX(I + 3) DY(I + 4) = DX(I + 4) DY(I + 5) = DX(I + 5) DY(I + 6) = DX(I + 6) 50 CONTINUE RETURN C CODE FOR EQUAL, POSITIVE, NONUNIT INCREMENTS. 60 CONTINUE NS=N*INCX DO 70 I=1,NS,INCX DY(I) = DX(I) 70 CONTINUE RETURN END *DDOT DOUBLE PRECISION FUNCTION DDOT(N,DX,INCX,DY,INCY) C***BEGIN PROLOGUE DDOT C***DATE WRITTEN 791001 (YYMMDD) C***REVISION DATE 820801 (YYMMDD) C***CATEGORY NO. D1A4 C***KEYWORDS BLAS,DOUBLE PRECISION,INNER PRODUCT,LINEAR ALGEBRA,VECTOR C***AUTHOR LAWSON, C. L., (JPL) C HANSON, R. J., (SNLA) C KINCAID, D. R., (U. OF TEXAS) C KROGH, F. T., (JPL) C***PURPOSE D.P. INNER PRODUCT OF D.P. VECTORS C***DESCRIPTION C B L A S SUBPROGRAM C DESCRIPTION OF PARAMETERS C --INPUT-- C N NUMBER OF ELEMENTS IN INPUT VECTOR(S) C DX DOUBLE PRECISION VECTOR WITH N ELEMENTS C INCX STORAGE SPACING BETWEEN ELEMENTS OF DX C DY DOUBLE PRECISION VECTOR WITH N ELEMENTS C INCY STORAGE SPACING BETWEEN ELEMENTS OF DY C --OUTPUT-- C DDOT DOUBLE PRECISION DOT PRODUCT (ZERO IF N .LE. 0) C RETURNS THE DOT PRODUCT OF DOUBLE PRECISION DX AND DY. C DDOT = SUM FOR I = 0 TO N-1 OF DX(LX+I*INCX) * DY(LY+I*INCY) C WHERE LX = 1 IF INCX .GE. 0, ELSE LX = (-INCX)*N, AND LY IS C DEFINED IN A SIMILAR WAY USING INCY. C***REFERENCES LAWSON C.L., HANSON R.J., KINCAID D.R., KROGH F.T., C *BASIC LINEAR ALGEBRA SUBPROGRAMS FOR FORTRAN USAGE*, C ALGORITHM NO. 539, TRANSACTIONS ON MATHEMATICAL C SOFTWARE, VOLUME 5, NUMBER 3, SEPTEMBER 1979, 308-323 C***ROUTINES CALLED (NONE) C***END PROLOGUE DDOT C...SCALAR ARGUMENTS INTEGER + INCX,INCY,N C...ARRAY ARGUMENTS DOUBLE PRECISION + DX(*),DY(*) C...LOCAL SCALARS INTEGER + I,IX,IY,M,MP1,NS C...INTRINSIC FUNCTIONS INTRINSIC + MOD C***FIRST EXECUTABLE STATEMENT DDOT DDOT = 0.D0 IF(N.LE.0)RETURN IF(INCX.EQ.INCY) IF(INCX-1) 5,20,60 5 CONTINUE C CODE FOR UNEQUAL OR NONPOSITIVE INCREMENTS. 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 DDOT = DDOT + DX(IX)*DY(IY) IX = IX + INCX IY = IY + INCY 10 CONTINUE RETURN C CODE FOR BOTH INCREMENTS EQUAL TO 1. C CLEAN-UP LOOP SO REMAINING VECTOR LENGTH IS A MULTIPLE OF 5. 20 M = MOD(N,5) IF( M .EQ. 0 ) GO TO 40 DO 30 I = 1,M DDOT = DDOT + DX(I)*DY(I) 30 CONTINUE IF( N .LT. 5 ) RETURN 40 MP1 = M + 1 DO 50 I = MP1,N,5 DDOT = DDOT + DX(I)*DY(I) + DX(I+1)*DY(I+1) + 1 DX(I + 2)*DY(I + 2) + DX(I + 3)*DY(I + 3) + DX(I + 4)*DY(I + 4) 50 CONTINUE RETURN C CODE FOR POSITIVE EQUAL INCREMENTS .NE.1. 60 CONTINUE NS = N*INCX DO 70 I=1,NS,INCX DDOT = DDOT + DX(I)*DY(I) 70 CONTINUE RETURN END *DNRM2 DOUBLE PRECISION FUNCTION DNRM2(N,DX,INCX) C***BEGIN PROLOGUE DNRM2 C***DATE WRITTEN 791001 (YYMMDD) C***REVISION DATE 820801 (YYMMDD) C***CATEGORY NO. D1A3B C***KEYWORDS BLAS,DOUBLE PRECISION,EUCLIDEAN,L2,LENGTH,LINEAR ALGEBRA, C NORM,VECTOR C***AUTHOR LAWSON, C. L., (JPL) C HANSON, R. J., (SNLA) C KINCAID, D. R., (U. OF TEXAS) C KROGH, F. T., (JPL) C***PURPOSE EUCLIDEAN LENGTH (L2 NORM) OF D.P. VECTOR C***DESCRIPTION C B L A S SUBPROGRAM C DESCRIPTION OF PARAMETERS C --INPUT-- C N NUMBER OF ELEMENTS IN INPUT VECTOR(S) C DX DOUBLE PRECISION VECTOR WITH N ELEMENTS C INCX STORAGE SPACING BETWEEN ELEMENTS OF DX C --OUTPUT-- C DNRM2 DOUBLE PRECISION RESULT (ZERO IF N .LE. 0) C EUCLIDEAN NORM OF THE N-VECTOR STORED IN DX() WITH STORAGE C INCREMENT INCX . C IF N .LE. 0 RETURN WITH RESULT = 0. C IF N .GE. 1 THEN INCX MUST BE .GE. 1 C C.L. LAWSON, 1978 JAN 08 C FOUR PHASE METHOD USING TWO BUILT-IN CONSTANTS THAT ARE C HOPEFULLY APPLICABLE TO ALL MACHINES. C CUTLO = MAXIMUM OF DSQRT(U/EPS) OVER ALL KNOWN MACHINES. C CUTHI = MINIMUM OF DSQRT(V) OVER ALL KNOWN MACHINES. C WHERE C EPS = SMALLEST NO. SUCH THAT EPS + 1. .GT. 1. C U = SMALLEST POSITIVE NO. (UNDERFLOW LIMIT) C V = LARGEST NO. (OVERFLOW LIMIT) C BRIEF OUTLINE OF ALGORITHM.. C PHASE 1 SCANS ZERO COMPONENTS. C MOVE TO PHASE 2 WHEN A COMPONENT IS NONZERO AND .LE. CUTLO C MOVE TO PHASE 3 WHEN A COMPONENT IS .GT. CUTLO C MOVE TO PHASE 4 WHEN A COMPONENT IS .GE. CUTHI/M C WHERE M = N FOR X() REAL AND M = 2*N FOR COMPLEX. C VALUES FOR CUTLO AND CUTHI.. C FROM THE ENVIRONMENTAL PARAMETERS LISTED IN THE IMSL CONVERTER C DOCUMENT THE LIMITING VALUES ARE AS FOLLOWS.. C CUTLO, S.P. U/EPS = 2**(-102) FOR HONEYWELL. CLOSE SECONDS ARE C UNIVAC AND DEC AT 2**(-103) C THUS CUTLO = 2**(-51) = 4.44089E-16 C CUTHI, S.P. V = 2**127 FOR UNIVAC, HONEYWELL, AND DEC. C THUS CUTHI = 2**(63.5) = 1.30438E19 C CUTLO, D.P. U/EPS = 2**(-67) FOR HONEYWELL AND DEC. C THUS CUTLO = 2**(-33.5) = 8.23181D-11 C CUTHI, D.P. SAME AS S.P. CUTHI = 1.30438D19 C DATA CUTLO, CUTHI / 8.232D-11, 1.304D19 / C DATA CUTLO, CUTHI / 4.441E-16, 1.304E19 / C***REFERENCES LAWSON C.L., HANSON R.J., KINCAID D.R., KROGH F.T., C *BASIC LINEAR ALGEBRA SUBPROGRAMS FOR FORTRAN USAGE*, C ALGORITHM NO. 539, TRANSACTIONS ON MATHEMATICAL C SOFTWARE, VOLUME 5, NUMBER 3, SEPTEMBER 1979, 308-323 C***ROUTINES CALLED (NONE) C***END PROLOGUE DNRM2 C...SCALAR ARGUMENTS INTEGER + INCX,N C...ARRAY ARGUMENTS DOUBLE PRECISION + DX(*) C...LOCAL SCALARS DOUBLE PRECISION + CUTHI,CUTLO,HITEST,ONE,SUM,XMAX,ZERO INTEGER + I,J,NEXT,NN C...INTRINSIC FUNCTIONS INTRINSIC + DABS,DSQRT,FLOAT C...DATA STATEMENTS DATA + ZERO,ONE/0.0D0,1.0D0/ DATA + CUTLO,CUTHI/8.232D-11,1.304D19/ C***FIRST EXECUTABLE STATEMENT DNRM2 XMAX = ZERO IF(N .GT. 0) GO TO 10 DNRM2 = ZERO GO TO 300 10 ASSIGN 30 TO NEXT SUM = ZERO NN = N * INCX C BEGIN MAIN LOOP I = 1 C 20 GO TO NEXT,(30, 50, 70, 110) 20 GO TO NEXT 30 IF( DABS(DX(I)) .GT. CUTLO) GO TO 85 ASSIGN 50 TO NEXT XMAX = ZERO C PHASE 1. SUM IS ZERO 50 IF( DX(I) .EQ. ZERO) GO TO 200 IF( DABS(DX(I)) .GT. CUTLO) GO TO 85 C PREPARE FOR PHASE 2. ASSIGN 70 TO NEXT GO TO 105 C PREPARE FOR PHASE 4. 100 I = J ASSIGN 110 TO NEXT SUM = (SUM / DX(I)) / DX(I) 105 XMAX = DABS(DX(I)) GO TO 115 C PHASE 2. SUM IS SMALL. C SCALE TO AVOID DESTRUCTIVE UNDERFLOW. 70 IF( DABS(DX(I)) .GT. CUTLO ) GO TO 75 C COMMON CODE FOR PHASES 2 AND 4. C IN PHASE 4 SUM IS LARGE. SCALE TO AVOID OVERFLOW. 110 IF( DABS(DX(I)) .LE. XMAX ) GO TO 115 SUM = ONE + SUM * (XMAX / DX(I))**2 XMAX = DABS(DX(I)) GO TO 200 115 SUM = SUM + (DX(I)/XMAX)**2 GO TO 200 C PREPARE FOR PHASE 3. 75 SUM = (SUM * XMAX) * XMAX C FOR REAL OR D.P. SET HITEST = CUTHI/N C FOR COMPLEX SET HITEST = CUTHI/(2*N) 85 HITEST = CUTHI/FLOAT( N ) C PHASE 3. SUM IS MID-RANGE. NO SCALING. DO 95 J =I,NN,INCX IF(DABS(DX(J)) .GE. HITEST) GO TO 100 95 SUM = SUM + DX(J)**2 DNRM2 = DSQRT( SUM ) GO TO 300 200 CONTINUE I = I + INCX IF ( I .LE. NN ) GO TO 20 C END OF MAIN LOOP. C COMPUTE SQUARE ROOT AND ADJUST FOR SCALING. DNRM2 = XMAX * DSQRT(SUM) 300 CONTINUE RETURN END *DPODI SUBROUTINE DPODI(A,LDA,N,DET,JOB) C***BEGIN PROLOGUE DPODI C***DATE WRITTEN 780814 (YYMMDD) C***REVISION DATE 820801 (YYMMDD) C***CATEGORY NO. D2B1B,D3B1B C***KEYWORDS DETERMINANT,DOUBLE PRECISION,FACTOR,INVERSE, C LINEAR ALGEBRA,LINPACK,MATRIX,POSITIVE DEFINITE C***AUTHOR MOLER, C. B., (U. OF NEW MEXICO) C***PURPOSE COMPUTES THE DETERMINANT AND INVERSE OF A CERTAIN DOUBLE C PRECISION SYMMETRIC POSITIVE DEFINITE MATRIX (SEE ABSTRACT) C USING THE FACTORS COMPUTED BY DPOCO, DPOFA OR DQRDC. C***DESCRIPTION C DPODI COMPUTES THE DETERMINANT AND INVERSE OF A CERTAIN C DOUBLE PRECISION SYMMETRIC POSITIVE DEFINITE MATRIX (SEE BELOW) C USING THE FACTORS COMPUTED BY DPOCO, DPOFA OR DQRDC. C ON ENTRY C A DOUBLE PRECISION(LDA, N) C THE OUTPUT A FROM DPOCO OR DPOFA C OR THE OUTPUT X FROM DQRDC. C LDA INTEGER C THE LEADING DIMENSION OF THE ARRAY A . C N INTEGER C THE ORDER OF THE MATRIX A . C JOB INTEGER C = 11 BOTH DETERMINANT AND INVERSE. C = 01 INVERSE ONLY. C = 10 DETERMINANT ONLY. C ON RETURN C A IF DPOCO OR DPOFA WAS USED TO FACTOR A , THEN C DPODI PRODUCES THE UPPER HALF OF INVERSE(A) . C IF DQRDC WAS USED TO DECOMPOSE X , THEN C DPODI PRODUCES THE UPPER HALF OF INVERSE(TRANS(X)*X) C WHERE TRANS(X) IS THE TRANSPOSE. C ELEMENTS OF A BELOW THE DIAGONAL ARE UNCHANGED. C IF THE UNITS DIGIT OF JOB IS ZERO, A IS UNCHANGED. C DET DOUBLE PRECISION(2) C DETERMINANT OF A OR OF TRANS(X)*X IF REQUESTED. C OTHERWISE NOT REFERENCED. C DETERMINANT = DET(1) * 10.0**DET(2) C WITH 1.0 .LE. DET(1) .LT. 10.0 C OR DET(1) .EQ. 0.0 . C ERROR CONDITION 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 DPOCO OR DPOFA HAS SET INFO .EQ. 0 . C LINPACK. THIS VERSION DATED 08/14/78 . C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB. C***REFERENCES DONGARRA J.J., BUNCH J.R., MOLER C.B., STEWART G.W., C *LINPACK USERS GUIDE*, SIAM, 1979. C***ROUTINES CALLED DAXPY,DSCAL C***END PROLOGUE DPODI C...SCALAR ARGUMENTS INTEGER JOB,LDA,N C...ARRAY ARGUMENTS DOUBLE PRECISION A(LDA,*),DET(*) C...LOCAL SCALARS DOUBLE PRECISION S,T INTEGER I,J,JM1,K,KP1 C...EXTERNAL SUBROUTINES EXTERNAL DAXPY,DSCAL C...INTRINSIC FUNCTIONS INTRINSIC MOD C***FIRST EXECUTABLE STATEMENT DPODI IF (JOB/10 .EQ. 0) GO TO 70 DET(1) = 1.0D0 DET(2) = 0.0D0 S = 10.0D0 DO 50 I = 1, N DET(1) = A(I,I)**2*DET(1) C ...EXIT IF (DET(1) .EQ. 0.0D0) GO TO 60 10 IF (DET(1) .GE. 1.0D0) GO TO 20 DET(1) = S*DET(1) DET(2) = DET(2) - 1.0D0 GO TO 10 20 CONTINUE 30 IF (DET(1) .LT. S) GO TO 40 DET(1) = DET(1)/S DET(2) = DET(2) + 1.0D0 GO TO 30 40 CONTINUE 50 CONTINUE 60 CONTINUE 70 CONTINUE C COMPUTE INVERSE(R) IF (MOD(JOB,10) .EQ. 0) GO TO 140 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 FORM INVERSE(R) * TRANS(INVERSE(R)) DO 130 J = 1, N JM1 = J - 1 IF (JM1 .LT. 1) GO TO 120 DO 110 K = 1, JM1 T = A(K,J) CALL DAXPY(K,T,A(1,J),1,A(1,K),1) 110 CONTINUE 120 CONTINUE T = A(J,J) CALL DSCAL(J,T,A(1,J),1) 130 CONTINUE 140 CONTINUE RETURN END *DQRDC SUBROUTINE DQRDC(X,LDX,N,P,QRAUX,JPVT,WORK,JOB) C***BEGIN PROLOGUE DQRDC C***DATE WRITTEN 780814 (YYMMDD) C***REVISION DATE 820801 (YYMMDD) C***CATEGORY NO. D5 C***KEYWORDS DECOMPOSITION,DOUBLE PRECISION,LINEAR ALGEBRA,LINPACK, C MATRIX,ORTHOGONAL TRIANGULAR C***AUTHOR STEWART, G. W., (U. OF MARYLAND) C***PURPOSE USES HOUSEHOLDER TRANSFORMATIONS TO COMPUTE THE QR FACTORI- C ZATION OF N BY P MATRIX X. COLUMN PIVOTING IS OPTIONAL. C***DESCRIPTION C DQRDC 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 USER'S OPTION. C ON ENTRY C X DOUBLE PRECISION(LDX,P), WHERE LDX .GE. N. C X CONTAINS THE MATRIX WHOSE DECOMPOSITION IS TO BE C COMPUTED. C LDX INTEGER. C LDX IS THE LEADING DIMENSION OF THE ARRAY X. C N INTEGER. C N IS THE NUMBER OF ROWS OF THE MATRIX X. C P INTEGER. C P IS THE NUMBER OF COLUMNS OF THE MATRIX X. 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 IF JPVT(K) .GT. 0, THEN X(K) IS AN INITIAL C COLUMN. C IF JPVT(K) .EQ. 0, THEN X(K) IS A FREE COLUMN. C IF JPVT(K) .LT. 0, THEN X(K) IS A FINAL COLUMN. 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 WORK DOUBLE PRECISION(P). C WORK IS A WORK ARRAY. WORK IS NOT REFERENCED IF C JOB .EQ. 0. 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 ON RETURN 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 ORTHOGONAL 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 QRAUX DOUBLE PRECISION(P). C QRAUX CONTAINS FURTHER INFORMATION REQUIRED TO RECOVER C THE ORTHOGONAL PART OF THE DECOMPOSITION. 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 LINPACK. THIS VERSION DATED 08/14/78 . C G. W. STEWART, UNIVERSITY OF MARYLAND, ARGONNE NATIONAL LAB. C***REFERENCES DONGARRA J.J., BUNCH J.R., MOLER C.B., STEWART G.W., C *LINPACK USERS GUIDE*, SIAM, 1979. C***ROUTINES CALLED DAXPY,DDOT,DNRM2,DSCAL,DSWAP C***END PROLOGUE DQRDC C...SCALAR ARGUMENTS INTEGER + JOB,LDX,N,P C...ARRAY ARGUMENTS DOUBLE PRECISION + QRAUX(*),WORK(*),X(LDX,*) INTEGER + JPVT(*) C...LOCAL SCALARS DOUBLE PRECISION + MAXNRM,NRMXL,T,TT INTEGER + J,JJ,JP,L,LP1,LUP,MAXJ,PL,PU LOGICAL + NEGJ,SWAPJ C...EXTERNAL FUNCTIONS DOUBLE PRECISION + DDOT,DNRM2 EXTERNAL + DDOT,DNRM2 C...EXTERNAL SUBROUTINES EXTERNAL + DAXPY,DSCAL,DSWAP C...INTRINSIC FUNCTIONS INTRINSIC + DABS,DMAX1,DSIGN,DSQRT,MIN0 C***FIRST EXECUTABLE STATEMENT DQRDC PL = 1 PU = 0 IF (JOB .EQ. 0) GO TO 60 C PIVOTING HAS BEEN REQUESTED. REARRANGE THE COLUMNS C ACCORDING TO JPVT. 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 DSWAP(N,X(1,PL),1,X(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 DSWAP(N,X(1,PU),1,X(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 COMPUTE THE NORMS OF THE FREE COLUMNS. IF (PU .LT. PL) GO TO 80 DO 70 J = PL, PU QRAUX(J) = DNRM2(N,X(1,J),1) WORK(J) = QRAUX(J) 70 CONTINUE 80 CONTINUE C PERFORM THE HOUSEHOLDER REDUCTION OF X. LUP = MIN0(N,P) DO 200 L = 1, LUP IF (L .LT. PL .OR. L .GE. PU) GO TO 120 C LOCATE THE COLUMN OF LARGEST NORM AND BRING IT C INTO THE PIVOT POSITION. MAXNRM = 0.0D0 MAXJ = L DO 100 J = L, PU IF (QRAUX(J) .LE. MAXNRM) GO TO 90 MAXNRM = QRAUX(J) MAXJ = J 90 CONTINUE 100 CONTINUE IF (MAXJ .EQ. L) GO TO 110 CALL DSWAP(N,X(1,L),1,X(1,MAXJ),1) QRAUX(MAXJ) = QRAUX(L) WORK(MAXJ) = WORK(L) JP = JPVT(MAXJ) JPVT(MAXJ) = JPVT(L) JPVT(L) = JP 110 CONTINUE 120 CONTINUE QRAUX(L) = 0.0D0 IF (L .EQ. N) GO TO 190 C COMPUTE THE HOUSEHOLDER TRANSFORMATION FOR COLUMN L. NRMXL = DNRM2(N-L+1,X(L,L),1) IF (NRMXL .EQ. 0.0D0) GO TO 180 IF (X(L,L) .NE. 0.0D0) NRMXL = DSIGN(NRMXL,X(L,L)) CALL DSCAL(N-L+1,1.0D0/NRMXL,X(L,L),1) X(L,L) = 1.0D0 + X(L,L) C APPLY THE TRANSFORMATION TO THE REMAINING COLUMNS, C UPDATING THE NORMS. LP1 = L + 1 IF (P .LT. LP1) GO TO 170 DO 160 J = LP1, P T = -DDOT(N-L+1,X(L,L),1,X(L,J),1)/X(L,L) CALL DAXPY(N-L+1,T,X(L,L),1,X(L,J),1) IF (J .LT. PL .OR. J .GT. PU) GO TO 150 IF (QRAUX(J) .EQ. 0.0D0) GO TO 150 TT = 1.0D0 - (DABS(X(L,J))/QRAUX(J))**2 TT = DMAX1(TT,0.0D0) T = TT TT = 1.0D0 + 0.05D0*TT*(QRAUX(J)/WORK(J))**2 IF (TT .EQ. 1.0D0) GO TO 130 QRAUX(J) = QRAUX(J)*DSQRT(T) GO TO 140 130 CONTINUE QRAUX(J) = DNRM2(N-L,X(L+1,J),1) WORK(J) = QRAUX(J) 140 CONTINUE 150 CONTINUE 160 CONTINUE 170 CONTINUE C SAVE THE TRANSFORMATION. QRAUX(L) = X(L,L) X(L,L) = -NRMXL 180 CONTINUE 190 CONTINUE 200 CONTINUE RETURN END *DQRSL SUBROUTINE DQRSL(X,LDX,N,K,QRAUX,Y,QY,QTY,B,RSD,XB,JOB,INFO) C***BEGIN PROLOGUE DQRSL C***DATE WRITTEN 780814 (YYMMDD) C***REVISION DATE 820801 (YYMMDD) C***CATEGORY NO. D9,D2A1 C***KEYWORDS DOUBLE PRECISION,LINEAR ALGEBRA,LINPACK,MATRIX, C ORTHOGONAL TRIANGULAR,SOLVE C***AUTHOR STEWART, G. W., (U. OF MARYLAND) C***PURPOSE APPLIES THE OUTPUT OF DQRDC TO COMPUTE COORDINATE C TRANSFORMATIONS, PROJECTIONS, AND LEAST SQUARES SOLUTIONS. C***DESCRIPTION C DQRSL APPLIES THE OUTPUT OF DQRDC TO COMPUTE COORDINATE C TRANSFORMATIONS, PROJECTIONS, AND LEAST SQUARES SOLUTIONS. C FOR K .LE. MIN(N,P), LET XK BE THE MATRIX C XK = (X(JPVT(1)),X(JPVT(2)), ... ,X(JPVT(K))) C FORMED FROM COLUMNNS JPVT(1), ... ,JPVT(K) OF THE ORIGINAL C N X P MATRIX X THAT WAS INPUT TO DQRDC (IF NO PIVOTING WAS C DONE, XK CONSISTS OF THE FIRST K COLUMNS OF X IN THEIR C ORIGINAL ORDER). DQRDC PRODUCES A FACTORED ORTHOGONAL MATRIX Q C AND AN UPPER TRIANGULAR MATRIX R SUCH THAT C XK = Q * (R) C (0) C THIS INFORMATION IS CONTAINED IN CODED FORM IN THE ARRAYS C X AND QRAUX. C ON ENTRY C X DOUBLE PRECISION(LDX,P). C X CONTAINS THE OUTPUT OF DQRDC. C LDX INTEGER. C LDX IS THE LEADING DIMENSION OF THE ARRAY X. C N INTEGER. C N IS THE NUMBER OF ROWS OF THE MATRIX XK. IT MUST C HAVE THE SAME VALUE AS N IN DQRDC. C K INTEGER. C K IS THE NUMBER OF COLUMNS OF THE MATRIX XK. K C MUST NOT BE GREATER THAN MIN(N,P), WHERE P IS THE C SAME AS IN THE CALLING SEQUENCE TO DQRDC. C QRAUX DOUBLE PRECISION(P). C QRAUX CONTAINS THE AUXILIARY OUTPUT FROM DQRDC. C Y DOUBLE PRECISION(N) C Y CONTAINS AN N-VECTOR THAT IS TO BE MANIPULATED C BY DQRSL. C JOB INTEGER. C JOB SPECIFIES WHAT IS TO BE COMPUTED. JOB HAS C THE DECIMAL EXPANSION ABCDE, WITH THE FOLLOWING C MEANING. 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 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 ON RETURN C QY DOUBLE PRECISION(N). C QY CONTAINS Q*Y, IF ITS COMPUTATION HAS BEEN C REQUESTED. C QTY DOUBLE PRECISION(N). C QTY CONTAINS TRANS(Q)*Y, IF ITS COMPUTATION HAS C BEEN REQUESTED. HERE TRANS(Q) IS THE C TRANSPOSE OF THE MATRIX Q. C B DOUBLE PRECISION(K) C B CONTAINS THE SOLUTION OF THE LEAST SQUARES PROBLEM C MINIMIZE NORM2(Y - XK*B), C IF ITS COMPUTATION HAS BEEN REQUESTED. (NOTE THAT C IF PIVOTING WAS REQUESTED IN DQRDC, THE J-TH C COMPONENT OF B WILL BE ASSOCIATED WITH COLUMN JPVT(J) C OF THE ORIGINAL MATRIX X THAT WAS INPUT INTO DQRDC.) C RSD DOUBLE PRECISION(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 XB DOUBLE PRECISION(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 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 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 CALL DQRSL(X,LDX,N,K,QRAUX,Y,DUM,Y,B,Y,DUM,110,INFO) 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 CALLING SEQUENCE. C 1. (Y,QTY,B) (RSD) (XB) (QY) C 2. (Y,QTY,RSD) (B) (XB) (QY) C 3. (Y,QTY,XB) (B) (RSD) (QY) C 4. (Y,QY) (QTY,B) (RSD) (XB) C 5. (Y,QY) (QTY,RSD) (B) (XB) C 6. (Y,QY) (QTY,XB) (B) (RSD) C IN ANY GROUP THE VALUE RETURNED IN THE ARRAY ALLOCATED TO C THE GROUP CORRESPONDS TO THE LAST MEMBER OF THE GROUP. C LINPACK. THIS VERSION DATED 08/14/78 . C G. W. STEWART, UNIVERSITY OF MARYLAND, ARGONNE NATIONAL LAB. C***REFERENCES DONGARRA J.J., BUNCH J.R., MOLER C.B., STEWART G.W., C *LINPACK USERS GUIDE*, SIAM, 1979. C***ROUTINES CALLED DAXPY,DCOPY,DDOT C***END PROLOGUE DQRSL C...SCALAR ARGUMENTS INTEGER + INFO,JOB,K,LDX,N C...ARRAY ARGUMENTS DOUBLE PRECISION + B(*),QRAUX(*),QTY(*),QY(*),RSD(*),X(LDX,*),XB(*), + Y(*) C...LOCAL SCALARS DOUBLE PRECISION + T,TEMP INTEGER + I,J,JJ,JU,KP1 LOGICAL + CB,CQTY,CQY,CR,CXB C...EXTERNAL FUNCTIONS DOUBLE PRECISION + DDOT EXTERNAL + DDOT C...EXTERNAL SUBROUTINES EXTERNAL + DAXPY,DCOPY C...INTRINSIC FUNCTIONS INTRINSIC + MIN0,MOD C***FIRST EXECUTABLE STATEMENT DQRSL INFO = 0 C DETERMINE WHAT IS TO BE COMPUTED. 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 SPECIAL ACTION WHEN N=1. IF (JU .NE. 0) GO TO 40 IF (CQY) QY(1) = Y(1) IF (CQTY) QTY(1) = Y(1) IF (CXB) XB(1) = Y(1) IF (.NOT.CB) GO TO 30 IF (X(1,1) .NE. 0.0D0) GO TO 10 INFO = 1 GO TO 20 10 CONTINUE B(1) = Y(1)/X(1,1) 20 CONTINUE 30 CONTINUE IF (CR) RSD(1) = 0.0D0 GO TO 250 40 CONTINUE C SET UP TO COMPUTE QY OR QTY. IF (CQY) CALL DCOPY(N,Y,1,QY,1) IF (CQTY) CALL DCOPY(N,Y,1,QTY,1) IF (.NOT.CQY) GO TO 70 C COMPUTE QY. DO 60 JJ = 1, JU J = JU - JJ + 1 IF (QRAUX(J) .EQ. 0.0D0) GO TO 50 TEMP = X(J,J) X(J,J) = QRAUX(J) T = -DDOT(N-J+1,X(J,J),1,QY(J),1)/X(J,J) CALL DAXPY(N-J+1,T,X(J,J),1,QY(J),1) X(J,J) = TEMP 50 CONTINUE 60 CONTINUE 70 CONTINUE IF (.NOT.CQTY) GO TO 100 C COMPUTE TRANS(Q)*Y. DO 90 J = 1, JU IF (QRAUX(J) .EQ. 0.0D0) GO TO 80 TEMP = X(J,J) X(J,J) = QRAUX(J) T = -DDOT(N-J+1,X(J,J),1,QTY(J),1)/X(J,J) CALL DAXPY(N-J+1,T,X(J,J),1,QTY(J),1) X(J,J) = TEMP 80 CONTINUE 90 CONTINUE 100 CONTINUE C SET UP TO COMPUTE B, RSD, OR XB. IF (CB) CALL DCOPY(K,QTY,1,B,1) KP1 = K + 1 IF (CXB) CALL DCOPY(K,QTY,1,XB,1) IF (CR .AND. K .LT. N) CALL DCOPY(N-K,QTY(KP1),1,RSD(KP1),1) IF (.NOT.CXB .OR. KP1 .GT. N) GO TO 120 DO 110 I = KP1, N XB(I) = 0.0D0 110 CONTINUE 120 CONTINUE IF (.NOT.CR) GO TO 140 DO 130 I = 1, K RSD(I) = 0.0D0 130 CONTINUE 140 CONTINUE IF (.NOT.CB) GO TO 190 C COMPUTE B. DO 170 JJ = 1, K J = K - JJ + 1 IF (X(J,J) .NE. 0.0D0) GO TO 150 INFO = J C ......EXIT GO TO 180 150 CONTINUE B(J) = B(J)/X(J,J) IF (J .EQ. 1) GO TO 160 T = -B(J) CALL DAXPY(J-1,T,X(1,J),1,B,1) 160 CONTINUE 170 CONTINUE 180 CONTINUE 190 CONTINUE IF (.NOT.CR .AND. .NOT.CXB) GO TO 240 C COMPUTE RSD OR XB AS REQUIRED. DO 230 JJ = 1, JU J = JU - JJ + 1 IF (QRAUX(J) .EQ. 0.0D0) GO TO 220 TEMP = X(J,J) X(J,J) = QRAUX(J) IF (.NOT.CR) GO TO 200 T = -DDOT(N-J+1,X(J,J),1,RSD(J),1)/X(J,J) CALL DAXPY(N-J+1,T,X(J,J),1,RSD(J),1) 200 CONTINUE IF (.NOT.CXB) GO TO 210 T = -DDOT(N-J+1,X(J,J),1,XB(J),1)/X(J,J) CALL DAXPY(N-J+1,T,X(J,J),1,XB(J),1) 210 CONTINUE X(J,J) = TEMP 220 CONTINUE 230 CONTINUE 240 CONTINUE 250 CONTINUE RETURN END *DROT SUBROUTINE DROT(N,DX,INCX,DY,INCY,DC,DS) C***BEGIN PROLOGUE DROT C***DATE WRITTEN 791001 (YYMMDD) C***REVISION DATE 820801 (YYMMDD) C***CATEGORY NO. D1A8 C***KEYWORDS BLAS,GIVENS ROTATION,LINEAR ALGEBRA,VECTOR C***AUTHOR LAWSON, C. L., (JPL) C HANSON, R. J., (SNLA) C KINCAID, D. R., (U. OF TEXAS) C KROGH, F. T., (JPL) C***PURPOSE APPLY D.P. GIVENS ROTATION C***DESCRIPTION C B L A S SUBPROGRAM C DESCRIPTION OF PARAMETERS C --INPUT-- C N NUMBER OF ELEMENTS IN INPUT VECTOR(S) C DX DOUBLE PRECISION VECTOR WITH N ELEMENTS C INCX STORAGE SPACING BETWEEN ELEMENTS OF DX C DY DOUBLE PRECISION VECTOR WITH N ELEMENTS C INCY STORAGE SPACING BETWEEN ELEMENTS OF DY C DC D.P. ELEMENT OF ROTATION MATRIX C DS D.P. ELEMENT OF ROTATION MATRIX C --OUTPUT-- C DX ROTATED VECTOR (UNCHANGED IF N .LE. 0) C DY ROTATED VECTOR (UNCHANGED IF N .LE. 0) C MULTIPLY THE 2 X 2 MATRIX ( DC DS) TIMES THE 2 X N MATRIX (DX**T) C (-DS DC) (DY**T) C WHERE **T INDICATES TRANSPOSE. THE ELEMENTS OF DX ARE IN C DX(LX+I*INCX), I = 0 TO N-1, WHERE LX = 1 IF INCX .GE. 0, ELSE C LX = (-INCX)*N, AND SIMILARLY FOR DY USING LY AND INCY. C***REFERENCES LAWSON C.L., HANSON R.J., KINCAID D.R., KROGH F.T., C *BASIC LINEAR ALGEBRA SUBPROGRAMS FOR FORTRAN USAGE*, C ALGORITHM NO. 539, TRANSACTIONS ON MATHEMATICAL C SOFTWARE, VOLUME 5, NUMBER 3, SEPTEMBER 1979, 308-323 C***ROUTINES CALLED (NONE) C***END PROLOGUE DROT C...SCALAR ARGUMENTS DOUBLE PRECISION + DC,DS INTEGER + INCX,INCY,N C...ARRAY ARGUMENTS DOUBLE PRECISION + DX(*),DY(*) C...LOCAL SCALARS DOUBLE PRECISION + ONE,W,Z,ZERO INTEGER + I,KX,KY,NSTEPS C...DATA STATEMENTS DATA + ZERO,ONE/0.D0,1.D0/ C***FIRST EXECUTABLE STATEMENT DROT IF(N .LE. 0 .OR. (DS .EQ. ZERO .AND. DC .EQ. ONE)) GO TO 40 IF(.NOT. (INCX .EQ. INCY .AND. INCX .GT. 0)) GO TO 20 NSTEPS=INCX*N DO 10 I=1,NSTEPS,INCX W=DX(I) Z=DY(I) DX(I)=DC*W+DS*Z DY(I)=-DS*W+DC*Z 10 CONTINUE GO TO 40 20 CONTINUE KX=1 KY=1 IF(INCX .LT. 0) KX=1-(N-1)*INCX IF(INCY .LT. 0) KY=1-(N-1)*INCY DO 30 I=1,N W=DX(KX) Z=DY(KY) DX(KX)=DC*W+DS*Z DY(KY)=-DS*W+DC*Z KX=KX+INCX KY=KY+INCY 30 CONTINUE 40 CONTINUE RETURN END *DROTG SUBROUTINE DROTG(DA,DB,DC,DS) C***BEGIN PROLOGUE DROTG C***DATE WRITTEN 791001 (YYMMDD) C***REVISION DATE 820801 (YYMMDD) C***CATEGORY NO. D1B10 C***KEYWORDS BLAS,GIVENS ROTATION,LINEAR ALGEBRA,VECTOR C***AUTHOR LAWSON, C. L., (JPL) C HANSON, R. J., (SNLA) C KINCAID, D. R., (U. OF TEXAS) C KROGH, F. T., (JPL) C***PURPOSE CONSTRUCT D.P. PLANE GIVENS ROTATION C***DESCRIPTION C B L A S SUBPROGRAM C DESCRIPTION OF PARAMETERS C --INPUT-- C DA DOUBLE PRECISION SCALAR C DB DOUBLE PRECISION SCALAR C --OUTPUT-- C DA DOUBLE PRECISION RESULT R C DB DOUBLE PRECISION RESULT Z C DC DOUBLE PRECISION RESULT C DS DOUBLE PRECISION RESULT C DESIGNED BY C. L. LAWSON, JPL, 1977 SEPT 08 C CONSTRUCT THE GIVENS TRANSFORMATION C ( DC DS ) C G = ( ) , DC**2 + DS**2 = 1 , C (-DS DC ) C WHICH ZEROS THE SECOND ENTRY OF THE 2-VECTOR (DA,DB)**T . C THE QUANTITY R = (+/-)DSQRT(DA**2 + DB**2) OVERWRITES DA IN C STORAGE. THE VALUE OF DB IS OVERWRITTEN BY A VALUE Z WHICH C ALLOWS DC AND DS TO BE RECOVERED BY THE FOLLOWING ALGORITHM. C IF Z=1 SET DC=0.D0 AND DS=1.D0 C IF DABS(Z) .LT. 1 SET DC=DSQRT(1-Z**2) AND DS=Z C IF DABS(Z) .GT. 1 SET DC=1/Z AND DS=DSQRT(1-DC**2) C NORMALLY, THE SUBPROGRAM DROT(N,DX,INCX,DY,INCY,DC,DS) WILL C NEXT BE CALLED TO APPLY THE TRANSFORMATION TO A 2 BY N MATRIX. C***REFERENCES LAWSON C.L., HANSON R.J., KINCAID D.R., KROGH F.T., C *BASIC LINEAR ALGEBRA SUBPROGRAMS FOR FORTRAN USAGE*, C ALGORITHM NO. 539, TRANSACTIONS ON MATHEMATICAL C SOFTWARE, VOLUME 5, NUMBER 3, SEPTEMBER 1979, 308-323 C***ROUTINES CALLED (NONE) C***END PROLOGUE DROTG C...SCALAR ARGUMENTS DOUBLE PRECISION + DA,DB,DC,DS C...LOCAL SCALARS DOUBLE PRECISION + R,U,V C...INTRINSIC FUNCTIONS INTRINSIC + DABS,DSQRT C***FIRST EXECUTABLE STATEMENT DROTG IF (DABS(DA) .LE. DABS(DB)) GO TO 10 C *** HERE DABS(DA) .GT. DABS(DB) *** U = DA + DA V = DB / U C NOTE THAT U AND R HAVE THE SIGN OF DA R = DSQRT(.25D0 + V**2) * U C NOTE THAT DC IS POSITIVE DC = DA / R DS = V * (DC + DC) DB = DS DA = R RETURN C *** HERE DABS(DA) .LE. DABS(DB) *** 10 IF (DB .EQ. 0.D0) GO TO 20 U = DB + DB V = DA / U C NOTE THAT U AND R HAVE THE SIGN OF DB C (R IS IMMEDIATELY STORED IN DA) DA = DSQRT(.25D0 + V**2) * U C NOTE THAT DS IS POSITIVE DS = DB / DA DC = V * (DS + DS) IF (DC .EQ. 0.D0) GO TO 15 DB = 1.D0 / DC RETURN 15 DB = 1.D0 RETURN C *** HERE DA = DB = 0.D0 *** 20 DC = 1.D0 DS = 0.D0 RETURN END *DSCAL SUBROUTINE DSCAL(N,DA,DX,INCX) C***BEGIN PROLOGUE DSCAL C***DATE WRITTEN 791001 (YYMMDD) C***REVISION DATE 820801 (YYMMDD) C***CATEGORY NO. D1A6 C***KEYWORDS BLAS,LINEAR ALGEBRA,SCALE,VECTOR C***AUTHOR LAWSON, C. L., (JPL) C HANSON, R. J., (SNLA) C KINCAID, D. R., (U. OF TEXAS) C KROGH, F. T., (JPL) C***PURPOSE D.P. VECTOR SCALE X = A*X C***DESCRIPTION C B L A S SUBPROGRAM C DESCRIPTION OF PARAMETERS C --INPUT-- C N NUMBER OF ELEMENTS IN INPUT VECTOR(S) C DA DOUBLE PRECISION SCALE FACTOR C DX DOUBLE PRECISION VECTOR WITH N ELEMENTS C INCX STORAGE SPACING BETWEEN ELEMENTS OF DX C --OUTPUT-- C DX DOUBLE PRECISION RESULT (UNCHANGED IF N.LE.0) C REPLACE DOUBLE PRECISION DX BY DOUBLE PRECISION DA*DX. C FOR I = 0 TO N-1, REPLACE DX(1+I*INCX) WITH DA * DX(1+I*INCX) C***REFERENCES LAWSON C.L., HANSON R.J., KINCAID D.R., KROGH F.T., C *BASIC LINEAR ALGEBRA SUBPROGRAMS FOR FORTRAN USAGE*, C ALGORITHM NO. 539, TRANSACTIONS ON MATHEMATICAL C SOFTWARE, VOLUME 5, NUMBER 3, SEPTEMBER 1979, 308-323 C***ROUTINES CALLED (NONE) C***END PROLOGUE DSCAL C...SCALAR ARGUMENTS DOUBLE PRECISION + DA INTEGER + INCX,N C...ARRAY ARGUMENTS DOUBLE PRECISION + DX(*) C...LOCAL SCALARS INTEGER + I,M,MP1,NS C...INTRINSIC FUNCTIONS INTRINSIC + MOD C***FIRST EXECUTABLE STATEMENT DSCAL IF(N.LE.0)RETURN IF(INCX.EQ.1)GOTO 20 C CODE FOR INCREMENTS NOT EQUAL TO 1. NS = N*INCX DO 10 I = 1,NS,INCX DX(I) = DA*DX(I) 10 CONTINUE RETURN C CODE FOR INCREMENTS EQUAL TO 1. C CLEAN-UP LOOP SO REMAINING VECTOR LENGTH IS A MULTIPLE OF 5. 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 *DSWAP SUBROUTINE DSWAP(N,DX,INCX,DY,INCY) C***BEGIN PROLOGUE DSWAP C***DATE WRITTEN 791001 (YYMMDD) C***REVISION DATE 820801 (YYMMDD) C***CATEGORY NO. D1A5 C***KEYWORDS BLAS,DOUBLE PRECISION,INTERCHANGE,LINEAR ALGEBRA,VECTOR C***AUTHOR LAWSON, C. L., (JPL) C HANSON, R. J., (SNLA) C KINCAID, D. R., (U. OF TEXAS) C KROGH, F. T., (JPL) C***PURPOSE INTERCHANGE D.P. VECTORS C***DESCRIPTION C B L A S SUBPROGRAM C DESCRIPTION OF PARAMETERS C --INPUT-- C N NUMBER OF ELEMENTS IN INPUT VECTOR(S) C DX DOUBLE PRECISION VECTOR WITH N ELEMENTS C INCX STORAGE SPACING BETWEEN ELEMENTS OF DX C DY DOUBLE PRECISION VECTOR WITH N ELEMENTS C INCY STORAGE SPACING BETWEEN ELEMENTS OF DY C --OUTPUT-- C DX INPUT VECTOR DY (UNCHANGED IF N .LE. 0) C DY INPUT VECTOR DX (UNCHANGED IF N .LE. 0) C INTERCHANGE DOUBLE PRECISION DX AND DOUBLE PRECISION DY. C FOR I = 0 TO N-1, INTERCHANGE DX(LX+I*INCX) AND DY(LY+I*INCY), C WHERE LX = 1 IF INCX .GE. 0, ELSE LX = (-INCX)*N, AND LY IS C DEFINED IN A SIMILAR WAY USING INCY. C***REFERENCES LAWSON C.L., HANSON R.J., KINCAID D.R., KROGH F.T., C *BASIC LINEAR ALGEBRA SUBPROGRAMS FOR FORTRAN USAGE*, C ALGORITHM NO. 539, TRANSACTIONS ON MATHEMATICAL C SOFTWARE, VOLUME 5, NUMBER 3, SEPTEMBER 1979, 308-323 C***ROUTINES CALLED (NONE) C***END PROLOGUE DSWAP C...SCALAR ARGUMENTS INTEGER + INCX,INCY,N C...ARRAY ARGUMENTS DOUBLE PRECISION + DX(*),DY(*) C...LOCAL SCALARS DOUBLE PRECISION + DTEMP1,DTEMP2,DTEMP3 INTEGER + I,IX,IY,M,MP1,NS C...INTRINSIC FUNCTIONS INTRINSIC + MOD C***FIRST EXECUTABLE STATEMENT DSWAP IF(N.LE.0)RETURN IF(INCX.EQ.INCY) IF(INCX-1) 5,20,60 5 CONTINUE C CODE FOR UNEQUAL OR NONPOSITIVE INCREMENTS. 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 DTEMP1 = DX(IX) DX(IX) = DY(IY) DY(IY) = DTEMP1 IX = IX + INCX IY = IY + INCY 10 CONTINUE RETURN C CODE FOR BOTH INCREMENTS EQUAL TO 1 C CLEAN-UP LOOP SO REMAINING VECTOR LENGTH IS A MULTIPLE OF 3. 20 M = MOD(N,3) IF( M .EQ. 0 ) GO TO 40 DO 30 I = 1,M DTEMP1 = DX(I) DX(I) = DY(I) DY(I) = DTEMP1 30 CONTINUE IF( N .LT. 3 ) RETURN 40 MP1 = M + 1 DO 50 I = MP1,N,3 DTEMP1 = DX(I) DTEMP2 = DX(I+1) DTEMP3 = DX(I+2) DX(I) = DY(I) DX(I+1) = DY(I+1) DX(I+2) = DY(I+2) DY(I) = DTEMP1 DY(I+1) = DTEMP2 DY(I+2) = DTEMP3 50 CONTINUE RETURN 60 CONTINUE C CODE FOR EQUAL, POSITIVE, NONUNIT INCREMENTS. NS = N*INCX DO 70 I=1,NS,INCX DTEMP1 = DX(I) DX(I) = DY(I) DY(I) = DTEMP1 70 CONTINUE RETURN END *DTRCO SUBROUTINE DTRCO(T,LDT,N,RCOND,Z,JOB) C***BEGIN PROLOGUE DTRCO C***DATE WRITTEN 780814 (YYMMDD) C***REVISION DATE 820801 (YYMMDD) C***CATEGORY NO. D2A3 C***KEYWORDS CONDITION,DOUBLE PRECISION,FACTOR,LINEAR ALGEBRA,LINPACK, C MATRIX,TRIANGULAR C***AUTHOR MOLER, C. B., (U. OF NEW MEXICO) C***PURPOSE ESTIMATES THE CONDITION OF A DOUBLE PRECISION TRIANGULAR C MATRIX. C***DESCRIPTION C DTRCO ESTIMATES THE CONDITION OF A DOUBLE PRECISION TRIANGULAR C MATRIX. C ON ENTRY C T DOUBLE PRECISION(LDT,N) C T CONTAINS THE TRIANGULAR MATRIX. THE ZERO C ELEMENTS OF THE MATRIX ARE NOT REFERENCED, AND C THE CORRESPONDING ELEMENTS OF THE ARRAY CAN BE C USED TO STORE OTHER INFORMATION. C LDT INTEGER C LDT IS THE LEADING DIMENSION OF THE ARRAY T. C N INTEGER C N IS THE ORDER OF THE SYSTEM. C JOB INTEGER C = 0 T IS LOWER TRIANGULAR. C = NONZERO T IS UPPER TRIANGULAR. C ON RETURN C RCOND DOUBLE PRECISION C AN ESTIMATE OF THE RECIPROCAL CONDITION OF T . C FOR THE SYSTEM T*X = B , RELATIVE PERTURBATIONS C IN T 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 T MAY BE SINGULAR TO WORKING C PRECISION. IN PARTICULAR, RCOND IS ZERO IF C EXACT SINGULARITY IS DETECTED OR THE ESTIMATE C UNDERFLOWS. C Z DOUBLE PRECISION(N) C A WORK VECTOR WHOSE CONTENTS ARE USUALLY UNIMPORTANT. C IF T 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 LINPACK. THIS VERSION DATED 08/14/78 . C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB. C***REFERENCES DONGARRA J.J., BUNCH J.R., MOLER C.B., STEWART G.W., C *LINPACK USERS GUIDE*, SIAM, 1979. C***ROUTINES CALLED DASUM,DAXPY,DSCAL C***END PROLOGUE DTRCO C...SCALAR ARGUMENTS DOUBLE PRECISION + RCOND INTEGER + JOB,LDT,N C...ARRAY ARGUMENTS DOUBLE PRECISION + T(LDT,*),Z(*) C...LOCAL SCALARS DOUBLE PRECISION + EK,S,SM,TNORM,W,WK,WKM,YNORM INTEGER + I1,J,J1,J2,K,KK,L LOGICAL + LOWER C...EXTERNAL FUNCTIONS DOUBLE PRECISION + DASUM EXTERNAL + DASUM C...EXTERNAL SUBROUTINES EXTERNAL + DAXPY,DSCAL C...INTRINSIC FUNCTIONS INTRINSIC + DABS,DMAX1,DSIGN C***FIRST EXECUTABLE STATEMENT DTRCO LOWER = JOB .EQ. 0 C COMPUTE 1-NORM OF T TNORM = 0.0D0 DO 10 J = 1, N L = J IF (LOWER) L = N + 1 - J I1 = 1 IF (LOWER) I1 = J TNORM = DMAX1(TNORM,DASUM(L,T(I1,J),1)) 10 CONTINUE C RCOND = 1/(NORM(T)*(ESTIMATE OF NORM(INVERSE(T)))) . C ESTIMATE = NORM(Z)/NORM(Y) WHERE T*Z = Y AND TRANS(T)*Y = E . C TRANS(T) IS THE TRANSPOSE OF T . C THE COMPONENTS OF E ARE CHOSEN TO CAUSE MAXIMUM LOCAL C GROWTH IN THE ELEMENTS OF Y . C THE VECTORS ARE FREQUENTLY RESCALED TO AVOID OVERFLOW. C SOLVE TRANS(T)*Y = E EK = 1.0D0 DO 20 J = 1, N Z(J) = 0.0D0 20 CONTINUE DO 100 KK = 1, N K = KK IF (LOWER) K = N + 1 - KK IF (Z(K) .NE. 0.0D0) EK = DSIGN(EK,-Z(K)) IF (DABS(EK-Z(K)) .LE. DABS(T(K,K))) GO TO 30 S = DABS(T(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 (T(K,K) .EQ. 0.0D0) GO TO 40 WK = WK/T(K,K) WKM = WKM/T(K,K) GO TO 50 40 CONTINUE WK = 1.0D0 WKM = 1.0D0 50 CONTINUE IF (KK .EQ. N) GO TO 90 J1 = K + 1 IF (LOWER) J1 = 1 J2 = N IF (LOWER) J2 = K - 1 DO 60 J = J1, J2 SM = SM + DABS(Z(J)+WKM*T(K,J)) Z(J) = Z(J) + WK*T(K,J) S = S + DABS(Z(J)) 60 CONTINUE IF (S .GE. SM) GO TO 80 W = WKM - WK WK = WKM DO 70 J = J1, J2 Z(J) = Z(J) + W*T(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) YNORM = 1.0D0 C SOLVE T*Z = Y DO 130 KK = 1, N K = N + 1 - KK IF (LOWER) K = KK IF (DABS(Z(K)) .LE. DABS(T(K,K))) GO TO 110 S = DABS(T(K,K))/DABS(Z(K)) CALL DSCAL(N,S,Z,1) YNORM = S*YNORM 110 CONTINUE IF (T(K,K) .NE. 0.0D0) Z(K) = Z(K)/T(K,K) IF (T(K,K) .EQ. 0.0D0) Z(K) = 1.0D0 I1 = 1 IF (LOWER) I1 = K + 1 IF (KK .GE. N) GO TO 120 W = -Z(K) CALL DAXPY(N-KK,W,T(I1,K),1,Z(I1),1) 120 CONTINUE 130 CONTINUE C MAKE ZNORM = 1.0 S = 1.0D0/DASUM(N,Z,1) CALL DSCAL(N,S,Z,1) YNORM = S*YNORM IF (TNORM .NE. 0.0D0) RCOND = YNORM/TNORM IF (TNORM .EQ. 0.0D0) RCOND = 0.0D0 RETURN END *DTRSL SUBROUTINE DTRSL(T,LDT,N,B,JOB,INFO) C***BEGIN PROLOGUE DTRSL C***DATE WRITTEN 780814 (YYMMDD) C***REVISION DATE 820801 (YYMMDD) C***CATEGORY NO. D2A3 C***KEYWORDS DOUBLE PRECISION,LINEAR ALGEBRA,LINPACK,MATRIX,SOLVE, C TRIANGULAR C***AUTHOR STEWART, G. W., (U. OF MARYLAND) C***PURPOSE SOLVES SYSTEMS OF THE FORM T*X=B OR TRANS(T)*X=B WHERE T C IS A TRIANGULAR MATRIX OF ORDER N. C***DESCRIPTION C DTRSL SOLVES SYSTEMS OF THE FORM C T * X = B C OR C TRANS(T) * X = B C WHERE T IS A TRIANGULAR MATRIX OF ORDER N. HERE TRANS(T) C DENOTES THE TRANSPOSE OF THE MATRIX T. C ON ENTRY C T DOUBLE PRECISION(LDT,N) C T CONTAINS THE MATRIX OF THE SYSTEM. THE ZERO C ELEMENTS OF THE MATRIX ARE NOT REFERENCED, AND C THE CORRESPONDING ELEMENTS OF THE ARRAY CAN BE C USED TO STORE OTHER INFORMATION. C LDT INTEGER C LDT IS THE LEADING DIMENSION OF THE ARRAY T. C N INTEGER C N IS THE ORDER OF THE SYSTEM. C B DOUBLE PRECISION(N). C B CONTAINS THE RIGHT HAND SIDE OF THE SYSTEM. C JOB INTEGER C JOB SPECIFIES WHAT KIND OF SYSTEM IS TO BE SOLVED. C IF JOB IS C 00 SOLVE T*X=B, T LOWER TRIANGULAR, C 01 SOLVE T*X=B, T UPPER TRIANGULAR, C 10 SOLVE TRANS(T)*X=B, T LOWER TRIANGULAR, C 11 SOLVE TRANS(T)*X=B, T UPPER TRIANGULAR. C ON RETURN C B B CONTAINS THE SOLUTION, IF INFO .EQ. 0. C OTHERWISE B IS UNALTERED. C INFO INTEGER C INFO CONTAINS ZERO IF THE SYSTEM IS NONSINGULAR. C OTHERWISE INFO CONTAINS THE INDEX OF C THE FIRST ZERO DIAGONAL ELEMENT OF T. C LINPACK. THIS VERSION DATED 08/14/78 . C G. W. STEWART, UNIVERSITY OF MARYLAND, ARGONNE NATIONAL LAB. C***REFERENCES DONGARRA J.J., BUNCH J.R., MOLER C.B., STEWART G.W., C *LINPACK USERS GUIDE*, SIAM, 1979. C***ROUTINES CALLED DAXPY,DDOT C***END PROLOGUE DTRSL C...SCALAR ARGUMENTS INTEGER + INFO,JOB,LDT,N C...ARRAY ARGUMENTS DOUBLE PRECISION + B(*),T(LDT,*) C...LOCAL SCALARS DOUBLE PRECISION + TEMP INTEGER + CASE,J,JJ C...EXTERNAL FUNCTIONS DOUBLE PRECISION + DDOT EXTERNAL + DDOT C...EXTERNAL SUBROUTINES EXTERNAL + DAXPY C...INTRINSIC FUNCTIONS INTRINSIC + MOD C***FIRST EXECUTABLE STATEMENT DTRSL C BEGIN BLOCK PERMITTING ...EXITS TO 150 C CHECK FOR ZERO DIAGONAL ELEMENTS. DO 10 INFO = 1, N C ......EXIT IF (T(INFO,INFO) .EQ. 0.0D0) GO TO 150 10 CONTINUE INFO = 0 C DETERMINE THE TASK AND GO TO IT. CASE = 1 IF (MOD(JOB,10) .NE. 0) CASE = 2 IF (MOD(JOB,100)/10 .NE. 0) CASE = CASE + 2 GO TO (20,50,80,110), CASE C SOLVE T*X=B FOR T LOWER TRIANGULAR 20 CONTINUE B(1) = B(1)/T(1,1) IF (N .LT. 2) GO TO 40 DO 30 J = 2, N TEMP = -B(J-1) CALL DAXPY(N-J+1,TEMP,T(J,J-1),1,B(J),1) B(J) = B(J)/T(J,J) 30 CONTINUE 40 CONTINUE GO TO 140 C SOLVE T*X=B FOR T UPPER TRIANGULAR. 50 CONTINUE B(N) = B(N)/T(N,N) IF (N .LT. 2) GO TO 70 DO 60 JJ = 2, N J = N - JJ + 1 TEMP = -B(J+1) CALL DAXPY(J,TEMP,T(1,J+1),1,B(1),1) B(J) = B(J)/T(J,J) 60 CONTINUE 70 CONTINUE GO TO 140 C SOLVE TRANS(T)*X=B FOR T LOWER TRIANGULAR. 80 CONTINUE B(N) = B(N)/T(N,N) IF (N .LT. 2) GO TO 100 DO 90 JJ = 2, N J = N - JJ + 1 B(J) = B(J) - DDOT(JJ-1,T(J+1,J),1,B(J+1),1) B(J) = B(J)/T(J,J) 90 CONTINUE 100 CONTINUE GO TO 140 C SOLVE TRANS(T)*X=B FOR T UPPER TRIANGULAR. 110 CONTINUE B(1) = B(1)/T(1,1) IF (N .LT. 2) GO TO 130 DO 120 J = 2, N B(J) = B(J) - DDOT(J-1,T(1,J),1,B(1),1) B(J) = B(J)/T(J,J) 120 CONTINUE 130 CONTINUE 140 CONTINUE 150 CONTINUE RETURN END *IDAMAX INTEGER FUNCTION IDAMAX(N,DX,INCX) C***BEGIN PROLOGUE IDAMAX C***DATE WRITTEN 791001 (YYMMDD) C***REVISION DATE 820801 (YYMMDD) C***CATEGORY NO. D1A2 C***KEYWORDS BLAS,DOUBLE PRECISION,LINEAR ALGEBRA,MAXIMUM COMPONENT, C VECTOR C***AUTHOR LAWSON, C. L., (JPL) C HANSON, R. J., (SNLA) C KINCAID, D. R., (U. OF TEXAS) C KROGH, F. T., (JPL) C***PURPOSE FIND LARGEST COMPONENT OF D.P. VECTOR C***DESCRIPTION C B L A S SUBPROGRAM C DESCRIPTION OF PARAMETERS C --INPUT-- C N NUMBER OF ELEMENTS IN INPUT VECTOR(S) C DX DOUBLE PRECISION VECTOR WITH N ELEMENTS C INCX STORAGE SPACING BETWEEN ELEMENTS OF DX C --OUTPUT-- C IDAMAX SMALLEST INDEX (ZERO IF N .LE. 0) C FIND SMALLEST INDEX OF MAXIMUM MAGNITUDE OF DOUBLE PRECISION DX. C IDAMAX = FIRST I, I = 1 TO N, TO MINIMIZE ABS(DX(1-INCX+I*INCX) C***REFERENCES LAWSON C.L., HANSON R.J., KINCAID D.R., KROGH F.T., C *BASIC LINEAR ALGEBRA SUBPROGRAMS FOR FORTRAN USAGE*, C ALGORITHM NO. 539, TRANSACTIONS ON MATHEMATICAL C SOFTWARE, VOLUME 5, NUMBER 3, SEPTEMBER 1979, 308-323 C***ROUTINES CALLED (NONE) C***END PROLOGUE IDAMAX C...SCALAR ARGUMENTS INTEGER + INCX,N C...ARRAY ARGUMENTS DOUBLE PRECISION + DX(*) C...LOCAL SCALARS DOUBLE PRECISION + DMAX,XMAG INTEGER + I,II,NS C...INTRINSIC FUNCTIONS INTRINSIC + DABS C***FIRST EXECUTABLE STATEMENT IDAMAX IDAMAX = 0 IF(N.LE.0) RETURN IDAMAX = 1 IF(N.LE.1)RETURN IF(INCX.EQ.1)GOTO 20 C CODE FOR INCREMENTS NOT EQUAL TO 1. DMAX = DABS(DX(1)) NS = N*INCX II = 1 DO 10 I = 1,NS,INCX XMAG = DABS(DX(I)) IF(XMAG.LE.DMAX) GO TO 5 IDAMAX = II DMAX = XMAG 5 II = II + 1 10 CONTINUE RETURN C CODE FOR INCREMENTS EQUAL TO 1. 20 DMAX = DABS(DX(1)) DO 30 I = 2,N XMAG = DABS(DX(I)) IF(XMAG.LE.DMAX) GO TO 30 IDAMAX = I DMAX = XMAG 30 CONTINUE RETURN END