program main !*****************************************************************************80 ! !! MAIN is the main program for MATALG. ! ! Discussion: ! ! MATALG is the main program for an interactive matrix algebra program. ! ! Journal: ! ! 10 January 2004 ! ! Replaced all > signs by < signs. ! ! 18 March 2003 ! ! Added HUGE_I, HUGE_R, the largest integer and real. ! ! 15 March 2001 ! ! Ditching IOUNIT, CHRWRT, paging, input files. ! ! 29 October 2000 ! ! Dropped routines to read from and write to example files. ! Easing out the IOUNIT variable, and the nice but cumbersome ! I/O handling. ! Dropped MAXFIX argument to COMRPN. ! ! 02 October 2000 ! ! Drop COSD, SIND, TAND. ! Add COT, SEC, CSC. ! Add internal variable "ANGLE_TO_RADIAN" and commands ! RADIAN and DEGREES that allow you to set the current angle measure. ! ! 20 June 2000 ! ! Converted to FORTRAN90 format. ! ! 30 April 2000 ! ! Renamed LUFACT to SGE_PLU. Deleted LUFAC. ! ! For some reason, I seem to need to make the first argument of ! COMRPN a variable rather than passing it in as a constant. ! I fear some variable is running out of bounds, and the inadequate ! bounds checking compiler options aren't helping me. ! ! If I try to enter "3,3" for the row and column dimensions, ! the routine gets confused, apparently trying to read the comma. ! ! 15 April 2000 ! ! FORTRAN90 changes: logical operator names, S_LENGTH=>LEN_TRIM, ! PARAMETER statement. ! ! 31 December 1999 ! ! A DO loop in RPNSET was zeroing out IRPN from 1 to MAXRPN. But ! when more than one formula is compiled, RPNSET only gets a portion ! of IRPN, but thinks it has it all, and thus the DO loop goes out ! of bounds. ! ! 15 April 1999 ! ! Renamed "RANDOM" to "UNIFORM_01_SAMPLE". ! ! 04 December 1998 ! ! Cleaned up R4MAT_POLY_CHAR, because I want to extract a copy for elsewhere. ! ! Added double exclamation markers to all routines. ! ! Version next ! ! Could you allow "max(V)" meaning the maximum element in the ! vector V? Could you allow "max(A)", and mean either the ! maximum in the matrix A, or a vector of row or column maxes, ! as in MATLAB? ! ! We need a better version of CHRREL, which prints the ! most information in the given space. ! ! Store input/output unit information in INDATA, rather than ! passing it all over. ! ! Write a FUNVEC routine to isolate vector and matrix functions, ! as was done for FUNSCL. ! ! Drop one letter commands, move to full names of commands. ! Then you can drop the "=" from a formula. ! ! Add rounding commands: ! ! ROUND2(A,ndigits) rounds A to so many digits. ! ROUND3(A,tol) drops all entries of magnitude less than TOL. ! ! Version 1.40 14 February 1996 ! ! I accidentally assumed INTSYM was of length MAXSYM rather than ! 80 in INICOM, causing illegal memory accesses. ! ! Had to shorten the header, which was getting to be more than ! 22 lines. ! ! Searching for the source of bugs on the Mac, which cause a ! simple command like "x = 1.5" to be mis-compiled! Instead of ! being converted to "X 1.5 = $" it gets converted to "X = 1.5 $". ! The DATA statements in COMPRN may not be properly treated by ! the Mac FORTRAN compiler, so I replaced them by assignment ! statements in INICOM, which fixed the problem. ! ! Added the "DEBUG" command, to force COMRPN to print internal data. ! ! Dropped the "A" command that merely announced the MATALG version. ! ! Added comments to every routine to list the arguments. ! ! Increased MAXLEN to 20. ! ! I corrected a problem with the naming of large constants, ! which caused MATALG to think that numbers like 1346.7 and ! 1329.3 were the same number since their representation was ! "13E4 in both cases. I increased the number of significant ! digits in the constant name (and that's why I had to increase ! MAXLEN to 20). ! ! Corrected an error in RELREA which referred to string entry ! LCHAR+1 without checking whether this was a legal reference. ! ! I corrected the rounding stuff, where the user could change ! RTOL but it wouldn't be used in the rounding. ! ! Version 1.39 12 October 1995 ! ! Corrected MAXLEN declarations. ! ! Version 1.38 09 February 1994 ! ! Adding the variable MAXLEN to control the maximum length ! of variable names, in case I want to increase them later. ! ! Added "KRYLOV" command to compute a given number of Krylov ! vectors. KRYLOV(A,X,N). ! ! Added "E" as a constant. ! ! Added code to stop users from altering constants. ! Original code was not effective in preventing this! ! ! I added a simple rounding command, ROUND(A) replaces all ! entries smaller than RTOL by zero. ! ! Trying to change LUFAC so that it uses SGEFA instead of my ! own routines. This is to avoid duplication, and to make it ! possible to get condition numbers, too. ! ! Since SGEFA/SGESL assume square matrix, you might find it ! useful to convert them to rectangular codes, so that you can ! dump your copies. ! ! Passed EPS to FUNVAL, so that the smallness of the determinant ! can be compared to EPS, rather than to 1.0e-6. ! ! Added the "A" command, so that the INPUT file can print ! out the MATALG version. ! ! Added source code for separate Frobenius and EISPACK matrix ! norm routines. Added "NORME" to MATALG list. ! ! Included notes in MATALG.HLP about making a double precision version. ! ! Dropped labels from all DO loops. ! ! Added GCD and LCM functions. ! ! Version 1.37 25 July 1993 ! ! Added "HILBERT" and "HILBINV". I had a bit of trouble ! with implicit dimensioning of the form "a=hilbert(5)", ! but it seems to have cleared up. ! ! I added "RCOND", the reciprocal condition number as computed ! by the LINPACK routine SGE_CO. ! ! I added routine NORM2, to allow computation of the two norm ! of a matrix. I haven't checked it yet. ! ! Added "NORMF", Frobenius norm of a matrix. ! ! Made "SEED" a variable that the user can access. ! ! Added "NORMS", spectral "norm" of a matrix. ! ! Documented the $ and % commands in the HELP program and ! in MATALG.HLP. ! ! Removed the SINE function since the SIN function is all we need. ! ! I changed "WRITE(*,..." to "WRITE(6,..." and "READ(*,..." to ! READ(5,..." so that input and output redirection could be used ! on PC's. As far as I know, the two statements are equivalent ! on other machines, except possibly on the Macintosh, so I'll ! have to check that. ! ! I removed the "N" and "C" commands, since the "#" command ! and the assignment operator are preferable. ! ! Added the new operator "\", which is used as in "A\b". This ! computes the LU factorization of A, and uses it to solve ! A*x=b, returning x. B may be a vector or matrix, that is, ! it may have one or several columns. ! This is apparently an operator available in MATLAB, and is ! preferable, because more accurate and stable and efficient, to ! computing INVERSE(A)*b. ! ! Inserted code for "<" input of file, but didn't check it. ! ! Reversed order of arguments in CHRWRT. ! ! Version 1.36 27 April 1993 ! ! Added (commented out) call to Language Systems routine ! "MoveOutWindow" which might help this output problem I've ! been seeing. ! ! Increased filename lengths to 80 characters. ! ! Renamed the variable NINT to NINTS to avoid conflict with the ! FORTRAN function of that name. ! ! Removed documentation for the "C" command. ! ! Documented the "#" command, which begins a comment. I ! hope to replace the "N" command by this, next version. ! ! CGC reported problems with ID, RAN and ZEROS functions. ! Fixed them. ! ! NORM0 and NORM1 can now handle matrices as well as vectors. ! ! Version 1.35 20 March 1993 ! ! CGC requested the ability to specify a particular label for ! data written to the example file. I added that in WRITEX. ! ! I corrected a problem in CHRREL, which previously stored a ! real number as a ratio of two integers. If those integers ! got too large, "wrap around" occurred, making nonsense of the ! results. ! ! MS FORTRAN had a fit, and wouldn't let me pass a real vector ! and use it as an integer, going from FUNVAL to EIGEN. So ! I stuck in a temporary vector called IWORK, of size 100. ! This is a minor sin, but fix it sometime! ! ! I thought the LS FORTRAN call "Call OutWindowScroll(9999)" ! would fix the scrolling problems I had been seeing, but I ! could not get it to work. However, the "-vax" compiler ! option fixed it right up. ! ! I made a Mac II version using the "-mc68020 -mc68881" ! compiler option. ! ! Reading from the example file is a real ugly piece of code. ! I corrected a mistake, and got it to at least print out ! the example once it's read in. But does it really have ! to share data with MATMAN? If not, it would be easier ! to type in examples by hand. I could even have the ! example file data have the same format as what you would ! type if you were entering the data, that is, ! ! e a 3,3 ! 1,2,3 ! 4,5,6 ! 7,8,9 ! ! The calculation of EPS cannot be trusted. It's really ! difficult to get a correct answer. I'm just setting it ! to 2**-23 for now! ! ! Got program to print "variable-name=" rather than "formula=' ! if formula has form "variable-name=formula". I did NOT do ! this for the case where the formula is re-evaluated via the ! V command. ! ! Language Systems FORTRAN for Macintosh seems to want a ! carriage control character, so I had to modify CHRWRT. ! ! Increased legal filename lengths to 60 characters. ! ! Added paging functions $ and %, using models from ! MATMAN. These are essential for doing big tests. ! ! Added # to make simple one line comments. ! ! INT function now calls FORTRAN INT function, rather than ! ANINT. The separate NINT function calls FORTRAN ANINT. ! ! Removed "ICOMP" argument from CHRWRT. ! ! Added LOG2 function. ! ! Extracted "FUNSCL" from FUNVAL, so that scalar code is separate. ! ! Routine names in error messages were suppressed, after ! complaints by reviewers that these were "frightening". ! ! Added COSD, SIND, TAND functions for degree-based ! calculations. ! ! Modified EIGEN so that, if only eigenvalues are wanted, then ! only HQR will be called. ! ! Inserted better random number generator. ! ! Changed C and E commands to print out matrix afterwards. ! ! I tried to compute 128*128*62, and got the answer ! 0.101581E+07, which seems pointless to me, since I had ! 14 places to print out. So I modified MATWRT so that ! if a quantity is less than 1.0e13 in magnitude, it's ! printed in I14 format. ! ! It was reported that complex eigenvectors are not computed ! correctly. This was caused by trying to back out the ! eigenvector transformations by hand. Instead, I now call ! RG and let it take care of everything. ! ! Changed the printout of the disk file name so that excess ! blanks were removed. ! ! Fixed the CHANGE command. It used to be that if you asked ! to change a nonexistent variable, the program would tell you ! it didn't exist and then ask you for its new value. ! ! The program will print out a warning if the inverse of a ! matrix with very small determinant is to be computed. ! ! Added HOUSE for householder matrix. ! ! Using "LEQI" for case insensitivity. ! ! It was requested that the "I" command be harder to invoke ! accidentally. So it's been changed to "INIT". In fact ! all the commands are now four or more letter keywords. ! This will make it possible to drop the "=" sign from ! the formula command shortly. ! ! Moved most of the main program into separate modules. ! ! Version 1.34 16 February 1992 ! ! Added a diagnostic to SYMADD, trying to track down why the ! Language Systems FORTRAN version of MATALG on the Macintosh ! thinks any new variable has already been entered. ! ! Version 1.33 08 December 1990 ! ! Made FACTOR more tolerant of input. It was failing if the ! user typed all the input on one line, for instance. ! ! Sign error in QRFACT, introduced in version 1.28, corrected. ! ! Version 1.32 12 October 1990 ! ! MAXROW increased to 12 from 10. ! ! CHRREA changed to allow caller to control capitalization of ! input. This means program no longer automatically ! capitalizes filenames, which was a bad idea! ! ! Version 1.31, 26 August 1990 ! ! The storage of the singular value decomposition factors was ! corrected. ! ! Version 1.31 08 August 1990 ! ! Argument EPS was not used by FUNVAL, MAXFIX not used by TOKENS. ! ! Storage of SVD items corrected. ! U is M by N, S and V are N by N. ! ! Version 1.30 10 July 1990 ! ! Internal coding changes. character*1 infix(maxfix) changed to ! character ( len = maxfix ) infix. ! ! CHRLOC replaced by call to INDEX. ! ! Version 1.29 10 June 1990 ! ! Before, you could say ! ! A(2,2)=3.0E+00 ! ! for instance, but not ! ! F(2,2)=3.0E+00 ! ! because the F confused the program. It thought the F was a ! command. Now, I've given up on the F command, since the "=" ! can do everything that it can do, and avoids the ambiguity. ! ! Version 1.28 16 May 1989 ! ! Modified QRFACT so that the R matrix output always has ! nonnegative diagonal. ! ! Added SVD decomposition A=U*S*TRANS(V) ! !*****************************************************************************80 ! ! Credits and Acknowledgements and Copyright: ! ! MATALG includes the text of routines from the LINPACK and ! EISPACK libraries. In some cases, these routines have ! been modified. ! ! This program was developed by ! ! John Burkardt, ! Department of Mathematics, ! Virginia Tech, ! Blacksburg, Virginia, 24060. ! ! Charles Cullen, ! Department of Mathematics and Statistics, ! University of Pittsburgh, ! Pittsburgh, Pennsylvania, 15260. ! ! EMAIL: ! ! burkardt@scs.fsu.edu ! ! This program is distributed free of charge with: ! ! Charles G. Cullen, ! An Introduction to Linear Algebra, Second Edition ! Harper Collins College Publishers, ! Glenview, Illinois, 1996 ! ! Development of this program was partially supported by a ! courseware development grant from the College of General ! Studies of the University of Pittsburgh. ! !*****************************************************************************80 ! ! MAXFIX specifies the maximum length of a formula that a user ! can enter, in characters. ! ! MAXROW specifies the maximum number of rows and columns that ! a variable can have. ! ! MAXRPN specifies the maximum number of symbols allowed in a compiled ! RPN formula. ! implicit none integer, parameter :: maxfix = 80 integer, parameter :: maxrow = 12 integer, parameter :: maxrpn = 80 integer, parameter :: maxlen = 20 character com character ( len = 3 ) command integer ierror integer ifrm integer iloc character ( len = maxfix ) infix integer irpn(maxrpn) logical s_eqi character ( len = 80 ) line character ( len = 80 ) line2 character ( len = maxlen ) namtmp character ( len = maxlen ) namvar integer ncol integer nform integer nrow character ( len = 80 ) prompt real value(maxrow,maxrow) ! ! Initialize ! call timestamp ( ) command = ' ' ierror = 0 ifrm = 0 ncol = 0 nrow = 0 nform = 0 namvar = ' ' com = 'I' call comrpn ( com, ierror, ifrm, infix, irpn, maxrow, maxrpn, namvar, & ncol, nrow, value ) ! ! Say hello. ! call hello ( maxrow ) ! ! Read next command. ! do line = ' ' if ( command(1:1) == '#' ) then command = ' ' prompt = ' ' else write ( *, '(a)' ) ' ' prompt = 'command (H for menu)' end if call chrrea ( line2, line, prompt, ierror, 0 ) ! ! If the input line was blank, don't bother processing it. ! Just go get the next line. ! if ( len_trim ( line2 ) <= 0 ) then command(1:1) = '#' cycle end if ! ! Check for "=" sign on a non-comment line. ! if ( line2(1:1) == '#' ) then iloc = 0 else iloc = index ( line2, '=' ) end if line = line2 if ( 0 < iloc ) then command = 'F' if ( iloc == 1 ) then line(1:1) = ' ' end if else prompt = 'command (H for menu)' call chrrea ( command, line, prompt, ierror, 1 ) end if if ( ierror /= 0 ) then ierror = 0 command = 'Q' end if ! ! DEBUG = print out current contents of COMRPN arrays. ! if ( s_eqi ( command(1:3), 'DEB' ) ) then com = 'L' namvar = 'DEBUG' call comrpn ( com, ierror, ifrm, infix, irpn, maxrow, maxrpn, namvar, & ncol, nrow, value ) ! ! DEGREES ! else if ( s_eqi ( command(1:3), 'DEG' ) ) then line = 'ANGLE_TO_RADIAN = PI / 180.0' call inform ( ierror, infix, irpn, line, & maxrow, maxrpn, ncol, nform, nrow, value ) ! ! ENTER = Enter variable, immediately request value too. ! else if ( s_eqi ( command(1:1), 'E' ) ) then call enter ( ierror, line, maxrow, maxrpn, namvar ) if ( ierror == 0 ) then com = 'L' call comrpn ( com, ierror, ifrm, infix, irpn, maxrow, maxrpn, namvar, & ncol, nrow, value ) end if ! ! F=Enter formula, interpret and evaluate. ! else if ( s_eqi ( command(1:1), 'F' ) ) then call inform ( ierror, infix, irpn, line, & maxrow, maxrpn, ncol, nform, nrow, value ) ! ! HELP = Help. ! else if ( s_eqi ( command(1:1), 'H' ) ) then call help ! ! INIT = Initialize. ! else if ( s_eqi ( command(1:1), 'I' ) ) then prompt = '"YES" to confirm you want to re-initialize.' line = ' ' call chrrea ( command, line, prompt, ierror, 1 ) if ( ierror /= 0 ) then cycle end if if ( s_eqi ( command(1:1), 'Y' ) ) then com = 'I' nform = 0 call comrpn ( com, ierror, ifrm, infix, irpn, maxrow, maxrpn, namvar, & ncol, nrow, value ) end if ! ! KILL = Kill variable. ! else if ( s_eqi ( command(1:1), 'K' ) ) then prompt = 'variable(s) to kill.' do call chrrea ( namvar, line, prompt, ierror, 3 ) if ( ierror /= 0 ) then exit end if if ( len_trim(namvar) <= 0 ) then exit end if com = 'D' call comrpn ( com, ierror, ifrm, infix, irpn, maxrow, maxrpn, namvar, & ncol, nrow, value ) if ( ierror /= 0 ) then exit end if if ( len_trim ( line ) <= 0 ) then exit end if end do ! ! Matrix factorization, 'LU', 'QR' or 'SVD'. ! else if ( s_eqi ( command(1:1), 'M' ) ) then call factor ( ierror, line, maxrow, maxrpn ) ! ! PART = partition. ! else if ( s_eqi ( command(1:1), 'U' ) ) then call part ( ierror, ifrm, line, maxrow, maxrpn, value ) ! ! POLY = Construct polynomial from root vector, ! or compute characteristic polynomial of matrix. ! else if ( s_eqi ( command(1:1), 'P' ) ) then call poly ( ierror, line, maxrow, maxrpn ) ! ! Q = Quit. ! QUit = QUIT NOW! ! else if ( s_eqi ( command(1:1), 'Q' ) ) then if ( s_eqi ( command(2:2), 'U' ) ) then command(1:1) = 'Y' else prompt = '"YES" to confirm you want to quit.' line = ' ' call chrrea ( command, line, prompt, ierror, 1 ) if ( ierror /= 0 ) then command(1:1) = 'Y' end if end if if ( s_eqi ( command(1:1), 'Y' ) ) then exit end if ! ! RADIANS ! else if ( s_eqi ( command(1:3), 'RAD' ) ) then line = 'ANGLE_TO_RADIAN = 1' call inform ( ierror, infix, irpn, line, & maxrow, maxrpn, ncol, nform, nrow, value ) ! ! TABLE = Tabulate. ! else if ( s_eqi ( command(1:1), 'L' ) ) then call table ( ierror, infix, irpn, line, maxrow, & maxrpn, ncol, nform, nrow, value ) ! ! TYPE = Type variables. ! else if ( s_eqi ( command(1:1), 'T' ) ) then if ( len_trim(line) == 0 ) then namvar = ' ' else prompt = 'variable to be listed.' call chrrea ( namvar, line, prompt, ierror, 3 ) if ( ierror /= 0 ) then cycle end if end if com = 'L' call comrpn ( com, ierror, ifrm, infix, irpn, maxrow, maxrpn, namvar, & ncol, nrow, value ) ! ! VALUE = Evaluate formula. ! else if ( s_eqi ( command(1:1), 'V' ) ) then if ( nform <= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'MATALG - Error!' write ( *, '(a)' ) ' Please enter a formula first!' ierror = 1 cycle end if com = 'E' ifrm = 1 call comrpn ( com, ierror, ifrm, infix, irpn, maxrow, maxrpn, namvar, & ncol, nrow, value ) if ( ierror /= 0 ) then cycle end if namtmp = 'Formula' call r4mat_print ( maxrow, nrow, ncol, value, namtmp ) ! ! # Comment. ! else if ( command(1:1) == '#' ) then ! ! Unknown command. ! else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'MATALG - Error!' write ( *, '(a)' ) ' Did not recognize the command:' write ( *, '(a)' ) '"' // trim ( command ) // '"' ierror = 1 end if ! ! If an error occurred, point that out. ! if ( ierror /= 0 ) then ierror = 0 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'MATALG noticed an error or mistake.' write ( *, '(a)' ) 'Your command was not carried out.' write ( *, '(a)' ) 'We are back at the main command menu.' end if end do write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'MATALG:' write ( *, '(a)' ) ' Normal end of execution.' stop end subroutine balanc ( nm, n, a, low, igh, scale ) !*****************************************************************************80 ! !! BALANC balances a real matrix before eigenvalue calculations. ! ! Discussion: ! ! This subroutine balances a real matrix and isolates eigenvalues ! whenever possible. ! ! Suppose that the principal submatrix in rows LOW through IGH ! has been balanced, that P(J) denotes the index interchanged ! with J during the permutation step, and that the elements ! of the diagonal matrix used are denoted by D(I,J). Then ! ! SCALE(J) = P(J), J = 1,...,LOW-1, ! = D(J,J), J = LOW,...,IGH, ! = P(J) J = IGH+1,...,N. ! ! The order in which the interchanges are made is N to IGH+1, ! then 1 to LOW-1. ! ! Note that 1 is returned for LOW if IGH is zero formally. ! ! Reference: ! ! J H Wilkinson and C Reinsch, ! Handbook for Automatic Computation, ! Volume II, Linear Algebra, Part 2, ! Springer Verlag, 1971. ! ! B Smith, J Boyle, J Dongarra, B Garbow, Y Ikebe, V Klema, C Moler, ! Matrix Eigensystem Routines, EISPACK Guide, ! Lecture Notes in Computer Science, Volume 6, ! Springer Verlag, 1976. ! ! Modified: ! ! 15 March 2001 ! ! Parameters: ! ! Input, integer NM, the leading dimension of A, which must ! be at least N. ! ! Input, integer N, the order of the matrix. ! ! Input/output, real A(NM,N), the N by N matrix. On output, ! the matrix has been balanced. ! ! Output, integer LOW, IGH, indicate that A(I,J) is equal to zero if ! (1) I is greater than J and ! (2) J=1,...,LOW-1 or I=IGH+1,...,N. ! ! Output, real SCALE(N), contains information determining the ! permutations and scaling factors used. ! implicit none integer nm integer n real a(nm,n) real b2 real c real f real g integer i integer iexc integer igh integer j integer k integer l integer low integer m logical noconv real r real, parameter :: radix = 16.0E+00 real s real scale(n) iexc = 0 j = 0 m = 0 b2 = radix**2 k = 1 l = n go to 100 20 continue scale(m) = j if ( j /= m ) then do i = 1, l call r4_swap ( a(i,j), a(i,m) ) end do do i = k, n call r4_swap ( a(j,i), a(m,i) ) end do end if if ( iexc == 2 ) go to 130 ! ! Search for rows isolating an eigenvalue and push them down. ! if ( l == 1 ) then low = k igh = l return end if l = l - 1 100 continue do j = l, 1, -1 do i = 1, l if ( i /= j ) then if ( a(j,i) /= 0.0E+00 ) then go to 120 end if end if end do m = l iexc = 1 go to 20 120 continue end do go to 140 ! ! Search for columns isolating an eigenvalue and push them left. ! 130 continue k = k + 1 140 continue do j = k, l do i = k, l if ( i /= j ) then if ( a(i,j) /= 0.0E+00 ) then go to 170 end if end if end do m = k iexc = 2 go to 20 170 continue end do ! ! Balance the submatrix in rows K to L. ! scale(k:l) = 1.0E+00 ! ! Iterative loop for norm reduction. ! noconv = .true. do while ( noconv ) noconv = .false. do i = k, l c = 0.0E+00 r = 0.0E+00 do j = k, l if ( j /= i ) then c = c + abs ( a(j,i) ) r = r + abs ( a(i,j) ) end if end do ! ! Guard against zero C or R due to underflow. ! if ( c /= 0.0E+00 .and. r /= 0.0E+00 ) then g = r / radix f = 1.0E+00 s = c + r do while ( c < g ) f = f * radix c = c * b2 end do g = r * radix do while ( g <= c ) f = f / radix c = c / b2 end do ! ! Balance. ! if ( ( c + r ) / f < 0.95E+00 * s ) then g = 1.0E+00 / f scale(i) = scale(i) * f noconv = .true. a(i,k:n) = a(i,k:n) * g a(1:l,i) = a(1:l,i) * f end if end if end do end do low = k igh = l return end subroutine balbak ( lda, n, low, igh, scale, m, z ) !*****************************************************************************80 ! !! BALBAK back transforms eigenvectors to undo the effect of BALANC. ! ! Discussion: ! ! This subroutine forms the eigenvectors of a real general ! matrix by back transforming those of the corresponding ! balanced matrix determined by BALANC. ! ! Reference: ! ! Parlett and Reinsch, ! Numerische Mathematik, ! Volume 13, pages 293-304, 1969. ! ! J H Wilkinson and C Reinsch, ! Handbook for Automatic Computation, ! Volume II, Linear Algebra, Part 2, ! Springer Verlag, 1971. ! ! B Smith, J Boyle, J Dongarra, B Garbow, Y Ikebe, V Klema, C Moler, ! Matrix Eigensystem Routines, EISPACK Guide, ! Lecture Notes in Computer Science, Volume 6, ! Springer Verlag, 1976. ! ! Modified: ! ! 18 February 2001 ! ! Parameters: ! ! Input, integer LDA, the leading dimension of Z. ! ! Input, integer N, the order of the matrix. ! ! Input, integer LOW, IGH, column indices determined by BALANC. ! ! Input, real SCALE(N), contains information determining ! the permutations and scaling factors used by BALANC. ! ! Input, integer M, the number of columns of Z to be back-transformed. ! ! Input/output, real Z(LDA,M), contains the real and imaginary parts ! of the eigenvectors, which, on return, have been back-transformed. ! implicit none integer lda integer m integer n integer i integer igh integer ii integer j integer k integer low real scale(n) real z(lda,m) if ( m <= 0 ) then return end if if ( igh /= low ) then do i = low, igh z(i,1:m) = scale(i) * z(i,1:m) end do end if do ii = 1, n i = ii if ( i < low .or. igh < i ) then if ( i < low ) then i = low - ii end if k = int ( scale(i) ) if ( k /= i ) then do j = 1, m call r4_swap ( z(i,j), z(k,j) ) end do end if end if end do return end subroutine cdiv ( ar, ai, br, bi, cr, ci ) !*****************************************************************************80 ! !! CDIV carries out complex division. ! ! Discussion: ! ! CDIV computes: ! ! (CR,CI) = (AR,AI) / (BR,BI) ! ! using real arithmetic. ! ! Modified: ! ! 15 March 2001 ! ! Parameters: ! ! Input, real AR, AI, the real and imaginary parts of the ! number to be divided. ! ! Input, real BR, BI, the real and imaginary parts of the divisor. ! ! Output, real CR, CI, the real and imaginary parts of the resultant. ! implicit none real ai real ais real ar real ars real bi real bis real br real brs real ci real cr real s s = abs ( br ) + abs ( bi ) ars = ar / s ais = ai / s brs = br / s bis = bi / s s = brs**2 + bis**2 cr = ( ars * brs + ais * bis ) / s ci = ( ais * brs - ars * bis ) / s return end subroutine ch_cap ( c ) !*****************************************************************************80 ! !! CH_CAP capitalizes a single character. ! ! Modified: ! ! 19 July 1998 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, character C, the character to capitalize. ! implicit none character c integer itemp itemp = ichar ( c ) if ( 97 <= itemp .and. itemp <= 122 ) then c = char ( itemp - 32 ) end if return end function ch_eqi ( c1, c2 ) !*****************************************************************************80 ! !! CH_EQI is a case insensitive comparison of two characters for equality. ! ! Examples: ! ! C_EQI ( 'A', 'a' ) is .TRUE. ! ! Modified: ! ! 14 August 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character C1, C2, the characters to compare. ! ! Output, logical CH_EQI, the result of the comparison. ! implicit none logical ch_eqi character c1 character c2 character cc1 character cc2 cc1 = c1 cc2 = c2 call ch_cap ( cc1 ) call ch_cap ( cc2 ) if ( cc1 == cc2 ) then ch_eqi = .true. else ch_eqi = .false. end if return end function ch_is_alpha ( c ) !*****************************************************************************80 ! !! CH_IS_ALPHA returns TRUE if C is an alphabetic character. ! ! ! Modified: ! ! 05 August 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character C, a character to check. ! ! Output, logical CH_IS_ALPHA is TRUE if C is an alphabetic character. ! implicit none character c logical ch_is_alpha if ( ( lle ( 'a', c ) .and. lle ( c, 'z' ) ) .or. & ( lle ( 'A', c ) .and. lle ( c, 'Z' ) ) ) then ch_is_alpha = .true. else ch_is_alpha = .false. end if return end subroutine ch_to_digit ( c, digit ) !*****************************************************************************80 ! !! CH_TO_DIGIT returns the integer value of a base 10 digit. ! ! Example: ! ! C DIGIT ! --- ----- ! '0' 0 ! '1' 1 ! ... ... ! '9' 9 ! ' ' 0 ! 'X' -1 ! ! Modified: ! ! 04 August 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character C, the decimal digit, '0' through '9' or blank ! are legal. ! ! Output, integer DIGIT, the corresponding integer value. If C was ! 'illegal', then DIGIT is -1. ! implicit none character c integer digit if ( lge ( c, '0' ) .and. lle ( c, '9' ) ) then digit = ichar ( c ) - 48 else if ( c == ' ' ) then digit = 0 else digit = -1 end if return end subroutine chrinp ( ierror, line, prompt ) !*****************************************************************************80 ! !! CHRINP checks for more input from the buffer, or prints a prompt. ! ! Discussion: ! ! The routine checks to see whether there is any more information in ! the buffer array LINE. If so, it simply updates the prompt ! and returns. Otherwise, it prints the prompt string out, ! reads the input from the user, and reprints the prompt and ! the user input on those I/O units where it is appropriate. ! ! Modified: ! ! 15 March 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, integer IERROR. ! ! If IERROR is nonzero on input, CHRINP stops. It is the ! calling routine's responsibility to make sure IERROR is ! zero on input. This is because CHRINP signals problems ! to the calling routine using IERROR. If the routine ! does not take the trouble to reset IERROR, then it ! is likely not to have addressed the problem itself. ! These problems can include things like end of input, ! so a failure to act can be catastrophic. ! ! On output, IERROR = ! 0, if no errors were detected, ! 1, if there was an error in the read, ! 2, if there was an end-of-file in the read. ! ! Input/output, character ( len = 80 ) LINE. ! ! On input, LINE may contain information that the calling ! program can use, or LINE may be empty. ! ! On output, LINE is unchanged if it contained information ! on input. But if the input LINE was empty, then the ! output LINE contains whatever information the user typed. ! ! Input/output, character ( len = 80 ) PROMPT, used for the prompt string. ! implicit none integer icomma integer ierror integer ios character ( len = 80 ) line character ( len = 80 ) prompt ! ! Catch nasty errors in calling routines. ! if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'CHRINP - Error!' write ( *, '(a,i8)' ) ' Nonzero input value of IERROR = ', ierror stop end if ! ! If there is nothing in the LINE buffer, then: ! print the prompt line, ! read the input line, ! remove double blanks. ! if ( len_trim ( line ) <= 0 ) then write ( *, '(a)' ) 'Enter ' // trim ( prompt ) read ( *, '(a)', iostat = ios ) line if ( ios /= 0 ) then ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'CHRINP - Fatal Error!' write ( *, '(a)' ) ' Illegal input:' write ( *, '(a)' ) trim ( line ) stop end if call s_blanks_delete ( line ) end if ! ! If an item was read, remove that item from prompt list. ! if ( 0 < len_trim ( line ) ) then icomma = index ( prompt, ',' ) if ( 0 < icomma .and. icomma < 80 .and. & prompt(icomma+1:icomma+1) == ' ' ) then icomma = icomma + 1 end if call s_chop ( prompt, 1, icomma ) end if return end subroutine chrrea ( string, line, prompt, ierror, iterm ) !*****************************************************************************80 ! !! CHRREA controls the process of getting input from the user. ! ! Discussion: ! ! CHRREA accepts LINE, which is assumed to contain user ! input characters and a prompt line. ! ! If the line is blank, the prompt is printed and user input into LINE. ! ! Then characters are read from LINE into STRING, and removed ! from LINE. ! ! PROMPT is also updated. On satisfactory input of string, ! everything in PROMPT up to and including the first comma is ! removed. ! ! Modified: ! ! 15 March 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, character ( len = * ) STRING, the string extracted from the input. ! ! Input/output, character ( len = 80 ) LINE, a buffer containing ! unprocessed user input. ! ! Input/output, character ( len = 80 ) PROMPT, used for the prompt string. ! ! Output, integer IERROR, error flag. ! 0, no error occurred. ! 1, Format error during read. ! 2, End of file during read. ! ! Input, integer ITERM, tells the program what rules are used to ! decide if input has terminated. ! ! 0, no check for terminators. ! 1, Blank, slash, comma, and equal terminate input. ! 2, nonalphabetic terminates input ! 3, nonalphanumeric terminates input ! implicit none character chrtmp integer i integer ierror integer iterm integer lchar character ( len = 80 ) line integer nchar logical num character ( len = 80 ) prompt logical s_is_alpha character ( len = * ) string ierror = 0 string = ' ' call chrinp ( ierror, line, prompt ) if ( ierror /= 0 ) then return end if ! ! Remove double blanks. ! if ( iterm == 2 .or. iterm == 3 ) then call s_blanks_delete ( line ) end if lchar = 0 nchar = len ( string ) do i = 1, nchar if ( lchar /= 0 ) then exit end if chrtmp = line(i:i) if ( iterm == 1 ) then if ( chrtmp == ' ' .or. chrtmp == '/' .or. chrtmp == ',' .or. & chrtmp == '=' ) then lchar = i end if else if ( iterm == 2 ) then if ( .not. s_is_alpha ( chrtmp ) ) then lchar = i end if else if ( iterm == 3 ) then num = lge ( chrtmp, '0' ) .and. lle ( chrtmp, '9' ) if ( ( .not. s_is_alpha ( chrtmp ) ) .and. ( .not. num ) ) then lchar = i end if end if if ( lchar == 0 ) then string(i:i) = chrtmp line(i:i) = ' ' end if end do if ( lchar == 0 ) then lchar = nchar end if ! ! Chop the extracted string from the input line. ! call s_chop ( line, 1, lchar ) return end subroutine comrpn ( com, ierror, ifrm, infix, irpn, maxrow, maxrpn, & namvar, ncol, nrow, value ) !*****************************************************************************80 ! !! COMRPN can translate formulas you type in and evaluate them. ! ! Discussion: ! ! This means that you can design interactive FORTRAN programs ! which can input their formulas at run time, rather than ! recompiling every time you want a new formula. ! ! If you declare vectors or matrices, you may reference entries ! such as X(3) or A(3,2) or even X(2+3). ! ! Formulas are made from operators, constants, punctuation, ! variables, and functions. ! ! Arithmetic Operators: ! ! The list of arithmetic operators includes: ! ! + - * / ** ^ = \ ! ! * means multiplication, and is standard matrix multiplication if ! both arguments are (conformable) matrices. If one argument is a ! scalar, it multiplies all the entries of the other argument. ! If both arguments are row or column vectors, the dot product ! is taken. ! ! Thus, if we assume that ! ! X is a scalar equal to 2, ! Y is a scalar equal to 3, ! U is a vector equal to (1, 2, 3), ! V is a vector equal to (1, 0, 2), ! ! then ! ! 2*3 results in 6. ! ! X*Y results in 6. ! ! U*V results in 7. ! ! X*V would result in the vector (2, 0, 4). ! ! / means division, as in A/B, but for matrices the only allowed ! form requires that B be a scalar, in which case each element of ! A is divided by B. ! ! + and - are allowed for pairs of scalars, vectors or matrices of ! the same order, and also for a square matrix and a scalar, ! in which case A+2 is interpreted as A+2*I. ! ! ** or ^ means exponentiation. In the scalar case, any ! nonnegative scalar can be taken to any power, positive, zero, or ! negative. A negative scalar may only be taken to an integer power. ! ! ** is also supported for square matrices, which can be taken ! to any nonnegative integer power. Nonsingular matrices may ! also be taken to negative integer powers. ! ! = is used to assign values. For vectors or matrices, this may ! also be used to assign a single entry, as in A(3,2)=7. ! ! / is standard scalar division. It is used in a formula like ! ! X/Y ! ! or ! ! (A*B)/(X+1) ! ! Normally, X would be a scalar quantity. However, as NONSTANDARD ! shorthand, you may 'divide' a matrix by a scalar, in which case ! each entry of the matrix is divided by the scalar, e.g. ! ! A/2 = (1/2)*A. ! ! MATALG will NOT allow you to use the "/" operator to represent ! multiplication by an inverse matrix. Thus, if A is a matrix, ! the statement ! ! B/A ! ! is illegal. However, the statement ! ! INV(A)*B ! ! will compute the inverse of A and multiply by B, and the statement ! ! A \ B ! ! will compute the LU factorization of A, and use that to ! compute the inverse of A times B. ! ! Constants: ! ! Real and integer constants may appear in your formulas, as well ! as the symbol 'PI' and the machine unit roundoff number 'EPS' ! which is the power of two such that 1+EPS>1 but 1+(EPS/2) = 1. ! ! Punctuation: ! ! Blanks may be used anywhere. ! ! Commas separate arguments in a function, such as MAX(X,Y). ! ! Parentheses are used: ! ! to group quantities: (X+Y)*Z ! ! to reference an entry of a vector or matrix: A(2,2) ! ! to mark the argument of a function: SIN(X), MAX(X,Y) ! ! ! Named Functions: ! ! S, S1, S2: Arguments which may only be scalar. ! V : Arguments may only be a scalar or vector. ! M : Arguments may only be scalar or square matrix. ! MV : Arguments may only be scalar, vector, or square matrix. ! * : Arguments which may be scalar, vector, or matrix. ! I : Arguments which may only be positive integers. ! ! ABS(*) Absolute value of *. ! ! ACOS(S) The arc cosine of S. ! -1 <= S < = 1. ! ! ASIN(S) The arc sine of S. ! -1 <= S < = 1. ! ! ATAN(S) The arc tangent of S. ! ! ATAN2(S1,S2) The arc tangent of (S1/S2) ! Correctly computes ATAN2(1.0,0.0)=PI/2. ! ! COS(S) The cosine of angle S. ! ! COT(S) The cotangent of angle S. ! ! COSH(S) The hyperbolic cosine of angle S. ! ! CSC(S) The cosecant of angle S. ! ! DET(M) Determinant of square matrix M. ! ! DIAG(*) The diagonal entries of *, stored in a vector. ! ! E The base of the natural logarithm system. ! E = 2.71828... ! You may not change the value of E. ! ! EPS The machine epsilon, or unit roundoff number. ! EPS is the power of 2 such that ! ! 1 < 1+EPS and 1 = 1+(EPS/2). ! ! You may not change the value of EPS. ! ! EVAL(M) Real and imaginary parts of eigenvalues of matrix ! M, stored in a matrix of N rows and 2 columns, ! with the real parts in column 1, and the ! imaginary parts in column 2. ! ! EVEC(M) Eigenvector matrix of a square matrix M. ! ! (Eigenvectors corresponding to a complex pair: ! If eigenvalues j and j+1 are a complex pair, ! then the eigenvector for eigenvalue j is ! column j + i*column j+1, and the eigenvector ! for eigenvalue j+1 is column j - i*column j+1.) ! ! EXP(S) Exponential of S. ! ! GCD(I,J) Greatest common divisor of two integers. ! ! HOUSE(V) The Householder elementary reflector matrix H: ! H = I - 2*(V * TRANSPOSE(V)) / (TRANS(V)*V) ! ! HUGE_I The largest integer (not represented exactly, because ! it is saved as a real number!). ! ! HUGE_R The largest real number. ! ! ID(I) The matrix identity of order I. ! ID(3) is the 3 by 3 identity, for instance. ! ! INT(*) Truncates real values to their integer part. ! INT(1.0) = INT(1.1) = INT(1.9) = 1. ! ! INV(M) The inverse matrix of the square matrix M. ! If M is singular, there will be no inverse. ! ! IVAL(M) Imaginary parts of eigenvalues of square matrix M ! stored in a column vector. ! ! KRYLOV(M,V,I) Returns a rectangular matrix whose columns are ! a sequence of I Krylov vectors: ! ! V, M*V, M*M*V, ..., M**(I-1)*V. ! ! LCM(I,J) Least common multiple of two integers. ! ! LN(S) The natural logarithm of S. ! S must be greater than 0. ! ! LOG(S) The natural logarithm of S. ! S must be greater than 0. ! ! LOG10(S) The logarithm base 10 of S. ! S must be greater than 0. ! ! LOG2(S) The logarithm base 2 of S. ! S must be greater than 0. ! ! MAX(S1,S2) The maximum of S1 and S2. ! ! MIN(S1,S2) Minimum of S1 and S2. ! ! NEG(*) Changes sign of *. ! ! NINT(*) Nearest integer value to entries of *. ! ! NORM0(*) Maximum or infinity norm of a vector or matrix. ! NORM0(S) = ABS(S) ! NORM0(V) returns the maximum of the absolute ! values of the entries of V. ! NORM0(M) returns the maximum of the sum of the ! absolute values of the entries in each row. ! ! NORM1(*) The L1-norm of a vector or matrix. ! NORM1(S) = ABS(S). ! NORM1(V) returns the sum of the absolute values ! of the entries of V. ! NORM0(M) returns the maximum of the sum of the ! absolute values of the entries in each column. ! ! NORM2(MV) L2-norm, Euclidean norm or root-mean-square norm ! of a vector or a square matrix M. NORM2(V) returns ! the square root of the sum of the squares of the ! entries of V. NORM2(M) returns the maximum magnitude ! of the eigenvalues of Transpose(M)*M. ! ! NORMF(*) The Frobenius norm. NORMF(M) returns the square ! root of the sum of the squares of all the entries ! of a matrix M. NORMF may also be applied to a ! vector, giving the same results as NORM2(V). ! ! NORMS(MV) The spectral "norm". NORMS(M) returns the ! value of the maximum of the absolute values of ! all the eigenvalues of the square matrix M. ! ! PI 3.14159265358979... ! You may not change the value of PI. ! ! POLY(V,M) Polynomial evaluation. V contains the ! coefficients of the polynomial, with V(1) the ! coefficient of M**(N-1) and V(N) the constant term. ! M is the scalar or square matrix argument of the polynomial. ! ! RAN(*) Fills * with random numbers between 0 and 1. ! ! RCOND(M) RCOND(M) is the estimate of the reciprocal ! condition number of the matrix M, as computed ! by the LINPACK routine SGE_CO. ! ! The condition number of a matrix M is defined as ! ! COND(M) = NORM(M) * NORM(INVERSE(M)) ! ! RCOND assumes that the L-1 norm is used to define ! COND(M), and then RCOND(M) is defined by ! ! RCOND(M) = 1 / COND(M). ! ! If the matrix M is singular, then the condition ! number is infinite, and RCOND will be zero. ! ! ROUND(*) Replaces all "small" entries of * by 0. ! Here "small" means less than RTOL in absolute value. ! ! RTOL The tolerance for the ROUND command, initially set ! to EPS. To change the value of RTOL, simply ! reassign it with a statement like "RTOL = 0.1E-6". ! ! RVAL(M) Real parts of eigenvalues of square matrix M, ! stored in a column vector. ! ! SEC(S) The secant of angle S. ! ! SIN(S) The sine of angle S. ! ! SINH(S) The hyperbolic sine of angle S. ! ! SQRT(S) Square root of S. S must be non-negative. ! ! STEP(S) Heavyside step function. ! ! STEP = 0 if S < 0, ! STEP = 1 if 0 <= S. ! ! TAN(S) The tangent of angle S. ! ! TANH(S) The hyperbolic tangent of angle S. ! ! TRACE(*) Sum of diagonal elements of *. ! ! TRANS(*) Transpose of matrix or vector. ! ! A=TRANS(A) is a legal formula if A is square. ! ! ZEROS(I) Matrix zero of order I. ! ZEROS(4) is the 4 by 4 zero matrix. ! ! I! Factorial of I. I should be an integer from 0 to ! 25. 0! = 1, 1!=1, 2!=1*2 = 2, 3!=1*2*3=6, and so on. ! For large I, the exact value is not returned. ! ! Extensions: ! ! To add a new function to COMPRN, update the information in the ! IOPSYM, IPRSYM, and SYMBOL arrays, and then insert text into FUNVAL ! or FUNSCL which evaluates the function, and determines the dimensions ! of its result. ! ! Modified: ! ! 10 January 2004 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character COM, the command to be executed, which determines ! what other input must be supplied. ! ! COM='A' Add (=declare) a variable: ! ! Causes the program to add the variable named in NAMVAR ! to its list of symbols. From now on, it may be used in ! formulas, and values may be assigned to it. If you ! declare a name which already exists, then if that name ! was already declared by the program, your command is ! ignored. Otherwise, your new definition overrides your old one. ! Input for the "A" command includes NAMVAR, NROW, and NCOL. ! ! COM='D' delete a variable: ! ! This command can be used to free up memory. You ! must also give the name of the variable to be deleted. ! The variable must be one you have already defined. ! Input for the "D" command includes NAMVAR. ! ! COM='E' evaluate a formula: ! ! Using current values of variables, evaluate the formula. ! The formula number is given by IFRM. The value of the ! formula is returned in VALUE, and its dimensions in NROW and NCOL. ! Input to the "E" command includes IFRM. ! Output from the "E" command includes VALUE. ! ! COM='F', Formula is to be compiled: ! ! The user assigns an index, IFRM, to the formula, so that ! it can be referred to later. The formula is stored in ! INFIX on input. The formula may refer to variables which ! have not yet been declared. ! Output from the "F" command includes IFRM. ! ! COM='G', Formula is to be compiled, variables already declared: ! ! The same as 'F' except that the formula may not refer to ! any variables which have not been declared. ! Output from the "G" command includes IFRM. ! ! COM='I' Initialize: ! ! This must be the first call, to initialize COMRPN's internal ! data. Also, if old data is to be flushed out, this command ! can be used. ! ! COM='L' List symbols and values: ! ! This command prints out the internal symbols, variables ! and values, and the formulas you have supplied. ! The value of NAMVAR determines what is printed. ! If NAMVAR=' ', all user variables are printed. ! If NAMVAR='ALL', internal variables are also printed. ! If NAMVAR='DEBUG', all that and more is printed. ! Otherwise, NAMVAR is assumed to be the name of a ! particular variable and only that is printed. ! If NAMVAR='#', the contents of VALUE are listed. ! ! COM='R' Read value of symbol: ! ! The program returns the current value of the variable ! NAMVAR. VALUE will contain the value, and NROW and NCOL ! the dimensions. ! ! COM='V' assign value to symbol: ! ! This command assigns the value of the numbers in VALUE ! to the variable NAMVAR, with NROW and NCOL specifying ! the dimension of the variable. ! ! Output, integer IERROR, an error flag. ! 0 means no error occurred. ! nonzero means some kind of error. These errors are usually ! signaled by a printed message. ! ! Input/output, integer IFRM. ! ! COMRPN can store more than one formula internally, at ! one time. COMRPN refers to the formulas by an index ! number. When a formula is entered, with the "F" or "G" ! command, it is assigned an index number, whose value is ! returned to the user. Then, whenever the formula is to ! be evaluated, the user must input the value of IFRM, so ! that COMRPN knows which formula to evaluate. In many ! cases, only one formula is being used, so that IFRM is ! usually just 1. ! For the "E" command, IFRM is the index of the formula to evaluate. ! For the "F" or "G" commands, the index assigned to the formula. ! ! Input, character ( len = MAXFIX ) INFIX, for the F or G commands, ! contains the user's formula to be compiled and evaluated. ! ! Workspace, integer IRPN(MAXRPN), for the F, G or E commands, ! a vector used to store the compiled versions of the input formulas. ! ! Input, integer MAXROW, the maximum number of rows and columns. ! ! Input, integer MAXRPN, the maximum length of IRPN. 80 should be enough. ! ! Input, character ( len = * ) NAMVAR, used for the A, D, L, R and ! V commands, containing the name of the variable to be ! added, deleted, listed, read or valued. ! NAMVAR should not be a blank, nor should it be the name 'VOID'. ! For the L command only, using NAMVAR='DEBUG' produces extra output. ! ! Input/output, integer NCOL, the number of columns in the variable. ! ! Input/output, integer NROW, the number of rows in the variable. ! ! Input/output, real VALUE(MAXROW,MAXROW), used to pass values in and out ! of the program. The actual size of the object in VALUE ! is (NROW,NCOL). Thus to pass a scalar, set VALUE(1,1) ! to the value and set NROW=NCOL = 1. To pass a 2 by 3 ! matrix, set the values of VALUE and set NROW = 2, NCOL=3. ! For the V command, VALUE is the input value to be assigned to a variable. ! For the E command, VALUE is the output value of the formula. ! For the R command, VALUE is the output value of the variable. ! ! Internal variables: ! ! Internal, integer IOPSYM(MAXSYM), the number of operands per symbol. ! Constants have 0, unaries have 1, and so on. ! ! Internal, integer IPRSYM(MAXSYM), the priority of each symbol. ! ! Internal, integer IPVAL(MAXSYM), contains, for each ! symbol, the address in VALSYM where the value of the ! symbol is stored, if it is a scalar. If the symbol ! represents a vector or matrix, IPVAL(IARG) points to ! the location of the first entry. ! ! Internal, integer MAXIW, controls the size of the integer ! work vector IWORK, which is needed by EIGEN, SGE_CO, and SGE_FA. ! ! Internal, integer MAXSYM, controls the maximum number ! of symbols allowed. This includes permanent symbols ! defined by the program, as well as symbols the user ! declares while running the program. ! implicit none integer, parameter :: maxsym = 150 integer, parameter :: maxval = 3000 integer, parameter :: maxlen = 20 integer maxrpn character com integer ierror integer, save :: ifinis integer, save :: ifree integer ifreesv integer ifrm integer ihi integer ilo integer imin integer implic integer, save :: indx1 integer, save :: indx2 integer, save :: ineg character ( len = * ) infix integer, save, dimension ( 80 ) :: intsym integer, save, dimension ( maxsym ) :: iopsym integer, save, dimension ( maxsym ) :: iprsym integer, save, dimension ( maxsym ) :: ipval integer, save :: irad integer irpn(maxrpn) integer, save :: irtol integer, save, dimension ( maxsym ) :: istack integer maxfix integer maxfrm integer maxrow character ( len = * ) namvar character ( len = maxlen ) namvr integer ncol integer, save :: nints integer nrow integer nrpn integer, save :: nsym integer, save :: nsymp integer, save :: nsyms integer, save, dimension ( 2, maxsym ) :: numdim integer nuvar logical s_eqi logical s_paren_check character ( len = maxlen ), save, dimension ( maxsym ) :: symbol real value(maxrow,maxrow) real, save, dimension ( maxval ) :: valsym ierror = 0 implic = 0 namvr = namvar call s_blank_delete ( namvr ) maxfix = len ( infix ) if ( 0 < maxfix ) then maxfrm = ( maxrpn / maxfix ) else maxfrm = 0 end if if ( s_eqi ( com, 'E' ) .or. s_eqi ( com, 'G' ) ) then if ( ifrm <= 0 .or. maxfrm < ifrm ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'COMRPN - Error!' write ( *, '(a,i8)' ) ' Illegal formula index = ', ifrm ierror = 1 return end if end if ! ! COM=I initialize. ! if ( s_eqi ( com, 'I' ) ) then call inicom ( ifinis, ifree, indx1, indx2, ineg, infix, intsym, iopsym, & iprsym, ipval, irad, irpn, irtol, istack, maxrow, maxrpn, & maxsym, maxval, nints, nsym, nsymp, nsyms, numdim, symbol, valsym, value ) return ! ! COM=L, List the value of symbol NAMVR. ! else if ( s_eqi ( com, 'L' ) ) then if ( namvr == '#' ) then namvr = ' ' call r4mat_print ( maxrow, nrow, ncol, value, namvr ) else call symbol_list ( ifree, intsym, iopsym, iprsym, ipval, irpn, & maxfrm, maxrpn, maxsym, maxval, namvr, nints, nsym, nsymp, nsyms, & numdim, symbol, valsym ) end if return ! ! COM=A, Add variable. ! else if ( s_eqi ( com, 'A' ) ) then call symbol_add ( ierror, ifree, iopsym, iprsym, ipval, maxrow, & maxsym, maxval, namvr, ncol, nrow, nsym, nsymp, numdim, symbol, valsym ) if ( ierror == 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'COMRPN:' write ( *, '(a)' ) ' Added the variable ' // trim ( namvr ) end if return ! ! COM=D, Delete variable. ! COM=V, Set variable value. ! COM=R, Get variable value. ! else if ( s_eqi ( com, 'D' ) .or. s_eqi ( com, 'R' ) .or. & s_eqi ( com, 'V' ) ) then call symbol_value ( com, ierror, ifree, iopsym, iprsym, ipval, maxrow, & maxsym, maxval, namvr, ncol, nrow, nsym, nsymp, numdim, symbol, valsym, & value ) return ! ! COM=F or G, compile formula. ! else if ( s_eqi ( com, 'F' ) .or. s_eqi ( com, 'G' ) ) then if ( s_eqi ( com, 'F' ) ) then nuvar = 1 else if ( s_eqi ( com, 'G' ) ) then nuvar = 0 end if ! ! Remove all blanks from the formula. ! call s_blank_delete ( infix ) ! ! Reject the formula if it has no information in it. ! if ( len_trim ( infix ) <= 0 ) then ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'COMRPN - Error!' write ( *, '(a)' ) ' The formula has zero length.' return end if ! ! Do a simple check on parentheses. ! if ( .not. s_paren_check ( infix ) ) then ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'COMRPN - Error!' write ( *, '(a)' ) ' The formula did not pass the parentheses checks.' return end if ! ! Convert INFIX, the string of characters, into INTSYM, an array of ! integers which are indices of tokens. ! call tokens ( ierror, ifinis, ifree, implic, indx1, indx2, ineg, infix, & intsym, iopsym, iprsym, ipval, maxrow, maxsym, maxval, & nints, nsym, nsymp, numdim, nuvar, symbol, valsym ) if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'COMRPN - Error!' write ( *, '(a)' ) ' Could not convert formula into tokens.' if ( implic /= 0 ) then symbol(implic) = 'VOID' numdim(1,implic) = 1 numdim(2,implic) = 1 if ( implic == nsym ) then nsym = nsym - 1 end if end if return end if ! ! Convert INTSYM, which is an infix formula, into IRPN, which ! is an RPN formula. ! imin = ( ifrm - 1 ) * 80 call rpnset ( ierror, ifinis, imin, intsym, iopsym, iprsym, irpn, istack, & maxrpn, maxsym, nints, nrpn, symbol ) if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'COMRPN - Error!' write ( *, '(a)' ) ' Could not compile formula.' if ( implic /= 0 ) then symbol(implic) = 'VOID' numdim(1,implic) = 1 numdim(2,implic) = 1 if ( implic == nsym ) then nsym = nsym - 1 end if end if return end if ! ! Check that operators and arguments are balanced. ! ihi = imin + nrpn call rpnchk ( ierror, ihi, ilo, iopsym, irpn, maxrpn, maxsym ) if ( ierror /= 0 .or. ilo /= imin + 1 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'COMRPN - Error!' write ( *, '(a)' ) ' Parentheses mismatch.' ierror = 1 if ( implic /= 0 ) then symbol(implic) = 'VOID' numdim(1,implic) = 1 numdim(2,implic) = 1 if ( implic == nsym ) then nsym = nsym - 1 end if end if return end if ! ! For implicit definition via equality, evaluate formula to ! get dimensions. ! if ( implic == 0 ) then return end if numdim(1,implic) = 0 numdim(2,implic) = 0 else if ( s_eqi ( com, 'E' ) ) then else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'COMRPN - Error!' write ( *, '(a)' ) ' Unknown command = "' // trim ( com ) // '"' ierror = 1 return end if ! ! COM=E, Evaluate compiled formula IFRM. ! (You can also get here when setting up an implicitly defined variable) ! imin = ( ifrm - 1 ) * 80 nrpn = 80 ifreesv = ifree call rpnval ( ierror, ifree, imin, iopsym, iprsym, ipval, irad, irpn, & irtol, istack, maxrow, maxrpn, maxsym, maxval, ncol, nrow, nrpn, nsym, & nsyms, numdim, symbol, valsym, value ) ifree = ifreesv if ( ierror /= 0 ) then nrow = 1 ncol = 1 value(1,1) = 0.0E+00 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'COMRPN - Error!' write ( *, '(a)' ) ' The formula could not be evaluated!' if ( implic /= 0 ) then symbol(implic) = 'VOID' numdim(1,implic) = 1 numdim(2,implic) = 1 if ( implic == nsym ) then nsym = nsym - 1 end if end if end if return end subroutine digit_to_ch ( digit, c ) !*****************************************************************************80 ! !! DIGIT_TO_CH returns the character representation of a decimal digit. ! ! Example: ! ! DIGIT C ! ----- --- ! 0 '0' ! 1 '1' ! ... ... ! 9 '9' ! 17 '*' ! ! Modified: ! ! 04 August 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer DIGIT, the digit value between 0 and 9. ! ! Output, character C, the corresponding character, or '*' if DIGIT ! was illegal. ! implicit none character c integer digit if ( 0 <= digit .and. digit <= 9 ) then c = char ( digit + 48 ) else c = '*' end if return end subroutine eigen ( a, evali, evalr, evvect, ierror, matz, nrow ) !*****************************************************************************80 ! !! EIGEN computes the eigenvalues and eigenvectors for a matrix A. ! ! Discussion: ! ! It does this by calling the EISPACK routine RG. If eigenvectors ! are computed, EIGEN normalizes those eigenvectors to have ! Euclidean norm 1. Note that this is slightly complicated in ! the case where two columns of the eigenvector array are used ! to store the real and imaginary parts of a complex eigenvector. ! ! Modified: ! ! 15 March 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, real A(NROW,NROW). ! ! On input, the matrix to be analyzed. ! ! On output, A has been "destroyed", that is, it has been ! overwritten by information needed for the computation of ! the eigenvalues. ! ! Output, real EVALI(NROW), will contain the imaginary ! parts of the eigenvalues of A. ! ! Output, real EVALR(NROW), will contain the real parts of ! the eigenvalues of A. ! ! Output, real EVVECT(NROW,NROW), will contain the matrix ! of eigenvectors. But EVVECT is only computed if MATZ is nonzero. ! ! Output, integer IERROR, error flag. ! 0, no error. ! nonzero, an error occurred. Usually, this error ! indicates that the eigenvalue iteration did not converge. ! ! Input, integer MATZ. ! 0, Only eigenvalues are wanted. ! nonzero, Eigenvalues and eigenvectors are wanted. ! ! Input, integer NROW, the order of the matrix A. ! implicit none integer nrow real a(nrow, nrow) real colnrm real evali(nrow) real evalr(nrow) real evvect(nrow, nrow) integer ierror integer j integer matz ! ! Call the EISPACK routine RG to produce the eigenvalues and ! eigenvectors. ! call rg ( nrow, nrow, a, evalr, evali, matz, evvect, ierror ) ! ! Normalize each column of the eigenvector array so that the ! two-norm is 1. ! ! In cases where a complex conjugate pair has been computed, ! the eigenvalue with positive imaginary part comes first. ! This allows us to properly normalize both columns. ! if ( matz /= 0 ) then do j = 1, nrow if ( evali(j) == 0.0E+00 ) then colnrm = sqrt ( sum ( evvect(1:nrow,j)**2 ) ) evvect(1:nrow,j) = evvect(1:nrow,j) / colnrm else if ( 0.0E+00 < evali(j) ) then colnrm = sqrt ( sum ( evvect(1:nrow,j)**2 + evvect(1:nrow,j+1)**2) ) evvect(1:nrow,j) = evvect(1:nrow,j) / colnrm evvect(1:nrow,j+1) = evvect(1:nrow,j+1) / colnrm else if ( evali(j) < 0.0E+00 ) then end if end do end if return end subroutine elmhes ( lda, n, low, igh, a, jint ) !*****************************************************************************80 ! !! ELMHES reduces all, or a portion of a matrix, to upper Hessenberg form. ! ! Discussion: ! ! The routine uses stabilized elementary similarity transformations. ! ! Modified: ! ! 15 March 2001 ! ! Parameters: ! ! Input, integer LDA, the leading dimension of A. ! ! Input, integer N, the number of rows and columns in A. ! ! Input, integer LOW, IGH. If subroutine BALANC was called ! before ELMHES, then it will set LOW and IGH. If BALANC ! was not called, then the user should set LOW = 1 and IGH=N. ! ! Input/output, real A(LDA,N). ! ! On input, A contains the matrix to be transformed. ! ! On output, A contains the Hessenberg matrix. The multipliers ! which were used in the reduction are stored in the ! remaining triangle under the Hessenberg matrix. ! ! Output, integer JINT(IGH), contains information on the rows ! and columns interchanged in the reduction. ! Only elements LOW through IGH are used. ! implicit none integer igh integer lda integer n real a(lda,n) integer i integer j integer jint(igh) integer low integer m real x real y do m = low+1, igh-1 ! ! Look for the largest element in the column A(J,M-1), where ! J goes from M to IGH. Store the row number as I. ! x = 0.0E+00 i = m do j = m, igh if ( abs ( x ) < abs ( a(j,m-1) ) ) then x = a(j,m-1) i = j end if end do jint(m) = i ! ! If I is not M, interchange rows and columns I and M of A. ! if ( i /= m ) then do j = m-1, n call r4_swap ( a(i,j), a(m,j) ) end do do j = 1, igh call r4_swap ( a(j,i), a(j,m) ) end do end if if ( x /= 0.0E+00 ) then do i = m+1, igh y = a(i,m-1) if ( y /= 0.0E+00 ) then y = y / x a(i,m-1) = y a(i,m:n) = a(i,m:n) - y * a(m,m:n) do j = 1, igh a(j,m) = a(j,m) + y * a(j,i) end do end if end do end if end do return end subroutine eltran ( lda, n, low, igh, a, jint, z ) !*****************************************************************************80 ! !! ELTRAN accumulates transformations used by ELMHES. ! ! Modified: ! ! 15 March 2001 ! ! Parameters: ! ! Input, integer LDA, the leading dimension of A and Z. ! ! Input, integer N, the number of rows and columns in A. ! ! Input, integer LOW, IGH. If BALANC was called ! before ELMHES, then it will set LOW and IGH. If BALANC ! was not called, then the user should set LOW = 1 and IGH=N. ! ! Input, integer A(LDA,N), contains, in its lower triangle, ! the multipliers used by ELMHES in the reduction. ! ! Input, integer JINT(IGH), contains information on the rows ! and columns interchanged in the reduction. ! ! Output, real Z(LDA,N), contains the transformation matrix ! produced in the reduction by ELMHES. ! implicit none integer igh integer lda integer n real a(lda,igh) integer i integer j integer jint(igh) integer low integer mm integer mp real z(lda,n) ! ! Initialize Z to the identity matrix. ! call r4mat_identity ( lda, n, z ) do mm = 1, igh-low-1 mp = igh - mm do i = mp+1, igh z(i,mp) = a(i,mp-1) end do i = jint(mp) if ( i /= mp ) then do j = mp, igh z(mp,j) = z(i,j) z(i,j) = 0.0E+00 end do z(i,mp) = 1.0E+00 end if end do return end subroutine enter ( ierror, line, maxrow, maxrpn, namvar ) !*****************************************************************************80 ! !! ENTER allows a user to enter a new variable, and assign it a value. ! ! Modified: ! ! 01 May 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer IERROR, error flag. ! 0, no error occurred. ! nonzero, an error occurred. ! ! Input/output, character ( len = 80 ) LINE, a buffer containing ! unprocessed user input. ! ! Input, integer MAXROW, the maximum number of rows and columns. ! ! Input, integer MAXRPN, specifies the length of IRPN. ! ! Output, character ( len = MAXLEN ) NAMVAR, the name of the variable that ! has been entered. ! implicit none integer, parameter :: maxfix = 80 integer, parameter :: maxlen = 20 integer maxrow integer maxrpn character com integer i integer ierror integer ifrm character ( len = maxfix ) infix integer irpn(maxrpn) integer j character ( len = 80 ) line character ( len = 15 ) name1 character ( len = 15 ) name2 character ( len = maxlen ) namvar integer ncol integer nhi integer nrow character ( len = 80 ) prompt real rval real value(maxrow,maxrow) ! ! Get NAMVAR, the name of the variable. ! prompt = 'variable, number of rows, number of columns' call chrrea ( namvar, line, prompt, ierror, 3 ) if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'ENTER - Error!' write ( *, '(a)' ) ' Could not read the variable name.' return end if ! ! Get NROW, the number of rows. ! call intrea ( ierror, nrow, line, prompt ) if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'ENTER - Error!' write ( *, '(a)' ) ' Could not read the number of rows.' return end if ! ! Get NCOL, the number of columns. ! call intrea ( ierror, ncol, line, prompt ) if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'ENTER - Error!' write ( *, '(a)' ) ' Could not read the number of columns.' return end if ! ! Add NAMVAR to the list of variables. ! ifrm = 0 com = 'A' call comrpn ( com, ierror, ifrm, infix, irpn, maxrow, maxrpn, namvar, & ncol, nrow, value ) if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'ENTER - Warning!' write ( *, '(a)' ) ' Could not add variable name to symbol list.' return end if ! ! Get a scalar value. ! if ( nrow == 1 .and. ncol == 1 ) then prompt = namvar call relrea ( rval, line, prompt, ierror ) if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'ENTER - Error!' write ( *, '(a)' ) ' Could not read the variable value.' return end if value(1,1) = rval ! ! Get a vector value. ! else if ( ncol == 1 .or. nrow == 1 ) then if ( ncol == 1 ) then nhi = nrow else nhi = ncol end if line = ' ' do i = 1, nhi if ( i < nhi ) then write ( prompt, '(''entries '', i3, '' through '', i3)' ) i, nhi else write ( prompt, '(''entry '', i3 )' ) i end if call s_blanks_delete ( prompt ) call relrea ( rval, line, prompt, ierror ) if ( ierror /= 0 ) then return end if if ( ncol == 1 ) then value(i,1) = rval else value(1,i) = rval end if end do ! ! Get a matrix value. ! else line = ' ' do i = 1, nrow do j = 1, ncol write ( name1, '( ''('', i6, '','', i6, '')'' )' ) i, j call s_blank_delete ( name1 ) write ( name2, '( ''('', i6, '','', i6, '')'' )' ) i, ncol call s_blank_delete ( name2 ) if ( j < ncol ) then prompt = 'entries ' // trim ( name1 ) // ' through ' // trim ( name2 ) else prompt = 'entry ' // trim ( name1 ) end if call s_blanks_delete ( prompt ) call relrea ( rval, line, prompt, ierror ) if ( ierror /= 0 ) then return end if value(i,j) = rval end do end do end if ! ! Store the value of the variable. ! com = 'V' call comrpn ( com, ierror, ifrm, infix, irpn, maxrow, maxrpn, namvar, & ncol, nrow, value ) if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'ENTER - Error!' write ( *, '(a)' ) ' Could not store the variable value.' return end if return end subroutine factor ( ierror, line, maxrow, maxrpn ) !*****************************************************************************80 ! !! FACTOR computes the PLU, QR or USV factorization of a matrix. ! ! Discussion: ! ! A=P*L*U ! P is an M by M permuation matrix, ! L is an M by M lower triangular matrix, ! U is an M by N upper triangular matrix. ! ! A=Q*R ! Q orthogonal, R upper triangular ! A is M by N, and ! Q is M by M, R is M by N. ! ! A=U*S*TRANS(V) ! U and V orthogonal, S diagonal ! A is M by N, and ! U is M by N, S is N by N, and V is N by N. ! ! Modified: ! ! 01 May 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer IERROR, error flag. ! 0, no error. ! nonzero, an error occurred. ! ! Input/output, character ( len = 80 ) LINE, a buffer containing ! unprocessed user input. ! ! Input, integer MAXFIX, specifies the maximum length of a ! formula that a user can enter, in characters. ! ! Input, integer MAXROW, the maximum number of rows and columns. ! ! Input, integer MAXRPN, specifies the length of IRPN. ! implicit none integer, parameter :: maxfix = 80 integer, parameter :: maxlen = 20 integer maxrow integer maxrpn character com integer i integer ierror integer ifrm character ( len = maxfix ) infix integer irpn(maxrpn) character isay character ( len = 3 ) itype integer lcol logical s_eqi character ( len = 80 ) line integer lrow character ( len = maxlen ) namvar integer ncol integer nrow character ( len = 80 ) prompt real pvec(maxrow) real value(maxrow,maxrow) real value1(maxrow,maxrow) real value2(maxrow,maxrow) real value3(maxrow,maxrow) prompt = 'PLU or QR or SVD for factorization desired.' call chrrea ( itype, line, prompt, ierror, 0 ) if ( ierror /= 0 ) then return end if if ( s_eqi ( itype(1:1), 'Q' ) ) then itype = 'QR' else if ( s_eqi ( itype(1:1), 'P' ) ) then itype = 'PLU' else if ( s_eqi ( itype(1:1), 'S' ) ) then itype = 'SVD' else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FACTOR - Error!' write ( *, '(a)' ) ' Your choice was not acceptable!' ierror = 1 return end if prompt = 'matrix to factor.' call chrrea ( namvar, line, prompt, ierror, 3 ) if ( ierror /= 0 ) then return end if com = 'R' ifrm = 1 ncol = 1 nrow = 1 call comrpn ( com, ierror, ifrm, infix, irpn, maxrow, maxrpn, namvar, & ncol, nrow, value ) if ( ierror /= 0 ) then return end if if ( s_eqi ( itype, 'PLU' ) ) then call sge_plu ( maxrow, nrow, ncol, value, value3, value1, value2 ) else if ( s_eqi ( itype, 'QR' ) ) then call qrfact ( value, nrow, maxrow, ncol, value1, value2 ) else call svd ( maxrow, nrow, ncol, value, pvec, .true., value1, .true., & value2, ierror ) value3(1:nrow,1:ncol) = 0.0E+00 do i = 1, nrow value3(i,i) = pvec(i) end do end if ! ! Print the factors. ! if ( s_eqi ( itype, 'QR' ) ) then namvar = 'Q Factor' lrow = nrow lcol = nrow else if ( s_eqi ( itype, 'PLU' ) ) then namvar = 'L Factor' lrow = nrow lcol = nrow else namvar = 'U Factor' lrow = nrow lcol = ncol end if call r4mat_print ( maxrow, lrow, lcol, value1, namvar ) if ( s_eqi ( itype, 'SVD' ) ) then namvar = 'Singular Values' call r4mat_print ( maxrow, ncol, ncol, value3, namvar ) end if if ( s_eqi ( itype, 'QR' ) ) then namvar = 'R Factor' lrow = nrow lcol = ncol else if ( s_eqi ( itype, 'PLU' ) ) then namvar = 'U Factor' lrow = nrow lcol = ncol else namvar = 'V Factor' lrow = ncol lcol = ncol end if call r4mat_print ( maxrow, lrow, lcol, value2, namvar ) if ( s_eqi ( itype, 'PLU' ) ) then namvar = 'P Factor' call r4mat_print ( maxrow, nrow, nrow, value3, namvar ) end if ! ! Save factors? ! line = ' ' if ( s_eqi ( itype, 'QR' ) ) then isay = 'Q' lrow = nrow lcol = nrow else if ( s_eqi ( itype, 'PLU' ) ) then isay = 'L' lrow = nrow lcol = nrow else isay = 'U' lrow = nrow lcol = ncol end if prompt = 'matrix to save ' // isay // ' in, or return' namvar = ' ' call chrrea ( namvar, line, prompt, ierror, 0 ) if ( ierror /= 0 ) then return end if if ( namvar /= ' ' ) then com = 'A' call comrpn ( com, ierror, ifrm, infix, irpn, maxrow, maxrpn, & namvar, lcol, lrow, value1 ) if ( ierror /= 0 ) then return end if com = 'V' call comrpn ( com, ierror, ifrm, infix, irpn, maxrow, maxrpn, & namvar, lcol, lrow, value1 ) if ( ierror /= 0 ) then return end if end if if ( s_eqi ( itype, 'SVD' ) ) then line = ' ' prompt = 'matrix to save S in, or return' namvar = ' ' call chrrea ( namvar, line, prompt, ierror, 0 ) if ( ierror /= 0 ) then return end if if ( namvar /= ' ' ) then com = 'A' call comrpn ( com, ierror, ifrm, infix, irpn, maxrow, maxrpn, & namvar, ncol, ncol, value3 ) if ( ierror /= 0 ) then return end if com = 'V' call comrpn ( com, ierror, ifrm, infix, irpn, maxrow, maxrpn, & namvar, ncol, ncol, value3 ) if ( ierror /= 0 ) then return end if end if end if if ( s_eqi ( itype, 'QR' ) ) then isay = 'R' lrow = nrow lcol = ncol else if ( s_eqi ( itype, 'PLU' ) ) then isay = 'U' lrow = nrow lcol = ncol else isay = 'V' lrow = ncol lcol = ncol end if line = ' ' prompt = 'matrix to save ' // isay // ' in, or return' namvar = ' ' call chrrea ( namvar, line, prompt, ierror, 0 ) if ( ierror /= 0 ) then return end if if ( namvar /= ' ' ) then com = 'A' call comrpn ( com, ierror, ifrm, infix, irpn, maxrow, maxrpn, namvar, & lcol, lrow, value2 ) if ( ierror /= 0 ) then return end if com = 'V' call comrpn ( com, ierror, ifrm, infix, irpn, maxrow, maxrpn, namvar, & lcol, lrow, value2 ) if ( ierror /= 0 ) then return end if end if if ( s_eqi ( itype, 'PLU' ) ) then isay = 'P' prompt = 'matrix to save ' // isay // ' in or return' namvar = ' ' call chrrea ( namvar, line, prompt, ierror, 0 ) if ( ierror /= 0 ) then return end if if ( namvar /= ' ' ) then com = 'A' call comrpn ( com, ierror, ifrm, infix, irpn, maxrow, maxrpn, & namvar, nrow, nrow, value3 ) if ( ierror /= 0 ) then return end if com = 'V' call comrpn ( com, ierror, ifrm, infix, irpn, maxrow, maxrpn, & namvar, nrow, nrow, value3 ) if ( ierror /= 0 ) then return end if end if end if return end subroutine funscl ( angle_to_radian, arg1, arg2, ierror, result, rtol, sym ) ! !*****************************************************************************80 ! !! FUNSCL evaluates a scalar function of one or more scalar arguments. ! ! ! Modified: ! ! 15 April 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ARG2, ARG2, the values of the arguments of the function. ! ! Output, integer IERROR, error flag. ! 0, no error occurred. ! nonzero, an error occurred. ! ! Output, real RESULT, the value of the scalar function. ! ! Input, real RTOL, the rounding tolerance. ! ! Input, character ( len = * ) SYM, the symbolic name of the ! function to be evaluated. ! implicit none real, parameter :: degrad = 3.14159265358979 / 180.0E+00 real angle_to_radian real arg1 real arg2 integer i4_lcm integer iarg integer iarg1 integer iarg2 integer ierror integer i4_gcd real result real rtol logical s_eqi real sarg character ( len = * ) sym real temp ierror = 0 result = 0.0E+00 ! ! Addition. ! if ( sym == '+' ) then result = arg1 + arg2 ! ! Subtraction. ! else if ( sym == '-' ) then result = arg1 - arg2 ! ! Division. ! else if ( sym == '/' ) then if ( arg2 == 0.0E+00 ) then ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FUNSCL - Error!' write ( *, '(a)' ) ' Attempt to divide by 0!' return end if result = arg1 / arg2 ! ! Multiplication. ! else if ( sym == '*' ) then result = arg1 * arg2 ! ! Exponentiation. ! else if ( sym == '**' .or. sym == '^' ) then if ( arg1 == 0.0E+00 .and. arg2 == 0.0E+00 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FUNSCL - Error!' write ( *, '(a)' ) ' Attempt to compute 0**0 !' ierror = 1 else if ( arg1 < 0.0E+00 .and. real(int(arg2)) /= arg2 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FUNSCL - Error!' write ( *, '(a)' ) ' Illegal exponentiation:' write ( *, '(a,g14.6)' ) ' Base = ', arg1 write ( *, '(a,g14.6)' ) ' Power = ', arg2 ierror = 1 else sarg = 1.0E+00 if ( arg1 < 0.0E+00 ) then iarg = int ( arg2 ) if ( 2 * ( iarg / 2 ) /= iarg ) then sarg = - 1.0E+00 end if end if result = sarg * abs ( arg1 )**arg2 end if ! ! Assignment. ! else if ( sym == '=' ) then result = arg2 ! ! Absolute value. ! else if ( s_eqi ( sym, 'ABS' ) ) then result = abs ( arg1 ) ! ! Arc Cosine. ! else if ( s_eqi ( sym, 'ACOS' ) ) then if ( arg1 < -1.0E+00 .or. 1.0E+00 < arg1 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FUNSCL - Error!' write ( *, '(a,g14.6)' ) ' Illegal inverse cosine of ', arg1 ierror = 1 return end if result = acos ( arg1 ) / angle_to_radian ! ! Arc Sine. ! else if ( s_eqi ( sym, 'ASIN' ) ) then if ( arg1 < - 1.0E+00 .or. 1.0E+00 < arg1 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FUNSCL - Error!' write ( *, '(a,g14.6)' ) ' Illegal inverse sine of ', arg1 ierror = 1 end if result = asin ( arg1 ) / angle_to_radian ! ! Arc Tangent. ! else if ( s_eqi ( sym, 'ATAN' ) ) then result = atan ( arg1 ) / angle_to_radian ! ! Arc Tangent of ratio. ! else if ( s_eqi ( sym, 'ATAN2' ) ) then if ( arg1 == 0.0E+00 .and. arg2 == 0.0E+00 ) then result = 0.0E+00 else result = atan2 ( arg1, arg2 ) / angle_to_radian end if ! ! Cosine. ! else if ( s_eqi ( sym, 'COS' ) ) then result = cos ( angle_to_radian * arg1 ) ! ! Cotangent of angle. ! else if ( s_eqi ( sym, 'COT' ) ) then result = cos ( angle_to_radian * arg1 ) / sin ( angle_to_radian * arg1 ) ! ! Hyperbolic cosine. ! else if ( s_eqi ( sym, 'COSH' ) ) then result = cosh ( angle_to_radian * arg1 ) ! ! Cosecant of angle. ! else if ( s_eqi ( sym, 'CSC' ) ) then result = 1.0E+00 / sin ( angle_to_radian * arg1 ) ! ! Determinant (of a scalar). ! else if ( s_eqi ( sym, 'DET' ) ) then result = arg1 ! ! Diagonal (of a scalar). ! else if ( s_eqi ( sym, 'DIAG' ) ) then result = arg1 ! ! Exponential. ! else if ( s_eqi ( sym, 'EXP' ) ) then result = exp ( arg1 ) ! ! Greatest common factor. ! else if ( s_eqi ( sym, 'GCD' ) ) then iarg1 = nint ( arg1 ) iarg2 = nint ( arg2 ) result = real ( i4_gcd ( iarg1, iarg2 ) ) ! ! Identity. ! else if ( s_eqi ( sym, 'ID' ) ) then result = 1.0E+00 ! ! Integer value. ! else if ( s_eqi ( sym, 'INT' ) ) then result = real ( int ( arg1 ) ) ! ! Inverse (of a scalar). ! else if ( s_eqi ( sym, 'INV' ) ) then if ( arg1 == 0.0E+00 ) then ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FUNSCL - Error!' write ( *, '(a)' ) ' Attempt to compute 1/0!' return end if result = 1.0E+00 / arg1 ! ! Least common multiple. ! else if ( s_eqi ( sym, 'LCM' ) ) then iarg1 = nint ( arg1 ) iarg2 = nint ( arg2 ) result = real ( i4_lcm ( iarg1, iarg2 ) ) ! ! Natural logarithm. ! else if ( s_eqi ( sym, 'ALOG' ) .or. s_eqi ( sym, 'LN' ) .or. & s_eqi ( sym, 'LOG' ) ) then if ( arg1 <= 0.0E+00 ) then ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FUNSCL - Error!' write ( *, '(a,g14.6)' ) ' Illegal LOG of ', arg1 return end if result = log ( arg1 ) ! ! Logarithm base 10. ! else if ( s_eqi ( sym, 'ALOG10' ) .or. s_eqi ( sym, 'LOG10' ) ) then if ( arg1 <= 0.0E+00 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FUNSCL - Error!' write ( *, '(a,g14.6)' ) ' Illegal LOG10 of ', arg1 ierror = 1 return end if result = log10 ( arg1 ) ! ! Logarithm base 2. ! else if ( s_eqi ( sym, 'LOG2' ) ) then if ( arg1 <= 0.0E+00 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FUNSCL - Error!' write ( *, '(a,g14.6)' ) ' Illegal LOG2 of ', arg1 ierror = 1 return end if result = log ( arg1 ) / log ( 2.0E+00 ) ! ! Maximum of two values. ! else if ( s_eqi ( sym, 'MAX' ) ) then if ( arg2 < arg1 ) then result = arg1 else result = arg2 end if ! ! Minimum of two values. ! else if ( s_eqi ( sym, 'MIN' ) ) then if ( arg1 < arg2 ) then result = arg1 else result = arg2 end if ! ! Negation ! else if ( s_eqi ( sym, 'NEG' ) ) then result = - arg1 ! ! Nearest integer value ! else if ( s_eqi ( sym, 'NINT' ) ) then result = anint ( arg1 ) ! ! NORM0, NORM1, NORM2, NORME or NORMF of a scalar. ! else if ( s_eqi ( sym, 'NORM0' ) .or. s_eqi ( sym, 'NORM1' ) .or. & s_eqi ( sym, 'NORM2' ) .or. s_eqi ( sym, 'NORME' ) .or. & s_eqi ( sym, 'NORMF' ) ) then result = abs ( arg1 ) ! ! POLY (scalar coefficient array, means constant polynomial). ! else if ( s_eqi ( sym, 'POLY' ) ) then result = arg1 ! ! Random value. ! else if ( s_eqi ( sym, 'RAN' ) ) then call random_number ( result ) ! ! Rounding. ! else if ( s_eqi ( sym, 'ROUND' ) ) then if ( abs ( arg1 ) < rtol ) then result = 0.0E+00 else result = arg1 end if ! ! Secant. ! else if ( s_eqi ( sym, 'SEC' ) ) then result = 1.0E+00 / cos ( angle_to_radian * arg1 ) ! ! Sine. ! else if ( s_eqi ( sym, 'SIN' ) ) then result = sin ( angle_to_radian * arg1 ) ! ! Hyperbolic sine. ! else if ( s_eqi ( sym, 'SINH' ) ) then result = sinh ( angle_to_radian * arg1 ) ! ! Square root. ! else if ( s_eqi ( sym, 'SQRT' ) ) then if ( arg1 < 0.0E+00 ) then ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FUNSCL - Error!' write ( *, '(a,g14.6)' ) ' Illegal SQRT of ', arg1 return end if result = sqrt ( arg1 ) ! ! Step function. ! else if ( s_eqi ( sym, 'STEP' ) ) then if ( arg1 < 0.0E+00 ) then result = 0.0E+00 else result = 1.0E+00 end if ! ! Tangent. ! else if ( s_eqi ( sym, 'TAN' ) ) then result = tan ( angle_to_radian * arg1 ) ! ! Hyperbolic tangent. ! else if ( s_eqi ( sym, 'TANH' ) ) then result = tanh ( angle_to_radian * arg1 ) ! ! Trace (of a scalar) ! else if ( s_eqi ( sym, 'TRACE' ) ) then result = arg1 ! ! Transpose (of a scalar) ! else if ( s_eqi ( sym, 'TRANS' ) ) then result = arg1 ! ! Zero ! else if ( s_eqi ( sym, 'ZEROS' ) ) then result = 0.0E+00 ! ! Factorial function. ! else if ( sym == '!' ) then if ( arg1 < 0.0E+00 .or. 25.0E+00 < arg1 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FUNSCL - Error!' write ( *, '(a)' ) ' Argument out of range for factorial!' ierror = 1 return end if temp = arg1 result = 1.0E+00 do while ( 1.0E+00 < temp ) result = result * temp temp = temp - 1.0E+00 end do ! ! Unknown function. ! else result = 0.0E+00 ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FUNSCL - Error!' write ( *, '(a)' ) ' Unrecognized function name:' // trim ( sym ) end if return end subroutine funval ( iarg1, iarg2, iarg3, ierror, ifree, iopsym, iprsym, & ipset, ipval, irad, irtol, itemp, maxrow, maxsym, maxval, nsyms, numdim, & sym, symbol, valsym, value ) !*****************************************************************************80 ! !! FUNVAL evaluates a scalar, vector or matrix valued function. ! ! Modified: ! ! 17 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer IARG1, IARG2, IARG3, are the indices ! of the arguments to the current function. Actually, ! the function might have just one or two arguments, ! in which case IARG2 and IARG3 may be zero. ! ! Output, integer IERROR. ! ! IERROR is an error flag. If it is zero on return, ! then no error was detected, and the formula was ! evaluated successfully. If it is nonzero on return, ! then an error was found, and the evaluation of the ! formula could not be carried out. ! ! Input, integer IFREE, the next free memory location in VALSYM. ! ! Input, integer IOPSYM(MAXSYM), the number of operands per symbol. ! Constants have 0, unaries have 1, and so on. ! ! Input, integer IPRSYM(MAXSYM), the priority of each symbol. ! ! Output, integer IPSET. If the function was an 'INDX1' or ! 'INDX2', then IPSET is the relative offset of the value. ! ! Input, integer IPVAL(MAXSYM), contains, for each symbol, ! the address in VALSYM where the value of the symbol ! is stored, if it is a scalar. If the symbol represents ! a vector or matrix, IPVAL(IARG) points to location of ! the first entry. ! ! Input, integer ITEMP, the index of the constant. ! ! Input, integer MAXIW, the amount of space in IWORK, which ! should be at least MAXROW. ! ! Input, integer MAXROW, the maximum number of rows and columns. ! ! Input, integer MAXSYM, the maximum number of symbols allowed. ! ! Input, integer MAXVAL, the maximum number of values allowed. ! ! Input, integer NSYMS, the number of declared symbols. ! ! Input, integer NUMDIM(2,MAXSYM). ! For each symbol I, NUMDIM(1,I) is the number of rows in its ! value, and NUMDIM(2,I) is the number of columns. ! ! Input, character ( len = MAXLEN ) SYM, the symbolic name of the function. ! ! Input, character ( len = MAXLEN ) SYMBOL(MAXSYM), the symbolic names ! of the internally defined constants and functions. ! ! Input, real VALSYM(MAXVAL), the values of the symbolic variables. ! ! Workspace, real VALUE(MAXROW,MAXROW), space used to hold ! the value of a variable. ! implicit none integer, parameter :: maxlen = 20 integer maxrow integer maxsym integer maxval real adet real angle_to_radian real arg1 real arg2 real coef character ( len = 3 ) ctemp integer i integer iarg1 integer iarg2 integer iarg3 integer ierror integer ifree integer index1 integer index2 integer index3 integer index4 integer index5 integer index6 integer indx integer info integer iopsym(maxsym) integer iprsym(maxsym) integer ipset integer ipval(maxsym) integer irad integer irtol integer itemp integer iwork(maxrow) integer j integer job integer matz integer na integer nb integer ncoef integer ncol integer ncol1 integer ncol2 integer ndim integer ndim1 integer ndim2 integer npow integer nrow integer nrow1 integer nrow2 integer nsyms integer ntemp integer numdim(2,maxsym) real rcond real result real r4mat_trace real r4mat_norme real r4mat_normf real rtol logical s_eqi character ( len = maxlen ) sym character ( len = maxlen ) symbol(maxsym) real temp real temp1 real tnorm real v1 real v2 real valsym(maxval) real value(maxrow,maxrow) ! ipset = 0 ierror = 0 ! ! If the operator is assignment, then the user may be trying to ! reset the values of certain reserved variables, including ! PI and EPS. ! if ( sym == '=' ) then if ( s_eqi ( symbol(iarg1), 'E' ) .or. & s_eqi ( symbol(iarg1), 'EPS' ) .or. & s_eqi ( symbol(iarg1), 'HUGE_I' ) .or. & s_eqi ( symbol(iarg1), 'HUGE_R' ) .or. & s_eqi ( symbol(iarg1), 'PI' ) ) then ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FUNVAL - Error!' write ( *, '(a)' ) ' You may not change the value of ' // & trim ( symbol(iarg1) ) return end if end if ! ! Set dimensions for ID and ZEROS functions, ! which otherwise would look like scalar functions and get ! trapped in FUNSCL. ! if ( s_eqi ( sym, 'ID' ) .or. & s_eqi ( sym, 'ZEROS' ) ) then nrow1 = int ( valsym(ipval(iarg1)) ) ncol1 = nrow1 else nrow1 = numdim(1,iarg1) ncol1 = numdim(2,iarg1) end if ndim1 = nrow1 * ncol1 index1 = ipval(iarg1) nrow2 = 0 ncol2 = 0 ndim2 = 0 index2 = 0 if ( iarg2 /= 0 ) then nrow2 = numdim(1,iarg2) ncol2 = numdim(2,iarg2) ndim2 = nrow2 * ncol2 index2 = ipval(iarg2) end if index3 = 0 if ( iarg3 /= 0 ) then index3 = ipval(iarg3) end if if ( maxsym <= nsyms ) then ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FUNVAL - Error!' write ( *, '(a)' ) ' Not enough free memory left!' write ( *, '(a)' ) ' The KILL or INIT commands may help!' return end if nsyms = nsyms + 1 iopsym(nsyms) = 0 iprsym(nsyms) = 10 itemp = itemp + 1 ctemp = ' ' call i4_to_s_zero ( itemp, ctemp ) symbol(nsyms) = 'STK000' symbol(nsyms)(4:6) = ctemp ipval(nsyms) = ifree index4 = ipval(nsyms) ! ! Check for simple scalar arithmetic. ! if ( ( ndim1 == 1 .and. iarg2 == 0 ) .or. & ( ndim1 == 1 .and. ndim2 == 1 ) ) then numdim(1,nsyms) = 1 numdim(2,nsyms) = 1 arg1 = valsym(index1) if ( iarg2 == 0 ) then arg2 = 0.0E+00 else arg2 = valsym(ipval(iarg2)) end if angle_to_radian = valsym(irad) rtol = valsym(irtol) call funscl ( angle_to_radian, arg1, arg2, ierror, result, rtol, sym ) valsym(index4) = result if ( sym == '=' ) then arg1 = arg2 valsym(index1) = arg2 end if return end if ! ! Scalar+square matrix or vector+vector or matrix+matrix ! addition or subtraction. ! if ( sym == '+'.or.sym=='-' ) then if ( ndim1 == 0 .and. ndim2 /= 0 ) then ndim1 = ndim2 nrow1 = nrow2 ncol1 = ncol2 numdim(1,iarg1) = numdim(1,iarg2) numdim(2,iarg1) = numdim(2,iarg2) else if (ndim1 /= 0 .and. ndim2 == 0 ) then ndim2 = ndim1 nrow2 = nrow1 ncol2 = ncol1 numdim(1,iarg2) = numdim(1,iarg1) numdim(2,iarg2) = numdim(2,iarg1) end if if ( ( ndim1 /= 1 .or. nrow2 /= ncol2 ) .and. & ( ndim2 /= 1 .or. nrow1 /= ncol1 ) ) then if ( nrow1 /= nrow2 .or. ncol1 /= ncol2 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FUNVAL - Error!' write ( *, '(a)' ) ' The matrices are incompatible for addition' write ( *, '(a)' ) ' or subtraction!' ierror = 1 return end if end if nrow = max ( nrow1, nrow2 ) ncol = max ( ncol1, ncol2 ) numdim(1,nsyms) = nrow numdim(2,nsyms) = ncol do j = 1, ncol do i = 1, nrow if ( ndim1 == 1 .and. i /= j ) then v1 = 0.0E+00 else v1 = valsym(index1) end if if ( ndim2 == 1 .and. i /= j ) then v2 = 0.0E+00 else v2 = valsym(index2) end if if ( sym == '+' ) then valsym(index4) = v1+v2 else valsym(index4) = v1-v2 end if if ( ndim1 /= 1 ) then index1 = index1+1 end if if ( ndim2 /= 1 ) then index2 = index2+1 end if index4 = index4+1 end do end do ! ! Scalar*vector, scalar*matrix, ! Vector*scalar, vector*vector, vector*matrix, ! Matrix*scalar, matrix*vector, matrix*matrix multiplication. ! else if ( sym == '*' ) then if ( ncol1 == nrow2 ) then numdim(1,nsyms) = nrow1 numdim(2,nsyms) = ncol2 do j = 1, ncol2 do i = 1, nrow1 valsym(index4) = dot_product ( & valsym(index1+i-1:index1+i-1+(ncol1-1)*nrow1:nrow1), & valsym(index2+(j-1)*nrow2:index2+(j-1)*nrow2+ncol1-1) ) index4 = index4+1 end do end do else if (ndim1 == 1 ) then numdim(1,nsyms) = numdim(1,iarg2) numdim(2,nsyms) = numdim(2,iarg2) valsym(index4:index4+ndim2-1) = valsym(index2:index2+ndim2-1) do i = index4, index4+ndim2-1 valsym(i) = valsym(i) * valsym(index1) end do else if ( ndim2 == 1 ) then numdim(1,nsyms) = numdim(1,iarg1) numdim(2,nsyms) = numdim(2,iarg1) valsym(index4:index4+ndim1-1) = valsym(index1:index1+ndim1-1) do i = index4, index4+ndim1-1 valsym(i) = valsym(i) * valsym(index2) end do else if ( ( nrow1 == nrow2 .and. ncol1 == 1 .and. ncol2 == 1 ) .or. & ( ncol1 == ncol2 .and. nrow1 == 1 .and. nrow2 == 1 ) ) then valsym(index4) = dot_product ( valsym(index1:index1+ndim1-1), & valsym(index2:index2+ndim2-1) ) numdim(1,nsyms) = 1 numdim(2,nsyms) = 1 else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FUNVAL - Error!' write ( *, '(a)' ) ' The matrices are incompatible for multiplication.' ierror = 1 end if ! ! A \ x ! meaning, solve for y in A*y=x, using PLU factorization of A. ! ! We can't explicitly use the backslash character because ninny ! UNIX FORTRAN compilers interpret it as an escape! ! ! I could eliminate one of the temporary vectors used here. ! else if ( sym == char(92) ) then if ( nrow1 /= ncol1 ) then ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FUNVAL - Error!' write ( *, '(a)' ) ' The matrix must be square for "\"!' return end if if ( ncol1 /= nrow2 ) then ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FUNVAL - Error!' write ( *, '(a)' ) ' Vector/matrix incompatible in "\"' return end if if ( maxval < ifree + nrow1 * ncol1 + 2 * nrow1 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FUNVAL - Error!' write ( *, '(a)' ) ' There is not enough memory for the "\" operation.' ierror = 1 return end if numdim(1,nsyms) = nrow1 numdim(2,nsyms) = ncol2 index5 = index4+nrow1*ncol2 index6 = index5+nrow1*ncol1 ! ! Make a copy of the matrix A. ! valsym(index5:index5+nrow1*ncol1-1) = & valsym(index1:index1+nrow1*ncol1-1) ! ! Factor the copy. ! call sge_fa ( valsym(index5), nrow1, nrow1, iwork, info ) if ( info /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FUNVAL - Error!' write ( *, '(a)' ) ' The matrix in "\" expresssion is singular!' return end if job = 0 ! ! Copy each column of "b" into the temporary vector, solve ! A*x=b, and copy "x" into the appropriate column of the result. ! do i = 1, ncol2 valsym(index6:index6+nrow1-1) = valsym(index2+(i-1)*nrow1:index2+i*nrow1-1) call sge_sl ( valsym(index5), nrow1, nrow1, iwork, valsym(index6), job ) valsym(index4+(i-1)*nrow1:index4+i*nrow1-1) = valsym(index6:index6+nrow1-1) end do ! ! Matrix divided by scalar, A/9, for example. ! else if ( sym == '/' ) then if ( ndim2 == 1 ) then numdim(1,nsyms) = numdim(1,iarg1) numdim(2,nsyms) = numdim(2,iarg1) arg2 = valsym(index2) if ( arg2 == 0.0E+00 ) then ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FUNVAL - Error!' write ( *, '(a)' ) ' Attempt to divide by 0!' return end if valsym(index4:index4+ndim-1) = valsym(index1:index1+ndim1-1) / arg2 else ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FUNVAL - Error!' write ( *, '(a)' ) ' Attempt to divide a matrix by a nonscalar.' end if ! ! Matrix power, A**N or A^N. ! else if ( sym == '**' .or. sym == '^' ) then if ( nrow1 /= ncol1 ) then ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FUNVAL - Error!' write ( *, '(a)' ) ' Matrix must be square for exponentiation!' return end if npow = int ( valsym(index2) ) ! ! We may need to compute the inverse. ! if ( npow <= -1 ) then if ( maxval < ifree + nrow1 * ncol1 - 1 ) then ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FUNVAL - Error!' write ( *, '(a)' ) ' Not enough memory for matrix power!' return end if index5 = ifree + nrow1 * nrow1 call matrix_inverse ( nrow1, nrow1, valsym(index1), valsym(index5), ierror ) if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FUNVAL - Error!' write ( *, '(a)' ) ' The matrix is singular.' write ( *, '(a)' ) ' Negative powers may not be taken!' return end if npow = - npow call r4mat_power ( nrow1, nrow1, valsym(index5), npow, valsym(index4) ) numdim(1,nsyms) = nrow1 numdim(2,nsyms) = ncol1 return end if call r4mat_power ( nrow1, nrow1, valsym(index1), npow, valsym(index4) ) numdim(1,nsyms) = nrow1 numdim(2,nsyms) = nrow1 ! ! Vector or matrix equality, with possibility of implicit ! dimensioning. ! else if ( sym == '=' ) then if ( ndim1 == 0 ) then nrow1 = numdim(1,iarg2) ncol1 = numdim(2,iarg2) numdim(1,iarg1) = nrow1 numdim(2,iarg1) = ncol1 ndim1 = nrow1*ncol1 end if if ( nrow1 /= nrow2 .or. ncol1 /= ncol2 ) then ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FUNVAL - Error!' write ( *, '(a)' ) ' The matrices are incompatible for copying!' return end if call scopy ( ndim1,valsym(index2),1,valsym(index1),1) call scopy ( ndim1,valsym(index2),1,valsym(index4),1) numdim(1,nsyms) = nrow1 numdim(2,nsyms) = ncol1 ! ! ABS: Vector or matrix absolute value. ! else if ( s_eqi ( sym, 'ABS' ) ) then numdim(1,nsyms) = nrow1 numdim(2,nsyms) = ncol1 do j = 1, ncol1*nrow1 valsym(index4) = abs ( valsym(index1) ) index1 = index1 + 1 index4 = index4 + 1 end do ! ! DET: Square matrix determinant. ! else if ( s_eqi ( sym, 'DET' ) ) then if ( nrow1 /= ncol1 ) then ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FUNVAL - Error!' write ( *, '(a)' ) ' The matrix must be square to take the determinant.' return end if call matrix_det ( nrow1, nrow1, valsym(index1), adet ) numdim(1,nsyms) = 1 numdim(2,nsyms) = 1 valsym(index4) = adet ! ! DIAG: Vector or matrix diagonal. ! else if ( s_eqi ( sym, 'DIAG' ) ) then ntemp = min ( nrow1, ncol1 ) numdim(1,nsyms) = ntemp numdim(2,nsyms) = 1 valsym(index4:index4+ntemp-1) = & valsym(index1:index1+(ntemp-1)*(nrow1+1):nrow1+1) ! ! HOUSE: Householder matrix from vector. ! else if ( s_eqi ( sym, 'HOUSE' ) ) then numdim(1,nsyms) = nrow1 numdim(2,nsyms) = nrow1 call house ( valsym(index4), nrow1, nrow1, valsym(index1) ) ! ! ID: Matrix identity. ! else if ( s_eqi ( sym, 'ID' ) ) then call r4mat_identity ( nrow1, nrow1, valsym(index4) ) numdim(1,nsyms) = nrow1 numdim(2,nsyms) = nrow1 ! ! INDEX1: Vector reference. ! else if ( s_eqi ( sym, 'INDEX1' ) ) then i = int ( valsym(index2) ) ipset = index1 + i - 1 numdim(1,nsyms) = 1 numdim(2,nsyms) = 1 ! ! INDEX2: Matrix reference. ! else if ( s_eqi ( sym, 'INDEX2' ) ) then i = int ( valsym(index2) ) j = int ( valsym(index3) ) ipset = index1 + ( j - 1 ) * nrow1 + i - 1 numdim(1,nsyms) = 1 numdim(2,nsyms) = 1 ! ! INT: Vector or matrix integer part (truncation). ! else if ( s_eqi ( sym, 'INT' ) ) then numdim(1,nsyms) = nrow1 numdim(2,nsyms) = ncol1 do j = 1, ncol1*nrow1 valsym(index4) = real ( int ( valsym(index1) ) ) index1 = index1 + 1 index4 = index4 + 1 end do ! ! INV: Matrix inverse. ! else if ( s_eqi ( sym, 'INV' ) ) then if ( nrow1 /= ncol1 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FUNVAL - Error!' write ( *, '(a)' ) ' The matrix must be square to compute the inverse!' ierror = 1 return end if call matrix_inverse ( nrow1, nrow1, valsym(index1), valsym(index4), ierror ) numdim(1,nsyms) = nrow1 numdim(2,nsyms) = ncol1 if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FUNVAL - Error!' write ( *, '(a)' ) ' The matrix is numerically singular!' return end if ! ! KRYLOV: Compute a given number of Krylov vectors. ! ! Form of expression: KRYLOV(A,X,N) ! else if ( s_eqi ( sym, 'KRYLOV' ) ) then if ( numdim(1,index1) /= numdim(2,index1) ) then ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FUNVAL - Error!' write ( *, '(a)' ) ' The matrix A must be square for use with Krylov!' return end if na = numdim(1,index1) nb = int(valsym(index3)) numdim(1,nsyms) = na numdim(2,nsyms) = nb call krylov ( valsym(index1), valsym(index4), na, nb, valsym(index2) ) ! ! NEG: Vector or matrix change sign. ! else if ( s_eqi ( sym, 'NEG' ) ) then numdim(1,nsyms) = nrow1 numdim(2,nsyms) = ncol1 valsym(index4:index4+nrow1*ncol1-1) = - valsym(index1:index1+nrow1*ncol1-1) ! ! NINT: Vector or matrix nearest integer. ! else if ( s_eqi ( sym, 'NINT' ) ) then numdim(1,nsyms) = nrow1 numdim(2,nsyms) = ncol1 do j = 1, ncol1*nrow1 valsym(index4) = anint(valsym(index1)) index1 = index1+1 index4 = index4+1 end do ! ! NORM0: Infinity or maximum norm of a vector or matrix. ! else if ( s_eqi ( sym, 'NORM0' ) ) then numdim(1,nsyms) = 1 numdim(2,nsyms) = 1 temp = 0.0E+00 if ( nrow1 == 1 .or. ncol1 == 1 ) then do i = 1, max(nrow1, ncol1) if ( temp < abs ( valsym(index1) ) ) then temp = abs ( valsym(index1) ) end if index1 = index1+1 end do else do i = 1, nrow1 temp1 = 0 indx = index1+i-1 do j = 1, ncol1 temp1 = temp1 + abs ( valsym(indx) ) indx = indx + nrow1 end do temp = max ( temp, temp1 ) end do end if valsym(index4) = temp ! ! NORM1: L1 norm of a matrix or vector. ! else if ( s_eqi ( sym, 'NORM1' ) ) then numdim(1,nsyms) = 1 numdim(2,nsyms) = 1 temp = 0.0E+00 if ( nrow1 == 1 .or. ncol1 == 1 ) then do i = 1, max(nrow1, ncol1) temp = temp + abs ( valsym(index1) ) index1 = index1 + 1 end do else do j = 1, ncol1 temp1 = 0.0E+00 indx = index1 + (j-1) * nrow1 do i = 1, nrow1 temp1 = temp1 + abs ( valsym(indx) ) indx = indx + 1 end do temp = max ( temp, temp1 ) end do end if valsym(index4) = temp ! ! NORM2: L2 norm of a vector. ! else if ( s_eqi ( sym, 'NORM2' ) ) then if ( nrow1 == 1 .or. ncol1 == 1 ) then ndim=max(nrow1, ncol1) numdim(1,nsyms) = 1 numdim(2,nsyms) = 1 valsym(index4) = 0.0E+00 do i = 1, ndim valsym(index4) = valsym(index4) + valsym(index1+i-1)**2 end do valsym(index4) = sqrt ( valsym(index4) ) else if ( nrow1 == ncol1 ) then if ( maxval < ifree + ( nrow1 + 3 ) * nrow1 - 1 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FUNVAL - Error!' write ( *, '(a)' ) ' There is not enough real free memory' write ( *, '(a)' ) ' for the requested NORM2 operation.' write ( *, '(a)' ) ' The internal value of MAXVAL must be increased' write ( *, '(a)' ) ' in routine COMRPN.' ierror = 1 return end if index5 = index4+nrow1*ncol1 index6 = index5+nrow1 call scopy ( nrow1*nrow1, valsym(index1), 1, value, 1 ) call r4mat_norm2 ( value, nrow1, tnorm ) numdim(1,nsyms) = 1 numdim(2,nsyms) = 1 valsym(index4) = tnorm else ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FUNVAL - Error!' write ( *, '(a)' ) ' The matrix must be square to compute the Two-Norm!' end if ! ! NORME: EISPACK norm of a vector or matrix. ! else if ( s_eqi ( sym, 'NORME' ) ) then numdim(1,nsyms) = 1 numdim(2,nsyms) = 1 valsym(index4) = r4mat_norme ( valsym(index1), nrow1, nrow1, ncol1 ) ! ! NORMF: Frobenius norm of a vector or matrix. ! else if ( s_eqi ( sym, 'NORMF' ) ) then numdim(1,nsyms) = 1 numdim(2,nsyms) = 1 valsym(index4) = r4mat_normf ( valsym(index1), nrow1, nrow1, ncol1 ) ! ! POLY: Scalar or matrix polynomial. ! else if ( s_eqi ( sym, 'POLY' ) ) then if ( ndim1 == 0 .or. ndim2 == 0 .or. ( ncol1 /= 1 .and. nrow1 /= 1 ) .or. & nrow2 /= ncol2 ) then ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FUNVAL - Error!' write ( *, '(a)' ) ' Illegal use of POLY function!' return end if nrow = nrow2 ncol = ncol2 numdim(1,nsyms) = nrow numdim(2,nsyms) = ncol if ( nrow == 1 ) then temp = 0.0E+00 do i = 1, ndim1 coef = valsym(index1+ndim1-i) temp = temp * valsym(index2) + coef end do valsym(index4) = temp else ncoef = max ( nrow1, ncol1 ) call r4mat_poly_val ( valsym(index2), value, valsym(index4), & nrow, valsym(index1), ncoef ) end if ! ! RAN: Vector or matrix randomize. ! else if ( s_eqi ( sym, 'RAN' ) ) then numdim(1,nsyms) = nrow1 numdim(2,nsyms) = ncol1 call random_number ( valsym(index1:index1+nrow1*ncol1-1) ) valsym(index4:index4+nrow1*ncol1-1) = & valsym(index1:index1+nrow1*ncol1-1) ! ! RCOND: Square matrix reciprocal condition number, as estimated ! by LINPACK routine SGE_CO. ! else if ( s_eqi ( sym, 'RCOND' ) ) then if ( nrow1 /= ncol1 ) then ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FUNVAL - Error!' write ( *, '(a)' ) ' The matrix must be square to compute' write ( *, '(a)' ) ' its condition number!' return end if index5 = index4+1 ! ! Make a copy of the matrix A. ! valsym(index5:index5+nrow1*ncol1-1) = & valsym(index1:index1+nrow1*ncol1-1) ! ! Factor the copy and estimate the condition number. ! call sge_co ( valsym(index5), nrow1, nrow1, iwork, rcond ) numdim(1,nsyms) = 1 numdim(2,nsyms) = 1 valsym(index4) = rcond ! ! ROUND: Set to zero all entries that are less than EPS. ! else if ( s_eqi ( sym, 'ROUND' ) ) then valsym(index4:index4+nrow1*ncol1-1) = & valsym(index1:index1+nrow1*ncol1-1) rtol = valsym(irtol) call round ( valsym(index4), nrow1, nrow1, ncol1, rtol ) numdim(1,nsyms) = nrow1 numdim(2,nsyms) = ncol1 ! ! TRACE: Vector or matrix trace. ! else if ( s_eqi ( sym, 'TRACE' ) ) then ntemp = min ( nrow1, ncol1 ) index1 = ipval(iarg1) temp = r4mat_trace ( nrow1, ntemp, valsym(index1) ) numdim(1,nsyms) = 1 numdim(2,nsyms) = 1 valsym(index4) = temp ! ! TRANS: Vector or matrix transpose. ! else if ( s_eqi ( sym, 'TRANS' ) ) then do i = 1, ncol1 call scopy(nrow1,valsym(index1),1,valsym(index4), ncol1) index4 = index4+1 index1=index1+nrow1 end do numdim(1,nsyms) = ncol1 numdim(2,nsyms) = nrow1 ! ! ZEROS: Vector or matrix zero. ! else if ( s_eqi ( sym, 'ZEROS' ) ) then numdim(1,nsyms) = nrow1 numdim(2,nsyms) = nrow1 do i = index4, index4+nrow1*nrow1-1 valsym(i) = 0.0E+00 end do ! ! EVEC, ! EVAL, ! IVAL, ! NORMS, ! or RVAL: Matrix eigen operations. ! else if ( s_eqi ( sym, 'EVEC' ) .or. s_eqi ( sym, 'EVAL' ) .or. & s_eqi ( sym, 'IVAL' ) .or. s_eqi ( sym, 'NORMS' ) .or. & s_eqi ( sym, 'RVAL' ) ) then if ( nrow1 /= ncol1 ) then ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FUNVAL - Error!' write ( *, '(a)' ) ' The matrix must be square for the operation ' // & trim ( sym ) return end if if ( maxval < ifree + ( nrow1 + 3 ) * nrow1 - 1 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FUNVAL - Error!' write ( *, '(a)' ) ' There is not enough real free memory' write ( *, '(a)' ) ' for the requested eigenvalue operation.' write ( *, '(a)' ) ' The internal value of MAXVAL must be increased' write ( *, '(a)' ) ' in routine COMRPN.' ierror = 1 return end if if ( s_eqi ( sym, 'EVEC' ) ) then matz = 1 else matz = 0 end if index5 = index4+nrow1*ncol1 index6 = index5+nrow1 call scopy ( nrow1*nrow1,valsym(index1),1,value,1) call eigen ( value, valsym(index6), valsym(index5), & valsym(index4), ierror, matz, nrow1 ) if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FUNVAL - Warning!' write ( *, '(a)' ) ' Could not compute eigenvalues/vectors' write ( *, '(a,i8)' ) ' 1 through ', ierror ierror = 0 end if if ( s_eqi ( sym, 'EVEC' ) ) then numdim(1,nsyms) = nrow1 numdim(2,nsyms) = ncol1 else if ( s_eqi ( sym, 'EVAL' ) ) then numdim(1,nsyms) = nrow1 numdim(2,nsyms) = 2 call scopy(nrow1,valsym(index5),1,valsym(index4),1) call scopy(nrow1,valsym(index6),1,valsym(index4+nrow1),1) else if ( s_eqi ( sym, 'NORMS' ) ) then numdim(1,nsyms) = 1 numdim(2,nsyms) = 1 temp = 0.0E+00 do i = 1, nrow1 temp1 = sqrt(valsym(index5+i-1)**2+valsym(index6+i-1)**2) temp = max ( temp, temp1 ) end do valsym(index4) = temp else numdim(1,nsyms) = nrow1 numdim(2,nsyms) = 1 if ( s_eqi ( sym, 'RVAL' ) ) then call scopy ( nrow1, valsym(index5), 1, valsym(index4), 1 ) else call scopy ( nrow1, valsym(index6), 1, valsym(index4), 1 ) end if end if ! ! Unknown operation! ! else ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FUNVAL - Error!' write ( *, '(a)' ) ' Unknown operation: ' // trim ( sym ) end if return end subroutine hello ( maxrow ) !*****************************************************************************80 ! !! HELLO prints out an introductory message for MATALG. ! ! Modified: ! ! 18 March 2003 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer MAXROW, the maximum number of rows and columns. ! implicit none integer maxrow write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'MATALG' write ( *, '(a)' ) ' Version 1.43' write ( *, '(a)' ) ' Last modified on 18 March 2003' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' MATALG evaluates formulas involving matrices.' write ( *, '(a,i8)' ) ' Matrices may have order up to ', maxrow write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' It was developed by Charles Cullen and John Burkardt.' write ( *, '(a)' ) ' All rights reserved by the authors. This program may' write ( *, '(a)' ) ' not be reproduced without written permission.' return end subroutine help !*****************************************************************************80 ! !! HELP prints the help menu. ! ! Modified: ! ! 27 October 2000 ! ! Author: ! ! John Burkardt ! implicit none write ( *, '(a)' ) 'DEBUG Print debugging info.' write ( *, '(a)' ) 'E variable, nrow, ncol Enter a new variable' write ( *, '(a)' ) 'I Initialize.' write ( *, '(a)' ) 'K v Kill variable(s).' write ( *, '(a)' ) 'L v, xlo, xhi, nstep Tabulate a formula' write ( *, '(a)' ) 'M PLU/QR/SVD, variable Matrix factors.' write ( *, '(a)' ) 'P variable Characteristic polynomial.' write ( *, '(a)' ) 'Q Quit' write ( *, '(a)' ) 'T name/ALL/blank Type variable' write ( *, '(a)' ) 'U name Partition a matrix.' write ( *, '(a)' ) 'V Re-evaluate formula.' write ( *, '(a)' ) '# comment Begin a comment.' write ( *, '(a)' ) '$, % Turn paging on or off.' write ( *, '(a)' ) '< file Get input from file.' write ( *, '(a)' ) 'variable = formula Set variable to value of formula.' return end subroutine house ( a, lda, n, u ) !*****************************************************************************80 ! !! HOUSE constructs the Householder elementary reflector for the vector U. ! ! Formula: ! ! A = I - (2*U*transpose(U))/(tranpose(U)*U) ! ! Example: ! ! U = (1, 1, 1, 0, -1) ! ! 1/2 -1/2 -1/2 0 1/2 ! -1/2 1/2 -1/2 0 1/2 ! -1/2 -1/2 1/2 0 1/2 ! 0 0 0 1 0 ! 1/2 1/2 1/2 0 1/2 ! ! Properties: ! ! A is symmetric. ! A is orthogonal (Transpose(A)*A=I) ! If Y and Z are nonzero vectors of equal length, and ! U=(Y-Z)/NORM(Y-Z), then A*y=z. ! The matrix A represents a reflection through the plane which ! is perpendicular to the vector U. In particular, A*U=-U. ! ! Modified: ! ! 15 March 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, real A(LDA,N), the Householder matrix. ! ! Input, integer LDA, the leading dimension of A. ! ! Input, integer N, the order of A. ! ! Input, real U(N), the vector that defines the Householder matrix. ! implicit none integer lda integer n real a(lda,n) integer i integer j real u(n) real udot ! ! Compute U dot U, for normalization. ! udot = dot_product ( u(1:n), u(1:n) ) ! ! Set A to the identity. ! call r4mat_identity ( lda, n, a ) if ( udot /= 0.0E+00 ) then do i = 1, n do j = 1, n a(i,j) = a(i,j) - 2.0E+00 * u(i) * u(j) / udot end do end do end if return end subroutine hqr ( lda, n, low, igh, h, wr, wi, ierror ) !*****************************************************************************80 ! !! HQR finds the eigenvalues of a real upper Hessenberg matrix by the QR method. ! ! Modified: ! ! 15 March 2001 ! ! Parameters: ! ! Input, integer LDA, the leading dimension of H. LDA must be ! at least N. ! ! Input, integer N, the number of rows and columns in the matrix. ! ! Input, integer LOW, IGH, column indices determined by ! BALANC. If BALANC is not used, set LOW = 1, IGH=N. ! ! Input/output, real H(LDA,N). ! ! On input, H contains the upper Hessenberg matrix. Information ! about the transformations used in the reduction to Hessenberg ! form by ELMHES or ORTHES, if performed, is stored ! in the remaining triangle under the Hessenberg matrix. ! ! On output, the information that was in H has been destroyed. ! ! Output, real WR(N), WI(N), contain the real and imaginary parts, ! respectively, of the eigenvalues. The eigenvalues ! are unordered except that complex conjugate pairs ! of values appear consecutively with the eigenvalue ! having the positive imaginary part first. If an ! error exit is made, the eigenvalues should be correct ! for indices IERROR+1 through N. ! ! Output, integer IERROR, error flag. ! 0, no error occurred. ! J, if the limit of 30*N iterations is exhausted ! while the J-th eigenvalue is being sought. ! implicit none integer lda integer n integer en integer enm2 real h(lda,n) real hnorm integer i integer ierror integer igh integer itn integer its integer j integer k integer l integer ll integer low integer m integer mm integer na logical notlas real p real q real r real s real t real tst1 real tst2 real w real wi(n) real wr(n) real x real y real zz ierror = 0 ! ! Compute the norm of the upper Hessenberg matrix. ! hnorm = 0.0E+00 do i = 1, n do j = max ( i-1, 1 ), n hnorm = hnorm + abs ( h(i,j) ) end do end do ! ! Store roots isolated by BALANC. ! do i = 1, n if (i < low .or. igh < i ) then wr(i) = h(i,i) wi(i) = 0.0E+00 end if end do en = igh t = 0.0E+00 itn = 30 * n ! ! Search for next eigenvalues. ! 60 continue if ( en < low ) then return end if its = 0 na = en - 1 enm2 = na - 1 ! ! Look for single small sub-diagonal element. ! 70 continue do ll = low, en l = en + low - ll if ( l == low ) then exit end if s = abs ( h(l-1,l-1) ) + abs ( h(l,l) ) if ( s == 0.0E+00 ) then s = hnorm end if tst1 = s tst2 = tst1 + abs ( h(l,l-1) ) if ( tst2 == tst1 ) then exit end if end do ! ! Form shift. ! x = h(en,en) if ( l == en ) then wr(en) = x + t wi(en) = 0.0E+00 en = na go to 60 end if y = h(na,na) w = h(en,na) * h(na,en) if ( l == na) then go to 280 end if if ( itn == 0 ) then ierror = en return end if ! ! Form exceptional shift. ! if ( its == 10 .or. its == 20 ) then t = t + x do i = low, en h(i,i) = h(i,i) - x end do s = abs ( h(en,na) ) + abs ( h(na,enm2) ) x = 0.75E+00 * s y = x w = -0.4375E+00 * s * s end if its = its + 1 itn = itn - 1 ! ! Look for two consecutive small sub-diagonal elements. ! do mm = l, enm2 m = enm2 + l - mm zz = h(m,m) r = x - zz s = y - zz p = (r * s - w) / h(m+1,m) + h(m,m+1) q = h(m+1,m+1) - zz - r - s r = h(m+2,m+1) s = abs ( p ) + abs ( q ) + abs ( r ) p = p / s q = q / s r = r / s if ( m == l ) then exit end if tst1 = abs ( p ) * ( abs ( h(m-1,m-1) ) + abs ( zz ) + abs ( h(m+1,m+1) ) ) tst2 = tst1 + abs ( h(m,m-1) ) * ( abs ( q ) + abs ( r ) ) if ( tst2 == tst1 ) then exit end if end do do i = m+2, en h(i,i-2) = 0.0E+00 if ( i /= m+2 ) then h(i,i-3) = 0.0E+00 end if end do ! ! Double QR step involving rows l to en and columns m to en. ! do k = m, na notlas = k /= na if ( k /= m ) then p = h(k,k-1) q = h(k+1,k-1) r = 0.0E+00 if ( notlas ) then r = h(k+2,k-1) end if x = abs ( p ) + abs ( q ) + abs ( r ) if ( x == 0.0E+00 ) then cycle end if p = p / x q = q / x r = r / x end if s = sign ( sqrt ( p*p+q*q+r*r), p ) if ( k /= m ) then h(k,k-1) = -s * x else if ( l /= m ) then h(k,k-1) = -h(k,k-1) end if end if p = p + s x = p / s y = q / s zz = r / s q = q / p r = r / p if ( .not. notlas ) then ! ! Row modification. ! do j = k, n p = h(k,j) + q * h(k+1,j) h(k,j) = h(k,j) - p * x h(k+1,j) = h(k+1,j) - p * y end do j = min ( en, k+3 ) ! ! Column modification. ! do i = 1, j p = x * h(i,k) + y * h(i,k+1) h(i,k) = h(i,k) - p h(i,k+1) = h(i,k+1) - p * q end do else ! ! Row modification. ! do j = k, n p = h(k,j) + q * h(k+1,j) + r * h(k+2,j) h(k,j) = h(k,j) - p * x h(k+1,j) = h(k+1,j) - p * y h(k+2,j) = h(k+2,j) - p * zz end do j = min(en,k+3) ! ! Column modification. ! do i = 1, j p = x * h(i,k) + y * h(i,k+1) + zz * h(i,k+2) h(i,k) = h(i,k) - p h(i,k+1) = h(i,k+1) - p * q h(i,k+2) = h(i,k+2) - p * r end do end if end do go to 70 ! ! Two roots found. ! 280 continue p = ( y - x ) / 2.0E+00 q = p * p + w zz = sqrt ( abs ( q ) ) x = x + t ! ! Real pair. ! if ( 0.0E+00 <= q ) then zz = p + sign(zz,p) wr(na) = x + zz if ( zz == 0.0E+00 ) then wr(en) = wr(na) else wr(en) = x - w / zz end if wi(na) = 0.0E+00 wi(en) = 0.0E+00 ! ! Complex pair. ! else wr(na) = x + p wr(en) = x + p wi(na) = zz wi(en) = -zz end if ! ! Deduct the two eigenvalues we have found from the total to ! be found, and proceed. ! en = enm2 go to 60 end subroutine hqr2 ( lda, n, low, igh, h, wr, wi, z, ierror ) !*****************************************************************************80 ! !! HQR2 finds the eigenvalues and eigenvectors of a real upper Hessenberg matrix by the QR method. ! ! Modified: ! ! 15 March 2001 ! ! Parameters: ! ! Input, integer LDA, the leading dimension of H and Z. ! LDA must be at least N. ! ! Input, integer N, the number of rows and columns in the ! matrix H. ! ! Input, integer LOW, IGH, column indices determined by ! BALANC. If BALANC is not used, set LOW = 1, IGH=N. ! ! Input/output, real H(LDA,N). ! On input, H contains the upper Hessenberg matrix. ! On output, the information that was in H has been destroyed. ! ! Output, real WR(N), WI(N), contain the real and imaginary parts, ! respectively, of the eigenvalues. The eigenvalues ! are unordered except that complex conjugate pairs ! of values appear consecutively with the eigenvalue ! having the positive imaginary part first. If an ! error exit is made, the eigenvalues should be correct ! for indices IERROR+1 through N. ! ! Input/output, real Z(LDA,N). ! ! On input, Z contains the transformation matrix produced by ! ELTRAN after the reduction by ELMHES if performed. If the ! eigenvectors of the Hessenberg matrix are desired, Z must ! contain the identity matrix. ! ! On output, Z contains the real and imaginary parts of the ! eigenvectors. If the I-th eigenvalue is real, the I-th column ! of Z contains its eigenvector. If the I-th eigenvalue is complex ! with positive imaginary part, the I-th and (I+1)-th ! columns of Z contain the real and imaginary parts of its ! eigenvector. The eigenvectors are unnormalized. If an ! error exit is made, none of the eigenvectors has been found. ! ! Output, integer IERROR, error flag. ! 0, no error occurred. ! J, if the limit of 30*N iterations is exhausted ! while the J-th eigenvalue is being sought. ! implicit none integer lda integer n integer en integer enm2 real h(lda,n) real hnorm integer i integer ierror integer igh integer ii integer itn integer its integer j integer jj integer k integer l integer ll integer low integer m integer mm integer na integer nn logical notlas real p real q real r real ra real s real sa real t real temp real tst1 real tst2 real vi real vr real w real wi(n) real wr(n) real x real y real z(lda,n) real zz ierror = 0 ! ! Compute the norm of the upper Hessenberg matrix. ! hnorm = 0.0E+00 do i = 1, n do j = max ( i-1, 1 ), n hnorm = hnorm + abs ( h(i,j) ) end do end do ! ! Store roots isolated by BALANC. ! do i = 1, n if ( i < low .or. igh < i ) then wr(i) = h(i,i) wi(i) = 0.0E+00 end if end do en = igh t = 0.0E+00 itn = 30 * n ! ! Search for next eigenvalues ! 60 continue if ( en < low ) then go to 340 end if its = 0 na = int(en-1) enm2 = na - 1 ! ! Look for single small sub-diagonal element. ! 70 continue do ll = low, en l = en + low - ll if ( l == low ) then exit end if s = abs ( h(l-1,l-1) ) + abs ( h(l,l) ) if ( s == 0.0E+00 ) then s = hnorm end if tst1 = s tst2 = tst1 + abs ( h(l,l-1) ) if ( tst2 == tst1 ) then exit end if end do ! ! Form shift. ! x = h(en,en) if ( l == en ) then go to 270 end if y = h(na,na) w = h(en,na) * h(na,en) if ( l == na ) then go to 280 end if if ( itn == 0 ) then ierror = en return end if ! ! Form exceptional shift. ! if ( its == 10 .or. its == 20 ) then t = t + x do i = low, en h(i,i) = h(i,i) - x end do s = abs ( h(en,na) ) + abs ( h(na,enm2) ) x = 0.75E+00 * s y = x w = -0.4375E+00 * s * s end if its = its + 1 itn = itn - 1 ! ! Look for two consecutive small sub-diagonal elements. ! do mm = l, enm2 m = enm2 + l - mm zz = h(m,m) r = x - zz s = y - zz p = (r * s - w) / h(m+1,m) + h(m,m+1) q = h(m+1,m+1) - zz - r - s r = h(m+2,m+1) s = abs ( p ) + abs ( q ) + abs ( r ) p = p / s q = q / s r = r / s if ( m == l) then exit end if tst1 = abs ( p ) * ( abs ( h(m-1,m-1) ) + abs ( zz ) + abs ( h(m+1,m+1))) tst2 = tst1 + abs ( h(m,m-1) ) * ( abs ( q ) + abs ( r ) ) if ( tst2 == tst1 ) then exit end if end do do i = m+2, en h(i,i-2) = 0.0E+00 if ( i /= m+2 ) then h(i,i-3) = 0.0E+00 end if end do ! ! Double QR step involving rows L to EN and columns M to EN. ! do k = m, na notlas = k /= na if ( k /= m ) then p = h(k,k-1) q = h(k+1,k-1) r = 0.0E+00 if ( notlas ) r = h(k+2,k-1) x = abs ( p ) + abs ( q ) + abs ( r ) if ( x == 0.0E+00 ) then cycle end if p = p / x q = q / x r = r / x end if s = sign ( sqrt ( p*p + q*q + r*r ), p ) if ( k /= m ) then h(k,k-1) = -s * x else if ( l /= m ) h(k,k-1) = -h(k,k-1) end if p = p + s x = p / s y = q / s zz = r / s q = q / p r = r / p if ( .not. notlas ) then ! ! Row modification. ! do j = k, n p = h(k,j) + q * h(k+1,j) h(k,j) = h(k,j) - p * x h(k+1,j) = h(k+1,j) - p * y end do j = min ( en, k + 3 ) ! ! Column modification. ! do i = 1, j p = x * h(i,k) + y * h(i,k+1) h(i,k) = h(i,k) - p h(i,k+1) = h(i,k+1) - p * q end do ! ! Accumulate transformations. ! do i = low, igh p = x * z(i,k) + y * z(i,k+1) z(i,k) = z(i,k) - p z(i,k+1) = z(i,k+1) - p * q end do else ! ! Row modification. ! do j = k, n p = h(k,j) + q * h(k+1,j) + r * h(k+2,j) h(k,j) = h(k,j) - p * x h(k+1,j) = h(k+1,j) - p * y h(k+2,j) = h(k+2,j) - p * zz end do j=min(en,k+3) ! ! Column modification. ! do i = 1, j p = x * h(i,k) + y * h(i,k+1) + zz * h(i,k+2) h(i,k) = h(i,k) - p h(i,k+1) = h(i,k+1) - p * q h(i,k+2) = h(i,k+2) - p * r end do ! ! Accumulate transformations. ! do i = low, igh p = x * z(i,k) + y * z(i,k+1) + zz * z(i,k+2) z(i,k) = z(i,k) - p z(i,k+1) = z(i,k+1) - p * q z(i,k+2) = z(i,k+2) - p * r end do end if end do go to 70 ! ! One root found. ! 270 continue h(en,en) = x + t wr(en) = h(en,en) wi(en) = 0.0E+00 en = na go to 60 ! ! Two roots found. ! 280 p = ( y - x ) / 2.0E+00 q = p * p + w zz = sqrt ( abs ( q ) ) h(en,en) = x + t x = h(en,en) h(na,na) = y + t if ( q < 0.0E+00 ) then go to 320 end if ! ! Real pair. ! zz = p + sign(zz,p) wr(na) = x + zz wr(en)=wr(na) if ( zz /= 0.0E+00 ) wr(en) = x - w / zz wi(na) = 0.0E+00 wi(en) = 0.0E+00 x = h(en,na) s = abs ( x ) + abs ( zz ) p = x / s q = zz / s r = sqrt(p*p+q*q) p = p / r q = q / r ! ! Row modification. ! do j = na, n zz = h(na,j) h(na,j) = q * zz + p * h(en,j) h(en,j) = q * h(en,j) - p * zz end do ! ! Column modification. ! do i = 1, en zz = h(i,na) h(i,na) = q * zz + p * h(i,en) h(i,en) = q * h(i,en) - p * zz end do ! ! Accumulate transformations. ! do i = low, igh zz = z(i,na) z(i,na) = q * zz + p * z(i,en) z(i,en) = q * z(i,en) - p * zz end do go to 330 ! ! Complex pair. ! 320 continue wr(na) = x + p wr(en) = x + p wi(na) = zz wi(en) = -zz 330 continue en = enm2 go to 60 ! ! All roots found. Backsubstitute to find vectors of upper ! triangular form. ! 340 continue if ( hnorm == 0.0E+00 ) then return end if do nn = 1, n en = n + 1 - nn p=wr(en) q=wi(en) na = en - 1 if ( q < 0.0E+00 ) then go to 710 end if if ( 0.0E+00 < q ) then go to 800 end if ! ! Real vector. ! m = en h(en,en) = 1.0E+00 if ( na == 0 ) then go to 800 end if do ii = 1, na i = en - ii w = h(i,i) - p r = 0.0E+00 do j = m, en r = r + h(i,j) * h(j,en) end do if ( wi(i) < 0.0E+00 ) then zz = w s = r go to 700 end if m = i if ( wi(i) == 0.0E+00 ) then t = w if ( t == 0.0E+00 ) then tst1 = hnorm t = tst1 632 continue t = 0.01E+00 * t tst2 = hnorm + t if ( tst1 < tst2 ) then go to 632 end if end if h(i,en) = -r / t go to 680 end if ! ! Solve real equations. ! x = h(i,i+1) y = h(i+1,i) q = (wr(i) - p) * (wr(i) - p) + wi(i) * wi(i) t = (x * s - zz * r) / q h(i,en) = t if ( abs ( zz ) < abs ( x ) ) then h(i+1,en) = (-r - w * t) / x else h(i+1,en) = (-s - y * t) / zz end if ! ! Overflow control. ! 680 continue t = abs ( h(i,en) ) if ( t == 0.0E+00 ) then go to 700 end if tst1 = t tst2 = tst1 + 1.0E+00 / tst1 if ( tst2 <= tst1 ) then do j = i, en h(j,en)=h(j,en)/t end do end if 700 continue end do ! ! End real vector. ! go to 800 ! ! Complex vector. ! 710 continue m = na ! ! Last vector component chosen imaginary so that ! eigenvector matrix is triangular. ! if ( abs ( h(na,en) ) < abs ( h(en,na) ) ) then h(na,na) = q / h(en,na) h(na,en) = -(h(en,en) - p) / h(en,na) else temp = 0.0E+00 call cdiv(temp,-h(na,en),h(na,na)-p,q,h(na,na),h(na,en)) end if h(en,na) = 0.0E+00 h(en,en) = 1.0E+00 enm2 = na - 1 do ii = 1, enm2 i = na - ii w = h(i,i) - p ra = 0.0E+00 sa = 0.0E+00 do j = m, en ra = ra + h(i,j) * h(j,na) sa = sa + h(i,j) * h(j,en) end do if ( wi(i) < 0.0E+00 ) then zz=w r = ra s = sa go to 795 end if m = i if ( wi(i) == 0.0E+00 ) then call cdiv ( -ra, -sa, w, q, h(i,na), h(i,en) ) go to 790 end if ! ! Solve complex equations. ! x = h(i,i+1) y = h(i+1,i) vr = (wr(i) - p) * (wr(i) - p) + wi(i) * wi(i) - q * q vi = (wr(i) - p) * 2.0E+00 * q if ( vr == 0.0E+00 .and. vi == 0.0E+00 ) then tst1 = hnorm * ( abs ( w ) + abs ( q ) + abs ( x ) + abs ( y ) & + abs ( zz ) ) vr = tst1 783 continue vr = 0.01E+00 * vr tst2 = tst1 + vr if ( tst1 < tst2 ) then go to 783 end if end if call cdiv ( x*r-zz*ra+q*sa, x*s-zz*sa-q*ra, vr, vi, h(i,na), h(i,en) ) if ( abs ( zz ) + abs ( q ) < abs ( x ) ) then h(i+1,na) = (-ra - w * h(i,na) + q * h(i,en)) / x h(i+1,en) = (-sa - w * h(i,en) - q * h(i,na)) / x else call cdiv ( -r-y*h(i,na), -s-y*h(i,en), zz, q, h(i+1,na), h(i+1,en) ) end if ! ! Overflow control ! 790 continue t = max ( abs ( h(i,na) ), abs ( h(i,en) ) ) if ( t /= 0.0E+00 ) then tst1 = t tst2 = tst1 + 1.0E+00 / tst1 if ( tst2 <= tst1 ) then h(i:en,na) = h(i:en,na) / t h(i:en,en) = h(i:en,en) / t end if end if 795 continue end do ! ! End complex vector ! 800 continue end do ! ! End back substitution. ! ! Vectors of isolated roots ! do i = 1, n if ( i < low .or. igh < i ) then z(i,i:n) = h(i,i:n) end if end do ! ! Multiply by transformation matrix to give ! vectors of original full matrix. ! do jj = low, n j = n + low - jj m = min ( j, igh ) do i = low, igh zz = 0.0E+00 do k = low, m zz = zz + z(i,k) * h(k,j) end do z(i,j) = zz end do end do return end function i4_gcd ( i, j ) !*****************************************************************************80 ! !! I4_GCD finds the greatest common divisor of I and J. ! ! Modified: ! ! 03 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer I, J, two numbers whose greatest common divisor ! is desired. ! ! Output, integer I4_GCD, the greatest common divisor of I and J. ! ! Note that only the absolute values of I and J are ! considered, so that the result is always nonnegative. ! ! If I or J is 0, I4_GCD is returned as max ( 1, abs ( I ), abs ( J ) ). ! ! If I and J have no common factor, I4_GCD is returned as 1. ! ! Otherwise, using the Euclidean algorithm, I4_GCD is the ! largest common factor of I and J. ! implicit none integer i integer i4_gcd integer ip integer iq integer ir integer j i4_gcd = 1 ! ! Return immediately if either I or J is zero. ! if ( i == 0 ) then i4_gcd = max ( 1, abs ( j ) ) return else if ( j == 0 ) then i4_gcd = max ( 1, abs ( i ) ) return end if ! ! Set IP to the larger of I and J, IQ to the smaller. ! This way, we can alter IP and IQ as we go. ! ip = max ( abs ( i ), abs ( j ) ) iq = min ( abs ( i ), abs ( j ) ) ! ! Carry out the Euclidean algorithm. ! do ir = mod ( ip, iq ) if ( ir == 0 ) then exit end if ip = iq iq = ir end do i4_gcd = iq return end function i4_lcm ( i, j ) !*****************************************************************************80 ! !! I4_LCM computes the least common multiple of two integers. ! ! Definition: ! ! The least common multiple may be defined as ! ! LCM(I,J) = ABS( I * J ) / GCD(I,J) ! ! where GCD(I,J) is the greatest common factor of I and J. ! ! Modified: ! ! 01 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer I, J, the integers whose I4_LCM is desired. ! ! Output, integer I4_LCM, the least common multiple of I and J. ! I4_LCM is never negative. I4_LCM is 0 if either I or J is zero. ! implicit none integer i integer i4_gcd integer j integer i4_lcm i4_lcm = abs ( i * ( j / i4_gcd ( i, j ) ) ) return end subroutine i4_to_s_zero ( intval, s ) !*****************************************************************************80 ! !! I4_TO_S_ZERO converts an integer to a string, with zero padding. ! ! Examples: ! ! Assume that S is 6 characters long: ! ! INTVAL S ! ! 1 000001 ! -1 -00001 ! 0 000000 ! 1952 001952 ! 123456 123456 ! 1234567 ****** <-- Not enough room! ! ! Modified: ! ! 04 August 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer INTVAL, an integer to be converted. ! ! Output, character ( len = * ) S, the representation of the integer. ! The integer will be right justified, and zero padded. ! If there is not enough space, the string will be filled with stars. ! implicit none character c integer i integer idig integer ihi integer ilo integer intval integer ipos integer ival character ( len = * ) s s = ' ' ilo = 1 ihi = len ( s ) if ( ihi <= 0 ) then return end if ! ! Make a copy of the integer. ! ival = intval ! ! Handle the negative sign. ! if ( ival < 0 ) then if ( ihi <= 1 ) then s(1:1) = '*' return end if ival = - ival s(1:1) = '-' ilo = 2 end if ! ! Working from right to left, strip off the digits of the integer ! and place them into S(ILO:IHI). ! ipos = ihi do while ( ival /= 0 .or. ipos == ihi ) idig = mod ( ival, 10 ) ival = ival / 10 if ( ipos < ilo ) then do i = 1, ihi s(i:i) = '*' end do return end if call digit_to_ch ( idig, c ) s(ipos:ipos) = c ipos = ipos - 1 end do ! ! Fill the empties with zeroes. ! do i = ilo, ipos s(i:i) = '0' end do return end subroutine inform ( ierror, infix, irpn, line, maxrow, maxrpn, & ncol, nform, nrow, value ) !*****************************************************************************80 ! !! INFORM accepts a formula, compiles, evaluates, and displays the value. ! ! Modified: ! ! 17 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer IERROR, error flag. ! 0, no error occurred. ! nonzero, an error occurred. ! ! Workspace, integer IRPN(MAXRPN), used to store the compiled ! versions of user formulas. ! ! Input/output, character ( len = 80 ) LINE, a buffer containing ! unprocessed user input. ! ! Input, integer MAXROW, the maximum number of rows and columns. ! ! Input, integer MAXRPN, specifies the length of IRPN. ! ! Output, integer NCOL, the number of columns in the result. ! ! Output, integer NFORM. ! 0, if the formula could not be read or compiled. ! 1, if the formula was read and compiled. ! ! Output, integer NROW, the number of rows in the result. ! ! Workspace, real VALUE(MAXROW,MAXROW), space used to hold ! the value of the variable to be saved. ! implicit none integer, parameter :: maxlen = 20 integer maxrow integer maxrpn character com integer ierror integer ifrm integer iloc character ( len = * ) infix integer irpn(maxrpn) character ( len = 80 ) line character ( len = maxlen ) namtmp character ( len = maxlen ) namvar integer ncol integer nform integer nrow character ( len = 80 ) prompt real value(maxrow,maxrow) namvar = ' ' nform = 0 nrow = 0 ncol = 0 ! ! Get INFIX, the formula as the user types it. ! prompt = 'formula.' call chrrea ( infix, line, prompt, ierror, 0 ) if ( ierror /= 0 ) then return end if ! ! Remove all blanks from the formula. ! call s_blank_delete ( infix ) ! ! Compile the formula. ! com = 'F' ifrm = 1 call comrpn ( com, ierror, ifrm, infix, irpn, maxrow, maxrpn, & namvar, ncol, nrow, value ) if ( ierror /= 0 ) then return end if nform = 1 ! ! Evaluate the formula. ! com = 'E' call comrpn ( com, ierror, ifrm, infix, irpn, maxrow, maxrpn, & namvar, ncol, nrow, value ) if ( ierror /= 0 ) then return end if ! ! Print the formula's value. ! ! If the formula is an assignment, then use the assigned ! variable name as a label; otherwise, use a dummy label. ! iloc = index ( infix, '=' ) if ( 1 < iloc ) then namtmp = infix(1:iloc-1) else namtmp = 'Formula' end if call r4mat_print ( maxrow, nrow, ncol, value, namtmp ) return end subroutine inicom ( ifinis, ifree, indx1, indx2, ineg, infix, intsym, iopsym, & iprsym, ipval, irad, irpn, irtol, istack, maxrow, maxrpn, maxsym, & maxval, nints, nsym, nsymp, nsyms, numdim, symbol, valsym, value ) !*****************************************************************************80 ! !! INICOM initializes data for COMRPN. ! ! Modified: ! ! 02 October 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer IFINIS, the index in SYMBOL of the symbol ! '$', meaning the end of the formula. ! ! Output, integer IFREE, the index of the first free address in VALSYM. ! ! Output, integer INDX1, the index of the symbol 'INDEX1'. ! ! Output, integer INDX2, the index of the symbol 'INDEX2'. ! ! Output, integer INEG, the index of the symbol 'NEG'. ! ! Output, character ( len = MAXFIX ) INFIX, space to hold the text ! of an infix formula typed in by the user. ! ! Output, integer INTSYM(80), a set of integers which ! are the indices of symbols, representing an infix formula. ! ! Output, integer IPVAL(MAXSYM), contains, for each symbol, ! the address in VALSYM where the value of the symbol ! is stored, if it is a scalar. If the symbol represents ! a vector or matrix, IPVAL(IARG) points to location of ! the first entry. ! ! Output, integer IRAD, the index of the symbol 'ANGLE_TO_RADIAN'. ! ! Output, integer IRPN(MAXRPN), used to store the compiled ! versions of user formulas. ! ! Workspace, integer ISTACK(MAXSYM), workspace for interpreting ! the formula. ! ! Input, integer MAXROW, the maximum number of rows and columns. ! ! Input, integer MAXRPN, specifies the length of IRPN. ! ! Input, integer MAXSYM, the maximum number of symbols allowed. ! ! Input, integer MAXVAL, the maximum number of values allowed. ! ! Output, integer NINTS, the number of useful entries in INTSYM. ! ! Output, integer NSYM, the total number of symbols, including temporaries. ! ! Output, integer NSYMP, the number of permanent symbols. ! ! Output, integer NYMS, the number of declared symbols. ! ! Output, integer NUMDIM(2,MAXSYM). ! For each symbol I, NUMDIM(1,I) is the number of rows in its ! value, and NUMDIM(2,I) is the number of columns. ! ! Output, character ( len = MAXLEN ) SYMBOL(MAXSYM), the symbolic names ! of the internally defined constants and functions. ! ! Output, real VALSYM(MAXVAL), the values of the symbolic variables. ! ! Workspace, real VALUE(MAXROW,MAXROW), space used to hold ! the value of the variable to be saved. ! implicit none integer, parameter :: maxlen = 20 integer maxrow integer maxrpn integer maxsym integer maxval integer i integer ifinis integer ifree integer indx1 integer indx2 integer ineg character ( len = * ) infix integer intsym(80) integer iopsym(maxsym) integer iprsym(maxsym) integer ipval(maxsym) integer irad integer irpn(maxrpn) integer irtol integer istack(maxsym) integer nints integer nsym integer nsymp integer nsyms integer numdim(2,maxsym) logical s_eqi character ( len = maxlen ) symbol(maxsym) real valsym(maxval) real value(maxrow,maxrow) nsym = 0 nsym = nsym + 1 symbol(nsym) = 'ABS' iopsym(nsym) = 1 iprsym(nsym) = 9 numdim(1,nsym) = 1 numdim(2,nsym) = 1 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'ACOS' iopsym(nsym) = 1 iprsym(nsym) = 9 numdim(1,nsym) = 1 numdim(2,nsym) = 1 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'ALOG' iopsym(nsym) = 1 iprsym(nsym) = 9 numdim(1,nsym) = 1 numdim(2,nsym) = 1 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'ALOG10' iopsym(nsym) = 1 iprsym(nsym) = 9 numdim(1,nsym) = 1 numdim(2,nsym) = 1 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'ANGLE_TO_RADIAN' iopsym(nsym) = 0 iprsym(nsym) = 10 numdim(1,nsym) = 1 numdim(2,nsym) = 1 valsym(nsym) = 1.0E+00 nsym = nsym + 1 symbol(nsym) = 'ASIN' iopsym(nsym) = 1 iprsym(nsym) = 9 numdim(1,nsym) = 1 numdim(2,nsym) = 1 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'ATAN' iopsym(nsym) = 1 iprsym(nsym) = 9 numdim(1,nsym) = 1 numdim(2,nsym) = 1 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'ATAN2' iopsym(nsym) = 2 iprsym(nsym) = 9 numdim(1,nsym) = 1 numdim(2,nsym) = 1 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'COS' iopsym(nsym) = 1 iprsym(nsym) = 9 numdim(1,nsym) = 1 numdim(2,nsym) = 1 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'COSH' iopsym(nsym) = 1 iprsym(nsym) = 9 numdim(1,nsym) = 1 numdim(2,nsym) = 1 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'COT' iopsym(nsym) = 1 iprsym(nsym) = 9 numdim(1,nsym) = 1 numdim(2,nsym) = 1 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'CSC' iopsym(nsym) = 1 iprsym(nsym) = 9 numdim(1,nsym) = 1 numdim(2,nsym) = 1 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'DET' iopsym(nsym) = 1 iprsym(nsym) = 9 numdim(1,nsym) = 1 numdim(2,nsym) = 1 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'DIAG' iopsym(nsym) = 1 iprsym(nsym) = 9 numdim(1,nsym) = 1 numdim(2,nsym) = 1 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'E' iopsym(nsym) = 0 iprsym(nsym) = 10 numdim(1,nsym) = 1 numdim(2,nsym) = 1 valsym(nsym) = exp ( 1.0E+00 ) nsym = nsym + 1 symbol(nsym) = 'EPS' iopsym(nsym) = 0 iprsym(nsym) = 10 numdim(1,nsym) = 1 numdim(2,nsym) = 1 valsym(nsym) = epsilon ( 1.0E+00 ) nsym = nsym + 1 symbol(nsym) = 'EVAL' iopsym(nsym) = 1 iprsym(nsym) = 9 numdim(1,nsym) = 1 numdim(2,nsym) = 1 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'EVEC' iopsym(nsym) = 1 iprsym(nsym) = 9 numdim(1,nsym) = 1 numdim(2,nsym) = 1 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'EXP' iopsym(nsym) = 1 iprsym(nsym) = 9 numdim(1,nsym) = 1 numdim(2,nsym) = 1 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'GCD' iopsym(nsym) = 2 iprsym(nsym) = 9 numdim(1,nsym) = 1 numdim(2,nsym) = 1 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'HOUSE' iopsym(nsym) = 1 iprsym(nsym) = 9 numdim(1,nsym) = 1 numdim(2,nsym) = 1 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'HUGE_I' iopsym(nsym) = 0 iprsym(nsym) = 10 numdim(1,nsym) = 1 numdim(2,nsym) = 1 valsym(nsym) = real ( huge ( 1 ) ) nsym = nsym + 1 symbol(nsym) = 'HUGE_R' iopsym(nsym) = 0 iprsym(nsym) = 10 numdim(1,nsym) = 1 numdim(2,nsym) = 1 valsym(nsym) = huge ( 1.0E+00 ) nsym = nsym + 1 symbol(nsym) = 'ID' iopsym(nsym) = 1 iprsym(nsym) = 9 numdim(1,nsym) = 0 numdim(2,nsym) = 0 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'INT' iopsym(nsym) = 1 iprsym(nsym) = 9 numdim(1,nsym) = 1 numdim(2,nsym) = 1 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'INV' iopsym(nsym) = 1 iprsym(nsym) = 9 numdim(1,nsym) = 1 numdim(2,nsym) = 1 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'IVAL' iopsym(nsym) = 1 iprsym(nsym) = 9 numdim(1,nsym) = 1 numdim(2,nsym) = 1 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'KRYLOV' iopsym(nsym) = 3 iprsym(nsym) = 9 numdim(1,nsym) = 1 numdim(2,nsym) = 1 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'LCM' iopsym(nsym) = 2 iprsym(nsym) = 9 numdim(1,nsym) = 1 numdim(2,nsym) = 1 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'LN' iopsym(nsym) = 1 iprsym(nsym) = 9 numdim(1,nsym) = 1 numdim(2,nsym) = 1 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'LOG' iopsym(nsym) = 1 iprsym(nsym) = 9 numdim(1,nsym) = 1 numdim(2,nsym) = 1 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'LOG10' iopsym(nsym) = 1 iprsym(nsym) = 9 numdim(1,nsym) = 1 numdim(2,nsym) = 1 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'LOG2' iopsym(nsym) = 1 iprsym(nsym) = 9 numdim(1,nsym) = 1 numdim(2,nsym) = 1 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'MAX' iopsym(nsym) = 2 iprsym(nsym) = 9 numdim(1,nsym) = 1 numdim(2,nsym) = 1 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'MIN' iopsym(nsym) = 2 iprsym(nsym) = 9 numdim(1,nsym) = 1 numdim(2,nsym) = 1 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'NEG' iopsym(nsym) = 1 iprsym(nsym) = 5 numdim(1,nsym) = 1 numdim(2,nsym) = 1 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'NINT' iopsym(nsym) = 1 iprsym(nsym) = 9 numdim(1,nsym) = 1 numdim(2,nsym) = 1 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'NORM0' iopsym(nsym) = 1 iprsym(nsym) = 9 numdim(1,nsym) = 1 numdim(2,nsym) = 1 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'NORM1' iopsym(nsym) = 1 iprsym(nsym) = 9 numdim(1,nsym) = 1 numdim(2,nsym) = 1 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'NORM2' iopsym(nsym) = 1 iprsym(nsym) = 9 numdim(1,nsym) = 1 numdim(2,nsym) = 1 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'NORME' iopsym(nsym) = 1 iprsym(nsym) = 9 numdim(1,nsym) = 1 numdim(2,nsym) = 1 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'NORMF' iopsym(nsym) = 1 iprsym(nsym) = 9 numdim(1,nsym) = 1 numdim(2,nsym) = 1 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'NORMS' iopsym(nsym) = 1 iprsym(nsym) = 9 numdim(1,nsym) = 1 numdim(2,nsym) = 1 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'PI' iopsym(nsym) = 0 iprsym(nsym) = 10 numdim(1,nsym) = 1 numdim(2,nsym) = 1 valsym(nsym) = 4.0E+00 * atan2 ( 1.0, 1.0E+00 ) nsym = nsym + 1 symbol(nsym) = 'POLY' iopsym(nsym) = 2 iprsym(nsym) = 9 numdim(1,nsym) = 1 numdim(2,nsym) = 1 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'RAN' iopsym(nsym) = 1 iprsym(nsym) = 9 numdim(1,nsym) = 0 numdim(2,nsym) = 0 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'RCOND' iopsym(nsym) = 1 iprsym(nsym) = 9 numdim(1,nsym) = 0 numdim(2,nsym) = 0 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'ROUND' iopsym(nsym) = 1 iprsym(nsym) = 9 numdim(1,nsym) = 1 numdim(2,nsym) = 1 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'RTOL' iopsym(nsym) = 0 iprsym(nsym) = 10 numdim(1,nsym) = 1 numdim(2,nsym) = 1 valsym(nsym) = epsilon ( 1.0E+00 ) nsym = nsym + 1 symbol(nsym) = 'RVAL' iopsym(nsym) = 1 iprsym(nsym) = 9 numdim(1,nsym) = 1 numdim(2,nsym) = 1 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'SEC' iopsym(nsym) = 1 iprsym(nsym) = 9 numdim(1,nsym) = 1 numdim(2,nsym) = 1 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'SIN' iopsym(nsym) = 1 iprsym(nsym) = 9 numdim(1,nsym) = 1 numdim(2,nsym) = 1 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'SINH' iopsym(nsym) = 1 iprsym(nsym) = 9 numdim(1,nsym) = 1 numdim(2,nsym) = 1 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'SQRT' iopsym(nsym) = 1 iprsym(nsym) = 9 numdim(1,nsym) = 1 numdim(2,nsym) = 1 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'STEP' iopsym(nsym) = 1 iprsym(nsym) = 9 numdim(1,nsym) = 1 numdim(2,nsym) = 1 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'TAN' iopsym(nsym) = 1 iprsym(nsym) = 9 numdim(1,nsym) = 1 numdim(2,nsym) = 1 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'TANH' iopsym(nsym) = 1 iprsym(nsym) = 9 numdim(1,nsym) = 1 numdim(2,nsym) = 1 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'TRACE' iopsym(nsym) = 1 iprsym(nsym) = 9 numdim(1,nsym) = 1 numdim(2,nsym) = 1 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'TRANS' iopsym(nsym) = 1 iprsym(nsym) = 9 numdim(1,nsym) = 1 numdim(2,nsym) = 1 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'ZEROS' iopsym(nsym) = 1 iprsym(nsym) = 9 numdim(1,nsym) = 0 numdim(2,nsym) = 0 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = '**' iopsym(nsym) = 2 iprsym(nsym) = 8 numdim(1,nsym) = 1 numdim(2,nsym) = 1 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = '^' iopsym(nsym) = 2 iprsym(nsym) = 8 numdim(1,nsym) = 1 numdim(2,nsym) = 1 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = '+' iopsym(nsym) = 2 iprsym(nsym) = 5 numdim(1,nsym) = 1 numdim(2,nsym) = 1 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = '-' iopsym(nsym) = 2 iprsym(nsym) = 5 numdim(1,nsym) = 1 numdim(2,nsym) = 1 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = '/' iopsym(nsym) = 2 iprsym(nsym) = 7 numdim(1,nsym) = 1 numdim(2,nsym) = 1 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = '*' iopsym(nsym) = 2 iprsym(nsym) = 6 numdim(1,nsym) = 1 numdim(2,nsym) = 1 valsym(nsym) = 0.0E+00 ! ! Take care of fact that ninny UNIX FORTRAN compilers won't let us ! explicitly use a backslash character! ! nsym = nsym + 1 symbol(nsym) = char ( 92 ) iopsym(nsym) = 2 iprsym(nsym) = 9 numdim(1,nsym) = 1 numdim(2,nsym) = 1 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = '(' iopsym(nsym) = -1 iprsym(nsym) = 1 numdim(1,nsym) = 1 numdim(2,nsym) = 1 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = ')' iopsym(nsym) = -1 iprsym(nsym) = 2 numdim(1,nsym) = 1 numdim(2,nsym) = 1 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = '=' iopsym(nsym) = 2 iprsym(nsym) = 3 numdim(1,nsym) = 1 numdim(2,nsym) = 1 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = ',' iopsym(nsym) = -1 iprsym(nsym) = 9 numdim(1,nsym) = 1 numdim(2,nsym) = 1 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = '!' iopsym(nsym) = 1 iprsym(nsym) = 9 numdim(1,nsym) = 1 numdim(2,nsym) = 1 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'INDEX1' iopsym(nsym) = 2 iprsym(nsym) = 9 numdim(1,nsym) = 1 numdim(2,nsym) = 1 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'INDEX2' iopsym(nsym) = 3 iprsym(nsym) = 9 numdim(1,nsym) = 1 numdim(2,nsym) = 1 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = '$' iopsym(nsym) = -1 iprsym(nsym) = 0 numdim(1,nsym) = 1 numdim(2,nsym) = 1 valsym(nsym) = 0.0E+00 nsymp = nsym nsyms = nsymp ifree = nsymp + 1 infix = ' ' intsym(1:80) = 0 irpn(1:maxrpn) = 0 istack(1:maxsym) = 0 nints = 0 do i = 1, nsymp ipval(i) = i end do do i = 1, nsymp if ( s_eqi ( symbol(i), 'ANGLE_TO_RADIAN' ) ) then irad = i else if ( s_eqi ( symbol(i), 'NEG' ) ) then ineg = i else if ( symbol(i) == '$' ) then ifinis = i else if ( s_eqi ( symbol(i), 'INDEX1' ) ) then indx1 = i else if ( s_eqi ( symbol(i), 'INDEX2' ) ) then indx2 = i else if ( s_eqi ( symbol(i), 'RTOL' ) ) then irtol = i end if end do call random_seed ( ) value(1:maxrow,1:maxrow) = 0.0E+00 return end subroutine intrea ( ierror, intval, line, prompt ) !*****************************************************************************80 ! !! INTREA gets an integer from the user. ! ! Discussion: ! ! INTREA accepts LINE and a prompt line. If the line is blank, ! the prompt will be printed and line read from the input unit, ! ! In either case, the integer INTVAL will be read from LINE, ! beginning at character 1 and ending at the first comma, slash, ! blank, or the end of line. ! ! The prompt should consist of a string of names of data items, ! separated by commas, with the current one first. INTREA will ! print ! ! 'Enter' prompt ! ! and after reading LINE will strip the characters corresponding ! to INTVAL from LINE, and the part of PROMPT up to the first ! comma, leaving LINE and PROMPT ready for another call to ! CHRREA, INTREA, RATREA or RELREA. ! ! Modified: ! ! 01 May 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer IERROR, error flag. ! 0, no error occurred. ! nonzero, an error occurred. ! ! Output, integer INTVAL, the integer that was read. ! ! Input/output, character ( len = 80 ) LINE, a buffer containing ! unprocessed user input. ! ! Input/output, character ( len = 80 ) PROMPT, the prompt string. ! implicit none character c integer ierror integer intval integer lchar character ( len = 80 ) line integer line_length character ( len = 80 ) prompt intval = 0 ! ! If LINE is empty, fill it up. ! do call chrinp ( ierror, line, prompt ) if ( ierror /= 0 ) then return end if line_length = len_trim ( line ) if ( 0 < line_length ) then exit end if end do ! ! Convert the line into an integer. ! call s_to_i4 ( line, intval, ierror, lchar ) if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'INTREA - Warning!' write ( *, '(a)' ) ' Could not interpret input as an integer.' return end if ! ! Advance LCHAR through blanks and commas. ! line_length = len_trim ( line ) do if ( line_length < lchar + 1 ) then exit end if c = line(lchar+1:lchar+1) if ( c /= ' ' .and. c /= ',' ) then exit end if lchar = lchar + 1 end do ! ! Chop out the portion of LINE that we have used. ! call s_chop ( line, 1, lchar ) return end subroutine krylov ( a, b, na, nb, x ) !*****************************************************************************80 ! !! KRYLOV computes Krylov vectors for a given matrix A and vector X. ! ! Discussion: ! ! In particular, KRYLOV computes a matrix B, containing NB columns, ! where column J is equal to ! ! A**(J-1) * X ! ! Modified: ! ! 18 February 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real A(NA,NA), the multiplying matrix. ! ! Output, real B(NA,NB), the matrix of Krylov vectors. ! ! Input, integer NA, the number of rows and columns in A. ! ! Input, integer NB, the number of columns in B. ! ! Input, real X(NA), the starting vector. ! implicit none integer na integer nb real a(na,na) real b(na,nb) integer j real x(na) ! ! First column of B is X. ! b(1:na,1) = x(1:na) ! ! Column J of B is A * column J-1 of B. ! do j = 2, nb b(1:na,j) = matmul ( a(1:na,1:na), b(1:na,j-1) ) end do return end subroutine matrix_det ( lda, n, a, adet ) !*****************************************************************************80 ! !! MATRIX_DET computes the determinant of a real square matrix. ! ! Modified: ! ! 02 May 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer LDA, the leading dimension of A. ! ! Input, integer N, the number of rows and columns in A. ! ! Input, real A(LDA,N), the matrix whose inverse and ! determinant are desired. ! ! Output, real ADET, the determinant of A. ! implicit none integer lda integer n real a(lda,n) real acopy(n,n) real adet real amax integer j integer k integer kk integer l real s adet = 1.0E+00 ! ! Make a copy of A. ! acopy(1:n,1:n) = a(1:n,1:n) do k = 1, n-1 ! ! Search for pivot index L. ! l = k amax = abs ( acopy(l,l) ) do kk = k+1, n s = abs ( acopy(kk,k) ) if ( amax < s ) then amax = abs ( s ) l = kk end if end do if ( amax == 0.0E+00 ) then adet = 0.0E+00 return end if if ( l /= k ) then adet = - adet call r4row_swap ( n, n, n, acopy, l, k ) end if adet = adet * acopy(k,k) acopy(k+1:n,k) = acopy(k+1:n,k) / acopy(k,k) do j = k+1, n acopy(k:n,j) = acopy(k:n,j) - acopy(k,j) * acopy(k:n,k) end do end do adet = adet * acopy(n,n) return end subroutine matrix_inverse ( lda, n, a, a_inverse, ierror ) !*****************************************************************************80 ! !! MATRIX_INVERSE computes the inverse of a real square matrix. ! ! Modified: ! ! 28 February 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer LDA, the leading dimension of A. ! ! Input, integer N, the number of rows and columns in A. ! ! Input, real A(LDA,N), the matrix whose inverse and ! determinant are desired. ! ! Output, real A_INVERSE(LDA,N), the inverse of A. ! ! Output, integer IERROR, error flag. ! 0, no error occurred. ! nonzero, an error occurred. ! implicit none integer lda integer n real a(lda,n) real a_inverse(lda,n) integer ierror integer pivot(n) ierror = 0 a_inverse(1:lda,1:n) = a(1:lda,1:n) call sge_fa ( a_inverse, lda, n, pivot, ierror ) if ( ierror /= 0 ) then return end if call sge_inv ( lda, n, a_inverse, pivot ) return end subroutine part ( ierror, ifrm, line, maxrow, maxrpn, value ) !*****************************************************************************80 ! !! PART supervises the partitioning of a matrix. ! ! Modified: ! ! 17 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer IERROR, error flag. ! 0, no error occurred. ! nonzero, an error occurred. ! ! Input/output, character ( len = 80 ) LINE, a buffer containing ! unprocessed user input. ! ! Input, integer MAXROW, the maximum number of rows and columns. ! ! Input, integer MAXRPN, specifies the length of IRPN. ! ! Workspace, real VALUE(MAXROW,MAXROW), space used to hold ! the value of the variable. ! implicit none integer, parameter :: maxfix = 80 integer, parameter :: maxlen = 20 integer maxrow integer maxrpn character com integer i integer icol1 integer icola(maxrow) integer ierror integer ifrm character ( len = maxfix ) infix integer ipiece integer irow1 integer irowa(maxrow) integer irpn(maxrpn) character isay integer j character ( len = 80 ) line character ( len = maxlen ) namvar integer ncol integer ncol1 integer nrow integer nrow1 character ( len = 80 ) prompt logical s_eqi real value(maxrow,maxrow) real value1(maxrow,maxrow) ! ! Get the name of the matrix to be partitioned. ! prompt = 'matrix to partition.' call chrrea ( namvar, line, prompt, ierror, 3 ) if ( ierror /= 0 ) then return end if ! ! Get its value. ! nrow = 0 ncol = 0 com = 'R' call comrpn ( com, ierror, ifrm, infix, irpn, maxrow, maxrpn, & namvar, ncol, nrow, value ) if ( ierror /= 0 ) then return end if ! ! Get the rows to keep. ! irowa(1:nrow) = 0 line = ' ' do prompt = 'rows to keep, 0 to end.' call intrea ( ierror, i, line, prompt ) if ( ierror /= 0 ) then return end if if ( i < 1 .or. nrow < i ) then exit end if irowa(i) = 1 end do ! ! Get the columns to keep. ! icola(1:ncol) = 0 line = ' ' do prompt = 'columns to keep, 0 to end.' call intrea ( ierror, i, line, prompt ) if ( ierror /= 0 ) then return end if if ( i < 1 .or. ncol < i ) then exit end if icola(i) = 1 end do ! ! Compute the dimensions of the result. ! nrow1 = sum ( irowa(1:nrow) ) ncol1 = sum ( icola(1:ncol) ) ! ! Copy the submatrix out of the original matrix. ! ipiece = 0 do if ( nrow1 == 0 .or. ncol1 == 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Skipping null submatrix.' else irow1 = 0 do i = 1, nrow if ( 0 < irowa(i) ) then irow1 = irow1 + 1 icol1 = 0 do j = 1, ncol if ( 0 < icola(j) ) then icol1 = icol1 + 1 value1(irow1,icol1) = value(i,j) end if end do end if end do ! ! Get name of matrix in which result is to be stored. ! If matrix already exists, make sure it's big enough. ! If extra-large, zero out extra entries. ! line = ' ' prompt = 'matrix to save submatrix, or return.' call chrrea ( namvar, line, prompt, ierror, 3 ) if ( ierror /= 0 ) then return end if if ( namvar == ' ' ) then namvar = '#' else com = 'A' call comrpn ( com, ierror, ifrm, infix, irpn, maxrow, maxrpn, & namvar, ncol1, nrow1, value ) if ( ierror /= 0 ) then return end if com = 'V' call comrpn ( com, ierror, ifrm, infix, irpn, maxrow, maxrpn, & namvar, ncol1, nrow1, value1 ) if ( ierror /= 0 ) then return end if end if ! ! Print the submatrix. ! com = 'L' call comrpn ( com, ierror, ifrm, infix, irpn, maxrow, maxrpn, & namvar, ncol1, nrow1, value1 ) end if if ( 3 <= ipiece ) then return end if ! ! See if user wants to examine other portions of matrix. ! line = ' ' prompt = '"YES" to see other pieces of matrix.' call chrrea ( isay, line, prompt, ierror, 0 ) if ( .not. s_eqi ( isay, 'Y' ) ) then return end if ipiece = ipiece + 1 if ( ipiece == 1 .or. ipiece == 3 ) then nrow1 = nrow - nrow1 irowa(1:nrow) = 1 - irowa(1:nrow) end if if ( ipiece == 2 ) then ncol1 = ncol - ncol1 icola(1:ncol) = 1 - icola(1:ncol) end if end do return end subroutine poly ( ierror, line, maxrow, maxrpn ) !*****************************************************************************80 ! !! POLY computes a polynomial from its roots, or the characteristic polynomial of a matrix. ! ! Modified: ! ! 17 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer IERROR, error flag. ! 0, no error occurred. ! nonzero, an error occurred. ! ! Input/output, character ( len = 80 ) LINE, a buffer containing ! unprocessed user input. ! ! Input, integer MAXROW, the maximum number of rows and columns. ! ! Input, integer MAXRPN, specifies the length of IRPN. ! implicit none integer, parameter :: maxfix = 80 integer, parameter :: maxlen = 20 integer maxrow integer maxrpn character com integer ierror integer ifrm character ( len = maxfix ) infix integer irpn(maxrpn) character ( len = 80 ) line character ( len = maxlen ) namvar integer ncol integer nmat integer nrow character ( len = 80 ) prompt real pvec(maxrow) real value(maxrow,maxrow) real xvec(maxrow) ! ! Get NAMVAR, the name of the root vector or characteristic ! polynomial matrix. ! prompt = 'vector of roots, or matrix for characteristic polynomial.' call chrrea ( namvar, line, prompt, ierror, 3 ) if ( ierror /= 0 ) then return end if ! ! See if we know about this variable. ! ifrm = 0 nrow = 0 ncol = 0 com = 'R' call comrpn ( com, ierror, ifrm, infix, irpn, maxrow, maxrpn, & namvar, ncol, nrow, value ) if ( ierror /= 0 ) then return end if ! ! Make sure the variable is a vector or a square matrix. ! if ( nrow /= ncol .and. nrow /= 1 .and. ncol /= 1 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'POLY - Error!' write ( *, '(a)' ) ' You must enter a vector, or a square matrix!' ierror = 1 return end if nmat = max ( nrow, ncol ) ! ! If it's a vector, then we assume it's a vector of roots. ! if ( nrow == 1 .or. ncol == 1 ) then if ( nrow == 1 ) then xvec(1:nmat) = value(1,1:nmat) else xvec(1:nmat) = value(1:nmat,1) end if call roots_to_poly ( nmat, xvec, pvec ) ! ! Otherwise, it's a matrix and we want the characteristic polynomial. ! else call r4mat_poly_char ( maxrow, nmat, value, pvec ) end if nrow = nmat + 1 ncol = 1 ! ! Print the polynomial. ! call r4poly_print ( nmat, pvec, ' The polynomial:' ) ! ! Offer to save it somewhere. ! prompt = 'variable to save polynomial in, or RETURN to discard.' line = ' ' call chrrea ( namvar, line, prompt, ierror, 3 ) if ( ierror /= 0 ) then return end if if ( namvar == ' ' ) then return end if ! ! Add the name. ! com = 'A' call comrpn ( com, ierror, ifrm, infix, irpn, maxrow, maxrpn, & namvar, ncol, nrow, value ) if ( ierror /= 0 ) then return end if ! ! Add the value. ! value(1:nrow,1) = pvec(1:nrow) com = 'V' call comrpn ( com, ierror, ifrm, infix, irpn, maxrow, maxrpn, & namvar, ncol, nrow, value ) if ( ierror /= 0 ) then return end if write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'POLY:' write ( *, '(a)' ) ' Polynomial coefficients saved in ' // trim ( namvar ) return end function pythag ( a, b ) !*****************************************************************************80 ! !! PYTHAG finds SQRT(A**2+B**2) without overflow or destructive underflow. ! ! Modified: ! ! 15 March 2001 ! ! Parameters: ! ! Input, real A, B, the two numbers whose "hypotenuse" is ! to be computed. ! ! Output, real PYTHAG, the computed value of SQRT(A**2+B**2). ! implicit none real a real b real p real pythag real r real s real t real u p = max ( abs ( a ), abs ( b ) ) if ( p == 0.0E+00 ) then pythag = p return end if r = ( min ( abs ( a ), abs ( b ) ) / p )**2 do t = 4.0E+00 + r if ( t == 4.0E+00 ) then pythag = p exit end if s = r / t u = 1.0E+00 + 2.0E+00 * s p = u * p r = ( s / u )**2 * r end do return end subroutine qrfact ( a, m, lda, n, q, r ) !*****************************************************************************80 ! !! QRFACT computes the QR factorization of a matrix. ! ! Discussion: ! ! The M by N matrix A is factored as Q * R, where Q is an ! M by M orthogonal matrix, and R is an M by N upper triangular ! matrix. ! ! Modified: ! ! 17 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real A(LDA,N), the M by N matrix to be factored. ! ! Input, integer M, the number of rows in A. ! ! Input, integer LDA, the leading dimension of A. ! ! Input, integer N, the number of columns of A. ! ! Output, real Q(LDA,M), the M by M orthogonal matrix. ! ! Output, real R(LDA,N), the M by N upper triangular matrix. ! implicit none integer lda integer m integer n real a(lda,n) real c integer i integer ii integer j integer jj integer k real q(lda,m) real qi real qj real r(lda,n) real rad real ri real rj real s real scale ! ! Copy A into R. ! r = a ! ! Initialize Q to the identity matrix. ! call r4mat_identity ( lda, m, q ) do j = 1, min ( n, m - 1 ) do i = j + 1, m if ( r(i,j) /= 0.0E+00 ) then scale = abs ( r(i,j) ) + abs ( r (j,j) ) ri = r(i,j) / scale rj = r(j,j) / scale rad = sqrt ( rj**2 + ri**2 ) if ( r(i,j) < 0.0E+00 ) then rad = - rad end if c = rj / rad s = ri / rad do jj = 1, n ri = r(i,jj) rj = r(j,jj) r(j,jj) = c * rj + s * ri r(i,jj) = - s * rj + c * ri end do r(i,j) = 0.0E+00 do ii = 1, m qi = q(ii,i) qj = q(ii,j) q(ii,i) = c * qi - s * qj q(ii,j) = s * qi + c * qj end do end if end do end do ! ! Ensure that the diagonal of R is nonnegative. ! do k = 1, min ( m, n ) if ( r(k,k) < 0.0E+00 ) then do j = 1, n r(k,j) = - r(k,j) end do do i = 1, m q(i,k) = - q(i,k) end do end if end do return end subroutine r4_swap ( x, y ) !*****************************************************************************80 ! !! R4_SWAP switches two R4's. ! ! Modified: ! ! 30 November 1998 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, real X, Y. On output, the values of X and ! Y have been interchanged. ! implicit none real x real y real z z = x x = y y = z return end subroutine r4col_swap ( lda, m, n, a, i, j ) !*****************************************************************************80 ! !! R4COL_SWAP swaps columns I and J of an R4COL. ! ! Example: ! ! Input: ! ! M = 3, N = 4, I = 2, J = 4 ! ! A = ( ! 1. 2. 3. 4. ! 5. 6. 7. 8. ! 9. 10. 11. 12. ) ! ! Output: ! ! A = ( ! 1. 4. 3. 2. ! 5. 8. 7. 6. ! 9. 12. 11. 10. ) ! ! Modified: ! ! 22 July 1998 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer LDA, the leading dimension of the array. ! ! Input, integer M, N, the number of rows and columns. ! ! Input, real A(LDA,N), the M by N array. ! ! Input, integer I, J, the columns to be swapped. ! implicit none integer lda integer n real a(lda,n) integer i integer j integer k integer m if ( 1 <= i .and. i <= n .and. 1 <= j .and. j <= n ) then do k = 1, m call r4_swap ( a(k,i), a(k,j) ) end do else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'R4COL_SWAP - Fatal error!' write ( *, '(a)' ) ' I or J is out of bounds.' write ( *, '(a,i8)' ) ' I = ', i write ( *, '(a,i8)' ) ' J = ', j write ( *, '(a,i8)' ) ' NCOL = ', n stop end if return end subroutine relrea ( rval, line, prompt, ierror ) !*****************************************************************************80 ! !! RELREA extracts a real number from a string of input. ! ! Discussion: ! ! RELREA accepts a line of characters which may contain some ! user input. If not, it prints out the prompt and reads new ! information into line, seeking to find a real number RVAL to ! return. ! ! RELREA will accept integers, real numbers, and ratios of the ! form R1/R2. Real numbers may be in scientific notation, as ! +12.34E-56.78 ! ! Modified: ! ! 15 March 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, real RVAL, the real value that was read in. ! ! Input/output, character ( len = 80 ) LINE, a buffer containing ! unprocessed user input. ! ! Input/output, character ( len = 80 ) PROMPT, used for the prompt string. ! ! Output, integer IERROR, error flag. ! 0, no error occurred. ! nonzero, an error occurred. ! implicit none real bot integer ierror integer lchar character ( len = 80 ) line character ( len = 80 ) prompt real rval ! ! Set the default value. ! rval = 0.0E+00 ! ! Get the input LINE, but if it's blank, demand real input. ! do call chrinp ( ierror, line, prompt ) if ( ierror /= 0 ) then return end if if ( 0 < len_trim ( line ) ) then exit end if end do ! ! Convert the string to a real value. ! call s_to_r4 ( line, rval, ierror, lchar ) ! ! If there's leftover stuff, check for a division sign. ! if ( lchar+1 < len_trim ( line ) ) then if ( line(lchar+1:lchar+1) == '/' ) then lchar = lchar + 1 call s_chop ( line, 1, lchar ) call s_to_r4 ( line, bot, ierror, lchar ) if ( bot == 0.0E+00 ) then bot = 1.0E+00 end if rval = rval / bot end if end if ! ! Remove the characters you used from the input line. ! call s_chop ( line, 1, lchar ) return end subroutine rg ( lda, n, a, wr, wi, matz, z, ierror ) !*****************************************************************************80 ! !! RG finds the eigenvalues and eigenvectors of a real general matrix. ! ! Modified: ! ! 01 May 2000 ! ! Parameters: ! ! Input, integer LDA, the leading dimension of A and Z. ! LDA must be at least N. ! ! Input, integer N, the number of rows and columns of A. ! ! Input/output, real A(LDA,N). ! ! On input, A contains the N by N matrix whose eigenvalues and ! eigenvectors are desired. ! ! On output, A has been overwritten by other information. ! ! Output, real WR(N), WI(N), contain the real and imaginary parts, ! respectively, of the eigenvalues. Complex conjugate ! pairs of eigenvalues appear consecutively with the ! eigenvalue having the positive imaginary part first. ! ! Input, integer MATZ, zero if only eigenvalues are desired. ! non-zero, for both eigenvalues and eigenvectors. ! ! Output, real Z(LDA,N), contains the real and imaginary parts of ! the eigenvectors if MATZ is not zero. If the J-th eigenvalue ! is real, the J-th column of Z contains its eigenvector. If the ! J-th eigenvalue is complex with positive imaginary part, the ! J-th and (J+1)-th columns of Z contain the real and ! imaginary parts of its eigenvector. The conjugate of this ! vector is the eigenvector for the conjugate eigenvalue. ! ! Output, integer IERROR, error flag. ! 0, no error. ! nonzero, an error occurred. ! implicit none integer lda integer n real a(lda,n) real fv1(n) integer ierror integer is1 integer is2 integer iv1(n) integer matz real wi(n) real wr(n) real z(lda,n) ierror = 0 if ( lda < n ) then ierror = 10 * n return end if ! ! Balance the matrix. ! call balanc ( lda, n, a, is1, is2, fv1 ) ! ! Put the matrix into upper Hessenberg form. ! call elmhes ( lda, n, is1, is2, a, iv1 ) if ( matz == 0 ) then call hqr ( lda, n, is1, is2, a, wr, wi, ierror ) if ( ierror /= 0 ) then return end if else call eltran ( lda, n, is1, is2, a, iv1, z ) call hqr2 ( lda, n, is1, is2, a, wr, wi, z, ierror ) if ( ierror /= 0 ) then return end if call balbak ( lda, n, is1, is2, fv1, n, z ) end if return end subroutine r4mat_identity ( lda, n, a ) !*****************************************************************************80 ! !! R4MAT_IDENTITY sets the square matrix A to the identity. ! ! Modified: ! ! 24 March 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer LDA, the leading dimension of A. ! ! Input, integer N, the order of A. ! ! Output, real A(LDA,N), the matrix which has been ! set to the identity. ! implicit none integer lda integer n real a(lda,n) integer i integer j do i = 1, n do j = 1, n if ( i == j ) then a(i,j) = 1.0E+00 else a(i,j) = 0.0E+00 end if end do end do return end subroutine r4mat_norm2 ( a, nrow, tnorm ) !*****************************************************************************80 ! !! R4MAT_NORM2 computes the two-norm of a square matrix A. ! ! Discussion: ! ! The two-norm can be computed as the square root of the spectral ! radius of Transpose(A)*A. ! ! Modified: ! ! 01 May 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real A(NROW,NROW), the matrix whose two-norm is desired. ! ! Input, integer NROW, the number of rows and columns in A. ! ! Output, real TNORM, the two-norm of A. ! implicit none integer nrow real a(nrow,nrow) real b(nrow,nrow) real evali(nrow) real evalr(nrow) real evvect(1,1) integer i integer ierror integer j integer matz real tnorm ! ! Compute A' * A. ! do i = 1, nrow do j = 1, nrow b(i,j) = dot_product ( a(1:nrow,i), a(1:nrow,j) ) end do end do ! ! Compute the eigenvalues of A' * A. ! matz = 0 call eigen ( b, evali, evalr, evvect, ierror, matz, nrow ) if ( ierror /= 0 ) then tnorm = -1.0E+00 return end if ! ! Compute the maximum eigenvalue norm. ! tnorm = sqrt ( maxval ( evalr(1:nrow)**2 + evali(1:nrow)**2 ) ) return end function r4mat_norme ( a, lda, m, n ) !*****************************************************************************80 ! !! R4MAT_NORME returns the EISPACK norm of an M by N matrix. ! ! Definition: ! ! The EISPACK norm is defined as: ! ! norm = Sum i = 1 to M ( Sum j = 1 to N ( Abs ( A(I,J) ) ) ) ! ! Modified: ! ! 01 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real A(LDA,N), the matrix whose EISPACK norm is desired. ! ! Input, integer LDA, the leading dimension of A. ! ! Input, integer M, the number of rows in A. ! ! Input, integer N, the number of columns in A. ! ! Output, real R4MAT_NORME, the EISPACK norm of A. ! implicit none integer lda integer n real a(lda,n) integer m real r4mat_norme r4mat_norme = sum ( abs ( a(1:m,1:n) ) ) return end function r4mat_normf ( a, lda, m, n ) !*****************************************************************************80 ! !! R4MAT_NORMF returns the Frobenius norm of an M by N matrix. ! ! Definition: ! ! The Frobenius norm is defined as ! ! norm = Sqrt ( Sum i = 1 to M ( Sum j = 1 to N ( A(I,J)**2) ) ) ! ! Note: ! ! The matrix Frobenius-norm is not derived from a vector norm, but ! is compatible with the vector 2-norm, so that: ! ! Norm2 ( A*x ) <= Normf ( A ) * Norm2 ( x ). ! ! Modified: ! ! 01 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real A(LDA,N), the matrix whose Frobenius norm is desired. ! ! Input, integer LDA, the leading dimension of A. ! ! Input, integer M, the number of rows in A. ! ! Input, integer N, the number of columns in A. ! ! Output, real R4MAT_NORMF, the Frobenius norm of A. ! implicit none integer lda integer n real a(lda,n) integer m real r4mat_normf r4mat_normf = sqrt ( sum ( a(1:m,1:n)**2 ) ) return end subroutine r4mat_poly_char ( lda, n, a, p ) !*****************************************************************************80 ! !! R4MAT_POLY_CHAR computes the characteristic polynomial of a matrix. ! ! Modified: ! ! 15 March 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer LDA, the leading dimension of the matrix A. ! ! Input, integer N, the order of the matrix A. ! ! Input, real A(LDA,N), the N by N matrix. ! ! Output, real P(0:N), the coefficients of the characteristic ! polynomial of A. P(N) contains the coefficient of X**N ! (which will be 1), P(I) contains the coefficient of X**I, ! and P(0) contains the constant term. ! implicit none integer lda integer n real a(lda,n) integer i integer iorder real p(0:n) real r4mat_trace real trace real work1(n,n) real work2(n,n) ! ! Initialize WORK1 to the identity matrix. ! call r4mat_identity ( n, n, work1 ) p(n) = 1.0E+00 do iorder = n-1, 0, -1 ! ! Work2 = A * WORK1. ! work2(1:n,1:n) = matmul ( a(1:n,1:n), work1(1:n,1:n) ) ! ! Take the trace. ! trace = r4mat_trace ( n, n, work2 ) ! ! P(IORDER) = - Trace ( WORK2 ) / ( N - IORDER ) ! p(iorder) = - trace / real ( n - iorder ) ! ! WORK1 := WORK2 + P(IORDER) * Identity. ! work1(1:n,1:n) = work2(1:n,1:n) do i = 1, n work1(i,i) = work1(i,i) + p(iorder) end do end do return end subroutine r4mat_poly_val ( a, work, result, n, coef, npoly ) !*****************************************************************************80 ! !! R4MAT_POLY_VAL evaluates a polynomial with a square matrix as its argument. ! ! Modified: ! ! 17 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real A(N,N), the matrix which is the argument of the polynomial. ! ! Workspace, real WORK(N,N). ! ! Output, real RESULT(N,N), the value of the polynomial. ! ! Input, integer N, the order of A. ! ! Input, real COEF(0:NPOLY), the polynomial coefficients. ! COEF(NPOLY) is the coefficient of A**NPOLY ! COEF(I) is the coefficient of A**(I) ! COEF(0) is the constant term. ! ! Input, integer NPOLY, the order of the polynomial. ! implicit none integer n integer npoly real a(n,n) real coef(0:npoly) integer i integer j integer k integer l real result(n,n) real work(n,n) ! ! Set Result = CN * I ! do i = 1, n do j = 1, n if ( i == j ) then result(i,j) = coef(npoly) else result(i,j) = 0.0E+00 end if end do end do ! ! Set Result = (Result*A+Ck*I) ! do l = npoly-1, 0, -1 do i = 1, n do j = 1, n work(i,j) = 0.0E+00 do k = 1, n work(i,j) = work(i,j) + result(i,k) * a(k,j) end do end do end do result(1:n,1:n) = work(1:n,1:n) do i = 1, n result(i,i) = result(i,i) + coef(l) end do end do return end subroutine r4mat_power ( maxrow, n, a, npow, b ) !*****************************************************************************80 ! !! R4MAT_POWER computes A**NPOW, the nonnegative power of a real square matrix. ! ! Discussion: ! ! The algorithm is: ! ! B=I ! Do NPOW times: ! B=A*B ! End do ! ! Modified: ! ! 14 February 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer MAXROW, the leading dimension of A and B. ! ! Input, integer N, the order of A. ! ! Input, real A(MAXROW,N), the matrix to be raised to a power. ! ! Input, integer NPOW, the power to which A is to be raised. ! ! Output, real B(MAXROW,N), the value of A**NPOW. ! implicit none integer maxrow integer n real a(maxrow,n) real b(maxrow,n) integer ipow integer npow call r4mat_identity ( maxrow, n, b ) do ipow = 1, npow b(1:n,1:n) = matmul ( a(1:n,1:n), b(1:n,1:n) ) end do return end subroutine r4mat_print ( maxrow, nrow, ncol, value, namvr ) !*****************************************************************************80 ! !! R4MAT_PRINT prints out a scalar, vector or array. ! ! Discussion: ! ! It tries to print out integer valued reals in compact format. ! It transposes column vectors. It can print up to 5 columns ! of an array or vector per line. ! ! For an array, it will label the columns. ! ! Modified: ! ! 26 October 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer MAXROW, the maximum number of rows and columns. ! ! Input, character ( len = MAXLEN ) NAMVR, the name of the item to be ! printed. ! ! Input, integer NCOL, NROW, the number of columns and rows ! in the item to be printed. ! ! Input, real VALUE(MAXROW,MAXROW), space used to hold ! the value of the variable to be saved. ! implicit none integer, parameter :: maxlen = 20 ! ! You don't want to know why I had to use two billion! ! real, parameter :: big = 2000000000.0E+00 integer maxrow integer ncol character ( len = 20 ) chrtmp integer i integer i1 integer i2 integer ihi integer ii integer ilo integer inc integer iv integer j integer jhi integer jinc integer jlo character ( len = maxlen ) namvr integer ndim integer nrow character ( len = 100 ) output real v real value(maxrow, ncol) if ( nrow <= 0 ) then return end if if ( ncol <= 0 ) then return end if ! ! Print a scalar value. ! if ( nrow == 1 .and. ncol == 1 ) then v = value(1,1) if ( big < abs ( v ) ) then write ( *, '(a,a,g16.8)' ) trim ( namvr ), ' = ', v else if ( real ( int ( v ) ) == v ) then write ( *, '(a,a,i12)' ) trim ( namvr ), ' = ', int ( v ) else write ( *, '(a,a,g16.8)' ) trim ( namvr ), ' = ', v end if ! ! Print a row or column vector. ! else if ( nrow == 1 .or. ncol == 1 ) then write ( *, '(a)' ) ' ' if ( nrow == 1 ) then output = namvr else output = namvr // ' Transposed' call s_blanks_delete ( output ) end if write ( *, '(a)' ) trim ( output ) write ( *, '(a)' ) ' ' ndim = nrow * ncol inc = 4 do ilo = 1, ndim, inc ihi = min ( ilo+inc-1, ndim ) ii = 0 output = ' ' do i = ilo, ihi ii = ii + 1 if ( nrow == 1 ) then v = value(1,i) else v = value(i,1) end if write ( chrtmp, '(g20.8)' ) v if ( abs ( v ) <= big ) then iv = int ( v ) if ( real ( iv ) == v .and. abs ( iv ) <= 999999 ) then write ( chrtmp, '(i7)' ) iv end if end if i1 = ( ii - 1 ) * 20 + 1 i2 = i1 + 20 - 1 output(i1:i2) = chrtmp end do write ( *, '(a)' ) trim ( output ) end do write ( *, '(a)' ) ' ' ! ! Print an array. ! ! If the matrix is too "wide" to print all at once, declare which ! columns are being printed. ! else write ( *, '(a)' ) ' ' write ( *, '(a)' ) trim ( namvr ) write ( *, '(a)' ) ' ' jinc = 4 do jlo = 1, ncol, jinc jhi = min ( jlo+jinc-1, ncol ) if ( jlo /= 1 .or. jhi /= ncol ) then write ( *, '(a)' ) ' ' if ( jhi /= jlo ) then write ( *, '(a,i8,a,i8)' ) 'Columns ', jlo, ' to ', jhi else write ( *, '(a,i8)' ) 'Column ', jlo end if end if do i = 1, nrow ii = 0 output = ' ' do j = jlo, jhi ii = ii + 1 v = value(i,j) write ( chrtmp, '(g20.8)' ) v if ( abs ( v ) <= big ) then iv = int ( v ) if ( v == iv .and. abs ( iv ) <= 999999 ) then write ( chrtmp, '(i7)' ) iv end if end if i1 = ( ii - 1 ) * 20 + 1 i2 = i1 + 20 - 1 output(i1:i2) = chrtmp end do write ( *, '(a)' ) trim ( output ) end do end do end if return end function r4mat_trace ( lda, n, a ) !*****************************************************************************80 ! !! R4MAT_TRACE computes the trace of a real matrix. ! ! Definition: ! ! The trace of a square matrix is the sum of the diagonal elements. ! ! Modified: ! ! 20 July 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer LDA, the leading dimension of the matrix A. ! ! Input, integer N, the order of the matrix A. ! ! Input, real A(LDA,N), the N by N matrix whose trace is desired. ! ! Output, real R4MAT_TRACE, the trace of the matrix. ! implicit none integer lda integer n real a(lda,n) integer i real r4mat_trace r4mat_trace = 0.0E+00 do i = 1, n r4mat_trace = r4mat_trace + a(i,i) end do return end subroutine roots_to_poly ( n, x, c ) !*****************************************************************************80 ! !! ROOTS_TO_POLY converts polynomial roots to polynomial coefficients. ! ! Example: ! ! N = 3 ! X = (/ 1, 2, 3 /) ! ! P(X) = ( X - 1 ) * ( X - 2 ) * ( X - 3 ) ! = X**3 - 6 * X**2 + 11 * X - 6 ! ! C = (/ -6, 11, -6, 1 /) ! ! Modified: ! ! 17 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer N, the number of roots specified. ! ! Input, real X(N), the roots. ! ! Output, real C(0:N), the coefficients of the polynomial. ! implicit none integer n real c(0:n) integer i integer j real x(n) ! ! Initialize C to (0, 0, ..., 0, 1). ! Essentially, we are setting up a divided difference table. ! c(0:n-1) = 0.0E+00 c(n) = 1.0E+00 ! ! Convert to standard polynomial form by shifting the abscissas ! of the divided difference table to 0. ! do j = 1, n do i = 1, n+1-j c(n-i) = c(n-i) - x(n+1-i-j+1) * c(n-i+1) end do end do return end subroutine round ( a, lda, m, n, rtol ) !*****************************************************************************80 ! !! ROUND rounds the entries of the rectangular array A. ! ! Modified: ! ! 15 March 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, real A(LDA,N), the M by N array containing ! entries to be rounded. ! ! Input, integer LDA, the leading dimension of A, which must ! be at least M. ! ! Input, integer M, N, the number of rows and columns in A. ! ! Input, real RTOL, the rounding tolerance. ! implicit none integer lda integer n real a(lda,n) integer i integer j integer m real rtol do i = 1, m do j = 1, n if ( abs ( a(i,j) ) < rtol ) then a(i,j) = 0.0E+00 end if end do end do return end subroutine rpnchk ( ierror, ihi, ilo, iopsym, irpn, maxrpn, maxsym ) !*****************************************************************************80 ! !! RPNCHK examines an RPN formula, looking for a complete RPN expression. ! ! Discussion: ! ! The routine starts at location IHI, and finds the position ILO ! such that IRPN(ILO)...IRPN(IHI) represents a single argument. ! ! Modified: ! ! 24 October 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer IERROR, 0, no error; 1 an error. ! ! Input, integer IHI, the location in IRPN where the search begins. ! ! Output, integer ILO, the location in IRPN such that IRPN(ILO) through ! IRPN(IHI) represents a single argument, or IHI+1 if no such location ! was found. ! ! Input, integer IOPSYM(MAXSYM), the number of operands per symbol. ! Constants have 0, unaries have 1, and so on. ! ! Workspace, integer IRPN(MAXRPN), used to store the compiled ! versions of user formulas. ! ! Input, integer MAXRPN, specifies the length of IRPN. ! ! Input, integer MAXSYM, the maximum number of symbols allowed. ! implicit none integer maxrpn integer maxsym integer ierror integer ihi integer ilo integer iopsym(maxsym) integer irpn(maxrpn) integer isum isum = 0 ierror = 0 ilo = ihi + 1 do ilo = ilo - 1 if ( ilo <= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'RPNCHK - Fatal error!' write ( *, '(a)' ) ' Reached beginning of RPN formula,' write ( *, '(a)' ) ' but still had operators missing arguments.' ierror = 1 return end if if ( iopsym(irpn(ilo)) < 0 .and. isum == 0 ) then cycle end if isum = isum + 1 - iopsym(irpn(ilo)) if ( isum == 1 ) then exit end if end do return end subroutine rpnset ( ierror, ifinis, imin, intsym, iopsym, iprsym, irpn, & istack, maxrpn, maxsym, nints, nrpn, symbol ) !*****************************************************************************80 ! !! RPNSET converts the infix formula into an RPN formula. ! ! Modified: ! ! 25 October 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer IERROR, error flag. ! 0, no error occurred. ! nonzero, an error occurred. ! ! Input, integer IFINIS, the index in SYMBOL of the symbol ! '$', meaning the end of the formula. ! ! Input, integer IMIN, an offset for accessing IRPN. ! ! Input, integer INTSYM(80), a set of integers which ! are the indices of symbols, representing an infix formula. ! ! Input, integer IOPSYM(MAXSYM), the number of operands per symbol. ! Constants have 0, unaries have 1, and so on. ! ! Input, integer IPRSYM(MAXSYM), the priority of each symbol. ! ! Output, integer IRPN(MAXRPN), the RPN version of the infix formula. ! ! Workspace, integer ISTACK(MAXSYM). ! ! Input, integer MAXRPN, the maximum number of RPN symbols allowed. ! ! Input, integer MAXSYM, the maximum number of symbols allowed. ! ! Input, integer NINTS, the number of useful entries in INTSYM. ! ! Output, integer NRPN, the number of useful entries in IRPN. ! ! Input, character ( len = MAXLEN ) SYMBOL(MAXSYM), the symbolic names ! of the internally defined constants and functions. ! implicit none integer, parameter :: maxlen = 20 integer maxrpn integer maxsym integer ierror integer ifinis integer imin integer intsym(80) integer iopsym(maxsym) integer iprsym(maxsym) integer iread integer irpn(maxrpn) integer istack(maxsym) integer istak integer isym integer jsym integer nints integer nrpn character ( len = maxlen ) sym character ( len = maxlen ) sym2 character ( len = maxlen ) symbol(maxsym) ierror = 0 nrpn = 0 ! ! An error if the formula seems to have nothing in it. ! if ( nints <= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'RPNSET - Error!' write ( *, '(a)' ) ' This formula seems to be blank!' ierror = 1 return end if iread = 0 istak = 0 ! ! Read symbol number IREAD from INTSYM. ! do iread = iread + 1 if ( nints < iread ) then ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'RPNSET - Error!' write ( *, '(a)' ) ' This formula does not make sense!' return end if isym = intsym(iread) if ( isym == 0 ) then cycle end if sym = symbol(isym) if ( sym == ',' ) then cycle end if ! ! If the symbol is "$", then it's time to pop the stack. ! if ( sym == '$' ) then do if ( istak <= 0 ) then nrpn = nrpn + 1 irpn(nrpn+imin) = ifinis if ( iread < nints ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'RPNSET - Error!' write ( *, '(a)' ) ' Some of the formula is left over!' ierror = 1 end if return end if nrpn = nrpn + 1 irpn(nrpn+imin) = istack(istak) istak = istak - 1 end do end if ! ! A left parenthesis goes onto the stack. ! if ( sym == '(' ) then istak = istak + 1 istack(istak) = isym cycle end if ! ! A variable or constant goes immediately into IRPN. ! if ( iopsym(isym) == 0 ) then nrpn = nrpn + 1 irpn(nrpn+imin) = isym cycle end if ! ! Put operators onto the stack. ! do if ( istak <= 0 ) then istak = istak + 1 istack(istak) = isym exit end if jsym = istack(istak) if ( iprsym(jsym) <= iprsym(isym) ) then jsym = istack(istak) sym2 = symbol(jsym) if ( sym == ')' .and. sym2 == '(' ) then istak = istak - 1 else istak = istak + 1 istack(istak) = isym end if exit end if nrpn = nrpn + 1 irpn(nrpn+imin) = istack(istak) istak = istak - 1 end do end do return end subroutine rpnval ( ierror, ifree, imin, iopsym, iprsym, ipval, & irad, irpn, irtol, istack, maxrow, maxrpn, maxsym, maxval, ncol, nrow, & nrpn, nsym, nsyms, numdim, symbol, valsym, value ) !*****************************************************************************80 ! !! RPNVAL evaluates the symbolic functions in an RPN formula. ! ! Discussion: ! ! The routine determines the number of arguments, gathering their values, ! and calling FUNVAL to evaluate the given function. ! ! Modified: ! ! 01 November 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer IERROR. ! 0, no error occurred. ! nonzero, an error occurred. ! ! Input, integer IFREE, the next free location in VALSYM. ! ! Input, integer IMIN, an offset for indexing IRPN. ! ! Input, integer IOPSYM(MAXSYM), the number of operands per symbol. ! Constants have 0, unaries have 1, and so on. ! ! Input, integer IPRSYM(MAXSYM), the priority of each symbol. ! ! Input, integer IPVAL(MAXSYM), contains, for each symbol, ! the address in VALSYM where the value of the symbol ! is stored, if it is a scalar. If the symbol represents ! a vector or matrix, IPVAL(IARG) points to location of ! the first entry. ! ! Workspace, integer IRPN(MAXRPN), used to store the compiled ! versions of user formulas. ! ! Workspace, integer ISTACK(MAXSYM), workspace for interpreting the formula. ! ! Input, integer MAXROW, the maximum number of rows and columns. ! ! Input, integer MAXRPN, specifies the length of IRPN. ! ! Input, integer MAXSYM, the maximum number of symbols allowed. ! ! Input, integer MAXVAL, the maximum number of values allowed. ! ! Output, integer NCOL, the number of columns in the result. ! ! Output, integer NROW, the number of rows in the result. ! ! Input, integer NRPN, the number of useful entries in IRPN. ! ! Input, integer NSYM, the total number of symbols, including temporaries. ! ! Input, integer NSYMS, the number of declared symbols. ! ! Input, integer NUMDIM(2,MAXSYM). ! For each symbol I, NUMDIM(1,I) is the number of rows in its ! value, and NUMDIM(2,I) is the number of columns. ! ! Input, character ( len = MAXLEN ) SYMBOL(MAXSYM), the symbolic names ! of the internally defined constants and functions. ! ! Input, real VALSYM(MAXVAL), the values of the symbolic variables. ! ! Workspace, real VALUE(MAXROW,MAXROW), the value of the variable. ! implicit none integer, parameter :: maxlen = 20 integer maxrow integer maxrpn integer maxsym character ( len = 3 ) ctemp integer iarg1 integer iarg2 integer iarg3 integer ierror integer ifree integer imin integer index1 integer index4 integer iopsym(maxsym) integer iprsym(maxsym) integer ipset integer ipval(maxsym) integer irad integer iread integer irpn(maxrpn) integer irtol integer istack(maxsym) integer istak integer isym integer isymo integer itemp integer j integer maxval integer ncol integer nrow integer nrpn integer nsym integer nsyms integer numdim(2,maxsym) character ( len = maxlen ) sym character ( len = maxlen ) symbol(maxsym) real valsym(maxval) real value(maxrow,maxrow) nrow = 1 ncol = 1 value(1,1) = 0.0E+00 if ( nrpn <= 1 ) then return end if istak = 0 isymo = 0 nsyms = nsym itemp = 0 iread = 0 do iread = iread + 1 if ( nrpn < iread ) then ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'RPNVAL - Error!' write ( *, '(a)' ) ' The formula ended unexpectedly!' exit end if isym = irpn(iread+imin) if ( 0 < isym .and. isym <= maxsym ) then sym = symbol(isym) else sym = '$' end if if ( sym == '$' .or. nrpn < iread ) then nrow = numdim(1,isymo) ncol = numdim(2,isymo) index4 = ipval(isymo) do j = 1, ncol value(1:nrow,j) = valsym(index4:index4+nrow-1) index4 = index4 + nrow end do exit end if ! ! Constants and variables go into stack. ! if ( iopsym(isym) == 0 ) then if ( isym <= nsym ) then istak = istak+1 istack(istak) = isym isymo = isym else if ( maxsym <= nsyms ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'RPNVAL - Error!' write ( *, '(a)' ) ' There is no memory left for more symbols!' write ( *, '(a)' ) ' Perhaps the KILL or INIT command would help.' ierror = 1 exit end if nsyms = nsyms+1 isymo = nsyms istak = istak+1 istack(istak) = nsyms iprsym(nsyms) = 10 iopsym(nsyms) = 0 itemp = itemp+1 ctemp = ' ' call i4_to_s_zero ( itemp, ctemp ) symbol(nsyms) = 'STK000' symbol(nsyms)(4:6) = ctemp nrow = numdim(1,isym) ncol = numdim(2,isym) numdim(1,nsyms) = nrow numdim(2,nsyms) = ncol ipval(nsyms) = ifree ifree = ifree + numdim(1,nsyms) * numdim(2,nsyms) index1 = ipval(isym) index4 = ipval(nsyms) valsym(index4:index4+nrow*ncol-1) = & valsym(index1:index1+nrow*ncol-1) end if cycle end if ! ! Pull off arguments. ! iarg1 = istack(istak) iarg2 = 0 iarg3 = 0 if ( iopsym(isym) == 2 ) then iarg2 = istack(istak) istak = istak-1 iarg1 = istack(istak) else if (iopsym(isym) == 3 ) then iarg3 = istack(istak) istak = istak-1 iarg2 = istack(istak) istak = istak-1 iarg1 = istack(istak) end if sym = symbol(isym) call funval ( iarg1, iarg2, iarg3, ierror, ifree, iopsym, iprsym, & ipset, ipval, irad, irtol, itemp, maxrow, maxsym, maxval, nsyms, & numdim, sym, symbol, valsym, value ) isymo = nsyms if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'RPNVAL - Warning!' write ( *, '(a)' ) ' Evaluation of this formula is abandoned.' nrow = 1 ncol = 1 value(1,1) = 0.0E+00 exit end if if ( ipset == 0 ) then ipval(nsyms) = ifree ifree = ifree + numdim(1,nsyms) * numdim(2,nsyms) else ipval(nsyms) = ipset ipset = 0 end if istack(istak) = nsyms end do return end subroutine r4poly_print ( n, poly_cof, title ) !*****************************************************************************80 ! !! R4POLY_PRINT prints out a polynomial. ! ! Modified: ! ! 09 February 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer N, the degree of the polynomial. ! ! Input, real POLY_COF(0:N), the polynomial coefficients. ! POLY_COF(0) is the constant term and ! POLY_COF(N) is the coefficient of X**N. ! ! Input, character ( len = * ), an optional title. ! implicit none integer n integer i real mag integer n2 character plus_minus real poly_cof(0:n) character ( len = * ) title if ( 0 < len_trim ( title ) ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) trim ( title ) end if ! ! Search for the first nonzero coefficient. If none, we'll take ! the constant term. ! n2 = 0 do i = n, 1, -1 if ( poly_cof(i) /= 0.0E+00 ) then n2 = i exit end if end do write ( *, '(a)' ) ' ' if ( poly_cof(n2) < 0.0E+00 ) then plus_minus = '-' else plus_minus = ' ' end if mag = abs ( poly_cof(n2) ) if ( 2 <= n2 ) then write ( *, '( '' p(x) = '', a1, g14.6, '' * x ^ '', i3 )' ) & plus_minus, mag, n2 else if ( n2 == 1 ) then write ( *, '( '' p(x) = '', a1, g14.6, '' * x'' )' ) & plus_minus, mag else if ( n2 == 0 ) then write ( *, '( '' p(x) = '', a1, g14.6 )' ) plus_minus, mag end if do i = n2-1, 0, -1 if ( poly_cof(i) < 0.0E+00 ) then plus_minus = '-' else plus_minus = '+' end if mag = abs ( poly_cof(i) ) if ( mag /= 0.0E+00 ) then if ( 2 <= i ) then write ( *, ' ( '' '', a1, g14.6, '' * x ^ '', i3 )' ) & plus_minus, mag, i else if ( i == 1 ) then write ( *, ' ( '' '', a1, g14.6, '' * x'' )' ) plus_minus, mag else if ( i == 0 ) then write ( *, ' ( '' '', a1, g14.6 )' ) plus_minus, mag end if end if end do return end subroutine r4row_swap ( lda, m, n, a, irow1, irow2 ) !*****************************************************************************80 ! !! R4ROW_SWAP swaps two rows of an R4ROW. ! ! Modified: ! ! 14 February 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer LDA, the first dimension of A. ! ! Input, integer M, N, the number of rows and columns. ! ! Input/output, real A(LDA,N), the M by N array. ! ! Input, integer IROW1, IROW2, the two rows to swap. ! implicit none integer lda integer n real a(lda,n) integer irow1 integer irow2 integer j integer m if ( irow1 < 1 .or. m < irow1 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'R4ROW_SWAP - Fatal error!' write ( *, '(a)' ) ' IROW1 is out of range.' stop end if if ( irow2 < 1 .or. m < irow2 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'R4ROW_SWAP - Fatal error!' write ( *, '(a)' ) ' IROW2 is out of range.' stop end if if ( irow1 == irow2 ) then return end if do j = 1, n call r4_swap ( a(irow1,j), a(irow2,j) ) end do return end subroutine s_blank_delete ( s ) !*****************************************************************************80 ! !! S_BLANK_DELETE removes blanks from a string, left justifying the remainder. ! ! Comment: ! ! All TAB characters are also removed. ! ! Modified: ! ! 26 July 1998 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, character ( len = * ) S, the string to be transformed. ! implicit none character c integer iget integer iput character ( len = * ) s character, parameter :: TAB = char ( 9 ) iput = 0 do iget = 1, len ( s ) c = s(iget:iget) if ( c /= ' ' .and. c /= TAB ) then iput = iput + 1 s(iput:iput) = c end if end do s(iput+1:) = ' ' return end subroutine s_blanks_delete ( s ) !*****************************************************************************80 ! !! S_BLANKS_DELETE replaces consecutive blanks by one blank. ! ! Discussion: ! ! The remaining characters are left justified and right padded with blanks. ! TAB characters are converted to spaces. ! ! Modified: ! ! 26 July 1998 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, character ( len = * ) S, the string to be transformed. ! implicit none integer i integer j character newchr character oldchr character ( len = * ) s character, parameter :: TAB = char ( 9 ) j = 0 newchr = ' ' do i = 1, len ( s ) oldchr = newchr newchr = s(i:i) if ( newchr == TAB ) then newchr = ' ' end if s(i:i) = ' ' if ( oldchr /= ' ' .or. newchr /= ' ' ) then j = j + 1 s(j:j) = newchr end if end do return end subroutine s_chop ( string, ilo, ihi ) !*****************************************************************************80 ! !! S_CHOP "chops out" a portion of a string, and closes up the hole. ! ! Example: ! ! STRING = 'Fred is not a jerk!' ! ! call s_chop ( STRING, 9, 12 ) ! ! STRING = 'Fred is a jerk! ' ! ! Modified: ! ! 06 July 1998 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, character ( len = * ) STRING, the string to be transformed. ! ! Input, integer ILO, IHI, the locations of the first and last ! characters to be removed. ! implicit none integer ihi integer ihi2 integer ilo integer ilo2 integer lens character ( len = * ) string lens = len ( string ) ilo2 = max ( ilo, 1 ) ihi2 = min ( ihi, lens ) if ( ihi2 < ilo2 ) then return end if string(ilo2:lens+ilo2-ihi2-1) = string(ihi2+1:lens) string(lens+ilo2-ihi2:lens) = ' ' return end function s_eqi ( s1, s2 ) !*****************************************************************************80 ! !! S_EQI is a case insensitive comparison of two strings for equality. ! ! Examples: ! ! S_EQI ( 'Anjana', 'ANJANA' ) is .TRUE. ! ! Modified: ! ! 14 April 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) S1, S2, the strings to compare. ! ! Output, logical S_EQI, the result of the comparison. ! implicit none character c1 character c2 integer i integer len1 integer len2 integer lenc logical s_eqi character ( len = * ) s1 character ( len = * ) s2 len1 = len ( s1 ) len2 = len ( s2 ) lenc = min ( len1, len2 ) s_eqi = .false. do i = 1, lenc c1 = s1(i:i) c2 = s2(i:i) call ch_cap ( c1 ) call ch_cap ( c2 ) if ( c1 /= c2 ) then return end if end do do i = lenc + 1, len1 if ( s1(i:i) /= ' ' ) then return end if end do do i = lenc + 1, len2 if ( s2(i:i) /= ' ' ) then return end if end do s_eqi = .true. return end function s_indexi ( string, sub ) !*****************************************************************************80 ! !! S_INDEXI is a case-insensitive INDEX function. ! ! Discussion: ! ! The runction returns the location in STRING at which the substring SUB is ! first found, or 0 if the substring does not occur at all. ! ! The routine is also trailing blank insensitive. This is very ! important for those cases where you have stored information in ! larger variables. If STRING is of length 80, and SUB is of ! length 80, then if STRING = 'FRED' and SUB = 'RED', a match would ! not be reported by the standard FORTRAN INDEX, because it treats ! both variables as being 80 characters long! This routine assumes that ! trailing blanks represent garbage! ! ! Because of the suppression of trailing blanks, this routine cannot be ! used to find, say, the first occurrence of the two-character ! string 'A '. However, this routine treats as a special case the ! occurrence where STRING or SUB is entirely blank. Thus you can ! use this routine to search for occurrences of double or triple blanks ! in a string, for example, although INDEX itself would be just as ! suitable for that problem. ! ! Modified: ! ! 14 April 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) STRING, the string to be searched. ! ! Input, character ( len = * ) SUB, the substring to search for. ! ! Output, integer S_INDEXI. 0 if SUB does not occur in ! STRING. Otherwise STRING(S_INDEXI:S_INDEXI+LENS-1) = SUB, ! where LENS is the length of SUB, and is the first place ! this happens. However, note that S_INDEXI ignores case, ! unlike the standard FORTRAN INDEX function. ! implicit none integer i integer llen1 integer llen2 logical s_eqi integer s_indexi character ( len = * ) string character ( len = * ) sub s_indexi = 0 llen1 = len_trim ( string ) llen2 = len_trim ( sub ) ! ! In case STRING or SUB is blanks, use LEN. ! if ( llen1 == 0 ) then llen1 = len ( string ) end if if ( llen2 == 0 ) then llen2 = len ( sub ) end if if ( llen1 < llen2 ) then return end if do i = 1, llen1 + 1 - llen2 if ( s_eqi ( string(i:i+llen2-1), sub ) ) then s_indexi = i return end if end do return end function s_is_alpha ( s ) !*****************************************************************************80 ! !! S_IS_ALPHA returns .TRUE. if the string contains only alphabetic characters. ! ! Discussion: ! ! Here, alphabetic characters are 'A' through 'Z' and 'a' through 'z'. ! ! Modified: ! ! 14 April 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) S, the string to be checked. ! ! Output, logical S_IS_ALPHA, .TRUE. if the string contains only ! alphabetic characters, .FALSE. otherwise. ! implicit none logical ch_is_alpha integer i character ( len = * ) s logical s_is_alpha s_is_alpha = .false. do i = 1, len ( s ) if ( .not. ch_is_alpha ( s(i:i) ) ) then return end if end do s_is_alpha = .true. return end function s_paren_check ( s ) !*****************************************************************************80 ! !! S_PAREN_CHECK checks the parentheses in a string. ! ! Discussion: ! ! Blanks are removed from the string, and then the following checks ! are made: ! ! 1) as we read the string left to right, there must never be more ! right parentheses than left ones; ! 2) there must be an equal number of left and right parentheses; ! 3) there must be no occurrences of the abutting packages '...)(...'. ! 4) there must be no occurrences of the empty package '()'. ! ! Modified: ! ! 20 November 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) S, the string to check. ! ! Output, logical S_PAREN_CHECK is TRUE if the string passed the checks. ! implicit none integer i integer isum character ( len = * ) s character ( len = 256 ) s_copy integer s_len logical s_paren_check s_copy = s call s_blank_delete ( s_copy) s_len = len_trim ( s_copy ) ! ! 1) Letting '(' = +1 and ')' = -1, check that the running parentheses sum ! is always nonnegative. ! isum = 0 do i = 1, s_len if ( s_copy(i:i) == '(' ) then isum = isum + 1 end if if ( s_copy(i:i) == ')' ) then isum = isum - 1 if ( isum < 0 ) then s_paren_check = .false. return end if end if end do ! ! 2) Check that the final parentheses sum is zero. ! if ( isum /= 0 ) then s_paren_check = .false. return end if ! ! 3) Check that there are no ")(" pairs. ! do i = 2, s_len if ( s_copy(i-1:i) == ')(' ) then s_paren_check = .false. return end if end do ! ! 4) Check that there are no "()" pairs. ! do i = 2, s_len if ( s_copy(i-1:i) == '()' ) then s_paren_check = .false. return end if end do ! ! The checks were passed. ! s_paren_check = .true. return end subroutine s_to_i4 ( s, ival, ierror, last ) !*****************************************************************************80 ! !! S_TO_I4 reads an I4 from a string. ! ! Modified: ! ! 28 June 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) S, a string to be examined. ! ! Output, integer IVAL, the integer value read from the string. ! If the string is blank, then IVAL will be returned 0. ! ! Output, integer IERROR, an error flag. ! 0, no error. ! 1, an error occurred. ! ! Output, integer LAST, the last character of S used to make IVAL. ! implicit none character c integer i integer ierror integer isgn integer istate integer ival integer last character ( len = * ) s ierror = 0 istate = 0 isgn = 1 ival = 0 do i = 1, len_trim ( s ) c = s(i:i) ! ! Haven't read anything. ! if ( istate == 0 ) then if ( c == ' ' ) then else if ( c == '-' ) then istate = 1 isgn = -1 else if ( c == '+' ) then istate = 1 isgn = + 1 else if ( lle ( '0', c ) .and. lle ( c, '9' ) ) then istate = 2 ival = ichar ( c ) - ichar ( '0' ) else ierror = 1 return end if ! ! Have read the sign, expecting digits. ! else if ( istate == 1 ) then if ( c == ' ' ) then else if ( lle ( '0', c ) .and. lle ( c, '9' ) ) then istate = 2 ival = ichar ( c ) - ichar ( '0' ) else ierror = 1 return end if ! ! Have read at least one digit, expecting more. ! else if ( istate == 2 ) then if ( lle ( '0', c ) .and. lle ( c, '9' ) ) then ival = 10 * ival + ichar ( c ) - ichar ( '0' ) else ival = isgn * ival last = i - 1 return end if end if end do ! ! If we read all the characters in the string, see if we're OK. ! if ( istate == 2 ) then ival = isgn * ival last = len_trim ( s ) else ierror = 1 last = 0 end if return end subroutine s_to_r4 ( s, r, ierror, lchar ) !*****************************************************************************80 ! !! S_TO_R4 reads an R4 from a string. ! ! Discussion: ! ! This routine will read as many characters as possible until it reaches ! the end of the string, or encounters a character which cannot be ! part of the real number. ! ! Legal input is: ! ! 1 blanks, ! 2 '+' or '-' sign, ! 2.5 spaces ! 3 integer part, ! 4 decimal point, ! 5 fraction part, ! 6 'E' or 'e' or 'D' or 'd', exponent marker, ! 7 exponent sign, ! 8 exponent integer part, ! 9 exponent decimal point, ! 10 exponent fraction part, ! 11 blanks, ! 12 final comma or semicolon. ! ! with most quantities optional. ! ! Examples: ! ! S R ! ! '1' 1.0 ! ' 1 ' 1.0 ! '1A' 1.0 ! '12,34,56' 12.0 ! ' 34 7' 34.0 ! '-1E2ABCD' -100.0 ! '-1X2ABCD' -1.0 ! ' 2E-1' 0.2 ! '23.45' 23.45 ! '-4.2E+2' -420.0 ! '17d2' 1700.0 ! '-14e-2' -0.14 ! 'e2' 100.0 ! '-12.73e-9.23' -12.73 * 10.0**(-9.23) ! ! Modified: ! ! 12 February 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) S, the string containing the ! data to be read. Reading will begin at position 1 and ! terminate at the end of the string, or when no more ! characters can be read to form a legal real. Blanks, ! commas, or other nonnumeric data will, in particular, ! cause the conversion to halt. ! ! Output, real R, the real value that was read from the string. ! ! Output, integer IERROR, error flag. ! ! 0, no errors occurred. ! ! 1, 2, 6 or 7, the input number was garbled. The ! value of IERROR is the last type of input successfully ! read. For instance, 1 means initial blanks, 2 means ! a plus or minus sign, and so on. ! ! Output, integer LCHAR, the number of characters read from ! the string to form the number, including any terminating ! characters such as a trailing comma or blanks. ! implicit none logical ch_eqi character c integer ierror integer ihave integer isgn integer iterm integer jbot integer jsgn integer jtop integer lchar integer nchar integer ndig real r real rbot real rexp real rtop character ( len = * ) s character, parameter :: TAB = char ( 9 ) nchar = len_trim ( s ) ierror = 0 r = 0.0E+00 lchar = - 1 isgn = 1 rtop = 0.0E+00 rbot = 1.0E+00 jsgn = 1 jtop = 0 jbot = 1 ihave = 1 iterm = 0 do lchar = lchar + 1 c = s(lchar+1:lchar+1) ! ! Blank or TAB character. ! if ( c == ' ' .or. c == TAB ) then if ( ihave == 2 ) then else if ( ihave == 6 .or. ihave == 7 ) then iterm = 1 else if ( 1 < ihave ) then ihave = 11 end if ! ! Comma. ! else if ( c == ',' .or. c == ';' ) then if ( ihave /= 1 ) then iterm = 1 ihave = 12 lchar = lchar + 1 end if ! ! Minus sign. ! else if ( c == '-' ) then if ( ihave == 1 ) then ihave = 2 isgn = - 1 else if ( ihave == 6 ) then ihave = 7 jsgn = - 1 else iterm = 1 end if ! ! Plus sign. ! else if ( c == '+' ) then if ( ihave == 1 ) then ihave = 2 else if ( ihave == 6 ) then ihave = 7 else iterm = 1 end if ! ! Decimal point. ! else if ( c == '.' ) then if ( ihave < 4 ) then ihave = 4 else if ( 6 <= ihave .and. ihave <= 8 ) then ihave = 9 else iterm = 1 end if ! ! Exponent marker. ! else if ( ch_eqi ( c, 'E' ) .or. ch_eqi ( c, 'D' ) ) then if ( ihave < 6 ) then ihave = 6 else iterm = 1 end if ! ! Digit. ! else if ( ihave < 11 .and. lge ( c, '0' ) .and. lle ( c, '9' ) ) then if ( ihave <= 2 ) then ihave = 3 else if ( ihave == 4 ) then ihave = 5 else if ( ihave == 6 .or. ihave == 7 ) then ihave = 8 else if ( ihave == 9 ) then ihave = 10 end if call ch_to_digit ( c, ndig ) if ( ihave == 3 ) then rtop = 10.0E+00 * rtop + real ( ndig ) else if ( ihave == 5 ) then rtop = 10.0E+00 * rtop + real ( ndig ) rbot = 10.0E+00 * rbot else if ( ihave == 8 ) then jtop = 10 * jtop + ndig else if ( ihave == 10 ) then jtop = 10 * jtop + ndig jbot = 10 * jbot end if ! ! Anything else is regarded as a terminator. ! else iterm = 1 end if ! ! If we haven't seen a terminator, and we haven't examined the ! entire string, go get the next character. ! if ( iterm == 1 .or. nchar <= lchar + 1 ) then exit end if end do ! ! If we haven't seen a terminator, and we have examined the ! entire string, then we're done, and LCHAR is equal to NCHAR. ! if ( iterm /= 1 .and. lchar+1 == nchar ) then lchar = nchar end if ! ! Number seems to have terminated. Have we got a legal number? ! Not if we terminated in states 1, 2, 6 or 7! ! if ( ihave == 1 .or. ihave == 2 .or. ihave == 6 .or. ihave == 7 ) then ierror = ihave return end if ! ! Number seems OK. Form it. ! if ( jtop == 0 ) then rexp = 1.0E+00 else if ( jbot == 1 ) then rexp = 10.0E+00**( jsgn * jtop ) else rexp = jsgn * jtop rexp = rexp / jbot rexp = 10.0E+00**rexp end if end if r = isgn * rexp * rtop / rbot return end subroutine scopy ( n, sx, incx, sy, incy ) !*****************************************************************************80 ! !! SCOPY copies a vector, X, to a vector, Y. ! ! Modified: ! ! 15 March 2001 ! ! Parameters: ! ! Input, integer N, the number of entries to copy. ! ! Input, real SX(*), the vector to be copied. ! ! Input, integer INCX, the increment between successive ! entries of SX. ! ! Output, real SY(*), a copy of SX. ! ! Input, integer INCY, the increment between successive ! entries of SY. ! implicit none integer i integer incx integer incy integer ix integer iy integer n real sx(*) real sy(*) if ( n<=0 ) then return end if ix = 1 iy = 1 if ( incx<0)ix = (-n+1)*incx + 1 if ( incy<0)iy = (-n+1)*incy + 1 do i = 1, n sy(iy) = sx(ix) ix = ix + incx iy = iy + incy end do return end subroutine sge_check ( lda, m, n, ierror ) !*****************************************************************************80 ! !! SGE_CHECK checks the dimensions of a general matrix. ! ! Modified: ! ! 11 January 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer LDA, the leading dimension of the array. ! LDA must be at least M. ! ! Input, integer M, the number of rows of the matrix. ! M must be positive. ! ! Input, integer N, the number of columns of the matrix. ! N must be positive. ! ! Output, integer IERROR, reports whether any errors were detected. ! IERROR is set to 0 before the checks are made, and then: ! IERROR = IERROR + 1 if LDA is illegal; ! IERROR = IERROR + 2 if M is illegal; ! IERROR = IERROR + 4 if N is illegal. ! implicit none integer ierror integer lda integer m integer n ierror = 0 if ( lda < m ) then ierror = ierror + 1 write ( *, '(a)' ) ' ' write ( *, '(a,i8)' ) 'SGE_CHECK - Illegal LDA = ', lda end if if ( m < 1 ) then ierror = ierror + 2 write ( *, '(a)' ) ' ' write ( *, '(a,i8)' ) 'SGE_CHECK - Illegal M = ', m end if if ( n < 1 ) then ierror = ierror + 4 write ( *, '(a)' ) ' ' write ( *, '(a,i8)') 'SGE_CHECK - Illegal N = ', n end if return end subroutine sge_co ( a, lda, n, ipvt, rcond ) !*****************************************************************************80 ! !! SGE_CO factors a real matrix and estimates its condition number. ! ! Discussion: ! ! If the condition number is not needed, then SGE_FA is slightly ! faster that SGE_CO. ! ! Modified: ! ! 15 March 2001 ! ! Parameters: ! ! Input/output, real A(LDA,N). ! On input, the matrix to be factored. ! On output, A has been overwritten by information that ! defines the LU factors of A. ! ! Input, integer LDA, the leading dimension of A, which ! must be at least N. ! ! Input, integer N, the number of rows and columns of A. ! ! Output, integer IPVT(N), the pivot row indices. ! ! Output, real RCOND, an estimate of the reciprocal of the ! condition number of A. For the system A*X = B, relative ! perturbations in A and B of size EPSILON may cause ! relative perturbations in X of size EPSILON/RCOND. ! If RCOND is so small that the logical expression ! 1.0E+00 + RCOND == 1.0E+00 ! is true, then A may be singular to working precision. In ! particular, RCOND is zero if exact singularity is detected ! or the estimate underflows. ! implicit none integer lda integer n real a(lda,n) real anorm real ek integer info integer ipvt(n) integer j integer k integer l real rcond real s real sm real t real wk real wkm real ynorm real z(n) ! ! Compute the L1-norm of A. ! anorm = 0.0E+00 do j = 1, n anorm = max ( anorm, sum ( abs ( a(1:n,j) ) ) ) end do ! ! Factor the matrix A. ! call sge_fa ( a, lda, n, ipvt, info ) ! ! RCOND = 1 / ( norm(a) * (estimate of norm(inverse(a)))) . ! ! Estimate = norm(z)/norm(y) where a*z = y and trans(a)*y = e . ! trans(a) is the transpose of a . the components of e are ! chosen to cause maximum local growth in the elements of w where ! trans(u)*w = e . ! ! The vectors are frequently rescaled to avoid overflow. ! ! solve trans(u)*w = e ! ek = 1.0E+00 z(1:n) = 0.0E+00 do k = 1, n if ( z(k) /= 0.0E+00 ) then ek = sign ( ek, - z(k) ) end if if ( abs ( a(k,k) ) < abs ( ek - z(k) ) ) then s = abs ( a(k,k) ) / abs ( ek - z(k) ) z(1:n) = s * z(1:n) ek = s * ek end if wk = ek - z(k) wkm = -ek - z(k) s = abs ( wk ) sm = abs ( wkm ) if ( a(k,k) /= 0.0E+00 ) then wk = wk / a(k,k) wkm = wkm / a(k,k) else wk = 1.0E+00 wkm = 1.0E+00 end if sm = sm + sum ( abs ( z(k+1:n) + wkm * a(k,k+1:n) ) ) z(k+1:n) = z(k+1:n) + wk * a(k,k+1:n) s = s + sum ( abs ( z(k+1:n) ) ) if ( s < sm ) then t = wkm - wk wk = wkm z(k+1:n) = z(k+1:n) + t * a(k,k+1:n) end if z(k) = wk end do s = sum ( abs ( z(1:n) ) ) z(1:n) = z(1:n) / s ! ! Solve trans(l)*y = w. ! do k = n, 1, -1 if ( k < n ) then z(k) = z(k) + dot_product ( a(k+1:n,k), z(k+1:n) ) end if if ( 1.0E+00 < abs ( z(k) ) ) then z(1:n) = z(1:n) / abs ( z(k) ) end if l = ipvt(k) call r4_swap ( z(l), z(k) ) end do s = sum ( abs ( z(1:n) ) ) z(1:n) = z(1:n) / s ynorm = 1.0E+00 ! ! Solve l*v = y ! do k = 1, n l = ipvt(k) call r4_swap ( z(l), z(k) ) if ( k < n ) then z(k+1:n) = z(k+1:n) + z(k) * a(k+1:n,k) end if if ( 1.0E+00 < abs ( z(k) ) ) then z(1:n) = z(1:n) / abs ( z(k) ) ynorm = s * ynorm end if end do s = sum ( abs ( z(1:n) ) ) z(1:n) = z(1:n) / s ynorm = ynorm / s ! ! Solve u*z = v. ! do k = n, 1, -1 if ( abs ( a(k,k) ) < abs ( z(k) ) ) then s = abs ( a(k,k) ) / abs ( z(k) ) z(1:n) = z(1:n) * s ynorm = s * ynorm end if if ( a(k,k) /= 0.0E+00 ) then z(k) = z(k) / a(k,k) else z(k) = 1.0E+00 end if z(1:k-1) = z(1:k-1) - z(k) * a(1:k-1,k) end do ! ! Make ZNORM = 1. ! s = sum ( abs ( z(1:n) ) ) z(1:n) = z(1:n) / s ynorm = ynorm / s if ( anorm /= 0.0E+00 ) then rcond = ynorm / anorm else rcond = 0.0E+00 end if return end subroutine sge_fa ( a, lda, n, ipivot, info ) !*****************************************************************************80 ! !! SGE_FA factors a general matrix. ! ! Note: ! ! SGE_FA is a simplified version of the LINPACK routine SGEFA. ! ! Modified: ! ! 26 February 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer LDA, the leading dimension of the array. ! LDA must be at least N. ! ! Input, integer N, the order of the matrix. ! N must be positive. ! ! Input/output, real A(LDA,N), the matrix to be factored. ! On output, A contains an upper triangular matrix and the multipliers ! which were used to obtain it. The factorization can be written ! A = L * U, where L is a product of permutation and unit lower ! triangular matrices and U is upper triangular. ! ! Output, integer IPIVOT(N), a vector of pivot indices. ! ! Output, integer INFO, singularity flag. ! 0, no singularity detected. ! nonzero, the factorization failed on the INFO-th step. ! implicit none integer lda integer n real a(lda,n) integer i integer ierror integer info integer ipivot(n) integer j integer k integer l ! ! Check the dimensions. ! call sge_check ( lda, n, n, ierror ) if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SGE_FA - Fatal error!' write ( *, '(a)' ) ' Illegal dimensions.' return end if info = 0 do k = 1, n-1 ! ! Find L, the index of the pivot row. ! l = k do i = k+1, n if ( abs ( a(l,k) ) < abs ( a(i,k) ) ) then l = i end if end do ipivot(k) = l ! ! If the pivot index is zero, the algorithm has failed. ! if ( a(l,k) == 0.0E+00 ) then info = k write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SGE_FA - Fatal error!' write ( *, '(a,i8)' ) ' Zero pivot on step ', info return end if ! ! Interchange rows L and K if necessary. ! if ( l /= k ) then call r4_swap ( a(l,k), a(k,k) ) end if ! ! Normalize the values that lie below the pivot entry A(K,K). ! a(k+1:n,k) = -a(k+1:n,k) / a(k,k) ! ! Row elimination with column indexing. ! do j = k+1, n if ( l /= k ) then call r4_swap ( a(l,j), a(k,j) ) end if a(k+1:n,j) = a(k+1:n,j) + a(k+1:n,k) * a(k,j) end do end do ipivot(n) = n if ( a(n,n) == 0.0E+00 ) then info = n write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SGE_FA - Fatal error!' write ( *, '(a,i8)' ) ' Zero pivot on step ', info end if return end subroutine sge_inv ( lda, n, a, ipivot ) !*****************************************************************************80 ! !! SGE_INV computes the inverse of a matrix factored by SGE_FA. ! ! Note: ! ! SGE_INV is a simplified standalone version of the LINPACK routine ! SGEDI. ! ! Modified: ! ! 16 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer LDA, the leading dimension of the array A, ! which must be at least N. ! ! Input, integer N, the order of the matrix A. ! ! Input/output, real A(LDA,N). ! On input, the factor information computed by SGE_FA. ! On output, the inverse matrix. ! ! Input, integer IPIVOT(N), the pivot vector from SGE_FA. ! implicit none integer lda integer n real a(lda,n) integer i integer ierror integer ipivot(n) integer j integer k real temp real work(n) ! ! Check the dimensions. ! call sge_check ( lda, n, n, ierror ) if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SGE_INV - Fatal error!' write ( *, '(a)' ) ' Illegal dimensions.' return end if ! ! Compute Inverse(U). ! do k = 1, n a(k,k) = 1.0E+00 / a(k,k) a(1:k-1,k) = -a(1:k-1,k) * a(k,k) do j = k + 1, n temp = a(k,j) a(k,j) = 0.0E+00 a(1:k,j) = a(1:k,j) + temp * a(1:k,k) end do end do ! ! Form Inverse(U) * Inverse(L). ! do k = n - 1, 1, -1 work(k+1:n) = a(k+1:n,k) a(k+1:n,k) = 0.0E+00 do j = k + 1, n a(1:n,k) = a(1:n,k) + a(1:n,j) * work(j) end do if ( ipivot(k) /= k ) then do i = 1, n call r4_swap ( a(i,k), a(i,ipivot(k)) ) end do end if end do return end subroutine sge_plu ( lda, m, n, a, p, l, u ) !*****************************************************************************80 ! !! SGE_PLU produces the PLU factors of a real rectangular matrix. ! ! Discussion: ! ! The PLU factors of the M by N matrix A are: ! ! P, an M by M permutation matrix P, ! L, an M by M unit lower triangular matrix, ! U, an M by N upper triangular matrix. ! ! Modified: ! ! 30 April 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer LDA, the leading dimension of A. ! ! Input, integer M, the number of rows in A. ! ! Input, integer N, the number of columns in A. ! ! Input, real A(LDA,N), the M by N matrix to be factored. ! ! Output, real P(LDA,M), the M by M permutation factor. ! ! Output, real L(LDA,M), the M by M unit lower triangular factor. ! ! Output, real U(LDA,N), the M by N upper triangular factor. ! implicit none integer lda integer m integer n real a(lda,n) integer i integer j integer k real l(lda,m) real p(lda,m) integer pivot_row real pivot_value real u(lda,n) ! ! Initialize: ! ! P:=Identity ! L:=Identity ! U:=A ! call r4mat_identity ( lda, m, p ) l = p u = a ! ! On step J, find the pivot row and the pivot value. ! do j = 1, min ( m-1, n ) pivot_value = 0.0E+00 pivot_row = 0 do i = j, m if ( pivot_value < abs ( u(i,j) ) ) then pivot_value = abs ( u(i,j) ) pivot_row = i end if end do ! ! If the pivot row is nonzero, swap rows J and PIVOT_ROW. ! if ( pivot_row /= 0 ) then call r4row_swap ( lda, m, n, u, j, pivot_row ) call r4row_swap ( lda, m, m, l, j, pivot_row ) call r4col_swap ( lda, m, m, l, j, pivot_row ) call r4col_swap ( lda, m, m, p, j, pivot_row ) ! ! Zero out the entries in column J, from row J+1 to M. ! do i = j+1, m if ( u(i,j) /= 0.0E+00 ) then l(i,j) = u(i,j) / u(j,j) u(i,j) = 0.0E+00 do k = j+1, n u(i,k) = u(i,k) - l(i,j) * u(j,k) end do end if end do end if end do return end subroutine sge_sl ( a, lda, n, ipivot, b, job ) !*****************************************************************************80 ! !! SGE_SL solves a system factored by SGE_FA. ! ! Discussion: ! ! SGE_SL is a simplified version of the LINPACK routine SGESL. ! ! Modified: ! ! 04 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer LDA, the leading dimension of the array. ! LDA must be at least N. ! ! Input, integer N, the order of the matrix. ! N must be positive. ! ! Input, real A(LDA,N), the LU factors from SGE_FA. ! ! Input, integer IPIVOT(N), the pivot vector from SGE_FA. ! ! Input/output, real B(N). ! On input, the right hand side vector. ! On output, the solution vector. ! ! Input, integer JOB, specifies the operation. ! 0, solve A * x = b. ! nonzero, solve transpose ( A ) * x = b. ! implicit none integer lda integer n real a(lda,n) real b(n) integer ierror integer ipivot(n) integer job integer k integer l ! ! Check the dimensions. ! call sge_check ( lda, n, n, ierror ) if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SGE_SL - Fatal error!' write ( *, '(a)' ) ' Illegal dimensions.' return end if ! ! Solve A * x = b. ! if ( job == 0 ) then ! ! Solve PL * Y = B. ! do k = 1, n-1 l = ipivot(k) if ( l /= k ) then call r4_swap ( b(l), b(k) ) end if b(k+1:n) = b(k+1:n) + a(k+1:n,k) * b(k) end do ! ! Solve U * X = Y. ! do k = n, 1, -1 b(k) = b(k) / a(k,k) b(1:k-1) = b(1:k-1) - a(1:k-1,k) * b(k) end do ! ! Solve transpose ( A ) * X = B. ! else ! ! Solve transpose ( U ) * Y = B. ! do k = 1, n b(k) = ( b(k) - dot_product ( b(1:k-1), a(1:k-1,k) ) ) / a(k,k) end do ! ! Solve transpose ( PL ) * X = Y. ! do k = n-1, 1, -1 b(k) = b(k) + dot_product ( b(k+1:n), a(k+1:n,k) ) l = ipivot(k) if ( l /= k ) then call r4_swap ( b(l), b(k) ) end if end do end if return end subroutine svd ( lda, m, n, a, w, matu, u, matv, v, ierror ) !*****************************************************************************80 ! !! SVD determines the singular value decomposition of a matrix. ! ! Discussion: ! ! The decomposition has the form ! ! A = U * S * V' ! ! where U is an M by N matrix whose columns are orthogonal, and V is an ! N by N orthogonal matrix, and S is an N by N diagonal matrix, whose ! zero entries, if any, come last. ! ! Householder bidiagonalization and a variant of the QR algorithm are used. ! ! Reference: ! ! Golub and Reinsch, ! Numerische Mathematik, ! Volume 14, pages 403-420, 1970, ! ! Handbook for Automatic Computation, ! Volume II, Linear Algebra, pages 134-151, 1971. ! ! Modified: ! ! 15 March 2001 ! ! Parameters: ! ! Input, integer LDA, the leading dimension of the arrays ! A, U and V. LDA must be at least as large as the maximum ! of M and N. ! ! Input, integer M, the number of rows of A and U. ! ! Input, integer N, the number of columns of A and U, ! and the number of rows and columns of V. ! ! Input, real A(LDA,M), the M by N matrix to be decomposed. ! ! Output, real W(N), contains the N (non-negative) singular ! values of A, that is, the diagonal elements of S. They are ! unordered. If an error exit is made, the singular values ! should be correct for indices IERROR+1, IERROR+2,...,N. ! ! Input, logical MATU, set to .TRUE. if the U matrix in the ! decomposition is desired, and .FALSE. otherwise. ! ! Output, real U(NM,N), contains the matrix U, of orthogonal ! column vectors, of the decomposition if MATU has been set to ! .TRUE. Otherwise, U is used as a temporary array. U may ! coincide with A. If an error exit is made, the columns of ! U corresponding to indices of correct singular values should ! be correct. ! ! Input, logical MATV, set to .TRUE. if the V matrix in the ! decomposition is desired, and .FALSE. otherwise. ! ! Output, real V(LDA,N), contains the orthogonal matrix V of the ! decomposition if MATV has been set to .TRUE. Otherwise V is ! not referenced. V may also coincide with A if U is not needed. ! If an error exit is made, the columns of V corresponding to ! indices of correct singular values should be correct. ! ! Output, integer IERROR, error flag. ! 0, no error. ! K, if the K-th singular value has not been determined after ! 30 iterations. ! implicit none integer lda integer n real a(lda,n) real c real f real g real h integer i integer i1 integer ierror integer ii integer its integer j integer k integer k1 integer kk integer l integer l1 integer ll integer m logical matu logical matv integer mn real pythag real rv1(n) real s real scale real temp real tst1 real tst2 real u(lda,n) real v(lda,n) real w(n) real x real y real z ierror = 0 ! ! Copy A into U. ! u(1:m,1:n) = a(1:m,1:n) ! ! Householder reduction to bidiagonal form. ! g = 0.0E+00 scale = 0.0E+00 x = 0.0E+00 do i = 1, n l = i + 1 rv1(i) = scale * g g = 0.0E+00 s = 0.0E+00 scale = 0.0E+00 if ( i <= m ) then do k = i, m scale = scale + abs ( u(k,i) ) end do if ( scale == 0.0E+00 ) then go to 210 end if do k = i, m u(k,i) = u(k,i) / scale s = s + u(k,i)**2 end do f = u(i,i) g = - sign ( sqrt ( s ), f ) h = f * g - s u(i,i) = f - g if ( i /= n ) then do j = l, n s = 0.0E+00 do k = i, m s = s + u(k,i) * u(k,j) end do f = s / h do k = i, m u(k,j) = u(k,j) + f * u(k,i) end do end do end if u(i:m,i) = scale * u(i:m,i) end if 210 continue w(i) = scale * g g = 0.0E+00 s = 0.0E+00 scale = 0.0E+00 if ( m < i .or. i == n ) then go to 290 end if do k = l, n scale = scale + abs ( u(i,k) ) end do if ( scale == 0.0E+00 ) then go to 290 end if do k = l, n u(i,k) = u(i,k) / scale s = s + u(i,k)**2 end do f = u(i,l) g = - sign ( sqrt ( s ), f ) h = f * g - s u(i,l) = f - g do k = l, n rv1(k) = u(i,k) / h end do if ( i /= m ) then do j = l, m s = 0.0E+00 do k = l, n s = s + u(j,k) * u(i,k) end do do k = l, n u(j,k) = u(j,k) + s * rv1(k) end do end do end if u(i,l:n) = scale * u(i,l:n) 290 continue x = max ( x, abs ( w(i) ) + abs ( rv1(i) ) ) end do ! ! Accumulation of right-hand transformations. ! if (.not. matv ) then go to 410 end if do ii = 1, n i = n + 1 - ii if ( i == n ) then go to 390 end if if ( g /= 0.0E+00 ) then do j = l, n v(j,i) = (u(i,j) / u(i,l)) / g end do do j = l, n s = 0.0E+00 do k = l, n s = s + u(i,k) * v(k,j) end do do k = l, n v(k,j) = v(k,j) + s * v(k,i) end do end do end if do j = l, n v(i,j) = 0.0E+00 v(j,i) = 0.0E+00 end do 390 continue v(i,i) = 1.0E+00 g = rv1(i) l = i end do ! ! Accumulation of left-hand transformations. ! 410 continue if ( .not. matu ) then go to 510 end if mn = n if ( m < n ) mn = m do ii = 1, mn i = mn + 1 - ii l = i + 1 g = w(i) if ( i /= n ) then do j = l, n u(i,j) = 0.0E+00 end do end if if ( g == 0.0E+00 ) then go to 475 end if if ( i /= mn ) then do j = l, n s = 0.0E+00 do k = l, m s = s + u(k,i) * u(k,j) end do f = ( s / u(i,i)) / g do k = i, m u(k,j) = u(k,j) + f * u(k,i) end do end do end if do j = i, m u(j,i) = u(j,i) / g end do go to 490 475 continue do j = i, m u(j,i) = 0.0E+00 end do 490 continue u(i,i) = u(i,i) + 1.0E+00 end do ! ! Diagonalization of the bidiagonal form. ! 510 continue tst1 = x do kk = 1, n k1 = n - kk k = k1 + 1 its = 0 ! ! Test for splitting. ! 520 continue do ll = 1, k l1 = k - ll l = l1 + 1 tst2 = tst1 + abs ( rv1(l) ) if ( tst2 == tst1 ) then go to 565 end if ! ! rv1(1) is always zero, so there is no exit through the bottom ! of the loop. ! tst2 = tst1 + abs ( w(l1) ) if ( tst2 == tst1 ) then go to 540 end if end do ! ! Cancellation of rv1(l) if l greater than 1. ! 540 continue c = 0.0E+00 s = 1.0E+00 do i = l, k f = s * rv1(i) rv1(i) = c * rv1(i) tst2 = tst1 + abs ( f ) if ( tst2 == tst1 ) then go to 565 end if g = w(i) h = pythag ( f, g ) w(i) = h c = g / h s = -f / h if (matu ) then do j = 1, m y = u(j,l1) z = u(j,i) u(j,l1) = y * c + z * s u(j,i) = -y * s + z * c end do end if end do ! ! Test for convergence. ! 565 continue z = w(k) if ( l == k ) then go to 650 end if ! ! Shift from bottom 2 by 2 minor. ! if ( 30 <= its ) then ierror = k return end if its = its + 1 x = w(l) y = w(k1) g = rv1(k1) h = rv1(k) f = 0.5E+00 * (((g + z) / h) * ((g - z) / y) + y / h - h / y) temp = 1.0E+00 g = pythag ( f, temp ) f = x - (z / x) * z + (h / x) * (y / (f + sign(g,f)) - h) ! ! Next QR transformation. ! c = 1.0E+00 s = 1.0E+00 do i1 = l, k1 i = i1 + 1 g = rv1(i) y = w(i) h = s * g g = c * g z = pythag(f,h) rv1(i1) = z c = f / z s = h / z f = x * c + g * s g = -x * s + g * c h = y * s y = y * c if ( matv ) then do j = 1, n x = v(j,i1) z = v(j,i) v(j,i1) = x * c + z * s v(j,i) = -x * s + z * c end do end if z = pythag ( f, h ) w(i1) = z ! ! Rotation can be arbitrary if z is zero. ! if ( z /= 0.0E+00 ) then c = f / z s = h / z end if f = c * g + s * y x = -s * g + c * y if ( matu ) then do j = 1, m y = u(j,i1) z = u(j,i) u(j,i1) = y * c + z * s u(j,i) = -y * s + z * c end do end if end do rv1(l) = 0.0E+00 rv1(k) = f w(k) = x go to 520 ! ! Convergence. ! 650 continue ! ! W(K) is made non-negative. ! if ( z < 0.0E+00 ) then w(k) = -z if ( matv ) then v(1:n,k) = -v(1:n,k) end if end if end do return end subroutine symbol_add ( ierror, ifree, iopsym, iprsym, ipval, maxrow, & maxsym, maxval, namvr, ncol, nrow, nsym, nsymp, numdim, symbol, valsym ) !*****************************************************************************80 ! !! SYMBOL_ADD adds a symbol name to the list of symbolic names. ! ! Modified: ! ! 10 January 2004 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer IERROR, error flag. ! 0, no error occurred, variable was added to list. ! 1, error occurred, variable was not added to list. ! ! Input/output, integer IFREE, the next free location in VALSYM. ! ! Input, integer IOPSYM(MAXSYM), the number of operands per symbol. ! Constants have 0, unaries have 1, and so on. ! ! Input, integer IPRSYM(MAXSYM), the priority of each symbol. ! ! Input, integer IPVAL(MAXSYM), contains, for each symbol, ! the address in VALSYM where the value of the symbol ! is stored, if it is a scalar. If the symbol represents ! a vector or matrix, IPVAL(IARG) points to location of ! the first entry. ! ! Input, integer MAXROW, the maximum number of rows and columns. ! ! Input, integer MAXSYM, the maximum number of symbols allowed. ! ! Input, integer MAXVAL, the maximum number of values. ! ! Input, character ( len = MAXLEN ) NAMVR, the name of the variable. ! ! Input, integer NCOL, the number of columns in the variable. ! ! Input, integer NROW, the number of rows in the variable. ! ! Input, integer NSYM, the total number of symbols, including temporaries. ! ! Input, integer NSYMP, the number of permanent symbols. ! ! Input, integer NUMDIM(2,MAXSYM). ! For each symbol I, NUMDIM(1,I) is the number of rows in its ! value, and NUMDIM(2,I) is the number of columns. ! ! Input, character ( len = MAXLEN ) SYMBOL(MAXSYM), the symbolic names ! of the internally defined constants and functions. ! ! Input, real VALSYM(MAXVAL), the values of the symbolic variables. ! implicit none integer, parameter :: maxlen = 20 integer maxsym integer maxval integer i integer ierror integer ifree integer indx integer iopsym(maxsym) integer iprsym(maxsym) integer ipval(maxsym) integer maxrow character ( len = maxlen ) namvr integer ncol integer nrow integer nsym integer nsymp integer numdim(2,maxsym) logical s_eqi character ( len = maxlen ) symbol(maxsym) real valsym(maxval) ! ! Get the length of the variable name. ! ierror = 0 if ( len_trim ( namvr ) <= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SYMBOL_ADD - Error!' write ( *, '(a)' ) ' This variable name has zero length!' ierror = 1 return end if if ( ncol <= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SYMBOL_ADD - Error!' write ( *, '(a)' ) ' This variable name has a nonpositive NCOL!' ierror = 1 return end if if ( nrow <= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SYMBOL_ADD - Error!' write ( *, '(a)' ) ' This variable name has a nonpositive NROW!' ierror = 1 return end if ! ! Check if the name is already in use. ! do i = 1, nsymp if ( s_eqi ( namvr, symbol(i) ) ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SYMBOL_ADD - Error!' write ( *, '(a)' ) ' The name ' // trim ( namvr ) // ' is already in use.' ierror = 1 return end if end do ! ! Check for user overriding earlier use of same name. ! do i = nsymp+1, nsym if ( s_eqi ( namvr, symbol(i) ) ) then if ( numdim(1,i) == nrow .and. numdim(2,i) == ncol ) then return end if if ( numdim(1,i) * numdim(2,i) == nrow * ncol ) then numdim(1,i) = nrow numdim(2,i) = ncol return end if symbol(i) = 'VOID' end if end do ! ! Is there enough space in SYMBOL for another symbol name? ! if ( maxsym <= nsym ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SYMBOL_ADD - Error!' write ( *, '(a)' ) ' There is not enough memory for ' // trim ( namvr ) write ( *, '(a)' ) ' Perhaps the KILL or INIT commands would help.' ierror = 1 return end if ! ! Check dimensions. ! if ( nrow < 1 ) then ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SYMBOL_ADD - Error!' write ( *, '(a)' ) ' Nonpositive row dimension specified.' return end if if ( ncol < 1 ) then ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SYMBOL_ADD - Error!' write ( *, '(a)' ) ' Nonpositive column dimension specified.' return end if ! ! Is either dimension too large? ! if ( maxrow < nrow ) then ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SYMBOL_ADD - Error!' write ( *, '(a)' ) ' The variable has too many rows!' write ( *, '(a,i8)' ) ' The maximum allowed is ', maxrow write ( *, '(a,i8)' ) ' Your variable requests ', nrow return end if if ( maxrow < ncol ) then ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SYMBOL_ADD - Error!' write ( *, '(a)' ) ' The variable has too many columns!' write ( *, '(a,i8)' ) ' The maximum allowed is ', maxrow write ( *, '(a,i8)' ) ' Your variable requests ', ncol return end if ! ! Is there enough space in VALSYM for the symbol's value? ! if ( maxval < ifree + nrow * ncol - 1 ) then ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SYMBOL_ADD - Error!' write ( *, '(a)' ) ' There is not enough memory for ' // trim ( namvr ) write ( *, '(a)' ) ' Perhaps the KILL or INIT commands would help.' return end if ! ! Insert information about the symbol into various tables. ! nsym = nsym + 1 symbol(nsym) = namvr iopsym(nsym) = 0 iprsym(nsym) = 10 numdim(1,nsym) = nrow numdim(2,nsym) = ncol ipval(nsym) = ifree ifree = ifree + nrow * ncol indx = ipval(nsym) valsym(indx:indx+nrow*ncol-1) = 0.0E+00 return end subroutine symbol_value ( com, ierror, ifree, iopsym, iprsym, ipval, maxrow, & maxsym, maxval, namvar, ncol, nrow, nsym, nsymp, numdim, symbol, valsym, & value ) !*****************************************************************************80 ! !! SYMBOL_VALUE sets, evaluates, or deletes a variable. ! ! Modified: ! ! 17 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character COM, determines what task is to be done: ! 'R', "reads" the variable, returning its value and dimensions. ! 'V', "sets" the value and dimensions of the variable. ! 'D', "deletes" the variable from memory. ! ! Output, integer IERROR, error flag. ! 0, no error occurred. ! nonzero, an error occurred. ! ! Input/output, integer IFREE, the next free memory location in VALSYM. ! ! Input, integer IOPSYM(MAXSYM), the number of operands per symbol. ! Constants have 0, unaries have 1, and so on. ! ! Input, integer IPRSYM(MAXSYM), priority of each symbol. ! ! Input, integer IPVAL(MAXSYM), contains, for each symbol, ! the address in VALSYM where the value of the symbol ! is stored, if it is a scalar. If the symbol represents ! a vector or matrix, IPVAL(IARG) points to location of ! the first entry. ! ! Input, integer MAXROW, the maximum number of rows and columns. ! ! Input, integer MAXSYM, the maximum number of symbols allowed. ! ! Input, integer MAXVAL, the maximum number of values allowed. ! ! Input, character ( len = MAXLEN ) NAMVAR, the name of the variable. ! ! Input/output, integer NCOL, the number of columns in the variable. ! If COM is 'V', then NCOL is an input quantity. ! If COM is 'R', then NCOL is an output quantity. ! ! Input/output, integer NROW, the number of rows in the variable. ! If COM is 'V', then NROW is an input quantity. ! If COM is 'R', then NROW is an output quantity. ! ! Input, integer NSYM, the total number of symbols, including temporaries. ! ! Input, integer NSYMP, the number of permanent symbols. ! ! Input, integer NUMDIM(2,MAXSYM). ! For each symbol I, NUMDIM(1,I) is the number of rows in its ! value, and NUMDIM(2,I) is the number of columns. ! ! Input, character ( len = MAXLEN ) SYMBOL(MAXSYM), the symbolic names ! of the internally defined constants and functions. ! ! Input, real VALSYM(MAXVAL), the values of the symbolic variables. ! ! Workspace, real VALUE(MAXROW,MAXROW), the value of the variable to ! be saved. ! implicit none integer, parameter :: maxlen = 20 integer maxrow integer maxsym integer maxval character com integer i integer ierror integer ifree integer ihi integer ilo integer indx integer iopsym(maxsym) integer iprsym(maxsym) integer ipval(maxsym) integer j integer match character ( len = maxlen ) namvar integer ncol integer nfree integer nrow integer nsym integer nsymp integer numdim(2,maxsym) logical s_eqi character ( len = maxlen ) symbol(maxsym) real valsym(maxval) real value(maxrow,maxrow) ierror = 0 if ( len_trim ( namvar ) <= 0 ) then ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SYMBOL_VALUE - Error!' write ( *, '(a)' ) ' Cannot understand the variable ' // trim ( namvar ) return end if ! ! Search for a match between NAMVAR and any predefined symbol. ! match = 0 do i = 1, nsym if ( s_eqi ( namvar, symbol(i) ) ) then match = i exit end if end do if ( match == 0 ) then ! ! If no such variable seems to exist, but we have been asked ! to DELETE or READ such a variable, then exit. ! if ( s_eqi ( com, 'D' ) .or. s_eqi ( com, 'R' ) ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SYMBOL_VALUE - Error!' write ( *, '(a)' ) ' Could not find a variable named ' // trim ( namvar ) if ( s_eqi ( com, 'R' ) ) then ierror = 1 return end if end if ! ! Add the name of the new variable to the internal list. ! call symbol_add ( ierror, ifree, iopsym, iprsym, ipval, maxrow, maxsym, & maxval, namvar, ncol, nrow, nsym, nsymp, numdim, symbol, valsym ) if ( ierror /= 0 ) then return end if match = nsym else ! ! The matched symbol has index MATCH. ! ! If the symbol is not the name of a variable, but ! rather the name of an operator, then this is an error. ! if ( iopsym(match) /= 0 ) then ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SYMBOL_VALUE - Error!' write ( *, '(a)' ) ' You are assigning a value to a function or operator:' write ( *, '(a)' ) trim ( symbol(match) ) return end if ! ! Some variables may not be revalued. ! if ( s_eqi ( com, 'V' ) ) then if ( s_eqi ( symbol(match), 'E' ) ) then ierror = 1 else if ( s_eqi ( symbol(match), 'EPS' ) ) then ierror = 1 else if ( s_eqi ( symbol(match), 'HUGE_I' ) ) then ierror = 1 else if ( s_eqi ( symbol(match), 'HUGE_R' ) ) then ierror = 1 else if ( s_eqi ( symbol(match), 'PI' ) ) then ierror = 1 else if ( symbol(match)(1:1) == '"' ) then ierror = 1 end if end if ! ! Some variables may not be deleted. ! if ( s_eqi ( com, 'D' ) ) then if ( s_eqi ( symbol(match), 'E' ) ) then ierror = 1 else if ( s_eqi ( symbol(match), 'EPS' ) ) then ierror = 1 else if ( s_eqi ( symbol(match), 'HUGE_I' ) ) then ierror = 1 else if ( s_eqi ( symbol(match), 'HUGE_R' ) ) then ierror = 1 else if ( s_eqi ( symbol(match), 'PI' ) ) then ierror = 1 else if ( s_eqi ( symbol(match), 'RTOL' ) ) then ierror = 1 else if ( symbol(match)(1:1) == '"' ) then ierror = 1 end if end if if ( ierror /= 0 ) then ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SYMBOL_VALUE - Error!' write ( *, '(a)' ) ' You are trying to change a special variable:' write ( *, '(a)' ) trim ( symbol(match) ) return end if nrow = numdim(1,match) ncol = numdim(2,match) end if ! ! Now carry out the requested operation. ! indx = ipval(match) if ( s_eqi ( com, 'V' ) ) then do j = 1, ncol valsym(indx:indx+nrow-1) = value(1:nrow,j) indx = indx + nrow end do else if ( s_eqi ( com, 'R' ) ) then do j = 1, ncol value(1:nrow,j) = valsym(indx:indx+nrow-1) indx = indx + nrow end do else if ( s_eqi ( com, 'D' ) ) then ilo = indx + nrow * ncol - 1 ihi = ifree - 1 call scopy ( ihi+1-ilo, valsym(ilo), 1, valsym(ilo-nrow*ncol), 1 ) ifree = ifree - nrow * ncol do i = match+1, nsym iopsym(i-1) = iopsym(i) iprsym(i-1) = iprsym(i) ipval(i-1) = ipval(i) - nrow * ncol numdim(1,i-1) = numdim(1,i) numdim(2,i-1) = numdim(2,i) symbol(i-1) = symbol(i) end do nsym = nsym - 1 nfree = maxval + 1 - ifree write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SYMBOL_VALUE - Note:' write ( *, '(a)' ) ' Removed the variable ' // trim ( namvar ) write ( *, '(a,i8)' ) ' The amount of free memory is now ', nfree end if return end subroutine symbol_list ( ifree, intsym, iopsym, iprsym, ipval, irpn, maxfrm, & maxrpn, maxsym, maxval, namvr, nints, nsym, nsymp, nsyms, numdim, symbol, & valsym ) !*****************************************************************************80 ! !! SYMBOL_LIST prints out information about a variable. ! ! Modified: ! ! 17 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer IFREE, the next free memory location in VALSYM. ! ! Input, integer INTSYM(80), the indices of the symbols in INFIX. ! ! Input, integer IOPSYM(MAXSYM), the number of operands per symbol. ! Constants have 0, unaries have 1, and so on. ! ! Input, integer IPRSYM(MAXSYM), priority of each symbol. ! ! Input, integer IPVAL(MAXSYM), contains, for each symbol, ! the address in VALSYM where the value of the symbol ! is stored, if it is a scalar. If the symbol represents ! a vector or matrix, IPVAL(IARG) points to location of ! the first entry. ! ! Workspace, integer IRPN(MAXRPN), the compiled versions of user formulas. ! ! Input, integer MAXFRM, the maximum number of formulas. ! ! Input, integer MAXRPN, specifies the length of IRPN. ! ! Input, integer MAXSYM, the maximum number of symbols allowed. ! ! Input, integer MAXVAL, the maximum number of values allowed. ! ! Input, character ( len = MAXLEN ) NAMVR, the name of the variable ! for which information is desired. ! ! There are a few "special" values for NAMVR: ! ! blank, prints values of all user defined symbols. ! 'ALL', same as blank, plus values of all built in symbols. ! 'DEBUG', same as 'ALL', plus a table of built in symbols. ! ! Input, integer NINTS, the number of useful entries in INTSYM. ! ! Input, integer NSYM, the total number of symbols, including temporaries. ! ! Input, integer NSYMP, the number of permanent symbols. ! ! Input, integer NSYMS, the number of declared symbols. ! ! Input, integer NUMDIM(2,MAXSYM). ! For each symbol I, NUMDIM(1,I) is the number of rows in its ! value, and NUMDIM(2,I) is the number of columns. ! ! Input, character ( len = MAXLEN ) SYMBOL(MAXSYM), the symbolic names ! of the internally defined constants and functions. ! ! Input, real VALSYM(MAXVAL), the values of the symbolic variables. ! implicit none integer, parameter :: maxlen = 20 integer maxrpn integer maxsym integer maxval integer i integer ifree integer ifrm integer ihi integer ilo integer inc integer intsym(80) integer iopsym(maxsym) integer ipoint integer iprsym(maxsym) integer ipval(maxsym) integer irpn(maxrpn) integer itemp integer jhi integer jinc integer jlo integer khi integer klo integer lens integer maxfrm character ( len = maxlen ) namtmp character ( len = maxlen ) namvr integer ncol integer nints integer nrow integer nsym integer nsymp integer nsyms integer numdim(2,maxsym) character ( len = 100 ) output logical s_eqi character ( len = maxlen ) sym character ( len = maxlen ) symbol(maxsym) real valsym(maxval) ! ! Special output for DEBUG. ! if ( s_eqi ( namvr, 'DEBUG' ) ) then nsyms = max ( nsym, nsyms ) write ( *, '(a)' ) ' Symbol Loc Prior Args Nrow Ncol' write ( *, '(a)' ) ' ' do i = 1, nsyms write ( *, '(a20,1x,5i6)' ) symbol(i), ipval(i), iprsym(i), & iopsym(i), numdim(1,i), numdim(2,i) end do write ( *, '(a)' ) ' ' write ( *, '(a,i8)' ) 'Next free memory location at ', ifree write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'Symbol Value' do i = 1, nsyms jlo = ipval(i) jhi = jlo + numdim(1,i) * numdim(2,i)-1 jinc = 3 do klo = jlo, jhi, jinc khi = min(klo+jinc-1,jhi) write ( *, '(i4,3g20.8)' ) i, valsym(klo:khi) end do end do end if ! ! Special output for DEBUG or ALL. ! if ( s_eqi ( namvr, 'DEBUG' ) .or. s_eqi ( namvr, 'ALL' ) ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'Predefined symbols:' write ( *, '(a)' ) ' ' inc = 3 ihi = ( ( nsymp - 1 ) / inc ) + 1 do i = 1, ihi jhi = i + (inc-1) * ihi jhi = min ( jhi, nsymp ) write ( *, '(3(a20,1x))' ) symbol(i:jhi:ihi) end do end if ! ! Special output for DEBUG or ALL or ' '. ! if ( s_eqi ( namvr, 'DEBUG' ) .or. s_eqi ( namvr, 'ALL' ) .or. & namvr == ' ' ) then if ( nsymp < nsym ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'User-defined variables:' do i = nsymp+1, nsym nrow = numdim(1,i) ncol = numdim(2,i) ipoint = ipval(i) namtmp = symbol(i) call r4mat_print ( nrow, nrow, ncol, valsym(ipoint), namtmp ) end do else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'There are NO user defined variables now.' end if write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'User RPN formulas:' write ( *, '(a)' ) ' ' do ifrm = 1, maxfrm+1 if ( ifrm == maxfrm+1 ) then if ( nints <= 0 ) then cycle end if write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'Most recent user infix formula:' write ( *, '(a)' ) ' ' ilo = 1 ihi = nints else ilo = ( ifrm - 1 ) * 80 + 1 ihi = ilo + 79 end if jlo = 0 output = ' ' do i = ilo, ihi if ( ifrm == maxfrm+1 ) then itemp = intsym(i) else itemp = irpn(i) end if if ( 1 <= itemp .and. itemp <= maxsym ) then sym = symbol(itemp) lens = len_trim ( symbol(itemp) ) if ( 80 < jlo + lens ) then write ( *, '(a)' ) trim ( output ) output = ' ' jlo = 0 end if jhi = jlo + lens jlo = jlo + 1 output(jlo:jhi) = sym(1:lens) jlo = jhi + 1 if ( sym == '$' ) then exit end if end if end do write ( *, '(a)' ) trim ( output ) end do write ( *, '(a)' ) ' ' return end if ! ! List a single variable. ! do i = 1, nsym if ( s_eqi ( namvr, symbol(i) ) ) then nrow = numdim(1,i) ncol = numdim(2,i) ipoint = ipval(i) call r4mat_print ( nrow, nrow, ncol, valsym(ipoint), namvr ) return end if end do write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SYMBOL_LIST - Warning!' write ( *, '(a)' ) ' There is no variable named ' // trim ( namvr ) return end subroutine table ( ierror, infix, irpn, line, maxrow, maxrpn, & ncol, nform, nrow, value ) !*****************************************************************************80 ! !! TABLE prints a table of a formula for regularly spaced argument values. ! ! Modified: ! ! 17 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer IERROR, error flag. ! 0, no error occurred. ! nonzero, an error occurred. ! ! Input/output, character ( len = 80 ) LINE, a buffer containing ! unprocessed user input. ! ! Input, integer MAXROW, the maximum number of rows and columns. ! ! Input, integer MAXRPN, specifies the length of IRPN. ! ! Input, integer NFORM, the number of formulas that have been entered. ! ! Workspace, real VALUE(MAXROW,MAXROW), space used to hold ! the value of the variable to be saved. ! implicit none integer, parameter :: maxlen = 20 integer maxrow integer maxrpn character com integer ierror integer ifrm character ( len = * ) infix integer irpn(maxrpn) integer istep character ( len = 80 ) line character ( len = maxlen ) namtmp character ( len = maxlen ) namvar integer ncol integer nform integer nrow integer nstep character ( len = 80 ) prompt real value(maxrow,maxrow) real xhi real xlo real xval if ( nform <= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TABLE - Error!' write ( *, '(a)' ) ' Please enter a formula first!' ierror = 1 return end if ! ! Get the variable name, NAMVAR. ! namvar = ' ' prompt = 'variable, xlo, xhi, nsteps.' call chrrea ( namvar, line, prompt, ierror, 3 ) if ( ierror /= 0 ) then return end if ! ! Get the dimensions NROW and NCOL, and the value VALUE of NAMVAR. ! com = 'R' ifrm = 0 nrow = 0 ncol = 0 call comrpn ( com, ierror, ifrm, infix, irpn, maxrow, maxrpn, & namvar, ncol, nrow, value ) if ( ierror /= 0 ) then return end if if ( nrow /= 1 .or. ncol /= 1 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TABLE - Error!' write ( *, '(a)' ) ' The variable must be a scalar!' ierror = 1 return end if ! ! Get the lower and upper limits for X. ! call relrea ( xlo, line, prompt, ierror ) if ( ierror /= 0 ) then return end if call relrea ( xhi, line, prompt, ierror ) if ( ierror /= 0 ) then return end if ! ! Get the number of steps to take. ! call intrea ( ierror, nstep, line, prompt ) if ( ierror /= 0 ) then return end if if ( nstep <= 0 ) then ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TABLE - Error!' write ( *, '(a)' ) ' The number of steps must be positive!' return end if write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'Table of X, F(X)' write ( *, '(a)' ) ' ' do istep = 1, nstep + 1 com = 'V' xval = ( real ( nstep + 1 - istep ) * xlo + real ( istep - 1 ) * xhi ) & / real ( nstep ) value(1,1) = xval call comrpn ( com, ierror, ifrm, infix, irpn, maxrow, maxrpn, & namvar, ncol, nrow, value ) if ( ierror /= 0 ) then return end if ! ! Evaluate the formula at XVAL. ! com = 'E' ifrm = 1 call comrpn ( com, ierror, ifrm, infix, irpn, maxrow, maxrpn, & namvar, ncol, nrow, value ) if ( ierror /= 0 ) then return end if ! ! Print out the value of XVAL versus the formula F(XVAL). ! if ( nrow == 1 .and. ncol == 1 ) then write ( *, '(2g14.6)' ) xval, value(1,1) else write ( *, '(a)' ) ' ' write ( *, '(g14.6)' ) xval namtmp = ' ' call r4mat_print ( maxrow, nrow, ncol, value, namtmp ) end if end do return end subroutine timestamp ( ) !*****************************************************************************80 ! !! TIMESTAMP prints the current YMDHMS date as a time stamp. ! ! Example: ! ! May 31 2001 9:45:54.872 AM ! ! Modified: ! ! 31 May 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! None ! implicit none character ( len = 8 ) ampm integer d character ( len = 8 ) date integer h integer m integer mm character ( len = 9 ), parameter, dimension(12) :: month = (/ & 'January ', 'February ', 'March ', 'April ', & 'May ', 'June ', 'July ', 'August ', & 'September', 'October ', 'November ', 'December ' /) integer n integer s character ( len = 10 ) time integer values(8) integer y character ( len = 5 ) zone call date_and_time ( date, time, zone, values ) y = values(1) m = values(2) d = values(3) h = values(5) n = values(6) s = values(7) mm = values(8) if ( h < 12 ) then ampm = 'AM' else if ( h == 12 ) then if ( n == 0 .and. s == 0 ) then ampm = 'Noon' else ampm = 'PM' end if else h = h - 12 if ( h < 12 ) then ampm = 'PM' else if ( h == 12 ) then if ( n == 0 .and. s == 0 ) then ampm = 'Midnight' else ampm = 'AM' end if end if end if write ( *, '(a,1x,i2,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) & trim ( month(m) ), d, y, h, ':', n, ':', s, '.', mm, trim ( ampm ) return end subroutine tokens ( ierror, ifinis, ifree, implic, indx1, indx2, ineg, infix, & intsym, iopsym, iprsym, ipval, maxrow, maxsym, maxval, nints, nsym, nsymp, & numdim, nuvar, symbol, valsym ) !*****************************************************************************80 ! !! TOKENS parses a character string into recognized symbols and constants. ! ! Discussion: ! ! The routine produces the output array of integers INTSYM, containing ! whose entries stand for recognized symbols, variable names or ! constants. ! ! Modified: ! ! 02 April 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer IERROR, error flag. ! 0, no error occurred. ! nonzero, an error occurred. ! ! Input, integer IFINIS, the index in SYMBOL of the "end" symbol '$'. ! ! Input, integer IFREE, the next free memory location in VALSYM. ! ! Output, integer IMPLIC, implicit variable flag. ! 0, the assignment variable was not implicit. ! nonzero, the assignment variable was implicit, and has been ! added to the symbol table, and assigned the symbol index of IMPLIC. ! ! Input, integer INDX1, the index of the symbol 'INDEX1'. ! ! Input, integer INDX2, the index of the symbol 'INDEX2'. ! ! Input, integer INEG, the index of the symbol 'NEG'. ! ! Input, character ( len = * ) INFIX, an infix formula. ! ! Output, integer INTSYM(80), a sequence of integers which are ! the indices of the symbols used in INFIX. ! ! Input, integer IOPSYM(MAXSYM), the number of operands per symbol. ! Constants have 0, unaries have 1, and so on. ! ! Input, integer IPRSYM(MAXSYM), the priority of each symbol. ! ! Input, integer IPVAL(MAXSYM), contains, for each symbol, ! the address in VALSYM where the value of the symbol ! is stored, if it is a scalar. If the symbol represents ! a vector or matrix, IPVAL(IARG) points to location of ! the first entry. ! ! Input, integer MAXROW, the maximum number of rows and columns. ! ! Input, integer MAXSYM, the maximum number of symbols allowed. ! ! Input, integer MAXVAL, the maximum number of values allowed. ! ! Output, integer NINTS, the number of useful entries in INTSYM. ! ! Input, integer NSYM, the total number of symbols, including temporaries. ! ! Input, integer NSYMP, the number of permanent symbols. ! ! Input, integer NUMDIM(2,MAXSYM). ! For each symbol I, NUMDIM(1,I) is the number of rows in its ! value, and NUMDIM(2,I) is the number of columns. ! ! Input, integer NUVAR. ! If NUVAR is 0, then the formula may refer to variables which ! have not yet been declared. ! Otherwise, all variables on the right hand side of the "=" ! sign must have been previously declared. ! ! Input, character ( len = MAXLEN ) SYMBOL(MAXSYM), the symbolic names ! of the internally defined constants and functions. ! ! Input, real VALSYM(MAXVAL), the values of the symbolic variables. ! implicit none integer, parameter :: maxlen = 20 integer maxsym character ctemp integer i integer iback integer icom1 integer icom2 integer ierror integer ifinis integer ifree integer ilpr integer implic integer indx1 integer indx2 integer ineg character ( len = * ) infix integer intsym(80) integer iopsym(maxsym) integer ipoint integer ipos integer iprsym(maxsym) integer ipval(maxsym) integer irpr integer j integer jback integer lchar integer lenfix integer lenmat integer lens integer loc integer loc1 integer loc2 integer match integer matcho integer maxrow integer maxval character ( len = maxlen ) namvr integer ncol integer nints logical notnum integer nrow integer nsym integer nsymp integer numdim(2,maxsym) integer nuvar real rval logical s_eqi integer s_indexi logical s_is_alpha character sym character ( len = maxlen ) sym1 character ( len = maxlen ) sym2 character ( len = maxlen ) sym3 character ( len = maxlen ) symbol(maxsym) real valsym(maxval) implic = 0 ierror = 0 nints = 0 match = 0 ipos = 1 lenfix = len_trim ( infix ) ! ! 1: Does this formula begin with a variable name followed by an equals sign? ! sym = '=' loc1 = s_indexi ( infix, sym ) if ( 1 < loc1 ) then ! ! The formula begins with a variable name followed by an equals sign. ! ! Pick off the name of the variable, but first make sure ! you stop before any left hand parenthesis, in case the ! formula is something like "alfred(7)=5". ! sym = '(' loc2 = s_indexi ( infix, sym ) if ( loc2 <= 0 ) then loc = loc1 else loc = min ( loc1, loc2 ) end if if ( 1 < loc .and. loc <= maxlen+1 ) then loc = loc - 1 namvr = infix(1:loc) ! ! Has the assignment variable name already been defined? ! match = 0 do i = 1, nsym if ( s_eqi ( namvr, symbol(i) ) ) then match = i exit end if end do ! ! The assignment variable has not previously been defined. ! Prepare to add it to the table of symbols. ! if ( match == 0 ) then ipos = loc + 1 nrow = maxrow ncol = maxrow write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TOKENS:' write ( *, '(a)' ) ' Implicit definition of ', trim ( namvr ) call symbol_add ( ierror, ifree, iopsym, iprsym, ipval, maxrow, & maxsym, maxval, namvr, ncol, nrow, nsym, nsymp, numdim, symbol, & valsym ) if ( ierror /= 0 ) then return end if nints = nints + 1 intsym(nints) = nsym match = nsym implic = nsym end if end if end if ! ! Process the next chunk of the formula. ! do ! ! If we have reached the end of the input string, we're done. ! Append an "End-of-formula" marker "$" and return. ! if ( lenfix < ipos ) then nints = nints + 1 intsym(nints) = ifinis exit end if ! ! We're not done yet. Remember the last match in MATCHO, and ! look for the next match. Go through all the symbols, and ! try to find the longest match possible. ! matcho = match match = 0 lenmat = 0 do i = 1, nsym lens = len_trim ( symbol(i) ) if ( lenmat < lens .and. ipos + lens - 1 <= lenfix ) then if ( s_eqi ( infix(ipos:ipos+lens-1), symbol(i) ) ) then lenmat = lens match = i end if end if end do if ( match == 0 ) then go to 110 end if ! ! We found a match. ! ! But watch out! The user might be implicitly defining a name, ! but because the first part matches an earlier symbol, you don't ! realize it! For example, in a formula involving SINGLE, you ! might think you see "SIN". So make sure that the match you've ! got is not followed by an alphabetic character. ! if ( iopsym(match) == 0 ) then if ( ipos+lenmat-1 < lenfix ) then ctemp = infix(ipos+lenmat:ipos+lenmat) if ( s_is_alpha ( ctemp ) ) then go to 140 end if end if end if ! ! We have matched an old symbol. ! sym1 = symbol(match) if ( matcho == 0 ) then sym2 = ' ' else sym2 = symbol(matcho) end if ! ! Check for unary minus or plus. ! if ( sym1 == '-' ) then if ( sym2 == '**' .or. sym2 == ',' .or. & sym2 == '(' .or. sym2 == '=' .or. & ipos == 1 ) then sym = infix(ipos+1:ipos+1) if ( lge ( sym, '0' ) .and. lle ( sym, '9' ) ) then go to 110 end if if ( sym == '.' ) then go to 110 end if end if end if if ( sym1 == '-' ) then if ( ipos == 1 .or. sym2 == '(' .or. sym2 == ',' .or. sym2 == '=' ) then nints = nints + 1 intsym(nints) = ineg ipos = ipos + 1 cycle end if end if if ( sym1 == '+' ) then if ( ipos == 1 .or. sym2 == '(' .or. sym2 == '=' ) then ipos = ipos + 1 cycle end if end if ! ! If we've hit a right parenthesis, then we assume we're dealing with: ! ! a vector index: A(3) ! a matrix index: A(1,2) ! a two place operator: MAX(A,B), or MIN or POLY ! a three place operator: KRYLOV(A,X,N) ! ! We need to find the matching left parenthesis, count the ! arguments by assuming that commas separate them, and read ! off the name preceding the left parenthesis. ! if ( sym1 == ')' ) then icom1 = 0 icom2 = 0 ilpr = 0 irpr = 1 do iback = 1, nints i = nints + 1 - iback if ( intsym(i) <= 0 ) then go to 90 end if sym3 = symbol ( intsym(i) ) if ( sym3 == ',' .and. irpr-ilpr == 1 ) then if ( icom1 == 0 ) then icom1 = i else if ( icom2 == 0 ) then icom2 = i else ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TOKENS - Error!' write ( *, '(a)' ) ' Too many commas in formula.' return end if else if ( sym3 == '(' ) then ilpr = ilpr + 1 else if ( sym3 == ')' ) then irpr = irpr + 1 end if ! ! If ILPR=IRPR, we've reached the matching left parenthesis. ! if ( ilpr == irpr ) then if ( i <= 1 ) then go to 100 end if if ( intsym(i-1) <= 0 ) then go to 100 end if sym3 = symbol ( intsym(i-1) ) if ( iopsym(intsym(i-1)) /= 0 .and. & .not. s_eqi ( sym3, 'MAX' ) .and. & .not. s_eqi ( sym3, 'MIN' ) .and. & .not. s_eqi ( sym3, 'ATAN2' ) .and. & .not. s_eqi ( sym3, 'POLY' ) .and. & .not. s_eqi ( sym3, 'KRYLOV' ) ) then go to 100 end if ! ! Copy ) into INTSYM. ! nints = nints + 1 intsym(nints) = match ! ! INDX1 is our name for a vector reference. ! if ( icom1 == 0 ) then match = indx1 ! ! INDX2 is our name for a matrix reference. ! else if ( icom2 == 0 ) then intsym(icom1) = match do jback = 1, nints-icom1 j = nints + 1 - jback intsym(j+1) = intsym(j) end do nints = nints + 1 intsym(icom1+1) = match - 1 if ( iopsym(intsym(i-1)) /= 0 ) then ipos = ipos + 1 go to 11 end if match = indx2 ! ! Insert )( to divide a trio of arguments. ! else intsym(icom1) = match do jback = 1, nints-icom1 j = nints + 1 - jback intsym(j+1) = intsym(j) end do intsym(icom1+1) = match - 1 nints = nints + 1 intsym(icom2) = match do jback = 1, nints-icom2 j = nints + 1 - jback intsym(j+1) = intsym(j) end do intsym(icom2+1) = match - 1 nints = nints + 1 if ( iopsym(intsym(i-1)) /= 0 ) then ipos = ipos + 1 go to 11 end if match = indx2 end if ipos = ipos + 1 nints = nints + 1 intsym(nints) = match go to 11 end if 90 continue end do end if 100 continue lens = len_trim ( symbol(match) ) nints = nints + 1 intsym(nints) = match ipos = ipos + lens cycle ! ! Check for a constant. ! 110 continue call s_to_r4 ( infix(ipos:), rval, ierror, lchar ) if ( lchar <= 0 ) then go to 140 end if if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TOKENS - Error!' write ( *, '(a)' ) ' Illegal real number.' exit end if ! ! Look for a constant already defined with the same value. ! do i = nsymp+1, nsym if ( iprsym(i) == 10 ) then if ( symbol(i)(1:1) == '"' ) then ipoint = ipval(i) if ( valsym(ipoint) == rval ) then match = i go to 130 end if end if end if end do ! ! This is a new constant. Make up a name for it. ! if ( real ( int ( rval ) ) == rval ) then write ( namvr, '(''"'',i19)' ) int ( rval ) else write ( namvr, '(''"'', g19.10)' ) rval end if call s_blank_delete ( namvr ) ! ! Add the constant to the symbol list. ! nrow = 1 ncol = 1 call symbol_add ( ierror, ifree, iopsym, iprsym, ipval, maxrow, maxsym, & maxval, namvr, ncol, nrow, nsym, nsymp, numdim, symbol, valsym ) if ( ierror /= 0 ) then exit end if ipoint = ipval(nsym) valsym(ipoint) = rval match = nsym 130 continue ! ! Advance LCHAR characters, except that CHRCTR will read in a ! comma as the end of the number, and we have to back up in ! such a case. ! if ( infix(ipos+lchar-1:ipos+lchar-1) /= ',' ) then ipos = ipos + lchar else ipos = ipos + lchar - 1 end if nints = nints + 1 intsym(nints) = match cycle ! ! Consider implicit variable declaration on right hand side. ! 140 continue if ( nuvar == 0 ) then ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TOKENS - Error!' write ( *, '(a)' ) ' Undeclared variable in formula.' exit end if namvr = ' ' do i = 1, maxlen sym = infix(ipos-1+i:ipos-1+i) notnum = llt ( sym, '0' ) .or. lgt ( sym, '0' ) if ( ( .not. s_is_alpha(sym) ) .and. notnum ) then exit end if namvr(i:i) = sym end do if ( len_trim ( namvr ) <= 0 ) then ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TOKENS - Error!' write ( *, '(a)' ) ' The string starting at location ', ipos, & ' is unreadable.' exit end if ipos = ipos + lchar nrow = 1 ncol = 1 call symbol_add ( ierror, ifree, iopsym, iprsym, ipval, maxrow, maxsym, & maxval, namvr, ncol, nrow, nsym, nsymp, numdim, symbol, valsym ) nints = nints + 1 intsym(nints) = nsym match = nsym numdim(1,nsym) = 1 numdim(2,nsym) = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TOKENS:' write ( *, '(a)' ) ' The undeclared variable ', trim ( namvr ), & ' will be assumed to be a scalar.' 11 continue end do return end