program matman ! !******************************************************************************* ! !! MATMAN is a program for interactive linear algebra demonstrations. ! ! ! Discussion: ! ! This is version 1.63. ! ! Journal: ! ! 09 October 2000: ! CHANGE calls VALUE_READ so I can add a column easily too. ! COL_APP allows me to append a column. ! ! 29 September 2000: ! Adding IDENTITY operation which appends I. ! Added modular arithmetic option. ! Added RANDOM. ! ! 26 September 2000: ! Might almost be working. Now either implement ! modular arithmetic, or saving of pre and post multipliers. ! ! 25 September 2000: ! Add integer arithmetic option, decide what to do about division later. ! ! 24 September 2000: ! Throw out superfluous items! ! ! Modified: ! ! 09 October 2000 ! ! Author: ! ! John Burkardt ! implicit none ! integer, parameter :: maxcol = 30 integer, parameter :: maxrow = 16 ! real a(maxrow,maxcol) integer a_int(maxrow,maxcol) integer base logical ch_is_digit character ( len = 20 ) command integer i_temp integer iabot(maxrow,maxcol) integer iatop(maxrow,maxcol) integer icol integer icol1 integer icol2 integer ictop(maxrow,maxcol) integer idbot integer idtop integer ierror integer iform integer ihi integer ilo integer iprint integer irow integer irow1 integer irow2 character isay integer isbot integer istop integer iterm integer jhi integer jform integer jlo character ( len = 80 ) line character ( len = 80 ) line2 integer ncol integer nrow character ( len = 80 ) prompt logical s_eqi real sval ! call timestamp ( ) ! ! Initializations. ! call init ( a, a_int, base, iabot, iatop, ierror, iform, & line, maxrow, maxcol, nrow, ncol ) ! ! Say hello ! call hello ! ! Get the next command from the user. ! do iprint = 0 if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'MATMAN could not carry out your command:' write ( *, '(a)' ) '"' // trim ( command ) // '"' ierror = 0 end if write ( *, '(a)' ) ' ' line = ' ' prompt = 'command ( H for help )' ! ! No check for terminators. ! iterm = 0 call s_read ( line2, line, prompt, ierror, iterm ) line = line2 if ( line2 == ' ' ) then cycle end if if ( ierror /= 0 ) then ierror = 0 command = 'Q' end if ! ! Check to see if the command is an ERO, in a special format. ! call s_blank_delete ( line2 ) if ( s_eqi ( line2, 'R' ) ) then else if ( s_eqi ( line2(1:1), 'R' ) .and. & ( line2(2:2) == ' ' .or. ch_is_digit ( line2(2:2) ) ) ) then call row_op_check ( command, ierror, line2 ) if ( ierror /= 0 ) then cycle end if line = line2 else if ( s_eqi ( line2(1:3), 'ROW' ) .and. & ( line2(4:4) == ' ' .or. ch_is_digit ( line2(4:4) ) ) ) then call row_op_check ( command, ierror, line2 ) if ( ierror /= 0 ) then cycle end if line = line2 ! ! Check to see if the command is an ECO, in a special format. ! else if ( s_eqi ( line2, 'C' ) ) then else if ( s_eqi ( line2(1:1), 'C' ) .and. & ( line2(2:2) == ' ' .or. ch_is_digit ( line2(2:2) ) ) ) then call col_op_check ( command, ierror, line2 ) if ( ierror /= 0 ) then cycle end if line = line2 else if ( s_eqi ( line2(1:3), 'COL' ) .and. & ( line2(4:4) == ' ' .or. ch_is_digit ( line2(4:4) ) ) ) then call col_op_check ( command, ierror, line2 ) if ( ierror /= 0 ) then cycle end if line = line2 else if ( s_eqi ( line2(1:6), 'COLUMN' ) .and. & ( line2(7:7) == ' ' .or. ch_is_digit ( line2(7:7) ) ) ) then call col_op_check ( command, ierror, line2 ) if ( ierror /= 0 ) then cycle end if line = line2 ! ! Maybe this is all I need in order to get my change command. ! else if ( s_eqi(line2(1:2), 'A(' ) ) then command = 'A(' line = line2(3:) else command = ' ' end if ! ! If the command was not an ERO that had to be translated, ! read it the regular way. ! ! Blank, slash, comma, semicolon, equals terminate the command. ! The rest of the user input stays in LINE. ! if ( command == ' ' ) then iterm = 1 call s_read ( command, line, prompt, ierror, iterm ) if ( ierror /= 0 ) then ierror = 0 command = 'Q' end if end if ! ! A(I,J)=X ! if ( s_eqi ( command(1:2), 'A(' ) ) then call change ( a, a_int, iatop, iabot, base, ierror, iform, line, & nrow, ncol ) iprint = 1 ! ! B=Set up sample problem. ! else if ( s_eqi ( command, 'B' ) ) then call size_set ( ierror, line, maxrow, nrow, 'rows' ) call size_set ( ierror, line, maxcol, ncol, 'columns' ) if ( s_eqi ( isay,'S' ) ) then if ( ncol + 1 > maxcol ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Your problem is too big!' cycle end if ncol = ncol + 1 call problem_solve_random ( a, a_int, iatop, iabot, base, iform, & nrow, ncol ) end if ! ! BASE = set base for modular arithmetic. ! else if ( s_eqi ( command, 'BASE' ) ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'Enter a positive integer to use as the base' write ( *, '(a)' ) 'for modular arithmetic.' read ( *, * ) base jform = 4 call form ( a, a_int, iatop, iabot, base, iform, jform, nrow, ncol ) iprint = 1 ! ! CHECK = Echelon Form. ! else if ( s_eqi ( command, 'CHECK' ) ) then call la_opt ( a, a_int, iabot, iatop, base, iform, nrow, ncol ) ! ! COL_ADD=Add a multiple of one column to another. ! else if ( s_eqi ( command, 'COL_ADD' ) ) then call col_add_param ( ierror, base, iform, icol1, icol2, istop, isbot, & line, ncol, sval ) if ( ierror /= 0 ) then cycle end if call col_add ( a, a_int, iatop, iabot, ierror, base, iform, icol1, icol2, & nrow, ncol, sval, istop, isbot ) iprint = 1 ! ! COL_APP = Append column ! else if ( s_eqi ( command(1:7), 'COL_APP' ) ) then ncol = ncol + 1 call col_app ( a, a_int, iatop, iabot, base, ierror, iform, line, & nrow, ncol ) if ( ierror /= 0 ) then ncol = ncol - 1 else iprint = 1 end if ! ! COL_AUTO=Automatic column reduction. ! else if ( s_eqi ( command, 'COL_AUTO' ) ) then call col_auto ( a, a_int, iatop, iabot, ierror, base, iform, nrow, ncol ) iprint = 1 ! ! COL_DIV = Divide column by scalar. ! else if ( s_eqi ( command, 'COL_DIV' ) ) then call col_div_param ( icol, ierror, base, iform, isbot, istop, line, sval ) call col_div ( a, a_int, iatop, iabot, ierror, base, iform, icol, & nrow, ncol, sval, istop, isbot ) if ( ierror == 0 ) then iprint = 1 end if ! ! COL_MUL=Multiply column by scalar. ! else if ( s_eqi ( command, 'COL_MUL' ) ) then call col_mul_param ( ierror, base, iform, icol, istop, isbot, line, sval ) if ( ierror /= 0 ) then cycle end if call col_mul ( a, a_int, iatop, iabot, ierror, base, iform, icol, & nrow, ncol, sval, istop, isbot ) iprint = 1 ! ! COL_SWAP = Swap columns I and J. ! else if ( s_eqi ( command, 'COL_SWAP' ) ) then prompt = 'column I, column J.' call i_read ( icol1, line, prompt, ierror ) if ( ierror /= 0 ) then cycle end if call i_read ( icol2, line, prompt, ierror ) if ( ierror /= 0 ) then cycle end if call col_swap ( a, a_int, iatop, iabot, ierror, iform, icol1, icol2, & nrow, ncol ) iprint = 1 ! ! DEC_DIGIT=Set number of digits. ! else if ( s_eqi ( command, 'DEC_DIGIT' ) ) then call dec_digit_set ( ierror, line ) ! ! DECimal = use decimal arithmetic ! else if ( s_eqi ( command(1:3), 'DEC' ) ) then jform = 2 call form ( a, a_int, iatop, iabot, base, iform, jform, nrow, ncol ) iprint = 1 ! ! E=Enter problem definition ! else if ( s_eqi ( command, 'E' ) ) then call size_set ( ierror, line, maxrow, nrow, 'rows' ) call size_set ( ierror, line, maxcol, ncol, 'columns' ) call la_inp1 ( nrow, ncol, a, a_int, iabot, iatop, ierror, base, iform, & line, 1, nrow, 1, ncol ) if ( ierror /= 0 ) then cycle end if iprint = 1 ! ! H=Help. ! else if ( s_eqi ( command, 'H' ) ) then call help ! ! I_BIG = Set maximum integer. ! else if ( s_eqi ( command, 'I_BIG' ) ) then prompt = 'maximum integer for rational representations.' call i_read ( i_temp, line, prompt, ierror ) call i_data ( 'SET', 'I_BIG', i_temp ) ! ! IDENTITY ! else if ( s_eqi ( command(1:2), 'ID' ) ) then if ( nrow <= 0 ) then call size_set ( ierror, line, maxrow, nrow, 'rows' ) ncol = 0 if ( ierror /= 0 ) then cycle end if end if if ( ncol + nrow > maxcol ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Not enough room to append the identity matrix.' cycle end if call mat_append_identity ( nrow, ncol, a, a_int, iatop, iabot, iform ) iprint = 1 ! ! INIT = Initializations. ! else if ( s_eqi ( command(1:3), 'INIT' ) ) then call init ( a, a_int, base, iabot, iatop, ierror, iform, & line, maxrow, maxcol, nrow, ncol ) ! ! INT = Integer arithmetic ! else if ( s_eqi ( command(1:3), 'INT' ) ) then jform = 3 call form ( a, a_int, iatop, iabot, base, iform, jform, nrow, ncol ) iprint = 1 ! ! Q = Quit (after confirmation). ! QUit = QUIT NOW! ! else if ( s_eqi ( command(1:1), 'Q' ) ) then if ( s_eqi ( command(2:2), 'U' ) ) then exit end if line = ' ' prompt = '"Y" to confirm you want to quit.' iterm = 0 call s_read ( isay, line, prompt, ierror, iterm ) if ( ierror /= 0 ) then exit end if if ( s_eqi ( isay, 'Y' ) ) then exit end if ! ! RANDOM = Enter problem definition ! else if ( s_eqi ( command(1:3), 'RAN' ) ) then call size_set ( ierror, line, maxrow, nrow, 'rows' ) call size_set ( ierror, line, maxcol, ncol, 'columns' ) call mat_random ( nrow, ncol, base, iform, a, a_int, iatop, iabot ) iprint = 1 ! ! RATional = use rational arithmetic ! else if ( s_eqi ( command(1:3), 'RAT' ) ) then jform = 0 call form ( a, a_int, iatop, iabot, base, iform, jform, nrow, ncol ) iprint = 1 ! ! REAl = use real arithmetic ! else if ( s_eqi ( command(1:3), 'REA' ) ) then jform = 1 call form ( a, a_int, iatop, iabot, base, iform, jform, nrow, ncol ) iprint = 1 ! ! ROW_ADD=Add a multiple of one row to another. ! else if ( s_eqi ( command, 'ROW_ADD' ) ) then call row_add_param ( ierror, base, iform, irow1, irow2, istop, isbot, & line, nrow, sval ) if ( ierror /= 0 ) then cycle end if call row_add ( a, a_int, iatop, iabot, ierror, base, iform, irow1, irow2, & nrow, ncol, sval, istop, isbot ) iprint = 1 ! ! ROW_AUTO=Automatic row reduction. ! else if ( s_eqi ( command, 'ROW_AUTO' ) ) then call row_auto ( a, a_int, iatop, iabot, ierror, base, iform, nrow, ncol ) iprint = 1 ! ! ROW_DIV = Divide row by scalar. ! else if ( s_eqi ( command, 'ROW_DIV' ) ) then call row_div_param ( ierror, base, iform, irow, isbot, istop, line, sval ) call row_div ( a, a_int, iatop, iabot, ierror, base, iform, irow, & nrow, ncol, sval, istop, isbot ) if ( ierror == 0 ) then iprint = 1 end if ! ! ROW_MUL=Multiply row by scalar. ! else if ( s_eqi ( command, 'ROW_MUL' ) ) then call row_mul_param ( ierror, base, iform, irow, istop, isbot, line, sval ) if ( ierror /= 0 ) then cycle end if call row_mul ( a, a_int, iatop, iabot, ierror, base, iform, irow, & nrow, ncol, sval, istop, isbot ) iprint = 1 ! ! ROW_SWAP = Interchange rows I and J. ! else if ( s_eqi ( command, 'ROW_SWAP' ) ) then prompt = 'row I, row J.' call i_read ( irow1, line, prompt, ierror ) if ( ierror /= 0 ) then cycle end if call i_read ( irow2, line, prompt, ierror ) if ( ierror /= 0 ) then cycle end if call row_swap ( a, a_int, iatop, iabot, ierror, iform, & irow1, irow2, nrow, ncol ) iprint = 1 ! ! Otherwise, type out the current matrix ! else iprint = 1 end if ! ! After certain operations, print out the matrix. ! if ( ierror == 0 .and. iprint == 1 ) then call mat_print ( a, a_int, iabot, iatop, iform, nrow, ncol, & 'The current matrix:' ) iprint = 0 end if end do write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'MATMAN' write ( *, '(a)' ) ' Normal end of execution.' stop end subroutine ch_cap ( c ) ! !******************************************************************************* ! !! 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 ) ! !******************************************************************************* ! !! CH_EQI is a case insensitive comparison of two characters for equality. ! ! ! Examples: ! ! CH_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 ) ! !******************************************************************************* ! !! 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 function ch_is_digit ( c ) ! !******************************************************************************* ! !! CH_IS_DIGIT returns .TRUE. if a character is a decimal digit. ! ! ! Modified: ! ! 09 August 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character C, the character to be analyzed. ! ! Output, logical CH_IS_DIGIT, .TRUE. if C is a digit, .FALSE. otherwise. ! implicit none ! character c logical ch_is_digit ! if ( lge ( c, '0' ) .and. lle ( c, '9' ) ) then ch_is_digit = .true. else ch_is_digit = .false. end if return end subroutine ch_to_digit ( c, digit ) ! !******************************************************************************* ! !! 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 change ( a, a_int, iatop, iabot, base, ierror, iform, line, & nrow, ncol ) ! !******************************************************************************* ! !! CHANGE allows the user to change an entry in the array. ! ! ! Discussion: ! ! Expect an input line of the form: ! ! A(5,12) = 17 ! ! Modified: ! ! 12 July 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, real A(NROW,NCOL), the current matrix. ! ! Input/output, integer IATOP(NROW,NCOL), IABOT(NROW,NCOL). ! IATOP and IABOT represent the rational or decimal matrix ! whose entry is to be changed. ! ! Output, integer IERROR, 0 for no error, 1 for an error. ! ! Input, integer IFORM, 0/1/2/3/4, rat/real/dec/int/mod arithmetic. ! ! Workspace, character ( len = 80 ) LINE, used to hold the user's input. ! ! Input, integer NROW, NCOL, the number of rows and columns in the matrix. ! implicit none ! integer ncol integer nrow ! real a(nrow,ncol) integer a_int(nrow,ncol) integer base character ( len = 10 ) chrtmp1 character ( len = 10 ) chrtmp2 character ( len = 22 ) chrtmp3 integer i_gcd integer iabot(nrow,ncol) integer iatop(nrow,ncol) integer icol integer ierror integer iform integer irow integer isbot integer istop integer itemp integer ival character ( len = 80 ) line character ( len = 80 ) prompt real rval ! prompt = 'row I, column J, new value S.' ! ! Get the row number. ! call i_read ( irow, line, prompt, ierror ) if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'CHANGE - Error!' write ( *, '(a)' ) ' I_READ returned error flag!' return end if if ( irow < 1 .or. irow > nrow ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'CHANGE - Error!' write ( *, '(a)' ) ' Illegal row value!' ierror = 1 return end if ! ! Get the column number. ! if ( line(1:1) == ',' ) then line(1:1) = ' ' end if call i_read ( icol, line, prompt, ierror ) if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'CHANGE - Error!' write ( *, '(a)' ) ' I_READ returned error flag!' return end if if ( icol < 1 .or. icol > ncol ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'CHANGE - Error!' write ( *, '(a)' ) ' Illegal column value!' ierror = 1 return end if ! ! Read the value. ! if ( line(1:2) == ')=' ) then line(1:2) = ' ' end if call value_read ( irow, icol, iform, base, nrow, ncol, line, a, a_int, & iatop, iabot, ierror ) call i_to_s_left ( irow, chrtmp1 ) call i_to_s_left ( icol, chrtmp2 ) if ( iform == 0 ) then call rat_to_s_left ( iatop(irow,icol), iabot(irow,icol), chrtmp3 ) else if ( iform == 1 ) then call r_to_s_left ( a(irow,icol), chrtmp3 ) else if ( iform == 2 ) then call dec_to_s_left ( iatop(irow,icol), iabot(irow,icol), chrtmp3 ) else if ( iform == 3 ) then call i_to_s_left ( a_int(irow,icol), chrtmp3 ) else if ( iform == 4 ) then call i_to_s_left ( a_int(irow,icol), chrtmp3 ) end if write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' A(' // trim ( chrtmp1 ) // ',' // trim ( chrtmp2 ) // & ') = ' // trim ( chrtmp3 ) return end subroutine chrctf ( s, itop, ibot, ierror, lchar ) ! !******************************************************************************* ! !! CHRCTF reads an integer or rational fraction from a string. ! ! ! Discussion: ! ! The integer may be in real format, for example '2.25'. The routine ! returns ITOP and IBOT. If the input number is an integer, ITOP ! equals that integer, and IBOT is 1. But in the case of 2.25, ! the program would return ITOP = 225, IBOT = 100. ! ! Legal input is: ! ! blanks, ! initial sign, ! blanks, ! integer part, ! decimal point, ! fraction part, ! 'E' or 'e' or 'D' or 'd', exponent marker, ! exponent sign, ! exponent integer part, ! blanks, ! final comma or semicolon. ! ! with most quantities optional. ! ! Examples: ! ! S ITOP IBOT ! ! '1' 1 1 ! ' 1 ' 1 1 ! '1A' 1 1 ! '12,34,56' 12 1 ! ' 34 7' 34 1 ! '-1E2ABCD' -100 1 ! '-1X2ABCD' -1 1 ! ' 2E-1' 2 10 ! '23.45' 2345 100 ! ! Modified: ! ! 07 December 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) S, the string containing the ! data to be read. Reading will begin at position 1 and ! terminate when no more characters ! can be read to form a legal integer. Blanks, commas, ! or other nonnumeric data will, in particular, cause ! the conversion to halt. ! ! Output, integer ITOP, the integer read from the string, ! assuming that no negative exponents or fractional parts ! were used. Otherwise, the 'integer' is ITOP/IBOT. ! ! Output, integer IBOT, the integer divisor required to ! represent numbers which are in real format or have a ! negative exponent. ! ! Output, integer IERROR, error flag. ! 0 if no errors, ! Value of IHAVE when error occurred otherwise. ! ! Output, integer LCHAR, number of characters read from ! the string to form the number. ! implicit none ! logical ch_eqi character c integer ibot integer ierror integer ihave integer isgn integer iterm integer itop integer jsgn integer jtop integer lchar integer nchar integer ndig character ( len = * ) s ! nchar = len_trim ( s ) ierror = 0 lchar = - 1 isgn = 1 itop = 0 ibot = 1 jsgn = 1 jtop = 0 ihave = 1 iterm = 0 do while ( lchar < nchar ) lchar = lchar + 1 c = s(lchar+1:lchar+1) ! ! Blank. ! if ( c == ' ' ) then if ( ihave == 2 ) then else if ( ihave == 6 .or. ihave == 7 ) then iterm = 1 else if ( ihave > 1 ) 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 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 ( lge ( c, '0' ) .and. lle ( c, '9' ) .and. ihave < 11 ) then if ( ihave <= 2 ) then ihave = 3 else if ( ihave == 4 ) then ihave = 5 else if ( ihave == 6 .or. ihave == 7 ) then ihave = 8 end if call ch_to_digit ( c, ndig ) if ( ihave == 3 ) then itop = 10 * itop + ndig else if ( ihave == 5 ) then itop = 10 * itop + ndig ibot = 10 * ibot else if ( ihave == 8 ) then jtop = 10 * jtop + ndig end if ! ! Anything else is regarded as a terminator. ! else iterm = 1 end if if ( iterm == 1 ) then exit end if end do if ( iterm /= 1 .and. lchar+1 == nchar ) then lchar = nchar end if ! ! Number seems to have terminated. Have we got a legal number? ! if ( ihave == 1 .or. ihave == 2 .or. ihave == 6 .or. ihave == 7 ) then ierror = ihave write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'CHRCTF - Serious error!' write ( *, '(a,a)' ) ' Illegal input:', trim ( s ) return end if ! ! Number seems OK. Form it. ! if ( jsgn == 1 ) then itop = itop * 10**jtop else ibot = ibot * 10**jtop end if itop = isgn * itop return end subroutine chrctg ( string, itop, ibot, ierror, lchar ) ! !******************************************************************************* ! !! CHRCTG reads an integer, decimal fraction or a ratio from a string. ! ! ! Discussion: ! ! CHRCTG returns an equivalent ratio (ITOP/IBOT). ! ! If the input number is an integer, ITOP equals that integer, and ! IBOT is 1. But in the case of 2.25, the program would return ! ITOP = 225, IBOT = 100. ! ! A ratio is either ! a number ! or ! a number, "/", a number. ! ! A "number" is defined as: ! ! blanks, ! initial sign, ! integer part, ! decimal point, ! fraction part, ! E, ! exponent sign, ! exponent integer part, ! blanks, ! final comma or semicolon, ! ! with most quantities optional. ! ! Examples: ! ! STRING ITOP IBOT ! ! '1' 1 1 ! ' 1 ' 1 1 ! '1A' 1 1 ! '12,34,56' 12 1 ! ' 34 7' 34 1 ! '-1E2ABCD' -100 1 ! '-1X2ABCD' -1 1 ! ' 2E-1' 2 10 ! '23.45' 2345 100 ! ! Modified: ! ! 19 July 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) STRING, the string containing the ! data to be read. Reading will begin at position 1 and ! terminate when no more characters ! can be read to form a legal integer. Blanks, commas, ! or other nonnumeric data will, in particular, cause ! the conversion to halt. ! ! Output, integer ITOP, the integer read from the string, ! assuming that no negative exponents or fractional parts ! were used. Otherwise, the 'integer' is ITOP/IBOT. ! ! Output, integer IBOT, the integer divisor required to ! represent numbers which are in decimal format or have a ! negative exponent. ! ! Output, integer IERROR, 0 for no error, 1 for an error. ! ! Output, integer LCHAR, the number of characters read from ! STRING to form the number. ! implicit none ! integer i integer i_gcd integer ibot integer ibotb integer ierror integer itemp integer itop integer itopb integer lchar integer lchar2 integer nchar character ( len = * ) string ! itop = 0 ibot = 1 lchar = 0 call chrctf ( string, itop, ibot, ierror, lchar ) if ( ierror /= 0 ) then return end if ! ! The number is represented as a fraction. ! If the next nonblank character is "/", then read another number. ! nchar = len_trim ( string ) do i = lchar+1, nchar-1 if ( string(i:i) == '/' ) then call chrctf ( string(i+1:), itopb, ibotb, ierror, lchar2 ) if ( ierror /= 0 ) then return end if itop = itop * ibotb ibot = ibot * itopb itemp = i_gcd ( itop, ibot ) itop = itop / itemp ibot = ibot / itemp lchar = i + lchar2 return else if ( string(i:i) /= ' ' ) then return end if end do return end subroutine chrinp ( ierror, line, prompt ) ! !******************************************************************************* ! !! CHRINP requests new input if the LINE buffer is empty. ! ! ! Discussion: ! ! CHRINP 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: ! ! 16 May 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer IERROR, 0 for no error, 1 for an error. ! ! 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. ! On input, the prompt string to be printed. ! On output, PROMPT has been blanked out, up to the first comma. ! implicit none ! integer i integer icomma integer ierror integer ios integer lchar character ( len = 80 ) line character ( len = 80 ) prompt ! ierror = 0 if ( line == ' ' ) then write ( *, '(a)' ) 'Enter ' // trim ( prompt ) read ( *, '(a)', iostat = ios ) line if ( ios /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'CHRINP - Warning!' write ( *, '(a)' ) ' End of input!' stop end if call s_blanks_delete ( line ) end if ! ! If item was read, remove item from PROMPT list. ! if ( line /= ' ' ) then icomma = index ( prompt, ',' ) if ( icomma > 0 .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 col_add ( a, a_int, iatop, iabot, ierror, base, iform, icol1, & icol2, nrow, ncol, sval, istop, isbot ) ! !******************************************************************************* ! !! COL_ADD adds a multiple of one column to another. ! ! ! Modified: ! ! 20 February 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, real A(NROW,NCOL), the current matrix. ! ! Input/output, integer IATOP(NROW,NCOL), IABOT(NROW,NCOL). ! IATOP and IABOT represent the current rational or decimal matrix. ! ! Output, integer IERROR, 0 = no error, 1 = an error. ! ! Input, integer IFORM, 0/1/2/3/4, rat/real/dec/int/mod arithmetic. ! ! Input, integer ICOL1, the column which is to be modified. ! ! Input, integer ICOL2, the column which is to be multiplied by ! a given value and added to row ICOL1. ! ! Input, integer NROW, NCOL, the number of rows and columns in the matrix. ! ! Input, real SVAL, the multiplier to use if real arithmetic is employed. ! ! Input, integer ISTOP, ISBOT, the fractional or decimal ! multiplier to use. ! implicit none ! integer ncol integer nrow ! real a(nrow,ncol) integer a_int(nrow,ncol) integer base character ( len = 10 ) chrtmp1 character ( len = 10 ) chrtmp2 character ( len = 22 ) chrtmp3 integer i integer iabot(nrow,ncol) integer iatop(nrow,ncol) integer ibot integer icol1 integer icol2 integer ierror integer iform integer isbot integer isbot2 integer istop integer istop2 integer itop real sval ! ! Return immediately if the multiplier is zero. ! if ( ( iform == 0 .and. istop == 0 ) .or. & ( iform == 1 .and. sval == 0.0E+00 ) .or. & ( iform == 2 .and. istop == 0 ) .or. & ( iform == 3 .and. istop == 0 ) .or. & ( iform == 4 .and. istop == 0 ) ) then return end if ! ! Carry out the operation. ! if ( iform == 0 ) then do i = 1, nrow call rat_mul ( isbot2, isbot, iabot(i,icol2), istop2, istop, & iatop(i,icol2), ierror ) if ( ierror /= 0 ) then return end if call rat_add ( ibot, iabot(i,icol1), isbot2, itop, iatop(i,icol1), & istop2 ) iatop(i,icol1) = itop iabot(i,icol1) = ibot end do else if ( iform == 1 ) then a(1:nrow,icol1) = a(1:nrow,icol1) + sval * a(1:nrow,icol2) else if ( iform == 2 ) then do i = 1, nrow call dec_mul ( isbot2, isbot, iabot(i,icol2), istop2, istop, & iatop(i,icol2) ) call dec_add ( ibot, iabot(i,icol1), isbot2, itop, iatop(i,icol1), & istop2 ) iatop(i,icol1) = itop iabot(i,icol1) = ibot end do else if ( iform == 3 ) then a_int(1:nrow,icol1) = a_int(1:nrow,icol1) + istop * a_int(1:nrow,icol2) else if ( iform == 4 ) then a_int(1:nrow,icol1) = & mod ( a_int(1:nrow,icol1) + istop * a_int(1:nrow,icol2), base ) end if ! ! Print out a message. ! if ( iform == 0 ) then if ( istop == isbot ) then chrtmp3 = '+' else if ( istop == - isbot ) then chrtmp3 = '-' else call rat_to_s_left ( istop, isbot, chrtmp3 ) end if else if ( iform == 1 ) then if ( sval == 1.0E+00 ) then chrtmp3 = '+' else if ( sval == - 1.0E+00 ) then chrtmp3 = '-' else call r_to_s_left ( sval, chrtmp3 ) end if else if ( iform == 2 ) then if ( istop == 1 .and. isbot == 0 ) then chrtmp3 = '+' else if ( istop == - 1 .and. isbot == 0 ) then chrtmp3 = '-' else call dec_to_s_left ( istop, isbot, chrtmp3 ) end if else if ( iform == 3 ) then if ( istop == 1 ) then chrtmp3 = '+' else if ( istop == - 1 ) then chrtmp3 = '-' else call i_to_s_left ( istop, chrtmp3 ) end if else if ( iform == 4 ) then if ( istop == 1 ) then chrtmp3 = '+' else if ( istop == - 1 ) then chrtmp3 = '-' else call i_to_s_left ( istop, chrtmp3 ) end if end if call i_to_s_left ( icol1, chrtmp1 ) call i_to_s_left ( icol2, chrtmp2 ) if ( chrtmp3 == '-' .or. chrtmp3 == '+' ) then write ( *, '(a)' ) ' ECO: Col ' // trim ( chrtmp1 ) // ' <= ' // 'Col ' // & trim ( chrtmp1 ) // ' ' // chrtmp3(1:1) // ' Col ' // trim ( chrtmp2 ) else if ( chrtmp3(1:1) == '-' .or. chrtmp3(1:1) == '+' ) then write ( *, '(a)' ) ' ECO: Col ' // trim ( chrtmp1 ) // ' <= ' // 'Col ' // & trim ( chrtmp1 ) // ' ' // chrtmp3(1:1) // ' ' // trim ( chrtmp3(2: ) ) // & ' Col ' // trim ( chrtmp2 ) else write ( *, '(a)' ) ' ECO: Col ' // trim ( chrtmp1 ) // ' <= ' // 'Col ' // & trim ( chrtmp1 ) // ' + ' // trim ( chrtmp3 ) // ' Col ' // trim ( chrtmp2 ) end if return end subroutine col_add_param ( ierror, base, iform, icol1, icol2, istop, isbot, & line, ncol, sval ) ! !******************************************************************************* ! !! COL_ADD_PARAM gets and checks the column add parameters. ! ! ! Modified: ! ! 20 February 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer IERROR, 0 = no error, 1 = an error. ! ! Input, integer IFORM, 0/1/2/3/4, rat/real/dec/int/mod arithmetic. ! ! Input, integer ICOL1, the column to which the multiple is to be added. ! ! Input, integer ICOL2, the column which is to be multiplied and ! added to another row. ! ! Output, integer ISTOP, ISBOT, the parts of the rational ! or decimal fraction of the multiplier, if that is the ! arithmetic being used. ! ! Workspace, character ( len = 80 ) LINE, used to hold the user's input. ! ! Input, integer NCOL, the number of columns in the matrix. ! ! Output, real SVAL, the multiplier, if real arithmetic is used. ! implicit none ! integer base integer icol1 integer icol2 integer ierror integer iform integer isbot integer istop character ( len = 80 ) line integer ncol character ( len = 80 ) prompt real sval ! prompt = 'multiplier S, column I to add, target column J.' ! ! Get the multiplier, SVAL or ISTOP/ISBOT. ! if ( iform == 0 ) then call rat_read ( istop, isbot, line, prompt, ierror ) else if ( iform == 1 ) then call r_read ( sval, line, prompt, ierror ) else if ( iform == 2 ) then call dec_read ( istop, isbot, line, prompt, ierror ) if ( ierror == 0 ) then call dec_round ( istop, isbot ) end if else if ( iform == 3 ) then call i_read ( istop, line, prompt, ierror ) else if ( iform == 4 ) then call i_read ( istop, line, prompt, ierror ) istop = mod ( istop, base ) end if if ( ierror /= 0 ) then return end if ! ! Get the column to add. ! call i_read ( icol2, line, prompt, ierror ) if ( ierror /= 0 ) then ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Error reading column index from line:' write ( *, '(a)' ) '"' // trim ( line ) // '"' return end if if ( icol2 < 1 .or. icol2 > ncol ) then ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'Error! Column index was not acceptable!' return end if ! ! Get the column to which we are adding. ! call i_read ( icol1, line, prompt, ierror ) if ( ierror /= 0 ) return if ( icol1 < 1 .or. icol1 > ncol ) then ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'Error! The column index was not acceptable!' return end if ! ! Make sure the columns are different. ! if ( icol1 == icol2 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'Error! The columns should not be the same!' ierror = 1 return end if ierror = 0 return end subroutine col_app ( a, a_int, iatop, iabot, base, ierror, iform, line, & nrow, ncol ) ! !******************************************************************************* ! !! COL_APP allows the user to append a column to the matrix. ! ! ! Modified: ! ! 09 October 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, real A(NROW,NCOL), the real matrix. ! ! Input/output, real A_INT(NROW,NCOL), the int/mod matrix. ! ! Input/output, integer IATOP(NROW,NCOL), IABOT(NROW,NCOL), the rat/dec matrix. ! ! Input, integer BASE, the base of modular arithmetic. ! ! Output, integer IERROR, 0 = no error, 1 = an error. ! ! Input, integer IFORM, 0/1/2/3/4, rat/real/dec/int/mod arithmetic. ! ! Workspace, character ( len = 80 ) LINE, used to hold the user's input. ! ! Input, integer NROW, NCOL, the number of rows and columns in the matrix. ! implicit none ! integer ncol integer nrow ! real a(nrow,ncol) integer a_int(nrow,ncol) integer base integer iabot(nrow,ncol) integer iatop(nrow,ncol) integer icol integer ierror integer iform integer irow character ( len = 80 ) line ! icol = ncol do irow = 1, nrow call value_read ( irow, icol, iform, base, nrow, ncol, line, a, a_int, & iatop, iabot, ierror ) if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'COL_APP - Error!' write ( *, '(a)' ) ' VALUE_READ returned an error.' return end if end do return end subroutine col_auto ( a, a_int, iatop, iabot, ierror, base, iform, nrow, ncol ) ! !******************************************************************************* ! !! COL_AUTO automatically column reduces the current matrix. ! ! ! Modified: ! ! 29 March 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, real A(NROW,NCOL), the current matrix. ! ! Input/output, integer IATOP(NROW,NCOL), IABOT(NROW,NCOL). ! IATOP and IABOT represent the rational or decimal matrix ! to which elementary row operations will be applied. ! ! Output, integer IERROR, 0 an error occurred, 1 no error occurred. ! ! Input, integer IFORM, 0/1/2/3/4, rat/real/dec/int/mod arithmetic. ! ! Input, integer NROW, NCOL, the number of rows and columns in the matrix. ! implicit none ! integer ncol integer nrow ! real a(nrow,ncol) integer a_int(nrow,ncol) real amax real atemp integer base integer i integer iabot(nrow,ncol) integer iatop(nrow,ncol) integer ierror integer iform integer imax integer irow integer isbot integer istop integer j integer jcol integer kcol integer l integer lcol real sval ! ierror = 0 kcol = 0 do j = 1, ncol jcol = j do i = 1, nrow irow = i ! ! In row IROW, seek the column between ICOL and NCOL with ! maximum nonzero entry AMAX. ! imax = 0 amax = 0.0E+00 do kcol = jcol, ncol if ( iform == 0 ) then call rat_to_r ( iatop(irow,kcol), iabot(irow,kcol), atemp ) else if ( iform == 1 ) then atemp = a(irow,kcol) else if ( iform == 2 ) then call dec_to_r ( iatop(irow,kcol), iabot(irow,kcol), atemp ) else if ( iform == 3 ) then atemp = real ( a_int(irow,kcol) ) else if ( iform == 4 ) then atemp = real ( a_int(irow,kcol) ) end if atemp = abs ( atemp ) if ( atemp > amax ) then amax = atemp imax = kcol end if end do if ( imax /= 0 ) then kcol = imax exit end if end do if ( kcol == 0 ) then return end if ! ! Interchange the JCOL-th and the pivot columns. ! if ( kcol /= jcol ) then call col_swap ( a, a_int, iatop, iabot, ierror, iform, kcol, jcol, & nrow, ncol ) end if ! ! Divide the pivot column by A(IROW,JCOL) so that A(IROW,JCOL) = 1. ! if ( iform == 0 ) then istop = iatop(irow,jcol) isbot = iabot(irow,jcol) else if ( iform == 1 ) then sval = a(irow,jcol) else if ( iform == 2 ) then istop = iatop(irow,jcol) isbot = iabot(irow,jcol) else if ( iform == 3 ) then istop = a_int(irow,jcol) else if ( iform == 4 ) then istop = a_int(irow,jcol) end if call col_div ( a, a_int, iatop, iabot, ierror, base, iform, jcol, & nrow, ncol, sval, istop, isbot ) ! ! Annihilate A(IROW,L) for L not equal to JCOL. ! do l = 1, nrow lcol = l if ( lcol /= jcol ) then if ( iform == 0 ) then if ( iatop(irow,lcol) /= 0 ) then istop = - iatop(irow,lcol) isbot = iabot(irow,lcol) call col_add ( a, a_int, iatop, iabot, ierror, base, iform, & lcol, jcol, nrow, ncol, sval, istop, isbot ) iatop(irow,lcol) = 0 iabot(irow,lcol) = 1 end if else if ( iform == 1 ) then if ( a(irow,lcol) /= 0.0E+00 ) then sval = - a(irow,lcol) call col_add ( a, a_int, iatop, iabot, ierror, base, iform, & lcol, jcol, nrow, ncol, sval, istop, isbot ) a(irow,lcol) = 0.0E+00 end if else if ( iform == 2 ) then if ( iatop(irow,lcol) /= 0 ) then istop = - iatop(irow,lcol) isbot = iabot(irow,lcol) call col_add ( a, a_int, iatop, iabot, ierror, base, iform, & lcol, jcol, nrow, ncol, sval, istop, isbot ) iatop(irow,lcol) = 0 iabot(irow,lcol) = 0 end if else if ( iform == 3 ) then if ( a_int(irow,lcol) /= 0 ) then istop = - a_int(irow,lcol) call col_add ( a, a_int, iatop, iabot, ierror, base, iform, & lcol, jcol, nrow, ncol, sval, istop, isbot ) a_int(irow,lcol) = 0 end if else if ( iform == 4 ) then if ( a_int(irow,lcol) /= 0 ) then istop = - a_int(irow,lcol) call col_add ( a, a_int, iatop, iabot, ierror, base, iform, & lcol, jcol, nrow, ncol, sval, istop, isbot ) a_int(irow,lcol) = 0 end if end if end if end do end do return end subroutine col_div ( a, a_int, iatop, iabot, ierror, base, iform, icol, & nrow, ncol, sval, istop, isbot ) ! !******************************************************************************* ! !! COL_DIV divides a column of the matrix by a scale factor. ! ! ! Modified: ! ! 20 February 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, real A(NROW,NCOL), the current matrix. ! ! Input/output, integer IATOP(NROW,NCOL), IABOT(NROW,NCOL). ! IATOP and IABOT represent the current rational or decimal ! matrix. ! ! Output, integer IERROR, 0 an error occurred, 1 no error occurred. ! ! Input, integer IFORM, 0/1/2/3/4, rat/real/dec/int/mod arithmetic. ! ! Input, integer ICOL, the column to be divided. ! ! Input, integer NROW, NCOL, the number of rows and columns in the matrix. ! ! Input, real SVAL, the real divisor. ! ! Input, integer ISTOP, ISBOT, the fractional or decimal divisor. ! implicit none ! integer ncol integer nrow ! real a(nrow,ncol) integer a_int(nrow,ncol) integer base character ( len = 10 ) chrtmp1 character ( len = 24 ) chrtmp2 integer i integer iabot(nrow,ncol) integer iatop(nrow,ncol) integer ibot integer icol integer ierror integer iform integer isbot integer istop integer itop character ( len = 3 ) op real sval ! ! Make sure that the column number is legal. ! if ( icol < 1 .or. icol > ncol ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Error! The column number is out of range!' ierror = 1 return end if ! ! Check for an illegal divisor of 0. ! if ( ( iform == 0 .and. istop == 0 ) .or. & ( iform == 1 .and. sval == 0.0E+00 ) .or. & ( iform == 2 .and. istop == 0 ) .or. & ( iform == 3 .and. istop == 0 ) .or. & ( iform == 4 .and. istop == 0 ) ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'COL_DIV - Error!' write ( *, '(a)' ) ' It is illegal to divide by 0!' ierror = 1 return end if ! ! Check for a pointless division by 1. ! if ( ( iform == 0 .and. istop == isbot ) .or. & ( iform == 1 .and. sval == 1.0E+00 ) .or. & ( iform == 2 .and. istop == 1 .and. isbot == 0 ) .or. & ( iform == 3 .and. istop == 1 ) .or. & ( iform == 4 .and. istop == 1 ) ) then return end if ! ! Carry out the division. ! if ( iform == 0 ) then do i = 1, nrow call rat_mul ( ibot, iabot(i,icol), istop, itop, iatop(i,icol), isbot, & ierror ) if ( ierror /= 0 ) then return end if iatop(i,icol) = itop iabot(i,icol) = ibot end do else if ( iform == 1 ) then a(1:nrow,icol) = a(1:nrow,icol) / sval else if ( iform == 2 ) then do i = 1, nrow call dec_div ( iatop(i,icol), iabot(i,icol), istop, isbot, & iatop(i,icol), iabot(i,icol), ierror ) end do else if ( iform == 3 ) then a_int(1:nrow,icol) = a_int(1:nrow,icol) / istop else if ( iform == 4 ) then a_int(1:nrow,icol) = a_int(1:nrow,icol) / istop end if ! ! Print out a statement about what has been done. ! if ( iform == 0 ) then if ( isbot == 1 ) then call i_to_s_left ( istop, chrtmp2 ) op = ' / ' else call rat_to_s_left ( isbot, istop, chrtmp2 ) op = ' * ' end if else if ( iform == 1 ) then call r_to_s_left ( sval, chrtmp2 ) op = ' / ' else if ( iform == 2 ) then call dec_to_s_left ( istop, isbot, chrtmp2 ) op = ' / ' else if ( iform == 3 ) then call i_to_s_left ( istop, chrtmp2 ) op = ' / ' else if ( iform == 4 ) then call i_to_s_left ( istop, chrtmp2 ) op = ' / ' end if call i_to_s_left ( icol, chrtmp1 ) write ( *, '(a)' ) ' ECO: Col ' // trim ( chrtmp1 ) // ' <= Col ' // & trim ( chrtmp1 ) // op // trim ( chrtmp2 ) return end subroutine col_div_param ( icol, ierror, base, iform, isbot, istop, line, sval ) ! !******************************************************************************* ! !! COL_DIV_PARAM gets and checks the column divide parameters. ! ! ! Modified: ! ! 12 July 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer IERROR, 0 an error occurred, 1 no error occurred. ! ! Input, integer IFORM, 0/1/2/3/4, rat/real/dec/int/mod arithmetic. ! ! Workspace, character ( len = 80 ) LINE, used to hold the user's input. ! implicit none ! integer base integer icol integer ierror integer iform integer isbot integer istop character ( len = 80 ) line character ( len = 80 ) prompt real sval ! prompt = 'column I, divisor S.' ! ! Read the column number to be divided. ! call i_read ( icol, line, prompt, ierror ) if ( ierror /= 0 ) then return end if ! ! Read the divisor, either RVAL or ISTOP/ISBOT. ! if ( iform == 0 ) then call rat_read ( istop, isbot, line, prompt, ierror ) else if ( iform == 1 ) then call r_read ( sval, line, prompt, ierror ) else if ( iform == 2 ) then call dec_read ( istop, isbot, line, prompt, ierror ) if ( ierror /= 0 ) then call dec_round ( istop, isbot ) end if else if ( iform == 3 ) then call i_read ( istop, line, prompt, ierror ) else if ( iform == 4 ) then call i_read ( istop, line, prompt, ierror ) istop = mod ( istop, base ) end if return end subroutine col_mul ( a, a_int, iatop, iabot, ierror, base, iform, icol, & nrow, ncol, sval, istop, isbot ) ! !******************************************************************************* ! !! COL_MUL multiplies a column of the matrix by a scale factor. ! ! ! Modified: ! ! 20 February 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, real A(NROW,NCOL), the current matrix. ! ! Input/output, integer IATOP(NROW,NCOL), IABOT(NROW,NCOL). ! IATOP and IABOT represent the current rational or decimal matrix. ! ! Output, integer IERROR, 0 an error occurred, 1 no error occurred. ! ! Input, integer IFORM, 0/1/2/3/4, rat/real/dec/int/mod arithmetic. ! ! Input, integer ICOL, the column that is to be multiplied. ! ! Input, integer NROW, NCOL, the number of rows and columns in the matrix. ! ! Input, real SVAL, the real row multiplier. ! ! Input, integer ISTOP, ISBOT, the decimal or fractional row multiplier. ! implicit none ! integer ncol integer nrow ! real a(nrow,ncol) integer a_int(nrow,ncol) integer base character ( len = 10 ) chrtmp1 character ( len = 22 ) chrtmp integer i integer iabot(nrow,ncol) integer iatop(nrow,ncol) integer ibot integer icol integer ierror integer iform integer isbot integer istop integer itop real sval ! ! Make sure column number is OK. ! if ( icol < 1 .or. icol > ncol ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Error! The column number is out of range!' ierror = 1 return end if ! ! For rational arithmetic, make sure bottom of scale factor ! is not 0. ! if ( iform == 0 .and. isbot == 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Error! Illegal 0 divisor in multiplier!' ierror = 1 return end if ! ! Check for multiplication by 0. ! if ( ( iform == 0 .and. istop == 0 ) .or. & ( iform == 1 .and. sval == 0.0E+00 ) .or. & ( iform == 2 .and. istop == 0 ) .or. & ( iform == 3 .and. istop == 0 ) .or. & ( iform == 4 .and. istop == 0 ) ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Warning - Multiplication by zero is not an ERO.' ierror = 1 return end if ! ! Check for multiplication by 1. ! if ( ( iform == 0 .and. istop == isbot ) .or. & ( iform == 1 .and. sval == 1.0E+00 ) .or. & ( iform == 2 .and. istop == 1 .and. isbot == 0 ) .or. & ( iform == 3 .and. istop == 1 ) .or. & ( iform == 4 .and. istop == 1 ) ) then return end if ! ! Carry out the multiplication. ! if ( iform == 0 ) then do i = 1, nrow call rat_mul ( ibot, iabot(i,icol), isbot, itop, iatop(i,icol), istop, & ierror ) if ( ierror /= 0 ) then return end if iatop(i,icol) = itop iabot(i,icol) = ibot end do else if ( iform == 1 ) then a(1:nrow,icol) = sval * a(1:nrow,icol) else if ( iform == 2 ) then do i = 1, nrow call dec_mul ( ibot, iabot(i,icol), isbot, itop, iatop(i,icol), istop ) if ( ierror /= 0 ) return iatop(i,icol) = itop iabot(i,icol) = ibot end do else if ( iform == 3 ) then a_int(1:nrow,icol) = istop * a_int(1:nrow,icol) else if ( iform == 3 ) then a_int(1:nrow,icol) = mod ( istop * a_int(1:nrow,icol), base ) end if ! ! Confirm the operation. ! if ( iform == 0 ) then if ( istop == - isbot ) then chrtmp = '-' else call rat_to_s_left ( istop, isbot, chrtmp ) end if else if ( iform == 1 ) then if ( sval == - 1.0E+00 ) then chrtmp = '-' else call r_to_s_left ( sval, chrtmp ) end if else if ( iform == 2 ) then if ( istop == - 1 .and. isbot == 0 ) then chrtmp = '-' else call dec_to_s_left ( istop, isbot, chrtmp ) end if else if ( iform == 3 ) then if ( istop == - 1 ) then chrtmp = '-' else call i_to_s_left ( istop, chrtmp ) end if else if ( iform == 4 ) then if ( istop == - 1 ) then chrtmp = '-' else call i_to_s_left ( istop, chrtmp ) end if end if call i_to_s_left ( icol, chrtmp1 ) write ( *, '(a)' ) ' ECO: Col ' // trim ( chrtmp1 ) // ' <= ' // trim ( chrtmp ) // & ' Col ' // trim ( chrtmp1 ) return end subroutine col_mul_param ( ierror, base, iform, icol, istop, isbot, line, rval ) ! !******************************************************************************* ! !! COL_MUL_PARAM gets and checks the column multiply parameters. ! ! ! Modified: ! ! 20 February 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer IERROR, 0 an error occurred, 1 no error occurred. ! ! Input, integer IFORM, 0/1/2/3/4, rat/real/dec/int/mod arithmetic. ! ! Output, integer ICOL, the column to be multiplied. ! ! Output, integer ISTOP, ISBOT, the multiplier to use for ! fractional or decimal arithmetic. ! ! Workspace, character ( len = 80 ) LINE, used to hold the user's input. ! ! Output, real RVAL, the multiplier to use for real arithmetic. ! implicit none ! integer base integer icol integer ierror integer iform integer isbot integer istop character ( len = 80 ) line character ( len = 80 ) prompt real rval ! prompt = 'column I, multiplier S.' ! ! Read the column number to be multiplied. ! call i_read ( icol, line, prompt, ierror ) if ( ierror /= 0 ) return ! ! Read the multiplier, either RVAL or ISTOP/ISBOT. ! if ( iform == 0 ) then call rat_read ( istop, isbot, line, prompt, ierror ) else if ( iform == 1 ) then call r_read ( rval, line, prompt, ierror ) else if ( iform == 2 ) then call dec_read ( istop, isbot, line, prompt, ierror ) if ( ierror == 0 ) then call dec_round ( istop, isbot ) end if else if ( iform == 3 ) then call i_read ( istop, line, prompt, ierror ) else if ( iform == 4 ) then call i_read ( istop, line, prompt, ierror ) istop = mod ( istop, base ) end if if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'COL_MUL_PARAM - Error!' write ( *, '(a)' ) ' Your input could not be understood.' return end if return end subroutine col_op_check ( command, ierror, line2 ) ! !******************************************************************************* ! !! COL_OP_CHECK checks for commands given in the form of ECO's. ! ! ! Discussion: ! ! The form of the elementary column operation commands includes: ! ! The column interchange command: ! CI1 <=> CI2 ! Note that this will fail if user types "C I1 <=> C I2" ! ! The scalar multiply command: ! CI1 <= S * CI1 ! with or without the "*". ! ! The scalar divide command: ! CI1 <= CI1 / S ! ! The add row command: ! CI1 <= CI1 + S *CI2 ! or ! CI1 <= S * CI2 + CI1 ! ! Modified: ! ! 20 February 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, character ( len = 4 ) COMMAND. ! If the routine decides that the user has input an ERO in the ! natural format, then COMMAND contains the necessary ! one letter MATMAN command to carry out the ERO. ! ! Output, integer IERROR, 0 an error occurred, 1 no error occurred. ! ! Workspace, character ( len = 80 ) LINE2, a copy of the user input in LINE. ! implicit none ! character ( len = 22 ) chrtmp character ( len = 10 ) chrtmp1 character ( len = 10 ) chrtmp2 character ( len = 20 ) command integer icol1 integer icol2 integer icol3 integer idbot2 integer idbot3 integer idtop2 integer idtop3 integer ierror integer isbot2 integer isbot3 integer istop2 integer istop3 integer lchar logical ldiv character ( len = 80 ) line2 character ( len = 80 ) string ! command = ' ' ! ! 1. Remove all blanks from the line, and capitalize it. ! call s_blank_delete ( line2 ) call s_cap ( line2 ) ! ! 2. Is the first character a "C" or "COL" or "COLUMN"? ! if ( line2(1:1) /= 'C' ) then return end if if ( line2(1:6) == 'COLUMN' ) then call s_chop ( line2, 1, 6 ) else if ( line2(1:3) == 'COL' ) then call s_chop ( line2, 1, 3 ) else call s_chop ( line2, 1, 1 ) end if ! ! 3. The next item should be a column number, ICOL1. ! call s_to_i ( line2, icol1, ierror, lchar ) if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Your ECO command could not be understood.' write ( *, '(a)' ) ' The first column number "R1" did not make sense.' return end if call s_chop ( line2, 1, lchar ) ! ! 4. Check for the column interchange string "=", "<>", "<=>" or "<->". ! if ( line2(1:2) == '<>' ) then string = '<>' else if ( line2(1:3) == '<=>' ) then string = '<=>' else if ( line2(1:3) == '<->' ) then string = '<->' else if ( line2(1:2) == '<=' ) then string = '<=' else if ( line2(1:2) == '<-' ) then string = '<-' else if ( line2(1:2) == '=>' ) then string = '=>' else if ( line2(1:2) == '->' ) then string = '->' else if ( line2(1:1) == '=' ) then string = '=' else if ( line2(1:2) == ':=' ) then string = ':=' else ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Your ECO command could not be understood.' write ( *, '(a)' ) ' The assignment symbol <=> was missing.' return end if lchar = len_trim ( string ) call s_chop ( line2, 1, lchar ) ! ! 5. The next quantity could be an explicit signed scalar, S2, ! or an implicit +-1. ! if ( line2(1:1) == 'C' ) then istop2 = 1.0E+00 isbot2 = 1.0E+00 else if ( line2(1:2) == '+C' ) then istop2 = 1.0E+00 isbot2 = 1.0E+00 call s_chop ( line2, 1, 1 ) else if ( line2(1:2) == '-C' ) then istop2 = - 1.0E+00 isbot2 = 1.0E+00 call s_chop ( line2, 1, 1 ) else call chrctg ( line2, istop2, isbot2, ierror, lchar ) call s_chop ( line2, 1, lchar ) if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Your ECO command could not be understood.' write ( *, '(a)' ) ' The multiplier S2 did not make sense.' ierror = 1 return end if end if end if ! ! 6. Is the next character an optional "*"? ! if ( line2(1:1) == '*' ) then call s_chop ( line2, 1, 1 ) end if ! ! 7. Is the next character a "C"? ! if ( line2(1:6) == 'COLUMN' ) then call s_chop ( line2, 1, 6 ) else if ( line2(1:3) == 'COL' ) then call s_chop ( line2, 1, 3 ) else if ( line2(1:1) == 'C' ) then call s_chop ( line2, 1, 1 ) else ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Your ECO command could not be understood.' write ( *, '(a)' ) ' Could not find the second column index.' return end if ! ! 8. The next item should be a column number, ICOL2. ! call s_to_i ( line2, icol2, ierror, lchar ) if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Your ECO command could not be understood.' write ( *, '(a)' ) ' The index of the second column "C2" did not make sense.' ierror = 1 return end if call s_chop ( line2, 1, lchar ) ! ! 9. If there's nothing more, this must be an interchange ! or a scaling. Form the equivalent MATMAN command. ! if ( line2 == ' ' ) then if ( icol1 == icol2 ) then command = 'COL_MUL' if ( isbot2 < 0 ) then isbot2 = - isbot2 istop2 = - istop2 end if call rat_to_s_left ( istop2, isbot2, chrtmp ) call i_to_s_left ( icol1, chrtmp1 ) line2 = chrtmp1 // ' ' // chrtmp call s_blanks_delete ( line2 ) return end if if ( istop2 == 1 .and. isbot2 == 1 ) then command = 'COL_SWAP' call i_to_s_left ( icol1, chrtmp1 ) call i_to_s_left ( icol2, chrtmp2 ) line2 = chrtmp1 // ' ' // chrtmp2 call s_blanks_delete ( line2 ) return end if ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Your ECO command could not be understood.' write ( *, '(a)' ) ' A MULTIPLY command must have C1 and C2 the same.' write ( *, '(a)' ) ' An INTERCHANGE command cannot have a multiplier.' return end if ! ! 10. Is the next quantity a '/', or perhaps a '*'? ! ldiv = .false. if ( line2(1:1) == '/' ) then ldiv = .true. call s_chop ( line2, 1, 1 ) call chrctg ( line2, idtop2, idbot2, ierror, lchar ) if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Your ECO command could not be understood.' write ( *, '(a)' ) ' The divisor of column 2 did not make sense.' return end if istop2 = istop2 * idbot2 isbot2 = isbot2 * idtop2 if ( icol1 == icol2 ) then if ( ldiv ) then command = 'COL_DIV' call i_swap ( istop2, isbot2 ) else command = 'COL_MUL' end if if ( isbot2 < 0 ) then isbot2 = - isbot2 istop2 = - istop2 end if call i_to_s_left ( icol1, chrtmp1 ) call rat_to_s_left ( istop2, isbot2, chrtmp ) line2 = chrtmp1 // ' ' // chrtmp call s_blanks_delete ( line2 ) return end if else if ( line2(1:1) == '*' ) then call s_chop ( line2, 1, 1 ) call chrctg ( line2, idtop2, idbot2, ierror, lchar ) if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Your ECO command could not be understood.' write ( *, '(a)' ) ' The multiplier of column 2 did not make sense.' return end if istop2 = istop2 * idtop2 isbot2 = isbot2 * idbot2 if ( icol1 == icol2 ) then command = 'COL_MUL' if ( isbot2 < 0 ) then isbot2 = - isbot2 istop2 = - istop2 end if call i_to_s_left ( icol1, chrtmp1 ) call rat_to_s_left ( istop2, isbot2, chrtmp ) line2 = chrtmp1 // ' ' // chrtmp call s_blanks_delete ( line2 ) return end if end if ! ! 11. Is the next quantity a scalar, S3? ! if ( line2(1:2) == '+C' ) then istop3 = 1.0E+00 isbot3 = 1.0E+00 call s_chop ( line2, 1, 1) else if ( line2(1:2) == '-C' ) then istop3 = - 1.0E+00 isbot3 = 1.0E+00 call s_chop ( line2, 1, 1 ) else call chrctg ( line2, istop3, isbot3, ierror, lchar ) if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Your ECO command could not be understood.' write ( *, '(a)' ) ' The multiplier S2 did not make sense.' ierror = 1 return end if call s_chop ( line2, 1, lchar ) end if ! ! 12. Is the next quantity an optional "*"? ! if ( line2(1:1) == '*' ) then call s_chop ( line2, 1, 1 ) end if ! ! 13. Is the next quantity a "C"? ! if ( line2(1:6) == 'COLUMN' ) then call s_chop ( line2, 1, 6 ) else if ( line2(1:3) == 'COL' ) then call s_chop ( line2, 1, 3 ) else if ( line2(1:1) == 'C' ) then call s_chop ( line2, 1, 1) else ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Your ECO command could not be understood.' write ( *, '(a)' ) ' The "C" marking the third column was misplaced.' return end if ! ! 14. The next item should be a column number, ICOL3. ! call s_to_i ( line2, icol3, ierror, lchar ) if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Your ECO command could not be understood.' write ( *, '(a)' ) ' The third column number "C3" did not make sense.' ierror = 1 return end if call s_chop ( line2, 1, lchar ) ! ! 15. Is the next quantity a '/', or perhaps a '*'? ! if ( line2(1:1) == '/' ) then call s_chop ( line2, 1, 1 ) call chrctg ( line2, idtop3, idbot3, ierror, lchar ) if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Your ECO command could not be understood.' write ( *, '(a)' ) ' The divisor of column 3 did not make sense.' return end if istop3 = istop3 * idbot3 isbot3 = isbot3 * idtop3 else if ( line2(1:1) == '*' ) then call s_chop ( line2, 1, 1 ) call chrctg ( line2, idtop3, idbot3, ierror, lchar ) if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Your ECO command could not be understood.' write ( *, '(a)' ) ' The multiplier of column 3 did not make sense.' return end if istop3 = istop3 * idtop3 isbot3 = isbot3 * idbot3 end if ! ! 16. Form the equivalent MATMAN command. ! if ( icol1 == icol2 ) then command = 'COL_ADD' if ( isbot3 < 0 ) then isbot3 = - isbot3 istop3 = - istop3 end if call i_to_s_left ( icol3, chrtmp1 ) call i_to_s_left ( icol1, chrtmp2 ) call rat_to_s_left ( istop3, isbot3, chrtmp ) line2 = chrtmp // ' ' // chrtmp1 // ' ' // chrtmp2 call s_blanks_delete ( line2 ) else if ( icol1 == icol3 ) then command = 'COL_ADD' if ( isbot2 < 0 ) then isbot2 = - isbot2 istop2 = - istop2 end if call rat_to_s_left ( istop2, isbot2, chrtmp ) call i_to_s_left ( icol2, chrtmp1 ) call i_to_s_left ( icol1, chrtmp2 ) line2 = chrtmp // ' ' // chrtmp1 // ' ' // chrtmp2 call s_blanks_delete ( line2 ) else ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Your ECO command could not be understood.' write ( *, '(a)' ) ' C2 or C3 must equal C1 in an ECO command.' end if return end subroutine col_swap ( a, a_int, iatop, iabot, ierror, iform, icol1, icol2, & nrow, ncol ) ! !******************************************************************************* ! !! COL_SWAP swaps two columns of a matrix. ! ! ! Modified: ! ! 20 February 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, real A(NROW,NCOL), the current matrix. ! ! Input/output, integer IATOP(NROW,NCOL), IABOT(NROW,NCOL). ! IATOP and IABOT represent the current rational or decimal matrix. ! ! Output, integer IERROR, 0 an error occurred, 1 no error occurred. ! ! Input, integer IFORM, 0/1/2/3/4, rat/real/dec/int/mod arithmetic. ! ! Input, integer ICOL1, ICOL2, the rows to swap. ! ! Input, integer NROW, NCOL, the number of rows and columns in the matrix. ! implicit none ! integer ncol integer nrow ! real a(nrow,ncol) integer a_int(nrow,ncol) character ( len = 10 ) chrtmp1 character ( len = 10 ) chrtmp2 integer iabot(nrow,ncol) integer iatop(nrow,ncol) integer i integer icol1 integer icol2 integer ierror integer iform ! ! Skip out if the two columns are the same. ! if ( icol1 == icol2 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'COL_SWAP - Error!' write ( *, '(a)' ) ' You have asked to swap a column with itself!' return end if ! ! Refuse to continue if a row is out of bounds. ! if ( ( icol1 < 1 .or. icol1 > ncol ) .or. & ( icol2 < 1 .or. icol2 > ncol ) ) then ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'COL_SWAP - Error!' write ( *, '(a)' ) ' One of the columns is illegal!' return end if ! ! Swap the columns. ! do i = 1, nrow if ( iform == 0 ) then call i_swap ( iatop(i,icol1), iatop(i,icol2) ) call i_swap ( iabot(i,icol1), iabot(i,icol2) ) else if ( iform == 1 ) then call r_swap ( a(i,icol1), a(i,icol2) ) else if ( iform == 2 ) then call i_swap ( iatop(i,icol1), iatop(i,icol2) ) call i_swap ( iabot(i,icol1), iabot(i,icol2) ) else if ( iform == 3 ) then call i_swap ( a_int(i,icol1), a_int(i,icol2) ) end if end do call i_to_s_left ( icol1, chrtmp1 ) call i_to_s_left ( icol2, chrtmp2 ) write ( *, '(a)' ) ' ECO: Col ' // trim ( chrtmp1 ) // ' <=> Col ' & // trim ( chrtmp2 ) return end subroutine d_to_dec ( dval, itop, ibot ) ! !******************************************************************************* ! !! D_TO_DEC converts a double precision quantity to a decimal representation. ! ! ! Discussion: ! ! Given the double precision value DVAL, DBLDEC computes integers ITOP and ! IBOT so that it is approximately true that: ! ! DVAL = ITOP * 10 ** IBOT ! ! Only DEC_DIGIT digits of DVAL are used in constructing the ! representation, where DEC_DIGIT is the number of digits in the ! decimal representation. ! ! Modified: ! ! 01 July 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, double precision DVAL, the value whose decimal representation ! is desired. ! ! Output, integer ITOP, IBOT, the approximate decimal representation of DVAL. ! ITOP is an integer, strictly between -10**DEC_DIGIT and 10**DEC_DIGIT. ! IBOT is an integer exponent of 10. ! implicit none ! integer dec_digit integer ibot integer itop double precision dtop double precision dval real ten1 real ten2 ! ! Special cases. ! if ( dval == 0.0D+00 ) then itop = 0 ibot = 0 return end if dec_digit = 0 call i_data ( 'GET', 'DEC_DIGIT', dec_digit ) ! ! Factor DVAL = DTOP * 10**IBOT ! dtop = dval ibot = 0 ! ! Now normalize so that 10**(DEC_DIGIT-1) <= ABS(DTOP) < 10**(DEC_DIGIT) ! ten1 = 10.0E+00**( dec_digit - 1 ) ten2 = 10.0E+00**dec_digit do while ( abs ( dtop ) < ten1 ) dtop = dtop * 10.0D+00 ibot = ibot - 1 end do do while ( abs ( dtop ) >= ten2 ) dtop = dtop / 10.0D+00 ibot = ibot + 1 end do ! ! ITOP is the integer part of DTOP, rounded. ! itop = nint ( dtop ) ! ! Now divide out any factors of ten from ITOP. ! if ( itop /= 0 ) then do while ( 10 * ( itop / 10 ) == itop ) itop = itop / 10 ibot = ibot + 1 end do end if return end subroutine dec_add ( ibot, ibot1, ibot2, itop, itop1, itop2 ) ! !******************************************************************************* ! !! DEC_ADD adds two decimal quantities. ! ! ! Discussion: ! ! The routine computes ! ! ITOP * 10**IBOT = ITOP1 * 10**IBOT1 + ITOP2 * 10**IBOT2 ! ! while trying to avoid integer overflow. ! ! Modified: ! ! 12 July 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer IBOT, the exponent of the result. ! ! Input, integer IBOT1, IBOT2, the exponents of the numbers to be added. ! ! Output, integer ITOP, the coefficient of the result. ! ! Input, integer ITOP1, ITOP2, the coefficients of the numbers to be added. ! implicit none ! integer ibot integer ibot1 integer ibot2 integer itop integer itop1 integer itop2 integer jtop1 integer jtop2 ! if ( itop1 == 0 ) then itop = itop2 ibot = ibot2 return else if ( itop2 == 0 ) then itop = itop1 ibot = ibot1 return else if ( ibot1 == ibot2 ) then itop = itop1 + itop2 ibot = ibot1 call dec_round ( itop, ibot ) return end if ! ! Line up the exponents. ! jtop1 = itop1 jtop2 = itop2 if ( ibot1 < ibot2 ) then jtop2 = jtop2 * 10**(ibot2-ibot1) else jtop1 = jtop1 * 10**(ibot1-ibot2) end if ! ! Add the coefficients. ! itop = jtop1 + jtop2 ibot = min ( ibot1, ibot2 ) ! ! Clean up the result. ! call dec_round ( itop, ibot ) return end subroutine dec_digit_set ( ierror, line ) ! !******************************************************************************* ! !! DEC_DIGIT_SET allows the user to specify the number of decimal digits. ! ! ! Discussion: ! ! DEC_DIGIT is ! ! * the number of digits used when converting a real number ! to a fraction using the "FI" or "FD" command; ! ! * the maximum number of digits in a decimal. ! ! Modified: ! ! 12 July 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer IERROR, 0 an error occurred, 1 no error occurred. ! ! Workspace, character ( len = 80 ) LINE, used to hold the user's input. ! implicit none ! character ( len = 10 ) chrtmp1 integer dec_digit_max integer dec_digit integer ierror integer itemp character ( len = 80 ) line character ( len = 80 ) prompt ! dec_digit_max = 0 call i_data ( 'GET', 'DEC_DIGIT_MAX', dec_digit_max ) write ( *, '(a)' ) 'How many decimal places should be used in ' write ( *, '(a)' ) 'converting real results to a decimal?' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' 1 means 123.45 becomes 1 * 10**2' write ( *, '(a)' ) ' 2 means 123.45 becomes 12 * 10**1' write ( *, '(a)' ) ' 3 means 123.45 becomes 123' write ( *, '(a)' ) 'and so on.' call i_to_s_left ( dec_digit_max, chrtmp1 ) prompt = 'number of decimals (1 to ' // chrtmp1 // ').' call s_blanks_delete ( prompt ) call i_read ( itemp, line, prompt, ierror ) if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'Your choice was not acceptable!' return end if ! ! Absolutely do not let DEC_DIGIT be less than 1. ! if ( itemp < 1 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' The number of decimals must be positive!' ierror = 1 return end if ! ! Allow user to exceed the maximum, with a warning. ! if ( itemp > dec_digit_max ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Warning!' write ( *, '(a)' ) ' Your choice is larger than the recommended maximum!' write ( *, '(a)' ) ' which is ' // trim ( chrtmp1 ) write ( *, '(a)' ) ' It is possible that calculations will break down' write ( *, '(a)' ) ' at any time! Be careful!' end if dec_digit = itemp call i_data ( 'SET', 'DEC_DIGIT', dec_digit ) call i_to_s_left ( dec_digit, chrtmp1 ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'The number of decimal digits will now be ' // chrtmp1 return end subroutine dec_div ( itop1, ibot1, itop2, ibot2, itop, ibot, ierror ) ! !******************************************************************************* ! !! DEC_DIV divides two decimal values. ! ! ! Discussion: ! ! A decimal quantity is stored as ! ! (ITOP,IBOT) ! ! representing the value ! ! ITOP * 10 ** IBOT. ! ! The routine computes ! ! ITOP * 10**IBOT = (ITOP1 * 10**IBOT1) / (ITOP2 * 10**IBOT2) ! ! = (ITOP1/ITOP2) * 10**(IBOT1-IBOT2) ! ! while avoiding integer overflow. ! ! Modified: ! ! 23 January 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ITOP1, IBOT1, the numerator. ! ! Input, integer ITOP2, IBOT2, the denominator. ! ! Output, integer ITOP, IBOT, the result. ! ! Output, integer IERROR, 0 an error occurred, 1 no error occurred. ! implicit none ! double precision dval integer ibot integer ibot1 integer ibot2 integer ibot3 integer ierror integer itop integer itop1 integer itop2 integer itop3 ! ! First special case, top fraction is 0. ! if ( itop1 == 0 ) then itop = 0 ibot = 0 return end if ! ! First error, bottom of fraction is 0. ! if ( itop2 == 0 ) then ierror = 1 itop = 0 ibot = 0 return end if ! ! Second special case, result is 1. ! if ( itop1 == itop2 .and. ibot1 == ibot2 ) then itop = 1 ibot = 0 return end if ! ! Third special case, result is power of 10. ! if ( itop1 == itop2 ) then itop = 1 ibot = ibot1 - ibot2 return end if ! ! Fourth special case: ITOP1/ITOP2 is exact. ! if ( ( itop1 / itop2 ) * itop2 == itop1 ) then itop = itop1 / itop2 ibot = ibot1 - ibot2 return end if ! ! General case. ! dval = dble ( itop1 ) / dble ( itop2 ) call d_to_dec ( dval, itop3, ibot3 ) itop = itop3 ibot = ibot3 + ibot1 - ibot2 return end subroutine dec_mul ( ibot, ibot1, ibot2, itop, itop1, itop2 ) ! !******************************************************************************* ! !! DEC_MUL multiplies two decimals. ! ! ! Discussion: ! ! The routine computes ! ! ITOP * 10**IBOT = (ITOP1 * 10**IBOT1) * (ITOP2 * 10**IBOT2) ! = (ITOP1*ITOP2) * 10**(IBOT1+IBOT2) ! ! while avoiding integer overflow. ! ! Modified: ! ! 25 January 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer IBOT, the exponent of the result. ! ! Input, integer IBOT1, the exponent of the first factor. ! ! Input, integer IBOT2, the exponent of the second factor. ! ! Output, integer ITOP, the coefficient of the result. ! ! Input, integer ITOP1, the coefficient of the first factor. ! ! Input, integer ITOP2, the coefficient of the second factor. ! implicit none ! integer i_big double precision dval integer ibot integer ibot1 integer ibot2 integer ibot3 integer itop integer itop1 integer itop2 integer itop3 real rmax real temp ! i_big = 0 call i_data ( 'GET', 'I_BIG', i_big ) rmax = real ( i_big ) ! ! The result is zero if either ITOP1 or ITOP2 is zero. ! if ( itop1 == 0 .or. itop2 == 0 ) then itop = 0 ibot = 0 return end if ! ! The result is simple if either ITOP1 or ITOP2 is one. ! if ( itop1 == 1 .or. itop2 == 1 ) then itop = itop1 * itop2 ibot = ibot1 + ibot2 return end if temp = log ( real ( abs ( itop1 ) ) ) + log ( real ( abs ( itop2 ) ) ) if ( temp < log ( rmax ) ) then itop = itop1 * itop2 ibot = ibot1 + ibot2 else dval = dble ( itop1 ) * dble ( itop2 ) call d_to_dec ( dval, itop3, ibot3 ) itop = itop3 ibot = ibot3 + ( ibot1 + ibot2 ) end if ! ! Clean up the result. ! call dec_round ( itop, ibot ) return end subroutine dec_print ( iatop, iabot, nrow, ncol, title ) ! !******************************************************************************* ! !! DEC_PRINT prints out decimal vectors and matrices. ! ! ! Modified: ! ! 12 July 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer IATOP(NROW,NCOL), IABOT(NROW,NCOL). ! IATOP and IABOT represent the decimal matrix. ! ! Input, integer NROW, NCOL, the number of rows and columns in the matrix. ! ! Input, character ( len = * ) TITLE, a label for the object being printed. ! implicit none ! integer ncol integer nrow ! character ( len = 22 ) chrtmp character ( len = 10 ) chrtmp1 character ( len = 10 ) chrtmp2 character ( len = 10 ) chrtmp3 character ( len = 40 ) format2 integer i integer iabot(nrow,ncol) integer iatop(nrow,ncol) integer ichi integer iclo integer imax integer imin integer izhi integer izlo integer j integer jmax integer jmin integer khi integer klo integer kmax integer lenc integer, parameter :: ncolum = 80 integer npline character ( len = 100 ) output character ( len = * ) title ! ! Figure out how wide we must make each column. ! imax = 0 jmax = 0 do i = 1, nrow do j = 1, ncol call dec_to_s_left ( iatop(i,j), iabot(i,j), chrtmp ) lenc = len_trim ( chrtmp ) jmax = max ( jmax, lenc ) end do end do kmax = 2 + imax + 1 + jmax npline = ncolum / kmax ! ! Set up the format for the heading. ! call i_to_s_left ( npline, chrtmp2 ) call i_to_s_left ( kmax, chrtmp3 ) format2 = '(' // chrtmp2 // 'i' // chrtmp3 // ')' call s_blank_delete ( format2 ) do jmin = 1, ncol, npline jmax = min ( jmin+npline-1, ncol ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) trim ( title ) write ( *, '(a)' ) ' ' if ( jmin > 1 .or. jmax < ncol ) then write ( output, * ) 'Columns ', jmin, ' to ', jmax call s_blanks_delete ( output ) write ( *, '(a)' ) trim ( output ) write ( *, '(a)' ) ' ' end if do i = 1, nrow output = ' ' do j = jmin, jmax klo = 4 + (j-jmin) * kmax + 1 khi = 4 + (j-jmin) * kmax + kmax call dec_to_s_left ( iatop(i,j), iabot(i,j), chrtmp ) output(klo:khi) = adjustr ( chrtmp(1:kmax) ) end do write ( *, '(a)' ) trim ( output ) end do end do return end subroutine dec_read ( itop, ibot, line, prompt, ierror ) ! !******************************************************************************* ! !! DEC_READ reads a decimal, rational or integer, and returns a decimal fraction. ! ! ! Modified: ! ! 23 January 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer ITOP, IBOT, represents the decimal fraction. ! ! Workspace, character ( len = 80 ) LINE, used to hold the user's input. ! ! Input/output, character ( len = 80 ) PROMPT, the prompt string. ! ! Output, integer IERROR, 0 an error occurred, 1 no error occurred. ! implicit none ! integer ibot integer ibot1 integer ibot2 integer ierror integer itop integer itop1 integer itop2 integer lchar character ( len = 80 ) line character ( len = 80 ) prompt ! ierror = 0 itop = 0 ibot = 0 do call chrinp ( ierror, line, prompt ) if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'DEC_READ - Fatal error!' write ( *, '(a)' ) ' CHRINP returned error flag!' return end if if ( line /= ' ' ) then exit end if end do call s_to_dec ( line, itop1, ibot1, lchar ) if ( lchar >= len ( line ) ) then itop = itop1 ibot = ibot1 else if ( line(lchar+1:lchar+1) /= '/' ) then itop = itop1 ibot = ibot1 else call s_chop ( line, 1, lchar+1 ) call s_to_dec ( line, itop2, ibot2, lchar ) call dec_div ( itop1, ibot1, itop2, ibot2, itop, ibot, ierror ) end if call s_chop ( line, 1, lchar ) return end subroutine dec_round ( itop, ibot ) ! !******************************************************************************* ! !! DEC_ROUND rounds a decimal fraction to a given number of digits. ! ! ! Discussion: ! ! The routine takes an arbitrary decimal fraction represented by ! ! ITOP * 10**IBOT ! ! and makes sure that ITOP has no more than the allowed number of ! decimal digits. ! ! Modified: ! ! 12 July 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, integer ITOP, IBOT, the coefficient and exponent ! of a decimal fraction. On return, ITOP has no more than ! the allowed number of decimal digits. ! implicit none ! integer dec_digit integer ibot integer itop ! if ( itop == 0 ) then ibot = 0 return end if dec_digit = 0 call i_data ( 'GET', 'DEC_DIGIT', dec_digit ) do while ( abs ( itop ) >= 10**dec_digit ) itop = itop / 10 ibot = ibot + 1 end do do while ( ( itop / 10 ) * 10 == itop ) itop = itop / 10 ibot = ibot + 1 end do return end subroutine dec_to_r ( itop, ibot, a ) ! !***************************************************************************** ! !! DEC_TO_R converts a decimal ITOP * 10**IBOT to a real value. ! ! ! Modified: ! ! 12 July 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ITOP, IBOT, the coefficient and exponent ! of the decimal value. ! ! Output, real A, the equivalent real value. ! implicit none ! real a integer ibot integer itop ! a = itop * 10.0E+00**ibot return end subroutine dec_to_rat ( iatop, iabot ) ! !******************************************************************************* ! !! DEC_TO_RAT converts a decimal to a rational representation. ! ! ! Discussion: ! ! On input, a value is represented as IATOP * 10**IABOT. ! ! On output, approximately the same value is represented as IATOP / IABOT. ! ! Modified: ! ! 12 July 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, integer IATOP, IABOT. ! On input, these quantities represent the value IATOP * 10 ** IABOT. ! On output, these quantities represent the value IATOP / IABOT. ! implicit none ! integer i_gcd integer iabot integer iatop integer itmp ! if ( iabot >= 0 ) then iatop = iatop * 10**iabot iabot = 1 else iabot = 10**(-iabot) itmp = i_gcd ( iatop, iabot ) iatop = iatop / itmp iabot = iabot / itmp end if return end subroutine dec_to_s_left ( ival, jval, s ) ! !******************************************************************************* ! !! DEC_TO_S_LEFT returns a left-justified representation of IVAL * 10**JVAL. ! ! ! Examples: ! ! IVAL JVAL S ! ---- ---- ------ ! 0 0 0 ! 21 3 21000 ! -3 0 -3 ! 147 -2 14.7 ! 16 -5 0.00016 ! 34 30 Inf ! 123 -21 0.0000000000000000012 ! 34 -30 0.0E+00 ! ! Modified: ! ! 13 September 1998 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer IVAL, JVAL, integers which represent the decimal. ! ! Output, character ( len = * ) S, the representation of the value. ! The string is 'Inf' or '0.0' if the value was too large ! or small to represent with a fixed point format. ! implicit none ! character ( len = 22 ) chrrep integer i integer iget1 integer iget2 integer iput1 integer iput2 integer ival integer jval integer maxdigit integer ndigit integer nleft character ( len = * ) s ! s = ' ' if ( ival == 0 ) then s = '0' return end if maxdigit = len ( s ) ! ! Store a representation of IVAL in CHRREP. ! write ( chrrep, '(i22)' ) ival call s_blank_delete ( chrrep ) ndigit = len_trim ( chrrep ) ! ! Overflow if JVAL is positive, and NDIGIT + JVAL > MAXDIGIT. ! if ( jval > 0 ) then if ( ndigit + jval > maxdigit ) then s = 'Inf' return end if end if ! ! Underflow if JVAL is negative, and 3 + NDIGIT - JVAL > MAXDIGIT. ! if ( jval < 0 ) then if ( ival > 0 ) then if ( 3 - ndigit - jval > maxdigit ) then s = '0.0' return end if else if ( 5 - ndigit - jval > maxdigit ) then s = '0.0' return end if end if end if ! ! If JVAL is nonnegative, insert trailing zeros. ! if ( jval >= 0 ) then s(1:ndigit) = chrrep(1:ndigit) do i = ndigit+1, ndigit+jval s(i:i) = '0' end do else if ( jval < 0 ) then iput2 = 0 iget2 = 0 ! ! Sign. ! if ( ival < 0 ) then iput1 = 1 iput2 = 1 iget2 = 1 s(iput1:iput2) = '-' ndigit = ndigit - 1 end if ! ! Digits of the integral part. ! if ( ndigit + jval > 0 ) then iput1 = iput2 + 1 iput2 = iput1 + ndigit + jval -1 iget1 = iget2 + 1 iget2 = iget1 + ndigit+jval - 1 s(iput1:iput2) = chrrep(iget1:iget2) else iput1 = iput2 + 1 iput2 = iput1 s(iput1:iput2) = '0' end if ! ! Decimal point. ! iput1 = iput2 + 1 iput2 = iput1 s(iput1:iput2) = '.' ! ! Leading zeroes. ! do i = 1, - jval - ndigit iput1 = iput2 + 1 iput2 = iput1 s(iput1:iput2) = '0' end do nleft = min ( -jval, ndigit ) nleft = min ( nleft, maxdigit - iput2 ) iput1 = iput2 + 1 iput2 = iput1 + nleft - 1 iget1 = iget2 + 1 iget2 = iget1 + nleft - 1 s(iput1:iput2) = chrrep(iget1:iget2) end if return end subroutine digit_to_ch ( digit, c ) ! !******************************************************************************* ! !! 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 form ( a, a_int, iatop, iabot, base, iform, jform, nrow, ncol ) ! !******************************************************************************* ! !! FORM converts from one arithmetic form to another. ! ! ! Discussion: ! ! On input, IFORM contains a code for the current arithmetic ! form, and JFORM contains the code for the new form. ! ! Modified: ! ! 12 July 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, real A(NROW,NCOL), the current real matrix. ! ! Input/output, integer IATOP(NROW,NCOL), IABOT(NROW,NCOL), ! the current fractional or decimal matrix. ! ! Input/output, integer IFORM, 0/1/2/3/4, rat/real/dec/int/mod arithmetic. ! ! Input, integer JFORM, the arithmetic to be converted to. ! ! Input, integer NROW, NCOL, the number of rows and columns in the matrix. ! implicit none ! integer ncol integer nrow ! real a(nrow,ncol) integer a_int(nrow,ncol) integer base integer i integer iabot(nrow,ncol) integer iatop(nrow,ncol) integer ierror integer iform integer j integer jform real r ! if ( iform == jform ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' You are already using the arithmetic type that' write ( *, '(a)' ) ' you have requested.' return end if ! ! Convert the matrix data. ! do i = 1, nrow do j = 1, ncol if ( jform == 0 ) then if ( iform == 1 ) then call r_to_rat ( a(i,j), iatop(i,j), iabot(i,j) ) else if ( iform == 2 ) then call dec_to_rat ( iatop(i,j), iabot(i,j) ) else if ( iform == 3 ) then iatop(i,j) = a_int(i,j) iabot(i,j) = 1 else if ( iform == 4 ) then iatop(i,j) = a_int(i,j) iabot(i,j) = 1 end if else if ( jform == 1 ) then if ( iform == 0 ) then call rat_to_r ( iatop(i,j), iabot(i,j), a(i,j) ) else if ( iform == 2 ) then call dec_to_r ( iatop(i,j), iabot(i,j), a(i,j) ) else if ( iform == 3 ) then a(i,j) = real ( a_int(i,j) ) else if ( iform == 4 ) then a(i,j) = real ( a_int(i,j) ) end if else if ( jform == 2 ) then if ( iform == 0 ) then call rat_to_dec ( iatop(i,j), iabot(i,j), ierror ) else if ( iform == 1 ) then call r_to_dec ( a(i,j), iatop(i,j), iabot(i,j) ) else if ( iform == 3 ) then iatop(i,j) = a_int(i,j) iabot(i,j) = 0 else if ( iform == 4 ) then iatop(i,j) = a_int(i,j) iabot(i,j) = 0 end if else if ( jform == 3 ) then if ( iform == 0 ) then call rat_to_r ( iatop(i,j), iabot(i,j), r ) a_int(i,j) = nint ( r ) else if ( iform == 1 ) then a_int(i,j) = nint ( a(i,j) ) else if ( iform == 2 ) then call dec_to_r ( iatop(i,j), iabot(i,j), r ) a_int(i,j) = nint ( r ) else if ( iform == 4 ) then a_int(i,j) = a_int(i,j) end if else if ( jform == 4 ) then if ( iform == 0 ) then call rat_to_r ( iatop(i,j), iabot(i,j), r ) a_int(i,j) = mod ( nint ( r ), base ) else if ( iform == 1 ) then a_int(i,j) = mod ( nint ( a(i,j) ), base ) else if ( iform == 2 ) then call dec_to_r ( iatop(i,j), iabot(i,j), r ) a_int(i,j) = mod ( nint ( r ), base ) else if ( iform == 3 ) then a_int(i,j) = mod ( a_int(i,j), base ) end if end if end do end do ! ! Update the arithmetic form. ! iform = jform return end subroutine hello ! !******************************************************************************* ! !! HELLO greets the user on startup. ! ! ! Modified: ! ! 28 September 2000 ! ! Author: ! ! John Burkardt ! implicit none ! write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'MATMAN, version 1.63' write ( *, '(a)' ) ' Last modified on 13 January 2002.' write ( *, '(a)' ) ' An interactive program to perform elementary' write ( *, '(a)' ) ' row and column operations on a matrix.' write ( *, '(a)' ) ' Developed by Charles Cullen and John Burkardt.' write ( *, '(a)' ) ' Send comments to burkardt@iastate.edu.' return end subroutine help ! !******************************************************************************* ! !! HELP prints out a brief list of the available commands. ! ! ! Modified: ! ! 25 September 2000 ! ! Author: ! ! John Burkardt ! implicit none ! write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'MATMAN commands:' write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'A(I,J)=S Set matrix entry to S.' write ( *, '(a)' ) 'CHECK Check matrix for reduced row echelon form.' write ( *, '(a)' ) 'COL_APP Append a column to the matrix.' write ( *, '(a)' ) 'COL_AUTO Automatic column reduction.' write ( *, '(a)' ) 'DEC Use decimal arithmetic.' write ( *, '(a)' ) 'DEC_DIGIT Set the number of decimal digits.' write ( *, '(a)' ) 'E Enter matrix with I rows and J columns.' write ( *, '(a)' ) 'H for help.' write ( *, '(a)' ) 'I_BIG Set size of largest integer for fractions.' write ( *, '(a)' ) 'ID Append the identity matrix.' write ( *, '(a)' ) 'INIT Initialize data.' write ( *, '(a)' ) 'INT Use integer arithmetic.' write ( *, '(a)' ) 'Q Quit.' write ( *, '(a)' ) 'RAN Set up a problem with random data.' write ( *, '(a)' ) 'RAT Use rational arithmetic.' write ( *, '(a)' ) 'REAL Use real arithmetic.' write ( *, '(a)' ) 'ROW_AUTO Automatic row reduction.' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' ERO''s and ECO''s:' write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'R1 <=> R2 C1 <=> C2' write ( *, '(a)' ) 'R1 <= S R1 C1 <= S C1' write ( *, '(a)' ) 'R1 <= R1 + S R2 C1 <= C1 + S C2' return end subroutine i_data ( op, var, ival ) ! !******************************************************************************* ! !! I_DATA stores and retrieves common data items. ! ! ! Discussion: ! ! This routine works like a sort of COMMON block. It stores or returns ! the values of certain variables. Thus, it allows routines ! to "communicate" without having to have data passed up and ! down the calling tree in argument lists. ! ! The variables stored by this version of the routine are: ! ! 'DEC_DIGIT', the number of digits stored for decimals; ! 'DEC_DIGIT_MAX', the maximum number of digits stored for decimals; ! 'I_BIG', the biggest integer to use in calculations. ! ! Modified: ! ! 20 January 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) OP, describes the operation to be done. ! 'SET' means set a value. ! 'INC' means increment a value (and return its new value) ! 'GET' means get a value. ! ! Input, character ( len = * ) VAR, the name of the variable. ! ! Input/output, integer IVAL. ! If OP is 'SET', then the variable named in VAR is set to the ! value IVAL. ! If OP is 'GET', then the value of IVAL is set to the value of ! the variable named in VAR. ! If OP is 'INC', then the value of IVAL is incremented by 1, ! and its new value is returned in VAR. ! implicit none ! integer, save :: dec_digit = 7 integer, save :: dec_digit_max = 7 integer, save :: i_big = 2147483647 integer ival character ( len = * ) op logical s_eqi character ( len = * ) var ! if ( s_eqi ( op, 'SET' ) ) then if ( s_eqi ( var, 'DEC_DIGIT' ) ) then dec_digit = ival else if ( s_eqi ( var, 'DEC_DIGIT_MAX' ) ) then dec_digit_max = ival else if ( s_eqi ( var, 'I_BIG' ) ) then i_big = ival end if else if ( s_eqi ( op, 'GET' ) ) then if ( s_eqi ( var, 'DEC_DIGIT' ) ) then ival = dec_digit else if ( s_eqi ( var, 'DEC_DIGIT_MAX' ) ) then ival = dec_digit_max else if ( s_eqi ( var, 'I_BIG' ) ) then ival = i_big end if else if ( s_eqi ( op, 'INC' ) ) then if ( s_eqi ( var, 'DEC_DIGIT' ) ) then dec_digit = dec_digit + 1 ival = dec_digit else if ( s_eqi ( var, 'DEC_DIGIT_MAX' ) ) then dec_digit_max = dec_digit_max + 1 ival = dec_digit_max else if ( s_eqi ( var, 'I_BIG' ) ) then i_big = i_big + 1 ival = i_big end if end if return end function i_gcd ( i, j ) ! !******************************************************************************* ! !! I_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 I_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, I_GCD is returned as max ( 1, abs ( I ), abs ( J ) ). ! ! If I and J have no common factor, I_GCD is returned as 1. ! ! Otherwise, using the Euclidean algorithm, I_GCD is the ! largest common factor of I and J. ! implicit none ! integer i integer i_gcd integer ip integer iq integer ir integer j ! i_gcd = 1 ! ! Return immediately if either I or J is zero. ! if ( i == 0 ) then i_gcd = max ( 1, abs ( j ) ) return else if ( j == 0 ) then i_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 i_gcd = iq return end subroutine i_print ( a_int, nrow, ncol, title ) ! !******************************************************************************* ! !! I_PRINT prints out integer vectors and matrices. ! ! ! Modified: ! ! 26 September 1998 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer A_INT(NROW,NCOL), the current matrix. ! ! Input, integer NROW, NCOL, the number of rows and columns in the matrix. ! ! Input, character ( len = * ) TITLE, a label for the object being printed. ! implicit none ! integer ncol integer nrow ! integer a_int(nrow,ncol) logical allint character ( len = 10 ) chrtmp1 character ( len = 10 ) chrtmp2 character ( len = 10 ) chrtmp3 character ( len = 40 ) format1 character ( len = 40 ) format2 integer i integer ichi integer iclo integer ihi integer ilo integer imax integer imin integer izhi integer izlo integer j integer jhi integer jlo integer jmax integer jmin integer kmax integer kmin integer, parameter :: ncolum = 80 integer npline character ( len = 100 ) output character ( len = * ) title ! ! Figure out how many numbers we can fit in NCOLUM columns. ! kmin = 11 kmax = kmin do i = 1, nrow do j = 1, ncol do while ( abs ( a_int(i,j) ) >= 10.0E+00**(kmax-kmin) ) kmax = kmax + 1 end do end do end do npline = ncolum / kmax ! ! If all integers, cut down KMAX, the width of each number, ! and update NPLINE, the number of numbers we can print on one line. ! kmax = kmax - 7 npline = ncolum / kmax call i_to_s_left ( npline, chrtmp2 ) call i_to_s_left ( kmax, chrtmp3 ) format1 = '(' // chrtmp2 // 'i' // chrtmp3 // ')' call s_blank_delete ( format1 ) ! ! The second format is for ... ! format2 = '(' // chrtmp2 // 'i' // chrtmp3 // ')' call s_blank_delete ( format2 ) ! ! Now print the data. ! do jmin = 1, ncol, npline jmax = min ( jmin + npline - 1, ncol ) write ( *, '(a)' ) ' ' if ( jmin == 1 ) then write ( *, '(a)' ) trim ( title ) write ( *, '(a)' ) ' ' end if if ( jmin > 1 .or. jmax < ncol ) then write ( output, * ) 'Columns ', jmin, ' to ', jmax call s_blanks_delete ( output ) write ( *, '(a)' ) trim ( output ) write ( *, '(a)' ) ' ' end if do i = 1, nrow write ( *, format1 ) a_int(i,jmin:jmax) end do end do return end subroutine i_random ( ilo, ihi, i ) ! !******************************************************************************* ! !! I_RANDOM returns a random integer in a given range. ! ! ! Modified: ! ! 23 September 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ILO, IHI, the minimum and maximum acceptable values. ! ! Output, integer I, the randomly chosen integer. ! implicit none ! logical, save :: seed = .false. integer i integer ihi integer ilo real r real rhi real rlo ! if ( .not. seed ) then call random_seed seed = .true. end if ! ! Pick a random number in (0,1). ! call random_number ( harvest = r ) ! ! Set a real interval [RLO,RHI] which contains the integers [ILO,IHI], ! each with a "neighborhood" of width 1. ! rlo = real ( ilo ) - 0.5E+00 rhi = real ( ihi ) + 0.5E+00 ! ! Set I to the integer that is nearest the scaled value of R. ! i = nint ( ( 1.0E+00 - r ) * rlo + r * rhi ) ! ! In case of oddball events at the boundary, enforce the limits. ! i = max ( i, ilo ) i = min ( i, ihi ) return end subroutine i_read ( intval, line, prompt, ierror ) ! !******************************************************************************* ! !! I_READ reads an integer from the input buffer. ! ! ! Discussion: ! ! The routine accepts LINE which contains input and a PROMPT line. ! If LINE is empy, the PROMPT will be printed and LINE read from the ! 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. ! ! The program 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 I_READ, S_READ, RAT_READ or ! R_READ. ! ! Modified: ! ! 12 July 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer INTVAL, the integer that was read. ! ! Workspace, character ( len = 80 ) LINE, used to hold the user's input. ! ! Input, character ( len = 80 ) PROMPT, the prompt string. ! ! Output, integer IERROR, 0 an error occurred, 1 no error occurred. ! implicit none ! integer ierror integer intval integer lchar character ( len = 80 ) line character ( len = 80 ) prompt ! ierror = 0 intval = 0 ! ! Retrieve a likely character string from input. ! do call chrinp ( ierror, line, prompt ) if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'I_READ - Fatal error!' write ( *, '(a)' ) ' CHRINP returned error flag!' return end if if ( line /= ' ' ) then exit end if end do ! ! Convert the character string to an integer. ! call s_to_i ( line, intval, ierror, lchar ) if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'I_READ - Fatal error!' write ( *, '(a)' ) ' S_TO_I returned error flag!' return end if ! ! Remove the character string from the input line. ! if ( lchar < len ( line ) ) then if ( line(lchar+1:lchar+1) == ',' ) then lchar = lchar + 1 end if end if call s_chop ( line, 1, lchar ) return end subroutine i_swap ( i, j ) ! !******************************************************************************* ! !! I_SWAP switches two integer values. ! ! ! Modified: ! ! 30 November 1998 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, integer I, J. On output, the values of I and ! J have been interchanged. ! implicit none ! integer i integer j integer k ! k = i i = j j = k return end subroutine i_to_s_left ( intval, s ) ! !******************************************************************************* ! !! I_TO_S_LEFT converts an integer to a left-justified string. ! ! ! Examples: ! ! Assume that S is 6 characters long: ! ! INTVAL S ! ! 1 1 ! -1 -1 ! 0 0 ! 1952 1952 ! 123456 123456 ! 1234567 ****** <-- Not enough room! ! ! Modified: ! ! 28 July 2000 ! ! 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 left-justified. 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 ! ! The absolute value of the integer goes into S(ILO:IHI). ! ipos = ihi ! ! Find the last digit of IVAL, strip it off, and stick it into the string. ! do 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 if ( ival == 0 ) then exit end if end do ! ! Shift the string to the left. ! s(ilo:ilo+ihi-ipos-1) = s(ipos+1:ihi) s(ilo+ihi-ipos:ihi) = ' ' return end subroutine init ( a, a_int, base, iabot, iatop, ierror, iform, & line, maxrow, maxcol, nrow, ncol ) ! !******************************************************************************* ! !! INIT initializes the data. ! ! ! Modified: ! ! 20 January 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, real A(MAXROW,MAXCOL), the current matrix. ! ! Output, integer IABOT(MAXROW,MAXCOL), IATOP(MAXROW,MAXCOL). ! IATOP and IABOT represent the current rational or decimal ! matrix. ! ! Output, integer IERROR. ! The error flag, which is initialized to zero by this routine. ! ! Output, integer IFORM, 0/1/2/3/4, rat/real/dec/int/mod arithmetic. ! ! Output, character ( len = 80 ) LINE, ! a buffer used to hold the user's input. ! ! Input, integer MAXCOL, the maximum number of matrix columns. ! ! Input, integer MAXROW, the maximum number of matrix rows. ! ! Output, integer NROW, NCOL, the number of rows and columns in the matrix. ! implicit none ! integer maxcol integer maxrow ! real a(maxrow,maxcol) integer a_int(maxrow,maxcol) integer base integer dec_digit integer dec_digit_max integer i integer i_big integer iabot(maxrow,maxcol) integer iatop(maxrow,maxcol) integer ierror integer iform character ( len = * ) line integer ncol integer nrow ! a(1,1) = 1.0E+00 a_int(1,1) = 1 iatop(1,1) = 1 iabot(1,1) = 1 base = 0 ierror = 0 iform = 0 line = ' ' dec_digit_max = 7 call i_data ( 'SET', 'DEC_DIGIT_MAX', dec_digit_max ) dec_digit = 4 call i_data ( 'SET', 'DEC_DIGIT', dec_digit ) i_big = huge ( i_big ) call i_data ( 'SET', 'I_BIG', i_big ) ncol = 0 nrow = 0 return end subroutine la_inp1 ( nrow, ncol, a, a_int, iabot, iatop, ierror, base, iform, & line, row1, row2, col1, col2 ) ! !******************************************************************************* ! !! LA_INP1 accepts the values of the entries of a matrix from the user. ! ! ! Modified: ! ! 28 September 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer NROW, NCOL, the number of rows and columns in the matrix. ! ! Output, real A(NROW,NCOL), the current matrix. ! ! Output, integer IATOP(NROW,NCOL), IABOT(NROW,NCOL). ! IATOP and IABOT represent the current rational or decimal ! matrix. ! ! Output, integer IERROR, 0 an error occurred, 1 no error occurred. ! ! Input, integer IFORM, 0/1/2/3/4, rat/real/dec/int/mod arithmetic. ! ! Workspace, character ( len = 80 ) LINE, used to hold the user's input. ! implicit none ! integer ncol integer nrow ! real a(nrow,ncol) integer a_int(nrow,ncol) integer base character ( len = 10 ) chrtmp1 character ( len = 10 ) chrtmp2 character ( len = 10 ) chrtmp3 integer col1 integer col2 integer i integer iabot(nrow,ncol) integer iatop(nrow,ncol) integer ierror integer iform integer j character ( len = 80 ) line character ( len = 80 ) prompt integer row1 integer row2 ! ! Enter a single value ! if ( row1 == row2 .and. col1 == col2 ) then call i_to_s_left ( row1, chrtmp1 ) call i_to_s_left ( col1, chrtmp2 ) prompt = 'A(' // trim ( chrtmp1 ) // ',' // trim ( chrtmp2 ) //')' if ( iform == 0 ) then call rat_read ( iatop(row1,col1), iabot(row1,col1), line, prompt, ierror ) else if ( iform == 1 ) then call r_read ( a(row1,col1), line, prompt, ierror ) else if ( iform == 2 ) then call dec_read ( iatop(row1,col1), iabot(row1,col1), line, prompt, ierror ) call dec_round ( iatop(row1,col1), iabot(row1,col1) ) else if ( iform == 3 ) then call i_read ( a_int(row1,col1), line, prompt, ierror ) else if ( iform == 4 ) then call i_read ( a_int(row1,col1), line, prompt, ierror ) a_int(row1,col1) = mod ( a_int(row1,col1), base ) end if if ( ierror /= 0 ) then return end if ! ! Enter a single column. ! else if ( col1 == col2 ) then do i = row1, row2 call i_to_s_left ( i, chrtmp1 ) call i_to_s_left ( row2, chrtmp2 ) call i_to_s_left ( col1, chrtmp3 ) prompt = 'entries ' // chrtmp1 // ' to ' // chrtmp2 // ' of column ' // & chrtmp3 call s_blanks_delete ( prompt ) if ( iform == 0 ) then call rat_read ( iatop(i,col1), iabot(i,col1), line, prompt, ierror ) else if ( iform == 1 ) then call r_read ( a(i,col1), line, prompt, ierror ) else if ( iform == 2 ) then call dec_read ( iatop(i,col1), iabot(i,col1), line, prompt, ierror ) call dec_round ( iatop(i,col1), iabot(i,col1) ) else if ( iform == 3 ) then call i_read ( a_int(i,col1), line, prompt, ierror ) else if ( iform == 4 ) then call i_read ( a_int(i,col1), line, prompt, ierror ) a_int(i,col1) = mod ( a_int(i,col1), base ) end if if ( ierror /= 0 ) then return end if end do ! ! Enter one or more partial rows of at least 2 entries. ! else do i = row1, row2 do j = col1, col2 call i_to_s_left ( j, chrtmp1 ) call i_to_s_left ( col2, chrtmp2 ) call i_to_s_left ( i, chrtmp3 ) prompt = 'entries ' // chrtmp1 // ' to ' // chrtmp2 // ' of row ' // & chrtmp3 call s_blanks_delete ( prompt ) if ( iform == 0 ) then call rat_read ( iatop(i,j), iabot(i,j), line, prompt, ierror ) else if ( iform == 1 ) then call r_read ( a(i,j), line, prompt, ierror ) else if ( iform == 2 ) then call dec_read ( iatop(i,j), iabot(i,j), line, prompt, ierror ) call dec_round ( iatop(i,j), iabot(i,j) ) else if ( iform == 3 ) then call i_read ( a_int(i,j), line, prompt, ierror ) else if ( iform == 3 ) then call i_read ( a_int(i,j), line, prompt, ierror ) a_int(i,j) = mod ( a_int(i,j), base ) end if if ( ierror /= 0 ) then return end if end do end do end if return end subroutine la_opt ( a, a_int, iabot, iatop, base, iform, nrow, ncol ) ! !******************************************************************************* ! !! LA_OPT checks for row echelon or reduced row echelon form. ! ! ! Discussion: ! ! A matrix is in row echelon form if: ! ! * The first nonzero entry in each row is 1. ! ! * The leading 1 in a given row occurs in a column to ! the right of the leading 1 in the previous row. ! ! * Rows which are entirely zero must occur last. ! ! The matrix is in reduced row echelon form if, in addition to ! the first three conditions, it also satisfies: ! ! * Each column containing a leading 1 has no other nonzero entries. ! ! Modified: ! ! 16 May 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real A(NROW,NCOL), the current matrix. ! ! Input, integer IATOP(NROW,NCOL), IABOT(NROW,NCOL). ! IATOP and IABOT represent the current rational or decimal ! matrix. ! ! Input, integer IFORM, 0/1/2/3/4, rat/real/dec/int/mod arithmetic. ! ! Input, integer NROW, NCOL, the number of rows and columns in the matrix. ! implicit none ! integer ncol integer nrow ! real a(nrow,ncol) integer a_int(nrow,ncol) integer base character ( len = 22 ) chrtmp character ( len = 10 ) chrtmp1 character ( len = 10 ) chrtmp2 integer i integer iabot(nrow,ncol) integer iatop(nrow,ncol) integer iform integer ii integer izer integer j integer lead integer leadp ! write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'Checking the matrix for row echelon form...' ! ! Check rule 1. ! do i = 1, nrow do j = 1, ncol if ( iform == 0 ) then if ( iatop(i,j) == 0 ) then cycle else if ( iatop(i,j) == iabot(i,j) ) then exit end if else if ( iform == 1 ) then if ( a(i,j) == 0.0E+00 ) then cycle else if ( a(i,j) == 1.0E+00 ) then exit end if else if ( iform == 2 ) then if ( iatop(i,j) == 0 ) then cycle else if ( iatop(i,j) == 1 .and. iabot(i,j) == 0 ) then exit end if else if ( iform == 3 ) then if ( a_int(i,j) == 0 ) then cycle else if ( a_int(i,j) == 1 ) then exit end if else if ( iform == 4 ) then if ( a_int(i,j) == 0 ) then cycle else if ( a_int(i,j) == 1 ) then exit end if end if write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' The matrix is NOT in row echelon form.' call i_to_s_left ( i, chrtmp1 ) write ( *, '(a)' ) ' The first nonzero entry in row ' // chrtmp1 call i_to_s_left ( j, chrtmp1 ) write ( *, '(a)' ) ' which occurs in column ' // chrtmp1 if ( iform == 0 ) then call rat_to_s_left ( iatop(i,j), iabot(i,j), chrtmp ) else if ( iform == 1 ) then call r_to_s_left ( a(i,j), chrtmp ) else if ( iform == 2 ) then call dec_to_s_left ( iatop(i,j), iabot(i,j), chrtmp ) else if ( iform == 3 ) then call i_to_s_left ( a_int(i,j), chrtmp ) else if ( iform == 4 ) then call i_to_s_left ( a_int(i,j), chrtmp ) end if write ( *, '(a)' ) ' is ' // trim ( chrtmp ) // ' rather than 1.' return end do end do ! ! Check rule 2. ! lead = 0 do i = 1, nrow do j = 1, ncol if ( iform == 0 ) then if ( iatop(i,j) == 0 ) then cycle else if ( iatop(i,j) == iabot(i,j) ) then leadp = lead lead = j if ( lead > leadp ) then exit end if end if else if ( iform == 1 ) then if ( a(i,j) == 0.0E+00 ) then cycle else if ( a(i,j) == 1.0E+00 ) then leadp = lead lead = j if ( lead > leadp ) then exit end if end if else if ( iform == 2 ) then if ( iatop(i,j) == 0 ) then cycle else if ( iatop(i,j) == 1 .and. iabot(i,j) == 0 ) then leadp = lead lead = j if ( lead > leadp ) then exit end if end if else if ( iform == 3 ) then if ( a_int(i,j) == 0 ) then cycle else if ( a_int(i,j) == 1 ) then leadp = lead lead = j if ( lead > leadp ) then exit end if end if else if ( iform == 4 ) then if ( a_int(i,j) == 0 ) then cycle else if ( a_int(i,j) == 1 ) then leadp = lead lead = j if ( lead > leadp ) then exit end if end if end if write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' The matrix is NOT in row echelon form.' call i_to_s_left ( i, chrtmp1 ) write ( *, '(a)' ) ' The first 1 in row ' // trim ( chrtmp1 ) // ' does' call i_to_s_left ( i-1, chrtmp1 ) write ( *, '(a)' ) ' NOT occur to the right of the first 1 in row ' // & trim ( chrtmp1 ) // '.' return end do end do ! ! Check rule 3. ! izer = 0 do i = 1, nrow if ( izer == 0 ) then izer = i do j = 1, ncol if ( ( iform == 0 .and. iatop(i,j) /= 0 ) .or. & ( iform == 1 .and. a(i,j) /= 0.0E+00 ) .or. & ( iform == 2 .and. iatop(i,j) /= 0 ) .or. & ( iform == 3 .and. a_int(i,j) /= 0 ) .or. & ( iform == 4 .and. a_int(i,j) /= 0 ) ) then izer = 0 exit end if end do else do j = 1, ncol if ( ( iform == 0 .and. iatop(i,j) /= 0 ) .or. & ( iform == 1 .and. a(i,j) /= 0.0E+00 ) .or. & ( iform == 2 .and. iatop(i,j) /= 0 ) .or. & ( iform == 3 .and. a_int(i,j) /= 0 ) .or. & ( iform == 4 .and. a_int(i,j) /= 0 ) ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' The matrix is NOT in row echelon form.' call i_to_s_left ( izer, chrtmp1 ) write ( *, '(a)' ) ' Row ' // trim ( chrtmp1 ) // ' is entirely zero.' call i_to_s_left ( i, chrtmp1 ) write ( *, '(a)' ) ' Row ' // trim ( chrtmp1 ) // ' occurs later, and has' write ( *, '(a)' ) ' nonzero entries in it!' return end if end do end if end do write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' The matrix is in row echelon form.' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Now checking to see if the matrix is in ' write ( *, '(a)' ) ' REDUCED row echelon form...' ! ! Check rule 4. ! do i = 1, nrow do j = 1, ncol ! ! We know first nonzero in this row will be 1. ! if ( ( iform == 0 .and. iatop(i,j) == 0 ) .or. & ( iform == 1 .and. a(i,j) == 0.0E+00 ) .or. & ( iform == 2 .and. iatop(i,j) == 0 ) .or. & ( iform == 3 .and. a_int(i,j) == 0 ) .or. & ( iform == 4 .and. a_int(i,j) == 0 ) ) then cycle end if ! ! The leading 1 of this row is entry (i,j). ! do ii = 1, nrow if ( ii /= i ) then if ( ( iform == 0 .and. iatop(ii,j) == 0 ) .or. & ( iform == 1 .and. a(ii,j) == 0.0E+00 ) .or. & ( iform == 2 .and. iatop(ii,j) == 0 ) .or. & ( iform == 3 .and. a_int(ii,j) == 0 ) .or. & ( iform == 4 .and. a_int(ii,j) == 0 ) ) then cycle end if write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' The matrix is NOT in reduced row echelon form.' call i_to_s_left ( i, chrtmp1 ) call i_to_s_left ( j, chrtmp2 ) write ( *, '(a)' ) ' Row ' // trim ( chrtmp1 ) // & ' has its leading 1 in column ' // trim ( chrtmp2 ) // '.' write ( *, '(a)' ) ' This means that all other entries of that ' // & 'column should be zero.' if ( iform == 0 ) then call rat_to_s_left ( iatop(ii,j), iabot(ii,j), chrtmp ) else if ( iform == 1 ) then call r_to_s_left ( a(ii,j), chrtmp ) else if ( iform == 2 ) then call dec_to_s_left ( iatop(ii,j), iabot(ii,j), chrtmp ) else if ( iform == 3 ) then call i_to_s_left ( a_int(ii,j), chrtmp ) else if ( iform == 4 ) then call i_to_s_left ( a_int(ii,j), chrtmp ) end if call i_to_s_left ( ii, chrtmp1 ) write ( *, '(a)' ) ' But the entry in row ' // trim ( chrtmp1 ) & // ' is ' // trim ( chrtmp ) return end if end do exit end do end do write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'The matrix is in reduced row echelon form.' return end subroutine mat_append_identity ( nrow, ncol, a, a_int, iatop, iabot, iform ) ! !******************************************************************************* ! !! MAT_APPEND_IDENTITY appends the identity matrix. ! ! ! Modified: ! ! 29 September 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, integer NROW, NCOL, the number of rows and columns of ! the matrix. On output, NCOL has been increased to NCOL+NROW. ! ! Input/output, real A(NROW,NCOL), the current matrix. ! ! Input/output, integer IATOP(NROW,NCOL), IABOT(NROW,NCOL). ! IATOP and IABOT represent the current rational or decimal ! matrix. ! ! Input, integer IFORM, 0/1/2/3/4, rat/real/dec/int/mod arithmetic. ! implicit none ! integer nrow ! real a(nrow,*) integer a_int(nrow,*) integer i integer iabot(nrow,*) integer iatop(nrow,*) integer iform integer j integer jform integer ncol ! do i = 1, nrow do j = 1, nrow if ( i == j ) then if ( iform == 0 ) then iatop(i,ncol+j) = 1 iabot(i,ncol+j) = 1 else if ( iform == 1 ) then a(i,ncol+j) = 1.0E+00 else if ( iform == 2 ) then iatop(i,ncol+j) = 1 iabot(i,ncol+j) = 0 else if ( iform == 3 ) then a_int(i,ncol+j) = 1 end if else if ( iform == 0 ) then iatop(i,ncol+j) = 0 iabot(i,ncol+j) = 1 else if ( iform == 1 ) then a(i,ncol+j) = 0.0E+00 else if ( iform == 2 ) then iatop(i,ncol+j) = 0 iabot(i,ncol+j) = 0 else if ( iform == 3 ) then a_int(i,ncol+j) = 0 end if end if end do end do ncol = ncol + nrow return end subroutine mat_print ( a, a_int, iabot, iatop, iform, nrow, ncol, title ) ! !******************************************************************************* ! !! MAT_PRINT prints out the matrix. ! ! ! Modified: ! ! 09 October 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real A(NROW,NCOL), the current matrix. ! ! Input, integer IABOT(NROW,NCOL), IATOP(NROW,NCOL). ! IATOP and IABOT represent the current rational or decimal matrix. ! ! Input, integer IFORM, 0/1/2/3/4, rat/real/dec/int/mod arithmetic. ! ! Input, integer NROW, NCOL, the number of rows and columns in the matrix. ! ! Input, character ( len = * ) TITLE, the title of the object to be printed. ! implicit none ! integer ncol integer nrow ! real a(nrow,ncol) integer a_int(nrow,ncol) integer iabot(nrow,ncol) integer iatop(nrow,ncol) integer iform character ( len = * ) title ! if ( iform == 0 ) then call rat_print ( iatop, iabot, nrow, ncol, title ) else if ( iform == 1 ) then call r_print ( a, nrow, ncol, title ) else if ( iform == 2 ) then call dec_print ( iatop, iabot, nrow, ncol, title ) else if ( iform == 3 ) then call i_print ( a_int, nrow, ncol, title ) else if ( iform == 4 ) then call i_print ( a_int, nrow, ncol, title ) end if return end subroutine r_print ( a, nrow, ncol, title ) ! !******************************************************************************* ! !! R_PRINT prints out real vectors and matrices. ! ! ! Modified: ! ! 26 September 1998 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real A(NROW,NCOL), the current matrix. ! ! Input, integer NROW, NCOL, the number of rows and columns in the matrix. ! ! Input, character ( len = * ) TITLE, a label for the object being printed. ! implicit none ! integer ncol integer nrow ! real a(nrow,ncol) logical allint character ( len = 10 ) chrtmp1 character ( len = 10 ) chrtmp2 character ( len = 10 ) chrtmp3 character ( len = 40 ) format1 character ( len = 40 ) format2 integer i integer ichi integer iclo integer ihi integer ilo integer imax integer imin integer izhi integer izlo integer j integer jhi integer jlo integer jmax integer jmin integer kmax integer kmin integer, parameter :: ncolum = 80 integer npline character ( len = 100 ) output character ( len = * ) title ! ! Figure out how many numbers we can fit in NCOLUM columns. ! kmin = 11 kmax = kmin do i = 1, nrow do j = 1, ncol do while ( abs ( a(i,j) ) >= 10.0E+00**(kmax-kmin) ) kmax = kmax + 1 end do end do end do npline = ncolum / kmax ! ! Check to see if the matrix entries are all integers. ! allint = .true. do i = 1, nrow do j = 1, ncol if ( a(i,j) /= real ( int ( a(i,j) ) ) ) then allint = .false. end if end do end do ! ! If all integers, cut down KMAX, the width of each number, ! and update NPLINE, the number of numbers we can print on one line. ! if ( allint ) then kmax = kmax - 7 npline = ncolum / kmax call i_to_s_left ( npline, chrtmp2 ) call i_to_s_left ( kmax, chrtmp3 ) format1 = '(' // chrtmp2 // 'f' // chrtmp3 // '.0)' ! ! If nonintegral entries, print 7 decimals. ! else call i_to_s_left ( npline, chrtmp2 ) call i_to_s_left ( kmax, chrtmp3 ) format1 = '(' // chrtmp2 // 'f' // chrtmp3 // '.7)' end if call s_blank_delete ( format1 ) ! ! The second format is for ... ! format2 = '(' // chrtmp2 // 'i' // chrtmp3 // ')' call s_blank_delete ( format2 ) ! ! Now print the data. ! do jmin = 1, ncol, npline jmax = min ( jmin + npline - 1, ncol ) write ( *, '(a)' ) ' ' if ( jmin == 1 ) then write ( *, '(a)' ) trim ( title ) write ( *, '(a)' ) ' ' end if if ( jmin > 1 .or. jmax < ncol ) then write ( output, * ) 'Columns ', jmin, ' to ', jmax call s_blanks_delete ( output ) write ( *, '(a)' ) trim ( output ) write ( *, '(a)' ) ' ' end if do i = 1, nrow write ( *, format1 ) a(i,jmin:jmax) end do end do return end subroutine r_read ( rval, line, prompt, ierror ) ! !******************************************************************************* ! !! R_READ "reads" a real value from a line of text. ! ! ! Discussion: ! ! The routine 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. ! ! The routine will accept integers, decimals, and ratios of the ! form R1/R2. Real numbers may be in scientific notation, as ! +12.34E-56.78 ! ! Modified: ! ! 12 July 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, real RVAL, the real value found in LINE. ! ! Workspace, character ( len = 80 ) LINE, used to hold the user's input. ! ! Workspace, character ( len = 80 ) PROMPT, the prompt string. ! ! Output, integer IERROR, 0 an error occurred, 1 no error occurred. ! implicit none ! real bot integer ierror integer lchar character ( len = 80 ) line character ( len = 80 ) prompt real rval real top ! ierror = 0 rval = 0.0E+00 top = 0.0E+00 bot = 1.0E+00 ! ! Read a character string. ! do call chrinp ( ierror, line, prompt ) if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'R_READ - Fatal error!' write ( *, '(a)' ) ' CHRINP returned error flag!' return end if if ( line /= ' ' ) then exit end if end do ! ! Convert the character string to a decimal value, TOP. ! call s_to_r ( line, top, ierror, lchar ) ! ! If we haven't used up all our characters, ! and if the next character is '/', ! then the user means to input the value as a ratio, ! so prepare to read BOT as well. ! if ( lchar+1 < len ( line ) ) then if ( line(lchar+1:lchar+1) == '/' ) then lchar = lchar + 1 call s_chop ( line, 1, lchar ) call s_to_r ( line, bot, ierror, lchar ) if ( bot == 0.0E+00 ) then bot = 1.0E+00 end if end if end if ! ! Set the value of RVAL. ! rval = top / bot ! ! Chop out the characters that were used. ! if ( lchar < len ( line ) ) then if ( line(lchar+1:lchar+1) == ',' ) then lchar = lchar + 1 end if end if call s_chop ( line, 1, lchar ) return end subroutine r_swap ( x, y ) ! !******************************************************************************* ! !! R_SWAP switches two real values. ! ! ! 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 z real x real y ! z = x x = y y = z return end subroutine r_to_dec ( rval, itop, ibot ) ! !******************************************************************************* ! !! R_TO_DEC converts a real value to a decimal fraction form. ! ! ! Discussion: ! ! The routine is given RVAL, and computes ITOP and IBOT, so that ! approximately: ! ! RVAL = ITOP * 10 ** IBOT ! ! However, only DEC_DIGIT digits of RVAL are used in constructing the ! representation, where DEC_DIGIT is the maximum number of decimal ! digits used in the decimal representation. ! ! Modified: ! ! 12 July 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real RVAL, the real number to be converted. ! ! Output, integer ITOP, IBOT, the approximate decimal representation. ! ! ITOP is an integer, strictly between -10**DEC_DIGIT and 10**DEC_DIGIT. ! IBOT is an integer exponent of 10. ! implicit none ! integer dec_digit integer ibot integer itop real rtop real rval real ten1 real ten2 ! ! Special cases. ! if ( rval == 0.0E+00 ) then itop = 0 ibot = 0 return end if dec_digit = 0 call i_data ( 'GET', 'DEC_DIGIT', dec_digit ) ! ! Factor RVAL = RTOP * 10**IBOT ! rtop = rval ibot = 0 ! ! Now normalize so that 10**(DEC_DIGIT-1) <= ABS(RTOP) < 10**(DEC_DIGIT) ! ten1 = 10.0E+00**( dec_digit - 1 ) ten2 = 10.0E+00**dec_digit do while ( abs ( rtop ) < ten1 ) rtop = rtop * 10.0E+00 ibot = ibot - 1 end do do while ( abs ( rtop ) >= ten2 ) rtop = rtop / 10.0E+00 ibot = ibot + 1 end do ! ! ITOP is the integer part of RTOP, rounded. ! itop = nint ( rtop ) ! ! Now divide out any factors of ten from ITOP. ! if ( itop /= 0 ) then do while ( mod ( itop, 10 ) == 0 ) itop = itop / 10 ibot = ibot + 1 end do end if return end subroutine r_to_rat ( a, iatop, iabot ) ! !******************************************************************************* ! !! R_TO_RAT converts a real value to a rational value. ! ! ! Discussion: ! ! The rational value (IATOP/IABOT) is essentially computed by truncating ! the decimal representation of the real value after a given number of ! decimal digits. ! ! Modified: ! ! 12 July 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real A, the real value to be converted. ! ! Output, integer IATOP, IABOT, the numerator and denominator ! of the rational value that approximates A. ! implicit none ! real a integer dec_digit real factor integer i_gcd integer iabot integer iatop integer ibot integer ifac integer itemp integer itop integer jfac ! dec_digit = 0 call i_data ( 'GET', 'DEC_DIGIT', dec_digit ) factor = 10.0E+00**dec_digit if ( dec_digit > 0 ) then ifac = 10**dec_digit jfac = 1 else ifac = 1 jfac = 10**( - dec_digit ) end if itop = nint ( a * factor ) * jfac ibot = ifac ! ! Factor out the greatest common factor. ! itemp = i_gcd ( itop, ibot ) iatop = itop / itemp iabot = ibot / itemp return end subroutine r_to_s_left ( rval, s ) ! !******************************************************************************* ! !! R_TO_S_LEFT represents a real using 14 left_justified characters. ! ! ! Modified: ! ! 22 January 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real RVAL, a real number. ! ! Output, character ( len = * ) S, a left-justified character variable ! containing the representation of RVAL. ! implicit none ! character ( len = 14 ) chrtmp integer i real rval character ( len = * ) s ! ! We can't seem to write directly into the string because of compiler ! quibbles. ! if ( real ( int ( rval ) ) == rval .and. abs ( rval ) < 1.0E+13 ) then write ( chrtmp, '(i14)' ) int ( rval ) else write ( chrtmp, '(g14.6)' ) rval end if do i = 1, len ( chrtmp ) if ( chrtmp(i:i) /= ' ' ) then s = chrtmp(i:) return end if end do s = ' ' return end subroutine rat_add ( ibot, ibot1, ibot2, itop, itop1, itop2 ) ! !******************************************************************************* ! !! RAT_ADD adds two rational values. ! ! ! Formula: ! ! ITOP / IBOT = ( ITOP1 / IBOT1 ) + ( ITOP2 / IBOT2 ) ! ! Note: ! ! If numeric overflow would occur, the computation is done in ! real arithmetic and then converted back to a rational value. ! ! Modified: ! ! 12 September 1998 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer IBOT, the denominator of the result. ! ! Input, integer IBOT1, IBOT2, the denominators of the ! two rational values to be added. ! ! Output, integer ITOP, the numerator of the result. ! ! Input, integer ITOP1, ITOP2, the numerators of the ! two rational values to be added. ! implicit none ! integer i_big integer i_gcd integer ibot integer ibot1 integer ibot2 integer ierror integer itemp integer itop integer itop1 integer itop2 integer jbot1 integer jbot2 integer jbot3 integer jtop1 integer jtop2 real rbot real rmax real rtop1 real rtop2 real rtop3 real rval real rval1 real rval2 ! ierror = 0 i_big = 0 call i_data ( 'GET', 'I_BIG', i_big ) rmax = real ( i_big ) if ( itop1 == 0 ) then itop = itop2 ibot = ibot2 return else if ( itop2 == 0 ) then itop = itop1 ibot = ibot1 return end if ! ! Make copies of the input arguments, since we will change them. ! jbot1 = ibot1 jbot2 = ibot2 jtop1 = itop1 jtop2 = itop2 ! ! Compute the greatest common factor of the two denominators, ! and factor it out. ! jbot3 = i_gcd ( jbot1, jbot2 ) jbot1 = jbot1 / jbot3 jbot2 = jbot2 / jbot3 ! ! The fraction may now be formally written as: ! ! (jtop1*jbot2 + jtop2*jbot1) / (jbot1*jbot2*jbot3) ! ! Check the tops for overflow. ! rtop1 = real ( jtop1 ) * real ( jbot2 ) if ( abs ( rtop1 ) > rmax ) then ierror = 1 itop = 0 else jtop1 = jtop1 * jbot2 end if rtop2 = real ( jtop2 ) * real ( jbot1 ) if ( abs ( rtop2 ) > rmax ) then ierror = 2 itop = 0 else jtop2 = jtop2 * jbot1 end if rtop3 = real ( jtop1 ) + real ( jtop2 ) if ( abs ( rtop3 ) > rmax ) then ierror = 3 itop = 0 else itop = jtop1 + jtop2 end if ! ! Check the bottom for overflow. ! rbot = real ( jbot1 ) * real ( jbot2 ) * real ( jbot3 ) if ( abs ( rbot ) > rmax ) then ierror = 4 ibot = 1 else ibot = jbot1 * jbot2 * jbot3 end if ! ! If there was potential overflow, then do the computation in ! real arithmetic and convert back. ! if ( ierror /= 0 ) then ierror = 0 rval1 = real ( itop1 ) / real ( ibot1 ) rval2 = real ( itop2 ) / real ( ibot2 ) rval = rval1 + rval2 call r_to_rat ( rval, itop, ibot ) end if ! ! Put the fraction in lowest terms. ! itemp = i_gcd ( itop, ibot ) itop = itop / itemp ibot = ibot / itemp ! ! Sign of bottom should be positive. ! if ( ibot < 0 ) then ibot = - ibot itop = - itop end if return end subroutine rat_mul ( ibot, ibot1, ibot2, itop, itop1, itop2, ierror ) ! !******************************************************************************* ! !! RAT_MUL multiplies two fractions. ! ! ! Formula: ! ! ITOP / IBOT = ( ITOP1 / IBOT1 ) * ( ITOP2 / IBOT2 ). ! ! Note: ! ! If numeric overflow would occur, the computation is done in ! real arithmetic and then converted back to a rational value. ! ! Modified: ! ! 12 September 1998 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer IBOT, the denominator of the result. ! ! Input, integer IBOT1, IBOT2, the denominators of the ! two rational values to be multiplied. ! ! Output, integer ITOP, the numerator of the result. ! ! Input, integer ITOP1, ITOP2, the numerators of the ! two rational values to be multiplied. ! ! Output, integer IERROR, 0 an error occurred, 1 no error occurred. ! implicit none ! integer i_big integer i_gcd integer ibot integer ibot1 integer ibot2 integer ierror integer itemp integer itop integer itop1 integer itop2 integer jbot1 integer jbot2 integer jtop1 integer jtop2 real rbot real rmax real rtop real temp ! ierror = 0 i_big = 0 call i_data ( 'GET', 'I_BIG', i_big ) rmax = real ( i_big ) if ( itop1 == 0 .or. itop2 == 0 ) then itop = 0 ibot = 1 return end if if ( ibot1 == 0 .or. ibot2 == 0 ) then ierror = 1 itop = 0 ibot = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'RAT_MUL - Fatal error!' write ( *, '(a)' ) ' A rational fraction has a zero denominator!' return end if ! ! Make copies of the input arguments, since we will change them. ! jbot1 = ibot1 jbot2 = ibot2 jtop1 = itop1 jtop2 = itop2 ! ! Get rid of all common factors in top and bottom. ! itemp = i_gcd ( jtop1, jbot1 ) jtop1 = jtop1 / itemp jbot1 = jbot1 / itemp itemp = i_gcd ( jtop1, jbot2 ) jtop1 = jtop1 / itemp jbot2 = jbot2 / itemp itemp = i_gcd ( jtop2, jbot1 ) jtop2 = jtop2 / itemp jbot1 = jbot1 / itemp itemp = i_gcd ( jtop2, jbot2 ) jtop2 = jtop2 / itemp jbot2 = jbot2 / itemp ! ! The fraction (ITOP1*ITOP2)/(IBOT1*IBOT2) is in lowest terms. ! ! Check the top ITOP1*ITOP2 for overflow. ! rtop = real ( jtop1 ) * real ( jtop2 ) if ( abs ( rtop ) > rmax ) then ierror = 1 else itop = jtop1 * jtop2 end if ! ! Check the bottom IBOT1*IBOT2 for overflow. ! rbot = real ( jbot1 ) * real ( jbot2 ) if ( abs ( rbot ) > rmax ) then ierror = 2 else ibot = jbot1 * jbot2 end if ! ! if there was an overflow, then compute RTOP / RBOT and convert ! back to rational. ! if ( ierror == 1 .or. ierror == 2 ) then ierror = 0 temp = rtop / rbot call r_to_rat ( temp, itop, ibot ) end if ! ! Sign of bottom should be positive. ! if ( ibot < 0 ) then ibot = - ibot itop = - itop end if ! ! The fraction is ITOP/IBOT with no loss of accuracy (unless there ! was overflow). ! return end subroutine rat_print ( iatop, iabot, nrow, ncol, title ) ! !******************************************************************************* ! !! RAT_PRINT prints out rational vectors or matrices. ! ! ! Modified: ! ! 12 July 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer IATOP(NROW,NCOL), IABOT(NROW,NCOL). ! IATOP and IABOT represent the current rational or decimal ! matrix. ! ! Input, integer NROW, NCOL, the number of rows and columns in the matrix. ! ! Input, character ( len = * ) TITLE, a label for the object being printed. ! implicit none ! integer ncol integer nrow ! character ( len = 10 ) chrtmp1 character ( len = 10 ) chrtmp2 character ( len = 10 ) chrtmp3 character ( len = 40 ) format1 character ( len = 40 ) format2 integer i integer iabot(nrow,ncol) integer iatop(nrow,ncol) integer ichi integer iclo integer imax integer imin integer ione integer itemp integer izhi integer izlo integer j integer jmax integer jmin integer kmax integer, parameter :: ncolum = 80 integer none integer npline character ( len = 100 ) output character ( len = * ) title ! ! Figure out how many rationals we can get in NCOLUM columns. ! kmax = 3 do i = 1, nrow do j = 1, ncol itemp = abs ( iatop(i,j) ) do while ( itemp >= 10**(kmax-2) ) kmax = kmax + 1 end do itemp = abs ( iabot(i,j) ) do while ( itemp > 10**(kmax-2) ) kmax = kmax + 1 end do end do end do kmax = kmax + 1 npline = ncolum / kmax ! ! Create the formats. ! call i_to_s_left ( npline, chrtmp2 ) call i_to_s_left ( kmax, chrtmp3 ) format1 = '(' // chrtmp2 // 'i' // chrtmp3 // ')' call s_blank_delete ( format1 ) format2 = '(' // chrtmp2 // 'i' // chrtmp3 // ')' call s_blank_delete ( format2 ) do jmin = 1, ncol, npline jmax = min ( jmin+npline-1, ncol ) write ( *, '(a)' ) ' ' if ( jmin == 1 ) then write ( *, '(a)' ) trim ( title ) write ( *, '(a)' ) ' ' end if if ( jmin > 1 .or. jmax < ncol ) then write ( output, * ) 'Columns ', jmin, ' to ', jmax call s_blanks_delete ( output ) write ( *, '(a)' ) trim ( output ) write ( *, '(a)' ) ' ' end if do i = 1, nrow write ( *, format1 ) iatop(i,jmin:jmax) write ( output, format1 ) iabot(i,jmin:jmax) ! ! Delete each denominator that is 1. If all are 1, don't ! even print out the line. ! none = 0 do j = jmin, jmax if ( iabot(i,j) == 1 ) then ione = (j-jmin+1) * kmax output(ione:ione) = ' ' else none = 1 end if end do write ( *, '(a)' ) trim ( output ) if ( jmax == ncol .and. i == nrow ) then else write ( *, '(a)' ) ' ' end if end do end do return end subroutine rat_read ( itop, ibot, line, prompt, ierror ) ! !******************************************************************************* ! !! RAT_READ reads a rational value, expressed as integer, decimal or fraction. ! ! ! Modified: ! ! 22 January 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer ITOP, IBOT, the top and bottom of the ! fraction that was read. ! ! Workspace, character ( len = 80 ) LINE, used to hold the user's input. ! ! Workspace, character ( len = 80 ) PROMPT. ! ! Output, integer IERROR, 0 an error occurred, 1 no error occurred. ! implicit none ! integer i_gcd integer ibot integer ibot1 integer ibot2 integer ierror integer itemp integer itop integer itop1 integer itop2 integer lchar character ( len = 80 ) line character ( len = 80 ) prompt ! ierror = 0 itop = 0 ibot = 1 do call chrinp ( ierror, line, prompt ) if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'RAT_READ - Fatal error!' write ( *, '(a)' ) ' CHRINP returned error flag!' return end if if ( line /= ' ' ) then exit end if end do call chrctf ( line, itop1, ibot1, ierror, lchar ) if ( lchar >= len ( line ) ) then itop = itop1 ibot = ibot1 else if ( line(lchar+1:lchar+1) /= '/' ) then itop = itop1 ibot = ibot1 else lchar = lchar + 1 call s_chop ( line, 1, lchar ) call chrctf ( line, itop2, ibot2, ierror, lchar ) itop = itop1 * ibot2 ibot = ibot1 * itop2 end if call s_chop ( line, 1, lchar ) ! ! Make sure fraction is in lowest terms. ! itemp = i_gcd ( itop, ibot ) itop = itop / itemp ibot = ibot / itemp return end subroutine rat_to_dec ( iatop, iabot, ierror ) ! !******************************************************************************* ! !! RAT_TO_DEC converts a rational value to a decimal value. ! ! ! Modified: ! ! 12 September 1998 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, integer IATOP, IABOT. ! ! On input, the rational value (IATOP/IABOT) to be converted. ! ! On output, the rational decimal value IATOP * 10**IABOT. ! ! Output, integer IERROR, 0 an error occurred, 1 no error occurred. ! implicit none ! double precision dval integer iabot integer iatop integer ierror ! if ( iabot /= 0 ) then ierror = 0 dval = dble ( iatop ) / dble ( iabot ) call d_to_dec ( dval, iatop, iabot ) else ierror = 1 iatop = 0 iabot = 0 end if return end subroutine rat_to_r ( iatop, iabot, a ) ! !******************************************************************************* ! !! RAT_TO_R converts rational values to real values. ! ! ! Modified: ! ! 12 July 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer IATOP, IABOT, the rational quantity ! (IATOP/IABOT) that is to be converted. ! ! Output, real A, the value of the rational quantity. ! implicit none ! real a integer iabot integer iatop ! if ( iabot == 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'RATREL - Warning!' write ( *, '(a)' ) ' The input fraction had a zero denominator.' a = 0.0E+00 else a = real ( iatop ) / real ( iabot ) end if return end subroutine rat_to_s_left ( ival, jval, string ) ! !******************************************************************************* ! !! RAT_TO_S_LEFT returns a left-justified representation of IVAL/JVAL. ! ! ! Discussion: ! ! If the ratio is negative, a minus sign precedes IVAL. ! A slash separates IVAL and JVAL. ! ! Modified: ! ! 13 September 1998 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer IVAL, JVAL, the two integers whose ! ratio IVAL/JVAL is to be represented. ! ! Note that if IVAL is nonzero and JVAL is 0, STRING will ! be returned as "Inf" or "-Inf" (Infinity), and if both ! IVAL and JVAL are zero, STRING will be returned as "NaN" ! (Not-a-Number). ! ! Output, character ( len = 22 ) STRING, a left-justified string ! containing the representation of IVAL/JVAL. ! implicit none ! integer ival integer ival2 integer jval integer jval2 character ( len = * ) string ! ! Take care of simple cases right away. ! if ( ival == 0 ) then if ( jval /= 0 ) then string = '0' else string = 'NaN' end if else if ( jval == 0 ) then if ( ival > 0 ) then string = 'Inf' else string = '-Inf' end if ! ! Make copies of IVAL and JVAL. ! else ival2 = ival jval2 = jval if ( jval2 == 1 ) then write ( string, '(i11)' ) ival2 else write ( string, '(i11, ''/'', i10)' ) ival2, jval2 end if call s_blank_delete ( string ) end if return end subroutine row_add ( a, a_int, iatop, iabot, ierror, base, iform, irow1, & irow2, nrow, ncol, sval, istop, isbot ) ! !******************************************************************************* ! !! ROW_ADD adds a multiple of one row to another. ! ! ! Modified: ! ! 26 September 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, real A(NROW,NCOL), the current matrix. ! ! Input/output, integer IATOP(NROW,NCOL), IABOT(NROW,NCOL). ! IATOP and IABOT represent the current rational or decimal matrix. ! ! Output, integer IERROR, 0 for no error, 1 for an error. ! ! Input, integer IFORM, 0/1/2/3/4, rat/real/dec/int/mod arithmetic. ! ! Input, integer IROW1, the row which is to be modified. ! ! Input, integer IROW2, the row which is to be multiplied by ! a given value and added to row IROW1. ! ! Input, integer NROW, NCOL, the number of rows and columns in the matrix. ! ! Input, real SVAL, the multiplier to use if real arithmetic is employed. ! ! Input, integer ISTOP, ISBOT, the fractional or decimal ! multiplier to use. ! implicit none ! integer ncol integer nrow ! real a(nrow,ncol) integer a_int(nrow,ncol) integer base character ( len = 10 ) chrtmp1 character ( len = 10 ) chrtmp2 character ( len = 22 ) chrtmp3 integer iabot(nrow,ncol) integer iatop(nrow,ncol) integer ibot integer ierror integer iform integer irow1 integer irow2 integer isbot integer isbot2 integer istop integer istop2 integer itop integer j real sval ! ! Return immediately if the multiplier is zero. ! if ( ( iform == 0 .and. istop == 0 ) .or. & ( iform == 1 .and. sval == 0.0E+00 ) .or. & ( iform == 2 .and. istop == 0 ) .or. & ( iform == 3 .and. istop == 0 ) .or. & ( iform == 4 .and. istop == 0 ) ) then return end if ! ! Carry out the operation. ! if ( iform == 0 ) then do j = 1, ncol call rat_mul ( isbot2, isbot, iabot(irow2,j), istop2, istop, & iatop(irow2,j), ierror ) if ( ierror /= 0 ) then return end if call rat_add ( ibot, iabot(irow1,j), isbot2, itop, iatop(irow1,j), istop2 ) iatop(irow1,j) = itop iabot(irow1,j) = ibot end do else if ( iform == 1 ) then a(irow1,1:ncol) = a(irow1,1:ncol) + sval * a(irow2,1:ncol) else if ( iform == 2 ) then do j = 1, ncol call dec_mul ( isbot2, isbot, iabot(irow2,j), istop2, istop, & iatop(irow2,j) ) call dec_add ( ibot, iabot(irow1,j), isbot2, itop, & iatop(irow1,j), istop2 ) iatop(irow1,j) = itop iabot(irow1,j) = ibot end do else if ( iform == 3 ) then a_int(irow1,1:ncol) = a_int(irow1,1:ncol) + istop * a_int(irow2,1:ncol) else if ( iform == 4 ) then a_int(irow1,1:ncol) = & mod ( a_int(irow1,1:ncol) + istop * a_int(irow2,1:ncol), base ) end if ! ! Print out a message. ! if ( iform == 0 ) then if ( istop == isbot ) then chrtmp3 = '+' else if ( istop == - isbot ) then chrtmp3 = '-' else call rat_to_s_left ( istop, isbot, chrtmp3 ) end if else if ( iform == 1 ) then if ( sval == 1.0E+00 ) then chrtmp3 = '+' else if ( sval == - 1.0E+00 ) then chrtmp3 = '-' else call r_to_s_left ( sval, chrtmp3 ) end if else if ( iform == 2 ) then if ( istop == 1 .and. isbot == 0 ) then chrtmp3 = '+' else if ( istop == - 1 .and. isbot == 0 ) then chrtmp3 = '-' else call dec_to_s_left ( istop, isbot, chrtmp3 ) end if else if ( iform == 3 ) then if ( istop == 1 ) then chrtmp3 = '+' else if ( istop == - 1 ) then chrtmp3 = '-' else call i_to_s_left ( istop, chrtmp3 ) end if else if ( iform == 4 ) then if ( istop == 1 ) then chrtmp3 = '+' else if ( istop == - 1 ) then chrtmp3 = '-' else call i_to_s_left ( istop, chrtmp3 ) end if end if call i_to_s_left ( irow1, chrtmp1 ) call i_to_s_left ( irow2, chrtmp2 ) if ( chrtmp3 == '-' .or. chrtmp3 == '+' ) then write ( *, '(a)' ) ' ERO: Row ' // trim ( chrtmp1 ) // ' <= ' // 'Row ' // & trim ( chrtmp1 ) // ' ' // chrtmp3(1:1) // ' Row ' // trim ( chrtmp2 ) else if ( chrtmp3(1:1) == '-' .or. chrtmp3(1:1) == '+' ) then write ( *, '(a)' ) ' ERO: Row ' // trim ( chrtmp1 ) // ' <= ' // 'Row ' // & trim ( chrtmp1 ) // ' ' // chrtmp3(1:1) // ' ' // trim ( chrtmp3 ) // & ' Row ' // trim ( chrtmp2 ) else write ( *, '(a)' ) ' ERO: Row ' // trim ( chrtmp1 ) // ' <= ' // 'Row ' // & trim ( chrtmp1 ) // ' + ' // trim ( chrtmp3 ) // ' Row ' // & trim ( chrtmp2 ) end if return end subroutine row_add_param ( ierror, base, iform, irow1, irow2, istop, isbot, & line, nrow, sval ) ! !******************************************************************************* ! !! ROW_ADD_PARAM gets and checks the row add parameters. ! ! ! Modified: ! ! 26 September 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer IERROR, 0 for no error, 1 for an error. ! ! Input, integer IFORM, 0/1/2/3/4, rat/real/dec/int/mod arithmetic. ! ! Input, integer IROW1, the row to which the multiple is to be added. ! ! Input, integer IROW2, the row which is to be multiplied and ! added to another row. ! ! Output, integer ISTOP, ISBOT, the parts of the rational ! or decimal fraction of the multiplier, if that is the ! arithmetic being used. ! ! Workspace, character ( len = 80 ) LINE, used to hold the user's input. ! ! Input, integer NROW, the number of rows in the matrix. ! ! Output, real SVAL, the multiplier, if real arithmetic is used. ! implicit none ! integer base integer ierror integer iform integer irow1 integer irow2 integer isbot integer istop character ( len = 80 ) line integer nrow character ( len = 80 ) prompt real sval ! prompt = 'multiplier S, row I to add, target row J.' ! ! Get the multiplier, SVAL or ISTOP/ISBOT. ! if ( iform == 0 ) then call rat_read ( istop, isbot, line, prompt, ierror ) else if ( iform == 1 ) then call r_read ( sval, line, prompt, ierror ) else if ( iform == 2 ) then call dec_read ( istop, isbot, line, prompt, ierror ) if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'ROW_ADD_PARAM - Error!' write ( *, '(a)' ) ' DEC_READ returns error code.' return end if call dec_round ( istop, isbot ) else if ( iform == 3 ) then call i_read ( istop, line, prompt, ierror ) else if ( iform == 4 ) then call i_read ( istop, line, prompt, ierror ) istop = mod ( istop, base ) end if if ( ierror /= 0 ) then return end if ! ! Get the row to add, IROW2. ! call i_read ( irow2, line, prompt, ierror ) if ( ierror /= 0 ) then ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Error reading row index from line:' write ( *, '(a)' ) '"' // trim ( line ) // '"' return end if if ( irow2 < 1 .or. irow2 > nrow ) then ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Error! Row index was not acceptable!' return end if ! ! Get the row to which we are adding, IROW1. ! call i_read ( irow1, line, prompt, ierror ) if ( ierror /= 0 ) return if ( irow1 < 1 .or. irow1 > nrow ) then ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Error! The row index was not acceptable!' return end if ! ! Make sure the rows are different. ! if ( irow1 == irow2 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Error! The rows should not be the same!' ierror = 1 return end if ierror = 0 return end subroutine row_auto ( a, a_int, iatop, iabot, ierror, base, iform, nrow, ncol ) ! !******************************************************************************* ! !! ROW_AUTO automatically row reduces the current matrix. ! ! ! Modified: ! ! 29 March 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, real A(NROW,NCOL), the current matrix. ! ! Input/output, integer IATOP(NROW,NCOL), IABOT(NROW,NCOL). ! IATOP and IABOT represent the rational or decimal matrix ! to which elementary row operations will be applied. ! ! Output, integer IERROR, 0 for no error, 1 for an error. ! ! Input, integer IFORM, 0/1/2/3/4, rat/real/dec/int/mod arithmetic. ! ! Input, integer NROW, NCOL, the number of rows and columns in the matrix. ! implicit none ! integer ncol integer nrow ! real a(nrow,ncol) integer a_int(nrow,ncol) real amax real atemp integer base integer i integer iabot(nrow,ncol) integer iatop(nrow,ncol) integer ierror integer iform integer imax integer irow integer isbot integer istop integer j integer jcol integer krow integer l integer lrow real sval ! ierror = 0 krow = 0 do i = 1, nrow irow = i do j = 1, ncol jcol = j ! ! In column JCOL, seek the row between IROW and NROW with ! maximum nonzero entry AMAX. ! imax = 0 amax = 0.0E+00 do krow = irow, nrow if ( iform == 0 ) then call rat_to_r ( iatop(krow,jcol), iabot(krow,jcol), atemp ) else if ( iform == 1 ) then atemp = a(krow,jcol) else if ( iform == 2 ) then call dec_to_r ( iatop(krow,jcol), iabot(krow,jcol), atemp ) else if ( iform == 3 ) then atemp = real ( a_int(krow,jcol) ) else if ( iform == 4 ) then atemp = real ( a_int(krow,jcol) ) end if atemp = abs ( atemp ) if ( atemp > amax ) then amax = atemp imax = krow end if end do if ( imax /= 0 ) then krow = imax exit end if end do if ( krow == 0 ) then return end if ! ! Interchange the IROW-th and the pivot rows. ! if ( krow /= irow ) then call row_swap ( a, a_int, iatop, iabot, ierror, iform, & krow, irow, nrow, ncol ) end if ! ! Divide the pivot row by A(IROW,JCOL) so that A(IROW,JCOL) = 1. ! if ( iform == 0 ) then istop = iatop(irow,jcol) isbot = iabot(irow,jcol) else if ( iform == 1 ) then sval = a(irow,jcol) else if ( iform == 2 ) then istop = iatop(irow,jcol) isbot = iabot(irow,jcol) else if ( iform == 3 ) then istop = a_int(irow,jcol) else if ( iform == 4 ) then istop = a_int(irow,jcol) end if call row_div ( a, a_int, iatop, iabot, ierror, base, iform, irow, & nrow, ncol, sval, istop, isbot ) ! ! Annihilate A(L,JCOL) for L not equal to IROW. ! do l = 1, nrow lrow = l if ( lrow /= irow ) then if ( iform == 0 ) then if ( iatop(lrow,jcol) /= 0 ) then istop = - iatop(lrow,jcol) isbot = iabot(lrow,jcol) call row_add ( a, a_int, iatop, iabot, ierror, base, iform, & lrow, irow, nrow, ncol, sval, istop, isbot ) iatop(lrow,jcol) = 0 iabot(lrow,jcol) = 1 end if else if ( iform == 1 ) then if ( a(lrow,jcol) /= 0.0E+00 ) then sval = - a(lrow,jcol) call row_add ( a, a_int, iatop, iabot, ierror, base, iform, & lrow, irow, nrow, ncol, sval, istop, isbot ) a(lrow,jcol) = 0.0E+00 end if else if ( iform == 2 ) then if ( iatop(lrow,jcol) /= 0 ) then istop = - iatop(lrow,jcol) isbot = iabot(lrow,jcol) call row_add ( a, a_int, iatop, iabot, ierror, base, iform, & lrow, irow, nrow, ncol, sval, istop, isbot ) iatop(lrow,jcol) = 0 iabot(lrow,jcol) = 0 end if else if ( iform == 3 ) then if ( a_int(lrow,jcol) /= 0 ) then istop = - a_int(lrow,jcol) call row_add ( a, a_int, iatop, iabot, ierror, base, iform, & lrow, irow, nrow, ncol, sval, istop, isbot ) a_int(lrow,jcol) = 0 end if else if ( iform == 4 ) then if ( a_int(lrow,jcol) /= 0 ) then istop = - a_int(lrow,jcol) call row_add ( a, a_int, iatop, iabot, ierror, base, iform, & lrow, irow, nrow, ncol, sval, istop, isbot ) a_int(lrow,jcol) = 0 end if end if end if end do end do return end subroutine row_div ( a, a_int, iatop, iabot, ierror, base, iform, irow, & nrow, ncol, sval, istop, isbot ) ! !******************************************************************************* ! !! ROW_DIV divides row IROW of the A matrix by a scale factor. ! ! ! Modified: ! ! 23 September 1998 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, real A(NROW,NCOL), the current matrix. ! ! Input/output, integer IATOP(NROW,NCOL), IABOT(NROW,NCOL). ! IATOP and IABOT represent the current rational or decimal ! matrix. ! ! Output, integer IERROR, 0 for no error, 1 for an error. ! ! Input, integer IFORM, 0/1/2/3/4, rat/real/dec/int/mod arithmetic. ! ! Input, integer IROW, the row to be divided. ! ! Input, integer NROW, NCOL, the number of rows and columns in the matrix. ! ! Input, real SVAL, the real divisor. ! ! Input, integer ISTOP, ISBOT, the fractional or decimal divisor. ! implicit none ! integer ncol integer nrow ! real a(nrow,ncol) integer a_int(nrow,ncol) integer base character ( len = 10 ) chrtmp1 character ( len = 24 ) chrtmp2 integer iabot(nrow,ncol) integer iatop(nrow,ncol) integer ibot integer ierror integer iform integer irow integer isbot integer istop integer itop integer j character ( len = 3 ) op real sval ! ! Make sure that the row number is legal. ! if ( irow < 1 .or. irow > nrow ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Error! The row number is out of range!' ierror = 1 return end if ! ! Check for an illegal divisor of 0. ! if ( ( iform == 0 .and. istop == 0 ) .or. & ( iform == 1 .and. sval == 0.0E+00 ) .or. & ( iform == 2 .and. istop == 0 ) .or. & ( iform == 3 .and. istop == 0 ) .or. & ( iform == 4 .and. istop == 0 ) ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'ROW_DIV - Error!' write ( *, '(a)' ) ' Requested division by zero.' ierror = 1 return end if ! ! Check for a pointless divisor of 1. ! if ( ( iform == 0 .and. istop == isbot ) .or. & ( iform == 1 .and. sval == 1.0E+00 ) .or.& ( iform == 2 .and. istop == 1 .and. isbot == 0 ) .or. & ( iform == 3 .and. istop == 1 ) .or. & ( iform == 4 .and. istop == 1 ) ) then return end if ! ! Carry out the division. ! if ( iform == 0 ) then do j = 1, ncol call rat_mul ( ibot, iabot(irow,j), istop, itop, iatop(irow,j), isbot, & ierror ) if ( ierror /= 0 ) then return end if iatop(irow,j) = itop iabot(irow,j) = ibot end do else if ( iform == 1 ) then a(irow,1:ncol) = a(irow,1:ncol) / sval else if ( iform == 2 ) then do j = 1, ncol call dec_div ( iatop(irow,j), iabot(irow,j), istop, isbot, & iatop(irow,j), iabot(irow,j), ierror ) end do else if ( iform == 3 ) then a_int(irow,1:ncol) = a_int(irow,1:ncol) / istop else if ( iform == 4 ) then a_int(irow,1:ncol) = a_int(irow,1:ncol) / istop end if ! ! Print out a statement about what has been done. ! if ( iform == 0 ) then if ( isbot == 1 ) then call i_to_s_left ( istop, chrtmp2 ) op = ' / ' else call rat_to_s_left ( isbot, istop, chrtmp2 ) op = ' * ' end if else if ( iform == 1 ) then call r_to_s_left ( sval, chrtmp2 ) op = ' / ' else if ( iform == 2 ) then call dec_to_s_left ( istop, isbot, chrtmp2 ) op = ' / ' else if ( iform == 3 ) then call i_to_s_left ( istop, chrtmp2 ) op = ' / ' else if ( iform == 4 ) then call i_to_s_left ( istop, chrtmp2 ) op = ' / ' end if call i_to_s_left ( irow, chrtmp1 ) write ( *, '(a)' ) ' ERO: Row ' // trim ( chrtmp1 ) // ' <= Row ' // & trim ( chrtmp1 ) // op // trim ( chrtmp2 ) return end subroutine row_div_param ( ierror, base, iform, irow, isbot, istop, line, sval ) ! !******************************************************************************* ! !! ROW_DIV_PARAM gets and checks the row divide parameters. ! ! ! Modified: ! ! 20 February 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer IERROR, 0 for no error, 1 for an error. ! ! Input, integer IFORM, 0/1/2/3/4, rat/real/dec/int/mod arithmetic. ! ! Output, integer IROW, the row to be divided. ! ! Workspace, character ( len = 80 ) LINE, used to hold the user's input. ! ! Output, real SVAL, the real divisor. ! implicit none ! integer base integer ierror integer iform integer irow integer isbot integer istop character ( len = 80 ) line character ( len = 80 ) prompt real sval ! prompt = 'row I, divisor S.' ! ! Read the row number to be divided. ! call i_read ( irow, line, prompt, ierror ) if ( ierror /= 0 ) then return end if ! ! Read the divisor, either RVAL or ISTOP/ISBOT. ! if ( iform == 0 ) then call rat_read ( istop, isbot, line, prompt, ierror ) else if ( iform == 1 ) then call r_read ( sval, line, prompt, ierror ) else if ( iform == 2 ) then call dec_read ( istop, isbot, line, prompt, ierror ) if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'ROW_DIV_PARAM - Fatal error!' write ( *, '(a)' ) ' DEC_READ returned error flag.' return end if call dec_round ( istop, isbot ) else if ( iform == 3 ) then call i_read ( istop, line, prompt, ierror ) else if ( iform == 4 ) then call i_read ( istop, line, prompt, ierror ) end if if ( ierror /= 0 ) then return end if return end subroutine row_mul ( a, a_int, iatop, iabot, ierror, base, iform, irow, & nrow, ncol, sval, istop, isbot ) ! !******************************************************************************* ! !! ROW_MUL multiplies a row of the matrix by a scale factor. ! ! ! Modified: ! ! 12 July 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, real A(NROW,NCOL), the current matrix. ! ! Input/output, integer IATOP(NROW,NCOL), IABOT(NROW,NCOL). ! IATOP and IABOT represent the current rational or decimal matrix. ! ! Output, integer IERROR, 0 for no error, 1 for an error. ! ! Input, integer IFORM, 0/1/2/3/4, rat/real/dec/int/mod arithmetic. ! ! Input, integer IROW, the row that is to be multiplied. ! ! Input, integer NROW, NCOL, the number of rows and columns in the matrix. ! ! Input, real SVAL, the real row multiplier. ! ! Input, integer ISTOP, ISBOT, the decimal or fractional row multiplier. ! implicit none ! integer ncol integer nrow ! real a(nrow,ncol) integer a_int(nrow,ncol) integer base character ( len = 10 ) chrtmp1 character ( len = 22 ) chrtmp integer iabot(nrow,ncol) integer iatop(nrow,ncol) integer ibot integer ierror integer iform integer irow integer isbot integer istop integer itop integer j real sval ! ! Make sure row number is OK. ! if ( irow < 1 .or. irow > nrow ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Error! The row number is out of range!' ierror = 1 return end if ! ! For rational arithmetic, make sure bottom of scale factor ! is not 0. ! if ( iform == 0 .and. isbot == 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Error! Illegal 0 divisor in multiplier!' ierror = 1 return end if ! ! Check for multiplication by 0. ! if ( ( iform == 0 .and. istop == 0 ) .or. & ( iform == 1 .and. sval == 0.0E+00 ) .or. & ( iform == 2 .and. istop == 0 ) .or. & ( iform == 3 .and. istop == 0 ) .or. & ( iform == 4 .and. istop == 0 ) ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Warning - Multiplication by zero is not an ERO.' ierror = 1 return end if ! ! Check for multiplication by 1. ! if ( ( iform == 0 .and. istop == isbot ) .or. & ( iform == 1 .and. sval == 1.0E+00 ) .or. & ( iform == 2 .and. istop == 1 .and. isbot == 0 ) .or. & ( iform == 3 .and. istop == 1 ) .or. & ( iform == 4 .and. istop == 1 ) ) then return end if ! ! Carry out the multiplication. ! if ( iform == 0 ) then do j = 1, ncol call rat_mul ( ibot, iabot(irow,j), isbot, itop, iatop(irow,j), istop, & ierror ) if ( ierror /= 0 ) then return end if iatop(irow,j) = itop iabot(irow,j) = ibot end do else if ( iform == 1 ) then a(irow,1:ncol) = sval * a(irow,1:ncol) else if ( iform == 2 ) then do j = 1, ncol call dec_mul ( ibot, iabot(irow,j), isbot, itop, iatop(irow,j), istop ) if ( ierror /= 0 ) return iatop(irow,j) = itop iabot(irow,j) = ibot end do else if ( iform == 3 ) then a_int(irow,1:ncol) = istop * a_int(irow,1:ncol) else if ( iform == 4 ) then a_int(irow,1:ncol) = mod ( istop * a_int(irow,1:ncol), base ) end if ! ! Confirm the operation. ! if ( iform == 0 ) then if ( istop == - isbot ) then chrtmp = '-' else call rat_to_s_left ( istop, isbot, chrtmp ) end if else if ( iform == 1 ) then if ( sval == - 1.0E+00 ) then chrtmp = '-' else call r_to_s_left ( sval, chrtmp ) end if else if ( iform == 2 ) then if ( istop == - 1 .and. isbot == 0 ) then chrtmp = '-' else call dec_to_s_left ( istop, isbot, chrtmp ) end if else if ( iform == 3 ) then if ( istop == - 1 ) then chrtmp = '-' else call i_to_s_left ( istop, chrtmp ) end if else if ( iform == 4 ) then if ( istop == - 1 ) then chrtmp = '-' else call i_to_s_left ( istop, chrtmp ) end if end if call i_to_s_left ( irow, chrtmp1 ) write ( *, '(a)' ) ' ERO: Row ' // trim ( chrtmp1 ) // ' <= ' // & trim ( chrtmp ) // ' Row ' // trim ( chrtmp1 ) return end subroutine row_mul_param ( ierror, base, iform, irow, istop, isbot, line, rval ) ! !******************************************************************************* ! !! ROW_MUL_PARAM gets and checks the row multiply parameters. ! ! ! Modified: ! ! 12 July 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer IERROR, 0 for no error, 1 for an error. ! ! Input, integer IFORM, 0/1/2/3/4, rat/real/dec/int/mod arithmetic. ! ! Output, integer IROW, the row to be multiplied. ! ! Output, integer ISTOP, ISBOT, the multiplier to use for ! fractional or decimal arithmetic. ! ! Workspace, character ( len = 80 ) LINE, used to hold the user's input. ! ! Output, real RVAL, the multiplier to use for real arithmetic. ! implicit none ! integer base integer ierror integer iform integer irow integer isbot integer istop character ( len = 80 ) line character ( len = 80 ) prompt real rval ! prompt = 'row I, multiplier S.' ! ! Read the row number to be multiplied. ! call i_read ( irow, line, prompt, ierror ) if ( ierror /= 0 ) return ! ! Read the multiplier, either RVAL or ISTOP/ISBOT. ! if ( iform == 0 ) then call rat_read ( istop, isbot, line, prompt, ierror ) else if ( iform == 1 ) then call r_read ( rval, line, prompt, ierror ) else if ( iform == 2 ) then call dec_read ( istop, isbot, line, prompt, ierror ) call dec_round ( istop, isbot ) else if ( iform == 3 ) then call i_read ( istop, line, prompt, ierror ) else if ( iform == 4 ) then call i_read ( istop, line, prompt, ierror ) istop = mod ( istop, base ) end if return end subroutine row_op_check ( command, ierror, line2 ) ! !******************************************************************************* ! !! ROW_OP_CHECK checks for commands given in the form of ERO's. ! ! ! Discussion: ! ! The form of the elementary row operation commands includes: ! ! The row interchange command: ! RI1 <=> RI2 ! Note that this will fail if user types "R I1 <=> R I2" ! ! The scalar multiply command: ! RI1 <= S * RI1 ! with or without the "*". ! ! The scalar divide command: ! RI1 <= RI1 / S ! ! The add row command: ! RI1 <= RI1 + S * RI2 ! or ! RI1 <= S * RI2 + RI1 ! ! Modified: ! ! 12 July 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, character ( len = 4 ) COMMAND. ! If the routine decides that the user has input an ERO in the ! natural format, then COMMAND contains the necessary ! one letter MATMAN command to carry out the ERO. ! ! Output, integer IERROR, 0 for no error, 1 for an error. ! ! Workspace, character ( len = 80 ) LINE2, a copy of the user input in LINE. ! implicit none ! character ( len = 22 ) chrtmp character ( len = 10 ) chrtmp1 character ( len = 10 ) chrtmp2 character ( len = 20 ) command integer idbot2 integer idbot3 integer idtop2 integer idtop3 integer ierror integer irow1 integer irow2 integer irow3 integer isbot2 integer isbot3 integer istop2 integer istop3 integer lchar logical ldiv character ( len = 80 ) line2 character ( len = 80 ) string ! command = ' ' ! ! 1. Remove all blanks from the line, and capitalize it. ! call s_blank_delete ( line2 ) call s_cap ( line2 ) ! ! 2. Is the first character an "R" or "ROW"? ! if ( line2(1:1) /= 'R' ) then return end if if ( line2(1:3) == 'ROW' ) then call s_chop ( line2, 1, 3 ) else call s_chop ( line2, 1, 1 ) end if ! ! 3. The next item should be a row number, IROW1. ! call s_to_i ( line2, irow1, ierror, lchar ) if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Your ERO command could not be understood.' write ( *, '(a)' ) ' The first row number "R1" did not make sense.' return end if call s_chop ( line2, 1, lchar ) ! ! 4. Check for the row interchange string "=", "<>", "<=>" or "<->". ! if ( line2(1:2) == '<>' ) then string = '<>' else if ( line2(1:3) == '<=>' ) then string = '<=>' else if ( line2(1:3) == '<->' ) then string = '<->' else if ( line2(1:2) == '<=' ) then string = '<=' else if ( line2(1:2) == '<-' ) then string = '<-' else if ( line2(1:2) == '=>' ) then string = '=>' else if ( line2(1:2) == '->' ) then string = '->' else if ( line2(1:1) == '=' ) then string = '=' else if ( line2(1:2) == ':=' ) then string = ':=' else ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Your ERO command could not be understood.' write ( *, '(a)' ) ' The assignment symbol <=> was missing.' return end if lchar = len_trim ( string ) call s_chop ( line2, 1, lchar ) ! ! 5. The next quantity could be an explicit signed scalar, S2, ! or an implicit +-1. ! if ( line2(1:1) == 'R' ) then istop2 = 1.0E+00 isbot2 = 1.0E+00 else if ( line2(1:2) == '+R' ) then istop2 = 1.0E+00 isbot2 = 1.0E+00 call s_chop ( line2, 1, 1 ) else if ( line2(1:2) == '-R' ) then istop2 = - 1.0E+00 isbot2 = 1.0E+00 call s_chop ( line2, 1, 1 ) else call chrctg ( line2, istop2, isbot2, ierror, lchar ) call s_chop ( line2, 1, lchar ) if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Your ERO command could not be understood.' write ( *, '(a)' ) ' The multiplier S2 did not make sense.' ierror = 1 return end if end if end if ! ! 6. Is the next character an optional "*"? ! if ( line2(1:1) == '*' ) then call s_chop ( line2, 1, 1 ) end if ! ! 7. Is the next character an "R"? ! if ( line2(1:3) == 'ROW' ) then call s_chop ( line2, 1, 3 ) else if ( line2(1:1) == 'R' ) then call s_chop ( line2, 1, 1 ) else ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Your ERO command could not be understood.' write ( *, '(a)' ) ' Could not find the second row index.' return end if ! ! 8. The next item should be a row number, IROW2. ! call s_to_i ( line2, irow2, ierror, lchar ) if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Your ERO command could not be understood.' write ( *, '(a)' ) ' The second row number "R2" did not make sense.' ierror = 1 return else call s_chop ( line2, 1, lchar ) end if ! ! 9. If there's nothing more, this must be an interchange ! or a scaling. Form the equivalent MATMAN command. ! if ( line2 == ' ' ) then if ( irow1 == irow2 ) then command = 'ROW_MUL' if ( isbot2 < 0 ) then isbot2 = - isbot2 istop2 = - istop2 end if call rat_to_s_left ( istop2, isbot2, chrtmp ) call i_to_s_left ( irow1, chrtmp1 ) line2 = chrtmp1 // ' ' // chrtmp call s_blanks_delete ( line2 ) return end if if ( istop2 == 1 .and. isbot2 == 1 ) then command = 'ROW_SWAP' call i_to_s_left ( irow1, chrtmp1 ) call i_to_s_left ( irow2, chrtmp2 ) line2 = chrtmp1 // ' ' // chrtmp2 call s_blanks_delete ( line2 ) return end if ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Your ERO command could not be understood.' write ( *, '(a)' ) ' A MULTIPLY command must have R1 and R2 the same.' write ( *, '(a)' ) ' An INTERCHANGE command cannot have a multiplier.' return end if ! ! 10. Is the next quantity a '/', or perhaps a '*'? ! ldiv = .false. if ( line2(1:1) == '/' ) then ldiv = .true. call s_chop ( line2, 1, 1 ) call chrctg ( line2, idtop2, idbot2, ierror, lchar ) if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Your ERO command could not be understood.' write ( *, '(a)' ) ' The divisor of row 2 did not make sense.' return end if istop2 = istop2 * idbot2 isbot2 = isbot2 * idtop2 if ( irow1 == irow2 ) then if ( ldiv ) then command = 'ROW_DIV' call i_swap ( istop2, isbot2 ) else command = 'ROW_MUL' end if if ( isbot2 < 0 ) then isbot2 = - isbot2 istop2 = - istop2 end if call i_to_s_left ( irow1, chrtmp1 ) call rat_to_s_left ( istop2, isbot2, chrtmp ) line2 = chrtmp1 // ' ' // chrtmp call s_blanks_delete ( line2 ) return end if else if ( line2(1:1) == '*' ) then call s_chop ( line2, 1, 1 ) call chrctg ( line2, idtop2, idbot2, ierror, lchar ) if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Your ERO command could not be understood.' write ( *, '(a)' ) ' The multiplier of row 2 did not make sense.' return end if istop2 = istop2 * idtop2 isbot2 = isbot2 * idbot2 if ( irow1 == irow2 ) then command = 'ROW_MUL' if ( isbot2 < 0 ) then isbot2 = - isbot2 istop2 = - istop2 end if call i_to_s_left ( irow1, chrtmp1 ) call rat_to_s_left ( istop2, isbot2, chrtmp ) line2 = chrtmp1 // ' ' // chrtmp call s_blanks_delete ( line2 ) return end if end if ! ! 11. Is the next quantity a scalar, S3? ! if ( line2(1:2) == '+R' ) then istop3 = 1.0E+00 isbot3 = 1.0E+00 call s_chop ( line2, 1, 1) else if ( line2(1:2) == '-R' ) then istop3 = - 1.0E+00 isbot3 = 1.0E+00 call s_chop ( line2, 1, 1 ) else call chrctg ( line2, istop3, isbot3, ierror, lchar ) if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Your ERO command could not be understood.' write ( *, '(a)' ) ' The multiplier S2 did not make sense.' ierror = 1 return end if call s_chop ( line2, 1, lchar ) end if ! ! 12. Is the next quantity an optional "*"? ! if ( line2(1:1) == '*' ) then call s_chop ( line2, 1, 1 ) end if ! ! 13. Is the next quantity an "R" or ROW? ! if ( line2(1:3) == 'ROW' ) then call s_chop ( line2, 1, 3 ) else if ( line2(1:1) == 'R' ) then call s_chop ( line2, 1, 1) else ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Your ERO command could not be understood.' write ( *, '(a)' ) ' The "R" marking the third row was misplaced.' return end if ! ! 14. The next item should be a row number, IROW3. ! call s_to_i ( line2, irow3, ierror, lchar ) if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Your ERO command could not be understood.' write ( *, '(a)' ) ' The third row number "R3" did not make sense.' ierror = 1 return end if call s_chop ( line2, 1, lchar ) ! ! 15. Is the next quantity a '/', or perhaps a '*'? ! if ( line2(1:1) == '/' ) then call s_chop ( line2, 1, 1 ) call chrctg ( line2, idtop3, idbot3, ierror, lchar ) if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Your ERO command could not be understood.' write ( *, '(a)' ) ' The divisor of R3 did not make sense.' return end if istop3 = istop3 * idbot3 isbot3 = isbot3 * idtop3 else if ( line2(1:1) == '*' ) then call s_chop ( line2, 1, 1 ) call chrctg ( line2, idtop3, idbot3, ierror, lchar ) if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Your ERO command could not be understood.' write ( *, '(a)' ) ' The multiplier of R3 did not make sense.' return end if istop3 = istop3 * idtop3 isbot3 = isbot3 * idbot3 end if ! ! 16. Form the equivalent MATMAN ADD command. ! if ( irow1 == irow2 ) then command = 'ROW_ADD' if ( isbot3 < 0 ) then isbot3 = - isbot3 istop3 = - istop3 end if call i_to_s_left ( irow3, chrtmp1 ) call i_to_s_left ( irow1, chrtmp2 ) call rat_to_s_left ( istop3, isbot3, chrtmp ) line2 = chrtmp // ' ' // chrtmp1 // ' ' // chrtmp2 call s_blanks_delete ( line2 ) else if ( irow1 == irow3 ) then command = 'ROW_ADD' if ( isbot2 < 0 ) then isbot2 = - isbot2 istop2 = - istop2 end if call rat_to_s_left ( istop2, isbot2, chrtmp ) call i_to_s_left ( irow2, chrtmp1 ) call i_to_s_left ( irow1, chrtmp2 ) line2 = chrtmp // ' ' // chrtmp1 // ' ' // chrtmp2 call s_blanks_delete ( line2 ) else ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Your ERO command could not be understood.' write ( *, '(a)' ) ' R2 or R3 must equal R1 in an ERO command.' end if return end subroutine row_swap ( a, a_int, iatop, iabot, ierror, iform, irow1, & irow2, nrow, ncol ) ! !******************************************************************************* ! !! ROW_SWAP swaps two rows of a matrix. ! ! ! Modified: ! ! 12 July 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, real A(NROW,NCOL), the current matrix. ! ! Input/output, integer IATOP(NROW,NCOL), IABOT(NROW,NCOL). ! IATOP and IABOT represent the current rational or decimal matrix. ! ! Output, integer IERROR, 0 for no error, 1 for an error. ! ! Input, integer IFORM, 0/1/2/3/4, rat/real/dec/int/mod arithmetic. ! ! Input, integer IROW1, IROW2, the rows to be swapped. ! ! Input, integer NROW, NCOL, the number of rows and columns in the matrix. ! implicit none ! integer ncol integer nrow ! real a(nrow,ncol) integer a_int(nrow,ncol) integer iabot(nrow,ncol) integer iatop(nrow,ncol) integer ierror integer iform integer irow1 integer irow2 integer j character ( len = 100 ) output ! ! Skip out if the two rows are the same. ! if ( irow1 == irow2 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'ROW_SWAP - Note:' write ( *, '(a)' ) ' You have asked to swap a row with itself!' return end if ! ! Refuse to continue if a row is out of bounds. ! if ( ( irow1 < 1 .or. irow1 > nrow ) .or. & ( irow2 < 1 .or. irow2 > nrow ) ) then ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'ROW_SWAP - Error!' write ( *, '(a)' ) ' One of the rows is illegal!' return end if ! ! Swap the rows. ! do j = 1, ncol if ( iform == 0 ) then call i_swap ( iatop(irow1,j), iatop(irow2,j) ) call i_swap ( iabot(irow1,j), iabot(irow2,j) ) else if ( iform == 1 ) then call r_swap ( a(irow1,j), a(irow2,j) ) else if ( iform == 2 ) then call i_swap ( iatop(irow1,j), iatop(irow2,j) ) call i_swap ( iabot(irow1,j), iabot(irow2,j) ) else if ( iform == 3 ) then call i_swap ( a_int(irow1,j), a_int(irow2,j) ) end if end do write ( output, * ) ' ERO: Row ', irow1, ' <=> Row ', irow2 call s_blanks_delete ( output ) write ( *, '(a)' ) trim ( output ) return end subroutine s_blank_delete ( s ) ! !******************************************************************************* ! !! 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 ( string ) ! !******************************************************************************* ! !! 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 = * ) STRING, the string to be transformed. ! implicit none ! integer i integer j integer nchar character newchr character oldchr character ( len = * ) string character, parameter :: TAB = char ( 9 ) ! nchar = len ( string ) j = 0 newchr = ' ' do i = 1, nchar oldchr = newchr newchr = string(i:i) if ( newchr == TAB ) then newchr = ' ' end if string(i:i) = ' ' if ( oldchr /= ' ' .or. newchr /= ' ' ) then j = j + 1 string(j:j) = newchr end if end do return end subroutine s_cap ( string ) ! !******************************************************************************* ! !! S_CAP replaces any lowercase letters by uppercase ones in a string. ! ! ! Modified: ! ! 16 May 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, character ( len = * ) STRING, the string to be transformed. ! implicit none ! character c integer i integer nchar character ( len = * ) string ! nchar = len ( string ) do i = 1, nchar c = string(i:i) call ch_cap ( c ) string(i:i) = c end do return end subroutine s_chop ( s, ilo, ihi ) ! !******************************************************************************* ! !! S_CHOP "chops out" a portion of a string, and closes up the hole. ! ! ! Example: ! ! S = 'Fred is not a jerk!' ! ! call s_chop ( S, 9, 12 ) ! ! S = 'Fred is a jerk! ' ! ! Modified: ! ! 06 July 1998 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, character ( len = * ) S, 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 = * ) s ! lens = len ( s ) ilo2 = max ( ilo, 1 ) ihi2 = min ( ihi, lens ) if ( ilo2 > ihi2 ) then return end if s(ilo2:lens+ilo2-ihi2-1) = s(ihi2+1:lens) s(lens+ilo2-ihi2:lens) = ' ' return end function s_eqi ( strng1, strng2 ) ! !******************************************************************************* ! !! 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 = * ) STRNG1, STRNG2, the strings to compare. ! ! Output, logical S_EQI, the result of the comparison. ! implicit none ! integer i integer len1 integer len2 integer lenc logical s_eqi character s1 character s2 character ( len = * ) strng1 character ( len = * ) strng2 ! len1 = len ( strng1 ) len2 = len ( strng2 ) lenc = min ( len1, len2 ) s_eqi = .false. do i = 1, lenc s1 = strng1(i:i) s2 = strng2(i:i) call ch_cap ( s1 ) call ch_cap ( s2 ) if ( s1 /= s2 ) then return end if end do do i = lenc + 1, len1 if ( strng1(i:i) /= ' ' ) then return end if end do do i = lenc + 1, len2 if ( strng2(i:i) /= ' ' ) then return end if end do s_eqi = .true. return end subroutine s_read ( string, line, prompt, ierror, iterm ) ! !******************************************************************************* ! !! S_READ extracts a character string from the input buffer. ! ! ! Discussion: ! ! The routine accepts an input LINE and a PROMPT line. ! ! If the LINE is empty, the PROMPT is printed and user input ! ! In either case, enough characters are read from LINE to fill ! STRING and the positions read are removed. ! ! PROMPT is also updated. On satisfactory input of STRING, ! everything in PROMPT up to and including the first comma is removed. ! ! Modified: ! ! 12 July 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, character ( len = * ) STRING. ! The user's response to the PROMPT, as read from LINE. ! ! Input/output, character ( len = 80 ) LINE. ! A buffer containing the user's input. ! ! Input/output, character ( len = 80 ) PROMPT. ! On input, a prompt string that will be printed if necessary. ! ! On output, if STRING has been read, then PROMPT is cleared out ! up to, and including, the first comma. ! ! Output, integer IERROR, 0 for no error, 1 for an error. ! ! Input, integer ITERM, ! 0, No check for terminators. ! 1, Blank, slash, comma, semicolon, equals, greater or ! lesser signs terminate input. ! implicit none ! character c integer i integer ierror integer iterm integer lchar character ( len = 80 ) line integer nchar character null character ( len = 80 ) prompt character ( len = * ) string ! ierror = 0 null = char(0) string = ' ' call chrinp ( ierror, line, prompt ) if ( ierror /= 0 ) then return end if ! ! Null input acceptable for character input only. ! if ( line == ' ' ) return lchar = 0 nchar = len ( string ) do i = 1, nchar if ( lchar /= 0 ) then exit end if c = line(i:i) if ( iterm == 1 ) then if ( c == ' ' .or. c == null .or. c == '/' .or. c == ',' .or. & c == ';' .or. c == '=' ) then lchar = i end if end if if ( lchar == 0 ) then string(i:i) = c end if end do ! ! Chop out the character positions that have been used. ! if ( lchar == 0 ) then lchar = nchar end if call s_chop ( line, 1, lchar ) ! ! Force the string to be flush left by removing leading blanks. ! line = adjustl ( line ) return end subroutine s_to_dec ( s, itop, ibot, length ) ! !******************************************************************************* ! !! S_TO_DEC reads a number from a string, returning a decimal result. ! ! ! Discussion: ! ! The integer may be in real format, for example '2.25'. It ! returns ITOP and IBOT. If the input number is an integer, ITOP ! equals that integer, and IBOT is 1. But in the case of 2.25, ! the program would return ITOP = 225, IBOT = 100. ! ! Legal input is ! ! blanks, ! 2 initial sign, ! blanks, ! 3 whole number, ! 4 decimal point, ! 5 fraction, ! 6 'E' or 'e' or 'D' or 'd', exponent marker, ! 7 exponent sign, ! 8 exponent, ! blanks ! 9 comma or semicolon ! 10end of information ! ! Examples: ! ! S ITOP IBOT Length Meaning ! ! '1' 1 0 1 1 ! ' 1 ' 1 0 6 1 ! '1A' 1 0 1 1 ! '12,34,56' 12 0 3 12 ! ' 34 7' 34 0 4 34 ! '-1E2ABCD' -1 2 4 -100 ! '-1X2ABCD' -1 0 2 -1 ! ' 2E-1' 2 -1 5 0.2 ! '23.45' 2345 -2 5 23.45 ! ! Modified: ! ! 04 February 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) S, the string containing the ! data to be read. Reading begins at position 1 and ! terminate when no more characters ! can be read to form a legal integer. Blanks, commas, ! or other nonnumeric data will, in particular, cause ! the conversion to halt. ! ! Output, integer ITOP, the integer read from the string, ! assuming that no negative exponents or fractional parts ! were used. Otherwise, the 'integer' is ITOP/IBOT. ! ! Output, integer IBOT, the integer divisor required to ! represent numbers which are in real format or have a ! negative exponent. ! ! Output, integer LENGTH, the number of characters used. ! implicit none ! logical ch_is_digit character c integer digit integer exponent integer exponent_sign integer ibot integer ihave integer iterm integer itop integer length integer mantissa_sign character ( len = * ) s logical s_eqi ! itop = 0 ibot = 0 if ( len ( s ) <= 0 ) then length = 0 return end if length = - 1 exponent_sign = 0 mantissa_sign = 1 exponent = 0 ihave = 1 iterm = 0 ! ! Consider the next character in the string. ! do length = length + 1 c = s(length+1:length+1) ! ! Blank. ! if ( c == ' ' ) then if ( ihave == 1 ) then else if ( ihave == 2 ) then else iterm = 1 end if ! ! Comma or semicolon. ! else if ( c == ',' .or. c == ';' ) then if ( ihave /= 1 ) then iterm = 1 ihave = 9 length = length + 1 end if ! ! Minus sign. ! else if ( c == '-' ) then if ( ihave == 1 ) then ihave = 2 mantissa_sign = - 1 else if ( ihave == 6 ) then ihave = 7 exponent_sign = -1 else iterm = 1 end if ! ! Plus sign. ! else if ( c == '+' ) then if ( ihave == 1 ) then ihave = 2 mantissa_sign = +1 else if ( ihave == 6 ) then ihave = 7 exponent_sign = +1 else iterm = 1 end if ! ! Decimal point. ! else if ( c == '.' ) then if ( ihave < 4 ) then ihave = 4 else iterm = 1 end if ! ! Exponent marker. ! else if ( s_eqi ( c, 'E' ) .or. s_eqi ( c, 'D' ) ) then if ( ihave < 6 ) then ihave = 6 else iterm = 1 end if ! ! Digit. ! else if ( ch_is_digit ( c ) ) then if ( ihave <= 3 ) then ihave = 3 call ch_to_digit ( c, digit ) itop = 10 * itop + digit else if ( ihave <= 5 ) then ihave = 5 call ch_to_digit ( c, digit ) itop = 10 * itop + digit ibot = ibot - 1 else if ( ihave <= 8 ) then ihave = 8 call ch_to_digit ( c, digit ) exponent = 10 * exponent + digit else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'S_TO_DEC: Fatal error!' write ( *, '(a,i6)' ) ' IHAVE = ', ihave stop end if ! ! Anything else is regarded as a terminator. ! else iterm = 1 end if if ( iterm == 1 ) then exit end if if ( length + 1 >= len ( s ) ) then length = len ( s ) exit end if end do ! ! Number seems to have terminated. ! Have we got a legal number? ! if ( ihave == 1 ) then return else if ( ihave == 2 .or. ihave == 6 .or. ihave == 7 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'S_TO_DEC - Serious error!' write ( *, '(a)' ) ' Illegal or nonnumeric input:' write ( *, '(a)' ) trim ( s ) return end if ! ! Normalize. ! if ( itop > 0 ) then do while ( mod ( itop, 10 ) == 0 ) itop = itop / 10 ibot = ibot + 1 end do end if ! ! Consolidate the number in the form ITOP * 10**IBOT. ! ibot = ibot + exponent_sign * exponent itop = mantissa_sign * itop if ( itop == 0 ) then ibot = 0 end if return end subroutine s_to_i ( s, ival, ierror, last ) ! !******************************************************************************* ! !! S_TO_I reads an integer value 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_r ( s, r, ierror, lchar ) ! !******************************************************************************* ! !! S_TO_R reads a real number 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 ( ihave > 1 ) 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 ( ihave >= 6 .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. lchar+1 >= nchar ) 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 problem_solve_random ( a, a_int, iatop, iabot, base, iform, & nrow, ncol ) ! !******************************************************************************* ! !! PROBLEM_SOLVE_RANDOM sets up a random linear system problem. ! ! ! Modified: ! ! 28 September 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, real A(NROW,NCOL), the current matrix. ! ! Output, integer IATOP(NROW,NCOL), IABOT(NROW,NCOL). ! IATOP and IABOT represent the current rational or decimal ! matrix. ! ! Input, integer IFORM, 0/1/2/3/4, rat/real/dec/int/mod arithmetic. ! ! Input, integer NROW, NCOL, the number of rows and columns in the matrix. ! implicit none ! integer ncol integer nrow ! real a(nrow,ncol) integer a_int(nrow,ncol) integer base integer i integer iabot(nrow,ncol) integer iatop(nrow,ncol) integer iform integer ihi integer ilo integer ival integer ival2 integer j integer jform integer sol(ncol-1) ! ilo = -10 ihi = 10 do i = 1, nrow do j = 1, ncol-1 call i_random ( ilo, ihi, ival ) a_int(i,j) = ival end do end do do j = 1, ncol-1 call i_random ( ilo, ihi, sol(j) ) end do do i = 1, nrow ival = 0 do j = 1, ncol-1 ival = ival + iatop(i,j) * sol(j) end do iatop(i,ncol) = ival end do jform = iform iform = 3 call form ( a, a_int, iatop, iabot, base, iform, jform, nrow, ncol ) return end subroutine size_set ( ierror, line, maxrow, nrow, title ) ! !******************************************************************************* ! !! SIZE_SET allows the user to specify the size of a new problem. ! ! ! Modified: ! ! 28 September 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer IERROR, 0 for no error, 1 for an error. ! ! Workspace, character ( len = 80 ) LINE, used to hold the user's input. ! ! Input, integer MAXROW, the maximum number of matrix rows. ! ! Output, integer NROW, the number of rows in the matrix. ! ! Input, character ( len = * ) TITLE, contains 'rows' or 'columns'. ! implicit none ! integer maxrow ! integer ierror character ( len = 80 ) line integer nrow character ( len = 80 ) prompt character ( len = * ) title ! ierror = 0 nrow = 0 prompt = 'number of ' // trim ( title ) call i_read ( nrow, line, prompt, ierror ) if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SIZE_SET - Fatal error!' write ( *, '(a)' ) ' I_READ returned error flag.' return end if if ( nrow < 1 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SIZE_SET - Error!' write ( *, '(a)' ) ' Your value is less than 1!' ierror = 1 return end if if ( nrow > maxrow ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SIZE_SET - Error!' write ( *, '(a,i6)' ) ' Your value is greater than ', maxrow ierror = 1 return end if write ( *, '(a,i6)' ) 'Set number of ' // trim ( title ) // ' to ', nrow return end subroutine mat_random ( nrow, ncol, base, iform, a, a_int, iatop, iabot ) ! !******************************************************************************* ! !! MAT_RANDOM sets up a random problem. ! ! ! Modified: ! ! 29 September 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer NROW, NCOL, the number of rows and columns in the matrix. ! ! Input, integer IFORM, 0/1/2/3/4, rat/real/dec/int/mod arithmetic. ! ! Output, real A(NROW,NCOL), the current matrix. ! ! Output, integer IATOP(NROW,NCOL), IABOT(NROW,NCOL). ! IATOP and IABOT represent the current rational or decimal ! matrix. ! implicit none ! integer ncol integer nrow ! real a(nrow,ncol) integer a_int(nrow,ncol) integer base integer i integer iabot(nrow,ncol) integer iatop(nrow,ncol) integer iform integer ihi integer ilo integer ival integer j ! if ( iform /= 4 ) then ilo = -10 ihi = 10 else ilo = 0 ihi = base - 1 end if do i = 1, nrow do j = 1, ncol call i_random ( ilo, ihi, ival ) if ( iform == 0 ) then iatop(i,j) = ival iabot(i,j) = 1 else if ( iform == 1 ) then a(i,j) = real ( ival ) else if ( iform == 2 ) then iatop(i,j) = ival iabot(i,j) = 0 else if ( iform == 3 ) then a_int(i,j) = ival else if ( iform == 4 ) then a_int(i,j) = ival end if end do end do return end subroutine value_read ( irow, icol, iform, base, nrow, ncol, line, a, a_int, & iatop, iabot, ierror ) ! !******************************************************************************* ! !! VALUE_READ allows the user to specify one entry in the array. ! ! ! Modified: ! ! 09 October 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer IROW, ICOL, the row and column of the entry. ! ! Input, integer IFORM, 0/1/2/3/4, rat/real/dec/int/mod arithmetic. ! ! Input, integer BASE, the base of modular arithmetic. ! ! Input, integer NROW, NCOL, the number of rows and columns in the matrix. ! ! Workspace, character ( len = 80 ) LINE, used to hold the user's input. ! ! Input/output, real A(NROW,NCOL), the real matrix. ! ! Input/output, real A_INT(NROW,NCOL), the int/mod matrix. ! ! Input/output, integer IATOP(NROW,NCOL), IABOT(NROW,NCOL), the ! rat/dec matrix. ! ! Output, integer IERROR, 0 for no error, 1 for an error. ! implicit none ! integer ncol integer nrow ! real a(nrow,ncol) integer a_int(nrow,ncol) integer base character ( len = 10 ) chrtmp1 character ( len = 10 ) chrtmp2 integer i_gcd integer iabot(nrow,ncol) integer iatop(nrow,ncol) integer icol integer ierror integer iform integer irow integer isbot integer istop integer itemp integer ival character ( len = 80 ) line character ( len = 80 ) prompt real rval ! call i_to_s_left ( irow, chrtmp1 ) call i_to_s_left ( icol, chrtmp2 ) prompt = 'value of A(' // trim ( chrtmp1 ) // ',' // trim ( chrtmp2 ) // ')' if ( iform == 0 ) then call rat_read ( istop, isbot, line, prompt, ierror ) if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'VALUE_READ - Error!' write ( *, '(a)' ) ' RAT_READ returned error flag!' return end if itemp = i_gcd ( istop, isbot ) iatop(irow,icol) = istop / itemp iabot(irow,icol) = isbot / itemp else if ( iform == 1 ) then call r_read ( rval, line, prompt, ierror ) if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'VALUE_READ - Error!' write ( *, '(a)' ) ' R_READ returned error flag!' return end if a(irow,icol) = rval else if ( iform == 2 ) then call dec_read ( istop, isbot, line, prompt, ierror ) if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'VALUE_READ - Error!' write ( *, '(a)' ) ' DEC_READ returned error flag!' return end if call dec_round ( istop, isbot ) iatop(irow,icol) = istop iabot(irow,icol) = isbot else if ( iform == 3 ) then call i_read ( ival, line, prompt, ierror ) if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'VALUE_READ - Error!' write ( *, '(a)' ) ' I_READ returned error flag!' return end if a_int(irow,icol) = ival else if ( iform == 4 ) then call i_read ( ival, line, prompt, ierror ) if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'VALUE_READ - Error!' write ( *, '(a)' ) ' I_READ returned error flag!' return end if ival = mod ( ival, base ) a_int(irow,icol) = ival end if return end subroutine timestamp ( ) ! !******************************************************************************* ! !! 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