*DECK LSSODS SUBROUTINE LSSODS (A, X, B, M, N, NRDA, IFLAG, IRANK, ISCALE, Q, + DIAG, KPIVOT, ITER, RESNRM, XNORM, Z, R, DIV, TD, SCALES) C***BEGIN PROLOGUE LSSODS C***SUBSIDIARY C***PURPOSE Subsidiary to BVSUP C***LIBRARY SLATEC C***TYPE SINGLE PRECISION (LSSODS-S) C***AUTHOR (UNKNOWN) C***DESCRIPTION C C LSSODS solves the same problem as SODS (in fact, it is called by C SODS) but is somewhat more flexible in its use. In particular, C LSSODS allows for iterative refinement of the solution, makes the C transformation and triangular reduction information more C accessible, and enables the user to avoid destruction of the C original matrix A. C C Modeled after the ALGOL codes in the articles in the REFERENCES C section. C C ********************************************************************** C INPUT C ********************************************************************** C C A -- Contains the matrix of M equations in N unknowns and must C be dimensioned NRDA by N. A remains unchanged C X -- Solution array of length at least N C B -- Given constant vector of length M, B remains unchanged C M -- Number of equations, M greater or equal to 1 C N -- Number of unknowns, N not larger than M C NRDA -- Row dimension of A, NRDA greater or equal to M C IFLAG -- Status indicator C = 0 for the first call (and for each new problem defined by C a new matrix A) when the matrix data is treated as exact C =-K for the first call (and for each new problem defined by C a new matrix A) when the matrix data is assumed to be C accurate to about K digits C = 1 for subsequent calls whenever the matrix A has already C been decomposed (problems with new vectors B but C same matrix a can be handled efficiently) C ISCALE -- Scaling indicator C =-1 if the matrix A is to be pre-scaled by C columns when appropriate C If the scaling indicator is not equal to -1 C no scaling will be attempted C For most problems scaling will probably not be necessary C ITER -- Maximum number of iterative improvement steps to be C performed, 0 .LE. ITER .LE. 10 (SODS uses ITER=0) C Q -- Matrix used for the transformation, must be dimensioned C NRDA by N (SODS puts A in the Q location which conserves C storage but destroys A) C When iterative improvement of the solution is requested, C ITER .GT. 0, this additional storage for Q must be C made available C DIAG,KPIVOT,Z,R, -- Arrays of length N (except for R which is M) C DIV,TD,SCALES used for internal storage C C ********************************************************************** C OUTPUT C ********************************************************************** C C IFLAG -- Status indicator C =1 if solution was obtained C =2 if improper input is detected C =3 if rank of matrix is less than N C if the minimal length least squares solution is C desired, simply reset IFLAG=1 and call the code again C C The next three IFLAG values can occur only when C the iterative improvement mode is being used. C =4 if the problem is ill-conditioned and maximal C machine accuracy is not achievable C =5 if the problem is very ill-conditioned and the solution C IS likely to have no correct digits C =6 if the allowable number of iterative improvement steps C has been completed without getting convergence C X -- Least squares solution of A X = B C IRANK -- Contains the numerically determined matrix rank C the user must not alter this value on succeeding calls C with input values of IFLAG=1 C Q -- Contains the strictly upper triangular part of the reduced C matrix and the transformation information in the lower C triangular part C DIAG -- Contains the diagonal elements of the triangular reduced C matrix C KPIVOT -- Contains the pivotal information. The column interchanges C performed on the original matrix are recorded here C ITER -- The actual number of iterative corrections used C RESNRM -- The Euclidean norm of the residual vector B - A X C XNORM -- The Euclidean norm of the solution vector C DIV,TD -- Contains transformation information for rank C deficient problems C SCALES -- Contains the column scaling parameters C C ********************************************************************** C C***SEE ALSO BVSUP C***REFERENCES G. Golub, Numerical methods for solving linear least C squares problems, Numerische Mathematik 7, (1965), C pp. 206-216. C P. Businger and G. Golub, Linear least squares C solutions by Householder transformations, Numerische C Mathematik 7, (1965), pp. 269-276. C***ROUTINES CALLED J4SAVE, OHTROR, ORTHOL, R1MACH, SDOT, SDSDOT, C XERMAX, XERMSG, XGETF, XSETF C***REVISION HISTORY (YYMMDD) C 750601 DATE WRITTEN C 890531 Changed all specific intrinsics to generic. (WRB) C 890831 Modified array declarations. (WRB) C 891214 Prologue converted to Version 4.0 format. (BAB) C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) C 900402 Added TYPE section. (WRB) C 910408 Updated the REFERENCES section. (WRB) C 920501 Reformatted the REFERENCES section. (WRB) C***END PROLOGUE LSSODS DIMENSION A(NRDA,*),X(*),B(*),Q(NRDA,*),DIAG(*), 1 Z(*),KPIVOT(*),R(*),DIV(*),TD(*),SCALES(*) C C ********************************************************************** C C MACHINE PRECISION (COMPUTER UNIT ROUNDOFF VALUE) IS DEFINED C THE FUNCTION R1MACH. C C***FIRST EXECUTABLE STATEMENT LSSODS URO = R1MACH(3) C C ********************************************************************** C IF (N .LT. 1 .OR. M .LT. N .OR. NRDA .LT. M) GO TO 1 IF (ITER .LT. 0) GO TO 1 IF (IFLAG .LE. 0) GO TO 5 IF (IFLAG .EQ. 1) GO TO 15 C C INVALID INPUT FOR LSSODS 1 IFLAG=2 CALL XERMSG ('SLATEC', 'LSSODS', 'INVALID INPUT PARAMETERS.', 2, + 1) RETURN C 5 CALL XGETF (NFATAL) MAXMES = J4SAVE (4,0,.FALSE.) IF (IFLAG .EQ. 0) GO TO 7 NFAT = -1 IF(NFATAL .EQ. 0) NFAT=0 CALL XSETF (NFAT) CALL XERMAX (1) C C COPY MATRIX A INTO MATRIX Q C 7 DO 10 J=1,N DO 10 K=1,M 10 Q(K,J)=A(K,J) C C USE ORTHOGONAL TRANSFORMATIONS TO REDUCE Q TO C UPPER TRIANGULAR FORM C CALL ORTHOL(Q,M,N,NRDA,IFLAG,IRANK,ISCALE,DIAG,KPIVOT,SCALES,Z,TD) C CALL XSETF (NFATAL) CALL XERMAX (MAXMES) IF (IRANK .EQ. N) GO TO 12 C C FOR RANK DEFICIENT PROBLEMS USE ADDITIONAL ORTHOGONAL C TRANSFORMATIONS TO FURTHER REDUCE Q C IF (IRANK .NE. 0) CALL OHTROR(Q,N,NRDA,DIAG,IRANK,DIV,TD) RETURN C C STORE DIVISORS FOR THE TRIANGULAR SOLUTION C 12 DO 13 K=1,N 13 DIV(K)=DIAG(K) C 15 IRM=IRANK-1 IRP=IRANK+1 ITERP=MIN(ITER+1,11) ACC=10.*URO C C ZERO OUT SOLUTION ARRAY C DO 20 K=1,N 20 X(K)=0. C IF (IRANK .GT. 0) GO TO 25 C C SPECIAL CASE FOR THE NULL MATRIX ITER=0 XNORM=0. RESNRM=SQRT(SDOT(M,B(1),1,B(1),1)) RETURN C C COPY CONSTANT VECTOR INTO R C 25 DO 30 K=1,M 30 R(K)=B(K) C C ********************************************************************** C SOLUTION SECTION C ITERATIVE REFINEMENT OF THE RESIDUAL VECTOR C ********************************************************************** C DO 100 IT=1,ITERP ITER=IT-1 C C APPLY ORTHOGONAL TRANSFORMATION TO R C DO 35 J=1,IRANK MJ=M-J+1 GAMMA=SDOT(MJ,Q(J,J),1,R(J),1)/(DIAG(J)*Q(J,J)) DO 35 K=J,M 35 R(K)=R(K)+GAMMA*Q(K,J) C C BACKWARD SUBSTITUTION FOR TRIANGULAR SYSTEM SOLUTION C Z(IRANK)=R(IRANK)/DIV(IRANK) IF (IRM .EQ. 0) GO TO 45 DO 40 L=1,IRM K=IRANK-L KP=K+1 40 Z(K)=(R(K)-SDOT(L,Q(K,KP),NRDA,Z(KP),1))/DIV(K) C 45 IF (IRANK .EQ. N) GO TO 60 C C FOR RANK DEFICIENT PROBLEMS OBTAIN THE C MINIMAL LENGTH SOLUTION C NMIR=N-IRANK DO 50 K=IRP,N 50 Z(K)=0. DO 55 K=1,IRANK GAM=((TD(K)*Z(K))+SDOT(NMIR,Q(K,IRP),NRDA,Z(IRP),1))/ 1 (TD(K)*DIV(K)) Z(K)=Z(K)+GAM*TD(K) DO 55 J=IRP,N 55 Z(J)=Z(J)+GAM*Q(K,J) C C REORDER SOLUTION COMPONENTS ACCORDING TO PIVOTAL POINTS C AND RESCALE ANSWERS AS DICTATED C 60 DO 65 K=1,N Z(K)=Z(K)*SCALES(K) L=KPIVOT(K) 65 X(L)=X(L)+Z(K) C C COMPUTE CORRECTION VECTOR NORM (SOLUTION NORM) C ZNORM=SQRT(SDOT(N,Z(1),1,Z(1),1)) IF (IT .EQ. 1) XNORM=ZNORM IF (ITERP .GT. 1) GO TO 80 C C NO ITERATIVE CORRECTIONS TO BE PERFORMED, SO COMPUTE C THE APPROXIMATE RESIDUAL NORM DEFINED BY THE EQUATIONS C WHICH ARE NOT SATISFIED BY THE SOLUTION C THEN WE ARE DONE C MMIR=M-IRANK IF (MMIR .EQ. 0) GO TO 70 RESNRM=SQRT(SDOT(MMIR,R(IRP),1,R(IRP),1)) RETURN 70 RESNRM=0. RETURN C C COMPUTE RESIDUAL VECTOR FOR THE ITERATIVE IMPROVEMENT PROCESS C 80 DO 85 K=1,M 85 R(K)=-SDSDOT(N,-B(K),A(K,1),NRDA,X(1),1) RESNRM=SQRT(SDOT(M,R(1),1,R(1),1)) IF (IT .EQ. 1) GO TO 100 C C TEST FOR CONVERGENCE C IF (ZNORM .LE. ACC*XNORM) RETURN C C COMPARE SUCCESSIVE REFINEMENT VECTOR NORMS C FOR LOOP TERMINATION CRITERIA C IF (ZNORM .LE. 0.25*ZNRM0) GO TO 100 IF (IT .EQ. 2) GO TO 90 C IFLAG=4 CALL XERMSG ('SLATEC', 'LSSODS', + 'PROBLEM MAY BE ILL-CONDITIONED. MAXIMAL MACHINE ACCURACY ' // + 'IS NOT ACHIEVABLE.', 3, 1) RETURN C 90 IFLAG=5 CALL XERMSG ('SLATEC', 'LSSODS', + 'PROBLEM IS VERY ILL-CONDITIONED. ITERATIVE ' // + 'IMPROVEMENT IS INEFFECTIVE.', 8, 1) RETURN C 100 ZNRM0=ZNORM C ********************************************************************** C C ********************************************************************** IFLAG=6 CALL XERMSG ('SLATEC', 'LSSODS', + 'CONVERGENCE HAS NOT BEEN OBTAINED WITH ALLOWABLE ' // + 'NUMBER OF ITERATIVE IMPROVEMENT STEPS.', 8, 1) C RETURN END