PROGRAM hilbert ! Demonstration program 1 using the 'quad' package. ! Inversion of 10x10 Hilbert matrix by Cholesky factorization. ! Matrix A has elements a(i,j) = 1/(i+j-1). ! To obtain full accuracy, the elements are all multiplied by ! 9 x 5 x 7 x 11 x 13 x 17 x 19 to make them into integers so ! that the input matrix is exact. The elements of the inverse ! will be multiplied by this integer at the end. ! The elements of the inverse should all be integers. ! Programmer: Alan Miller ! Latest Fortran 77 revision - 13 September 1986 ! Fortran 90 version - 3 December 1996 ! Latest revision - 3 October 1999 USE quadruple_precision IMPLICIT NONE TYPE (quad) :: a(55), total, exact_start, exact, relerr REAL (dp) :: det, scale_factor = 14549535.d0 INTEGER :: i, j, nrows = 10, ier, ipos, ik, k, kj, kjstart INTERFACE SUBROUTINE qchol(a, nrows, ier) USE quadruple_precision IMPLICIT NONE TYPE (quad), DIMENSION(:), INTENT(IN OUT) :: a INTEGER, INTENT(IN) :: nrows INTEGER, INTENT(OUT) :: ier END SUBROUTINE qchol END INTERFACE ! Enter the elements of the Hilbert matrix, storing them in a ! 1-dimensional array up to the diagonal in each row. ipos = 1 DO i = 1, nrows DO j = 1, i a(ipos) = quad(scale_factor/(i+j-1), 0._dp) ipos = ipos + 1 END DO END DO ! Get the Cholesky factorization. CALL qchol(a, nrows, ier) IF (ier /= 0) THEN WRITE(*, *) 'Error in qchol, number:', ier STOP END IF ! Calculate the determinant by multiplying together the diagonal elements ! of the Cholesky factor, squaring the result, and remembering that each ! element was multiplied by scale_factor. ipos = 1 det = 1.d0 DO i = 1, nrows det = det * a(ipos)%hi ipos = ipos + i + 1 END DO det = det * det / scale_factor**nrows WRITE(*, 900) det 900 FORMAT(' Determinant = ', g20.12/) ! Now invert the triangular matrix. ipos = 1 DO i = 1, nrows kjstart = 1 DO j = 1, i kj = kjstart ik = ipos IF (j < i) THEN total = quad(0._dp, 0._dp) DO k = j, i-1 total = total - a(ik) * a(kj) ik = ik + 1 kj = kj + k END DO a(ipos) = total / a(ik) ipos = ipos + 1 ELSE a(ipos) = quad(1._dp, 0._dp) / a(ipos) ipos = ipos + 1 END IF kjstart = kjstart + j + 1 END DO END DO ! We formed A = LL', so the inverse is L**(-T).L**(-1), ! where L**(-T) denotes the transpose of the inverse. ipos = 1 kjstart = 1 DO i = 1, nrows DO j = 1, i ik = ipos total = quad(0._dp, 0._dp) kj = kjstart DO k = i, nrows total = total + a(ik) * a(kj) ik = ik + k kj = kj + k END DO a(ipos) = total * scale_factor ipos = ipos + 1 END DO kjstart = kjstart + i + 1 END DO ! Output the result. ! The exact values for an n x n matrix are:- ! (-1)**(i+j).(n+i-1)!(n+j-1)! ! ------------------------------------- ! (i+j-1)[(i-1)!(j-1)!]**2.(n-i)!(n-j)! WRITE(*, 910) 910 FORMAT(' row col calculated exact Rel. error') ipos = 1 exact_start = quad(100._dp, 0._dp) DO i = 1, nrows exact = exact_start DO j = 1, i relerr = (a(ipos) - exact) / exact WRITE(*, 920) i, j, a(ipos)%hi, exact%hi, relerr%hi 920 FORMAT(2I4, 2F20.4, g14.5) exact = -exact * (10+j)*(10-j)*(i+j-1) / REAL((i+j)*j*j) ipos = ipos + 1 END DO exact_start = -exact_start * (10+i)*(10-i) / REAL((i+1)*i) END DO STOP END PROGRAM hilbert SUBROUTINE qchol(a, nrows, ier) ! Quadruple-precision Cholesky factorization of a +ve definite symmetric ! matrix, stored by rows up to the diagonal. The input matrix is ! overwritten by the lower-triangular factor L, where A = LL'. ! Arguments: ! a(:) INPUT. Array containing the lower triangle of matrix A. ! nrows INPUT. The number of rows (and columns) in A. ! ier OUTPUT. Error indicator ! = 0 if no error detected ! = 1 if nrows < 1 ! = 2 if matrix not +ve definite ! Programmer: Alan Miller ! Latest revision - 12 September 1986 ! Fortran 90 version - 3 December 1996 USE quadruple_precision IMPLICIT NONE TYPE (quad), DIMENSION(:), INTENT(IN OUT) :: a INTEGER, INTENT(IN) :: nrows INTEGER, INTENT(OUT) :: ier ! Local variables INTEGER :: i, j, ij, ikstart, jk, ik, k TYPE (quad) :: total IF (nrows <= 0) THEN ier = 1 RETURN END IF ij = 1 ikstart = 1 DO i = 1, nrows jk = 1 DO j = 1, i ik = ikstart total = a(ij) ! Form the sum a(i,j) - l(i,k).l(j,k) DO k = 1, j-1 total = total - a(ik) * a(jk) ik = ik + 1 jk = jk + 1 END DO ! On the diagonal l(i,i) = sqrt(total) ! Off diagonal l(i,j) = total/l(j,j) IF (i == j) THEN IF (total%hi <= 0._dp) THEN ier = 2 RETURN END IF a(ij) = SQRT(total) ELSE a(ij) = total / a(jk) jk = jk + 1 END IF ij = ij + 1 END DO ! ikstart = location of first element in row i. ikstart = ij END DO ier = 0 RETURN END SUBROUTINE qchol