program main !*****************************************************************************80 ! !! MAIN is the main program of FEM2D_POISSON_CG. ! ! Discussion: ! ! This program is a variant of FEM2D_POISSON. That program is ! particularly limited because of its use of banded matrix storage and ! solving routines. ! ! This program discards the banded approach. Instead, it uses a ! sparse matrix storage format and an iterative solver, ! which allow this program to solve larger problems faster. ! ! This program solves the Poisson equation ! ! -DEL H(X,Y) DEL U(X,Y) + K(X,Y) * U(X,Y) = F(X,Y) ! ! in a triangulated region in the plane. ! ! Along the boundary of the region, Dirichlet conditions ! are imposed: ! ! U(X,Y) = G(X,Y) ! ! The code uses continuous piecewise linear basis functions on ! triangles. ! ! Problem specification: ! ! The user defines the geometry by supplying two data files ! which list the node coordinates, and list the nodes that make up ! each element. ! ! The user specifies the right hand side of the Dirichlet boundary ! conditions by supplying a function ! ! subroutine dirichlet_condition ( node_num, node_xy, node_bc ) ! ! The user specifies the coefficient function H(X,Y) of the Poisson ! equation by supplying a routine of the form ! ! subroutine h_coef ( node_num, node_xy, node_h ) ! ! The user specifies the coefficient function K(X,Y) of the Poisson ! equation by supplying a routine of the form ! ! subroutine k_coef ( node_num, node_xy, node_k ) ! ! The user specifies the right hand side of the Poisson equation ! by supplying a routine of the form ! ! subroutine rhs ( node_num, node_xy, node_f ) ! ! Usage: ! ! fem2d_poisson_cg prefix ! ! where ! ! * prefix_nodes.txt is the file containing the coordinates of the nodes; ! ! * prefix_elements.txt is the file containing the indices of nodes ! that make up each element. ! ! Files created include: ! ! * prefix_values.txt, the value of the solution at every node. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 26 January 2013 ! ! Author: ! ! John Burkardt ! ! Local parameters: ! ! Local, real ( kind = 8 ) A(NZ_NUM), the nonzero entries of the coefficient ! matrix. ! ! Local, integer ( kind = 4 ) ELEMENT_NODE(3,ELEMENT_NUM); ! ELEMENT_NODE(I,J) is the global index of local node I in element J. ! ! Local, integer ( kind = 4 ) ELEMENT_NUM, the number of elements. ! ! Local, integer ( kind = 4 ) ELEMENT_ORDER, the element order. ! ! Local, real ( kind = 8 ) F(NODE_NUM), the right hand side. ! ! Local, integer ( kind = 4 ) IA(NZ_NUM), the row indices of the nonzero ! entries of the coefficient matrix. ! ! Local, integer ( kind = 4 ) JA(NZ_NUM), the column indices of the nonzero ! entries of the coefficient matrix. ! ! Local, logical NODE_BOUNDARY(NODE_NUM), is TRUE if the node is ! found to lie on the boundary of the region. ! ! Local, integer ( kind = 4 ) NODE_CONDITION(NODE_NUM), ! indicates the condition used to determine the variable at a node. ! 0, there is no condition (and no variable) at this node. ! 1, a finite element equation is used; ! 2, a Dirichlet condition is used. ! 3, a Neumann condition is used. ! ! Local, integer ( kind = 4 ) NODE_NUM, the number of nodes. ! ! Local, real ( kind = 8 ) NODE_U(NODE_NUM), the finite element coefficients. ! ! Local, real ( kind = 8 ) NODE_XY(2,NODE_NUM), the coordinates of nodes. ! ! Local, integer ( kind = 4 ) NZ_NUM, the number of nonzero entries ! in the coefficient matrix. ! ! Local, integer ( kind = 4 ) QUAD_NUM, the number of quadrature points ! used for assembly. This is currently set to 3, the lowest reasonable ! value. Legal values are 1, 3, 4, 6, 7, 9, 13, and for some problems, ! a value of QUAD_NUM greater than 3 may be appropriate. ! implicit none real ( kind = 8 ), allocatable, dimension (:) :: a integer ( kind = 4 ), allocatable, dimension (:) :: adj_col logical, parameter :: debug = .false. integer ( kind = 4 ), allocatable :: diag(:) integer ( kind = 4 ) dim_num character ( len = 255 ) element_filename integer ( kind = 4 ), allocatable, dimension(:,:) :: element_neighbor integer ( kind = 4 ), allocatable, dimension(:,:) :: element_node integer ( kind = 4 ) element_num integer ( kind = 4 ) element_order real ( kind = 8 ), allocatable, dimension (:) :: f integer ( kind = 4 ), allocatable, dimension (:) :: ia integer ( kind = 4 ) iarg integer ( kind = 4 ) iargc integer ( kind = 4 ) ierr integer ( kind = 4 ) ios integer ( kind = 4 ), allocatable, dimension (:) :: ja integer ( kind = 4 ) k integer ( kind = 4 ) node logical, allocatable, dimension(:) :: node_boundary integer ( kind = 4 ), allocatable, dimension(:) :: node_condition character ( len = 255 ) node_filename integer ( kind = 4 ) node_num real ( kind = 8 ), allocatable, dimension (:) :: node_u real ( kind = 8 ), allocatable, dimension(:,:) :: node_xy integer ( kind = 4 ) num_arg integer ( kind = 4 ) nz_num integer ( kind = 4 ), parameter :: quad_num = 3 character ( len = 255 ) prefix integer ( kind = 4 ) seed character ( len = 255 ) solution_filename call timestamp ( ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FEM2D_POISSON_CG' write ( *, '(a)' ) ' FORTRAN90 version:' write ( *, '(a)' ) ' A version of FEM2D_POISSON using sparse storage' write ( *, '(a)' ) ' and a conjugate gradient solver.' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Solution of the Poisson equation in an arbitrary region' write ( *, '(a)' ) ' in 2 dimensions.' write ( *, '(a)' ) ' ' write ( *, '(a)' ) & ' - DEL H(x,y) DEL U(x,y) + K(x,y) * U(x,y) = F(x,y) in the region' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' U(x,y) = G(x,y) on the boundary.' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' The finite element method is used,' write ( *, '(a)' ) ' with triangular elements,' write ( *, '(a)' ) ' which must be a 3 node linear triangle.' ! ! Get the number of command line arguments. ! num_arg = iargc ( ) if ( num_arg < 1 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'Enter the filename prefix:' read ( *, '(a)', iostat = ios ) prefix if ( ios /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FEM2D_POISSON_CG - Fatal error!' write ( *, '(a)' ) ' Unexpected read error!' stop end if else iarg = 1 call getarg ( iarg, prefix ) end if ! ! Create the filenames. ! node_filename = trim ( prefix ) // '_nodes.txt' element_filename = trim ( prefix ) // '_elements.txt' solution_filename = trim ( prefix ) // '_values.txt' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Node file is "' // trim ( node_filename ) // '".' write ( *, '(a)' ) ' Element file is "' // trim ( element_filename ) & // '".' ! ! Read the node coordinate file. ! call r8mat_header_read ( node_filename, dim_num, node_num ) write ( *, '(a,i8)' ) ' Number of nodes = ', node_num allocate ( node_boundary(node_num) ) allocate ( node_condition(node_num) ) allocate ( node_xy(dim_num,node_num) ) call r8mat_data_read ( node_filename, dim_num, node_num, node_xy ) call r8mat_transpose_print_some ( dim_num, node_num, node_xy, 1, 1, & dim_num, 10, ' First 10 nodes' ) ! ! Read the element description file. ! call i4mat_header_read ( element_filename, element_order, element_num ) write ( *, '(a)' ) ' ' write ( *, '(a,i8)' ) ' Element order = ', element_order write ( *, '(a,i8)' ) ' Number of elements = ', element_num if ( element_order /= 3 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FEM2D_POISSON_CG - Fatal error!' write ( *, '(a,i8)' ) ' The input triangulation has order ', element_order write ( *, '(a)' ) ' However, a triangulation of order 3 is required.' stop end if allocate ( element_node(3,element_num) ) call i4mat_data_read ( element_filename, element_order, element_num, & element_node ) call i4mat_transpose_print_some ( 3, element_num, & element_node, 1, 1, 3, 10, ' First 10 elements' ) write ( *, '(a)' ) ' ' write ( *, '(a,i8)' ) ' Quadrature order = ', quad_num ! ! Determine which nodes are boundary nodes and which have a ! finite element unknown. Then set the boundary values. ! call triangulation_order3_boundary_node ( node_num, element_num, & element_node, node_boundary ) ! ! Determine the node conditions. ! For now, we'll just assume all boundary nodes are Dirichlet. ! node_condition(1:node_num) = 1 do node = 1, node_num if ( node_boundary(node) ) then node_condition(node) = 2 end if end do ! ! Determine the element neighbor array, just so we can estimate ! the nonzeros. ! allocate ( element_neighbor(3,element_num) ) call triangulation_order3_neighbor_triangles ( element_num, element_node, & element_neighbor ) ! ! Count the number of nonzeros. ! allocate ( adj_col(1:node_num+1) ) call triangulation_order3_adj_count ( node_num, element_num, element_node, & element_neighbor, nz_num, adj_col ) write ( *, '(a)' ) ' ' write ( *, '(a,i8)' ) ' Number of nonzero coefficients NZ_NUM = ', nz_num ! ! Set up the sparse row and column index vectors. ! allocate ( ia(1:nz_num) ) allocate ( ja(1:nz_num) ) call triangulation_order3_adj_set2 ( node_num, element_num, element_node, & element_neighbor, nz_num, adj_col, ia, ja ) deallocate ( adj_col ) deallocate ( element_neighbor ) ! ! Index the diagonal elements for use by the CG solver. ! allocate ( diag(1:node_num) ) call diag_index ( nz_num, ia, ja, node_num, diag ) if ( debug ) then call i4vec_print ( node_num, diag, ' Diagonal adjacency vector:' ) end if ! ! Allocate space for the coefficient matrix A and right hand side F. ! allocate ( a(nz_num) ) allocate ( f(node_num) ) allocate ( node_u(node_num) ) ! ! Assemble the finite element coefficient matrix A and the right-hand side F. ! call assemble_poisson_dsp ( node_num, node_xy, element_num, & element_node, quad_num, nz_num, ia, ja, a, f ) if ( debug ) then call dsp_print_some ( node_num, node_num, nz_num, ia, ja, a, 1, 1, & 10, 10, ' Part of Finite Element matrix A:' ) call r8vec_print_some ( node_num, f, 1, 10, & ' Part of right hand side vector F:' ) end if ! ! Adjust the linear system to account for Dirichlet boundary conditions. ! call dirichlet_apply_dsp ( node_num, node_xy, node_condition, nz_num, & ia, ja, a, f ) if ( debug ) then call dsp_print_some ( node_num, node_num, nz_num, ia, ja, a, 1, 1, & 10, 10, ' Part of A after adjustment for Dirichlet condition:' ) call r8vec_print_some ( node_num, f, 1, 10, & ' Part of F after adjustment for Dirichlet condition:' ) end if ! ! Solve the linear system using the conjugate gradient method. ! do k = 1, nz_num if ( ia(k) < 1 .or. node_num < ia(k) ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FEM2D_POISSON_CG - Fatal error!' write ( *, '(a)' ) ' Illegal IA(K)' stop end if if ( ja(k) < 1 .or. node_num < ja(k) ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FEM2D_POISSON_CG - Fatal error!' write ( *, '(a)' ) ' Illegal JA(K)' stop end if end do call solve_cg ( node_num, diag, nz_num, ia, ja, a, f, node_u ) if ( debug .or. .true. ) then call r8vec_print_some ( node_num, node_u, 1, 10, & ' Part of the solution vector U:' ) end if ! ! Write an ASCII file that can be read into MATLAB. ! call r8mat_write ( solution_filename, 1, node_num, node_u ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FEM2D_POISSON_CG:' write ( *, '(a)' ) ' Wrote an ASCII file' write ( *, '(a)' ) ' "' // trim ( solution_filename ) // '"' write ( *, '(a)' ) ' of the form' write ( *, '(a)' ) ' U ( X(I), Y(I) )' write ( *, '(a)' ) ' which can be used for plotting.' ! ! Free memory. ! deallocate ( a ) deallocate ( diag ) deallocate ( f ) deallocate ( element_node ) deallocate ( ia ) deallocate ( ja ) deallocate ( node_boundary ) deallocate ( node_condition ) deallocate ( node_u ) deallocate ( node_xy ) ! ! Terminate. ! write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FEM2D_POISSON_CG:' write ( *, '(a)' ) ' Normal end of execution.' write ( *, '(a)' ) ' ' call timestamp ( ) stop end subroutine assemble_poisson_dsp ( node_num, node_xy, element_num, & element_node, quad_num, nz_num, ia, ja, a, f ) !*****************************************************************************80 ! !! ASSEMBLE_POISSON_DSP assembles the system for the Poisson equation. ! ! Discussion: ! ! The matrix is sparse, and stored in the DSP or "sparse triple" format. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 14 November 2006 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) NODE_NUM, the number of nodes. ! ! Input, real ( kind = 8 ) NODE_XY(2,NODE_NUM), the ! coordinates of nodes. ! ! Input, integer ( kind = 4 ) ELEMENT_NUM, the number of elements. ! ! Input, integer ( kind = 4 ) ELEMENT_NODE(3,ELEMENT_NUM); ! ELEMENT_NODE(I,J) is the global index of local node I in element J. ! ! Input, integer ( kind = 4 ) QUAD_NUM, the number of quadrature points ! used in assembly. ! ! Input, integer ( kind = 4 ) NZ_NUM, the number of nonzero entries. ! ! Input, integer ( kind = 4 ) IA(NZ_NUM), JA(NZ_NUM), the row and column ! indices of the nonzero entries. ! ! Output, real ( kind = 8 ) A(NZ_NUM), the nonzero entries of the matrix. ! ! Output, real ( kind = 8 ) F(NODE_NUM), the right hand side. ! ! Local parameters: ! ! Local, real ( kind = 8 ) BI, DBIDX, DBIDY, the value of some basis function ! and its first derivatives at a quadrature point. ! ! Local, real ( kind = 8 ) BJ, DBJDX, DBJDY, the value of another basis ! function and its first derivatives at a quadrature point. ! implicit none integer ( kind = 4 ) node_num integer ( kind = 4 ) quad_num integer ( kind = 4 ) element_num integer ( kind = 4 ) nz_num real ( kind = 8 ), dimension(nz_num) :: a real ( kind = 8 ) area integer ( kind = 4 ) basis real ( kind = 8 ) bi real ( kind = 8 ) bj real ( kind = 8 ) dbidx real ( kind = 8 ) dbidy real ( kind = 8 ) dbjdx real ( kind = 8 ) dbjdy integer ( kind = 4 ) element integer ( kind = 4 ), dimension(3,element_num) :: element_node real ( kind = 8 ), dimension(node_num) :: f integer ( kind = 4 ) i integer ( kind = 4 ) ia(nz_num) integer ( kind = 4 ) j integer ( kind = 4 ) ja(nz_num) integer ( kind = 4 ) k real ( kind = 8 ), dimension(2,node_num) :: node_xy real ( kind = 8 ), allocatable, dimension(:) :: phys_h real ( kind = 8 ), allocatable, dimension(:) :: phys_k real ( kind = 8 ), allocatable, dimension(:) :: phys_rhs real ( kind = 8 ) phys_xy(2,quad_num) integer ( kind = 4 ) quad real ( kind = 8 ), allocatable, dimension(:) :: quad_w real ( kind = 8 ), allocatable, dimension(:,:) :: quad_xy real ( kind = 8 ), dimension (2,3) :: t3 integer ( kind = 4 ) test real ( kind = 8 ) triangle_area_2d real ( kind = 8 ) w(quad_num) ! ! Initialize the arrays to zero. ! f(1:node_num) = 0.0D+00 a(1:nz_num) = 0.0D+00 allocate ( phys_h(1:quad_num) ) allocate ( phys_k(1:quad_num) ) allocate ( phys_rhs(1:quad_num) ) allocate ( quad_w(1:quad_num) ) allocate ( quad_xy(1:2,1:quad_num) ) ! ! Get the quadrature weights and nodes. ! call quad_rule ( quad_num, quad_w, quad_xy ) ! ! Add up all quantities associated with the ELEMENT-th element. ! do element = 1, element_num ! ! Make a copy of the triangle. ! t3(1:2,1:3) = node_xy(1:2,element_node(1:3,element)) ! ! Map the quadrature points QUAD_XY to points XY in the physical triangle. ! call reference_to_physical_t3 ( t3, quad_num, quad_xy, phys_xy ) area = abs ( triangle_area_2d ( t3 ) ) w(1:quad_num) = area * quad_w(1:quad_num) call rhs ( quad_num, phys_xy, phys_rhs ) call h_coef ( quad_num, phys_xy, phys_h ) call k_coef ( quad_num, phys_xy, phys_k ) ! ! Consider the QUAD-th quadrature point. ! do quad = 1, quad_num ! ! Consider the TEST-th test function. ! ! We generate an integral for every node associated with an unknown. ! But if a node is associated with a boundary condition, we do nothing. ! do test = 1, 3 i = element_node(test,element) call basis_one_t3 ( t3, test, phys_xy(1:2,quad), bi, dbidx, dbidy ) f(i) = f(i) + w(quad) * phys_rhs(quad) * bi ! ! Consider the BASIS-th basis function, which is used to form the ! value of the solution function. ! do basis = 1, 3 j = element_node(basis,element) call basis_one_t3 ( t3, basis, phys_xy(1:2,quad), bj, dbjdx, dbjdy ) call dsp_ij_to_k ( nz_num, ia, ja, i, j, k ) a(k) = a(k) + w(quad) * ( & phys_h(quad) * ( dbidx * dbjdx + dbidy * dbjdy ) & + phys_k(quad) * bj * bi ) end do end do end do end do deallocate ( phys_h ) deallocate ( phys_k ) deallocate ( phys_rhs ) deallocate ( quad_w ) deallocate ( quad_xy ) return end subroutine basis_one_t3 ( t, i, p, qi, dqidx, dqidy ) !*****************************************************************************80 ! !! BASIS_ONE_T3 evaluates a linear basis function. ! ! Discussion: ! ! The routine is given the coordinates of the nodes of a triangle. ! ! 3 ! / \ ! / \ ! / \ ! 1-------2 ! ! It evaluates the linear basis function Q(I)(X,Y) associated with ! node I, which has the property that it is a linear function ! which is 1 at node I and zero at the other two nodes. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 07 January 2006 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real T(2,3), the coordinates of the nodes. ! ! Input, integer ( kind = 4 ) I, the index of the desired basis function. ! I should be between 1 and 3. ! ! Input, real P(2), the coordinates of a point at which the basis ! function is to be evaluated. ! ! Output, real QI, DQIDX, DQIDY, the values of the basis function ! and its X and Y derivatives. ! implicit none real ( kind = 8 ) area real ( kind = 8 ) dqidx real ( kind = 8 ) dqidy integer ( kind = 4 ) i integer ( kind = 4 ) i4_wrap integer ( kind = 4 ) ip1 integer ( kind = 4 ) ip2 real ( kind = 8 ) p(2) real ( kind = 8 ) qi real ( kind = 8 ) t(2,3) area = t(1,1) * ( t(2,2) - t(2,3) ) & + t(1,2) * ( t(2,3) - t(2,1) ) & + t(1,3) * ( t(2,1) - t(2,2) ) if ( area == 0.0D+00 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'BASIS_ONE_T3 - Fatal error!' write ( *, '(a)' ) ' Element has zero area.' stop end if if ( i < 1 .or. 3 < i ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'BASIS_ONE_T3 - Fatal error!' write ( *, '(a)' ) ' Basis index I is not between 1 and 3.' write ( *, '(a,i8)' ) ' I = ', i stop end if ip1 = i4_wrap ( i + 1, 1, 3 ) ip2 = i4_wrap ( i + 2, 1, 3 ) qi = ( ( t(1,ip2) - t(1,ip1) ) * ( p(2) - t(2,ip1) ) & - ( t(2,ip2) - t(2,ip1) ) * ( p(1) - t(1,ip1) ) ) / area dqidx = - ( t(2,ip2) - t(2,ip1) ) / area dqidy = ( t(1,ip2) - t(1,ip1) ) / area return end subroutine cg_rc ( n, b, x, r, z, p, q, job ) !*****************************************************************************80 ! !! CG_RC is a reverse communication conjugate gradient routine. ! ! Discussion: ! ! This routine seeks a solution of the linear system A*x=b ! where b is a given right hand side vector, A is an n by n ! symmetric positive definite matrix, and x is an unknown vector ! to be determined. ! ! Under the assumptions that the matrix A is large and sparse, ! the conjugate gradient method may provide a solution when ! a direct approach would be impractical because of excessive ! requirements of storage or even of time. ! ! The conjugate gradient method presented here does not require the ! user to store the matrix A in a particular way. Instead, it only ! supposes that the user has a way of calculating ! y = alpha * A * x + b * y ! and of solving the preconditioned linear system ! M * x = b ! where M is some preconditioning matrix, which might be merely ! the identity matrix, or a diagonal matrix containing the ! diagonal entries of A. ! ! This routine was extracted from the "templates" package. ! There, it was not intended for direct access by a user; ! instead, a higher routine called "cg()" was called once by ! the user. The cg() routine then made repeated calls to ! cgrevcom() before returning the result to the user. ! ! The reverse communication feature of cgrevcom() makes it, by itself, ! a very powerful function. It allows the user to handle issues of ! storage and implementation that would otherwise have to be ! mediated in a fixed way by the function argument list. Therefore, ! this version of cgrecom() has been extracted from the templates ! library and documented as a stand-alone procedure. ! ! The user sets the value of JOB to 1 before the first call, ! indicating the beginning of the computation, and to the value of ! 2 thereafter, indicating a continuation call. ! The output value of JOB is set by cgrevcom(), which ! will return with an output value of JOB that requests a particular ! new action from the user. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 12 January 2013 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Richard Barrett, Michael Berry, Tony Chan, James Demmel, ! June Donato, Jack Dongarra, Victor Eijkhout, Roidan Pozo, ! Charles Romine, Henk van der Vorst, ! Templates for the Solution of Linear Systems: ! Building Blocks for Iterative Methods, ! SIAM, 1994, ! ISBN: 0898714710, ! LC: QA297.8.T45. ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the dimension of the matrix. ! ! Input, real ( kind = 8 ) B(N), the right hand side vector. ! ! Input/output, real ( kind = 8 ) X(N). On first call, the user ! should store an initial guess for the solution in X. On return with ! JOB = 4, X contains the latest solution estimate. ! ! Input/output, real ( kind = 8 ) R(N), Z(N), P(N), Q(N), ! information used by the program during the calculation. The user ! does not need to initialize these vectors. However, specific ! return values of JOB may require the user to carry out some computation ! using data in some of these vectors. ! ! Input/output, integer ( kind = 4 ) JOB, communicates the task to be done. ! The user needs to set the input value of JOB to 1, before the first call, ! and then to 2 for every subsequent call for the given problem. ! The output value of JOB indicates the requested user action. ! * JOB = 1, compute Q = A * P; ! * JOB = 2: solve M*Z=R, where M is the preconditioning matrix; ! * JOB = 3: compute R = R - A * X; ! * JOB = 4: check the residual R for convergence. ! If satisfactory, terminate the iteration. ! If too many iterations were taken, terminate the iteration. ! implicit none integer ( kind = 4 ) n real ( kind = 8 ) alpha real ( kind = 8 ) b(n) real ( kind = 8 ) beta integer ( kind = 4 ) iter integer ( kind = 4 ) job real ( kind = 8 ) p(n) real ( kind = 8 ) pdotq real ( kind = 8 ) q(n) real ( kind = 8 ) r(n) real ( kind = 8 ) rho real ( kind = 8 ) rho_old integer ( kind = 4 ) rlbl real ( kind = 8 ) x(n) real ( kind = 8 ) z(n) ! ! Some local variables must be preserved between calls. ! save iter save rho save rho_old save rlbl ! ! Initialization. ! Ask the user to compute the initial residual. ! if ( job == 1 ) then r(1:n) = b(1:n) job = 3 rlbl = 2 ! ! Begin first conjugate gradient loop. ! Ask the user for a preconditioner solve. ! else if ( rlbl == 2 ) then iter = 1 job = 2 rlbl = 3 ! ! Compute the direction. ! Ask the user to compute ALPHA. ! Save A*P to Q. ! else if ( rlbl == 3 ) then rho = dot_product ( r, z ) if ( 1 < iter ) then beta = rho / rho_old z(1:n) = z(1:n) + beta * p(1:n) end if p(1:n) = z(1:n) job = 1 rlbl = 4 ! ! Compute current solution vector. ! Ask the user to check the stopping criterion. ! else if ( rlbl == 4 ) then pdotq = dot_product ( p, q ) alpha = rho / pdotq x(1:n) = x(1:n) + alpha * p(1:n) r(1:n) = r(1:n) - alpha * q(1:n) job = 4 rlbl = 5 ! ! Begin the next step. ! Ask for a preconditioner solve. ! else if ( rlbl == 5 ) then rho_old = rho iter = iter + 1 job = 2 rlbl = 3 end if return end subroutine ch_cap ( c ) !*****************************************************************************80 ! !! CH_CAP capitalizes a single character. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 19 July 1998 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, character C, the character to capitalize. ! implicit none character c integer ( kind = 4 ) itemp itemp = ichar ( c ) if ( 97 <= itemp .and. itemp <= 122 ) then c = char ( itemp - 32 ) end if return end function ch_eqi ( c1, c2 ) !*****************************************************************************80 ! !! CH_EQI is a case insensitive comparison of two characters for equality. ! ! Example: ! ! CH_EQI ( 'A', 'a' ) is .TRUE. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 28 July 2000 ! ! 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 c1_cap character c2 character c2_cap c1_cap = c1 c2_cap = c2 call ch_cap ( c1_cap ) call ch_cap ( c2_cap ) if ( c1_cap == c2_cap ) then ch_eqi = .true. else ch_eqi = .false. end if return end subroutine ch_to_digit ( c, digit ) !*****************************************************************************80 ! !! CH_TO_DIGIT returns the value of a base 10 digit. ! ! Example: ! ! C DIGIT ! --- ----- ! '0' 0 ! '1' 1 ! ... ... ! '9' 9 ! ' ' 0 ! 'X' -1 ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 04 August 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character C, the decimal digit, '0' through '9' or blank ! are legal. ! ! Output, integer ( kind = 4 ) DIGIT, the corresponding value. ! If C was 'illegal', then DIGIT is -1. ! implicit none character c integer ( kind = 4 ) 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 diag_index ( m, ia, ja, n, diag ) !*****************************************************************************80 ! !! DIAG_INDEX determines where the diagonal matrix entries are stored. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 26 January 2013 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) M, the number of adjacencies. ! ! Input, integer ( kind = 4 ) IA(M), JA(M), the row and column indices ! of adjacencies. ! ! Input, integer ( kind = 4 ) N, the number of nodes. ! ! Output, integer ( kind = 4 ) DIAG(N), contains for each index 1 <= I <= N, ! the unique index J such that IA[J] = JA[J] = I. ! implicit none integer ( kind = 4 ) m integer ( kind = 4 ) n integer ( kind = 4 ) diag(n) integer ( kind = 4 ) i integer ( kind = 4 ) ia(m) integer ( kind = 4 ) j integer ( kind = 4 ) ja(m) diag(1:n) = -1 do j = 1, m if ( ia(j) == ja(j) ) then i = ia(j) if ( diag(i) /= -1 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'DIAG_INDEX - Fatal error!' write ( *, '(a)' ) ' Multiple occurrences of diagonal pairs.' write ( *, '(a,i4,a,i4,a,i4,a)' ) & ' IA(', j, ') = JA(', j, ') = ', ia(j), ' and' write ( *, '(a,i4,a,i4,a,i4)' ) & ' IA(', diag(i), ') = JA(', diag(i), ') = ', ia(j) stop end if diag(i) = j end if end do do i = 1, n if ( diag(i) == -1 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'DIAG_INDEX - Fatal error!' write ( *, '(a,i4,a)' ) ' DIAG(', i, ') = -1.' stop end if end do return end subroutine dirichlet_apply_dsp ( node_num, node_xy, node_condition, & nz_num, ia, ja, a, f ) !*****************************************************************************80 ! !! DIRICHLET_APPLY_DSP accounts for Dirichlet boundary conditions. ! ! Discussion: ! ! It is assumed that the matrix A and right hand side F have already been ! set up as though there were no boundary conditions. This routine ! then modifies A and F, essentially replacing the finite element equation ! at a boundary node NODE by a trivial equation of the form ! ! A(NODE,NODE) * U(NODE) = NODE_BC(NODE) ! ! where A(NODE,NODE) = 1. ! ! This routine assumes that the coefficient matrix is stored in a ! sparse triplet format. ! ! This routine implicitly assumes that the sparse matrix has a storage ! location for every diagonal element...or at least for those diagonal ! elements corresponding to boundary nodes. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 12 July 2007 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) NODE_NUM, the number of nodes. ! ! Input, real ( kind = 8 ) NODE_XY(2,NODE_NUM), the coordinates of nodes. ! ! Input, integer ( kind = 4 ) NODE_CONDITION(NODE_NUM), reports the ! condition used to set the unknown associated with the node. ! 0, unknown. ! 1, finite element equation. ! 2, Dirichlet condition; ! 3, Neumann condition. ! ! Input, integer ( kind = 4 ) NZ_NUM, the number of nonzero entries. ! ! Input, integer ( kind = 4 ) IA(NZ_NUM), JA(NZ_NUM), the row and column ! indices of the nonzero entries. ! ! Input/output, real ( kind = 8 ) A(NZ_NUM), the coefficient matrix, ! stored in sparse triplet format; on output, the matrix has been adjusted ! for Dirichlet boundary conditions. ! ! Input/output, real ( kind = 8 ) F(NODE_NUM), the right hand side. ! On output, the right hand side has been adjusted for Dirichlet ! boundary conditions. ! implicit none integer ( kind = 4 ) node_num integer ( kind = 4 ) nz_num real ( kind = 8 ), dimension(nz_num) :: a integer ( kind = 4 ) column integer ( kind = 4 ), parameter :: DIRICHLET = 2 real ( kind = 8 ), dimension(node_num) :: f integer ( kind = 4 ) ia(nz_num) integer ( kind = 4 ) ja(nz_num) integer ( kind = 4 ) node real ( kind = 8 ), dimension ( node_num ) :: node_bc integer ( kind = 4 ) node_condition(node_num) real ( kind = 8 ), dimension(2,node_num) :: node_xy integer ( kind = 4 ) nz ! ! Retrieve the Dirichlet boundary condition value at every node. ! call dirichlet_condition ( node_num, node_xy, node_bc ) ! ! Consider every matrix entry, NZ. ! ! If the row I corresponds to a boundary node, then ! zero out all off diagonal matrix entries, set the diagonal to 1, ! and the right hand side to the Dirichlet boundary condition value. ! do nz = 1, nz_num node = ia(nz) if ( node_condition(node) == DIRICHLET ) then column = ja(nz) if ( column == node ) then a(nz) = 1.0D+00 f(node) = node_bc(node) else a(nz) = 0.0D+00 end if end if end do return end subroutine dsp_ij_to_k ( nz_num, row, col, i, j, k ) !*****************************************************************************80 ! !! DSP_IJ_TO_K seeks the compressed index of the (I,J) entry of A. ! ! Discussion: ! ! If A(I,J) is nonzero, then its value is stored in location K. ! ! This routine searches the DSP storage structure for the index K ! corresponding to (I,J), returning -1 if no such entry was found. ! ! This routine assumes that the data structure has been sorted, ! so that the entries of ROW are ascending sorted, and that the ! entries of COL are ascending sorted, within the group of entries ! that have a common value of ROW. ! ! The DSP storage format stores the row, column and value of each nonzero ! entry of a sparse matrix. ! ! The DSP format is used by CSPARSE ("sparse triplet"), DLAP/SLAP ! ("nonsymmetric SLAP triad"), by MATLAB, and by SPARSEKIT ("COO" format). ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 11 July 2007 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) NZ_NUM, the number of nonzero elements in ! the matrix. ! ! Input, integer ( kind = 4 ) ROW(NZ_NUM), COL(NZ_NUM), the row and ! column indices of the nonzero elements. ! ! Input, integer ( kind = 4 ) I, J, the row and column indices of the ! matrix entry. ! ! Output, integer ( kind = 4 ) K, the DSP index of the (I,J) entry. ! implicit none integer ( kind = 4 ) nz_num integer ( kind = 4 ) col(nz_num) integer ( kind = 4 ) hi integer ( kind = 4 ) i integer ( kind = 4 ) j integer ( kind = 4 ) k integer ( kind = 4 ) lo integer ( kind = 4 ) md integer ( kind = 4 ) row(nz_num) lo = 1 hi = nz_num do if ( hi < lo ) then k = -1 exit end if md = ( lo + hi ) / 2 if ( row(md) < i .or. ( row(md) == i .and. col(md) < j ) ) then lo = md + 1 else if ( i < row(md) .or. ( row(md) == i .and. j < col(md) ) ) then hi = md - 1 else k = md exit end if end do return end subroutine dsp_print_some ( m, n, nz_num, row, col, a, ilo, jlo, & ihi, jhi, title ) !*****************************************************************************80 ! !! DSP_PRINT_SOME prints some of a DSP matrix. ! ! Discussion: ! ! This version of DSP_PRINT_SOME has been specifically modified to allow, ! and correctly handle, the case in which a single matrix location ! A(I,J) is referenced more than once by the sparse matrix structure. ! In such cases, the routine prints out the sum of all the values. ! ! The DSP storage format stores the row, column and value of each nonzero ! entry of a sparse matrix. ! ! It is possible that a pair of indices (I,J) may occur more than ! once. Presumably, in this case, the intent is that the actual value ! of A(I,J) is the sum of all such entries. This is not a good thing ! to do, but I seem to have come across this in MATLAB. ! ! The DSP format is used by CSPARSE ("sparse triplet"), DLAP/SLAP ! ("nonsymmetric SLAP triad"), by MATLAB, and by SPARSEKIT ("COO" format). ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 06 September 2006 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) M, N, the number of rows and columns ! of the matrix. ! ! Input, integer ( kind = 4 ) NZ_NUM, the number of nonzero elements ! in the matrix. ! ! Input, integer ( kind = 4 ) ROW(NZ_NUM), COL(NZ_NUM), the row and column ! indices of the nonzero elements. ! ! Input, real ( kind = 8 ) A(NZ_NUM), the nonzero elements of the matrix. ! ! Input, integer ( kind = 4 ) ILO, JLO, IHI, JHI, the first row and ! column, and the last row and column to be printed. ! ! Input, character ( len = * ) TITLE, a title to print. ! implicit none integer ( kind = 4 ), parameter :: incx = 5 integer ( kind = 4 ) nz_num real ( kind = 8 ) a(nz_num) real ( kind = 8 ) aij(incx) integer ( kind = 4 ) col(nz_num) integer ( kind = 4 ) i integer ( kind = 4 ) i2hi integer ( kind = 4 ) i2lo integer ( kind = 4 ) ihi integer ( kind = 4 ) ilo integer ( kind = 4 ) inc integer ( kind = 4 ) j integer ( kind = 4 ) j2 integer ( kind = 4 ) j2hi integer ( kind = 4 ) j2lo integer ( kind = 4 ) jhi integer ( kind = 4 ) jlo integer ( kind = 4 ) k integer ( kind = 4 ) m integer ( kind = 4 ) n integer ( kind = 4 ) row(nz_num) character ( len = * ) title if ( 0 < len_trim ( title ) ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) trim ( title ) end if ! ! Print the columns of the matrix, in strips of 5. ! do j2lo = jlo, jhi, incx j2hi = j2lo + incx - 1 j2hi = min ( j2hi, n ) j2hi = min ( j2hi, jhi ) inc = j2hi + 1 - j2lo write ( *, '(a)' ) ' ' write ( *, '('' Col: '',5(i7,7x))' ) ( j, j = j2lo, j2hi ) write ( *, '(a)' ) ' Row' write ( *, '(a)' ) ' ---' ! ! Determine the range of the rows in this strip. ! i2lo = max ( ilo, 1 ) i2hi = min ( ihi, m ) do i = i2lo, i2hi ! ! Print out (up to) 5 entries in row I, that lie in the current strip. ! aij(1:inc) = 0.0D+00 ! ! Is matrix entry K actually the value of A(I,J), with J2LO <= J <= J2HI? ! Because MATLAB seems to allow for multiple (I,J,A) entries, we have ! to sum up what we find. ! do k = 1, nz_num if ( i == row(k) .and. & j2lo <= col(k) .and. & col(k) <= j2hi ) then j2 = col(k) - j2lo + 1 aij(j2) = aij(j2) + a(k) end if end do if ( any ( aij(1:inc) /= 0.0D+00 ) ) then write ( *, '(i5,1x,5g14.6)' ) i, aij(1:inc) end if end do end do return end subroutine file_column_count ( input_file_name, column_num ) !*****************************************************************************80 ! !! FILE_COLUMN_COUNT counts the number of columns in the first line of a file. ! ! Discussion: ! ! The file is assumed to be a simple text file. ! ! Most lines of the file is presumed to consist of COLUMN_NUM words, ! separated by spaces. There may also be some blank lines, and some ! comment lines, ! which have a "#" in column 1. ! ! The routine tries to find the first non-comment non-blank line and ! counts the number of words in that line. ! ! If all lines are blanks or comments, it goes back and tries to analyze ! a comment line. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 21 June 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) INPUT_FILE_NAME, the name of the file. ! ! Output, integer ( kind = 4 ) COLUMN_NUM, the number of columns in the file. ! implicit none integer ( kind = 4 ) column_num logical got_one character ( len = * ) input_file_name integer ( kind = 4 ) input_status integer ( kind = 4 ) input_unit character ( len = 255 ) line ! ! Open the file. ! call get_unit ( input_unit ) open ( unit = input_unit, file = input_file_name, status = 'old', & form = 'formatted', access = 'sequential', iostat = input_status ) if ( input_status /= 0 ) then column_num = -1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FILE_COLUMN_COUNT - Fatal error!' write ( *, '(a,i8)' ) ' Could not open the input file "' & // trim ( input_file_name ) // '" on unit ', input_unit return end if ! ! Read one line, but skip blank lines and comment lines. ! got_one = .false. do read ( input_unit, '(a)', iostat = input_status ) line if ( input_status /= 0 ) then exit end if if ( len_trim ( line ) == 0 ) then cycle end if if ( line(1:1) == '#' ) then cycle end if got_one = .true. exit end do if ( .not. got_one ) then rewind ( input_unit ) do read ( input_unit, '(a)', iostat = input_status ) line if ( input_status /= 0 ) then exit end if if ( len_trim ( line ) == 0 ) then cycle end if got_one = .true. exit end do end if close ( unit = input_unit ) if ( .not. got_one ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FILE_COLUMN_COUNT - Warning!' write ( *, '(a)' ) ' The file does not seem to contain any data.' column_num = -1 return end if call s_word_count ( line, column_num ) return end subroutine file_name_specification ( node_file_name, element_file_name ) !*****************************************************************************80 ! !! FILE_NAME_SPECIFICATION determines the names of the input files. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 04 October 2006 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, character ( len = * ) NODE_FILE_NAME, the name of the node file. ! ! Output, character ( len = * ) ELEMENT_FILE_NAME, the name ! of the element file. ! implicit none integer ( kind = 4 ) arg_num character ( len = * ) :: element_file_name integer ( kind = 4 ) iarg integer ( kind = 4 ) iargc character ( len = * ) :: node_file_name ! ! Get the number of command line arguments. ! arg_num = iargc ( ) ! ! If at least one command line argument, it's the node file name. ! if ( 1 <= arg_num ) then iarg = 1 call getarg ( iarg, node_file_name ) else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FILE_NAME_SPECIFICATION:' write ( *, '(a)' ) ' Please enter the name of the node file.' read ( *, '(a)' ) node_file_name end if ! ! If at least two command line arguments, the second is the triangulation file. ! if ( 2 <= arg_num ) then iarg = 2 call getarg ( iarg, element_file_name ) else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FILE_NAME_SPECIFICATION:' write ( *, '(a)' ) ' Please enter the name of the triangulation file.' read ( *, '(a)' ) element_file_name end if return end subroutine file_row_count ( input_file_name, row_num ) !*****************************************************************************80 ! !! FILE_ROW_COUNT counts the number of row records in a file. ! ! Discussion: ! ! It does not count lines that are blank, or that begin with a ! comment symbol '#'. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 06 March 2003 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) INPUT_FILE_NAME, the name of the input file. ! ! Output, integer ( kind = 4 ) ROW_NUM, the number of rows found. ! implicit none integer ( kind = 4 ) bad_num integer ( kind = 4 ) comment_num integer ( kind = 4 ) ierror character ( len = * ) input_file_name integer ( kind = 4 ) input_status integer ( kind = 4 ) input_unit character ( len = 255 ) line integer ( kind = 4 ) record_num integer ( kind = 4 ) row_num call get_unit ( input_unit ) open ( unit = input_unit, file = input_file_name, status = 'old', & iostat = input_status ) if ( input_status /= 0 ) then row_num = -1; ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FILE_ROW_COUNT - Fatal error!' write ( *, '(a,i8)' ) ' Could not open the input file "' // & trim ( input_file_name ) // '" on unit ', input_unit stop end if comment_num = 0 row_num = 0 record_num = 0 bad_num = 0 do read ( input_unit, '(a)', iostat = input_status ) line if ( input_status /= 0 ) then ierror = record_num exit end if record_num = record_num + 1 if ( line(1:1) == '#' ) then comment_num = comment_num + 1 cycle end if if ( len_trim ( line ) == 0 ) then comment_num = comment_num + 1 cycle end if row_num = row_num + 1 end do close ( unit = input_unit ) return end subroutine get_unit ( iunit ) !*****************************************************************************80 ! !! GET_UNIT returns a free FORTRAN unit number. ! ! Discussion: ! ! A "free" FORTRAN unit number is an integer between 1 and 99 which ! is not currently associated with an I/O device. A free FORTRAN unit ! number is needed in order to open a file with the OPEN command. ! ! If IUNIT = 0, then no free FORTRAN unit could be found, although ! all 99 units were checked (except for units 5, 6 and 9, which ! are commonly reserved for console I/O). ! ! Otherwise, IUNIT is an integer between 1 and 99, representing a ! free FORTRAN unit. Note that GET_UNIT assumes that units 5 and 6 ! are special, and will never return those values. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 18 September 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer ( kind = 4 ) IUNIT, the free unit number. ! implicit none integer ( kind = 4 ) i integer ( kind = 4 ) ios integer ( kind = 4 ) iunit logical lopen iunit = 0 do i = 1, 99 if ( i /= 5 .and. i /= 6 .and. i /= 9 ) then inquire ( unit = i, opened = lopen, iostat = ios ) if ( ios == 0 ) then if ( .not. lopen ) then iunit = i return end if end if end if end do return end function i4_huge ( ) !*****************************************************************************80 ! !! I4_HUGE returns a "huge" I4. ! ! Discussion: ! ! On an IEEE 32 bit machine, I4_HUGE should be 2^31 - 1, and its ! bit pattern should be ! ! 01111111111111111111111111111111 ! ! In this case, its numerical value is 2147483647. ! ! Using the Dec/Compaq/HP Alpha FORTRAN compiler FORT, I could ! use I4_HUGE() and HUGE interchangeably. ! ! However, when using the G95, the values returned by HUGE were ! not equal to 2147483647, apparently, and were causing severe ! and obscure errors in my random number generator, which needs to ! add I4_HUGE to the seed whenever the seed is negative. So I ! am backing away from invoking HUGE, whereas I4_HUGE is under ! my control. ! ! Explanation: because under G95 the default integer type is 64 bits! ! So HUGE ( 1 ) = a very very huge integer indeed, whereas ! I4_HUGE ( ) = the same old 32 bit big value. ! ! An I4 is an integer ( kind = 4 ) value. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 26 January 2007 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer ( kind = 4 ) I4_HUGE, a "huge" I4. ! implicit none integer ( kind = 4 ) i4 integer ( kind = 4 ) i4_huge i4_huge = 2147483647 return end function i4_modp ( i, j ) !*****************************************************************************80 ! !! I4_MODP returns the nonnegative remainder of integer division. ! ! Discussion: ! ! If ! NREM = I4_MODP ( I, J ) ! NMULT = ( I - NREM ) / J ! then ! I = J * NMULT + NREM ! where NREM is always nonnegative. ! ! The MOD function computes a result with the same sign as the ! quantity being divided. Thus, suppose you had an angle A, ! and you wanted to ensure that it was between 0 and 360. ! Then mod(A,360) would do, if A was positive, but if A ! was negative, your result would be between -360 and 0. ! ! On the other hand, I4_MODP(A,360) is between 0 and 360, always. ! ! Example: ! ! I J MOD I4_MODP Factorization ! ! 107 50 7 7 107 = 2 * 50 + 7 ! 107 -50 7 7 107 = -2 * -50 + 7 ! -107 50 -7 43 -107 = -3 * 50 + 43 ! -107 -50 -7 43 -107 = 3 * -50 + 43 ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 02 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) I, the number to be divided. ! ! Input, integer ( kind = 4 ) J, the number that divides I. ! ! Output, integer ( kind = 4 ) I4_MODP, the nonnegative remainder when I is ! divided by J. ! implicit none integer ( kind = 4 ) i integer ( kind = 4 ) i4_modp integer ( kind = 4 ) j if ( j == 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'I4_MODP - Fatal error!' write ( *, '(a,i8)' ) ' I4_MODP ( I, J ) called with J = ', j stop end if i4_modp = mod ( i, j ) if ( i4_modp < 0 ) then i4_modp = i4_modp + abs ( j ) end if return end function i4_wrap ( ival, ilo, ihi ) !*****************************************************************************80 ! !! I4_WRAP forces an integer to lie between given limits by wrapping. ! ! Example: ! ! ILO = 4, IHI = 8 ! ! I I4_WRAP ! ! -2 8 ! -1 4 ! 0 5 ! 1 6 ! 2 7 ! 3 8 ! 4 4 ! 5 5 ! 6 6 ! 7 7 ! 8 8 ! 9 4 ! 10 5 ! 11 6 ! 12 7 ! 13 8 ! 14 4 ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 19 August 2003 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) IVAL, an integer value. ! ! Input, integer ( kind = 4 ) ILO, IHI, the desired bounds for the integer ! value. ! ! Output, integer ( kind = 4 ) I4_WRAP, a "wrapped" version of IVAL. ! implicit none integer ( kind = 4 ) i4_modp integer ( kind = 4 ) i4_wrap integer ( kind = 4 ) ihi integer ( kind = 4 ) ilo integer ( kind = 4 ) ival integer ( kind = 4 ) jhi integer ( kind = 4 ) jlo integer ( kind = 4 ) wide jlo = min ( ilo, ihi ) jhi = max ( ilo, ihi ) wide = jhi - jlo + 1 if ( wide == 1 ) then i4_wrap = jlo else i4_wrap = jlo + i4_modp ( ival - jlo, wide ) end if return end subroutine i4col_compare ( m, n, a, i, j, isgn ) !*****************************************************************************80 ! !! I4COL_COMPARE compares columns I and J of an I4COL. ! ! Example: ! ! Input: ! ! M = 3, N = 4, I = 2, J = 4 ! ! A = ( ! 1 2 3 4 ! 5 6 7 8 ! 9 10 11 12 ) ! ! Output: ! ! ISGN = -1 ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 30 June 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) M, N, the number of rows and columns. ! ! Input, integer ( kind = 4 ) A(M,N), an array of N columns of vectors ! of length M. ! ! Input, integer ( kind = 4 ) I, J, the columns to be compared. ! I and J must be between 1 and N. ! ! Output, integer ( kind = 4 ) ISGN, the results of the comparison: ! -1, column I < column J, ! 0, column I = column J, ! +1, column J < column I. ! implicit none integer ( kind = 4 ) m integer ( kind = 4 ) n integer ( kind = 4 ) a(m,n) integer ( kind = 4 ) i integer ( kind = 4 ) isgn integer ( kind = 4 ) j integer ( kind = 4 ) k ! ! Check. ! if ( i < 1 .or. n < i ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'I4COL_COMPARE - Fatal error!' write ( *, '(a)' ) ' Column index I is out of bounds.' write ( *, '(a,i8)' ) ' I = ', i stop end if if ( j < 1 .or. n < j ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'I4COL_COMPARE - Fatal error!' write ( *, '(a)' ) ' Column index J is out of bounds.' write ( *, '(a,i8)' ) ' J = ', j stop end if isgn = 0 if ( i == j ) then return end if k = 1 do while ( k <= m ) if ( a(k,i) < a(k,j) ) then isgn = -1 return else if ( a(k,j) < a(k,i) ) then isgn = +1 return end if k = k + 1 end do return end subroutine i4col_sort_a ( m, n, a ) !*****************************************************************************80 ! !! I4COL_SORT_A ascending sorts an I4COL. ! ! Discussion: ! ! In lexicographic order, the statement "X < Y", applied to two real ! vectors X and Y of length M, means that there is some index I, with ! 1 <= I <= M, with the property that ! ! X(J) = Y(J) for J < I, ! and ! X(I) < Y(I). ! ! In other words, the first time they differ, X is smaller. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 25 September 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) M, the number of rows of A, and the length of ! a vector of data. ! ! Input, integer ( kind = 4 ) N, the number of columns of A. ! ! Input/output, integer ( kind = 4 ) A(M,N). ! On input, the array of N columns of M-vectors. ! On output, the columns of A have been sorted in ascending ! lexicographic order. ! implicit none integer ( kind = 4 ) m integer ( kind = 4 ) n integer ( kind = 4 ) a(m,n) integer ( kind = 4 ) i integer ( kind = 4 ) indx integer ( kind = 4 ) isgn integer ( kind = 4 ) j if ( m <= 0 ) then return end if if ( n <= 1 ) then return end if ! ! Initialize. ! i = 0 indx = 0 isgn = 0 j = 0 ! ! Call the external heap sorter. ! do call sort_heap_external ( n, indx, i, j, isgn ) ! ! Interchange the I and J objects. ! if ( 0 < indx ) then call i4col_swap ( m, n, a, i, j ) ! ! Compare the I and J objects. ! else if ( indx < 0 ) then call i4col_compare ( m, n, a, i, j, isgn ) else if ( indx == 0 ) then exit end if end do return end subroutine i4col_swap ( m, n, a, i, j ) !*****************************************************************************80 ! !! I4COL_SWAP swaps columns I and J of an I4COL. ! ! Example: ! ! Input: ! ! M = 3, N = 4, I = 2, J = 4 ! ! A = ( ! 1 2 3 4 ! 5 6 7 8 ! 9 10 11 12 ) ! ! Output: ! ! A = ( ! 1 4 3 2 ! 5 8 7 6 ! 9 12 11 10 ) ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 04 April 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) M, N, the number of rows and columns in ! the array. ! ! Input/output, integer ( kind = 4 ) A(M,N), an array of N columns of ! length M. ! ! Input, integer ( kind = 4 ) I, J, the columns to be swapped. ! implicit none integer ( kind = 4 ) m integer ( kind = 4 ) n integer ( kind = 4 ) a(m,n) integer ( kind = 4 ) col(m) integer ( kind = 4 ) i integer ( kind = 4 ) j if ( i < 1 .or. n < i .or. j < 1 .or. n < j ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'I4COL_SWAP - Fatal error!' write ( *, '(a)' ) ' I or J is out of bounds.' write ( *, '(a,i8)' ) ' I = ', i write ( *, '(a,i8)' ) ' J = ', j write ( *, '(a,i8)' ) ' N = ', n stop end if if ( i == j ) then return end if col(1:m) = a(1:m,i) a(1:m,i) = a(1:m,j) a(1:m,j) = col(1:m) return end subroutine i4mat_transpose_print_some ( m, n, a, ilo, jlo, ihi, jhi, title ) !*****************************************************************************80 ! !! I4MAT_TRANSPOSE_PRINT_SOME prints some of the transpose of an I4mat. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 09 February 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) M, N, the number of rows and columns. ! ! Input, integer ( kind = 4 ) A(M,N), an M by N matrix to be printed. ! ! Input, integer ( kind = 4 ) ILO, JLO, the first row and column to print. ! ! Input, integer ( kind = 4 ) IHI, JHI, the last row and column to print. ! ! Input, character ( len = * ) TITLE, a title. ! implicit none integer ( kind = 4 ), parameter :: incx = 10 integer ( kind = 4 ) m integer ( kind = 4 ) n integer ( kind = 4 ) a(m,n) character ( len = 7 ) ctemp(incx) integer ( kind = 4 ) i integer ( kind = 4 ) i2 integer ( kind = 4 ) i2hi integer ( kind = 4 ) i2lo integer ( kind = 4 ) ihi integer ( kind = 4 ) ilo integer ( kind = 4 ) inc integer ( kind = 4 ) j integer ( kind = 4 ) j2hi integer ( kind = 4 ) j2lo integer ( kind = 4 ) jhi integer ( kind = 4 ) jlo character ( len = * ) title write ( *, '(a)' ) ' ' write ( *, '(a)' ) trim ( title ) do i2lo = max ( ilo, 1 ), min ( ihi, m ), incx i2hi = i2lo + incx - 1 i2hi = min ( i2hi, m ) i2hi = min ( i2hi, ihi ) inc = i2hi + 1 - i2lo write ( *, '(a)' ) ' ' do i = i2lo, i2hi i2 = i + 1 - i2lo write ( ctemp(i2), '(i7)') i end do write ( *, '('' Row '',10a7)' ) ctemp(1:inc) write ( *, '(a)' ) ' Col' write ( *, '(a)' ) ' ' j2lo = max ( jlo, 1 ) j2hi = min ( jhi, n ) do j = j2lo, j2hi do i2 = 1, inc i = i2lo - 1 + i2 write ( ctemp(i2), '(i7)' ) a(i,j) end do write ( *, '(i5,1x,10a7)' ) j, ( ctemp(i), i = 1, inc ) end do end do return end subroutine i4vec_print ( n, a, title ) !*****************************************************************************80 ! !! I4VEC_PRINT prints an I4VEC. ! ! Discussion: ! ! An I4VEC is a vector of I4's. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 02 May 2010 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the number of components of the vector. ! ! Input, integer ( kind = 4 ) A(N), the vector to be printed. ! ! Input, character ( len = * ) TITLE, a title. ! implicit none integer ( kind = 4 ) n integer ( kind = 4 ) a(n) integer ( kind = 4 ) i character ( len = * ) title write ( *, '(a)' ) ' ' write ( *, '(a)' ) trim ( title ) write ( *, '(a)' ) ' ' do i = 1, n write ( *, '(2x,i8,a,2x,i12)' ) i, ':', a(i) end do return end subroutine i4vec2_compare ( n, a1, a2, i, j, isgn ) !*****************************************************************************80 ! !! I4VEC2_COMPARE compares pairs of integers stored in two vectors. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 22 October 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the number of data items. ! ! Input, integer ( kind = 4 ) A1(N), A2(N), contain the two components ! of each item. ! ! Input, integer ( kind = 4 ) I, J, the items to be compared. ! ! Output, integer ( kind = 4 ) ISGN, the results of the comparison: ! -1, item I < item J, ! 0, item I = item J, ! +1, item J < item I. ! implicit none integer ( kind = 4 ) n integer ( kind = 4 ) a1(n) integer ( kind = 4 ) a2(n) integer ( kind = 4 ) i integer ( kind = 4 ) isgn integer ( kind = 4 ) j isgn = 0 if ( a1(i) < a1(j) ) then isgn = -1 else if ( a1(i) == a1(j) ) then if ( a2(i) < a2(j) ) then isgn = -1 else if ( a2(i) < a2(j) ) then isgn = 0 else if ( a2(j) < a2(i) ) then isgn = +1 end if else if ( a1(j) < a1(i) ) then isgn = +1 end if return end subroutine i4vec2_sort_a ( n, a1, a2 ) !*****************************************************************************80 ! !! I4VEC2_SORT_A ascending sorts a vector of pairs of integers. ! ! Discussion: ! ! Each item to be sorted is a pair of integers (I,J), with the I ! and J values stored in separate vectors A1 and A2. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 25 September 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the number of items of data. ! ! Input/output, integer ( kind = 4 ) A1(N), A2(N), the data to be sorted. ! implicit none integer ( kind = 4 ) n integer ( kind = 4 ) a1(n) integer ( kind = 4 ) a2(n) integer ( kind = 4 ) i integer ( kind = 4 ) indx integer ( kind = 4 ) isgn integer ( kind = 4 ) j integer ( kind = 4 ) temp if ( n <= 1 ) then return end if ! ! Initialize. ! i = 0 indx = 0 isgn = 0 j = 0 ! ! Call the external heap sorter. ! do call sort_heap_external ( n, indx, i, j, isgn ) ! ! Interchange the I and J objects. ! if ( 0 < indx ) then temp = a1(i) a1(i) = a1(j) a1(j) = temp temp = a2(i) a2(i) = a2(j) a2(j) = temp ! ! Compare the I and J objects. ! else if ( indx < 0 ) then call i4vec2_compare ( n, a1, a2, i, j, isgn ) else if ( indx == 0 ) then exit end if end do return end subroutine i4mat_data_read ( input_filename, m, n, table ) !*****************************************************************************80 ! !! I4MAT_DATA_READ reads data from an I4MAT file. ! ! Discussion: ! ! An I4MAT is an array of I4's. ! ! The file may contain more than N points, but this routine ! will return after reading N points. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 27 January 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) INPUT_FILENAME, the name of the input file. ! ! Input, integer ( kind = 4 ) M, the spatial dimension. ! ! Input, integer ( kind = 4 ) N, the number of points. ! ! Output, integer ( kind = 4 ) TABLE(M,N), the table data. ! implicit none integer ( kind = 4 ) m integer ( kind = 4 ) n integer ( kind = 4 ) ierror character ( len = * ) input_filename integer ( kind = 4 ) input_status integer ( kind = 4 ) input_unit integer ( kind = 4 ) j character ( len = 255 ) line integer ( kind = 4 ) table(m,n) integer ( kind = 4 ) x(m) ierror = 0 call get_unit ( input_unit ) open ( unit = input_unit, file = input_filename, status = 'old', & iostat = input_status ) if ( input_status /= 0 ) then ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'I4MAT_DATA_READ - Fatal error!' write ( *, '(a,i8)' ) ' Could not open the input file "' // & trim ( input_filename ) // '" on unit ', input_unit stop end if j = 0 do while ( j < n ) read ( input_unit, '(a)', iostat = input_status ) line if ( input_status /= 0 ) then ierror = 2 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'I4MAT_DATA_READ - Fatal error!' write ( *, '(a)' ) ' Error while reading lines of data.' write ( *, '(a,i8)' ) ' Number of values expected per line M = ', m write ( *, '(a,i8)' ) ' Number of data lines read, J = ', j write ( *, '(a,i8)' ) ' Number of data lines needed, N = ', n stop end if if ( line(1:1) == '#' .or. len_trim ( line ) == 0 ) then cycle end if call s_to_i4vec ( line, m, x, ierror ) if ( ierror /= 0 ) then cycle end if j = j + 1 table(1:m,j) = x(1:m) end do close ( unit = input_unit ) return end subroutine i4mat_header_read ( input_filename, m, n ) !*****************************************************************************80 ! !! I4MAT_HEADER_READ reads the header from an I4MAT. ! ! Discussion: ! ! An I4MAT is an array of I4's. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 04 June 2004 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) INPUT_FILENAME, the name of the input file. ! ! Output, integer ( kind = 4 ) M, spatial dimension. ! ! Output, integer ( kind = 4 ) N, the number of points. ! implicit none character ( len = * ) input_filename integer ( kind = 4 ) m integer ( kind = 4 ) n call file_column_count ( input_filename, m ) if ( m <= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'I4MAT_HEADER_READ - Fatal error!' write ( *, '(a)' ) ' There was some kind of I/O problem while trying' write ( *, '(a)' ) ' to count the number of data columns in' write ( *, '(a)' ) ' the file "' // trim ( input_filename ) // '".' stop end if call file_row_count ( input_filename, n ) if ( n <= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'I4MAT_HEADER_READ - Fatal error!' write ( *, '(a)' ) ' There was some kind of I/O problem while trying' write ( *, '(a)' ) ' to count the number of data rows in' write ( *, '(a)' ) ' the file "' // trim ( input_filename ) // '".' stop end if return end subroutine quad_rule ( quad_num, quad_w, quad_xy ) !*****************************************************************************80 ! !! QUAD_RULE sets the quadrature rule for assembly. ! ! Discussion: ! ! The quadrature rule is given for a reference element. ! ! 0 <= X, ! 0 <= Y, and ! X + Y <= 1. ! ! ^ ! 1 | * ! | |\ ! Y | | \ ! | | \ ! 0 | *---* ! +-------> ! 0 X 1 ! ! The rules have the following precision: ! ! QUAD_NUM Precision ! ! 1 1 ! 3 2 ! 4 3 ! 6 4 ! 7 5 ! 9 6 ! 13 7 ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 18 July 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) QUAD_NUM, the number of quadrature nodes. ! ! Output, real ( kind = 8 ) QUAD_W(QUAD_NUM), the quadrature weights. ! ! Output, real ( kind = 8 ) QUAD_XY(2,QUAD_NUM), ! the coordinates of the quadrature nodes. ! implicit none integer ( kind = 4 ), parameter :: dim_num = 2 integer ( kind = 4 ) quad_num real ( kind = 8 ) a real ( kind = 8 ) b real ( kind = 8 ) c real ( kind = 8 ) d real ( kind = 8 ) e real ( kind = 8 ) f real ( kind = 8 ) g real ( kind = 8 ) h real ( kind = 8 ), dimension(quad_num) :: quad_w real ( kind = 8 ), dimension(dim_num,quad_num) :: quad_xy real ( kind = 8 ) t real ( kind = 8 ) u real ( kind = 8 ) v real ( kind = 8 ) w if ( quad_num == 1 ) then quad_xy(1:dim_num,1:quad_num) = reshape ( (/ & 1.0D+00 / 3.0D+00, 1.0D+00 / 3.0D+00 /), (/ dim_num, quad_num /) ) quad_w(1:quad_num) = 1.0D+00 else if ( quad_num == 3 ) then quad_xy(1:dim_num,1:quad_num) = reshape ( (/ & 0.5D+00, 0.0D+00, & 0.5D+00, 0.5D+00, & 0.0D+00, 0.5D+00 /), (/ dim_num, quad_num /) ) quad_w(1:quad_num) = 1.0D+00 / 3.0D+00 else if ( quad_num == 4 ) then a = 6.0D+00 / 30.0D+00 b = 10.0D+00 / 30.0D+00 c = 18.0D+00 / 30.0D+00 d = 25.0D+00 / 48.0D+00 e = -27.0D+00 / 48.0D+00 quad_xy(1:dim_num,1:quad_num) = reshape ( (/ & b, b, & c, a, & a, c, & a, a /), (/ dim_num, quad_num /) ) quad_w(1:quad_num) = (/ e, d, d, d /) else if ( quad_num == 6 ) then a = 0.816847572980459D+00 b = 0.091576213509771D+00 c = 0.108103018168070D+00 d = 0.445948490915965D+00 v = 0.109951743655322D+00 w = 0.223381589678011D+00 quad_xy(1:dim_num,1:quad_num) = reshape ( (/ & a, b, & b, a, & b, b, & c, d, & d, c, & d, d /), (/ dim_num, quad_num /) ) quad_w(1:quad_num) = (/ v, v, v, w, w, w /) else if ( quad_num == 7 ) then a = 1.0D+00 / 3.0D+00 b = ( 9.0D+00 + 2.0D+00 * sqrt ( 15.0D+00 ) ) / 21.0D+00 c = ( 6.0D+00 - sqrt ( 15.0D+00 ) ) / 21.0D+00 d = ( 9.0D+00 - 2.0D+00 * sqrt ( 15.0D+00 ) ) / 21.0D+00 e = ( 6.0D+00 + sqrt ( 15.0D+00 ) ) / 21.0D+00 u = 0.225D+00 v = ( 155.0D+00 - sqrt ( 15.0D+00 ) ) / 1200.0D+00 w = ( 155.0D+00 + sqrt ( 15.0D+00 ) ) / 1200.0D+00 quad_xy(1:dim_num,1:quad_num) = reshape ( (/ & a, a, & b, c, & c, b, & c, c, & d, e, & e, d, & e, e /), (/ dim_num, quad_num /) ) quad_w(1:quad_num) = (/ u, v, v, v, w, w, w /) else if ( quad_num == 9 ) then a = 0.124949503233232D+00 b = 0.437525248383384D+00 c = 0.797112651860071D+00 d = 0.165409927389841D+00 e = 0.037477420750088D+00 u = 0.205950504760887D+00 v = 0.063691414286223D+00 quad_xy(1:dim_num,1:quad_num) = reshape ( (/ & a, b, & b, a, & b, b, & c, d, & c, e, & d, c, & d, e, & e, c, & e, d /), (/ dim_num, quad_num /) ) quad_w(1:quad_num) = (/ u, u, u, v, v, v, v, v, v /) else if ( quad_num == 13 ) then h = 1.0D+00 / 3.0D+00 a = 0.479308067841923D+00 b = 0.260345966079038D+00 c = 0.869739794195568D+00 d = 0.065130102902216D+00 e = 0.638444188569809D+00 f = 0.312865496004875D+00 g = 0.048690315425316D+00 w = -0.149570044467670D+00 t = 0.175615257433204D+00 u = 0.053347235608839D+00 v = 0.077113760890257D+00 quad_xy(1:dim_num,1:quad_num) = reshape ( (/ & h, h, & a, b, & b, a, & b, b, & c, d, & d, c, & d, d, & e, f, & e, g, & f, e, & f, g, & g, e, & g, f /), (/ dim_num, quad_num /) ) quad_w(1:quad_num) = (/ w, t, t, t, u, u, u, v, v, v, v, v, v /) else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'QUAD_RULE - Fatal error!' write ( *, '(a,i8)' ) ' No rule is available of order QUAD_NUM = ', & quad_num stop end if return end subroutine r8mat_data_read ( input_filename, m, n, table ) !*****************************************************************************80 ! !! R8MAT_DATA_READ reads data from an R8MAT file. ! ! Discussion: ! ! An R8MAT is an array of R8 values. ! ! Discussion: ! ! The file may contain more than N points, but this routine will ! return after reading N of them. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 18 October 2008 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) INPUT_FILENAME, the name of the input file. ! ! Input, integer ( kind = 4 ) M, the spatial dimension. ! ! Input, integer ( kind = 4 ) N, the number of points. ! ! Output, real ( kind = 8 ) TABLE(M,N), the table data. ! implicit none integer ( kind = 4 ) m integer ( kind = 4 ) n integer ( kind = 4 ) ierror character ( len = * ) input_filename integer ( kind = 4 ) input_status integer ( kind = 4 ) input_unit integer ( kind = 4 ) j character ( len = 255 ) line real ( kind = 8 ) table(m,n) real ( kind = 8 ) x(m) ierror = 0 call get_unit ( input_unit ) open ( unit = input_unit, file = input_filename, status = 'old', & iostat = input_status ) if ( input_status /= 0 ) then ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'R8MAT_DATA_READ - Fatal error!' write ( *, '(a,i8)' ) ' Could not open the input file "' // & trim ( input_filename ) // '" on unit ', input_unit stop end if j = 0 do while ( j < n ) read ( input_unit, '(a)', iostat = input_status ) line if ( input_status /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'R8MAT_DATA_READ - Fatal error!' write ( *, '(a)' ) ' Error while reading lines of data.' write ( *, '(a,i8)' ) ' Number of values expected per line M = ', m write ( *, '(a,i8)' ) ' Number of data lines read, J = ', j write ( *, '(a,i8)' ) ' Number of data lines needed, N = ', n stop end if if ( line(1:1) == '#' .or. len_trim ( line ) == 0 ) then cycle end if call s_to_r8vec ( line, m, x, ierror ) if ( ierror /= 0 ) then cycle end if j = j + 1 table(1:m,j) = x(1:m) end do close ( unit = input_unit ) return end subroutine r8mat_header_read ( input_filename, m, n ) !*****************************************************************************80 ! !! R8MAT_HEADER_READ reads the header from an R8MAT file. ! ! Discussion: ! ! An R8MAT is an array of R8 values. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 07 September 2004 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) INPUT_FILENAME, the name of the input file. ! ! Output, integer ( kind = 4 ) M, spatial dimension. ! ! Output, integer ( kind = 4 ) N, the number of points. ! implicit none character ( len = * ) input_filename integer ( kind = 4 ) m integer ( kind = 4 ) n call file_column_count ( input_filename, m ) if ( m <= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'R8MAT_HEADER_READ - Fatal error!' write ( *, '(a)' ) ' There was some kind of I/O problem while trying' write ( *, '(a)' ) ' to count the number of data columns in' write ( *, '(a)' ) ' the file "' // trim ( input_filename ) // '".' stop end if call file_row_count ( input_filename, n ) if ( n <= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'R8MAT_HEADER_READ - Fatal error!' write ( *, '(a)' ) ' There was some kind of I/O problem while trying' write ( *, '(a)' ) ' to count the number of data rows in' write ( *, '(a)' ) ' the file "' // trim ( input_filename ) // '".' stop end if return end subroutine r8mat_transpose_print_some ( m, n, a, ilo, jlo, ihi, jhi, title ) !*****************************************************************************80 ! !! R8MAT_TRANSPOSE_PRINT_SOME prints some of an R8MAT, transposed. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 14 June 2004 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) M, N, the number of rows and columns. ! ! Input, real ( kind = 8 ) A(M,N), an M by N matrix to be printed. ! ! Input, integer ( kind = 4 ) ILO, JLO, the first row and column to print. ! ! Input, integer ( kind = 4 ) IHI, JHI, the last row and column to print. ! ! Input, character ( len = * ) TITLE, a title. ! implicit none integer ( kind = 4 ), parameter :: incx = 5 integer ( kind = 4 ) m integer ( kind = 4 ) n real ( kind = 8 ) a(m,n) character ( len = 14 ) ctemp(incx) integer ( kind = 4 ) i integer ( kind = 4 ) i2 integer ( kind = 4 ) i2hi integer ( kind = 4 ) i2lo integer ( kind = 4 ) ihi integer ( kind = 4 ) ilo integer ( kind = 4 ) inc integer ( kind = 4 ) j integer ( kind = 4 ) j2hi integer ( kind = 4 ) j2lo integer ( kind = 4 ) jhi integer ( kind = 4 ) jlo character ( len = * ) title write ( *, '(a)' ) ' ' write ( *, '(a)' ) trim ( title ) do i2lo = max ( ilo, 1 ), min ( ihi, m ), incx i2hi = i2lo + incx - 1 i2hi = min ( i2hi, m ) i2hi = min ( i2hi, ihi ) inc = i2hi + 1 - i2lo write ( *, '(a)' ) ' ' do i = i2lo, i2hi i2 = i + 1 - i2lo write ( ctemp(i2), '(i7,7x)') i end do write ( *, '('' Row '',5a14)' ) ctemp(1:inc) write ( *, '(a)' ) ' Col' write ( *, '(a)' ) ' ' j2lo = max ( jlo, 1 ) j2hi = min ( jhi, n ) do j = j2lo, j2hi do i2 = 1, inc i = i2lo - 1 + i2 write ( ctemp(i2), '(g14.6)' ) a(i,j) end do write ( *, '(i5,1x,5a14)' ) j, ( ctemp(i), i = 1, inc ) end do end do return end subroutine r8mat_write ( output_filename, m, n, table ) !*****************************************************************************80 ! !! R8MAT_WRITE writes an R8MAT file. ! ! Discussion: ! ! An R8MAT is an array of R8 values. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 31 May 2009 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) OUTPUT_FILENAME, the output file name. ! ! Input, integer ( kind = 4 ) M, the spatial dimension. ! ! Input, integer ( kind = 4 ) N, the number of points. ! ! Input, real ( kind = 8 ) TABLE(M,N), the table data. ! implicit none integer ( kind = 4 ) m integer ( kind = 4 ) n integer ( kind = 4 ) j character ( len = * ) output_filename integer ( kind = 4 ) output_status integer ( kind = 4 ) output_unit character ( len = 30 ) string real ( kind = 8 ) table(m,n) ! ! Open the file. ! call get_unit ( output_unit ) open ( unit = output_unit, file = output_filename, & status = 'replace', iostat = output_status ) if ( output_status /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'R8MAT_WRITE - Fatal error!' write ( *, '(a,i8)' ) ' Could not open the output file "' // & trim ( output_filename ) // '" on unit ', output_unit output_unit = -1 stop end if ! ! Create a format string. ! ! For less precision in the output file, try: ! ! '(', m, 'g', 14, '.', 6, ')' ! if ( 0 < m .and. 0 < n ) then write ( string, '(a1,i8,a1,i8,a1,i8,a1)' ) '(', m, 'g', 24, '.', 16, ')' ! ! Write the data. ! do j = 1, n write ( output_unit, string ) table(1:m,j) end do end if ! ! Close the file. ! close ( unit = output_unit ) return end subroutine r8vec_print_some ( n, a, i_lo, i_hi, title ) !*****************************************************************************80 ! !! R8VEC_PRINT_SOME prints "some" of an R8VEC. ! ! Discussion: ! ! An R8VEC is a vector of R8 values. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 16 October 2006 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the number of entries of the vector. ! ! Input, real ( kind = 8 ) A(N), the vector to be printed. ! ! Input, integer ( kind = 4 ) I_LO, I_HI, the first and last indices ! to print. The routine expects 1 <= I_LO <= I_HI <= N. ! ! Input, character ( len = * ) TITLE, an optional title. ! implicit none integer ( kind = 4 ) n real ( kind = 8 ) a(n) integer ( kind = 4 ) i integer ( kind = 4 ) i_hi integer ( kind = 4 ) i_lo character ( len = * ) title if ( 0 < len_trim ( title ) ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) trim ( title ) end if write ( *, '(a)' ) ' ' do i = max ( i_lo, 1 ), min ( i_hi, n ) write ( *, '(2x,i8,2x,g14.8)' ) i, a(i) end do return end subroutine r8vec_uniform_01 ( n, seed, r ) !*****************************************************************************80 ! !! R8VEC_UNIFORM_01 returns a unit pseudorandom R8VEC. ! ! Discussion: ! ! An R8VEC is a vector of real ( kind = 8 ) values. ! ! For now, the input quantity SEED is an integer variable. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 31 May 2007 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Paul Bratley, Bennett Fox, Linus Schrage, ! A Guide to Simulation, ! Second Edition, ! Springer, 1987, ! ISBN: 0387964673, ! LC: QA76.9.C65.B73. ! ! Bennett Fox, ! Algorithm 647: ! Implementation and Relative Efficiency of Quasirandom ! Sequence Generators, ! ACM Transactions on Mathematical Software, ! Volume 12, Number 4, December 1986, pages 362-376. ! ! Pierre L'Ecuyer, ! Random Number Generation, ! in Handbook of Simulation, ! edited by Jerry Banks, ! Wiley, 1998, ! ISBN: 0471134031, ! LC: T57.62.H37. ! ! Peter Lewis, Allen Goodman, James Miller, ! A Pseudo-Random Number Generator for the System/360, ! IBM Systems Journal, ! Volume 8, 1969, pages 136-143. ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the number of entries in the vector. ! ! Input/output, integer ( kind = 4 ) SEED, the "seed" value, which ! should NOT be 0. On output, SEED has been updated. ! ! Output, real ( kind = 8 ) R(N), the vector of pseudorandom values. ! implicit none integer ( kind = 4 ) n integer ( kind = 4 ) i integer ( kind = 4 ) i4_huge integer ( kind = 4 ) k integer ( kind = 4 ) seed real ( kind = 8 ) r(n) if ( seed == 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'R8VEC_UNIFORM_01 - Fatal error!' write ( *, '(a)' ) ' Input value of SEED = 0.' stop end if do i = 1, n k = seed / 127773 seed = 16807 * ( seed - k * 127773 ) - k * 2836 if ( seed < 0 ) then seed = seed + i4_huge ( ) end if r(i) = real ( seed, kind = 8 ) * 4.656612875D-10 end do return end subroutine reference_to_physical_t3 ( t, n, ref, phy ) !*****************************************************************************80 ! !! REFERENCE_TO_PHYSICAL_T3 maps reference points to physical points. ! ! Discussion: ! ! Given the vertices of an order 3 physical triangle and a point ! (XSI,ETA) in the reference triangle, the routine computes the value ! of the corresponding image point (X,Y) in physical space. ! ! Note that this routine may also be appropriate for an order 6 ! triangle, if the mapping between reference and physical space ! is linear. This implies, in particular, that the sides of the ! image triangle are straight and that the "midside" nodes in the ! physical triangle are literally halfway along the sides of ! the physical triangle. ! ! Reference Element T3: ! ! | ! 1 3 ! | |\ ! | | \ ! S | \ ! | | \ ! | | \ ! 0 1-----2 ! | ! +--0--R--1--> ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 24 June 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) T(2,3), the coordinates of the vertices. ! The vertices are assumed to be the images of (0,0), (1,0) and ! (0,1) respectively. ! ! Input, integer ( kind = 4 ) N, the number of objects to transform. ! ! Input, real ( kind = 8 ) REF(2,N), points in the reference triangle. ! ! Output, real ( kind = 8 ) PHY(2,N), corresponding points in the ! physical triangle. ! implicit none integer ( kind = 4 ) n integer ( kind = 4 ) i real ( kind = 8 ) phy(2,n) real ( kind = 8 ) ref(2,n) real ( kind = 8 ) t(2,3) do i = 1, 2 phy(i,1:n) = t(i,1) * ( 1.0D+00 - ref(1,1:n) - ref(2,1:n) ) & + t(i,2) * ref(1,1:n) & + t(i,3) * ref(2,1:n) end do return end subroutine s_to_i4 ( s, ival, ierror, length ) !*****************************************************************************80 ! !! S_TO_I4 reads an I4 from a string. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 28 June 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) S, a string to be examined. ! ! Output, integer ( kind = 4 ) IVAL, the integer value read from the string. ! If the string is blank, then IVAL will be returned 0. ! ! Output, integer ( kind = 4 ) IERROR, an error flag. ! 0, no error. ! 1, an error occurred. ! ! Output, integer ( kind = 4 ) LENGTH, the number of characters of S ! used to make IVAL. ! implicit none character c integer ( kind = 4 ) i integer ( kind = 4 ) ierror integer ( kind = 4 ) isgn integer ( kind = 4 ) istate integer ( kind = 4 ) ival integer ( kind = 4 ) length 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 length = 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 length = len_trim ( s ) else ierror = 1 length = 0 end if return end subroutine s_to_i4vec ( s, n, ivec, ierror ) !*****************************************************************************80 ! !! S_TO_I4VEC reads an I4VEC from a string. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 08 October 2003 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) S, the string to be read. ! ! Input, integer ( kind = 4 ) N, the number of values expected. ! ! Output, integer ( kind = 4 ) IVEC(N), the values read from the string. ! ! Output, integer ( kind = 4 ) IERROR, error flag. ! 0, no errors occurred. ! -K, could not read data for entries -K through N. ! implicit none integer ( kind = 4 ) n integer ( kind = 4 ) i integer ( kind = 4 ) ierror integer ( kind = 4 ) ilo integer ( kind = 4 ) ivec(n) integer ( kind = 4 ) length character ( len = * ) s i = 0 ierror = 0 ilo = 1 do while ( i < n ) i = i + 1 call s_to_i4 ( s(ilo:), ivec(i), ierror, length ) if ( ierror /= 0 ) then ierror = -i exit end if ilo = ilo + length end do return end subroutine s_to_r8 ( s, dval, ierror, length ) !*****************************************************************************80 ! !! S_TO_R8 reads an R8 from a string. ! ! Discussion: ! ! The 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 number. ! ! Legal input is: ! ! 1 blanks, ! 2 '+' or '-' sign, ! 2.5 blanks ! 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. ! ! Example: ! ! S DVAL ! ! '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) ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 07 September 2004 ! ! 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 ( kind = 8 ) DVAL, the value read from the string. ! ! Output, integer ( kind = 4 ) 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 ( kind = 4 ) LENGTH, the number of characters read ! to form the number, including any terminating ! characters such as a trailing comma or blanks. ! implicit none logical ch_eqi character c real ( kind = 8 ) dval integer ( kind = 4 ) ierror integer ( kind = 4 ) ihave integer ( kind = 4 ) isgn integer ( kind = 4 ) iterm integer ( kind = 4 ) jbot integer ( kind = 4 ) jsgn integer ( kind = 4 ) jtop integer ( kind = 4 ) length integer ( kind = 4 ) nchar integer ( kind = 4 ) ndig real ( kind = 8 ) rbot real ( kind = 8 ) rexp real ( kind = 8 ) rtop character ( len = * ) s nchar = len_trim ( s ) ierror = 0 dval = 0.0D+00 length = -1 isgn = 1 rtop = 0 rbot = 1 jsgn = 1 jtop = 0 jbot = 1 ihave = 1 iterm = 0 do length = length + 1 if ( nchar < length+1 ) then exit end if c = s(length+1:length+1) ! ! Blank character. ! if ( c == ' ' ) then if ( ihave == 2 ) then else if ( ihave == 6 .or. ihave == 7 ) then iterm = 1 else if ( 1 < ihave ) then ihave = 11 end if ! ! Comma. ! else if ( c == ',' .or. c == ';' ) then if ( ihave /= 1 ) then iterm = 1 ihave = 12 length = length + 1 end if ! ! Minus sign. ! else if ( c == '-' ) then if ( ihave == 1 ) then ihave = 2 isgn = -1 else if ( ihave == 6 ) then ihave = 7 jsgn = -1 else iterm = 1 end if ! ! Plus sign. ! else if ( c == '+' ) then if ( ihave == 1 ) then ihave = 2 else if ( ihave == 6 ) then ihave = 7 else iterm = 1 end if ! ! Decimal point. ! else if ( c == '.' ) then if ( ihave < 4 ) then ihave = 4 else if ( 6 <= ihave .and. ihave <= 8 ) then ihave = 9 else iterm = 1 end if ! ! Scientific notation 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. lle ( '0', c ) .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.0D+00 * rtop + real ( ndig, kind = 8 ) else if ( ihave == 5 ) then rtop = 10.0D+00 * rtop + real ( ndig, kind = 8 ) rbot = 10.0D+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 ) 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 LENGTH is equal to NCHAR. ! if ( iterm /= 1 .and. length+1 == nchar ) then length = 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 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'S_TO_R8 - Serious error!' write ( *, '(a)' ) ' Illegal or nonnumeric input:' write ( *, '(a)' ) ' ' // trim ( s ) return end if ! ! Number seems OK. Form it. ! if ( jtop == 0 ) then rexp = 1.0D+00 else if ( jbot == 1 ) then rexp = 10.0D+00 ** ( jsgn * jtop ) else rexp = 10.0D+00 ** ( real ( jsgn * jtop, kind = 8 ) & / real ( jbot, kind = 8 ) ) end if end if dval = real ( isgn, kind = 8 ) * rexp * rtop / rbot return end subroutine s_to_r8vec ( s, n, rvec, ierror ) !*****************************************************************************80 ! !! S_TO_R8VEC reads an R8VEC from a string. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 07 September 2004 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) S, the string to be read. ! ! Input, integer ( kind = 4 ) N, the number of values expected. ! ! Output, real ( kind = 8 ) RVEC(N), the values read from the string. ! ! Output, integer ( kind = 4 ) IERROR, error flag. ! 0, no errors occurred. ! -K, could not read data for entries -K through N. ! implicit none integer ( kind = 4 ) n integer ( kind = 4 ) i integer ( kind = 4 ) ierror integer ( kind = 4 ) ilo integer ( kind = 4 ) lchar real ( kind = 8 ) rvec(n) character ( len = * ) s i = 0 ierror = 0 ilo = 1 do while ( i < n ) i = i + 1 call s_to_r8 ( s(ilo:), rvec(i), ierror, lchar ) if ( ierror /= 0 ) then ierror = -i exit end if ilo = ilo + lchar end do return end subroutine s_word_count ( s, nword ) !*****************************************************************************80 ! !! S_WORD_COUNT counts the number of "words" in a string. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 14 April 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) S, the string to be examined. ! ! Output, integer ( kind = 4 ) NWORD, the number of "words" in the string. ! Words are presumed to be separated by one or more blanks. ! implicit none logical blank integer ( kind = 4 ) i integer ( kind = 4 ) lens integer ( kind = 4 ) nword character ( len = * ) s nword = 0 lens = len ( s ) if ( lens <= 0 ) then return end if blank = .true. do i = 1, lens if ( s(i:i) == ' ' ) then blank = .true. else if ( blank ) then nword = nword + 1 blank = .false. end if end do return end subroutine solution_evaluate ( xy, t, node_u, u, dudx, dudy ) !*****************************************************************************80 ! !! SOLUTION_EVALUATE evaluates the solution at a point in an element. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 06 April 2004 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) XY(2), the point where the solution is ! to be evaluated. ! ! Input, real ( kind = 8 ) T(2,3), the coordinates of the vertices ! of the triangle which contains XY. ! ! Input, real ( kind = 8 ) NODE_U(3), the value of the solution ! at the nodes of the triangle. ! ! Output, real ( kind = 8 ) U, DUDX, DUDY, the solution and its X and ! Y derivatives at XY. ! implicit none real ( kind = 8 ) b real ( kind = 8 ) dbdx real ( kind = 8 ) dbdy real ( kind = 8 ) dudx real ( kind = 8 ) dudy integer ( kind = 4 ) i real ( kind = 8 ) node_u(3) real ( kind = 8 ) t(2,3) real ( kind = 8 ) u real ( kind = 8 ) xy(2) u = 0.0D+00 dudx = 0.0D+00 dudy = 0.0D+00 do i = 1, 3 call basis_one_t3 ( t, i, xy, b, dbdx, dbdy ) u = u + node_u(i) * b dudx = dudx + node_u(i) * dbdx dudy = dudy + node_u(i) * dbdy end do return end subroutine solve_cg ( n, diag, nz_num, ia, ja, a, b, x ) !*****************************************************************************80 ! !! SOLVE_CG solves a linear system using the conjugate gradient method. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 26 January 2013 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the number of nodes. ! ! Input, integer ( kind = 4 ) DIAG(N), contains for each index 1 <= I <= N, ! the unique index J such that IA(J) = JA(J) = I. ! ! Input, integer ( kind = 4 ) NZ_NUM, the number of nonzero entries. ! ! Input, integer ( kind = 4 ) IA(NZ_NUM), JA(NZ_NUM), the row and column ! indices of the nonzero entries. ! ! Input, real ( kind = 8 ) A(NZ_NUM), the nonzero entries of the matrix. ! ! Input, real ( kind = 8 ) B(N), the right hand side. ! ! Output, real ( kind = 8 ) X(N), the solution of the linear system. ! implicit none integer ( kind = 4 ) n integer ( kind = 4 ) nz_num real ( kind = 8 ) a(nz_num) real ( kind = 8 ) aii real ( kind = 8 ) b(n) real ( kind = 8 ) bnrm2 integer ( kind = 4 ) diag(n) integer ( kind = 4 ) i integer ( kind = 4 ) ia(nz_num) integer ( kind = 4 ) it integer ( kind = 4 ) it_max integer ( kind = 4 ) j integer ( kind = 4 ) ja(nz_num) integer ( kind = 4 ) job integer ( kind = 4 ) k real ( kind = 8 ) p(n) real ( kind = 8 ) q(n) real ( kind = 8 ) r(n) real ( kind = 8 ) rnrm2 real ( kind = 8 ) tol real ( kind = 8 ) x(n) real ( kind = 8 ) z(n) it = 0 it_max = 100 tol = 1.0D-08 bnrm2 = sqrt ( sum ( b(1:n)**2 ) ) x(1:n) = b(1:n) / a(diag(1:n)) write ( *, '(a)' ) '' write ( *, '(a)' ) ' Step Residual' write ( *, '(a)' ) '' job = 1 do call cg_rc ( n, b, x, r, z, p, q, job ) ! ! Compute q = A * p. ! if ( job == 1 ) then q(1:n) = 0.0D+00 do k = 1, nz_num i = ia(k) j = ja(k) q(i) = q(i) + a(k) * p(j) end do ! ! Solve M * z = r. ! else if ( job == 2 ) then z(1:n) = r(1:n) / a(diag(1:n)) ! ! Compute r = r - A * x. ! else if ( job == 3 ) then do k = 1, nz_num i = ia(k) j = ja(k) r(i) = r(i) - a(k) * x(j) end do ! ! Stopping test. ! else if ( job == 4 ) then rnrm2 = sqrt ( sum ( r(1:n)**2 ) ) if ( bnrm2 == 0.0D+00 ) then if ( rnrm2 <= tol ) then exit end if else if ( rnrm2 <= tol * bnrm2 ) then exit end if end if it = it + 1 write ( *, '(2x,i4,2x,g14.6)' ) it, rnrm2 if ( it_max <= it ) then write ( *, '(a)' ) '' write ( *, '(a)' ) ' Iteration limit exceeded.' write ( *, '(a)' ) ' Terminating early.' exit end if end if job = 2 end do write ( *, '(a)' ) '' write ( *, '(a,i4)' ) ' Number of iterations was ', it write ( *, '(a,g14.6)' ) ' Estimated error is ', rnrm2 return end subroutine sort_heap_external ( n, indx, i, j, isgn ) !*****************************************************************************80 ! !! SORT_HEAP_EXTERNAL externally sorts a list of items into ascending order. ! ! Discussion: ! ! The actual list of data is not passed to the routine. Hence this ! routine may be used to sort integers, reals, numbers, names, ! dates, shoe sizes, and so on. After each call, the routine asks ! the user to compare or interchange two items, until a special ! return value signals that the sorting is completed. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 05 February 2004 ! ! Author: ! ! Original FORTRAN77 version by Albert Nijenhuis, Herbert Wilf. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! A Nijenhuis and H Wilf, ! Combinatorial Algorithms, ! Academic Press, 1978, second edition, ! ISBN 0-12-519260-6. ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the number of items to be sorted. ! ! Input/output, integer ( kind = 4 ) INDX, the main communication signal. ! ! The user must set INDX to 0 before the first call. ! Thereafter, the user should not change the value of INDX until ! the sorting is done. ! ! On return, if INDX is ! ! greater than 0, ! * interchange items I and J; ! * call again. ! ! less than 0, ! * compare items I and J; ! * set ISGN = -1 if I < J, ISGN = +1 if J < I; ! * call again. ! ! equal to 0, the sorting is done. ! ! Output, integer ( kind = 4 ) I, J, the indices of two items. ! On return with INDX positive, elements I and J should be interchanged. ! On return with INDX negative, elements I and J should be compared, and ! the result reported in ISGN on the next call. ! ! Input, integer ( kind = 4 ) ISGN, results of comparison of elements ! I and J. (Used only when the previous call returned INDX less than 0). ! ISGN <= 0 means I is less than or equal to J; ! 0 <= ISGN means I is greater than or equal to J. ! implicit none integer ( kind = 4 ) i integer ( kind = 4 ), save :: i_save = 0 integer ( kind = 4 ) indx integer ( kind = 4 ) isgn integer ( kind = 4 ) j integer ( kind = 4 ), save :: j_save = 0 integer ( kind = 4 ), save :: k = 0 integer ( kind = 4 ), save :: k1 = 0 integer ( kind = 4 ) n integer ( kind = 4 ), save :: n1 = 0 ! ! INDX = 0: This is the first call. ! if ( indx == 0 ) then i_save = 0 j_save = 0 k = n / 2 k1 = k n1 = n ! ! INDX < 0: The user is returning the results of a comparison. ! else if ( indx < 0 ) then if ( indx == -2 ) then if ( isgn < 0 ) then i_save = i_save + 1 end if j_save = k1 k1 = i_save indx = -1 i = i_save j = j_save return end if if ( 0 < isgn ) then indx = 2 i = i_save j = j_save return end if if ( k <= 1 ) then if ( n1 == 1 ) then i_save = 0 j_save = 0 indx = 0 else i_save = n1 n1 = n1 - 1 j_save = 1 indx = 1 end if i = i_save j = j_save return end if k = k - 1 k1 = k ! ! 0 < INDX, the user was asked to make an interchange. ! else if ( indx == 1 ) then k1 = k end if do i_save = 2 * k1 if ( i_save == n1 ) then j_save = k1 k1 = i_save indx = -1 i = i_save j = j_save return else if ( i_save <= n1 ) then j_save = i_save + 1 indx = -2 i = i_save j = j_save return end if if ( k <= 1 ) then exit end if k = k - 1 k1 = k end do if ( n1 == 1 ) then i_save = 0 j_save = 0 indx = 0 i = i_save j = j_save else i_save = n1 n1 = n1 - 1 j_save = 1 indx = 1 i = i_save j = j_save end if return end subroutine timestamp ( ) !*****************************************************************************80 ! !! TIMESTAMP prints the current YMDHMS date as a time stamp. ! ! Example: ! ! May 31 2001 9:45:54.872 AM ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 31 May 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! None ! implicit none character ( len = 8 ) ampm integer ( kind = 4 ) d character ( len = 8 ) date integer ( kind = 4 ) h integer ( kind = 4 ) m integer ( kind = 4 ) mm character ( len = 9 ), parameter, dimension(12) :: month = (/ & 'January ', 'February ', 'March ', 'April ', & 'May ', 'June ', 'July ', 'August ', & 'September', 'October ', 'November ', 'December ' /) integer ( kind = 4 ) n integer ( kind = 4 ) s character ( len = 10 ) time integer ( kind = 4 ) values(8) integer ( kind = 4 ) 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 function triangle_area_2d ( t ) !*****************************************************************************80 ! !! TRIANGLE_AREA_2D computes the area of a triangle in 2D. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 17 December 2004 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) T(2,3), the triangle vertices. ! ! Output, real ( kind = 8 ) TRIANGLE_AREA_2D, the absolute area of ! the triangle. ! implicit none integer ( kind = 4 ), parameter :: dim_num = 2 real ( kind = 8 ) t(dim_num,3) real ( kind = 8 ) triangle_area_2d triangle_area_2d = 0.5D+00 * abs ( & t(1,1) * ( t(2,2) - t(2,3) ) & + t(1,2) * ( t(2,3) - t(2,1) ) & + t(1,3) * ( t(2,1) - t(2,2) ) ) return end subroutine triangulation_order3_adj_count ( node_num, triangle_num, & triangle_node, triangle_neighbor, adj_num, adj_col ) !*****************************************************************************80 ! !! TRIANGULATION_ORDER3_ADJ_COUNT counts adjacencies in a triangulation. ! ! Discussion: ! ! This routine is called to count the adjacencies, so that the ! appropriate amount of memory can be set aside for storage when ! the adjacency structure is created. ! ! The triangulation is assumed to involve 3-node triangles. ! ! Two nodes are "adjacent" if they are both nodes in some triangle. ! Also, a node is considered to be adjacent to itself. ! ! Diagram: ! ! 3 ! s |\ ! i | \ ! d | \ ! e | \ side 2 ! | \ ! 3 | \ ! | \ ! 1-------2 ! ! side 1 ! ! The local node numbering ! ! ! 21-22-23-24-25 ! |\ |\ |\ |\ | ! | \| \| \| \| ! 16-17-18-19-20 ! |\ |\ |\ |\ | ! | \| \| \| \| ! 11-12-13-14-15 ! |\ |\ |\ |\ | ! | \| \| \| \| ! 6--7--8--9-10 ! |\ |\ |\ |\ | ! | \| \| \| \| ! 1--2--3--4--5 ! ! A sample grid. ! ! ! Below, we have a chart that summarizes the adjacency relationships ! in the sample grid. On the left, we list the node, and its neighbors, ! with an asterisk to indicate the adjacency of the node to itself ! (in some cases, you want to count this self adjacency and in some ! you don't). On the right, we list the number of adjancencies to ! lower-indexed nodes, to the node itself, to higher-indexed nodes, ! the total number of adjacencies for this node, and the location ! of the first and last entries required to list this set of adjacencies ! in a single list of all the adjacencies. ! ! N Adjacencies Below Self Above Total First Last ! ! -- -- -- -- -- -- -- -- -- -- -- -- --- 0 ! 1: * 2 6 0 1 2 3 1 3 ! 2: 1 * 3 6 7 1 1 3 5 4 8 ! 3: 2 * 4 7 8 1 1 3 5 9 13 ! 4: 3 * 5 8 9 1 1 3 5 14 18 ! 5: 4 * 9 10 1 1 2 4 19 22 ! 6: 1 2 * 7 11 2 1 2 5 23 27 ! 7: 2 3 6 * 8 11 12 3 1 3 7 28 34 ! 8: 3 4 7 * 9 12 13 3 1 3 7 35 41 ! 9: 4 5 8 * 10 13 14 3 1 3 7 42 48 ! 10: 5 9 * 14 15 2 1 2 5 49 53 ! 11: 6 7 * 12 16 2 1 2 5 54 58 ! 12: 7 8 11 * 13 16 17 3 1 3 7 59 65 ! 13: 8 9 12 * 14 17 18 3 1 3 7 66 72 ! 14: 9 10 13 * 15 18 19 3 1 3 7 73 79 ! 15: 10 14 * 19 20 2 1 2 5 80 84 ! 16: 11 12 * 17 21 2 1 2 5 85 89 ! 17: 12 13 16 * 18 21 22 3 1 3 7 90 96 ! 18: 13 14 17 * 19 22 23 3 1 3 7 97 103 ! 19: 14 15 18 * 20 23 24 3 1 3 7 104 110 ! 20: 15 19 * 24 25 2 1 2 5 111 115 ! 21: 16 17 * 22 2 1 1 4 116 119 ! 22: 17 18 21 * 23 3 1 1 5 120 124 ! 23: 18 19 22 * 24 3 1 1 5 125 129 ! 24: 19 20 23 * 25 3 1 1 5 130 134 ! 25: 20 24 * 2 1 0 3 135 137 ! -- -- -- -- -- -- -- -- -- -- -- -- 138 --- ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 24 August 2006 ! ! Author: ! ! John Burkardt ! ! Parameters ! ! Input, integer ( kind = 4 ) NODE_NUM, the number of nodes. ! ! Input, integer ( kind = 4 ) TRIANGLE_NUM, the number of triangles. ! ! Input, integer ( kind = 4 ) TRIANGLE_NODE(3,TRIANGLE_NUM), lists the ! nodes that make up each triangle, in counterclockwise order. ! ! Input, integer ( kind = 4 ) TRIANGLE_NEIGHBOR(3,TRIANGLE_NUM), for each ! side of a triangle, lists the neighboring triangle, or -1 if there is ! no neighbor. ! ! Output, integer ( kind = 4 ) ADJ_NUM, the number of adjacencies. ! ! Output, integer ( kind = 4 ) ADJ_COL(NODE_NUM+1). Information about ! column J is stored in entries ADJ_COL(J) through ADJ_COL(J+1)-1 of ADJ. ! implicit none integer ( kind = 4 ) node_num integer ( kind = 4 ) triangle_num integer ( kind = 4 ), parameter :: triangle_order = 3 integer ( kind = 4 ) adj_num integer ( kind = 4 ) adj_col(node_num+1) integer ( kind = 4 ) i integer ( kind = 4 ) n1 integer ( kind = 4 ) n2 integer ( kind = 4 ) n3 integer ( kind = 4 ) triangle integer ( kind = 4 ) triangle2 integer ( kind = 4 ) triangle_neighbor(3,triangle_num) integer ( kind = 4 ) triangle_node(triangle_order,triangle_num) adj_num = 0 ! ! Set every node to be adjacent to itself. ! adj_col(1:node_num) = 1 ! ! Examine each triangle. ! do triangle = 1, triangle_num n1 = triangle_node(1,triangle) n2 = triangle_node(2,triangle) n3 = triangle_node(3,triangle) ! ! Add edge (1,2) if this is the first occurrence, ! that is, if the edge (1,2) is on a boundary (TRIANGLE2 <= 0) ! or if this triangle is the first of the pair in which the edge ! occurs (TRIANGLE < TRIANGLE2). ! triangle2 = triangle_neighbor(1,triangle) if ( triangle2 < 0 .or. triangle < triangle2 ) then adj_col(n1) = adj_col(n1) + 1 adj_col(n2) = adj_col(n2) + 1 end if ! ! Add edge (2,3). ! triangle2 = triangle_neighbor(2,triangle) if ( triangle2 < 0 .or. triangle < triangle2 ) then adj_col(n2) = adj_col(n2) + 1 adj_col(n3) = adj_col(n3) + 1 end if ! ! Add edge (3,1). ! triangle2 = triangle_neighbor(3,triangle) if ( triangle2 < 0 .or. triangle < triangle2 ) then adj_col(n1) = adj_col(n1) + 1 adj_col(n3) = adj_col(n3) + 1 end if end do ! ! We used ADJ_COL to count the number of entries in each column. ! Convert it to pointers into the ADJ array. ! adj_col(2:node_num+1) = adj_col(1:node_num) adj_col(1) = 1 do i = 2, node_num+1 adj_col(i) = adj_col(i-1) + adj_col(i) end do adj_num = adj_col(node_num+1) - 1 return end subroutine triangulation_order3_adj_set2 ( node_num, triangle_num, & triangle_node, triangle_neighbor, adj_num, adj_col, ia, ja ) !*****************************************************************************80 ! !! TRIANGULATION_ORDER3_ADJ_SET2 sets adjacencies in a triangulation. ! ! Discussion: ! ! This routine is called to set up the arrays IA and JA that ! record which nodes are adjacent in a triangulation. ! ! The triangulation is assumed to involve 3-node triangles. ! ! Two nodes are "adjacent" if they are both nodes in some triangle. ! Also, a node is considered to be adjacent to itself. ! ! This routine can be used to set up the sparse triplet storage ! for a linear triangle finite element discretization of Poisson's ! equation in two dimensions. ! ! Diagram: ! ! 3 ! s |\ ! i | \ ! d | \ ! e | \ side 2 ! | \ ! 3 | \ ! | \ ! 1-------2 ! ! side 1 ! ! The local node numbering ! ! ! 21-22-23-24-25 ! |\ |\ |\ |\ | ! | \| \| \| \| ! 16-17-18-19-20 ! |\ |\ |\ |\ | ! | \| \| \| \| ! 11-12-13-14-15 ! |\ |\ |\ |\ | ! | \| \| \| \| ! 6--7--8--9-10 ! |\ |\ |\ |\ | ! | \| \| \| \| ! 1--2--3--4--5 ! ! A sample grid ! ! ! Below, we have a chart that summarizes the adjacency relationships ! in the sample grid. On the left, we list the node, and its neighbors, ! with an asterisk to indicate the adjacency of the node to itself ! (in some cases, you want to count this self adjacency and in some ! you don't). On the right, we list the number of adjancencies to ! lower-indexed nodes, to the node itself, to higher-indexed nodes, ! the total number of adjacencies for this node, and the location ! of the first and last entries required to list this set of adjacencies ! in a single list of all the adjacencies. ! ! N Adjacencies Below Self Above Total First Last ! ! -- -- -- -- -- -- -- -- -- -- -- -- --- 0 ! 1: * 2 6 0 1 2 3 1 3 ! 2: 1 * 3 6 7 1 1 3 5 4 8 ! 3: 2 * 4 7 8 1 1 3 5 9 13 ! 4: 3 * 5 8 9 1 1 3 5 14 18 ! 5: 4 * 9 10 1 1 2 4 19 22 ! 6: 1 2 * 7 11 2 1 2 5 23 27 ! 7: 2 3 6 * 8 11 12 3 1 3 7 28 34 ! 8: 3 4 7 * 9 12 13 3 1 3 7 35 41 ! 9: 4 5 8 * 10 13 14 3 1 3 7 42 48 ! 10: 5 9 * 14 15 2 1 2 5 49 53 ! 11: 6 7 * 12 16 2 1 2 5 54 58 ! 12: 7 8 11 * 13 16 17 3 1 3 7 59 65 ! 13: 8 9 12 * 14 17 18 3 1 3 7 66 72 ! 14: 9 10 13 * 15 18 19 3 1 3 7 73 79 ! 15: 10 14 * 19 20 2 1 2 5 80 84 ! 16: 11 12 * 17 21 2 1 2 5 85 89 ! 17: 12 13 16 * 18 21 22 3 1 3 7 90 96 ! 18: 13 14 17 * 19 22 23 3 1 3 7 97 103 ! 19: 14 15 18 * 20 23 24 3 1 3 7 104 110 ! 20: 15 19 * 24 25 2 1 2 5 111 115 ! 21: 16 17 * 22 2 1 1 4 116 119 ! 22: 17 18 21 * 23 3 1 1 5 120 124 ! 23: 18 19 22 * 24 3 1 1 5 125 129 ! 24: 19 20 23 * 25 3 1 1 5 130 134 ! 25: 20 24 * 2 1 0 3 135 137 ! -- -- -- -- -- -- -- -- -- -- -- -- 138 --- ! ! For this example, the initial portion of the IA and JA arrays will be: ! ! (1,1), (1,2), (1,6), ! (2,1), (2,2), (2,3), (2,6), (2,7), ! (3,2), (3,3), (3,4), (3,7), (3,8), ! ... ! (25,20), (25,24), (25,25) ! ! for a total of 137 pairs of values. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 12 July 2007 ! ! Author: ! ! John Burkardt ! ! Parameters ! ! Input, integer ( kind = 4 ) NODE_NUM, the number of nodes. ! ! Input, integer ( kind = 4 ) TRIANGLE_NUM, the number of triangles. ! ! Input, integer ( kind = 4 ) TRIANGLE_NODE(3,TRIANGLE_NUM), lists the nodes ! that make up each triangle in counterclockwise order. ! ! Input, integer ( kind = 4 ) TRIANGLE_NEIGHBOR(3,TRIANGLE_NUM), for each ! side of a triangle, lists the neighboring triangle, or -1 if there is ! no neighbor. ! ! Input, integer ( kind = 4 ) ADJ_NUM, the number of adjacencies. ! ! Input, integer ( kind = 4 ) ADJ_COL(NODE_NUM+1). Information about ! column J is stored in entries ADJ_COL(J) through ADJ_COL(J+1)-1 of ADJ. ! ! Output, integer ( kind = 4 ) IA(ADJ_NUM), JA(ADJ_NUM), the adjacency ! information. ! implicit none integer ( kind = 4 ) adj_num integer ( kind = 4 ) node_num integer ( kind = 4 ) triangle_num integer ( kind = 4 ), parameter :: triangle_order = 3 integer ( kind = 4 ) adj_col(node_num+1) integer ( kind = 4 ) adj_copy(node_num) integer ( kind = 4 ) ia(adj_num) integer ( kind = 4 ) ja(adj_num) integer ( kind = 4 ) k1 integer ( kind = 4 ) k2 integer ( kind = 4 ) n1 integer ( kind = 4 ) n2 integer ( kind = 4 ) n3 integer ( kind = 4 ) node integer ( kind = 4 ) number integer ( kind = 4 ) triangle integer ( kind = 4 ) triangle2 integer ( kind = 4 ) triangle_neighbor(3,triangle_num) integer ( kind = 4 ) triangle_node(triangle_order,triangle_num) ia(1:adj_num) = -1 ja(1:adj_num) = -1 adj_copy(1:node_num) = adj_col(1:node_num) ! ! Set every node to be adjacent to itself. ! do node = 1, node_num ia(adj_copy(node)) = node ja(adj_copy(node)) = node adj_copy(node) = adj_copy(node) + 1 end do ! ! Examine each triangle. ! do triangle = 1, triangle_num n1 = triangle_node(1,triangle) n2 = triangle_node(2,triangle) n3 = triangle_node(3,triangle) ! ! Add edge (1,2) if this is the first occurrence, ! that is, if the edge (1,2) is on a boundary (TRIANGLE2 <= 0) ! or if this triangle is the first of the pair in which the edge ! occurs (TRIANGLE < TRIANGLE2). ! triangle2 = triangle_neighbor(1,triangle) if ( triangle2 < 0 .or. triangle < triangle2 ) then ia(adj_copy(n1)) = n1 ja(adj_copy(n1)) = n2 adj_copy(n1) = adj_copy(n1) + 1 ia(adj_copy(n2)) = n2 ja(adj_copy(n2)) = n1 adj_copy(n2) = adj_copy(n2) + 1 end if ! ! Add edge (2,3). ! triangle2 = triangle_neighbor(2,triangle) if ( triangle2 < 0 .or. triangle < triangle2 ) then ia(adj_copy(n2)) = n2 ja(adj_copy(n2)) = n3 adj_copy(n2) = adj_copy(n2) + 1 ia(adj_copy(n3)) = n3 ja(adj_copy(n3)) = n2 adj_copy(n3) = adj_copy(n3) + 1 end if ! ! Add edge (3,1). ! triangle2 = triangle_neighbor(3,triangle) if ( triangle2 < 0 .or. triangle < triangle2 ) then ia(adj_copy(n1)) = n1 ja(adj_copy(n1)) = n3 adj_copy(n1) = adj_copy(n1) + 1 ia(adj_copy(n3)) = n3 ja(adj_copy(n3)) = n1 adj_copy(n3) = adj_copy(n3) + 1 end if end do ! ! Lexically sort the IA, JA values. ! call i4vec2_sort_a ( adj_num, ia, ja ) return end subroutine triangulation_order3_boundary_node ( node_num, & element_num, element_node, node_boundary ) !*****************************************************************************80 ! !! TRIANGULATION_ORDER3_BOUNDARY_NODE indicates which nodes are on the boundary. ! ! Discussion: ! ! This routine is given a triangulation, an abstract list of sets of ! of nodes. It is assumed that the nodes in each triangle are listed ! in a counterclockwise order, although the routine should work ! if the nodes are consistently listed in a clockwise order as well. ! ! It is assumed that each edge of the triangulation is either ! * an INTERIOR edge, which is listed twice, once with positive ! orientation and once with negative orientation, or; ! * a BOUNDARY edge, which will occur only once. ! ! This routine should work even if the region has holes - as long ! as the boundary of the hole comprises more than 3 edges! ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 07 June 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) NODE_NUM, the number of nodes. ! ! Input, integer ( kind = 4 ) ELEMENT_NUM, the number of elements. ! ! Input, integer ( kind = 4 ) ELEMENT_NODE(3,ELEMENT_NUM), ! the nodes that make up the triangular elements. These should be listed ! in counterclockwise order. ! ! Output, logical NODE_BOUNDARY(NODE_NUM), is TRUE if the node ! is on a boundary edge. ! implicit none integer ( kind = 4 ) node_num integer ( kind = 4 ) element_num integer ( kind = 4 ) e1(3*element_num) integer ( kind = 4 ) e2(3*element_num) integer ( kind = 4 ) edge(2,3*element_num) integer ( kind = 4 ) i integer ( kind = 4 ) m integer ( kind = 4 ) n logical node_boundary(node_num) integer ( kind = 4 ) element_node(3,element_num) m = 2 n = 3 * element_num ! ! Set up the edge array. ! edge(1:2, 1: element_num) = element_node(1:2,1:element_num) edge(1:2, element_num+1:2*element_num) = element_node(2:3,1:element_num) edge(1, 2*element_num+1:3*element_num) = element_node(3, 1:element_num) edge(2, 2*element_num+1:3*element_num) = element_node(1, 1:element_num) ! ! In each column, force the smaller entry to appear first. ! e1(1:n) = minval ( edge(1:2,1:n), dim = 1 ) e2(1:n) = maxval ( edge(1:2,1:n), dim = 1 ) edge(1,1:n) = e1(1:n) edge(2,1:n) = e2(1:n) ! ! Ascending sort the column array. ! call i4col_sort_a ( m, n, edge ) ! ! Records which appear twice are internal edges and can be ignored. ! node_boundary(1:node_num) = .false. i = 0 do while ( i < 3 * element_num ) i = i + 1 if ( i == 3 * element_num ) then node_boundary(edge(1:m,i)) = .true. else if ( all ( edge(1:m,i) == edge(1:m,i+1) ) ) then i = i + 1 else node_boundary(edge(1:m,i)) = .true. end if end do return end subroutine triangulation_order3_neighbor_triangles ( triangle_num, & triangle_node, triangle_neighbor ) !*****************************************************************************80 ! !! TRIANGULATION_ORDER3_NEIGHBOR_TRIANGLES determines triangle neighbors. ! ! Discussion: ! ! A triangulation of a set of nodes can be completely described by ! the coordinates of the nodes, and the list of nodes that make up ! each triangle. However, in some cases, it is necessary to know ! triangle adjacency information, that is, which triangle, if any, ! is adjacent to a given triangle on a particular side. ! ! This routine creates a data structure recording this information. ! ! The primary amount of work occurs in sorting a list of 3 * TRIANGLE_NUM ! data items. ! ! Note that ROW is a work array allocated dynamically inside this ! routine. It is possible, for very large values of TRIANGLE_NUM, ! that the necessary amount of memory will not be accessible, and the ! routine will fail. This is a limitation of the implementation of ! dynamic arrays in FORTRAN90. One way to get around this would be ! to require the user to declare ROW in the calling routine ! as an allocatable array, get the necessary memory explicitly with ! an ALLOCATE statement, and then pass ROW into this routine. ! ! Of course, the point of dynamic arrays was to make it easy to ! hide these sorts of temporary work arrays from the poor user! ! ! This routine was revised to store the edge data in a column ! array rather than a row array. ! ! Example: ! ! The input information from TRIANGLE_NODE: ! ! Triangle Nodes ! -------- --------------- ! 1 3 4 1 ! 2 3 1 2 ! 3 3 2 8 ! 4 2 1 5 ! 5 8 2 13 ! 6 8 13 9 ! 7 3 8 9 ! 8 13 2 5 ! 9 9 13 7 ! 10 7 13 5 ! 11 6 7 5 ! 12 9 7 6 ! 13 10 9 6 ! 14 6 5 12 ! 15 11 6 12 ! 16 10 6 11 ! ! The output information in TRIANGLE_NEIGHBOR: ! ! Triangle Neighboring Triangles ! -------- --------------------- ! ! 1 -1 -1 2 ! 2 1 4 3 ! 3 2 5 7 ! 4 2 -1 8 ! 5 3 8 6 ! 6 5 9 7 ! 7 3 6 -1 ! 8 5 4 10 ! 9 6 10 12 ! 10 9 8 11 ! 11 12 10 14 ! 12 9 11 13 ! 13 -1 12 16 ! 14 11 -1 15 ! 15 16 14 -1 ! 16 13 15 -1 ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 14 February 2006 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) TRIANGLE_NUM, the number of triangles. ! ! Input, integer ( kind = 4 ) TRIANGLE_NODE(3,TRIANGLE_NUM), the nodes ! that make up each triangle. ! ! Output, integer ( kind = 4 ) TRIANGLE_NEIGHBOR(3,TRIANGLE_NUM), the three ! triangles that are direct neighbors of a given triangle. ! TRIANGLE_NEIGHBOR(1,I) is the index of the triangle which touches side 1, ! defined by nodes 2 and 3, and so on. TRIANGLE_NEIGHBOR(1,I) is negative ! if there is no neighbor on that side. In this case, that side of the ! triangle lies on the boundary of the triangulation. ! implicit none integer ( kind = 4 ) triangle_num integer ( kind = 4 ), parameter :: triangle_order = 3 integer ( kind = 4 ) col(4,3*triangle_num) integer ( kind = 4 ) i integer ( kind = 4 ) icol integer ( kind = 4 ) j integer ( kind = 4 ) k integer ( kind = 4 ) side1 integer ( kind = 4 ) side2 integer ( kind = 4 ) triangle_neighbor(3,triangle_num) integer ( kind = 4 ) tri integer ( kind = 4 ) triangle_node(triangle_order,triangle_num) integer ( kind = 4 ) tri1 integer ( kind = 4 ) tri2 ! ! Step 1. ! From the list of nodes for triangle T, of the form: (I,J,K) ! construct the three neighbor relations: ! ! (I,J,1,T) or (J,I,1,T), ! (J,K,2,T) or (K,J,2,T), ! (K,I,3,T) or (I,K,3,T) ! ! where we choose (I,J,1,T) if I < J, or else (J,I,1,T) ! do tri = 1, triangle_num i = triangle_node(1,tri) j = triangle_node(2,tri) k = triangle_node(3,tri) if ( i < j ) then col(1:4,3*(tri-1)+1) = (/ i, j, 1, tri /) else col(1:4,3*(tri-1)+1) = (/ j, i, 1, tri /) end if if ( j < k ) then col(1:4,3*(tri-1)+2) = (/ j, k, 2, tri /) else col(1:4,3*(tri-1)+2) = (/ k, j, 2, tri /) end if if ( k < i ) then col(1:4,3*(tri-1)+3) = (/ k, i, 3, tri /) else col(1:4,3*(tri-1)+3) = (/ i, k, 3, tri /) end if end do ! ! Step 2. Perform an ascending dictionary sort on the neighbor relations. ! We only intend to sort on rows 1 and 2; the routine we call here ! sorts on rows 1 through 4 but that won't hurt us. ! ! What we need is to find cases where two triangles share an edge. ! Say they share an edge defined by the nodes I and J. Then there are ! two columns of COL that start out ( I, J, ?, ? ). By sorting COL, ! we make sure that these two columns occur consecutively. That will ! make it easy to notice that the triangles are neighbors. ! call i4col_sort_a ( 4, 3*triangle_num, col ) ! ! Step 3. Neighboring triangles show up as consecutive columns with ! identical first two entries. Whenever you spot this happening, ! make the appropriate entries in TRIANGLE_NEIGHBOR. ! triangle_neighbor(1:3,1:triangle_num) = -1 icol = 1 do if ( 3 * triangle_num <= icol ) then exit end if if ( col(1,icol) /= col(1,icol+1) .or. col(2,icol) /= col(2,icol+1) ) then icol = icol + 1 cycle end if side1 = col(3,icol) tri1 = col(4,icol) side2 = col(3,icol+1) tri2 = col(4,icol+1) triangle_neighbor(side1,tri1) = tri2 triangle_neighbor(side2,tri2) = tri1 icol = icol + 2 end do return end