SUBROUTINE CL1(NEQNS,NEQC,NIQC,NVARS,NACT,IFL,MXS,PSW, * E,NER,X,F,EL1N,RES,INDX,W) C INTEGER IFL,INDX(1),MXS,NACT,NEQC,NEQNS,NER,NIQC,NVARS LOGICAL PSW REAL E(NER,1),EL1N,F(1),RES(1),W(1),X(1) C C C ********************************************************** C A PROGRAM FOR THE SOLUTION IN THE L1 SENSE C OF A LINEAR-EQUATION SYSTEM C (WITH OR WITHOUT LINEAR CONSTRAINTS). C RICHARD H. BARTELS AND ANDREW R. CONN. C LATEST UPDATE .... 13 APRIL, 1980. C C DEVELOPMENT OF THIS PROGRAM WAS SUPPORTED C IN PART BY U. S. NSF GRANT DCR75-07817 C AND BY FUNDS FROM THE NATIONAL BUREAU OF STANDARDS, C AND IN PART BY CANADIAN NRC GRANT A8639. C C +++++ PARAMETERS +++++ C ---------------------------------------------------------- C INPUT C NAME TYPE SUBSCRPT OUTPUT DESCRIPTION C SCRATCH C .......................................................... C NEQNS INT. NONE IN NUMBER OF EQUATIONS C (MAY BE ZERO) C C NEQC INT. NONE IN NUMBER OF EQUALITY C CONSTRAINTS C (MAY BE ZERO) C C NIQC INT. NONE IN NUMBER OF INEQUALITY C CONSTRAINTS C (MAY BE ZERO) C C NVARS INT. NONE IN NUMBER OF VARIABLES C C NACT INT. NONE OUT NUMBER OF ACTIVE C EQUATIONS/CONSTRAINTS C AT TERMINATION C (IF ANY, THEIR ASSOCIATED C COLUMN POSITIONS IN E WILL C BE LISTED IN INDX(1) C THROUGH INDX(NACT) ) C C IFL INT. NONE OUT TERMINATION CODE C (SEE BELOW) C C MXS INT. NONE IN MAXIMUM NUMBER OF STEPS C ALLOWED C C PSW LOGIC. NONE IN PRINT SWITCH C (SEE BELOW) C C E REAL 2 IN EQUATION/CONSTRAINT MATRIX C THE FIRST NEQNS COLUMNS C (SEE NOTE BELOW) SPECIFY C EQUATIONS, THE REMAINING C COLUMNS (IF ANY) SPECIFY C CONSTRAINTS. C C NER INT. NONE IN ROW DIMENSION OF E C C X REAL 1 IN STARTING VALUES FOR THE C UNKNOWNS (USE ZEROS IF NO C GUESS IS AVAILABLE) C OUT TERMINATION VALUES FOR C THE UNKNOWNS C C F REAL 1 IN EQUATION/CONSTRAINT C RIGHT-HAND SIDES C C EL1N REAL NONE OUT L1 NORM OF EQUATION C RESIDUALS AT TERMINATION C C RES REAL 1 OUT EQUATION/CONSTRAINT C RESIDUALS AT TERMINATION C C INDX INT. 1 OUT INDEX VECTOR USED TO RECORD C THE ORDER IN WHICH THE COLUMNS C OF E ARE BEING PROCESSED C C W REAL 1 SCR. WORKING STORAGE C ---------------------------------------------------------- C C +++++ PURPOSE +++++ C ---------------------------------------------------------- C THIS SUBROUTINE SOLVES THE NEQNS BY NVARS C SYSTEM OF EQUATIONS C C (A-TRANSPOSE) * X == B C C SUBJECT TO THE NEQC CONSTRAINTS C C (G-TRANSPOSE) * X .EQ. H C C AND THE NIQC INEQUALITY CONSTRAINTS C C (C-TRANSPOSE) * X .GE. D C C FOR THE UNKNOWNS X(1),...,X(NVARS). C C THE PROBLEM MUST BE WELL-POSED, NONTRIVIAL C AND OVERDETERMINED IN THE SENSE THAT C C NVARS .GE. 1 C NEQNS .GE. 0 C NEQC .GE. 0 C NIQC .GE. 0 C NEQNS+NEQC+NIQC .GE. NVARS. C C FURTHER, NO COLUMN OF A, G OR C SHOULD BE ZERO. C IF THESE CONDITIONS ARE NOT MET, THE PROGRAM C WILL TERMINATE WITHOUT PERFORMING ANY SUBSTANTIVE C COMPUTATIONS. C C A POINT X IS A SOLUTION IF IT MINIMIZES THE EQUATION C RESIDUALS FROM AMONG ALL POINTS WHICH SATISFY THE C CONSTRAINTS. AT ANY (NONDEGENERATE) SOLUTION C THERE WILL BE NACT EQUATIONS AND CONSTRAINTS C WHOSE RESIDUALS C C (A(I)-TRANSPOSE) * X - B(I) C C (G(I)-TRANSPOSE) * X - H(I) C C AND C C (C(I)-TRANSPOSE) * X - D(I) C C ARE ZERO. C C THE COLUMNS OF (A,G,C) CORRESPONDING TO THE ZERO RESIDUALS C ARE REFERRED TO AS ACTIVE COLUMNS THROUGHOUT THIS LISTING. C THE NUMBERS OF THE ACTIVE COLUMNS ARE MAINTAINED AS THE C ENTRIES 1,...,NACT OF THE ARRAY INDX. C C A SOLUTION X IS FOUND BY MINIMIZING A PIECEWISE C LINEAR PENALTY FUNCTION FORMED FROM THE L1 C NORM OF THE EQUATION RESIDUALS AND THE SUM OF THE C INFEASIBILITIES IN THE CONSTRAINTS. C THE MINIMIZATION PROCEEDS IN A STEP-BY-STEP C FASHION, TERMINATING AFTER A FINITE NUMBER OF STEPS. C C NOTE THAT A, G AND C APPEAR TRANSPOSED IN THE C PROBLEM FORMULATION. HENCE IT IS THE COLUMNS OF (A,G,C) C WHICH DEFINE THE EQUATIONS AND CONSTRAINTS RESPECTIVELY. C C THE ARRAY E IS A COMPOSITE OF A, G AND C C AND F IS A COMPOSITE OF B, H AND D. C E SHOULD CONTAIN A AS ITS FIRST NEQNS COLUMNS. C IT SHOULD CONTAIN G AS ITS NEXT NEQC COLUMNS AND C CONTAIN C AS ITS REMAINING NIQC COLUMNS. C SIMILARLY F SHOULD CONTAIN B AS ITS FIRST C NEQNS COMPONENTS, H AS ITS NEXT NEQC COMPONENTS C AND D AS ITS LAST NIQC COMPONENTS. C ---------------------------------------------------------- C C +++++ ARRAYS +++++ C ---------------------------------------------------------- C E IS TO BE DIMENSIONED AT LEAST N BY M, C X AT LEAST N, C F AT LEAST M, C RES AT LEAST M, C INDX AT LEAST M, C W AT LEAST ((3*N*N+11*N+2)/2) + (2*M). C C WHERE N = NVARS AND C M = NEQNS+NEQC+NIQC C ---------------------------------------------------------- C C +++++ INITIALIZATION +++++ C ---------------------------------------------------------- C THE USER MUST INITIALIZE C C NEQNS,NEQC,NIQC,NVARS,MXS,PSW,E,NER,X,F . C C THE FOLLOWING ARE SET BY CL1 C AND DO NOT REQUIRE INITIALIZATION C C NACT,INDX,RES . C C THE ARRAY W IS USED AS SCRATCH SPACE. C ---------------------------------------------------------- C C +++++ TERMINATION CODES AND INTERMEDIATE PRINTING +++++ C ---------------------------------------------------------- C MXS SETS A LIMIT ON THE NUMBER OF MINIMIZATION STEPS TO BE C TAKEN. C C UPON TERMINATION IFL WILL BE SET ACCORDING TO C THE FOLLOWING CODE ... C C IFL = 1 .... SUCCESSFUL TERMINATION. C C IFL = 2 .... UNSUCCESSFUL TERMINATION. C CONSTRAINTS CANNOT BE SATISFIED. C PROBLEM IS INFEASIBLE. C C IFL = 3 .... LIMIT IMPOSED BY MXS REACHED C WITHOUT FINDING A SOLUTION. C C IFL = 4 .... PROGRAM ABORTED. C NUMERICAL DIFFICULTIES C DUE TO ILL-CONDITIONING. C C IFL = 5 .... NEQNS, NVARS, NEQC AND/OR C NIQC HAVE IMPROPER VALUES C OR E CONTAINS A ZERO COLUMN. C C IN ALL CASES THE OUTPUT PARAMETERS X,EL1N AND RES C WILL CONTAIN THE VALUES WHICH THEY REACHED AT TERMINATION. C C INTERMEDIATE PRINTING WILL BE TURNED OFF IF PSW = .FALSE. . C ON THE OTHER HAND, DETAILS OF EACH MINIMIZATION CYCLE C WILL BE PRINTED IF PSW IS SET TO .TRUE. C ---------------------------------------------------------- C C +++++ REMARKS AND USER CAUTIONS +++++ C ---------------------------------------------------------- C 1. BEYOND SOME PRECAUTIONARY STEPS TAKEN IN C CERTAIN DIVISIONS, NO SPECIAL C OVERFLOW/UNDERFLOW PROTECTION IS PROVIDED. C 2. ALL TOLERANCES FOR CHECKING ZEROS AND LINEAR C DEPENDENCIES ARE DETERMINED FROM THE QUANTITY C EPS WHICH APPEARS IN DATA DECLARATIONS IN CL1 C AND SEVERAL OF ITS SUBROUTINES. EPS CAN BE SET TO THE C LEAST POSITIVE NUMBER SATISFYING (1.0 + EPS) .GT. 1.0 C IN THE PRECISION OF ARITHMETIC BEING USED. WITH THIS C SETTING, CL1 USES AN EXTREMELY STRICT ZERO TOLERANCE. C FOR A MORE FORGIVING VERSION OF CL1, EPS MAY BE REMOVED C FROM THE DATA DEFINITIONS AND INCLUDED AS A USER-SPECIFIED C ZERO-TESTING PARAMETER IN THE ARGUMENT LIST. IN SUCH AN C EVENT, IF THE PROBLEM DATA IS GIVEN TO NDIG SIGNIFICANT C DIGITS, THEN 10.0**(-NDIG) IS A REASONABLE CHOICE C FOR THE VALUE OF EPS. C OVERFLOW CHECKING PRIOR TO DIVISION IS C DONE USING THE QUANTITY BIG, ALSO SPECIFIED AS DATA. C BIG SHOULD BE THE LARGEST REPRESENTABLE FLOATING C POINT NUMBER. C 3. THIS IS A SINGLE PRECISION VERSION OF CL1. C TO CHANGE THIS CODE INTO DOUBLE PRECISION ... C A. CHANGE ALL OCCURRENCES OF - REAL - C DECLARATIONS TO - DOUBLE PRECISION - C B. CHANGE ALL OCCURRENCES OF - SYSTEM ROUTINES - C (AS LISTED IN THE HEADING OF EACH SUBROUTINE) C TO THEIR CORRESPONDING DOUBLE PRECISION VERSIONS C C. CHANGE ALL OCCURRENCES OF THE STRINGS E+ AND E- C TO D+ AND D- RESPECTIVELY C D. CHANGE ALL BASIC LINEAR ALGEBRA ROUTINES (BLAS) C TO THEIR DOUBLE-PRECISION EQUIVALENTS C E. BOTH EPS AND BIG WILL HAVE TO BE CHANGED C TO THEIR DOUBLE PRECISION EQUIVALENTS C F. THE REFERENCES TO - IFIX - IN SUBROUTINES C - RESID - AND - GETV - MUST BE CHANGED C FROM THE FORM C IFIX(FLOAT(K)*UNIF(...)) C TO THE FORM C IFIX(FLOAT(K)*SNGL(UNIF(...))) C G. THE REFERENCE TO - FLOAT - IN SUBROUTINE C - RESID - MUST BE CHANGED FROM THE FORM C FLOAT(...) C TO THE FORM C DBLE(FLOAT(...)) C H. REMOVE THESE COMMENT CARDS (3., 3.A.-3.H.). C ---------------------------------------------------------- C C +++++ PORTABILITY NOTE +++++ C ---------------------------------------------------------- C THE INTRINSIC FORTRAN FUNCTIONS WHICH HAVE BEEN USED C ARE NOT DECLARED IN TYPE STATEMENTS IN CL1 OR ANY C OF ITS SUBROUTINES. A VERY FEW FORTRAN COMPILERS DO C REQUIRE SUCH TYPING, AND THEY WILL FAIL QUITE EXPLICITLY C BECAUSE OF OUR OMISSION. ON THE OTHER HAND, SOME FEW C OTHER FORTRAN COMPILERS DO NOT EXPECT SUCH TYPING, AND C THEIR MODE OF FAILURE IS MORE SUBTILE. THESE COMPILERS C MAY ATTACH THE NAMES OF THE INTRINSIC FUNCTIONS TO THE C LIST OF SUBROUTINES TO BE OBTAINED FROM EXTERNAL SOURCES, C CAUSING A FAULT IN LIBRARY SEARCHING DURING PRE-EXECUTION C LOADING. SINCE THERE IS NO WAY TO WIN THE GAME, WE C HAVE OPTED TO LOSE IT, IN THOSE FEW CASES WHERE WE MUST, C IN THE MOST IMMEDIATE AND OBVIOUS FASHION. C NOTE -- INTRINSIC FUNCTIONS ARE LISTED AMONG THE C - SYSTEM ROUTINES - C AT THE START OF EACH SUBPROGRAM. C C THE QUANTITIES EPS AND BIG WHICH ARE TO BE FOUND C IN DATA STATEMENTS ARE MACHINE DEPENDENT. SEE 2. UNDER C + REMARKS AND USER CAUTIONS +. C C MANY OF THE COMPUTATIONS IN THIS PROGRAM ARE C DEFINED IN TERMS OF THE BASIC LINEAR ALGEBRA C SUBROUTINES (BLAS). REFER TO THE APPENDIX C OF THIS CODE FOR FURTHER INFORMATION. C ********************************************************** C INTEGER DDX,GRDX,ICYC,IADDC,IDELC INTEGER PX,PTEX,RRX,TOPX,ZZX REAL ALPHA,AMAG,CGMAG,PEN,PENPAR C C //////////////// BEGIN PROGRAM ///////////////////////// C CALL SETUP * (NEQNS,NEQC,NIQC,NVARS,DDX,GRDX,PX,PTEX,RRX, * TOPX,ZZX,ICYC,IFL,E,NER,AMAG,CGMAG,PENPAR) C C 10 CONTINUE C CALL NEWPEN * (IADDC,IDELC,NACT,NEQNS,NEQC, * NIQC,NVARS,IFL,E,NER,X,F,RES, * W(PTEX),ALPHA,PENPAR,INDX) C 20 CONTINUE C CALL UPDATE * (IADDC,IDELC,NACT,NEQNS,NEQC,NIQC,NVARS, * ICYC,IFL,MXS,E,NER,X,F,RES, * W(GRDX),EL1N,PEN,PENPAR, * INDX,W(ZZX),NVARS,W(DDX),W(RRX), * W(TOPX)) CALL MONIT * (NACT,NEQC,NIQC,NVARS,ICYC,PSW,X, * ALPHA,EL1N,PEN,PENPAR,INDX) CALL FINDP * (IDELC,NACT,NEQNS,NEQC,NIQC,NVARS,IFL,E,NER, * X,F,RES,W(GRDX),W(PX),EL1N,AMAG,CGMAG, * PENPAR,INDX,W(ZZX),NVARS,W(DDX), * W(RRX),W(TOPX)) CALL STEP * (IADDC,NACT,NEQNS,NEQC,NIQC,NVARS,IFL,E,NER, * X,RES,W(GRDX),W(PX),W(PTEX),ALPHA,PENPAR, * INDX,W(TOPX)) C IF (IFL .EQ. 0) GO TO 20 C IF (IFL .EQ. 2 * .AND. (CGMAG + PENPAR*AMAG) .NE. CGMAG) GO TO 10 C RETURN END C C --------------- C SECOND LEVEL SUBROUTINES -- C SETUP,NEWPEN,UPDATE,MONIT,FINDP,STEP,REFINE C --------------- C SUBROUTINE SETUP * (NEQNS,NEQC,NIQC,NVARS,DDX,GRDX,PX,PTEX,RRX, * TOPX,ZZX,ICYC,IFL,E,NER,AMAG,CGMAG,PENPAR) C INTEGER DDX,GRDX,ICYC,IFL,NEQC,NEQNS,NER INTEGER NIQC,NVARS,PX,PTEX,RRX,TOPX,ZZX REAL AMAG,CGMAG,E(NER,1),PENPAR C C *************** C CL1 VERSION. C C SET UP THE PROGRAM C PARAMETERS AND INDICES. C *************** C C +++++++++++++++ C SYSTEM ROUTINES ABS C +++++++++++++++ C INTEGER I,J,NCOLS,NQNP1 REAL OCT,TMP,ZERO C DATA OCT /8.0E+00/ DATA ZERO /0.0E+00/ C C ///////////////// BEGIN PROGRAM ////////////////// C C *************** C CHECK VALIDITY OF PROBLEM DIMENSIONS C *************** C NCOLS = NEQNS + NEQC + NIQC IF (NVARS .GE. 1 * .AND. NEQC .GE. 0 * .AND. NIQC .GE. 0 * .AND. NEQNS .GE. 0 * .AND. NCOLS .GE. NVARS * .AND. NER .GE. NVARS) GO TO 10 IFL = 5 GO TO 100 10 CONTINUE C C *************** C SET UP INDICES FOR THE TEMPORARY STORAGE VECTOR W. C *************** C NQNP1 = NEQNS + 1 GRDX = 1 PX = GRDX + NVARS PTEX = PX + NVARS DDX = PTEX + NCOLS RRX = DDX + NVARS ZZX = RRX + (((NVARS + 1)*(NVARS + 2))/2) TOPX = ZZX + NVARS*NVARS C C *************** C AMAG IS A ROUGH ESTIMATE OF THE NORM OF A. C CGMAG IS A ROUGH ESTIMATE OF THE NORM OF (G,C). C TOGETHER THEY ARE USED TO DETERMINE WHEN THE C PENALTY PARAMETER IS TOO SMALL AND WHEN THE C RESTRICTED GRADIENT IS ZERO. C *************** C AMAG = ZERO CGMAG = ZERO IF (1 .GT. NEQNS) GO TO 50 DO 40 J=1,NEQNS TMP = ZERO DO 20 I=1,NVARS TMP = TMP + ABS(E(I,J)) 20 CONTINUE IF (TMP .GT. ZERO) GO TO 30 IFL = 5 GO TO 100 30 CONTINUE IF (TMP .GT. AMAG) AMAG = TMP 40 CONTINUE 50 CONTINUE IF (NQNP1 .GT. NCOLS) GO TO 90 DO 80 J=NQNP1,NCOLS TMP = ZERO DO 60 I=1,NVARS TMP = TMP + ABS(E(I,J)) 60 CONTINUE IF (TMP .GT. ZERO) GO TO 70 IFL = 5 GO TO 100 70 CONTINUE IF (TMP .GT. CGMAG) CGMAG = TMP 80 CONTINUE 90 CONTINUE C C *************** C INITIALIZE IFL,ICYC,PENPAR C *************** C IFL = 2 ICYC = -1 PENPAR = OCT 100 CONTINUE RETURN END SUBROUTINE NEWPEN * (IADDC,IDELC,NACT,NEQNS,NEQC,NIQC, * NVARS,IFL,E,NER,X,F,RES,PTE, * ALPHA,PENPAR,INDX) C INTEGER IADDC,IDELC,IFL,INDX(1),NACT INTEGER NEQC,NEQNS,NER,NIQC,NVARS REAL ALPHA,E(NER,1),F(1),PENPAR,PTE(1),RES(1),X(1) C C *************** C CL1 VERSION. C C BEGIN A ROUND OF MINIMIZATION STEPS C WITH A NEW PENALTY PARAMETER VALUE. C *************** C C +++++++++++++++ C BLAS SDOT C +++++++++++++++ C INTEGER I,NCOLS REAL OCT,ONE,ZERO C REAL SDOT C DATA ZERO /0.0E+00/ DATA ONE /1.0E+00/ DATA OCT /8.0E+00/ C C ///////////////// BEGIN PROGRAM ////////////////// C C *************** C SET PENALTY PARAMETER VALUE. C ERASE RECORD OF ACTIVE EQUATION/CONSTRAINTS. C *************** C IF (IFL .NE. 2) GO TO 20 NCOLS = NEQNS + NEQC + NIQC IFL = 0 NACT = 0 IADDC = 0 IDELC = 0 ALPHA = ZERO PENPAR = PENPAR/OCT C C *************** C INITIALIZE INDX,RES,PTE,INDX C *************** C DO 10 I=1,NCOLS RES(I) = SDOT(NVARS,E(1,I),1,X,1) - F(I) PTE(I) = ZERO INDX(I) = I 10 CONTINUE 20 CONTINUE RETURN END SUBROUTINE UPDATE * (IADDC,IDELC,NACT,NEQNS,NEQC,NIQC,NVARS, * ICYC,IFL,MXS,E,NER,X,F,RES,GRD, * EL1N,PEN,PENPAR,INDX,ZZ, * NZZR,DD,RR,W) C INTEGER IADDC,IDELC,ICYC,IFL,INDX(1),MXS INTEGER NACT,NEQC,NEQNS,NER,NIQC,NVARS,NZZR REAL DD(1),E(NER,1),EL1N,F(1),GRD(1),PEN,PENPAR REAL RES(1),RR(1),W(1),X(1),ZZ(NZZR,1) C C *************** C CL1 VERSION. C C PREPARATION FOR NEXT MINIMIZATION STEP. C *************** C C +++++++++++++++ C SYSTEM ROUTINES ABS C +++++++++++++++ C INTEGER NALLQ,NCOLS REAL ONE,ZERO C DATA ONE /1.0E+00/ DATA ZERO /0.0E+00/ C C ///////////////// BEGIN PROGRAM ////////////////// C C *************** C DETERMINE THE ACTIVE EQUATIONS AND ACTIVE C CONSTRAINTS. COMPUTE RESIDUALS AND FUNCTION VALUE. C UPDATE THE Z*D*R DECOMPOSITION. C *************** C NALLQ = NEQNS + NEQC NCOLS = NALLQ + NIQC IF (IFL .NE. 0) GO TO 20 ICYC = ICYC + 1 IF (ICYC .LE. MXS) GO TO 10 IFL = 3 GO TO 20 10 CONTINUE CALL DELCOL(IADDC,IDELC,NACT,NVARS,ZZ,NZZR,DD,RR,INDX) CALL RESID(IADDC,NACT,NCOLS,NVARS,E,NER,X,F,RES,INDX) CALL ADDCOL(IADDC,IDELC,NACT,NVARS,ZZ,NZZR,DD,RR, * E,NER,INDX,W) CALL OBJECT(IADDC,NACT,NEQNS,NALLQ,NCOLS,NVARS,E,NER, * RES,GRD,EL1N,PEN,PENPAR,INDX) 20 CONTINUE RETURN END SUBROUTINE MONIT * (NACT,NEQC,NIQC,NVARS,ICYC,PSW,X, * ALPHA,EL1N,PEN,PENPAR,INDX) C INTEGER ICYC,INDX(1),NACT,NEQC,NIQC,NVARS LOGICAL PSW REAL ALPHA,EL1N,PEN,PENPAR,X(1) C C *************** C CL1 VERSION. C C PRINT OUT INFORMATION AT C EACH STEP IF PSW = .TRUE. C C TO CHANGE THE OUTPUT MEDIUM, C CHANGE THE DATA DECLARATION C APPEARING BELOW. C *************** C INTEGER I,NXOUTP C DATA NXOUTP /6/ C C ///////////////// BEGIN PROGRAM ////////////////// C IF (.NOT. PSW) GO TO 20 WRITE (NXOUTP,6000) ICYC,ALPHA WRITE (NXOUTP,6010) (X(I),I=1,NVARS) WRITE (NXOUTP,6020) NACT IF (NACT .LE. 0) GO TO 10 WRITE (NXOUTP,6030) WRITE (NXOUTP,6040) (INDX(I),I=1,NACT) 10 CONTINUE WRITE (NXOUTP,6050) EL1N IF (NEQC + NIQC .GT. 0) WRITE (NXOUTP,6060) PENPAR,PEN 20 CONTINUE RETURN C 6000 FORMAT(20H0***** CYCLE NUMBER ,I5/ * 6X,13HSTEP TAKEN = ,1PE15.7/ * 6X,11HX-VECTOR...) 6010 FORMAT(5X,3(1PE15.7)) 6020 FORMAT(6X,38HNUMBER OF ACTIVE EQUATIONS/CONSTRAINTS,I5) 6030 FORMAT(6X,31HACTIVE EQUATIONS/CONSTRAINTS...) 6040 FORMAT(5X,7I5) 6050 FORMAT(6X,19HRESIDUAL L1 NORM = ,1PE15.7) 6060 FORMAT(6X,20HPENALTY PARAMETER = ,1PE10.2/ * 6X,25HPENALTY FUNCTION VALUE = ,1PE15.7) END SUBROUTINE FINDP * (IDELC,NACT,NEQNS,NEQC,NIQC,NVARS,IFL, * E,NER,X,F,RES,GRD,P,EL1N,AMAG,CGMAG, * PENPAR,INDX,ZZ,NZZR,DD,RR,W) C INTEGER IDELC,IFL,INDX(1),NACT,NEQC,NEQNS,NER,NIQC,NVARS,NZZR REAL AMAG,CGMAG,DD(1),E(NER,1),EL1N,F(1),GRD(1),P(1),PENPAR REAL RES(1),RR(1),W(1),X(1),ZZ(NZZR,1) C C *************** C CL1 VERSION. C C DETERMINE DESCENT DIRECTION P C (OR DISCOVER OPTIMALITY) C *************** C C +++++++++++++++ C SYSTEM ROUTINES ABS C C BLAS SASUM,SCOPY,SSCAL C C EPS IS THE SMALLEST POSITIVE NUMBER WHICH C SATISFIES (1.0 + EPS) .GT. 1.0 IN THE C PRECISION OF THE ARITHMETIC BEING USED. C (ALTERNATIVELY, FOR LESS STRICT ZERO CHECKING, C EPS CAN BE SET TO A USER-SPECIFIED TOLERANCE.) C +++++++++++++++ C INTEGER COEFX,I,IX,NALLQ,NALQP1,NCOLS,NQNP1,TOPX LOGICAL FAIL REAL EPS,GRDNRM,ONE,PNRM,PROD,TEST,ZERO C REAL SASUM C C +++++++++++++++ C THE FOLLOWING DECLARATIONS ARE NECESSARY C FOR PORTABILITY WHEN SCOPY IS USED, AS C IT IS BELOW, TO FILL ARRAYS WITH A SINGLE VALUE C (ONE=UNITY AND ZERO=ZIP IN THIS CASE). C +++++++++++++++ C REAL UNITY(1),ZIP(1) EQUIVALENCE (ONE,UNITY(1)),(ZERO,ZIP(1)) C DATA EPS /1.19209290E-07/ DATA ONE /1.0E+00/ DATA ZERO /0.0E+00/ C C ///////////////// BEGIN PROGRAM ////////////////// C IDELC = 0 IF (IFL .NE. 0) GO TO 90 NALLQ = NEQNS + NEQC NALQP1 = NALLQ + 1 NCOLS = NALLQ + NIQC NQNP1 = NEQNS + 1 COEFX = 1 TOPX = COEFX + NVARS C C *************** C PROJECT THE NEGATIVE OF THE RESTRICTED GRADIENT C ONTO THE ORTHOGONAL COMPLEMENT OF THE SPACE C SPANNED BY THE ACTIVE COLUMNS. C *************** C CALL ZDRPOC(NVARS,NACT,ZZ,NZZR,DD,GRD,P,FAIL) IF (.NOT.FAIL) GO TO 5 IFL = 4 GO TO 90 5 CONTINUE CALL SSCAL(NVARS,-ONE,P,1) PNRM = SASUM(NVARS,P,1) GRDNRM = SASUM(NVARS,GRD,1) C C *************** C IF THE PROJECTION IS NOT ZERO, C IT WILL SERVE AS A DESCENT DIRECTION. C C OTHERWISE FIND THE REPRESENTATION OF C THE RESTRICTED GRADIENT AS A LINEAR C COMBINATION OF THE ACTIVE COLUMNS. C THE COEFFICIENTS OF THE LINEAR COMBINATION C ARE TO BE STORED IN THE ARRAY COEF C (THAT IS, IN W(COEFX),...,W(COEFX+NACT-1)). C *************** C IF (PNRM .GT. EPS*(AMAG*PENPAR + CGMAG)) GO TO 90 IF (NACT .EQ. 0) GO TO 50 CALL ZDRGNV(NVARS,NACT,ZZ,NZZR,RR,GRD,W(COEFX),FAIL) IF (.NOT. FAIL) GO TO 10 IFL = 4 GO TO 90 10 CONTINUE C C *************** C CONVERT THE COEFFICIENTS OF THE LINEAR C COMBINATION INTO A DESCENT DIRECTION P , C OR DETERMINE OPTIMALITY. C C IF THE OPTIMALITY TEST IS NOT SATISFIED, C GETV WILL INDICATE AN EQUATION/CONSTRAINT C TO BE DELETED FROM ACTIVITY BY THE VALUE C OF IDELC. FOR OPTIMALITY, IDELC=0. C *************** C CALL GETV(IDELC,NACT,NVARS,NEQNS,NALLQ,E,NER,GRD, * W(COEFX),PENPAR,INDX) PNRM = ZERO IF (IDELC .EQ. 0) GO TO 20 CALL ZDRGIT(NVARS,NACT,ZZ,NZZR,RR,W(COEFX),P,FAIL,W(TOPX)) IF (.NOT. FAIL) PNRM = SASUM(NVARS,P,1) IF (.NOT. FAIL) GO TO 20 IFL = 4 GO TO 90 20 CONTINUE C C *************** C IF A DESCENT DIRECTION P COULD HAVE BEEN FOUND, C IT HAS BEEN OBTAINED BY THIS POINT IN THE PROGRAM. C C CHECK FOR OPTIMALITY. C C PNRM HAS BEEN SET EXACTLY ZERO C AFTER THE CALL TO SUBROUTINE GETV C IF THE OPTIMALITY CONDITIONS ARE SATISFIED. C THE CHECK BELOW HAS BEEN MADE SOMEWHAT C COMPLICATED TO ALLOW FOR THE RARE EVENT THAT C THE RESTRICTED GRADIENT IS ZERO AND NO C COLUMNS ARE ACTIVE, OR THAT THE L1 NORM OF C (A-TRANSPOSE) * X - F C IS COMPUTATIONALLY ZERO. C (THE CALL TO THE SUBROUTINE REFINE C MAY BE OMITTED, IF DESIRED.) C *************** C IF (PNRM .LE. EPS*(AMAG*PENPAR + CGMAG)) GO TO 50 DO 40 I=1,NEQNS TEST = ABS(F(I)) DO 30 IX=1,NVARS PROD = ABS(E(IX,I)*X(IX)) IF (PROD .GT. TEST) TEST = PROD 30 CONTINUE IF (ABS(RES(I)) .GT. EPS*TEST) GO TO 90 40 CONTINUE 50 CONTINUE IFL = 1 CALL REFINE(NACT,NEQNS,NCOLS,NVARS,IFL,E,NER, * X,F,EL1N,RES,INDX,ZZ,NZZR,RR,W) IF (IFL .NE. 1) GO TO 90 C C *************** C IF THE PROBLEM HAS CONSTRAINTS, C CHECK FEASIBILITY. C *************** C IF (NQNP1 .GT. NCOLS) GO TO 90 DO 80 I=NQNP1,NCOLS TEST = ABS(F(I)) DO 60 IX=1,NVARS PROD = ABS(E(IX,I)*X(IX)) IF (PROD .GT. TEST) TEST = PROD 60 CONTINUE TEST = EPS*TEST IF (I .GT. NALLQ) GO TO 70 IF (ABS(RES(I)) .LE. TEST) GO TO 80 IFL = 2 GO TO 90 70 CONTINUE IF (RES(I) .GE. (-TEST)) GO TO 80 IFL = 2 GO TO 90 80 CONTINUE 90 CONTINUE RETURN END SUBROUTINE STEP * (IADDC,NACT,NEQNS,NEQC,NIQC,NVARS,IFL, * E,NER,X,RES,GRD,P,PTE,ALPHA,PENPAR,INDX,ALF) C INTEGER IADDC,IFL,INDX(1),NACT,NEQC,NEQNS,NER,NIQC,NVARS REAL ALPHA,ALF(1),E(NER,1),GRD(1),P(1) REAL PENPAR,PTE(1),RES(1),X(1) C C *************** C CL1 VERSION. C C PIECEWISE LINEAR LINE SEARCH. C *************** C C +++++++++++++++ C SYSTEM ROUTINES ABS,SIGN C C BLAS SASUM,SAXPY,SDOT C C EPS IS THE SMALLEST POSITIVE NUMBER WHICH C SATISFIES (1.0 + EPS) .GT. 1.0 IN THE C PRECISION OF THE ARITHMETIC BEING USED. C (ALTERNATIVELY, FOR LESS STRICT ZERO CHECKING, C EPS CAN BE SET TO A USER-SPECIFIED TOLERANCE.) C C BIG IS THE LARGEST POSITIVE NUMBER C WHICH CAN BE REPRESENTED IN THE C PRECISION OF THE ARITHMETIC BEING USED. C +++++++++++++++ C INTEGER I,IIN,IX,JX,NACTP1,NALLQ,NCOLS,NUM,NUMNAC REAL BIG,DEN,EPS,GRDNRM,ONE,PNRM,PTG REAL RATIO,RESID,TMP,TWO,ZERO C REAL SASUM,SDOT C DATA BIG /1.0E+30/ DATA EPS /1.19209290E-07/ DATA ONE /1.0E+00/ DATA TWO /2.0E+00/ DATA ZERO /0.0E+00/ C C ///////////////// BEGIN PROGRAM ////////////////// C C *************** C THIS ROUTINE DETERMINES ALL OF THE RATIOS ALF C OF THE FORM C -RES(I)/((E(.,I)-TRANSP)*P), C FOR I = K+1,...,MPL C WHICH ARE NONNEGATIVE AND HENCE INDICATE DISTANCES C FROM THE POINT X TO BREAKPOINTS WHICH WILL C BE ENCOUNTERED IN TRAVEL ALONG DIRECTION P. C THE INDEX VECTOR INDX IS REARRANGED SO THAT C ITS K+1 THROUGH NUM COMPONENTS CORRESPOND TO C THESE NONNEGATIVE RATIOS. C THE RESULTS ARE HEAPED SO THAT THE ALF VALUES CAN C BE INSPECTED IN ORDER FROM SMALLEST TO LARGEST. C THE BREAKPOINT ALPHA GIVING THE MINIMUM OBJECTIVE C FUNCTION VALUE IS FOUND, AND X IS C ADJUSTED TO X + ALPHA*P . C C THE INNER PRODUCTS (E(.,I)-TRANSPOSE)*P ARE SAVED C FOR LATER USE IN UPDATING THE RESIDUAL VALUES. C *************** C ALPHA = ZERO IF (IFL .NE. 0) GO TO 90 NALLQ = NEQNS + NEQC NCOLS = NALLQ + NIQC NACTP1 = NACT + 1 NUM = 0 IF (1 .GT. NACT) GO TO 20 DO 10 I=1,NACT IX = INDX(I) PTE(IX) = SDOT(NVARS,E(1,IX),1,P,1) 10 CONTINUE 20 CONTINUE IF (NACTP1 .LE. NCOLS) GO TO 30 IFL = 1 GO TO 90 30 CONTINUE DO 50 I=NACTP1,NCOLS IX = INDX(I) RESID = RES(IX) DEN = SDOT(NVARS,E(1,IX),1,P,1) PTE(IX) = DEN IF (SIGN(ONE,RESID) .EQ. SIGN(ONE,DEN) * .AND. RESID .NE. ZERO) GO TO 50 RESID = ABS(RESID) DEN = ABS(DEN) IF (DEN .GE. ONE) GO TO 40 IF (RESID .GE. DEN*BIG) GO TO 50 40 CONTINUE RATIO = RESID/DEN NUM = NUM + 1 NUMNAC = NUM + NACT JX = INDX(NUMNAC) INDX(NUMNAC) = IX INDX(I) = JX ALF(NUM) = RATIO 50 CONTINUE IF (NUM .GT. 0) GO TO 60 IFL = 2 GO TO 90 60 CONTINUE C C *************** C HEAP THE POSITIVE RATIOS C *************** C CALL DKHEAP(.TRUE.,NUM,INDX(NACTP1),ALF) C C *************** C TRAVEL ALONG P UNTIL NO FURTHER DECREASE IN THE C PENALTY FUNCTION IS POSSIBLE C *************** C IIN = NUM PTG = SDOT(NVARS,GRD,1,P,1) PNRM = SASUM(NVARS,P,1) GRDNRM = SASUM(NVARS,GRD,1) DO 70 I=1,NUM IX = INDX(NACTP1) TMP = -SIGN(ONE,RES(IX)) IF (IX .LE. NALLQ) TMP = TMP*TWO IF (IX .LE. NEQNS) TMP = TMP*PENPAR PTG = PTG + TMP*PTE(IX) IF (PTG .GE. (-EPS*GRDNRM*PNRM)) GO TO 80 CALL DKHEAP(.FALSE.,IIN,INDX(NACTP1),ALF) 70 CONTINUE IFL = 2 GO TO 90 80 CONTINUE IADDC = NACTP1 C C *************** C ADJUST X TO X + ALPHA*P C *************** C ALPHA = ALF(1) CALL SAXPY(NVARS,ALPHA,P,1,X,1) 90 CONTINUE RETURN END SUBROUTINE REFINE(NACT,NEQNS,NCOLS,NVARS,IFL,E,NER, * X,F,EL1N,RES,INDX,ZZ,NZZR,RR,W) C INTEGER IFL,INDX(1),NACT,NCOLS,NEQNS,NER,NVARS,NZZR REAL E(NER,1),EL1N,F(1),RES(1),RR(1),W(1),X(1),ZZ(NZZR,1) C C *************** C A ROUTINE FOR REFINING THE SOLUTION C PRODUCED BY CL1. C C (THIS ROUTINE MAY BE OMITTED IF DESIRED.) C *************** C C +++++++++++++++ C SYSTEM ROUTINES ABS C C BLAS SDOT C +++++++++++++++ C INTEGER I,IX LOGICAL FAIL REAL TMP,ZERO C REAL SDOT C DATA ZERO /0.0E+00/ C C /////////////// BEGIN PROGRAM /////////////// C IF (NACT .EQ. 0) GO TO 40 DO 10 I=1,NACT IX = INDX(I) RES(I) = F(IX) 10 CONTINUE CALL ZDRGIT(NVARS,NACT,ZZ,NZZR,RR,RES,X,FAIL,W) IF (.NOT. FAIL) GO TO 20 IFL = 4 GO TO 40 20 CONTINUE EL1N = ZERO DO 30 I=1,NCOLS TMP = SDOT(NVARS,E(1,I),1,X,1) - F(I) RES(I) = TMP IF (I .LE. NEQNS) EL1N = EL1N + ABS(TMP) 30 CONTINUE 40 CONTINUE RETURN END C C --------------- C THIRD LEVEL SUBROUTINES -- C DELCOL,RESID,ADDCOL,OBJECT,GETV C --------------- C SUBROUTINE DELCOL(IADDC,IDELC,NACT,NROW,ZZ,NZZR,DD,RR,INDX) C INTEGER INDX(1),NZZR,NACT,IDELC,IADDC,NROW REAL DD(1),RR(1),ZZ(NZZR,1) C C *************** C CL1 VERSION OF DELCOL. C C THIS ROUTINE ADMINISTERS THE DELETION OF THE COLUMN C INDICATED BY THE VALUE OF IDELC C FROM AN NROW BY NACT Z*D*R DECOMPOSITION. C NOTE THAT THE VALUE OF IDELC C IS THE NUMBER OF A COLUMN IN THE DECOMPOSITION C RATHER THAN A NUMBER WHICH REFERS TO C A COLUMN IN THE MATRIX E. C (THE E-COLUMN NUMBERS CORRESPONDING TO C THE COLUMNS OF THE FACTORIZATION ARE TO BE C FOUND IN INDX(1),...,INDX(NACT) . C THE CONTENTS OF INDX(NACT+1),...,INDX(IADDC) C INDICATE COLUMNS OF E WHICH ARE SLATED FOR C ADDITION TO THE DECOMPOSITION.) C THE VECTOR INDX IS REARRANGED BY C PERMUTING THE ELEMENT WHICH CORRESPONDS TO C THE DELETION OUT TO THE IADDC-TH POSITION. C NACT AND IADDC ARE DECREASED ACCORDINGLY. C *************** C INTEGER I,IDLP1,IXDLC LOGICAL FAIL C C ///////////////// BEGIN PROGRAM ////////////////// C IF (IDELC .EQ. 0) GO TO 20 IDLP1 = IDELC + 1 IXDLC = INDX(IDELC) DO 10 I=IDLP1,IADDC INDX(I - 1) = INDX(I) 10 CONTINUE INDX(IADDC) = IXDLC IADDC = IADDC - 1 CALL ZDRCOU(NROW,NACT,ZZ,NZZR,DD,RR,IDELC,FAIL) IDELC = IXDLC 20 CONTINUE RETURN END SUBROUTINE RESID(IADDC,NACT,NCOLS,NVARS,E,NER,X,F,RES,INDX) C INTEGER IADDC,INDX(1),NACT,NCOLS,NER,NVARS REAL E(NER,1),F(1),RES(1),X(1) C C C *************** C COMPUTE THE RESIDUALS C (E(.,IX)-TRANSP)*X - F(IX) . C THE RESIDUALS ARE STORED IN THE ARRAY RES. C INDX IS REARRANGED SO THAT THE ZERO RESIDUALS C CORRESPOND TO INDX(1),...,INDX(IADDC) . C *************** C C +++++++++++++++ C SYSTEM ROUTINES ABS,IFIX,FLOAT,SQRT C C BLAS SDOT C C EPS IS THE SMALLEST POSITIVE NUMBER WHICH C SATISFIES (1.0 + EPS) .GT. 1.0 IN THE C PRECISION OF THE ARITHMETIC BEING USED. C (ALTERNATIVELY, FOR LESS STRICT ZERO CHECKING, C EPS CAN BE SET TO A USER-SPECIFIED TOLERANCE.) C +++++++++++++++ C INTEGER I,IADP1,IDUMMY,IRAND,IX,J,NACTP1 REAL EPS,PROD,TEMP,TEST,TOL,ZERO C REAL SDOT,UNIF01 C DATA EPS /1.19209290E-07/ DATA ZERO /0.0E+00/ C C ///////////////// BEGIN PROGRAM ////////////////// C TOL = EPS*SQRT(FLOAT(NVARS)) NACTP1 = NACT + 1 IF (1 .GT. IADDC) GO TO 20 C C *************** C ZERO OUT ALL RESIDUALS KNOWN TO BE ZERO. C *************** C DO 10 I=1,IADDC IX = INDX(I) RES(IX) = ZERO 10 CONTINUE 20 CONTINUE C C *************** C COMPUTE THE REMAINING RESIDUALS. C DETECT ANY MORE RESIDUALS WHICH C ARE COMPUTATIONALLY ZERO, AND C SET THEM EXACTLY ZERO. THEIR C ASSOCIATED INDICES ARE PERMUTED C SO THAT THEY ARE STORED IN C INDX(NACT+1),...,NACT(IADDC). C C (A FAIRLY TIGHT ZERO CHECK IS USED. C IT IS FAR LESS EXPENSIVE IN RUNNING C TIME TO NEGLECT AN EXTRA ZERO C RESIDUAL THAN TO ACCEPT IT AND RISK C INVOKING THE ANTI-CYCLING C MECHANISMS IN THE PROGRAM. C THE ACCURACY OF THE SOLUTION AS C FINALLY DETERMINED IS NOT AFFECTED.) C *************** C IADP1 = IADDC + 1 IF (IADP1 .GT. NCOLS) GO TO 50 DO 40 I=IADP1,NCOLS IX = INDX(I) TEMP = SDOT(NVARS,E(1,IX),1,X,1) - F(IX) TEST = ABS(F(IX)) DO 25 J=1,NVARS PROD = ABS(E(J,IX)*X(J)) IF (PROD .GT. TEST) TEST = PROD 25 CONTINUE TEST = TOL*TEST IF (ABS(TEMP) .LE. TEST) GO TO 30 RES(IX) = TEMP GO TO 40 30 CONTINUE IADDC = IADDC + 1 INDX(I) = INDX(IADDC) INDX(IADDC) = IX RES(IX) = ZERO 40 CONTINUE 50 CONTINUE C C *************** C IF ANY NEW ZERO RESIDUALS HAVE C BEEN FOUND, RANDOMIZE THEIR C ORDERING AS AN ANTI-CYCLING C DEVICE FOR ADDCOL. C *************** C IF (IADDC .LE. NACTP1) GO TO 70 DO 60 I=NACTP1,IADDC IRAND = I + IFIX(FLOAT(IADDC - I + 1)*UNIF01(0,IDUMMY)) IX = INDX(IRAND) INDX(IRAND) = INDX(I) INDX(I) = IX 60 CONTINUE 70 CONTINUE RETURN END SUBROUTINE ADDCOL(IADDC,IDELC,NACT,NVARS,ZZ,NZZR,DD,RR, * E,NER,INDX,W) C INTEGER IADDC,IDELC,INDX(1),NACT,NER,NVARS,NZZR REAL DD(1),E(NER,1),RR(1) REAL W(1),ZZ(NZZR,1) C C *************** C CL1 VERSION OF ADDCOL. C C THIS ROUTINE ADMINISTERS THE ADJUSTMENT OF THE C Z*D*R DECOMPOSITION FOR ANY NEW ZERO RESIDUALS. C THE DATA CORRESPONDING TO THE ZERO RESIDUALS IS INDEXED C IN INDX(NACT+1),...,INDX(IADDC). C *************** C C +++++++++++++++ C BLAS SASUM C C EPS IS THE SMALLEST POSITIVE NUMBER WHICH C SATISFIES (1.0 + EPS) .GT. 1.0 IN THE C PRECISION OF THE ARITHMETIC BEING USED. C (ALTERNATIVELY, FOR LESS STRICT ZERO CHECKING, C EPS CAN BE SET TO A USER-SPECIFIED TOLERANCE.) C +++++++++++++++ C INTEGER I,ISTRT,IX,NACTP1,TOPX LOGICAL FAIL REAL COLNRM,EPS,PRJNRM C REAL SASUM C DATA EPS /1.19209290E-07/ C C ///////////////// BEGIN PROGRAM ////////////////// C TOPX = NVARS + 1 ISTRT = NACT + 1 IF (ISTRT .GT. IADDC) GO TO 20 C C *************** C CANDIDATES FOR ADDITION TO THE Z*D*R C FACTORIZATION ARE INSPECTED IN RANDOM C ORDER TO HINDER CYCLING. C THE RANDOMIZATION WAS CARRIED OUT BY RESID. C C IF A CANDIDATE HAS JUST BEEN RELEASED C FROM THE FACTORIZATION OR IS DEPENDENT UPON THE C COLUMNS IN THE FACTORIZATION, C THEN IT IS OMITTED FROM ADDITION. C C UPON EXIT, INDICES OF SUCH OMITTED C COLUMNS ARE TO BE FOUND IN C INDX(NACT+1),...,INDX(IADDC) . C *************** C DO 10 I=ISTRT,IADDC NACTP1 = NACT + 1 IX = INDX(I) CALL ZDRPOC(NVARS,NACT,ZZ,NZZR,DD,E(1,IX),W,FAIL) COLNRM = SASUM(NVARS,E(1,IX),1) PRJNRM = SASUM(NVARS,W,1) IF (PRJNRM .LE. EPS*COLNRM .OR. IX .EQ. IDELC) GO TO 10 INDX(I) = INDX(NACTP1) INDX(NACTP1) = IX CALL ZDRCIN(NVARS,NACT,ZZ,NZZR,DD,RR,E(1,IX),FAIL,W) 10 CONTINUE 20 CONTINUE RETURN END SUBROUTINE OBJECT(IADDC,NACT,NEQNS,NALLQ,NCOLS,NVARS, * E,NER,RES,GRD,EL1N,PEN,PENPAR,INDX) C INTEGER IADDC,INDX(1),NACT,NALLQ,NEQNS,NCOLS,NER,NVARS REAL E(NER,1),EL1N,GRD(1),PEN,PENPAR,RES(1) C C *************** C CL1 VERSION OF OBJECT. C C THIS ROUTINE ADMINISTERS THE EVALUATION OF THE C PENALTY (OBJECTIVE) FUNCTION GIVEN THE EQUATION C AND CONSTRAINT RESIDUALS. IT ALSO COMPUTES THE C RESTRICTED GRADIENT OF THE FUNCTION. C C COLUMNS WHICH ARE NOT IN THE Z*D*R FACTORIZATION C BUT WHICH ARE ASSOCIATED WITH ZERO RESIDUALS MUST C BE INCLUDED IN THE RESTRICTED GRADIENT WITH RANDOM C SIGNS AS AN ANTI-CYCLING DEVICE. C THE INDICES OF THESE COLUMNS ARE TO BE C FOUND IN INDX(NACT+1),...,INDX(IADDC) C *************** C C +++++++++++++++ C SYSTEM ROUTINES ABS,SIGN C C BLAS SAXPY,SCOPY C +++++++++++++++ C INTEGER I,IDUMMY,IX,NACTP1 REAL HALF,ONE,SGN,TMP,ZERO C REAL UNIF01 C C +++++++++++++++ C THE FOLLOWING DECLARATIONS ARE NECESSARY C FOR PORTABILITY WHEN SCOPY IS USED, AS C IT IS BELOW, TO FILL ARRAYS WITH A SINGLE C VALUE (ZERO=ZIP IN THIS CASE). C +++++++++++++++ C REAL ZIP(1) EQUIVALENCE (ZERO,ZIP(1)) C DATA HALF /0.5E+00/ DATA ONE /1.0E+00/ DATA ZERO /0.0E+00/ C C ///////////////// BEGIN PROGRAM ////////////////// C NACTP1 = NACT + 1 EL1N = ZERO PEN = ZERO CALL SCOPY(NVARS,ZIP,0,GRD,1) IF (NACTP1 .GT. NCOLS) GO TO 30 DO 20 I=NACTP1,NCOLS IX = INDX(I) TMP = RES(IX) SGN = SIGN(ONE,TMP) IF (I .LE. IADDC) SGN = SIGN(ONE,HALF - UNIF01(0,IDUMMY)) IF (IX .GT. NALLQ .AND. SGN .GT. ZERO) GO TO 20 TMP = ABS(TMP) IF (IX .GT. NEQNS) GO TO 10 EL1N = EL1N + TMP TMP = TMP*PENPAR SGN = SGN*PENPAR 10 CONTINUE PEN = PEN + TMP CALL SAXPY(NVARS,SGN,E(1,IX),1,GRD,1) 20 CONTINUE 30 CONTINUE RETURN END SUBROUTINE GETV(IDELC,NACT,NVARS,NEQNS,NALLQ,E,NER, * GRD,COEF,PENPAR,INDX) C INTEGER IDELC,INDX(1),NACT,NALLQ,NEQNS,NER,NVARS REAL COEF(1),E(NER,1),GRD(1),PENPAR C C *************** C CL1 VERSION. C C SET UP THE RIGHT-HAND-SIDE VECTOR C (AND STORE IN THE ARRAY COEF) C FOR THE LINEAR PROBLEM WHICH DETERMINES C A DESCENT DIRECTION P IN THE CASE WHERE C THE PROJECTION OF THE RESTRICTED GRADIENT IS ZERO. C *************** C C +++++++++++++++ C SYSTEM ROUTINES ABS,FLOAT,IFIX,SIGN C C BLAS SAXPY C C EPS IS THE SMALLEST POSITIVE NUMBER WHICH C SATISFIES (1.0 + EPS) .GT. 1.0 IN THE C PRECISION OF THE ARITHMETIC BEING USED. C (ALTERNATIVELY, FOR LESS STRICT ZERO CHECKING, C EPS CAN BE SET TO A USER-SPECIFIED TOLERANCE.) C +++++++++++++++ C INTEGER I,IDUMMY,IRAND,IX REAL CF,EPS,ONE,OPE,S,TMP,TMPSAV,ZERO C REAL UNIF01 C DATA EPS /1.19209290E-07/ DATA ONE /1.0E+00/ DATA ZERO /0.0E+00/ C C ///////////////// BEGIN PROGRAM ////////////////// C C *************** C FIND THE MOST OUT-OF-KILTER C COEFFICIENT. BEGIN INSPECTING C THE COEFFICIENTS AT A RANDOM INDEX C TO HINDER CYCLING. SET COEF C TO ZERO ON THE FLY. C *************** C OPE = ONE + EPS IDELC = 0 TMPSAV = ZERO IF (1 .GT. NACT) GO TO 40 IRAND = IFIX(FLOAT(NACT)*UNIF01(0,IDUMMY)) DO 30 I=1,NACT IRAND = IRAND + 1 IF (IRAND .GT. NACT) IRAND = 1 IX = INDX(IRAND) CF = COEF(IRAND) COEF(IRAND) = ZERO IF (IX .GT. NALLQ) GO TO 10 IF (IX .LE. NEQNS) CF = CF/PENPAR TMP = OPE - ABS(CF) GO TO 20 10 CONTINUE TMP = CF + EPS 20 CONTINUE IF (TMP .GE. TMPSAV) GO TO 30 IDELC = IRAND S = SIGN(ONE,CF) TMPSAV = TMP 30 CONTINUE C C *************** C IF NO COEFFICIENTS ARE OUT OF KILTER, C THEN RETURN. OTHERWISE SET A C VALUE IN AN APPROPRIATE COMPONENT C (INDICATED BY IDELC) OF COEF C AND ADJUST THE RESTRICTED GRADIENT C IF NECESSARY. C *************** C IF (IDELC .EQ. 0) GO TO 40 COEF(IDELC) = ONE IX = INDX(IDELC) IF (IX .GT. NALLQ) GO TO 40 COEF(IDELC) = -S TMP = -S IF (IX .LE. NEQNS) TMP = TMP*PENPAR CALL SAXPY(NVARS,TMP,E(1,IX),1,GRD,1) 40 CONTINUE RETURN END C C --------------- C FOURTH LEVEL SUBROUTINES -- C DKHEAP,UNIF01,ZDRCIN,ZDRCOU, C ZDRGIT,ZDRGNV,ZDRPOC C --------------- C SUBROUTINE DKHEAP(MAKE,IR,INDX,ARAY) C INTEGER INDX(1),IR LOGICAL MAKE REAL ARAY(1) C C *************** C AN ADAPTATION OF D. E. KNUTH,S HEAPING C ROUTINES (SEE VOLUME 3 OF C THE ART OF COMPUTER PROGRAMMING ). C IF MAKE IS .TRUE., THE FULL HEAP BUILDING C PROCESS IS CARRIED OUT ON C ARAY(1),...,ARAY(IR) , C AND THE VALUE OF IR IS UNCHANGED. C IF MAKE IS .FALSE., ONE STEP OF THE SORTING C PROCESS IS CARRIED OUT TO PROVIDE THE NEXT C ELEMENT OF ARAY IN ORDER, AND THE VARIABLE C IR IS DECREASED BY ONE. THE INTERRUPTION OF THE C SORTING PHASE IS BUILT IN VIA THE FLAG ONCE. C INDX IS AN INDEX VECTOR ASSOCIATED WITH C ARAY WHICH MUST BE REARRANGED IN PARALLEL C WITH IT. C *************** C INTEGER I,IL,IT,J LOGICAL ONCE REAL T C C ///////////////// BEGIN PROGRAM ////////////////// C IF (IR .GT. 1) GO TO 5 IF (.NOT.MAKE) IR = 0 RETURN 5 CONTINUE C C *************** C TEST WHETHER OR NOT THE INITIAL C HEAP IS TO BE BUILT C *************** C IL = 1 IF (MAKE) IL = (IR/2) + 1 ONCE = .FALSE. C C *************** C THE LOOP BEGINS HERE C *************** C 10 CONTINUE IF (IL .GT. 1) GO TO 20 C C *************** C THE SORTING PHASE USES THIS BRANCH C *************** C IF (MAKE .OR. ONCE) RETURN ONCE = .TRUE. IT = INDX(IR) T = ARAY(IR) INDX(IR) = INDX(1) ARAY(IR) = ARAY(1) IR = IR - 1 IF (IR .GT. 1) GO TO 30 INDX(1) = IT ARAY(1) = T RETURN 20 CONTINUE C C *************** C THE HEAP-BUILDING PHASE USES THIS BRANCH C *************** C IL = IL - 1 IT = INDX(IL) T = ARAY(IL) 30 CONTINUE C C *************** C THE REMAINING STATEMENTS ARE COMMON C TO BOTH PHASES AND EMBODY THE C HEAP-RECTIFYING (SIFTING) SECTION C *************** C J = IL 40 CONTINUE I = J J = 2*J IF (J - IR) 50, 60, 70 50 CONTINUE IF (ARAY(J) .LE. ARAY(J + 1)) GO TO 60 J = J + 1 60 CONTINUE IF (T .LE. ARAY(J)) GO TO 70 INDX(I) = INDX(J) ARAY(I) = ARAY(J) GO TO 40 70 CONTINUE INDX(I) = IT ARAY(I) = T GO TO 10 END REAL FUNCTION UNIF01(ISEED,IX) C INTEGER ISEED,IX,IX0 C DATA IX0/2/ C C +++++++++++++++ C SYSTEM ROUTINES FLOAT,MOD C +++++++++++++++ C C -------------------------------------------------------------- C -------------------------------------------------------------- C C *****PURPOSE- C THIS FUNCTION RETURNS A PSEUDO-RANDOM NUMBER DISTRIBUTED C UNIFORMLY IN THE INTERVAL (0,1). C C *****PARAMETER DESCRIPTION- C ON INPUT- C C ISEED, IF IT IS NONZERO MODULO 9973, BECOMES THE C NEW SEED, I.E. IT REPLACES THE INTERNALLY STORED C VALUE OF IX0. ON MACHINES WHERE FORTRAN VARIABLES C RETAIN THEIR VALUES BETWEEN CALLS, THE INTERNALLY C STORED VALUE IF IX0 IS THE VALUE ASSIGNED TO IX IN C THE PREVIOUS INVOCATION OF UNIF01. OTHERWISE -- AND C IN THE FIRST CALL TO UNIF01 -- IX0=2. C C ON OUTPUT- C C IX IS THE NEXT INTEGER IN A PSEUDO-RANDOM SEQUENCE OF C INTEGERS BETWEEN 1 AND 9972 AND IS GENERATED FROM ITS C PREDECESSOR IX0 (I.E. FROM ISEED, IF ISEED IS NONZERO C MODULO 9973). IX IS THE VALUE WHICH ISEED SHOULD HAVE C IN THE NEXT INVOCATION OF UNIF01 TO GET THE NEXT C PSEUDO-RANDOM NUMBER. THE CALLER WILL OFTEN PASS THE C SAME VARIABLE FOR ISEED AS FOR IX, C E.G. X = UNIF01(IX,IX). C C *****APPLICATION AND USAGE RESTRICTIONS- C UNIF01 SHOULD ONLY BE USED WHEN PORTABILITY IS IMPORTANT AND A C COURSE RANDOM NUMBER GENERATOR SUFFICES. APPLICATIONS REQUIRING C A FINE, HIGH PRECISON GENERATOR SHOULD USE ONE WITH A MUCH C LARGER MODULUS. C C *****ALGORITHM NOTES- C UNIF01 WILL RUN ON ANY MACHINE HAVING AT LEAST 20 BITS OF AC- C CURACY FOR FIXED-POINT ARITHMITIC. IT IS BASED ON A GENERATOR C RECOMMENDED IN (3), WHICH PASSES THE SPECTRAL TEST WITH FLYING C COLORS -- SEE (1) AND (2). C C REFERENCES- C (1) HOAGLIN, D.C. (1976), THEORETICAL PROPERTIES OF CONGRUENTIAL C RANDOM-NUMBER GENERATORS- AN EMPIRICAL VIEW, C MEMORANDUM NS-340, DEPT. OF STATISTICS, HARVARD UNIV. C C (2) KNUTH, D.E. (1969), THE ART OF COMPUTER PROGRAMMING, VOL. 2 C (SEMINUMERICAL ALGORITHMS), ADDISON-WESLEY, READING, MASS. C C (3) SMITH, C.S. (1971), MULTIPLICATIVE PSEUDO-RANDOM NUMBER C GENERATORS WITH PRIME MODULUS, J. ASSOC. COMPUT. MACH. 18, C PP. 586-593. C C *****GENERAL- C C THIS SUBROUTINE WAS WRITTEN IN CONNECTION WITH RESEARCH C SUPPORTED BY THE NATIONAL SCIENCE FOUNDATION UNDER GRANTS C MCS-7600324, DCR75-10143, 76-14311DSS, AND MCS76-11989. C C PERMISSION FOR THE USE OF UNIF01 IN CL1 WAS C GENEROUSLY GIVEN BY V. KLEMA AND D. HOAGLIN. C C -------------------------------------------------------------- C -------------------------------------------------------------- C IF (ISEED .EQ. 0) GO TO 10 IX = MOD(ISEED, 9973) IF (IX .NE. 0) IX0 = IX C *** C IN ORDER THAT ALL FIXED-POINT CALCULATIONS REQUIRE ONLY 20 BIT C ARITHMETIC, WE USE TWO CALLS TO MOD TO COMPUTE C IX0 = MOD(3432*IX0, 9973). C *** 10 IX0 = MOD(52*MOD(66*IX0, 9973), 9973) IX = IX0 UNIF01 = FLOAT(IX0)/9973.0E+00 RETURN END SUBROUTINE ZDRCIN(N,K,ZZ,NZZR,DD,RR,COL,FAIL,W) C INTEGER K,N,NZZR LOGICAL FAIL REAL COL(1),DD(1),RR(1),W(1),ZZ(NZZR,1) C C *************** C PREPARED BY RICHARD BARTELS C THE UNIVERSITY OF WATERLOO C COMPUTER SCIENCE DEPARTMENT C LATEST UPDATE .... 30 NOVEMBER, 1979. C C GIVEN THE FACTORIZATION C C ZZ*DD*RR C C OF SOME N BY K MATRIX C C (0 .LE. K .LT. N) C (N .GE. 1), C C WHERE C C (ZZ-TRANSP)*(ZZ) = (DD-INV), C DD IS DIAGONAL AND NONSINGULAR, C AND C RR HAS ZEROS BELOW THE DIAGONAL, C C AND GIVEN A (K+1)TH COLUMN C TO BE ADDEDTO THE ORIGINAL MATRIX, C THIS PROGRAM UPDATES ZZ,DD AND RR. C C THE VALUE OF K IS INCREASED BY ONE. C C W IS A SCRATCH ARRAY. C C USE IS MADE OF ROUTINES FROM THE LIBRARY C OF BASIC LINEAR ALGEBRA SUBROUTINES (BLAS). C C PARAMETERS... C C INPUT/ C NAME TYPE OUTPUT/ SUB- DESCRIPTION C SCRATCH SCRIPTS C ------------------------------------------- C N INT. I NUMBER OF ROWS C C K INT. I/O NUMBER OF COLUMNS C C ZZ REAL I/O 2 SCALED ORTHOGONAL C MATRIX C C NZZR INT. I ROW DIMENSION OF ZZ C C DD REAL I/O 1 DIAGONAL SCALING C MATRIX (DIAGONAL C ELEMENTS ONLY) C C RR REAL I/O 1 RIGHT-TRIANGULAR C MATRIX IN COMPACT FORM. C C COL REAL I 1 COLUMN TO BE C ADDED TO RR C C FAIL LOG. O .TRUE. IF K,N C ARE IMPROPER C C W REAL SCR 1 WORKSPACE C ------------------------------------------- C C THE I-TH SEGMENT OF THE ARRAY RR IS N-I+2 SPACES C LONG AND CONTAINS 1 WORK SPACE FOLLOWED BY THE C K-I+1 ELEMENTS OF ROW I FOLLOWED BY N-K C SCRATCH SPACES. C *************** C C +++++++++++++++ C BLAS SCOPY,SDOT,SROTM,SROTMG C +++++++++++++++ C INTEGER I,J,JDEL,KP1,KP2 REAL DI,ONE,PARAM(5),WI,ZERO C REAL SDOT C C +++++++++++++++ C THE FOLLOWING DECLARATIONS ARE NECESSARY C FOR PORTABILITY WHEN SCOPY IS USED, AS C IT IS BELOW, TO FILL ARRAYS WITH A SINGLE VALUE C (ONE=UNITY AND ZERO=ZIP IN THIS CASE). C +++++++++++++++ C REAL UNITY(1),ZIP(1) EQUIVALENCE (ONE,UNITY(1)),(ZERO,ZIP(1)) C DATA ONE /1.0E+00/ DATA ZERO /0.0E+00/ C C ///////////////// BEGIN PROGRAM ////////////////// C IF (K .GE. 0 .AND. K .LT. N .AND. N .LE. NZZR) GO TO 5 FAIL = .TRUE. RETURN 5 CONTINUE IF (K .GT. 0) GO TO 20 C C *************** C FOR THE SPECIAL CASE THAT THE C FACTORIZATION WAS VACUOUS, C RESET THE ARRAYS ZZ AND DD C TO REPRESENT THE IDENTITY. C *************** C K = 0 CALL SCOPY(N,UNITY,0,DD,1) CALL SCOPY(N*N,ZIP,0,ZZ,1) DO 10 I=1,N ZZ(I,I) = ONE 10 CONTINUE 20 CONTINUE KP1 = K + 1 KP2 = K + 2 C C *************** C TRANSFORM THE INCOMING COLUMN, C AND STORE THE RESULT IN W. C *************** C DO 30 I=1,N W(I) = SDOT(N,ZZ(1,I),1,COL,1) 30 CONTINUE C C *************** C ZERO OUT THE SPIKE WHICH WOULD RESULT FROM C STORING W IN RR. UPDATE ZZ AND DD. C *************** C IF (KP2 .GT. N) GO TO 50 DO 40 I=KP2,N DI = DD(I) WI = W(I) CALL SROTMG(DD(KP1),DI,W(KP1),WI,PARAM) W(I) = WI DD(I) = DI CALL SROTM(N,ZZ(1,KP1),1,ZZ(1,I),1,PARAM) 40 CONTINUE 50 CONTINUE C C *************** C STORE THE NEW COLUMN, WHICH IS STILL C IN THE ARRAY W, INTO RR. C *************** C J = KP2 JDEL = N DO 60 I=1,KP1 RR(J) = W(I) J = J + JDEL JDEL = JDEL - 1 60 CONTINUE K = KP1 FAIL = .FALSE. RETURN END SUBROUTINE ZDRCOU(N,K,ZZ,NZZR,DD,RR,IC,FAIL) C INTEGER IC,K,N,NZZR LOGICAL FAIL REAL DD(1),RR(1),ZZ(NZZR,1) C C *************** C PREPARED BY RICHARD BARTELS C THE UNIVERSITY OF WATERLOO C COMPUTER SCIENCE DEPARTMENT C LATEST UPDATE .... 30 NOVEMBER, 1979. C C GIVEN THE FACTORIZATION C C ZZ*DD*RR C C OF SOME N BY K MATRIX C C (1 .LE. K .LE. N) C (N .GE. 1), C C WHERE C C (ZZ-TRANSP)*(ZZ) = (DD-INV), C DD IS DIAGONAL AND NONSINGULAR, C AND C RR HAS ZEROS BELOW THE DIAGONAL, C C AND GIVEN THE INDEX IC OF A COLUMN C TO BE REMOVED (1 .LE. IC .LE. K), C THIS PROGRAM UPDATES ZZ,DD AND RR . C C THE VALUE OF K IS DECREASED BY ONE, AND C THE COLUMN ORDERING IN RR IS CHANGED. C C USE IS MADE OF ROUTINES FROM THE LIBRARY C OF BASIC LINEAR ALGEBRA SUBROUTINES (BLAS). C C PARAMETERS... C C NAME TYPE INPUT/ SUB- DESCRIPTION C OUTPUT SCRIPTS C ------------------------------------------- C N INT. I NUMBER OF ROWS C C K INT. I/O NUMBER OF COLUMNS C C ZZ REAL I/O 2 SCALED ORTHOGONAL C MATRIX C C NZZR INT. I ROW DIMENSION OF ZZ C C DD REAL I/O 1 DIAGONAL SCALING C MATRIX (DIAGONAL C ELEMENTS ONLY) C C RR REAL I/O 1 RIGHT-TRIANGULAR C MATRIX IN COMPACT FORM. C C IC INT. I INDEX OF COLUMN C TO BE REMOVED C C FAIL LOG. O .TRUE. IF K,N,IC C ARE IMPROPER C ------------------------------------------- C C THE I-TH SEGMENT OF THE ARRAY RR IS N-I+2 SPACES C LONG AND CONTAINS 1 WORK SPACE FOLLOWED BY THE C K-I+1 ELEMENTS OF ROW I FOLLOWED BY N-K C SCRATCH SPACES. C *************** C C +++++++++++++++ C BLAS SROTM,SROTMG C +++++++++++++++ C INTEGER I,IM1,J,JEND,JINC,JSTRT,KM1,LSTRT REAL DI,PARAM(5),RJ C C ///////////////// BEGIN PROGRAM ////////////////// C IF (K .GE. 1 .AND. K .LE. N .AND. N .LE. NZZR) GO TO 5 FAIL = .TRUE. RETURN 5 CONTINUE KM1 = K - 1 C C *************** C SPECIAL CASES ARE HANDLED FIRST. C 1. K=1 AND THE FACTORIZATION BECOMES NULL. C 2. IC=K AND THE UPDATING IS TRIVIAL. C *************** C IF (K .GT. 1) GO TO 10 K = 0 FAIL = .FALSE. RETURN 10 CONTINUE IF (IC .LT. K) GO TO 20 K = KM1 FAIL = .FALSE. RETURN 20 CONTINUE C C *************** C GENERAL UPDATING STEP. C THE COLUMN TO BE DELETED MUST BE PERMUTED C TO THE RIGHT, AND SUBDIAGONAL ELEMENTS C WHICH RESULT IN RR HAVE TO BE C TRANSFORMED TO ZERO. C *************** C JSTRT = IC + 1 JEND = K JINC = N DO 50 I=1,K C C *************** C PERMUTATION OF THE I-TH ROW OF RR. C *************** C DO 30 J=JSTRT,JEND RR(J) = RR(J + 1) 30 CONTINUE IF (I .LE. IC) GO TO 40 C C *************** C TRANSFORMATION OF THE CURRENT AND LAST C ROWS (I AND I-1) OF RR AS WELL AS C CORRESPONDING CHANGES TO ZZ AND DD. C C THE EXTRA VARIABLES DI AND RJ C ARE USED TO AVOID AN ERROR MESSAGE C FROM THE PFORT VERIFIER, AND THEY C MAY BE REMOVED, IF DESIRED, SO THAT C THE CALL TO SROTMG WOULD BE C C CALL SROTMG(DD(IM1),DD(I),RR(LSTRT),RR(JSTRT),PARAM) C C *************** C IM1 = I - 1 DI = DD(I) RJ = RR(JSTRT) CALL SROTMG(DD(IM1),DI,RR(LSTRT),RJ,PARAM) RR(JSTRT) = RJ DD(I) = DI CALL SROTM(JEND - JSTRT + 1,RR(LSTRT + 1),1, * RR(JSTRT + 1),1,PARAM) CALL SROTM(N,ZZ(1,IM1),1,ZZ(1,I),1,PARAM) JSTRT = JSTRT + 1 40 CONTINUE C C *************** C INDEX UPDATING C *************** C LSTRT = JSTRT JSTRT = JSTRT + JINC JEND = JEND + JINC JINC = JINC - 1 50 CONTINUE K = KM1 FAIL = .FALSE. RETURN END SUBROUTINE ZDRGIT(N,K,ZZ,NZZR,RR,GV,SOL,FAIL,W) C INTEGER K,N,NZZR LOGICAL FAIL REAL GV(1),RR(1),W(1),SOL(1),ZZ(NZZR,1) C C *************** C PREPARED BY RICHARD BARTELS C THE UNIVERSITY OF WATERLOO C COMPUTER SCIENCE DEPARTMENT C LATEST UPDATE .... 30 NOVEMBER, 1979. C C GIVEN THE FACTORIZATION C C ZZ*DD*RR C C OF SOME N BY K MATRIX C C (1 .LE. K .LE. N) C (N .GE. 1), C C WHERE C C (ZZ-TRANSP)*(ZZ) = (DD-INV), C DD IS DIAGONAL AND NONSINGULAR, C AND C RR HAS ZEROS BELOW THE DIAGONAL, C C AND GIVEN AN ARBITRARY VECTOR GV OF C APPROPRIATE DIMENSION, THIS ROUTINE FINDS THE C VECTOR SOL SATISFYING THE UNDERDETERMINED SYSTEM C C (ZZ*DD*RR-TRANSP.)*(SOL) = (GV). C C THAT IS, C C (SOL) = ((ZZ*DD*RR)-GEN.INV.-TRANSP.)*(GV). C C THE ARRAY DD IS NOT NEEDED BY ZDRGIT. C C USE IS MADE OF ROUTINES FROM THE LIBRARY C OF BASIC LINEAR ALGEBRA SUBROUTINES (BLAS). C C W IS A SCRATCH ARRAY. C C PARAMETERS... C C INPUT/ C NAME TYPE OUTPUT/ SUB- DESCRIPTION C SCRATCH SCRIPTS C ------------------------------------------- C N INT. I NUMBER OF ROWS C C K INT. I/O NUMBER OF COLUMNS C C ZZ REAL I/O 2 SCALED ORTHOGONAL C MATRIX C C NZZR INT. I ROW DIMENSION OF ZZ C C RR REAL I/O 1 RIGHT-TRIANGULAR C MATRIX IN COMPACT FORM. C C GV REAL I 1 GIVEN VECTOR C C SOL REAL O 1 SOLUTION C C FAIL LOG. O .TRUE. IF N,K C ARE IMPROPER, OR IF C RR IS SINGULAR C C W REAL SCR 1 WORKSPACE C ------------------------------------------- C C THE I-TH SEGMENT OF THE ARRAY RR IS N-I+2 SPACES C LONG AND CONTAINS 1 WORK SPACE FOLLOWED BY THE C K-I+1 ELEMENTS OF ROW I FOLLOWED BY N-K C SCRATCH SPACES. C C IF GV AND SOL ARE DIMENSIONED TO THE C MAXIMUM OF N AND K , THEN THE SAME C STORAGE ARRAY MAY BE USED FOR BOTH OF C THESE VECTORS. C *************** C C +++++++++++++++ C SYSTEM ROUTINES ABS C C BLAS SAXPY,SCOPY C C BIG IS THE LARGEST POSITIVE NUMBER C WHICH CAN BE REPRESENTED IN THE C PRECISION OF THE ARITHMETIC BEING USED. C +++++++++++++++ C INTEGER I,J,JDEL REAL BIG,WI,ONE,RRJ,ZERO C C +++++++++++++++ C THE FOLLOWING DECLARATIONS ARE NECESSARY C FOR PORTABILITY WHEN SCOPY IS USED, AS C IT IS BELOW, TO FILL ARRAYS WITH A SINGLE VALUE C (ZERO=ZIP IN THIS CASE). C +++++++++++++++ C REAL ZIP(1) EQUIVALENCE (ZERO,ZIP(1)) C DATA BIG /1.0E+30/ DATA ONE /1.0E+00/ DATA ZERO /0.0E+00/ C C ///////////////// BEGIN PROGRAM ////////////////// C IF (K .GE. 1 .AND. K .LE. N .AND. N .LE. NZZR) GO TO 10 FAIL = .TRUE. RETURN 10 CONTINUE C C *************** C FIRST SOLVE (RR-TRANSP.)*(W) = (GV) C *************** C CALL SCOPY(K,GV,1,W,1) J = 2 JDEL = N + 1 DO 30 I=1,K RRJ = RR(J) WI = W(I) IF (ABS(RRJ) .GE. ONE) GO TO 20 IF (ABS(WI) .LT. ABS(RRJ)*BIG) GO TO 20 FAIL = .TRUE. RETURN 20 CONTINUE W(I) = WI/RRJ IF (I.LT.K) CALL SAXPY(K-I,(-W(I)),RR(J+1),1,W(I+1),1) J = J + JDEL JDEL = JDEL - 1 30 CONTINUE C C *************** C NOW (SOL) = (ZZ)*(W) C *************** C CALL SCOPY(N,ZIP,0,SOL,1) DO 40 I=1,K CALL SAXPY(N,W(I),ZZ(1,I),1,SOL,1) 40 CONTINUE FAIL = .FALSE. RETURN END SUBROUTINE ZDRGNV(N,K,ZZ,NZZR,RR,GV,SOL,FAIL) C INTEGER K,N,NZZR LOGICAL FAIL REAL GV(1),RR(1),SOL(1),ZZ(NZZR,1) C C *************** C PREPARED BY RICHARD BARTELS C THE UNIVERSITY OF WATERLOO C COMPUTER SCIENCE DEPARTMENT C LATEST UPDATE .... 30 NOVEMBER, 1979. C C GIVEN THE FACTORIZATION C C ZZ*DD*RR C C OF SOME N BY K MATRIX C C (1 .LE. K .LE. N) C (N .GE. 1), C C WHERE C C (ZZ-TRANSP)*(ZZ) = (DD-INV), C DD IS DIAGONAL AND NONSINGULAR, C AND C RR HAS ZEROS BELOW THE DIAGONAL, C C AND GIVEN AN ARBITRARY VECTOR GV OF C APPROPRIATE DIMENSION, THIS ROUTINE FINDS THE C VECTOR SOL GIVEN BY C C (SOL) = ((ZZ*DD*RR)-GEN.INV.)*(GV), C C WHICH REPRESENTS THE LEAST SQUARES PROBLEM C C (ZZ*DD*RR)*(SOL) = (GV). C C THE ARRAY DD IS NOT NEEDED BY ZDRGNV. C C USE IS MADE OF ROUTINES FROM THE LIBRARY C OF BASIC LINEAR ALGEBRA SUBROUTINES (BLAS). C C PARAMETERS... C C NAME TYPE INPUT/ SUB- DESCRIPTION C OUTPUT/ SCRIPTS C ------------------------------------------- C N INT. I NUMBER OF ROWS C C K INT. I/O NUMBER OF COLUMNS C C ZZ REAL I/O 2 SCALED ORTHOGONAL C MATRIX C C NZZR INT. I ROW DIMENSION OF ZZ C C RR REAL I/O 1 RIGHT-TRIANGULAR C MATRIX IN COMPACT FORM. C C GV REAL I 1 GIVEN VECTOR C C SOL REAL O 1 SOLUTION C C FAIL LOG. O .TRUE. IF N,K C ARE IMPROPER, OR IF C RR IS SINGULAR C ------------------------------------------- C C THE I-TH SEGMENT OF THE ARRAY RR IS N-I+2 SPACES C LONG AND CONTAINS 1 WORK SPACE FOLLOWED BY THE C K-I+1 ELEMENTS OF ROW I FOLLOWED BY N-K C SCRATCH SPACES. C *************** C C +++++++++++++++ C SYSTEM ROUTINES ABS C C BLAS SDOT C C BIG IS THE LARGEST POSITIVE NUMBER C WHICH CAN BE REPRESENTED IN THE C PRECISION OF THE ARITHMETIC BEING USED. C +++++++++++++++ C INTEGER I,IX,J,JDEL REAL BIG,ONE,TDEN,TNUM C REAL SDOT C DATA BIG /1.0E+30/ DATA ONE /1.0E+00/ C C ///////////////// BEGIN PROGRAM ////////////////// C IF (K .GE. 1 .AND. K .LE. N .AND. N .LE. NZZR) GO TO 10 FAIL = .TRUE. RETURN 10 CONTINUE C C *************** C FORM (V) = (ZZ(1)-TRANSP)*(GV), WHERE ZZ(1) C IS THE MATRIX OF THE FIRST K COLUMNS OF ZZ C C V CAN BE STORED IN THE ARRAY SOL. C *************** C DO 20 I=1,K SOL(I) = SDOT(N,ZZ(1,I),1,GV,1) 20 CONTINUE C C *************** C BACKSOLVE THE SYSTEM C (RR)*(SOL) = (V) C FOR THE VECTOR SOL C C NOTE THAT SOL AND V C ARE STORED IN THE SAME ARRAY. C *************** C J = (((N + 1)*(N + 2) - (N - K + 3)*(N - K + 2))/2) + 2 JDEL = N - K + 3 DO 40 IX=1,K I = K - IX + 1 TDEN = RR(J) TNUM = SOL(I) IF (IX.GT.1) TNUM = TNUM-SDOT(IX-1,RR(J+1),1,SOL(I+1),1) IF (ABS(TDEN) .GE. ONE) GO TO 30 IF (ABS(TNUM) .LT. ABS(TDEN)*BIG) GO TO 30 FAIL = .TRUE. RETURN 30 CONTINUE SOL(I) = TNUM/TDEN J = J - JDEL JDEL = JDEL + 1 40 CONTINUE FAIL = .FALSE. RETURN END SUBROUTINE ZDRPOC(N,K,ZZ,NZZR,DD,GV,POC,FAIL) C INTEGER K,N,NZZR LOGICAL FAIL REAL DD(1),POC(1),GV(1),ZZ(NZZR,1) C C *************** C PREPARED BY RICHARD BARTELS C THE UNIVERSITY OF WATERLOO C COMPUTER SCIENCE DEPARTMENT C LATEST UPDATE .... 30 NOVEMBER, 1979. C C ZZ IS AN N BY N (N .GE. 1) SCALED C ORTHOGONAL MATRIX. DD CONTAINS THE C DIAGONAL ELEMENTS OF A DIAGONAL SCALING C MATRIX. GV IS A GIVEN VECTOR OF LENGTH N. C C WE HAVE C C (ZZ-TRANSP.)*(ZZ) = (DD-INV.) C C AND C C ZZ*DD*RR = MAT C C FOR SOME N BY K (0 .LE. K .LE. N) C MATRIX RR WITH ZEROS BELOW THE DIAGONAL C AND SOME GIVEN MATRIX MAT. (NIETHER RR C NOR MAT ARE NEEDED BY ZDRPOC.) C C THEN C C (PROJ(OC)) = (ZZ(2))*(DD(2))*(ZZ(2)-TRANSP.) C C IS THE (ORTHOGONAL) PROJECTOR ON THE C COMPLEMENT OF THE RANGE SPACE OF MAT, C WHERE ZZ(2) REPRESENTS THE LAST N-K C COLUMNS OF ZZ AND DD(2) REPRESENTS THE C LOWER-RIGHT-HAND N-K ORDER SUBMATRIX OF DD. C C ZDRPOC PRODUCES THE VECTOR C C POC = (PROJ(OC))*GV . C C USE IS MADE OF ROUTINES FROM THE LIBRARY C OF BASIC LINEAR ALGEBRA SUBROUTINES (BLAS). C C PARAMETERS... C C INPUT/ C NAME TYPE OUTPUT/ SUB- DESCRIPTION C SCRATCH SCRIPTS C ------------------------------------------- C N INT. I ORDER OF ZZ,DD C C K INT. I/O NUMBER OF COLUMNS C OF ZZ DEFINING C RANGE OF MAT C C ZZ REAL I/O 2 SCALED ORTHOGONAL C MATRIX C C NZZR INT. I ROW DIMENSION OF ZZ C C DD REAL I/O 1 DIAGONAL SCALING C MATRIX (DIAGONAL C ELEMENTS ONLY) C C GV REAL I 1 VECTOR TO BE PROJECTED C C POC REAL O 1 PROJECTION C C FAIL LOG. O .TRUE. IF N,K C ARE IMPROPER C C ------------------------------------------- C C *************** C C +++++++++++++++ C BLAS SAXPY,SCOPY,SDOT C +++++++++++++++ C INTEGER I,KP1 REAL WI,ZERO C REAL SDOT C C +++++++++++++++ C THE FOLLOWING DECLARATIONS ARE NECESSARY C FOR PORTABILITY WHEN SCOPY IS USED, AS C IT IS BELOW, TO FILL ARRAYS WITH A SINGLE VALUE C (ZERO=ZIP IN THIS CASE). C +++++++++++++++ C REAL ZIP(1) EQUIVALENCE (ZERO,ZIP(1)) C DATA ZERO /0.0E+00/ C C ///////////////// BEGIN PROGRAM ////////////////// C KP1 = K + 1 IF (K .GE. 0 .AND. K .LE. N * .AND. N .GE. 1 .AND. N .LE. NZZR) GO TO 5 FAIL = .TRUE. RETURN 5 CONTINUE IF (K .GT. 0) GO TO 10 C C *************** C CASE 1 ... ZZ(2)=ZZ (K=0) C *************** C CALL SCOPY(N,GV,1,POC,1) FAIL = .FALSE. RETURN 10 CONTINUE IF (K .LT. N) GO TO 20 C C *************** C CASE 2 ... ZZ(2) IS VACUOUS (K=N) C *************** C CALL SCOPY(N,ZIP,0,POC,1) FAIL = .FALSE. RETURN 20 CONTINUE C C *************** C CASE 3 ... ZZ(2) IS INTERMEDIATE C BETWEEN THE OTHER TWO CASES C (0 .LT. K .LT. N) C *************** C CALL SCOPY(N,ZIP,0,POC,1) DO 30 I=KP1,N WI = SDOT(N,ZZ(1,I),1,GV,1)*DD(I) CALL SAXPY(N,WI,ZZ(1,I),1,POC,1) 30 CONTINUE FAIL = .FALSE. RETURN END C C --------------------------------------------- C ** APPENDIX ** C C BASIC LINEAR ALGEBRA SUBROUTINES (BLAS) -- C SASUM,SAXPY,SCOPY,SDOT, C SROTM,SROTMG,SSCAL C C THE BLAS ARE DESCRIBED IN C C BASIC LINEAR ALGEBRA SUBPROGRAMS C FOR FORTRAN USAGE C C BY C C. L. LAWSON, R. J. HANSON, C D. R. KINCAID, F. T. KROGH. C ACM TRANS. ON MATH SOFTWARE C VOL. 5, NO. 3, 308-325 (1979) C C THE BLAS ARE TAKEN AS C PRIMITIVES. THEIR CODE IS NOT C PART OF THE DEFINITION OF CL1. C THEY ARE INCLUDED FOR REFERENCE C ONLY. PLEASE REFER TO THE ABOVE C PUBLICATION FOR CURRENT INFORMATION. C C IF THE BLAS ARE AVAILABLE INDEPENDANTLY C ON THE TARGET MACHINE, THIS APPENDIX SHOULD C BE REMOVED FROM THE CODE FOR CL1. C --------------------------------------------- C REAL FUNCTION SASUM(N,SX,INCX) C INTEGER INCX,N REAL SX(1) C C RETURNS SUM OF MAGNITUDES OF SINGLE PRECISION SX. C INTEGER I,M,MP1,NS SASUM = 0.0E+00 IF(N.LE.0)RETURN IF(INCX.EQ.1)GOTO 20 C C CODE FOR INCREMENTS NOT EQUAL TO 1. C NS = N*INCX DO 10 I=1,NS,INCX SASUM = SASUM + ABS(SX(I)) 10 CONTINUE RETURN C C CODE FOR INCREMENTS EQUAL TO 1. C C C CLEAN-UP LOOP SO REMAINING VECTOR LENGTH IS A MULTIPLE OF 6. C 20 M = N - (N/6)*6 IF( M .EQ. 0 ) GO TO 40 DO 30 I = 1,M SASUM = SASUM + ABS(SX(I)) 30 CONTINUE IF( N .LT. 6 ) RETURN 40 MP1 = M + 1 DO 50 I = MP1,N,6 SASUM = SASUM + ABS(SX(I)) + ABS(SX(I + 1)) + ABS(SX(I + 2)) * + ABS(SX(I + 3)) + ABS(SX(I + 4)) + ABS(SX(I + 5)) 50 CONTINUE RETURN END SUBROUTINE SAXPY(N,SA,SX,INCX,SY,INCY) C INTEGER INCX,INCY,N REAL SA,SX(1),SY(1) C C OVERWRITE SINGLE PRECISION SY WITH SINGLE PRECISION SA*SX +SY. C INTEGER I,IX,IY,M,MP1,NS C IF(N.LE.0.OR.SA.EQ.0.E0) RETURN IF(INCX.EQ.INCY) IF(INCX-1) 5,20,60 5 CONTINUE C C CODE FOR NONEQUAL OR NONPOSITIVE INCREMENTS. 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 SY(IY) = SY(IY) + SA*SX(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 SO REMAINING VECTOR LENGTH IS A MULTIPLE OF 4. C 20 M = N - (N/4)*4 IF( M .EQ. 0 ) GO TO 40 DO 30 I = 1,M SY(I) = SY(I) + SA*SX(I) 30 CONTINUE IF( N .LT. 4 ) RETURN 40 MP1 = M + 1 DO 50 I = MP1,N,4 SY(I) = SY(I) + SA*SX(I) SY(I + 1) = SY(I + 1) + SA*SX(I + 1) SY(I + 2) = SY(I + 2) + SA*SX(I + 2) SY(I + 3) = SY(I + 3) + SA*SX(I + 3) 50 CONTINUE RETURN C C CODE FOR EQUAL, POSITIVE, NONUNIT INCREMENTS. C 60 CONTINUE NS = N*INCX DO 70 I=1,NS,INCX SY(I) = SA*SX(I) + SY(I) 70 CONTINUE RETURN END SUBROUTINE SCOPY(N,SX,INCX,SY,INCY) C INTEGER INCX,INCY,N REAL SX(1),SY(1) C C COPY SINGLE PRECISION SX TO SINGLE PRECISION SY. C INTEGER I,IX,IY,M,MP1,NS C IF(N.LE.0)RETURN IF(INCX.EQ.INCY) IF(INCX-1) 5,20,60 5 CONTINUE C C CODE FOR UNEQUAL OR NONPOSITIVE INCREMENTS. 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 SY(IY) = SX(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 SO REMAINING VECTOR LENGTH IS A MULTIPLE OF 7. C 20 M = N - (N/7)*7 IF( M .EQ. 0 ) GO TO 40 DO 30 I = 1,M SY(I) = SX(I) 30 CONTINUE IF( N .LT. 7 ) RETURN 40 MP1 = M + 1 DO 50 I = MP1,N,7 SY(I) = SX(I) SY(I + 1) = SX(I + 1) SY(I + 2) = SX(I + 2) SY(I + 3) = SX(I + 3) SY(I + 4) = SX(I + 4) SY(I + 5) = SX(I + 5) SY(I + 6) = SX(I + 6) 50 CONTINUE RETURN C C CODE FOR EQUAL, POSITIVE, NONUNIT INCREMENTS. C 60 CONTINUE NS = N*INCX DO 70 I=1,NS,INCX SY(I) = SX(I) 70 CONTINUE RETURN END REAL FUNCTION SDOT(N,SX,INCX,SY,INCY) C INTEGER INCX,INCY,N REAL SX(1),SY(1) C C RETURNS THE DOT PRODUCT OF SINGLE PRECISION SX AND SY. C INTEGER I,IX,IY,M,MP1,NS C SDOT = 0.0E+00 IF(N.LE.0)RETURN IF(INCX.EQ.INCY) IF(INCX-1)5,20,60 5 CONTINUE C C CODE FOR UNEQUAL INCREMENTS OR NONPOSITIVE INCREMENTS. 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 SDOT = SDOT + SX(IX)*SY(IY) IX = IX + INCX IY = IY + INCY 10 CONTINUE RETURN C C CODE FOR BOTH INCREMENTS EQUAL TO 1 C C C CLEAN-UP LOOP SO REMAINING VECTOR LENGTH IS A MULTIPLE OF 5. C 20 M = N - (N/5)*5 IF( M .EQ. 0 ) GO TO 40 DO 30 I = 1,M SDOT = SDOT + SX(I)*SY(I) 30 CONTINUE IF( N .LT. 5 ) RETURN 40 MP1 = M + 1 DO 50 I = MP1,N,5 SDOT = SDOT + SX(I)*SY(I) + SX(I + 1)*SY(I + 1) + * SX(I + 2)*SY(I + 2) + SX(I + 3)*SY(I + 3) + SX(I + 4)*SY(I + 4) 50 CONTINUE RETURN C C CODE FOR POSITIVE EQUAL INCREMENTS .NE.1. C 60 CONTINUE NS=N*INCX DO 70 I=1,NS,INCX SDOT = SDOT + SX(I)*SY(I) 70 CONTINUE RETURN END SUBROUTINE SROTM (N,SX,INCX,SY,INCY,SPARAM) C INTEGER INCX,INCY,N REAL SPARAM(5),SX(1),SY(1) C C APPLY THE MODIFIED GIVENS TRANSFORMATION, H, TO THE 2 BY N MATRIX C C (SX(1) SX(N)) C ( ... ) C (SY(1) SY(N)) C C WITH SPARAM(1)=SFLAG, H HAS ONE OF THE FOLLOWING FORMS.. C C SFLAG=-1.E0 SFLAG=0.E0 SFLAG=1.E0 SFLAG=-2.E0 C C (SH11 SH12) (1.E0 SH12) (SH11 1.E0) (1.E0 0.E0) C H=( ) ( ) ( ) ( ) C (SH21 SH22), (SH21 1.E0), (-1.E0 SH22), (0.E0 1.E0). C INTEGER I,KX,KY,NSTEPS REAL SFLAG,SH11,SH12,SH21,SH22,TWO,W,Z,ZERO C DATA ZERO,TWO /0.0E+00,2.0E+00/ C SFLAG=SPARAM(1) IF(N .LE. 0 .OR.(SFLAG+TWO.EQ.ZERO)) GO TO 140 IF(.NOT.(INCX.EQ.INCY.AND. INCX .GT.0)) GO TO 70 C NSTEPS=N*INCX IF(SFLAG) 50,10,30 10 CONTINUE SH12=SPARAM(4) SH21=SPARAM(3) DO 20 I=1,NSTEPS,INCX W=SX(I) Z=SY(I) SX(I)=W+Z*SH12 SY(I)=W*SH21+Z 20 CONTINUE GO TO 140 30 CONTINUE SH11=SPARAM(2) SH22=SPARAM(5) DO 40 I=1,NSTEPS,INCX W=SX(I) Z=SY(I) SX(I)=W*SH11+Z SY(I)=-W+SH22*Z 40 CONTINUE GO TO 140 50 CONTINUE SH11=SPARAM(2) SH12=SPARAM(4) SH21=SPARAM(3) SH22=SPARAM(5) DO 60 I=1,NSTEPS,INCX W=SX(I) Z=SY(I) SX(I)=W*SH11+Z*SH12 SY(I)=W*SH21+Z*SH22 60 CONTINUE GO TO 140 70 CONTINUE KX=1 KY=1 IF(INCX .LT. 0) KX=1+(1-N)*INCX IF(INCY .LT. 0) KY=1+(1-N)*INCY C IF(SFLAG)120,80,100 80 CONTINUE SH12=SPARAM(4) SH21=SPARAM(3) DO 90 I=1,N W=SX(KX) Z=SY(KY) SX(KX)=W+Z*SH12 SY(KY)=W*SH21+Z KX=KX+INCX KY=KY+INCY 90 CONTINUE GO TO 140 100 CONTINUE SH11=SPARAM(2) SH22=SPARAM(5) DO 110 I=1,N W=SX(KX) Z=SY(KY) SX(KX)=W*SH11+Z SY(KY)=-W+SH22*Z KX=KX+INCX KY=KY+INCY 110 CONTINUE GO TO 140 120 CONTINUE SH11=SPARAM(2) SH12=SPARAM(4) SH21=SPARAM(3) SH22=SPARAM(5) DO 130 I=1,N W=SX(KX) Z=SY(KY) SX(KX)=W*SH11+Z*SH12 SY(KY)=W*SH21+Z*SH22 KX=KX+INCX KY=KY+INCY 130 CONTINUE 140 CONTINUE RETURN END SUBROUTINE SROTMG (SD1,SD2,SX1,SY1,SPARAM) C REAL SD1,SD2,SPARAM(5),SX1,SY1 C C CONSTRUCT THE MODIFIED GIVENS TRANSFORMATION MATRIX H C WHICH ZEROS THE SECOND COMPONENT OF THE 2-VECTOR C (SQRT(SD1)*SX1,SQRT(SD2)*SY2)**T. C WITH SPARAM(1)=SFLAG, H HAS ONE OF THE FOLLOWING FORMS.. C C SFLAG=-1.E0 SFLAG=0.E0 SFLAG=1.E0 SFLAG=-2.E0 C C (SH11 SH12) (1.E0 SH12) (SH11 1.E0) (1.E0 0.E0) C H=( ) ( ) ( ) ( ) C (SH21 SH22), (SH21 1.E0), (-1.E0 SH22), (0.E0 1.E0). C C +++++++++++++++ C SYSTEM ROUTINES ABS C +++++++++++++++ C INTEGER IFLAG,IGO REAL GAM,GAMSQ,ONE,RGAM,RGAMSQ,SFLAG,SH11,SH12 REAL SH21,SH22,SP1,SP2,SQ1,SQ2,STEMP,SU,TWO,ZERO C DATA ZERO,ONE,TWO /0.0E+00,1.0E+00,2.0E+00/,IFLAG/1/ DATA GAM/4096.0E+00/ DATA GAMSQ /1.678E+07/ DATA RGAM /2.441E-04/ DATA RGAMSQ /5.960E-08/ C IF(.NOT. SD1 .LT. ZERO) GO TO 10 C GO ZERO-H-D-AND-SX1.. GO TO 60 10 CONTINUE C CASE-SD1-NONNEGATIVE SP2=SD2*SY1 IF(.NOT. SP2 .EQ. ZERO) GO TO 20 SFLAG=-TWO GO TO 260 C REGULAR-CASE.. 20 CONTINUE SP1=SD1*SX1 SQ2=SP2*SY1 SQ1=SP1*SX1 C IF(.NOT. ABS(SQ1) .GT. ABS(SQ2)) GO TO 40 SH21=-SY1/SX1 SH12=SP2/SP1 C SU=ONE-SH12*SH21 C IF(.NOT. SU .LE. ZERO) GO TO 30 C GO ZERO-H-D-AND-SX1.. GO TO 60 30 CONTINUE SFLAG=ZERO SD1=SD1/SU SD2=SD2/SU SX1=SX1*SU C GO SCALE-CHECK.. GO TO 100 40 CONTINUE IF(.NOT. SQ2 .LT. ZERO) GO TO 50 C GO ZERO-H-D-AND-SX1.. GO TO 60 50 CONTINUE SFLAG=ONE SH11=SP1/SP2 SH22=SX1/SY1 SU=ONE+SH11*SH22 STEMP=SD2/SU SD2=SD1/SU SD1=STEMP SX1=SY1*SU C GO SCALE-CHECK GO TO 100 C PROCEDURE..ZERO-H-D-AND-SX1.. 60 CONTINUE SFLAG=-ONE SH11=ZERO SH12=ZERO SH21=ZERO SH22=ZERO C SD1=ZERO SD2=ZERO SX1=ZERO C RETURN.. GO TO 220 C PROCEDURE..FIX-H.. 70 CONTINUE IF(.NOT. SFLAG .GE. ZERO) GO TO 90 C IF(.NOT. SFLAG .EQ. ZERO) GO TO 80 SH11=ONE SH22=ONE SFLAG=-ONE GO TO 90 80 CONTINUE SH21=-ONE SH12=ONE SFLAG=-ONE 90 CONTINUE GO TO IGO,(120,150,180,210) C PROCEDURE..SCALE-CHECK 100 CONTINUE 110 CONTINUE IF(.NOT. SD1 .LE. RGAMSQ) GO TO 130 IF(.NOT. IFLAG.EQ.1) GO TO 105 C C RECOMPUTE RESCALING PARAMETERS C MORE ACCURATELY.. C RGAM = ONE/GAM GAMSQ = GAM**2 RGAMSQ = RGAM**2 IFLAG = 2 105 CONTINUE IF(SD1 .EQ. ZERO) GO TO 160 ASSIGN 120 TO IGO C FIX-H.. GO TO 70 120 CONTINUE SD1=SD1*GAMSQ SX1=SX1*RGAM SH11=SH11*RGAM SH12=SH12*RGAM GO TO 110 130 CONTINUE 140 CONTINUE IF(.NOT. SD1 .GE. GAMSQ) GO TO 160 ASSIGN 150 TO IGO C FIX-H.. GO TO 70 150 CONTINUE SD1=SD1*RGAMSQ SX1=SX1*GAM SH11=SH11*GAM SH12=SH12*GAM GO TO 140 160 CONTINUE 170 CONTINUE IF(.NOT. ABS(SD2) .LE. RGAMSQ) GO TO 190 IF(SD2 .EQ. ZERO) GO TO 220 ASSIGN 180 TO IGO C FIX-H.. GO TO 70 180 CONTINUE SD2=SD2*GAMSQ SH21=SH21*RGAM SH22=SH22*RGAM GO TO 170 190 CONTINUE 200 CONTINUE IF(.NOT. ABS(SD2) .GE. GAMSQ) GO TO 220 ASSIGN 210 TO IGO C FIX-H.. GO TO 70 210 CONTINUE SD2=SD2*RGAMSQ SH21=SH21*GAM SH22=SH22*GAM GO TO 200 220 CONTINUE IF(SFLAG)250,230,240 230 CONTINUE SPARAM(3)=SH21 SPARAM(4)=SH12 GO TO 260 240 CONTINUE SPARAM(2)=SH11 SPARAM(5)=SH22 GO TO 260 250 CONTINUE SPARAM(2)=SH11 SPARAM(3)=SH21 SPARAM(4)=SH12 SPARAM(5)=SH22 260 CONTINUE SPARAM(1)=SFLAG RETURN END SUBROUTINE SSCAL(N,SA,SX,INCX) C INTEGER INCX,N REAL SA,SX(1) C C REPLACE SINGLE PRECISION SX BY SINGLE PRECISION SA*SX. C INTEGER I,M,MP1,NS C IF(N.LE.0)RETURN IF(INCX.EQ.1)GOTO 20 C C CODE FOR INCREMENTS NOT EQUAL TO 1. C NS = N*INCX DO 10 I = 1,NS,INCX SX(I) = SA*SX(I) 10 CONTINUE RETURN C C CODE FOR INCREMENTS EQUAL TO 1. C C C CLEAN-UP LOOP SO REMAINING VECTOR LENGTH IS A MULTIPLE OF 5. C 20 M = N - (N/5)*5 IF( M .EQ. 0 ) GO TO 40 DO 30 I = 1,M SX(I) = SA*SX(I) 30 CONTINUE IF( N .LT. 5 ) RETURN 40 MP1 = M + 1 DO 50 I = MP1,N,5 SX(I) = SA*SX(I) SX(I + 1) = SA*SX(I + 1) SX(I + 2) = SA*SX(I + 2) SX(I + 3) = SA*SX(I + 3) SX(I + 4) = SA*SX(I + 4) 50 CONTINUE RETURN END