! Last change: DOS 29 Jul 2000 12:15 pm ! *** copyright 2000 *** ! *** filename runge1.f95 *** John F. Monahan ** ! ********************** program runge1 ! runge1 -- first demonstration of runge phenomenon ! approximate a smooth function with interpolating polynomial ! from equally spaced abscissas and watch it go weird ! implicit none real, dimension(10,10) :: a real, dimension(10) :: b real, dimension(256,3) :: y real xi,rcond,s integer, parameter :: npts = 100 integer, dimension(10) :: pi integer i,j,n,kn ! INTERFACE BLOCK INTERFACE subroutine gaucpp(a,n,pi,rcond) IMPLICIT NONE REAL, DIMENSION(:,:), INTENT(IN OUT) :: a REAL, INTENT(OUT) :: rcond INTEGER, INTENT (IN) :: n INTEGER, DIMENSION(:), INTENT(OUT) :: pi end subroutine gaucpp SUBROUTINE GAUSSE(A,N,PI,B) IMPLICIT NONE REAL, DIMENSION(:,:), INTENT(IN) :: A REAL, DIMENSION(:), INTENT(IN OUT) :: B INTEGER, INTENT (IN) :: N INTEGER, DIMENSION(:), INTENT(IN) :: PI END SUBROUTINE GAUSSE END INTERFACE ! 21 format(' n = ',i2,' rcond =',f9.6) 24 format(f9.6,2x,4f10.6) ! open( unit=6, file='runge1.out') open( unit=9, file='runge1.dat') ! do n = 5,9,2 ! three numbers of points kn = (n-3)/2 ! kn = 1, 2, 3 ! ! set up vandermonde matrix for interpolation do i = 1,n ! x's equally spaced on (-4,4) xi = 4.*real(2*i-1-n)/real(n-1) ! function interpolated is simple-looking b(i) = 1./(1.+xi*xi) ! make a copy a(i,1) = 1. do j = 2,n a(i,j) = xi*a(i,j-1) end do ! loop on j end do ! loop on i ! have matrix -- factor and get condition estimate call gaucpp(a,n,pi,rcond) write(6,21) n,rcond ! solve equations call gausse(a,n,pi,b) ! now have polynomial coefficients, ! compute points of interpolant do i = 1,npts xi = 4.1*real(2*i-1-npts)/real(npts-1) s = b(n) do j = 2,n s = xi*s + b(n-j+1) end do ! loop on j y(i,kn) = s end do ! loop on i end do ! loop on n ! write out results for 3 interpolants do i = 1,npts xi = 4.1*real(2*i-1-npts)/real(npts-1) s = 1./(1.+xi*xi) write(9,24) xi,s,(y(i,j),j=1,3) end do ! loop on i stop end program runge1 subroutine gaucpp(a,n,pi,rcond) ! computes the LU factorization of the matrix a using gaussian ! elimination with partial pivoting and computes condtion estimate IMPLICIT NONE REAL, DIMENSION(:,:), INTENT(IN OUT) :: a REAL, INTENT(OUT) :: rcond INTEGER, INTENT (IN) :: n INTEGER, DIMENSION(:), INTENT(OUT) :: pi ! REAL, DIMENSION(n) :: p,y real anorm,rnorm,znorm,ykp,ykm,sp,sm integer i,k,kk,kp1 ! ! a real input: matrix of real dimension n by n ! output: LU factorization, PA = LU ! n integer dimension of matrix ! pi integer vector of permutations -- forms P ! rcond real reciprocal of condition number of a using ! infinity norm ! ! first step is to compute the infinity norm of the matrix A ! then the LU factorization is computed, overwriting A, using GAUSPP ! then compute the vector y which leads to large norm estimate ! using Golub and van Loan Algorithm 4.5-1 ! then solve Az=y to finish using GAUSSE ! ! See Golub and van Loan, Matrix Computations, pp.76-78 ! Cline, Moler, Stewart, and Wilkinson (1979) An Estimate of the ! Condition Number of a Matrix, SIAM J Num Anal, 16, p368ff ! Dongarra, Bunch, Moler and Stewart (1979) Linpack User's Guide ! ! J F Monahan (June, 1992) Dept of Stat, NCSU, Raleigh, NC 27695-8203 ! revised June 1999, April 2000 for Fortran 95 ! ! INTERFACE BLOCK INTERFACE SUBROUTINE GAUSPP(A,N,PI) IMPLICIT NONE REAL, DIMENSION(:,:), INTENT(IN OUT) :: A INTEGER, INTENT (IN) :: N INTEGER, DIMENSION(:), INTENT(OUT) :: PI END SUBROUTINE GAUSPP SUBROUTINE GAUSSE(A,N,PI,B) IMPLICIT NONE REAL, DIMENSION(:,:), INTENT(IN) :: A REAL, DIMENSION(:), INTENT(IN OUT) :: B INTEGER, INTENT (IN) :: N INTEGER, DIMENSION(:), INTENT(IN) :: PI END SUBROUTINE GAUSSE END INTERFACE rcond = 0. ! compute the norm of A anorm = 0. do k = 1,n sm = 0. do i = 1,n sm = sm + abs( a(k,i) ) enddo ! loop on i anorm = max( anorm, sm ) enddo ! loop on k ! call gauspp to compute the factorization P A = L U call gauspp(a,n,pi) ! if singular, then return as is with rcond = 0. if( pi(n) .ne. n ) return ! have L and U factors stored in A, permutation in PI ! apply algorithm 4.5-1 to transpose(U) do i = 1,n p(i) = 0. enddo ! loop on i ! now the main loop rnorm = 1. do k = 1,n ykp = (1./rnorm - p(k) ) / a(k,k) ykm = - (1./rnorm + p(k) ) / a(k,k) sp = abs( ykp ) sm = abs( ykm ) kp1 = k + 1 if( kp1 .le. n ) then do i = kp1,n sp = sp + abs( (p(i) + a(i,k)*ykp)/a(i,i) ) sm = sm + abs( (p(i) + a(i,k)*ykm)/a(i,i) ) enddo ! loop on i endif ! have the best y(k) y(k) = ykp if( sp .ge. sm ) y(k) = ykm if( kp1 .le. n ) then do i = kp1,n p(i) = p(i) + a(i,k)*y(k) enddo ! loop on i endif ! normalize every time -- these grow fast! sm = abs( y(k) ) if( sm .gt. rnorm ) then rnorm = sm do i = 1,n y(i) = y(i) / rnorm p(i) = p(i) / rnorm enddo ! loop on i endif ! end of main loop enddo ! loop on k ! now solve transpose(L) r = y but overwrite do kk = 1,n k = n + 1 - kk if( k .ne. n ) then kp1 = k + 1 do i = kp1,n y(k) = y(k) - a(i,k)*y(i) enddo ! loop on i endif enddo ! loop on kk, k ! compute norm of r rnorm = 0. do i = 1,n rnorm = max( rnorm, abs( y(i) ) ) enddo ! loop on i ! solve a system in A call gausse(a,n,pi,y) ! compute norm of this solution (z) znorm = 0. do i = 1,n znorm = max( znorm, abs( y(i) ) ) enddo ! compute the reciprocal of condition number rcond = rnorm / (anorm*znorm) return end subroutine gaucpp SUBROUTINE GAUSPP(A,N,PI) ! SUBROUTINE TO DO THE LU FACTORIZATION OF GAUSSIAN ELIMINATION ! FACTORS P * A = L U WHERE P IS THE PERMUTATION MATRIX STORED IN PI ! A IS DESTROYED AS L AND U ARE OVERWRITTEN ON IT ! L IS UNIT LOWER AND U UPPER TRIANGULAR ! ! USE GAUSSE TO SOLVE A SYSTEM OF LINEAR EQUATIONS ! ! A IN MATRIX OF COEFFICIENTS N BY N ! OUT FACTORS L AND U OVERLAID ! N IN TRUE SIZE OF MATRIX A (NOT DECLARED SIZE) ! OUT UNCHANGED ! PI IN ANYTHING ! OUT PERMUTATION LIST OF N-1 PIVOT LOCATIONS ! PI(N) GIVES COMPUTED RANK, PI(N) = N IF NONSINGULAR ! ! J F MONAHAN, SEPT 1985, N C S U, RALEIGH, NC 27695-8203 USA ! REVISED JUNE 1999, APRIL 2000 FOR FORTRAN 95 ! IMPLICIT NONE REAL, DIMENSION(:,:), INTENT(IN OUT) :: A INTEGER, INTENT (IN) :: N INTEGER, DIMENSION(:), INTENT(OUT) :: PI REAL S,BIG INTEGER I,J,K,KP1,NM1,IBIG ! PI(N) = 0 IF( N .GT. 1 ) THEN ! DO NM1 STEPS TO FORM A TRIANGULAR SYSTEM NM1 = N - 1 ! DO K = 1,NM1 ! SEARCH FOR BEST PIVOT IN COLUMN K BIG = 0. IBIG = 0 DO I = K,N S = ABS( A(I,K) ) IF( S .GT. BIG ) THEN ! KEEP TRACK OF BEST I AND ITS SIZE BIG = S IBIG = I ENDIF ENDDO ! LOOP ON I ! TEST FOR SINGULARITY ! CHANGE 0. TO AN EPSILON IF DESIRED IF( BIG .EQ. 0. ) RETURN ! KEEP TRACK OF THE RANK IN PI(N) PI(N) = K ! STORE PIVOT LOCATION IN PERMUTATION LIST PI(K) = IBIG ! SWITCH ONLY IF NECESSARY IF( IBIG .NE. I ) THEN DO J = 1,N S = A(K,J) A(K,J) = A(IBIG,J) A(IBIG,J) = S ENDDO ! LOOP ON J ! ENDIF ! ELIMINATE REMAINDER OF COLUMN K AND CREATE ONE KP1 = K + 1 DO I = KP1,N A(I,K) = A(I,K) / A(K,K) DO J = KP1,N A(I,J) = A(I,J) - A(I,K) * A(K,J) ENDDO ! LOOP ON J ENDDO ! LOOP ON I ENDDO ! END OF MAIN LOOP ON K ENDIF ! TEST LAST (AND MAYBE ONLY) ONE FOR SINGULARITY IF ( ABS( A(N,N) ) .GT. 0. ) THEN ! ONLY (COMPUTATIONALLY) NONSINGULAR MATRICES ! REACH HERE PI(N) = N ENDIF RETURN END SUBROUTINE GAUSPP SUBROUTINE GAUSSE(A,N,PI,B) ! ! SUBROUTINE TO SOLVE A SYSTEM OF LINEAR EQUATIONS A * X = B ! ! A IN HOLDS FACTORS L AND U OVERLAID FROM GAUSPP ! OUT UNCHANGED ! N IN TRUE SIZE OF MATRIX A (NOT DECLARED SIZE) ! OUT UNCHANGED ! PI IN LIST OF PIVOTS -- HOLDS PERMUTATION MATRIX P ! IN FACTORED FORM ! OUT UNCHANGED ! ! *** GAUSPP MUST HAVE BEEN CALLED PREVIOUSLY *** ! ! GAUSPP FACTORS P * A = L U ! TO SOLVE A X = B, FIRST PERMUTE B TO GET P A X = L U X = P B ! THEN SOLVE TWO TRIANGULAR SYSTEMS OF EQUATIONS ! ! J F MONAHAN, SEPT 1985, N C S U, RALEIGH, NC 27695 - 8203 USA ! REVISED JUNE, 1999 FOR FORTRAN 95 ! ! DECLARATIONS FIRST IMPLICIT NONE REAL, DIMENSION(:,:), INTENT(IN OUT) :: A REAL, DIMENSION(:), INTENT(IN OUT) :: B INTEGER, INTENT (IN) :: N INTEGER, DIMENSION(:), INTENT(IN) :: PI REAL S INTEGER I,J,K,II,IP1,IM1,NM1 ! IF( N .GT. 1 ) THEN NM1 = N - 1 ! PERMUTE B FIRST DO K = 1,NM1 S = B(K) I = PI(K) B(K) = B(I) B(I) = S ENDDO ! LOOP ON K ! L IS UNIT TRIANGULAR MATRIX ! SOLVE L * Y = P * B FROM THE TOP DOWN DO I = 2,N IM1 = I - 1 DO J = 1,IM1 B(I) = B(I) - A(I,J) * B(J) ENDDO ! LOOP ON J ENDDO ! LOOP ON I ! U IS UPPER TRIANGULAR ! SOLVE U * X = Y ENDIF B(N) = B(N) / A(N,N) ! LAST ONE DONE FIRST IF( N .EQ. 1 ) RETURN ! NOW FROM THE BOTTOM UP, I = N-1, ... , 2, 1 DO II = 2,N I = N + 1 - II IP1 = I + 1 DO J = IP1,N B(I) = B(I) - A(I,J) * B(J) ENDDO ! LOOP ON J B(I) = B(I) / A(I,I) ENDDO ! LOOP ON II, I RETURN END SUBROUTINE GAUSSE ! *** end of filename runge1.f95 *********************