SUBROUTINE CHRCNT(XMESS,NMESS) C***BEGIN PROLOGUE CHRCNT C***REFER TO SQED C***ROUTINES CALLED (NONE) C***END PROLOGUE CHRCNT C COUNT THE NUMBER OF NON-BLANK CHARACTERS IN XMESS. C RETURN THE RESULT IN NMESS. CHARACTER * (*) XMESS INTEGER NMESS DO 10 I = LEN(XMESS),1,-1 IF (XMESS(I:I).NE.' ') THEN NMESS = I RETURN END IF 10 CONTINUE NMESS = 0 RETURN END function d1mach ( i ) c*********************************************************************72 c cc D1MACH returns double precision real machine-dependent constants. c c Discussion: c c D1MACH can be used to obtain machine-dependent parameters c for the local machine environment. It is a function c with one input argument, and can be called as follows: c c D = D1MACH ( I ) c c where I=1,...,5. The output value of D above is c determined by the input value of I:. c c D1MACH ( 1) = B**(EMIN-1), the smallest positive magnitude. c D1MACH ( 2) = B**EMAX*(1 - B**(-T)), the largest magnitude. c D1MACH ( 3) = B**(-T), the smallest relative spacing. c D1MACH ( 4) = B**(1-T), the largest relative spacing. c D1MACH ( 5) = LOG10(B) c c Modified: c c 25 April 2007 c c Author: c c Phyllis Fox, Andrew Hall, Norman Schryer c c Reference: c c Phyllis Fox, Andrew Hall, Norman Schryer, c Algorithm 528: c Framework for a Portable Library, c ACM Transactions on Mathematical Software, c Volume 4, Number 2, June 1978, page 176-188. c c Parameters: c c Input, integer I, the index of the desired constant. c c Output, double precision D1MACH, the value of the constant. c implicit none double precision d1mach integer i if ( i < 1 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'D1MACH - Fatal error!' write ( *, '(a)' ) ' The input argument I is out of bounds.' write ( *, '(a)' ) ' Legal values satisfy 1 <= I <= 5.' write ( *, '(a,i12)' ) ' I = ', i d1mach = 0.0D+00 stop else if ( i == 1 ) then d1mach = 4.450147717014403D-308 else if ( i == 2 ) then d1mach = 8.988465674311579D+307 else if ( i == 3 ) then d1mach = 1.110223024625157D-016 else if ( i == 4 ) then d1mach = 2.220446049250313D-016 else if ( i == 5 ) then d1mach = 0.301029995663981D+000 else if ( 5 < i ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'D1MACH - Fatal error!' write ( *, '(a)' ) ' The input argument I is out of bounds.' write ( *, '(a)' ) ' Legal values satisfy 1 <= I <= 5.' write ( *, '(a,i12)' ) ' I = ', i d1mach = 0.0D+00 stop end if return end function dasum ( n, dx, incx ) c*********************************************************************72 c cc DASUM takes the sum of the absolute values. c c Discussion: c c This routine uses double precision real arithmetic. c c Modified: c c 08 July 2007 c c Author: c c Jack Dongarra c c Reference: c c Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart, c LINPACK User's Guide, c SIAM, 1979, c ISBN13: 978-0-898711-72-1, c LC: QA214.L56. c c Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh, c Basic Linear Algebra Subprograms for FORTRAN usage, c ACM Transactions on Mathematical Software, c Volume 5, Number 3, pages 308-323, 1979. c c Parameters: c c Input, integer N, the number of entries in the vector. c c Input, double precision X(*), the vector to be examined. c c Input, integer INCX, the increment between successive entries of X. c INCX must not be negative. c c Output, double precision DASUM, the sum of the absolute values of X. c implicit none double precision dasum double precision dtemp double precision dx(*) integer i integer incx integer m integer n integer nincx dasum = 0.0d0 dtemp = 0.0d0 if( n.le.0 .or. incx.le.0 )return if(incx.eq.1)go to 20 c c code for increment not equal to 1 c nincx = n*incx do i = 1,nincx,incx dtemp = dtemp + dabs(dx(i)) end do dasum = dtemp return c c code for increment equal to 1 c c c clean-up loop c 20 m = mod(n,6) if( m .eq. 0 ) go to 40 do i = 1,m dtemp = dtemp + dabs(dx(i)) end do if( n .lt. 6 ) go to 60 40 continue do i = m + 1, n, 6 dtemp = dtemp + dabs(dx(i)) + dabs(dx(i + 1)) + dabs(dx(i + 2)) & + dabs(dx(i + 3)) + dabs(dx(i + 4)) + dabs(dx(i + 5)) end do 60 dasum = dtemp return end subroutine daxpy ( n, da, dx, incx, dy, incy ) c*********************************************************************72 c cc DAXPY computes constant times a vector plus a vector. c c Discussion: c c This routine uses double precision real arithmetic. c c This routine uses unrolled loops for increments equal to one. c c Modified: c c 08 July 2007 c c Author: c c Jack Dongarra c c Reference: c c Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart, c LINPACK User's Guide, c SIAM, 1979, c ISBN13: 978-0-898711-72-1, c LC: QA214.L56. c c Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh, c Basic Linear Algebra Subprograms for FORTRAN usage, c ACM Transactions on Mathematical Software, c Volume 5, Number 3, pages 308-323, 1979. c c Parameters: c c Input, integer N, the number of elements in DX and DY. c c Input, double precision DA, the multiplier of DX. c c Input, double precision DX(*), the first vector. c c Input, integer INCX, the increment between successive entries of DX. c c Input/output, double precision DY(*), the second vector. c On output, DY(*) has been replaced by DY(*) + DA * DX(*). c c Input, integer INCY, the increment between successive entries of DY. c implicit none double precision da double precision dx(*) double precision dy(*) integer i integer incx integer incy integer ix integer iy integer m integer n if ( n .le. 0 ) then return end if if ( da .eq. 0.0d0 ) then return end if if(incx.eq.1.and.incy.eq.1)go to 20 c c code for unequal increments or equal increments c not equal to 1 c ix = 1 iy = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 if(incy.lt.0)iy = (-n+1)*incy + 1 do i = 1,n dy(iy) = dy(iy) + da*dx(ix) ix = ix + incx iy = iy + incy end do return c c code for both increments equal to 1 c c c clean-up loop c 20 m = mod(n,4) if( m .eq. 0 ) go to 40 do i = 1,m dy(i) = dy(i) + da*dx(i) end do if( n .lt. 4 ) return 40 continue do i = m + 1, 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) end do return end SUBROUTINE DBOCLS(W,MDW,MCON,MROWS,NCOLS,BL,BU,IND,IOPT,X,RNORMC, * RNORM,MODE,RW,IW) C***BEGIN PROLOGUE DBOCLS C***DATE WRITTEN 821220 (YYMMDD) C***REVISION DATE 880722 (YYMMDD) C***CATEGORY NO. K1A2A,G2E,G2H1,G2H2 C***KEYWORDS BOUNDS,CONSTRAINTS,INEQUALITY,LEAST SQUARES,LINEAR C***AUTHOR HANSON, R. J., SNLA C***PURPOSE Solve the bounded and constrained least squares C problem consisting of solving the equation C E*X = F (in the least squares sense) C subject to the linear constraints C C*X = Y. C***DESCRIPTION C C This subprogram solves the bounded and constrained least squares C problem. The problem statement is: C C Solve E*X = F (least squares sense), subject to constraints C C*X=Y. C C In this formulation both X and Y are unknowns, and both may C have bounds on any of their components. This formulation C of the problem allows the user to have equality and inequality C constraints as well as simple bounds on the solution components. C C This constrained linear least squares subprogram solves E*X=F C subject to C*X=Y, where E is MROWS by NCOLS, C is MCON by NCOLS. C C The user must have dimension statements of the form C C DIMENSION W(MDW,NCOLS+MCON+1), BL(NCOLS+MCON), BU(NCOLS+MCON), C * X(2*(NCOLS+MCON)+2+NX), RW(6*NCOLS+5*MCON) C INTEGER IND(NCOLS+MCON), IOPT(17+NI), IW(2*(NCOLS+MCON)) C C (here NX=number of extra locations required for the options; NX=0 C if no options are in use. Also NI=number of extra locations C for options 1-9. C C INPUT C ----- C C ------------------------- C W(MDW,*),MCON,MROWS,NCOLS C ------------------------- C The array W contains the (possibly null) matrix [C:*] followed by C [E:F]. This must be placed in W as follows: C [C : *] C W = [ ] C [E : F] C The (*) after C indicates that this data can be undefined. The C matrix [E:F] has MROWS rows and NCOLS+1 columns. The matrix C is C placed in the first MCON rows of W(*,*) while [E:F] C follows in rows MCON+1 through MCON+MROWS of W(*,*). The vector F C is placed in rows MCON+1 through MCON+MROWS, column NCOLS+1. The C values of MDW and NCOLS must be positive; the value of MCON must C be nonnegative. An exception to this occurs when using option 1 C for accumulation of blocks of equations. In that case MROWS is an C OUTPUT variable only, and the matrix data for [E:F] is placed in C W(*,*), one block of rows at a time. See IOPT(*) contents, option C number 1, for further details. The row dimension, MDW, of the C array W(*,*) must satisfy the inequality: C C If using option 1, C MDW .ge. MCON + max(max. number of C rows accumulated, NCOLS) C C If using option 8, MDW .ge. MCON + MROWS. C Else, MDW .ge. MCON + max(MROWS, NCOLS). C C Other values are errors, but this is checked only when using C option=2. The value of MROWS is an output parameter when C using option number 1 for accumulating large blocks of least C squares equations before solving the problem. C See IOPT(*) contents for details about option 1. C C ------------------ C BL(*),BU(*),IND(*) C ------------------ C These arrays contain the information about the bounds that the C solution values are to satisfy. The value of IND(J) tells the C type of bound and BL(J) and BU(J) give the explicit values for C the respective upper and lower bounds on the unknowns X and Y. C The first NVARS entries of IND(*), BL(*) and BU(*) specify C bounds on X; the next MCON entries specify bounds on Y. C C 1. For IND(J)=1, require X(J) .ge. BL(J); C IF J.gt.NCOLS, Y(J-NCOLS) .ge. BL(J). C (the value of BU(J) is not used.) C 2. For IND(J)=2, require X(J) .le. BU(J); C IF J.gt.NCOLS, Y(J-NCOLS) .le. BU(J). C (the value of BL(J) is not used.) C 3. For IND(J)=3, require X(J) .ge. BL(J) and C X(J) .le. BU(J); C IF J.gt.NCOLS, Y(J-NCOLS) .ge. BL(J) and C Y(J-NCOLS) .le. BU(J). C (to impose equality constraints have BL(J)=BU(J)= C constraining value.) C 4. For IND(J)=4, no bounds on X(J) or Y(J-NCOLS) are required. C (the values of BL(J) and BU(J) are not used.) C C Values other than 1,2,3 or 4 for IND(J) are errors. In the case C IND(J)=3 (upper and lower bounds) the condition BL(J) .gt. BU(J) C is an error. The values BL(J), BU(J), J .gt. NCOLS, will be C changed. Significant changes mean that the constraints are C infeasible. (Users must make this decision themselves.) C The new values for BL(J), BU(J), J .gt. NCOLS, define a C region such that the perturbed problem is feasible. If users C know that their problem is feasible, this step can be skipped C by using option number 8 described below. C C ------- C IOPT(*) C ------- C This is the array where the user can specify nonstandard options C for DBOCLS( ). Most of the time this feature can be ignored by C setting the input value IOPT(1)=99. Occasionally users may have C needs that require use of the following subprogram options. For C details about how to use the options see below: IOPT(*) CONTENTS. C C Option Number Brief Statement of Purpose C ------ ------ ----- --------- -- ------- C 1 Return to user for accumulation of blocks C of least squares equations. The values C of IOPT(*) are changed with this option. C The changes are updates to pointers for C placing the rows of equations into position C for processing. C 2 Check lengths of all arrays used in the C subprogram. C 3 Column scaling of the data matrix, [C]. C [E] C 4 User provides column scaling for matrix [C]. C [E] C 5 Provide option array to the low-level C subprogram DBOLS( ). C {Provide option array to the low-level C subprogram DBOLSM( ) by imbedding an C option array within the option array to C DBOLS(). Option 6 is now disabled.} C 7 Move the IOPT(*) processing pointer. C 8 Do not preprocess the constraints to C resolve infeasibilities. C 9 Do not pretriangularize the least squares matrix. C 99 No more options to change. C C ---- C X(*) C ---- C This array is used to pass data associated with options 4,5 and C 6. Ignore this parameter (on input) if no options are used. C Otherwise see below: IOPT(*) CONTENTS. C C C OUTPUT C ------ C C ----------------- C X(*),RNORMC,RNORM C ----------------- C The array X(*) contains a solution (if MODE .ge.0 or .eq.-22) for C the constrained least squares problem. The value RNORMC is the C minimum residual vector length for the constraints C*X - Y = 0. C The value RNORM is the minimum residual vector length for the C least squares equations. Normally RNORMC=0, but in the case of C inconsistent constraints this value will be nonzero. C The values of X are returned in the first NVARS entries of X(*). C The values of Y are returned in the last MCON entries of X(*). C C ---- C MODE C ---- C The sign of MODE determines whether the subprogram has completed C normally, or encountered an error condition or abnormal status. A C value of MODE .ge. 0 signifies that the subprogram has completed C normally. The value of mode (.ge. 0) is the number of variables C in an active status: not at a bound nor at the value zero, for C the case of free variables. A negative value of MODE will be one C of the cases (-57)-(-41), (-37)-(-22), (-19)-(-2). Values .lt. -1 C correspond to an abnormal completion of the subprogram. These C error messages are in groups for the subprograms DBOCLS(), C DBOLSM(), and DBOLS(). An approximate solution will be returned C to the user only when max. iterations is reached, MODE=-22. C C ----------- C RW(*),IW(*) C ----------- C These are working arrays. (normally the user can ignore the C contents of these arrays.) C C IOPT(*) CONTENTS C ------- -------- C The option array allows a user to modify some internal variables C in the subprogram without recompiling the source code. A central C goal of the initial software design was to do a good job for most C people. Thus the use of options will be restricted to a select C group of users. The processing of the option array proceeds as C follows: a pointer, here called LP, is initially set to the value C 1. At the pointer position the option number is extracted and C used for locating other information that allows for options to be C changed. The portion of the array IOPT(*) that is used for each C option is fixed; the user and the subprogram both know how many C locations are needed for each option. The value of LP is updated C for each option based on the amount of storage in IOPT(*) that is C required. A great deal of error checking is done by the C subprogram on the contents of the option array. Nevertheless it C is still possible to give the subprogram optional input that is C meaningless. For example option 4 uses the locations C X(NCOLS+IOFF),...,X(NCOLS+IOFF+NCOLS-1) for passing scaling data. C The user must manage the allocation of these locations. C C 1 C - C This option allows the user to solve problems with a large number C of rows compared to the number of variables. The idea is that the C subprogram returns to the user (perhaps many times) and receives C new least squares equations from the calling program unit. C Eventually the user signals "that's all" and a solution is then C computed. The value of MROWS is an output variable when this C option is used. Its value is always in the range 0 .le. MROWS C .le. NCOLS+1. It is the number of rows after the C triangularization of the entire set of equations. If LP is the C processing pointer for IOPT(*), the usage for the sequential C processing of blocks of equations is C C C IOPT(LP)=1 C Move block of equations to W(*,*) starting at C the first row of W(*,*). C IOPT(LP+3)=# of rows in the block; user defined C C The user now calls DBOCLS( ) in a loop. The value of IOPT(LP+1) C directs the user's action. The value of IOPT(LP+2) points to C where the subsequent rows are to be placed in W(*,*). Both of C these values are first defined in the subprogram. The user C changes the value of IOPT(LP+1) (to 2) as a signal that all of C the rows have been processed. C C C . abs ( B ) c = sign ( B ) if abs ( A ) <= abs ( B ); c c R = SIGMA * ( A * A + B * B ); c c C = A / R if R is not 0 c = 1 if R is 0; c c S = B / R if R is not 0, c 0 if R is 0. c c The computed numbers then satisfy the equation c c ( C S ) ( A ) = ( R ) c ( -S C ) ( B ) = ( 0 ) c c The routine also computes c c Z = S if abs ( A ) > abs ( B ), c = 1 / C if abs ( A ) <= abs ( B ) and C is not 0, c = 1 if C is 0. c c The single value Z encodes C and S, and hence the rotation: c c If Z = 1, set C = 0 and S = 1; c If abs ( Z ) < 1, set C = sqrt ( 1 - Z * Z ) and S = Z; c if abs ( Z ) > 1, set C = 1/ Z and S = sqrt ( 1 - C * C ); c c Modified: c c 07 July 2007 c c Author: c c Jack Dongarra c c Reference: c c Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart, c LINPACK User's Guide, c SIAM, 1979, c ISBN13: 978-0-898711-72-1, c LC: QA214.L56. c c Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh, c Basic Linear Algebra Subprograms for FORTRAN usage, c ACM Transactions on Mathematical Software, c Volume 5, Number 3, pages 308-323, 1979. c c Parameters: c c Input/output, double precision SA, SB. On input, SA and SB are the values c A and B. On output, SA is overwritten with R, and SB is c overwritten with Z. c c Output, double precision C, S, the cosine and sine of the c Givens rotation. c implicit none double precision da,db,c,s,roe,scale,r,z roe = db if( dabs(da) .gt. dabs(db) ) roe = da scale = dabs(da) + dabs(db) if( scale .eq. 0.0d0 ) then c = 1.0d0 s = 0.0d0 r = 0.0d0 z = 0.0d0 else r = scale*dsqrt((da/scale)**2 + (db/scale)**2) r = dsign(1.0d0,roe)*r c = da/r s = db/r z = 1.0d0 if( dabs(da) .gt. dabs(db) ) z = s if( dabs(db) .ge. dabs(da) .and. c .ne. 0.0d0 ) z = 1.0d0/c end if da = r db = z return end subroutine dscal ( n, da, dx, incx ) c*********************************************************************72 c cc DSCAL scales a vector by a constant. c c Discussion: c c This routine uses double precision real arithmetic. c c Modified: c c 07 July 2007 c c Author: c c Jack Dongarra c c Reference: c c Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart, c LINPACK User's Guide, c SIAM, 1979, c ISBN13: 978-0-898711-72-1, c LC: QA214.L56. c c Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh, c Basic Linear Algebra Subprograms for FORTRAN usage, c ACM Transactions on Mathematical Software, c Volume 5, Number 3, pages 308-323, 1979. c c Parameters: c c Input, integer N, the number of entries in the vector. c c Input, double precision SA, the multiplier. c c Input/output, double precision X(*), the vector to be scaled. c c Input, integer INCX, the increment between successive entries of X. c implicit none double precision da,dx(*) integer i,incx,m,n,nincx if( n.le.0 .or. incx.le.0 )return if(incx.eq.1)go to 20 c c code for increment not equal to 1 c nincx = n*incx do i = 1,nincx,incx dx(i) = da*dx(i) end do return c c code for increment equal to 1 c c c clean-up loop c 20 m = mod(n,5) if( m .eq. 0 ) go to 40 do i = 1,m dx(i) = da*dx(i) end do if( n .lt. 5 ) return 40 continue do i = m+1, 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) end do return end subroutine dswap ( n, dx, incx, dy, incy ) c*********************************************************************72 c cc DSWAP interchanges two vectors. c c Discussion: c c This routine uses double precision real arithmetic. c c Modified: c c 07 July 2007 c c Author: c c Jack Dongarra c c Reference: c c Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart, c LINPACK User's Guide, c SIAM, 1979, c ISBN13: 978-0-898711-72-1, c LC: QA214.L56. c c Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh, c Basic Linear Algebra Subprograms for FORTRAN usage, c ACM Transactions on Mathematical Software, c Volume 5, Number 3, pages 308-323, 1979. c c Parameters: c c Input, integer N, the number of entries in the vectors. c c Input/output, double precision X(*), one of the vectors to swap. c c Input, integer INCX, the increment between successive entries of X. c c Input/output, double precision Y(*), one of the vectors to swap. c c Input, integer INCY, the increment between successive elements of Y. c implicit none double precision dx(*),dy(*),dtemp integer i,incx,incy,ix,iy,m,n c if(n.le.0)return if(incx.eq.1.and.incy.eq.1)go to 20 c c code for unequal increments or equal increments not equal c to 1 c ix = 1 iy = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 if(incy.lt.0)iy = (-n+1)*incy + 1 do i = 1,n dtemp = dx(ix) dx(ix) = dy(iy) dy(iy) = dtemp ix = ix + incx iy = iy + incy end do return c c code for both increments equal to 1 c c c clean-up loop c 20 m = mod(n,3) if( m .eq. 0 ) go to 40 do i = 1,m dtemp = dx(i) dx(i) = dy(i) dy(i) = dtemp end do if( n .lt. 3 ) return 40 continue do i = m+1, n, 3 dtemp = dx(i) dx(i) = dy(i) dy(i) = dtemp dtemp = dx(i + 1) dx(i + 1) = dy(i + 1) dy(i + 1) = dtemp dtemp = dx(i + 2) dx(i + 2) = dy(i + 2) dy(i + 2) = dtemp end do return end SUBROUTINE DVOUT(N,DX,IFMT,IDIGIT) C REVISED FEB. 27, 1981. DOUBLE PRECISION DX(*) CHARACTER IFMT*(*) C C DOUBLE PRECISION VECTOR OUTPUT ROUTINE. C C INPUT.. C C N,DX(*) PRINT THE DOUBLE PRECISION ARRAY DX(I),I=1,...,N, ON C OUTPUT UNIT LOUT. THE HEADING IN THE FORTRAN FORMAT C STATEMENT IFMT(*), DESCRIBED BELOW, IS PRINTED AS A FIRST C STEP. THE COMPONENTS DX(I) ARE INDEXED, ON OUTPUT, C IN A PLEASANT FORMAT. C IFMT(*) A FORTRAN FORMAT STATEMENT. THIS IS PRINTED ON OUTPUT C UNIT LOUT WITH THE VARIABLE FORMAT FORTRAN STATEMENT C WRITE(LOUT,IFMT) C IDIGIT PRINT AT LEAST IABS(IDIGIT) DECIMAL DIGITS PER NUMBER. C THE SUBPROGRAM WILL CHOOSE THAT INTEGER 6,14,20 OR 28 C WHICH WILL PRINT AT LEAST IABS(IDIGIT) NUMBER OF C PLACES. IF IDIGIT.LT.0, 72 PRINTING COLUMNS ARE UTILIZED C TO WRITE EACH LINE OF OUTPUT OF THE ARRAY DX(*). (THIS C CAN BE USED ON MOST TIME-SHARING TERMINALS). IF C IDIGIT.GE.0, 133 PRINTING COLUMNS ARE UTILIZED. (THIS CAN C BE USED ON MOST LINE PRINTERS). C C EXAMPLE.. C C PRINT AN ARRAY CALLED (COSTS OF PURCHASES) OF LENGTH 100 SHOWING C 6 DECIMAL DIGITS PER NUMBER. THE USER IS RUNNING ON A TIME-SHARING C SYSTEM WITH A 72 COLUMN OUTPUT DEVICE. C C DOUBLE PRECISION COSTS(100) C N = 100 C IDIGIT = -6 C CALL DVOUT(N,COSTS,'(''1COSTS OF PURCHASES'')',IDIGIT) C C C C AUTHORS JOHN A. WISNIEWSKI SANDIA LABS ALBUQUERQUE. C RICHARD J. HANSON SANDIA LABS ALBUQUERQUE. C DATE JULY 27,1978. C C C GET THE UNIT NUMBER WHERE OUTPUTWILL BE WRITTEN. J=2 LOUT=I1MACH(J) WRITE(LOUT,IFMT) IF(N.LE.0) RETURN NDIGIT = IDIGIT IF(IDIGIT.EQ.0) NDIGIT = 6 IF(IDIGIT.GE.0) GO TO 80 C NDIGIT = -IDIGIT IF(NDIGIT.GT.6) GO TO 20 C DO 10 K1=1,N,4 K2 = MIN0(N,K1+3) WRITE(LOUT,1000) K1,K2,(DX(I),I=K1,K2) 10 CONTINUE RETURN C 20 CONTINUE IF(NDIGIT.GT.14) GO TO 40 C DO 30 K1=1,N,2 K2 = MIN0(N,K1+1) WRITE(LOUT,1001) K1,K2,(DX(I),I=K1,K2) 30 CONTINUE RETURN C 40 CONTINUE IF(NDIGIT.GT.20) GO TO 60 C DO 50 K1=1,N,2 K2=MIN0(N,K1+1) WRITE(LOUT,1002) K1,K2,(DX(I),I=K1,K2) 50 CONTINUE RETURN C 60 CONTINUE DO 70 K1=1,N K2 = K1 WRITE(LOUT,1003) K1,K2,(DX(I),I=K1,K2) 70 CONTINUE RETURN C 80 CONTINUE IF(NDIGIT.GT.6) GO TO 100 C DO 90 K1=1,N,8 K2 = MIN0(N,K1+7) WRITE(LOUT,1000) K1,K2,(DX(I),I=K1,K2) 90 CONTINUE RETURN C 100 CONTINUE IF(NDIGIT.GT.14) GO TO 120 C DO 110 K1=1,N,5 K2 = MIN0(N,K1+4) WRITE(LOUT,1001) K1,K2,(DX(I),I=K1,K2) 110 CONTINUE RETURN C 120 CONTINUE IF(NDIGIT.GT.20) GO TO 140 C DO 130 K1=1,N,4 K2 = MIN0(N,K1+3) WRITE(LOUT,1002) K1,K2,(DX(I),I=K1,K2) 130 CONTINUE RETURN C 140 CONTINUE DO 150 K1=1,N,3 K2 = MIN0(N,K1+2) WRITE(LOUT,1003) K1,K2,(DX(I),I=K1,K2) 150 CONTINUE RETURN 1000 FORMAT(1X,I4,' - ',I4,1X,1P8D14.5) 1001 FORMAT(1X,I4,' - ',I4,1X,1P5D22.13) 1002 FORMAT(1X,I4,' - ',I4,1X,1P4D28.19) 1003 FORMAT(1X,I4,' - ',I4,1X,1P3D36.27) END function i1mach ( i ) c*********************************************************************72 c cc I1MACH returns integer machine dependent constants. c c Discussion: c c Input/output unit numbers. c c I1MACH(1) = the standard input unit. c I1MACH(2) = the standard output unit. c I1MACH(3) = the standard punch unit. c I1MACH(4) = the standard error message unit. c c Words. c c I1MACH(5) = the number of bits per integer storage unit. c I1MACH(6) = the number of characters per integer storage unit. c c Integers. c c Assume integers are represented in the S digit base A form: c c Sign * (X(S-1)*A**(S-1) + ... + X(1)*A + X(0)) c c where 0 <= X(1:S-1) < A. c c I1MACH(7) = A, the base. c I1MACH(8) = S, the number of base A digits. c I1MACH(9) = A**S-1, the largest integer. c c Floating point numbers c c Assume floating point numbers are represented in the T digit c base B form: c c Sign * (B**E) * ((X(1)/B) + ... + (X(T)/B**T) ) c c where 0 <= X(I) < B for I=1 to T, 0 < X(1) and EMIN <= E <= EMAX. c c I1MACH(10) = B, the base. c c Single precision c c I1MACH(11) = T, the number of base B digits. c I1MACH(12) = EMIN, the smallest exponent E. c I1MACH(13) = EMAX, the largest exponent E. c c Double precision c c I1MACH(14) = T, the number of base B digits. c I1MACH(15) = EMIN, the smallest exponent E. c I1MACH(16) = EMAX, the largest exponent E. c c Modified: c c 25 April 2007 c c Author: c c Phyllis Fox, Andrew Hall, Norman Schryer c c Reference: c c Phyllis Fox, Andrew Hall, Norman Schryer, c Algorithm 528, c Framework for a Portable Library, c ACM Transactions on Mathematical Software, c Volume 4, Number 2, June 1978, page 176-188. c c Parameters: c c Input, integer I, chooses the parameter to be returned. c 1 <= I <= 16. c c Output, integer I1MACH, the value of the chosen parameter. c implicit none integer i integer i1mach if ( i < 1 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'I1MACH - Fatal error!' write ( *, '(a)' ) ' The input argument I is out of bounds.' write ( *, '(a)' ) ' Legal values satisfy 1 <= I <= 16.' write ( *, '(a,i12)' ) ' I = ', i i1mach = 0 stop else if ( i == 1 ) then i1mach = 5 else if ( i == 2 ) then i1mach = 6 else if ( i == 3 ) then i1mach = 7 else if ( i == 4 ) then i1mach = 6 else if ( i == 5 ) then i1mach = 32 else if ( i == 6 ) then i1mach = 4 else if ( i == 7 ) then i1mach = 2 else if ( i == 8 ) then i1mach = 31 else if ( i == 9 ) then i1mach = 2147483647 else if ( i == 10 ) then i1mach = 2 else if ( i == 11 ) then i1mach = 24 else if ( i == 12 ) then i1mach = -125 else if ( i == 13 ) then i1mach = 128 else if ( i == 14 ) then i1mach = 53 else if ( i == 15 ) then i1mach = -1021 else if ( i == 16 ) then i1mach = 1024 else if ( 16 < i ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'I1MACH - Fatal error!' write ( *, '(a)' ) ' The input argument I is out of bounds.' write ( *, '(a)' ) ' Legal values satisfy 1 <= I <= 16.' write ( *, '(a,i12)' ) ' I = ', i i1mach = 0 stop end if return end function idamax ( n, dx, incx ) c*********************************************************************72 c cc IDAMAX finds the index of element having maximum absolute value. c c Discussion: c c This routine uses double precision real arithmetic. c c Modified: c c 07 July 2007 c c Author: c c Jack Dongarra c c Reference: c c Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart, c LINPACK User's Guide, c SIAM, 1979, c ISBN13: 978-0-898711-72-1, c LC: QA214.L56. c c Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh, c Basic Linear Algebra Subprograms for FORTRAN usage, c ACM Transactions on Mathematical Software, c Volume 5, Number 3, pages 308-323, 1979. c c Parameters: c c Input, integer N, the number of entries in the vector. c c Input, double precision X(*), the vector to be examined. c c Input, integer INCX, the increment between successive entries of SX. c c Output, integer IDAMAX, the index of the element of SX of maximum c absolute value. c implicit none double precision dx(*),dmax integer idamax integer i,incx,ix,n idamax = 0 if( n.lt.1 .or. incx.le.0 ) return idamax = 1 if(n.eq.1)return if(incx.eq.1)go to 20 c c code for increment not equal to 1 c ix = 1 dmax = dabs(dx(1)) ix = ix + incx do i = 2,n if(dabs(dx(ix)).le.dmax) go to 5 idamax = i dmax = dabs(dx(ix)) 5 ix = ix + incx end do return c c code for increment equal to 1 c 20 dmax = dabs(dx(1)) do i = 2,n if( dmax .lt. dabs(dx(i)) ) then idamax = i dmax = dabs(dx(i)) end if end do return end SUBROUTINE IVOUT(N,IX,IFMT,IDIGIT) C REVISED FEB. 27, 1981. DIMENSION IX(N) CHARACTER IFMT*(*) C C INTEGER VECTOR OUTPUT ROUTINE. C C INPUT.. C C N,IX(*) PRINT THE INTEGER ARRAY IX(I),I=1,...,N, ON OUTPUT C UNIT LOUT. THE HEADING IN THE FORTRAN FORMAT C STATEMENT IFMT(*), DESCRIBED BELOW, IS PRINTED AS A FIRST C STEP. THE COMPONENTS IX(I) ARE INDEXED, ON OUTPUT, C IN A PLEASANT FORMAT. C IFMT(*) A FORTRAN FORMAT STATEMENT. THIS IS PRINTED ON OUTPUT C UNIT LOUT WITH THE VARIABLE FORMAT FORTRAN STATEMENT C WRITE(LOUT,IFMT) C IDIGIT PRINT UP TO IABS(IDIGIT) DECIMAL DIGITS PER NUMBER. C THE SUBPROGRAM WILL CHOOSE THAT INTEGER 4,6,10 OR 14 C WHICH WILL PRINT AT LEAST IABS(IDIGIT) NUMBER OF C PLACES. IF IDIGIT.LT.0, 72 PRINTING COLUMNS ARE UTILIZED C TO WRITE EACH LINE OF OUTPUT OF THE ARRAY IX(*). (THIS C CAN BE USED ON MOST TIME-SHARING TERMINALS). IF C IDIGIT.GE.0, 133 PRINTING COLUMNS ARE UTILIZED. (THIS CAN C BE USED ON MOST LINE PRINTERS). C C EXAMPLE.. C C PRINT AN ARRAY CALLED (COSTS OF PURCHASES) OF LENGTH 100 SHOWING C 6 DECIMAL DIGITS PER NUMBER. THE USER IS RUNNING ON A TIME-SHARING C SYSTEM WITH A 72 COLUMN OUTPUT DEVICE. C C DIMENSION ICOSTS(100) C N = 100 C IDIGIT = -6 C CALL IVOUT(N,ICOSTS,'(''1COSTS OF PURCHASES'')',IDIGIT) C C C C AUTHORS JOHN A. WISNIEWSKI SANDIA LABS ALBUQUERQUE. C RICHARD J. HANSON SANDIA LABS ALBUQUERQUE. C DATE JULY 27,1978. C C C GET THE UNIT NUMBER WHERE OUTPUTWILL BE WRITTEN. J=2 LOUT=I1MACH(J) WRITE(LOUT,IFMT) IF(N.LE.0) RETURN NDIGIT = IDIGIT IF(IDIGIT.EQ.0) NDIGIT = 4 IF(IDIGIT.GE.0) GO TO 80 C NDIGIT = -IDIGIT IF(NDIGIT.GT.4) GO TO 20 C DO 10 K1=1,N,10 K2 = MIN(N,K1+9) WRITE(LOUT,1000) K1,K2,(IX(I),I=K1,K2) 10 CONTINUE RETURN C 20 CONTINUE IF(NDIGIT.GT.6) GO TO 40 C DO 30 K1=1,N,7 K2 = MIN(N,K1+6) WRITE(LOUT,1001) K1,K2,(IX(I),I=K1,K2) 30 CONTINUE RETURN C 40 CONTINUE IF(NDIGIT.GT.10) GO TO 60 C DO 50 K1=1,N,5 K2=MIN(N,K1+4) WRITE(LOUT,1002) K1,K2,(IX(I),I=K1,K2) 50 CONTINUE RETURN C 60 CONTINUE DO 70 K1=1,N,3 K2 = MIN(N,K1+2) WRITE(LOUT,1003) K1,K2,(IX(I),I=K1,K2) 70 CONTINUE RETURN C 80 CONTINUE IF(NDIGIT.GT.4) GO TO 100 C DO 90 K1=1,N,20 K2 = MIN(N,K1+19) WRITE(LOUT,1000) K1,K2,(IX(I),I=K1,K2) 90 CONTINUE RETURN C 100 CONTINUE IF(NDIGIT.GT.6) GO TO 120 C DO 110 K1=1,N,15 K2 = MIN(N,K1+14) WRITE(LOUT,1001) K1,K2,(IX(I),I=K1,K2) 110 CONTINUE RETURN C 120 CONTINUE IF(NDIGIT.GT.10) GO TO 140 C DO 130 K1=1,N,10 K2 = MIN(N,K1+9) WRITE(LOUT,1002) K1,K2,(IX(I),I=K1,K2) 130 CONTINUE RETURN C 140 CONTINUE DO 150 K1=1,N,7 K2 = MIN(N,K1+6) WRITE(LOUT,1003) K1,K2,(IX(I),I=K1,K2) 150 CONTINUE RETURN 1000 FORMAT(1X,I4,' - ',I4,20(1X,I5)) 1001 FORMAT(1X,I4,' - ',I4,15(1X,I7)) 1002 FORMAT(1X,I4,' - ',I4,10(1X,I11)) 1003 FORMAT(1X,I4,' - ',I4,7(1X,I15)) END SUBROUTINE XERROR (XMESS, NMESS, NERR, LEVEL) C C REPLACEMENT SUBROUTINE FOR THE SLATEC XERROR PACKAGE SUBROUTINE. CHARACTER*120 XMESS IF (LEVEL.GE.1) THEN IERR=I1MACH(3) WRITE(IERR,'(1X,A)') XMESS(1:NMESS) WRITE(IERR,'('' ERROR NUMBER = '',I5,'', MESSAGE LEVEL = '',I5)') . NERR,LEVEL END IF RETURN END SUBROUTINE XERRWV(XMESS,NMESS,NERR,LEVEL,NI,I1,I2, . NR,R1,R2) C C REPLACEMENT SUBROUTINE FOR THE SLATEC XERRWV PACKAGE SUBROUTINE. CHARACTER*120 XMESS IF(LEVEL.LT.1) RETURN IERR=I1MACH(3) WRITE(IERR,'(1X,A)') XMESS(1:NMESS) WRITE(IERR,'('' ERROR NUMBER = '',I5,'', MESSAGE LEVEL = '',I5)') . NERR,LEVEL ASSIGN 80 TO IGOIPR GO TO (10,20),NI GO TO 80 10 CONTINUE C C WRITE ONE INTEGER VALUE. WRITE(IERR,'('' I1 = '',I8)') I1 GO TO IGOIPR 20 ASSIGN 30 TO IGOIPR GO TO 10 30 CONTINUE C C WRITE ONE INTEGER VALUE. WRITE(IERR,'('' I2 = '',I8)') I2 80 CONTINUE ASSIGN 40 TO IGOIPR GO TO (50,60),NR 40 RETURN 50 CONTINUE C C WRITE ONE REAL VALUE. WRITE(IERR,'('' R1 = '',1PE15.7)') R1 GO TO IGOIPR 60 ASSIGN 70 TO IGOIPR GO TO 50 70 CONTINUE C C WRITE REAL VALUE. WRITE(IERR,'('' R2 = '',1PE15.7)') R2 GO TO 40 END