program main !*****************************************************************************80 ! !! vandermonde_approx_2d_test() tests vandermonde_approx_2d(). ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 23 September 2012 ! ! Author: ! ! John Burkardt ! implicit none integer, parameter :: m_test_num = 5 integer j integer m integer, dimension(m_test_num) :: m_test = (/ 0, 1, 2, 4, 8 /) integer grid integer prob integer prob_num call timestamp ( ) write ( *, '(a)' ) '' write ( *, '(a)' ) 'VANDERMONDE_APPROX_2D_TEST():' write ( *, '(a)' ) ' FORTRAN90 version' write ( *, '(a)' ) ' Test VANDERMONDE_APPROX_2D().' write ( *, '(a)' ) ' The QR_SOLVE library is needed.' write ( *, '(a)' ) ' The R8LIB library is needed.' write ( *, '(a)' ) ' This test also needs the TEST_INTERP_2D library.' call f00_num ( prob_num ) do prob = 1, prob_num grid = 1 do j = 1, m_test_num m = m_test(j) call test01 ( prob, grid, m ) end do end do ! ! Terminate. ! write ( *, '(a)' ) '' write ( *, '(a)' ) 'VANDERMONDE_APPROX_2D_TEST():' write ( *, '(a)' ) ' Normal end of execution.' write ( *, '(a)' ) '' call timestamp ( ) stop 0 end subroutine test01 ( prob, grd, m ) !*****************************************************************************80 ! !! VANDERMONDE_APPROX_2D_TEST01 tests VANDERMONDE_APPROX_2D_MATRIX. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 23 September 2012 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer PROB, the problem number. ! ! Input, integer GRD, the grid number. ! (Can't use GRID as the name because that's also a plotting function.) ! ! Input, integer M, the total polynomial degree. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) real ( kind = rk ), allocatable :: a(:,:) real ( kind = rk ) app_error real ( kind = rk ), allocatable :: c(:) integer grd integer m integer nd integer ni integer prob real ( kind = rk ) r8vec_norm_affine integer tm integer triangle_num real ( kind = rk ), allocatable :: xd(:) real ( kind = rk ), allocatable :: xi(:) real ( kind = rk ), allocatable :: yd(:) real ( kind = rk ), allocatable :: yi(:) real ( kind = rk ), allocatable :: zd(:) real ( kind = rk ), allocatable :: zi(:) write ( *, '(a)' ) '' write ( *, '(a)' ) 'TEST01:' write ( *, '(a,i4)' ) ' Approximate data from TEST_INTERP_2D problem #', prob write ( *, '(a,i4)' ) ' Use grid from TEST_INTERP_2D with index #', grd write ( *, '(a,i4)' ) ' Using polynomial approximant of total degree ', m call g00_size ( grd, nd ) write ( *, '(a,i6)' ) ' Number of data points = ', nd allocate ( xd(1:nd) ) allocate ( yd(1:nd) ) call g00_xy ( grd, nd, xd, yd ) allocate ( zd(1:nd) ) call f00_f0 ( prob, nd, xd, yd, zd ) if ( nd < 10 ) then call r8vec3_print ( nd, xd, yd, zd, ' X, Y, Z data:' ) end if ! ! Compute the Vandermonde matrix. ! tm = triangle_num ( m + 1 ); allocate ( a(1:nd,1:tm) ) call vandermonde_approx_2d_matrix ( nd, m, tm, xd, yd, a ) ! ! Solve linear system. ! allocate ( c(1:tm) ) call qr_solve ( nd, tm, a, zd, c ) ! ! #1: Does approximant match function at data points? ! ni = nd allocate ( xi(1:ni) ) allocate ( yi(1:ni) ) xi(1:ni) = xd(1:ni) yi(1:ni) = yd(1:ni) allocate ( zi(1:ni) ) call r8poly_value_2d ( m, c, ni, xi, yi, zi ) app_error = r8vec_norm_affine ( ni, zi, zd ) / real ( ni, kind = rk ) write ( *, '(a)' ) '' write ( *, '(a,g14.6)' ) ' L2 data approximation error = ', app_error deallocate ( a ) deallocate ( c ) deallocate ( xd ) deallocate ( xi ) deallocate ( yd ) deallocate ( yi ) deallocate ( zd ) deallocate ( zi ) return end subroutine r8mat_cholesky_factor ( n, a, c, flag ) !*****************************************************************************80 ! !! r8mat_cholesky_factor() computes the Cholesky factor of a symmetric matrix. ! ! Discussion: ! ! The matrix must be symmetric and positive semidefinite. ! ! For a positive semidefinite symmetric matrix A, the Cholesky factorization ! is a lower triangular matrix L such that: ! ! A = L * L' ! ! An R8MAT is an MxN array of R8's, stored by (I,J) -> [I+J*M]. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 08 April 2009 ! ! Author: ! ! John Burkardt ! ! Input: ! ! integer N, the number of rows and columns. ! ! real ( kind = rk8 ) A(N,N), the N by N matrix. ! ! Output: ! ! real ( kind = rk8 ) C(N,N), the lower triangular Cholesky factor. ! ! integer FLAG: ! 0, no error occurred. ! 1, the matrix is not positive definite. ! implicit none integer, parameter :: rk8 = kind ( 1.0D+00 ) integer n real ( kind = rk8 ) a(n,n) real ( kind = rk8 ) c(n,n) integer flag integer i integer j real ( kind = rk8 ) sum2 flag = 0 c(1:n,1:n) = a(1:n,1:n) do j = 1, n c(1:j-1,j) = 0.0D+00 do i = j, n sum2 = c(j,i) - dot_product ( c(j,1:j-1), c(i,1:j-1) ) if ( i == j ) then if ( sum2 <= 0.0D+00 ) then flag = 1 return else c(i,j) = sqrt ( sum2 ) end if else if ( c(j,j) /= 0.0D+00 ) then c(i,j) = sum2 / c(j,j) else c(i,j) = 0.0D+00 end if end if end do end do return end subroutine r8mat_cholesky_solve ( n, l, b, x ) !*****************************************************************************80 ! !! r8mat_cholesky_solve() solves a Cholesky factored linear system A * x = b. ! ! Discussion: ! ! An R8MAT is an MxN array of R8's, stored by (I,J) -> [I+J*M]. ! ! This routine works with the lower triangular Cholesky factor A = L * L'. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 07 December 2004 ! ! Author: ! ! John Burkardt ! ! Input: ! ! integer N, the number of rows and columns. ! ! real ( kind = rk8 ) L(N,N), the N by N lower Cholesky factor of the ! system matrix A. ! ! real ( kind = rk8 ) B(N), the right hand side of the linear system. ! ! Output: ! ! real ( kind = rk8 ) X(N), the solution of the linear system. ! implicit none integer, parameter :: rk8 = kind ( 1.0D+00 ) integer n real ( kind = rk8 ) b(n) real ( kind = rk8 ) l(n,n) real ( kind = rk8 ) x(n) ! ! Solve L * y = b. ! call r8mat_l_solve ( n, l, b, x ) ! ! Solve L' * x = y. ! call r8mat_lt_solve ( n, l, x, x ) return end subroutine r8mat_l_solve ( n, a, b, x ) !*****************************************************************************80 ! !! r8mat_l_solve() solves a lower triangular linear system. ! ! Discussion: ! ! An R8MAT is an MxN array of R8's, stored by (I,J) -> [I+J*M]. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 07 December 2004 ! ! Author: ! ! John Burkardt ! ! Input: ! ! integer N, the number of rows and columns. ! ! real ( kind = rk8 ) A(N,N), the N by N lower triangular matrix. ! ! real ( kind = rk8 ) B(N), the right hand side of the linear system. ! ! Output: ! ! real ( kind = rk8 ) X(N), the solution of the linear system. ! implicit none integer, parameter :: rk8 = kind ( 1.0D+00 ) integer n real ( kind = rk8 ) a(n,n) real ( kind = rk8 ) b(n) integer i real ( kind = rk8 ) x(n) ! ! Solve L * x = b. ! do i = 1, n x(i) = ( b(i) - dot_product ( a(i,1:i-1), x(1:i-1) ) ) / a(i,i) end do return end subroutine r8mat_lt_solve ( n, a, b, x ) !*****************************************************************************80 ! !! r8mat_lt_solve() solves a transposed lower triangular linear system. ! ! Discussion: ! ! An R8MAT is an MxN array of R8's, stored by (I,J) -> [I+J*M]. ! ! Given the lower triangular matrix A, the linear system to be solved is: ! ! A' * x = b ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 07 December 2004 ! ! Author: ! ! John Burkardt ! ! Input: ! ! integer N, the number of rows and columns. ! ! real ( kind = rk8 ) A(N,N), the N by N lower triangular matrix. ! ! real ( kind = rk8 ) B(N), the right hand side of the linear system. ! ! Output: ! ! real ( kind = rk8 ) X(N), the solution of the linear system. ! implicit none integer, parameter :: rk8 = kind ( 1.0D+00 ) integer n real ( kind = rk8 ) a(n,n) real ( kind = rk8 ) b(n) integer i real ( kind = rk8 ) x(n) ! ! Solve L'*x = b. ! do i = n, 1, -1 x(i) = ( b(i) - dot_product ( x(i+1:n), a(i+1:n,i) ) ) / a(i,i) end do return end subroutine r8poly_value_2d ( m, c, n, x, y, p ) !*****************************************************************************80 ! !! R8POLY_VALUE_2D evaluates a polynomial in 2 variables, X and Y. ! ! Discussion: ! ! We assume the polynomial is of total degree M, and has the form: ! ! p(x,y) = c00 ! + c10 * x + c01 * y ! + c20 * x^2 + c11 * xy + c02 * y^2 ! + ... ! + cm0 * x^(m) + ... + c0m * y^m. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 31 August 2012 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer M, the degree of the polynomial. ! ! Input, real ( kind = rk ) C(T(M+1)), the polynomial coefficients. ! C(1) is the constant term. T(M+1) is the M+1-th triangular number. ! The coefficients are stored consistent with the following ordering ! of monomials: 1, X, Y, X^2, XY, Y^2, X^3, X^2Y, XY^2, Y^3, X^4, ... ! ! Input, integer N, the number of evaluation points. ! ! Input, real ( kind = rk ) X(N), Y(N), the evaluation points. ! ! Output, real ( kind = rk ) P(N), the value of the polynomial at the ! evaluation points. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer n real ( kind = rk ) c(*) integer ex integer ey integer j integer m real ( kind = rk ) p(n) integer s real ( kind = rk ) x(n) real ( kind = rk ) y(n) p(1:n) = 0.0D+00 j = 0 do s = 0, m do ex = s, 0, -1 ey = s - ex j = j + 1 p(1:n) = p(1:n) + c(j) * x(1:n) ** ex * y(1:n) ** ey end do end do return end function r8vec_norm_affine ( n, v0, v1 ) !*****************************************************************************80 ! !! r8vec_norm_affine() returns the affine norm of an R8VEC. ! ! Discussion: ! ! An R8VEC is a vector of R8's. ! ! The affine vector L2 norm is defined as: ! ! R8VEC_NORM_AFFINE(V0,V1) ! = sqrt ( sum ( 1 <= I <= N ) ( V1(I) - V0(I) )^2 ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 27 October 2010 ! ! Author: ! ! John Burkardt ! ! Input: ! ! integer N, the order of the vectors. ! ! real ( kind = rk8 ) V0(N), the base vector. ! ! real ( kind = rk8 ) V1(N), the vector whose affine norm is desired. ! ! Output: ! ! real ( kind = rk8 ) R8VEC_NORM_AFFINE, the L2 norm of V1-V0. ! implicit none integer, parameter :: rk8 = kind ( 1.0D+00 ) integer n real ( kind = rk8 ) r8vec_norm_affine real ( kind = rk8 ) v0(n) real ( kind = rk8 ) v1(n) r8vec_norm_affine = sqrt ( sum ( ( v0(1:n) - v1(1:n) ) ** 2 ) ) return end subroutine r8vec3_print ( n, a1, a2, a3, title ) !*****************************************************************************80 ! !! r8vec3_print() prints an R8VEC3. ! ! Discussion: ! ! An R8VEC3 is a dataset consisting of N triples of R8's, stored ! as three separate vectors A1, A2, A3. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 22 September 2012 ! ! Author: ! ! John Burkardt ! ! Input: ! ! integer N, the number of components of the vector. ! ! real ( kind = rk8 ) A1(N), A2(N), A3(N), the vectors to be printed. ! ! character ( len = * ) TITLE, a title. ! implicit none integer, parameter :: rk8 = kind ( 1.0D+00 ) integer n real ( kind = rk8 ) a1(n) real ( kind = rk8 ) a2(n) real ( kind = rk8 ) a3(n) integer i character ( len = * ) title write ( *, '(a)' ) ' ' write ( *, '(a)' ) trim ( title ) write ( *, '(a)' ) ' ' do i = 1, n write ( *, '(i8,3g14.6)' ) i, a1(i), a2(i), a3(i) end do return end subroutine timestamp ( ) !*****************************************************************************80 ! !! timestamp() prints the current YMDHMS date as a time stamp. ! ! Example: ! ! 31 May 2001 9:45:54.872 AM ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 15 August 2021 ! ! Author: ! ! John Burkardt ! implicit none character ( len = 8 ) ampm integer d integer h integer m integer mm character ( len = 9 ), parameter, dimension(12) :: month = (/ & 'January ', 'February ', 'March ', 'April ', & 'May ', 'June ', 'July ', 'August ', & 'September', 'October ', 'November ', 'December ' /) integer n integer s integer values(8) integer y call date_and_time ( values = 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 ( *, '(i2.2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) & d, trim ( month(m) ), y, h, ':', n, ':', s, '.', mm, trim ( ampm ) return end