SUBROUTINE syminv(a, n, nn, c, nullty, ifault) ! Algorithm AS7, Applied Statistics, vol.17, 1968, p.198. ! N.B. Argument W has been removed. ! Forms in c( ) as lower triangle, a generalised inverse ! of the positive semi-definite symmetric matrix a( ) ! order n, stored as lower triangle. ! arguments:- ! a() = input, the symmetric matrix to be inverted, stored in ! lower triangular form ! n = input, order of the matrix ! nn = input, the size of the a and c arrays >= n*(n+1)/2 ! c() = output, the inverse of a (a generalized inverse if c is ! singular), also stored in lower triangular. ! c and a may occupy the same locations. ! nullty = output, the rank deficiency of a. ! ifault = output, error indicator ! = 1 if n < 1 ! = 2 if a is not +ve semi-definite ! = 3 if nn < n*(n+1)/2 ! = 0 otherwise !*************************************************************************** IMPLICIT NONE INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(15, 60) REAL (dp), INTENT(IN) :: a(:) INTEGER, INTENT(IN) :: n, nn REAL (dp), INTENT(OUT) :: c(:) INTEGER, INTENT(OUT) :: nullty, ifault INTERFACE SUBROUTINE chol (a, n, nn, u, nullty, ifault) IMPLICIT NONE INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(15, 60) REAL (dp), INTENT(IN) :: a(:) INTEGER, INTENT(IN) :: n, nn REAL (dp), INTENT(OUT) :: u(:) INTEGER, INTENT(OUT) :: nullty, ifault END SUBROUTINE chol END INTERFACE ! Local variables REAL (dp) :: w(n), x REAL (dp), PARAMETER :: zero = 0.0_dp, one = 1.0_dp INTEGER :: i, icol, irow, j, jcol, k, l, mdiag, ndiag ! Cholesky factorization of a, result in c CALL chol(a, n, nn, c, nullty, ifault) IF (ifault /= 0) RETURN ! Invert c & form the product (cinv)'*cinv, where cinv is the inverse ! of c, row by row starting with the last row. ! irow = the row number, ndiag = location of last element in the row. ndiag = nn DO irow = n, 1, -1 l = ndiag IF (c(ndiag) == zero) THEN DO j = irow, n c(l) = zero l = l + j END DO GO TO 80 END IF DO i = irow, n w(i) = c(l) l = l + i END DO icol = n jcol = nn mdiag = nn DO l = jcol x = zero IF(icol == irow) x = one / w(irow) k = n DO IF(k == irow) EXIT x = x - w(k)*c(l) k = k - 1 l = l - 1 IF(l > mdiag) l = l - k + 1 END DO c(l) = x/w(irow) IF(icol == irow) EXIT mdiag = mdiag - icol icol = icol - 1 jcol = jcol - 1 END DO 80 ndiag = ndiag - irow END DO RETURN END SUBROUTINE syminv SUBROUTINE subinv(a, nm, b, n, c, nullty, ifault, det) ! Remark asr 44.1 Appl. Statist. (1982) Vol.31, No.3 ! A revised and enhanced version of ! algorithm as7, applied statistics, vol.17, 1968. ! N.B. Arguments W & NDIM have been removed, and DET is OPTIONAL. ! Forms in c() as lower triangle, a generalised inverse ! of a sub-matrix of the positive semi-definite symmetric ! matrix a() of order n, stored as lower triangle. ! c() may co-incide with a(). nullty is returned as the nullity ! of a(). ifault is returned as 1 if n.lt.1, 2 if the ! submatrix is not positive semi-definite, 3 if the elements ! of b are inadmissible, otherwise zero. ! arguments:- ! a() = input, the symmetric matrix to be inverted, stored in ! lower triangular form ! nm = input, the order of a ! n = input, order of the sub-matrix to be inverted ! b = input, an array containing the numbers of rows and ! columns of a() to be included. ! c() = output, the inverse of a (a generalized inverse if c is ! singular), also stored in lower triangular. ! c and a may occupy the same locations. ! nullty = output, the rank deficiency of a. ! ifault = output, error indicator ! = 1 if n < 1 ! = 2 if the submatrix of a is not +ve semi-definite ! = 3 if elements of b are inadmissible ! = 0 otherwise ! Auxiliary routine required: SUBCHL from ASR44 (see AS6) !*************************************************************************** IMPLICIT NONE INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(15, 60) REAL (dp), INTENT(IN) :: a(:) INTEGER, INTENT(IN) :: nm, b(:), n REAL (dp), INTENT(OUT) :: c(:) REAL (dp), INTENT(OUT), OPTIONAL :: det INTEGER, INTENT(OUT) :: nullty, ifault INTERFACE SUBROUTINE subchl (a, b, n, u, nullty, ifault, det) IMPLICIT NONE INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(15, 60) REAL (dp), INTENT(IN) :: a(:) REAL (dp), INTENT(OUT) :: u(:) REAL (dp), INTENT(OUT), OPTIONAL :: det INTEGER, INTENT(IN) :: b(:), n INTEGER, INTENT(OUT) :: nullty, ifault END SUBROUTINE subchl END INTERFACE ! Local variables REAL (dp) :: w(n), x REAL (dp), PARAMETER :: zero = 0.0_dp, one = 1.0_dp INTEGER :: i, icol, irow, j, jcol, k, l, mdiag, ndiag, nn, nrow ifault = 3 IF (n > nm) RETURN IF (b(1) < 1 .OR. b(1) > nm - n + 1) RETURN IF (n == 1) GO TO 19 DO i = 2, n IF (b(i) <= b(i-1) .OR. b(i) > nm - n + i) RETURN END DO 19 nrow = n ifault = 1 IF(nrow <= 0) GO TO 100 ifault = 0 ! Cholesky factorization of a, result in c IF (PRESENT(det)) THEN CALL subchl(a, b, nrow, c, nullty, ifault, det) ELSE CALL subchl(a, b, nrow, c, nullty, ifault ) END IF IF(ifault /= 0) GO TO 100 ! invert c & form the product (cinv)'*cinv, where cinv is the inverse ! of c, row by row starting with the last row. ! irow = the row number, ndiag = location of last element in the row. nn = nrow*(nrow + 1)/2 irow = nrow ndiag = nn DO irow = nrow, 1, -1 IF(c(ndiag) == zero) GO TO 60 l = ndiag DO i = irow, nrow w(i) = c(l) l = l + i END DO icol = nrow jcol = nn mdiag = nn DO l = jcol x = zero IF(icol == irow) x = one / w(irow) k = nrow DO IF(k == irow) EXIT x = x - w(k)*c(l) k = k - 1 l = l - 1 IF(l > mdiag) l = l - k + 1 END DO c(l) = x / w(irow) IF(icol == irow) GO TO 80 mdiag = mdiag - icol icol = icol - 1 jcol = jcol - 1 END DO 60 l = ndiag DO j = irow, nrow c(l) = zero l = l + j END DO 80 ndiag = ndiag - irow END DO 100 RETURN END SUBROUTINE subinv PROGRAM test_as6_and_as7 IMPLICIT NONE INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(15, 60) REAL (dp), ALLOCATABLE :: a(:), c(:), aa(:), u(:) INTEGER, ALLOCATABLE :: b(:) INTEGER :: n, nn, nullty, ifault, status, row, col, pos1, pos2, & i, k, end_row, destn REAL (dp) :: total, det, maxerr INTERFACE SUBROUTINE syminv(a, n, nn, c, nullty, ifault) IMPLICIT NONE INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(15, 60) REAL (dp), INTENT(IN) :: a(:) INTEGER, INTENT(IN) :: n, nn REAL (dp), INTENT(OUT) :: c(:) INTEGER, INTENT(OUT) :: nullty, ifault END SUBROUTINE syminv SUBROUTINE subinv(a, nm, b, n, c, nullty, ifault, det) IMPLICIT NONE INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(15, 60) REAL (dp), INTENT(IN) :: a(:) INTEGER, INTENT(IN) :: nm, b(:), n REAL (dp), INTENT(OUT) :: c(:) REAL (dp), INTENT(OUT), OPTIONAL :: det INTEGER, INTENT(OUT) :: nullty, ifault END SUBROUTINE subinv END INTERFACE WRITE(*, *) 'Enter number of rows for matrix: ' READ(*, *) n nn = n*(n+1)/2 ALLOCATE( a(nn), c(nn), aa(nn), u(nn), b(n), STAT=status ) IF (status /= 0) THEN WRITE(*, *) 'Unable to allocate space for arrays' STOP END IF ! Generate random lower-triangular matrix. CALL RANDOM_NUMBER( c ) ! aa = c.c' ! This code could be used if aa & c occupied the same locations. det = 1.0_dp end_row = nn destn = nn DO row = n, 1, -1 pos1 = end_row DO col = row, 1, -1 pos2 = destn IF (col == row) det = det * c(destn) total = c(pos1)*c(pos2) DO k = col-1, 1, -1 pos1 = pos1 - 1 pos2 = pos2 - 1 total = total + c(pos1)*c(pos2) END DO aa(destn) = total destn = destn - 1 pos1 = pos1 - 1 END DO end_row = destn END DO WRITE(*, '(a, g13.5)') ' Determinant of original Cholesky factor = ', det WRITE(*, '(a, g13.5)') ' Expected determinant of inv(A''A) = ', & 1.0_dp / det**2 a = aa ! Form the Cholesky factorization, then invert & calculate inv(A). CALL syminv(a, n, nn, u, nullty, ifault) IF (nullty /= 0) WRITE(*, *) 'NULLTY = ', nullty IF (ifault /= 0) WRITE(*, *) 'IFAULT = ', ifault ! Reverse the process DO i = 1, n b(i) = i END DO ! Apart from signs, u should be identical to the original c matrix. CALL subinv(u, n, b, n, a, nullty, ifault, det) IF (nullty /= 0) WRITE(*, *) 'NULLTY = ', nullty IF (ifault /= 0) WRITE(*, *) 'IFAULT = ', ifault WRITE(*, '(a, g13.5)') ' Determinant now = ', det ! Calculate maximum absolute difference between elements of A & AA. maxerr = 0.0_dp DO i = 1, nn maxerr = MAX( maxerr, ABS( a(i) - aa(i) ) ) END DO WRITE(*, '(a, g12.4)') ' Max. abs. error = ', maxerr DEALLOCATE( a, c, aa, u, b ) STOP END PROGRAM test_as6_and_as7