subroutine i_swap ( i, j ) ! !******************************************************************************* ! !! I_SWAP swaps two integer values. ! ! ! Modified: ! ! 30 November 1998 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, integer I, J. On output, the values of I and ! J have been interchanged. ! implicit none ! integer i integer j integer k ! k = i i = j j = k return end function isamax ( n, x, incx ) ! !******************************************************************************* ! !! ISAMAX finds the index of the vector element of maximum absolute value. ! ! ! Modified: ! ! 08 April 1999 ! ! Reference: ! ! Lawson, Hanson, Kincaid, Krogh, ! Basic Linear Algebra Subprograms for Fortran Usage, ! Algorithm 539, ! ACM Transactions on Mathematical Software, ! Volume 5, Number 3, September 1979, pages 308-323. ! ! Dongarra, Moler, Bunch and Stewart, ! LINPACK User's Guide, ! SIAM, (Society for Industrial and Applied Mathematics), ! 3600 University City Science Center, ! Philadelphia, PA, 19104-2688. ! ISBN 0-89871-172-X ! ! Parameters: ! ! Input, integer N, the number of entries in the vector. ! ! Input, real X(*), the vector to be examined. ! ! Input, integer INCX, the increment between successive entries of SX. ! ! Output, integer ISAMAX, the index of the element of SX of maximum ! absolute value. ! implicit none ! integer i integer incx integer isamax integer ix integer n real samax real x(*) ! if ( n <= 0 ) then isamax = 0 else if ( n == 1 ) then isamax = 1 else if ( incx == 1 ) then isamax = 1 samax = abs ( x(1) ) do i = 2, n if ( abs ( x(i) ) > samax ) then isamax = i samax = abs ( x(i) ) end if end do else if ( incx >= 0 ) then ix = 1 else ix = ( - n + 1 ) * incx + 1 end if isamax = 1 samax = abs ( x(ix) ) ix = ix + incx do i = 2, n if ( abs ( x(ix) ) > samax ) then isamax = i samax = abs ( x(ix) ) end if ix = ix + incx end do end if return end subroutine r_swap ( x, y ) ! !******************************************************************************* ! !! R_SWAP swaps two real values. ! ! ! Modified: ! ! 01 May 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, real X, Y. On output, the values of X and ! Y have been interchanged. ! implicit none ! real x real y real z ! z = x x = y y = z return end function samax ( n, x, incx ) ! !******************************************************************************* ! !! SAMAX returns the maximum absolute value of the entries in a vector. ! ! ! Modified: ! ! 08 April 1999 ! ! Parameters: ! ! Input, integer N, the number of entries in the vector. ! ! Input, real X(*), the vector to be examined. ! ! Input, integer INCX, the increment between successive entries of X. ! ! Output, real SAMAX, the maximum absolute value of an element of X. ! implicit none ! integer i integer incx integer ix integer n real samax real x(*) ! if ( n <= 0 ) then samax = 0.0E+00 else if ( n == 1 ) then samax = abs ( x(1) ) else if ( incx == 1 ) then samax = abs ( x(1) ) do i = 2, n if ( abs ( x(i) ) > samax ) then samax = abs ( x(i) ) end if end do else if ( incx >= 0 ) then ix = 1 else ix = ( - n + 1 ) * incx + 1 end if samax = abs ( x(ix) ) ix = ix + incx do i = 2, n if ( abs ( x(ix) ) > samax ) then samax = abs ( x(ix) ) end if ix = ix + incx end do end if return end function sasum ( n, x, incx ) ! !******************************************************************************* ! !! SASUM sums the absolute values of the entries of a vector. ! ! ! Modified: ! ! 15 February 2001 ! ! Reference: ! ! Lawson, Hanson, Kincaid, Krogh, ! Basic Linear Algebra Subprograms for Fortran Usage, ! Algorithm 539, ! ACM Transactions on Mathematical Software, ! Volume 5, Number 3, September 1979, pages 308-323. ! ! Dongarra, Moler, Bunch and Stewart, ! LINPACK User's Guide, ! SIAM, (Society for Industrial and Applied Mathematics), ! 3600 University City Science Center, ! Philadelphia, PA, 19104-2688. ! ISBN 0-89871-172-X ! ! Parameters: ! ! Input, integer N, the number of entries in the vector. ! ! Input, real X(*), the vector to be examined. ! ! Input, integer INCX, the increment between successive entries of X. ! INCX must not be negative. ! ! Output, real SASUM, the sum of the absolute values of X. ! implicit none ! integer incx integer n real sasum real x(*) ! sasum = sum ( abs ( x(1:1+(n-1)*incx:incx) ) ) return end subroutine saxpy ( n, sa, x, incx, y, incy ) ! !******************************************************************************* ! !! SAXPY adds a constant times one vector to another. ! ! ! Modified: ! ! 08 April 1999 ! ! Reference: ! ! Lawson, Hanson, Kincaid, Krogh, ! Basic Linear Algebra Subprograms for Fortran Usage, ! Algorithm 539, ! ACM Transactions on Mathematical Software, ! Volume 5, Number 3, September 1979, pages 308-323. ! ! Dongarra, Moler, Bunch and Stewart, ! LINPACK User's Guide, ! SIAM, (Society for Industrial and Applied Mathematics), ! 3600 University City Science Center, ! Philadelphia, PA, 19104-2688. ! ISBN 0-89871-172-X ! ! Parameters: ! ! Input, integer N, the number of entries in the vector. ! ! Input, real SA, the multiplier. ! ! Input, real X(*), the vector to be scaled and added to Y. ! ! Input, integer INCX, the increment between successive entries of X. ! ! Input/output, real Y(*), the vector to which a multiple of X is to ! be added. ! ! Input, integer INCY, the increment between successive entries of Y. ! implicit none ! integer i integer incx integer incy integer ix integer iy integer n real sa real x(*) real y(*) ! if ( n <= 0 ) then else if ( sa == 0.0E+00 ) then else if ( incx == 1 .and. incy == 1 ) then y(1:n) = y(1:n) + sa * x(1:n) else if ( incx >= 0 ) then ix = 1 else ix = ( - n + 1 ) * incx + 1 end if if ( incy >= 0 ) then iy = 1 else iy = ( - n + 1 ) * incy + 1 end if do i = 1, n y(iy) = y(iy) + sa * x(ix) ix = ix + incx iy = iy + incy end do end if return end subroutine schdc ( a, lda, p, work, ipvt, job, info ) ! !******************************************************************************* ! !! SCHDC computes the Cholesky decomposition of a positive definite matrix. ! ! ! Discussion: ! ! A pivoting option allows the user to estimate the condition of a ! positive definite matrix or determine the rank of a positive ! semidefinite matrix. ! ! For positive definite matrices, INFO = P is the normal return. ! ! For pivoting with positive semidefinite matrices, INFO will ! in general be less than P. However, INFO may be greater than ! the rank of A, since rounding error can cause an otherwise zero ! element to be positive. Indefinite systems will always cause ! INFO to be less than P. ! ! Reference: ! ! Dongarra, Moler, Bunch and Stewart, ! LINPACK User's Guide, ! SIAM, (Society for Industrial and Applied Mathematics), ! 3600 University City Science Center, ! Philadelphia, PA, 19104-2688. ! ISBN 0-89871-172-X ! ! Author: ! ! J Dongarra and G Stewart, ! Argonne National Laboratory and University of Maryland. ! ! Parameters: ! ! Input/output, real A(LDA,P). ! On input, A contains the matrix whose decomposition is to ! be computed. Only the upper half of A need be stored. ! The lower part of the array a is not referenced. ! On output, A contains in its upper half the Cholesky factor ! of the input matrix, as it has been permuted by pivoting. ! ! Input, integer LDA, the leading dimension of the array A. ! ! Input, integer P, the order of the matrix. ! ! Input, real WORK(P) is a work array. ! ! Input/output, integer IPVT(P). ! On input, IPVT contains integers that control the selection ! of the pivot elements, if pivoting has been requested. ! Each diagonal element A(K,K) is placed in one of three classes ! according to the value of IPVT(K). ! ! > 0, then X(K) is an initial element. ! = 0, then X(K) is a free element. ! < 0, then X(K) is a final element. ! ! Before the decomposition is computed, initial elements are moved by ! symmetric row and column interchanges to the beginning of the array A ! and final elements to the end. Both initial and final elements are ! frozen in place during the computation and only free elements are moved. ! At the K-th stage of the reduction, if A(K,K) is occupied by a free ! element, it is interchanged with the largest free element A(L,L) with ! L >= K. IPVT is not referenced if JOB is 0. ! ! On output, IPVT(J) contains the index of the diagonal element ! of A that was moved into the J-th position, if pivoting was requested. ! ! Input, integer JOB, initiates column pivoting. ! 0, no pivoting is done. ! nonzero, pivoting is done. ! ! Output, integer INFO, contains the index of the last positive diagonal ! element of the Cholesky factor. ! implicit none ! integer lda integer p ! real a(lda,p) integer info integer j integer job integer jp integer ipvt(p) integer jt integer k integer l real maxdia integer maxl logical negk integer pl integer pu logical swapk real temp real work(*) ! pl = 1 pu = 0 info = p if ( job /= 0 ) then ! ! Pivoting has been requested. ! Rearrange the the elements according to IPVT. ! do k = 1, p swapk = ipvt(k) > 0 negk = ipvt(k) < 0 if ( negk ) then ipvt(k) = -k else ipvt(k) = k end if if ( swapk ) then if ( k /= pl ) then call sswap ( pl-1, a(1,k), 1, a(1,pl), 1 ) call r_swap ( a(k,k), a(pl,pl) ) do j = pl+1, p if ( j < k ) then call r_swap ( a(pl,j), a(j,k) ) else if ( j > k ) then call r_swap ( a(pl,j), a(k,j) ) end if end do ipvt(k) = ipvt(pl) ipvt(pl) = k end if pl = pl + 1 end if end do pu = p do k = p, pl, -1 if ( ipvt(k) < 0 ) then ipvt(k) = -ipvt(k) if ( pu /= k ) then call sswap ( k-1, a(1,k), 1, a(1,pu), 1 ) call r_swap ( a(k,k), a(pu,pu) ) do j = k+1, p if ( j < pu ) then call r_swap ( a(k,j), a(j,pu) ) else if ( j > pu ) then call r_swap ( a(k,j), a(pu,j) ) end if end do call i_swap ( ipvt(k), ipvt(pu) ) end if pu = pu - 1 end if end do end if do k = 1, p ! ! Reduction loop. ! maxdia = a(k,k) maxl = k ! ! Determine the pivot element. ! if ( k >= pl .and. k < pu ) then do l = k+1, pu if ( a(l,l) > maxdia ) then maxdia = a(l,l) maxl = l end if end do end if ! ! Quit if the pivot element is not positive. ! if ( maxdia <= 0.0E+00 ) then info = k - 1 return end if ! ! Start the pivoting and update IPVT. ! if ( k /= maxl ) then call sswap ( k-1, a(1,k), 1, a(1,maxl), 1 ) a(maxl,maxl) = a(k,k) a(k,k) = maxdia call i_swap ( ipvt(maxl), ipvt(k) ) end if ! ! Reduction step. ! Pivoting is contained across the rows. ! work(k) = sqrt ( a(k,k) ) a(k,k) = work(k) do j = k+1, p if ( k /= maxl ) then if ( j < maxl ) then call r_swap ( a(k,j), a(j,maxl) ) else if ( j > maxl ) then call r_swap ( a(k,j), a(maxl,j) ) end if end if a(k,j) = a(k,j) / work(k) work(j) = a(k,j) temp = -a(k,j) call saxpy ( j-k, temp, work(k+1), 1, a(k+1,j), 1 ) end do end do return end subroutine schdd ( r, ldr, p, x, z, ldz, nz, y, rho, c, s, info ) ! !******************************************************************************* ! !! SCHDD downdates an augmented Cholesky decomposition. ! ! ! Discussion: ! ! SCHDD can also downdate the triangular factor of an augmented QR ! decomposition. ! ! Specifically, given an upper triangular matrix R of order P, a ! row vector X, a column vector Z, and a scalar Y, SCHDD ! determines an orthogonal matrix U and a scalar ZETA such that ! ! ( R Z ) ( RR ZZ ) ! U * ( ) = ( ), ! ( 0 ZETA ) ( X Y ) ! ! where RR is upper triangular. ! ! If R and Z have been obtained from the factorization of a least squares ! problem, then RR and ZZ are the factors corresponding to the problem ! with the observation (X,Y) removed. In this case, if RHO ! is the norm of the residual vector, then the norm of ! the residual vector of the downdated problem is ! sqrt ( RHO**2 - ZETA**2 ). SCHDD will simultaneously downdate ! several triplets (Z, Y, RHO) along with R. ! ! For a less terse description of what SCHDD does and how ! it may be applied, see the LINPACK guide. ! ! The matrix U is determined as the product U(1)*...*U(P) ! where U(I) is a rotation in the (P+1,I)-plane of the form ! ! ( C(I) -S(I) ) ! ( ). ! ( S(I) C(I) ) ! ! The rotations are chosen so that C(I) is real. ! ! The user is warned that a given downdating problem may be impossible ! to accomplish or may produce inaccurate results. For example, this ! can happen if X is near a vector whose removal will reduce the ! rank of R. Beware. ! ! Modified: ! ! 17 April 2002 ! ! Reference: ! ! Dongarra, Moler, Bunch and Stewart, ! LINPACK User's Guide, ! SIAM, (Society for Industrial and Applied Mathematics), ! 3600 University City Science Center, ! Philadelphia, PA, 19104-2688. ! ISBN 0-89871-172-X ! ! Parameters: ! ! Input/output, real R(LDR,P), the upper triangular matrix that is to be ! downdated. The part of R below the diagonal is not referenced. ! ! Input, integer LDR, the leading dimension of the array R. ! LDR must be at least P. ! ! Input, integer P, the order of the matrix R. ! ! Input, real X(P), the row vector that is to be removed from R. ! ! Input/output, real Z(LDZ,NZ), an array of NZ P-vectors which are to ! be downdated along with R. ! ! Input, integer LDZ, the leading dimension of the array Z. ! LDZ must be at least P. ! ! Input, integer NZ, the number of vectors to be downdated. ! NZ may be zero, in which case Z, Y, and RHO are not referenced. ! ! Input, real Y(NZ), the scalars for the downdating of the vectors Z. ! ! Input/output, real RHO(NZ), the norms of the residual vectors. On ! output these have been changed along with R and Z. ! ! Output, real C(P), S(P), the cosines and sines of the transforming ! rotations. ! ! Output, integer INFO, return flag. ! 0, the entire downdating was successful. ! -1, if R could not be downdated. In this case, all quantities ! are left unaltered. ! 1, if some RHO could not be downdated. The offending RHO's are ! set to -1. ! implicit none ! integer ldr integer ldz integer nz integer p ! real a real alpha real azeta real b real c(p) integer i integer ii integer info integer j real norm real r(ldr,p) real rho(nz) real s(p) real scale real snrm2 real t real x(p) real xx real y(nz) real z(ldz,nz) real zeta ! ! Solve the system R'*A = X, placing the result in the array S. ! info = 0 s(1) = x(1) / r(1,1) do j = 2, p s(j) = x(j) - dot_product ( r(1:j-1,j), s(1:j-1) ) s(j) = s(j) / r(j,j) end do norm = snrm2 ( p, s, 1 ) if ( norm >= 1.0E+00 ) then info = -1 return end if alpha = sqrt ( 1.0E+00 - norm**2 ) ! ! Determine the transformations. ! do ii = 1, p i = p - ii + 1 scale = alpha + abs ( s(i) ) a = alpha / scale b = s(i) / scale norm = sqrt ( a**2 + b**2 ) c(i) = a / norm s(i) = b / norm alpha = scale * norm end do ! ! Apply the transformations to R. ! do j = 1, p xx = 0.0E+00 do ii = 1, j i = j - ii + 1 t = c(i) * xx + s(i) * r(i,j) r(i,j) = c(i) * r(i,j) - s(i) * xx xx = t end do end do ! ! If required, downdate Z and RHO. ! do j = 1, nz zeta = y(j) do i = 1, p z(i,j) = ( z(i,j) - s(i) * zeta ) / c(i) zeta = c(i) * zeta - s(i) * z(i,j) end do azeta = abs ( zeta ) if ( azeta > rho(j) ) then info = 1 rho(j) = -1.0E+00 else rho(j) = rho(j) * sqrt ( 1.0E+00 - ( azeta / rho(j) )**2 ) end if end do return end subroutine schex ( r, ldr, p, k, l, z, ldz, nz, c, s, job ) ! !******************************************************************************* ! !! SCHEX updates the Cholesky factorization of a positive definite matrix. ! ! ! Discussion: ! ! The factorization has the form ! ! A = R' * R ! ! where A is a positive definite matrix of order P. ! ! The updating involves diagonal permutations of the form ! ! E' * A * E ! ! where E is a permutation matrix. Specifically, given ! an upper triangular matrix R and a permutation matrix ! E (which is specified by K, L, and JOB), SCHEX determines ! a orthogonal matrix U such that ! ! U * R * E = RR, ! ! where RR is upper triangular. At the user's option, the ! transformation U will be multiplied into the array Z. ! If A = X'*X, so that R is the triangular part of the ! QR factorization of X, then RR is the triangular part of the ! QR factorization of X*E, i.e. X with its columns permuted. ! For a less terse description of what SCHEX does and how ! it may be applied, see the LINPACK guide. ! ! The matrix Q is determined as the product U(L-K)*...*U(1) ! of plane rotations of the form ! ! ( C(I) S(I) ) ! ( ), ! ( -S(I) C(I) ) ! ! where C(I) is real, the rows these rotations operate on ! are described below. ! ! There are two types of permutations, which are determined ! by the value of JOB. ! ! 1, right circular shift. The columns are rearranged in the order: ! ! 1,...,K-1,L,K,K+1,...,L-1,L+1,...,P. ! ! U is the product of L-K rotations U(I), where U(I) ! acts in the (L-I,L-I+1)-plane. ! ! 2, left circular shift: the columns are rearranged in the order ! ! 1,...,K-1,K+1,K+2,...,L,K,L+1,...,P. ! ! U is the product of L-K rotations U(I), where U(I) ! acts in the (K+I-1,K+I)-plane. ! ! Reference: ! ! Dongarra, Moler, Bunch and Stewart, ! LINPACK User's Guide, ! SIAM, (Society for Industrial and Applied Mathematics), ! 3600 University City Science Center, ! Philadelphia, PA, 19104-2688. ! ISBN 0-89871-172-X ! ! Parameters: ! ! Input/output, real R(LDR,P). On input, the upper triangular factor ! that is to be updated. Elements of R below the diagonal are not ! referenced. On output, R has been updated. ! ! Input, integer LDR, the leading dimension of the array R. ! LDR must be at least P. ! ! Input, integer P, the order of the matrix R. ! ! Input, integer K, the first column to be permuted. ! ! Input, integer L, the last column to be permuted. ! L must be strictly greater than K. ! ! Input/output real Z(LDZ,NZ), an array of NZ P-vectors into which the ! transformation U is multiplied. Z is not referenced if NZ = 0. ! On output, Z has been updated. ! ! Input, integer LDZ, the leading dimension of the array Z. ! LDZ must be at least P. ! ! Input, integer NZ, the number of columns of the matrix Z. ! ! Input, integer JOB, determines the type of permutation. ! 1, right circular shift. ! 2, left circular shift. ! ! Output, real C(P), S(P), the cosines and sines of the transforming ! rotations. ! implicit none ! integer ldr integer ldz integer p integer nz ! real c(p) integer job integer i integer ii integer il integer iu integer j integer jj integer k integer l integer lm1 integer lmk real r(ldr,p) real s(p) real t real z(ldz,nz) ! ! Initialize ! lmk = l - k lm1 = l - 1 ! ! Right circular shift. ! if ( job == 1 ) then ! ! Reorder the columns. ! do i = 1, l ii = l - i + 1 s(i) = r(ii,l) end do do jj = k, lm1 j = lm1 - jj + k do i = 1, j r(i,j+1) = r(i,j) end do r(j+1,j+1) = 0.0E+00 end do do i = 1, k-1 ii = l - i + 1 r(i,k) = s(ii) end do ! ! Calculate the rotations. ! t = s(1) do i = 1, lmk call srotg ( s(i+1), t, c(i), s(i) ) t = s(i+1) end do r(k,k) = t do j = k+1, p il = max (1,l-j+1) do ii = il, lmk i = l - ii t = c(ii) * r(i,j) + s(ii) * r(i+1,j) r(i+1,j) = c(ii) * r(i+1,j) - s(ii) * r(i,j) r(i,j) = t end do end do ! ! If required, apply the transformations to Z. ! do j = 1, nz do ii = 1, lmk i = l - ii t = c(ii) * z(i,j) + s(ii) * z(i+1,j) z(i+1,j) = c(ii) * z(i+1,j) - s(ii) * z(i,j) z(i,j) = t end do end do ! ! Left circular shift. ! else ! ! Reorder the columns. ! do i = 1, k ii = lmk + i s(ii) = r(i,k) end do do j = k, lm1 do i = 1, j r(i,j) = r(i,j+1) end do jj = j - k + 1 s(jj) = r(j+1,j+1) end do do i = 1, k ii = lmk + i r(i,l) = s(ii) end do r(k+1:l,l) = 0.0E+00 ! ! Reduction loop. ! do j = k, p ! ! Apply the rotations. ! if ( j /= k ) then iu = min ( j-1, l-1 ) do i = k, iu ii = i - k + 1 t = c(ii) * r(i,j) + s(ii) * r(i+1,j) r(i+1,j) = c(ii) * r(i+1,j) - s(ii) * r(i,j) r(i,j) = t end do end if if ( j < l ) then jj = j - k + 1 t = s(jj) call srotg ( r(j,j), t, c(jj), s(jj) ) end if end do ! ! Apply the rotations to Z. ! do j = 1, nz do i = k, lm1 ii = i - k + 1 t = c(ii) * z(i,j) + s(ii) * z(i+1,j) z(i+1,j) = c(ii) * z(i+1,j) - s(ii) * z(i,j) z(i,j) = t end do end do end if return end subroutine schud ( r, ldr, p, x, z, ldz, nz, y, rho, c, s ) ! !******************************************************************************* ! !! SCHUD updates an augmented Cholesky decomposition. ! ! ! Discussion: ! ! SCHUD can also update the triangular part of an augmented QR ! decomposition. ! ! Specifically, given an upper triangular matrix R of order P, a row vector ! X, a column vector Z, and a scalar Y, SCHUD determines a unitary matrix ! U and a scalar ZETA such that ! ! (R Z) (RR ZZ ) ! U * ( ) = ( ), ! (X Y) ( 0 ZETA) ! ! where RR is upper triangular. ! ! If R and Z have been obtained from the factorization of a least squares ! problem, then RR and ZZ are the factors corresponding to the problem ! with the observation (X,Y) appended. In this case, if RHO is the ! norm of the residual vector, then the norm of the residual vector of ! the updated problem is sqrt ( RHO**2 + ZETA**2 ). SCHUD will ! simultaneously update several triplets (Z, Y, RHO). ! ! For a less terse description of what SCHUD does and how ! it may be applied, see the LINPACK guide. ! ! The matrix U is determined as the product U(P)*...*U(1), ! where U(I) is a rotation in the (I,P+1) plane of the form ! ! ( C(I) S(I) ) ! ( ). ! ( -S(I) C(I) ) ! ! The rotations are chosen so that C(I) is real. ! ! Reference: ! ! Dongarra, Moler, Bunch and Stewart, ! LINPACK User's Guide, ! SIAM, (Society for Industrial and Applied Mathematics), ! 3600 University City Science Center, ! Philadelphia, PA, 19104-2688. ! ISBN 0-89871-172-X ! ! Parameters: ! ! Input, real R(LDR,P), the upper triangular matrix that is to be updated. ! the part of R below the diagonal is not referenced. ! ! Input, integer LDR, the leading dimension of the array R. ! LDR must be at least equal to P. ! ! Input, integer P, the order of the matrix R. ! ! Input, real X(P), the row to be added to R. ! ! Input/output, real Z(LDZ,NZ), contains NZ P-vectors to be updated with R. ! ! Input, integer LDZ, the leading dimension of the array Z. ! LDZ must be at least P. ! ! Input, integer NZ, the number of vectors to be updated. NZ may be ! zero, in which case Z, Y, and RHO are not referenced. ! ! Input, real Y(NZ), the scalars for updating the vectors Z. ! ! Input/output, real RHO(NZ). On input, the norms of the residual ! vectors that are to be updated. If RHO(J) is negative, it is left ! unaltered. ! ! Output, real C(P), S(P), the cosines and sines of the transforming ! rotations. ! implicit none ! integer ldr integer ldz integer nz integer p ! real azeta real c(p) integer i integer j real r(ldr,p) real rho(nz) real s(p) real scale real t real x(p) real xj real y(nz) real z(ldz,nz) real zeta ! ! Update R. ! do j = 1, p xj = x(j) ! ! Apply the previous rotations. ! do i = 1, j-1 t = c(i) * r(i,j) + s(i) * xj xj = c(i) * xj - s(i) * r(i,j) r(i,j) = t end do ! ! Compute the next rotation. ! call srotg ( r(j,j), xj, c(j), s(j) ) end do ! ! If required, update Z and RHO. ! do j = 1, nz zeta = y(j) do i = 1, p t = c(i) * z(i,j) + s(i) * zeta zeta = c(i) * zeta - s(i) * z(i,j) z(i,j) = t end do azeta = abs ( zeta ) if ( azeta /= 0.0E+00 .and. rho(j) >= 0.0 ) then scale = azeta + rho(j) rho(j) = scale * sqrt ( ( azeta / scale )**2 + ( rho(j) / scale )**2 ) end if end do return end subroutine scopy ( n, x, incx, y, incy ) ! !******************************************************************************* ! !! SCOPY copies one real vector into another. ! ! ! Modified: ! ! 08 April 1999 ! ! Reference: ! ! Lawson, Hanson, Kincaid, Krogh, ! Basic Linear Algebra Subprograms for Fortran Usage, ! Algorithm 539, ! ACM Transactions on Mathematical Software, ! Volume 5, Number 3, September 1979, pages 308-323. ! ! Dongarra, Moler, Bunch and Stewart, ! LINPACK User's Guide, ! SIAM, (Society for Industrial and Applied Mathematics), ! 3600 University City Science Center, ! Philadelphia, PA, 19104-2688. ! ISBN 0-89871-172-X ! ! Parameters: ! ! Input, integer N, the number of entries in the vector. ! ! Input, real X(*), the vector to be copied into Y. ! ! Input, integer INCX, the increment between successive entries of X. ! ! Output, real Y(*), the copy of X. ! ! Input, integer INCY, the increment between successive elements of Y. ! implicit none ! integer i integer incx integer incy integer ix integer iy integer n real x(*) real y(*) ! if ( n <= 0 ) then else if ( incx == 1 .and. incy == 1 ) then y(1:n) = x(1:n) else if ( incx >= 0 ) then ix = 1 else ix = ( - n + 1 ) * incx + 1 end if if ( incy >= 0 ) then iy = 1 else iy = ( - n + 1 ) * incy + 1 end if do i = 1, n y(iy) = x(ix) ix = ix + incx iy = iy + incy end do end if return end function sdot ( n, x, incx, y, incy ) ! !******************************************************************************* ! !! SDOT forms the dot product of two vectors. ! ! ! Discussion: ! ! For many purposes, SDOT can be replaced by the FORTRAN 90 ! intrinsic DOT_PRODUCT. ! ! Modified: ! ! 02 June 2000 ! ! Reference: ! ! Lawson, Hanson, Kincaid, Krogh, ! Basic Linear Algebra Subprograms for Fortran Usage, ! Algorithm 539, ! ACM Transactions on Mathematical Software, ! Volume 5, Number 3, September 1979, pages 308-323. ! ! Dongarra, Moler, Bunch and Stewart, ! LINPACK User's Guide, ! SIAM, (Society for Industrial and Applied Mathematics), ! 3600 University City Science Center, ! Philadelphia, PA, 19104-2688. ! ISBN 0-89871-172-X ! ! Parameters: ! ! Input, integer N, the number of entries in the vectors. ! ! Input, real X(*), one of the vectors to be multiplied. ! ! Input, integer INCX, the increment between successive entries of X. ! ! Input, real Y(*), one of the vectors to be multiplied. ! ! Input, integer INCY, the increment between successive elements of Y. ! ! Output, real SDOT, the dot product of X and Y. ! implicit none ! integer i integer incx integer incy integer ix integer iy integer n real sdot real stemp real x(*) real y(*) ! if ( n <= 0 ) then sdot = 0.0E+00 else if ( incx == 1 .and. incy == 1 ) then sdot = dot_product ( x(1:n), y(1:n) ) else if ( incx >= 0 ) then ix = 1 else ix = ( - n + 1 ) * incx + 1 end if if ( incy >= 0 ) then iy = 1 else iy = ( - n + 1 ) * incy + 1 end if stemp = 0.0E+00 do i = 1, n stemp = stemp + x(ix) * y(iy) ix = ix + incx iy = iy + incy end do sdot = stemp end if return end subroutine sgbco ( abd, lda, n, ml, mu, ipvt, rcond, z ) ! !******************************************************************************* ! !! SGBCO factors a real band matrix and estimates its condition. ! ! ! Discussion: ! ! If RCOND is not needed, SGBFA is slightly faster. ! ! To solve A*X = B, follow SGBCO by SGBSL. ! ! To compute inverse(A)*C, follow SGBCO by SGBSL. ! ! To compute determinant(A), follow SGBCO by SGBDI. ! ! Example: ! ! If the original matrix is ! ! 11 12 13 0 0 0 ! 21 22 23 24 0 0 ! 0 32 33 34 35 0 ! 0 0 43 44 45 46 ! 0 0 0 54 55 56 ! 0 0 0 0 65 66 ! ! then for proper band storage, ! ! N = 6, ML = 1, MU = 2, LDA >= 5 and ABD should contain ! ! * * * + + + * = not used ! * * 13 24 35 46 + = used for pivoting ! * 12 23 34 45 56 ! 11 22 33 44 55 66 ! 21 32 43 54 65 * ! ! Band storage: ! ! If A is a band matrix, the following program segment ! will set up the input. ! ! ml = (band width below the diagonal) ! mu = (band width above the diagonal) ! m = ml + mu + 1 ! ! do j = 1, n ! i1 = max ( 1, j-mu ) ! i2 = min ( n, j+ml ) ! do i = i1, i2 ! k = i - j + m ! abd(k,j) = a(i,j) ! end do ! end do ! ! This uses rows ML+1 through 2*ML+MU+1 of ABD. In addition, the first ! ML rows in ABD are used for elements generated during the ! triangularization. ! ! The total number of rows needed in ABD is 2*ML+MU+1. The ML+MU by ! ML+MU upper left triangle and the ML by ML lower right triangle are ! not referenced. ! ! Reference: ! ! Dongarra, Moler, Bunch and Stewart, ! LINPACK User's Guide, ! SIAM, (Society for Industrial and Applied Mathematics), ! 3600 University City Science Center, ! Philadelphia, PA, 19104-2688. ! ISBN 0-89871-172-X ! ! Parameters: ! ! Input/output, real ABD(LDA,N). On input, the matrix in band storage. ! The columns of the matrix are stored in the columns of ABD and ! the diagonals of the matrix are stored in rows ML+1 through 2*ML+MU+1 ! of ABD. On output, an upper triangular matrix in band storage and ! the multipliers which were used to obtain it. The factorization can ! be written A = L*U where L is a product of permutation and unit lower ! triangular matrices and U is upper triangular. ! ! Input, integer LDA, the leading dimension of the array ABD. ! LDA must be >= 2*ML + MU + 1. ! ! Input, integer N, the order of the matrix. ! ! Input, integer ML, MU, the number of diagonals below and above the ! main diagonal. 0 <= ML < N, 0 <= MU < N. ! ! Output, integer IPVT(N), the pivot indices. ! ! Output, real RCOND, an estimate of the reciprocal condition of A. ! For the system A*X = B, relative perturbations in A and B of size ! EPSILON may cause relative perturbations in X of size EPSILON/RCOND. ! If RCOND is so small that the logical expression ! 1.0E+00 + RCOND == 1.0E+00 ! is true, then A may be singular to working precision. In particular, ! RCOND is zero if exact singularity is detected or the estimate underflows. ! ! Workspace, real Z(N), a work vector whose contents are usually unimportant. ! If A is close to a singular matrix, then Z is an approximate null vector ! in the sense that ! norm(A*Z) = RCOND * norm(A) * norm(Z). ! implicit none ! integer lda integer n ! real abd(lda,n) real anorm real ek integer info integer ipvt(n) integer is integer j integer ju integer k integer l integer la integer lm integer lz integer m integer ml integer mm integer mu real rcond real s real sm real t real wk real wkm real ynorm real z(n) ! ! Compute the 1-norm of A. ! anorm = 0.0E+00 l = ml + 1 is = l + mu do j = 1, n anorm = max ( anorm, sum ( abs ( abd(is:is+l-1,j) ) ) ) if ( is > ml + 1 ) is = is - 1 if ( j <= mu ) l = l + 1 if ( j >= n - ml ) l = l - 1 end do ! ! Factor. ! call sgbfa ( abd, lda, n, ml, mu, ipvt, info ) ! ! RCOND = 1/(norm(A)*(estimate of norm(inverse(A)))). ! ! Estimate = norm(Z)/norm(Y) where a*z = y and A'*Y = E. ! ! A' is the transpose of A. The components of E are ! chosen to cause maximum local growth in the elements of W where ! U'*W = E. The vectors are frequently rescaled to avoid ! overflow. ! ! Solve U'*W = E. ! ek = 1.0E+00 z(1:n) = 0.0E+00 m = ml + mu + 1 ju = 0 do k = 1, n if ( z(k) /= 0.0E+00 ) then ek = sign ( ek, -z(k) ) end if if ( abs ( ek - z(k) ) > abs ( abd(m,k) ) ) then s = abs ( abd(m,k) ) / abs ( ek - z(k) ) z(1:n) = s * z(1:n) ek = s * ek end if wk = ek - z(k) wkm = -ek - z(k) s = abs ( wk ) sm = abs ( wkm ) if ( abd(m,k) /= 0.0E+00 ) then wk = wk / abd(m,k) wkm = wkm / abd(m,k) else wk = 1.0E+00 wkm = 1.0E+00 end if ju = min ( max ( ju, mu+ipvt(k) ), n ) mm = m if ( k+1 <= ju ) then do j = k+1, ju mm = mm - 1 sm = sm + abs ( z(j ) + wkm * abd(mm,j) ) z(j) = z(j) + wk * abd(mm,j) s = s + abs ( z(j) ) end do if ( s < sm ) then t = wkm - wk wk = wkm mm = m do j = k+1, ju mm = mm - 1 z(j) = z(j) + t * abd(mm,j) end do end if end if z(k) = wk end do z(1:n) = z(1:n) / sum ( abs ( z(1:n) ) ) ! ! Solve L'*Y = W. ! do k = n, 1, -1 lm = min ( ml, n-k ) if ( k < n ) then z(k) = z(k) + dot_product ( abd(m+1:m+lm,k), z(k+1:k+lm) ) end if if ( abs ( z(k) ) > 1.0E+00 ) then s = 1.0E+00 / abs ( z(k) ) z(1:n) = s * z(1:n) end if l = ipvt(k) t = z(l) z(l) = z(k) z(k) = t end do z(1:n) = z(1:n) / sum ( abs ( z(1:n) ) ) ynorm = 1.0E+00 ! ! Solve L*V = Y. ! do k = 1, n l = ipvt(k) t = z(l) z(l) = z(k) z(k) = t lm = min ( ml, n-k ) if ( k < n ) then call saxpy ( lm, t, abd(m+1,k), 1, z(k+1), 1 ) end if if ( abs ( z(k) ) > 1.0E+00 ) then s = 1.0E+00 / abs ( z(k) ) z(1:n) = s * z(1:n) ynorm = s * ynorm end if end do s = 1.0E+00 / sum ( abs ( z(1:n) ) ) z(1:n) = s * z(1:n) ynorm = s * ynorm ! ! Solve U*Z = W. ! do k = n, 1, -1 if ( abs ( z(k) ) > abs ( abd(m,k) ) ) then s = abs ( abd(m,k) ) / abs ( z(k) ) z(1:n) = s * z(1:n) ynorm = s * ynorm end if if ( abd(m,k) /= 0.0E+00 ) then z(k) = z(k) / abd(m,k) else z(k) = 1.0E+00 end if lm = min ( k, m ) - 1 la = m - lm lz = k - lm t = -z(k) call saxpy ( lm, t, abd(la,k), 1, z(lz), 1 ) end do ! ! Make ZNORM = 1.0. ! s = 1.0E+00 / sum ( abs ( z(1:n) ) ) z(1:n) = s * z(1:n) ynorm = s * ynorm if ( anorm /= 0.0E+00 ) then rcond = ynorm / anorm else rcond = 0.0E+00 end if return end subroutine sgbdi ( abd, lda, n, ml, mu, ipvt, det ) ! !******************************************************************************* ! !! SGBDI computes the determinant of a band matrix factored by SGBCO or SGBFA. ! ! ! Discussion: ! ! If the inverse is needed, use SGBSL N times. ! ! Reference: ! ! Dongarra, Moler, Bunch and Stewart, ! LINPACK User's Guide, ! SIAM, (Society for Industrial and Applied Mathematics), ! 3600 University City Science Center, ! Philadelphia, PA, 19104-2688. ! ISBN 0-89871-172-X ! ! Parameters: ! ! Input, real ABD(LDA,N), the output from SGBCO or SGBFA. ! ! Input, integer LDA, the leading dimension of the array ABD. ! ! Input, integer N, the order of the matrix. ! ! Input, integer ML, MU, the number of diagonals below and above the ! main diagonal. 0 <= ML < N, 0 <= MU < N. ! ! Input, integer IPVT(N), the pivot vector from SGBCO or SGBFA. ! ! Output, real DET(2), the determinant of the original matrix. ! determinant = DET(1) * 10.0**DET(2) ! with 1.0E+00 <= abs ( DET(1) ) < 10.0E+00 or DET(1) = 0.0E+00. ! implicit none ! integer lda integer n ! real abd(lda,n) real det(2) integer i integer ipvt(n) integer m integer ml integer mu real, parameter :: ten = 10.0E+00 ! m = ml + mu + 1 det(1) = 1.0E+00 det(2) = 0.0E+00 do i = 1, n if ( ipvt(i) /= i ) then det(1) = -det(1) end if det(1) = abd(m,i) * det(1) if ( det(1) == 0.0E+00 ) then return end if do while ( abs ( det(1) ) < 1.0E+00 ) det(1) = ten * det(1) det(2) = det(2) - 1.0E+00 end do do while ( abs ( det(1) ) >= ten ) det(1) = det(1) / ten det(2) = det(2) + 1.0E+00 end do end do return end subroutine sgbfa ( abd, lda, n, ml, mu, ipvt, info ) ! !******************************************************************************* ! !! SGBFA factors a real band matrix by elimination. ! ! ! Discussion: ! ! SGBFA is usually called by SGBCO, but it can be called ! directly with a saving in time if RCOND is not needed. ! ! Reference: ! ! Dongarra, Moler, Bunch and Stewart, ! LINPACK User's Guide, ! SIAM, (Society for Industrial and Applied Mathematics), ! 3600 University City Science Center, ! Philadelphia, PA, 19104-2688. ! ISBN 0-89871-172-X ! ! Parameters: ! ! Input/output, real ABD(LDA,N). On input, contains the matrix in band ! storage. The columns of the matrix are stored in the columns of ABD ! and the diagonals of the matrix are stored in rows ML+1 through ! 2*ML+MU+1 of ABD. On output, an upper triangular matrix in band storage ! and the multipliers which were used to obtain it. The factorization ! can be written A = L*U where L is a product of permutation and unit lower ! triangular matrices and U is upper triangular. ! ! Input, integer LDA, the leading dimension of the array ABD. ! LDA must be >= 2*ML + MU + 1. ! ! Input, integer N, the order of the matrix. ! ! Input, integer ML, MU, the number of diagonals below and above the ! main diagonal. 0 <= ML < N, 0 <= MU < N. ! ! Output, integer IPVT(N), the pivot indices. ! ! Output, integer INFO, error flag. ! 0, normal value. ! K, if U(K,K) == 0.0E+00. This is not an error condition for this ! subroutine, but it does indicate that SGBSL will divide by zero if ! called. Use RCOND in SGBCO for a reliable indication of singularity. ! implicit none ! integer lda integer n ! real abd(lda,n) integer i integer i0 integer info integer ipvt(n) integer isamax integer j integer j0 integer j1 integer ju integer jz integer k integer l integer lm integer m integer ml integer mm integer mu real t ! m = ml + mu + 1 info = 0 ! ! Zero initial fill-in columns. ! j0 = mu + 2 j1 = min ( n, m ) - 1 do jz = j0, j1 i0 = m + 1 - jz abd(i0:ml,jz) = 0.0E+00 end do jz = j1 ju = 0 ! ! Gaussian elimination with partial pivoting. ! do k = 1, n-1 ! ! Zero next fill-in column. ! jz = jz + 1 if ( jz <= n ) then abd(1:ml,jz) = 0.0E+00 end if ! ! Find L = pivot index. ! lm = min ( ml, n-k ) l = isamax ( lm+1, abd(m,k), 1 ) + m - 1 ipvt(k) = l + k - m ! ! Zero pivot implies this column already triangularized. ! if ( abd(l,k) == 0.0E+00 ) then info = k ! ! Interchange if necessary. ! else if ( l /= m ) then call r_swap ( abd(l,k), abd(m,k) ) end if ! ! Compute multipliers. ! t = -1.0E+00 / abd(m,k) call sscal ( lm, t, abd(m+1,k), 1 ) ! ! Row elimination with column indexing. ! ju = min ( max ( ju, mu+ipvt(k) ), n ) mm = m do j = k+1, ju l = l - 1 mm = mm - 1 t = abd(l,j) if ( l /= mm ) then abd(l,j) = abd(mm,j) abd(mm,j) = t end if call saxpy ( lm, t, abd(m+1,k), 1, abd(mm+1,j), 1 ) end do end if end do ipvt(n) = n if ( abd(m,n) == 0.0E+00 ) then info = n end if return end subroutine sgbsl ( abd, lda, n, ml, mu, ipvt, b, job ) ! !******************************************************************************* ! !! SGBSL solves a real banded system factored by SGBCO or SGBFA. ! ! ! Discussion: ! ! SGBSL can solve either A * X = B or A' * X = B. ! ! A division by zero will occur if the input factor contains a ! zero on the diagonal. Technically this indicates singularity ! but it is often caused by improper arguments or improper ! setting of LDA. It will not occur if the subroutines are ! called correctly and if SGBCO has set RCOND > 0.0E+00 ! or SGBFA has set INFO == 0. ! ! To compute inverse(A) * C where C is a matrix with P columns: ! ! call sgbco ( abd, lda, n, ml, mu, ipvt, rcond, z ) ! ! if ( rcond is too small ) then ! exit ! end if ! ! do j = 1, p ! call sgbsl ( abd, lda, n, ml, mu, ipvt, c(1,j), 0 ) ! end do ! ! Reference: ! ! Dongarra, Moler, Bunch and Stewart, ! LINPACK User's Guide, ! SIAM, (Society for Industrial and Applied Mathematics), ! 3600 University City Science Center, ! Philadelphia, PA, 19104-2688. ! ISBN 0-89871-172-X ! ! Parameters: ! ! Input, real ABD(LDA,N), the output from SGBCO or SGBFA. ! ! Input, integer LDA, the leading dimension of the array ABD. ! ! Input, integer N, the order of the matrix. ! ! Input, integer ML, MU, the number of diagonals below and above the ! main diagonal. 0 <= ML < N, 0 <= MU < N. ! ! Input, integer IPVT(N), the pivot vector from SGBCO or SGBFA. ! ! Input/output, real B(N). On input, the right hand side. ! On output, the solution. ! ! Input, integer JOB, job choice. ! 0, solve A*X=B. ! nonzero, solve A'*X=B. ! implicit none ! integer lda integer n ! real abd(lda,n) real b(n) integer ipvt(n) integer job integer k integer l integer la integer lb integer lm integer m integer ml integer mu real t ! m = mu + ml + 1 ! ! JOB = 0, Solve a * x = b. ! ! First solve l*y = b. ! if ( job == 0 ) then if ( ml > 0 ) then do k = 1, n-1 lm = min ( ml, n-k ) l = ipvt(k) t = b(l) if ( l /= k ) then b(l) = b(k) b(k) = t end if call saxpy ( lm, t, abd(m+1,k), 1, b(k+1), 1 ) end do end if ! ! Now solve u*x = y. ! do k = n, 1, -1 b(k) = b(k) / abd(m,k) lm = min ( k, m ) - 1 la = m - lm lb = k - lm t = -b(k) call saxpy ( lm, t, abd(la,k), 1, b(lb), 1 ) end do ! ! JOB nonzero, solve A' * X = B. ! ! First solve U'*Y = B. ! else do k = 1, n lm = min ( k, m ) - 1 t = dot_product ( abd(m-lm:m-1,k), b(k-lm:k-1) ) b(k) = ( b(k) - t ) / abd(m,k) end do ! ! Now solve L'*X = Y. ! if ( ml > 0 ) then do k = n-1, 1, -1 lm = min ( ml, n-k ) b(k) = b(k) + dot_product ( abd(m+1:m+lm,k), b(k+1:k+lm) ) l = ipvt(k) if ( l /= k ) then call r_swap ( b(l), b(k) ) end if end do end if end if return end subroutine sgeco ( a, lda, n, ipvt, rcond, z ) ! !******************************************************************************* ! !! SGECO factors a real matrix and estimates its condition number. ! ! ! Discussion: ! ! If RCOND is not needed, SGEFA is slightly faster. ! ! To solve A * X = B, follow SGECO by SGESL. ! ! To compute inverse ( A ) * C, follow SGECO by SGESL. ! ! To compute determinant ( A ), follow SGECO by SGEDI. ! ! To compute inverse ( A ), follow SGECO by SGEDI. ! ! For the system A * X = B, relative perturbations in A and B ! of size EPSILON may cause relative perturbations in X of size ! EPSILON/RCOND. ! ! If RCOND is so small that the logical expression ! 1.0E+00 + RCOND == 1.0E+00 ! is true, then A may be singular to working precision. In particular, ! RCOND is zero if exact singularity is detected or the estimate ! underflows. ! ! Reference: ! ! Dongarra, Moler, Bunch and Stewart, ! LINPACK User's Guide, ! SIAM, (Society for Industrial and Applied Mathematics), ! 3600 University City Science Center, ! Philadelphia, PA, 19104-2688. ! ISBN 0-89871-172-X ! ! Author: ! ! Cleve Moler, ! University of New Mexico / Argonne National Lab. ! ! Parameters: ! ! Input/output, real A(LDA,N). On input, a matrix to be factored. ! On output, the LU factorization of the matrix. ! ! Input, integer LDA, the leading dimension of the array A. ! ! Input, integer N, the order of the matrix A. ! ! Output, integer IPVT(N), the pivot indices. ! ! Output, real RCOND, an estimate of the reciprocal condition number of A. ! ! Output, real Z(N), a work vector whose contents are usually unimportant. ! If A is close to a singular matrix, then Z is an approximate null vector ! in the sense that ! norm ( A * Z ) = RCOND * norm ( A ) * norm ( Z ). ! implicit none ! integer lda integer n ! real a(lda,n) real anorm real ek integer info integer ipvt(n) integer j integer k integer l real rcond real s real sm real t real wk real wkm real ynorm real z(n) ! ! Compute the L1 norm of A. ! anorm = 0.0E+00 do j = 1, n anorm = max ( anorm, sum ( abs ( a(1:n,j) ) ) ) end do ! ! Compute the LU factorization. ! call sgefa ( a, lda, n, ipvt, info ) ! ! RCOND = 1 / ( norm(A) * (estimate of norm(inverse(A))) ) ! ! estimate of norm(inverse(A)) = norm(Z) / norm(Y) ! ! where ! A * Z = Y ! and ! A' * Y = E ! ! The components of E are chosen to cause maximum local growth in the ! elements of W, where U'*W = E. The vectors are frequently rescaled ! to avoid overflow. ! ! Solve U' * W = E. ! ek = 1.0E+00 z(1:n) = 0.0E+00 do k = 1, n if ( z(k) /= 0.0E+00 ) then ek = sign ( ek, -z(k) ) end if if ( abs ( ek - z(k) ) > abs ( a(k,k) ) ) then s = abs ( a(k,k) ) / abs ( ek - z(k) ) z(1:n) = s * z(1:n) ek = s * ek end if wk = ek - z(k) wkm = -ek - z(k) s = abs ( wk ) sm = abs ( wkm ) if ( a(k,k) /= 0.0E+00 ) then wk = wk / a(k,k) wkm = wkm / a(k,k) else wk = 1.0E+00 wkm = 1.0E+00 end if if ( k+1 <= n ) then do j = k+1, n sm = sm + abs ( z(j) + wkm * a(k,j) ) z(j) = z(j) + wk * a(k,j) s = s + abs ( z(j) ) end do if ( s < sm ) then t = wkm - wk wk = wkm z(k+1:n) = z(k+1:n) + t * a(k,k+1:n) end if end if z(k) = wk end do z(1:n) = z(1:n) / sum ( abs ( z(1:n) ) ) ! ! Solve L' * Y = W ! do k = n, 1, -1 z(k) = z(k) + dot_product ( a(k+1:n,k), z(k+1:n) ) if ( abs ( z(k) ) > 1.0E+00 ) then z(1:n) = z(1:n) / abs ( z(k) ) end if l = ipvt(k) t = z(l) z(l) = z(k) z(k) = t end do z(1:n) = z(1:n) / sum ( abs ( z(1:n) ) ) ynorm = 1.0E+00 ! ! Solve L * V = Y. ! do k = 1, n l = ipvt(k) t = z(l) z(l) = z(k) z(k) = t z(k+1:n) = z(k+1:n) + t * a(k+1:n,k) if ( abs ( z(k) ) > 1.0E+00 ) then ynorm = ynorm / abs ( z(k) ) z(1:n) = z(1:n) / abs ( z(k) ) end if end do s = sum ( abs ( z(1:n) ) ) z(1:n) = z(1:n) / s ynorm = ynorm / s ! ! Solve U * Z = V. ! do k = n, 1, -1 if ( abs ( z(k) ) > abs ( a(k,k) ) ) then s = abs ( a(k,k) ) / abs ( z(k) ) z(1:n) = s * z(1:n) ynorm = s * ynorm end if if ( a(k,k) /= 0.0E+00 ) then z(k) = z(k) / a(k,k) else z(k) = 1.0E+00 end if z(1:k-1) = z(1:k-1) - z(k) * a(1:k-1,k) end do ! ! Normalize Z in the L1 norm. ! s = 1.0E+00 / sum ( abs ( z(1:n) ) ) z(1:n) = s * z(1:n) ynorm = s * ynorm if ( anorm /= 0.0E+00 ) then rcond = ynorm / anorm else rcond = 0.0E+00 end if return end subroutine sgedi ( a, lda, n, ipvt, det, work, job ) ! !******************************************************************************* ! !! SGEDI computes the determinant and inverse of a matrix factored by SGECO or SGEFA. ! ! ! Discussion: ! ! A division by zero will occur if the input factor contains ! a zero on the diagonal and the inverse is requested. ! It will not occur if the subroutines are called correctly ! and if SGECO has set RCOND > 0.0E+00 or SGEFA has set INFO == 0. ! ! Reference: ! ! Dongarra, Moler, Bunch and Stewart, ! LINPACK User's Guide, ! SIAM, (Society for Industrial and Applied Mathematics), ! 3600 University City Science Center, ! Philadelphia, PA, 19104-2688. ! ISBN 0-89871-172-X ! ! Parameters: ! ! Input/output, real A(LDA,N), on input, the N by N factored matrix. ! as output by SGECO or SGEFA. On output, contains the inverse ! matrix if requested. ! ! Input, integer LDA, the leading dimension of the array A. ! ! Input, integer N, the order of the matrix A. ! ! Input, integer IPVT(N), the pivot vector from SGECO or SGEFA. ! ! Workspace, real WORK(N). ! ! Output, real DET(2), the determinant of original matrix if requested. ! determinant = DET(1) * 10.0**DET(2) ! with 1.0E+00 <= abs ( DET(1) ) < 10.0E+00 ! or DET(1) == 0.0E+00. ! ! Input, integer JOB, specifies what is to be computed. ! 11, both determinant and inverse. ! 01, inverse only. ! 10, determinant only. ! implicit none ! integer lda integer n ! real a(lda,n) real det(2) integer i integer ipvt(n) integer j integer job integer k integer kp1 integer l real t real, parameter :: ten = 10.0E+00 real work(n) ! ! Compute the determinant. ! if ( job / 10 /= 0 ) then det(1) = 1.0E+00 det(2) = 0.0E+00 do i = 1, n if ( ipvt(i) /= i ) then det(1) = - det(1) end if det(1) = a(i,i) * det(1) if ( det(1) == 0.0E+00 ) then exit end if do while ( abs ( det(1) ) < 1.0E+00 ) det(1) = ten * det(1) det(2) = det(2) - 1.0E+00 end do do while ( abs ( det(1) ) >= ten ) det(1) = det(1) / ten det(2) = det(2) + 1.0E+00 end do end do end if ! ! Compute inverse(U). ! if ( mod ( job, 10 ) /= 0 ) then do k = 1, n a(k,k) = 1.0E+00 / a(k,k) t = - a(k,k) call sscal ( k-1, t, a(1,k), 1 ) do j = k+1, n t = a(k,j) a(k,j) = 0.0E+00 call saxpy ( k, t, a(1,k), 1, a(1,j), 1 ) end do end do ! ! Form inverse(U) * inverse(L). ! do k = n-1, 1, -1 do i = k+1, n work(i) = a(i,k) a(i,k) = 0.0E+00 end do do j = k+1, n t = work(j) call saxpy ( n, t, a(1,j), 1, a(1,k), 1 ) end do l = ipvt(k) if ( l /= k ) then call sswap ( n, a(1,k), 1, a(1,l), 1 ) end if end do end if return end subroutine sgefa ( a, lda, n, ipvt, info ) ! !******************************************************************************* ! !! SGEFA factors a real matrix. ! ! ! Modified: ! ! 17 April 2002 ! ! Reference: ! ! Dongarra, Moler, Bunch and Stewart, ! LINPACK User's Guide, ! SIAM, (Society for Industrial and Applied Mathematics), ! 3600 University City Science Center, ! Philadelphia, PA, 19104-2688. ! ISBN 0-89871-172-X ! ! Parameters: ! ! Input/output, real A(LDA,N). ! On intput, the matrix to be factored. ! On output, an upper triangular matrix and the multipliers used to obtain ! it. The factorization can be written A=L*U, where L is a product of ! permutation and unit lower triangular matrices, and U is upper triangular. ! ! Input, integer LDA, the leading dimension of A. ! ! Input, integer N, the order of the matrix A. ! ! Output, integer IPVT(N), the pivot indices. ! ! Output, integer INFO, singularity indicator. ! 0, normal value. ! K, if U(K,K) == 0. This is not an error condition for this subroutine, ! but it does indicate that SGESL or SGEDI will divide by zero if called. ! Use RCOND in SGECO for a reliable indication of singularity. ! implicit none ! integer lda integer n ! real a(lda,n) integer info integer ipvt(n) integer isamax integer j integer k integer l real t ! ! Gaussian elimination with partial pivoting. ! info = 0 do k = 1, n - 1 ! ! Find L = pivot index. ! l = isamax ( n-k+1, a(k,k), 1 ) + k - 1 ipvt(k) = l ! ! Zero pivot implies this column already triangularized. ! if ( a(l,k) == 0.0E+00 ) then info = k cycle end if ! ! Interchange if necessary. ! if ( l /= k ) then call r_swap ( a(l,k), a(k,k) ) end if ! ! Compute multipliers. ! a(k+1:n,k) = - a(k+1:n,k) / a(k,k) ! ! Row elimination with column indexing. ! do j = k+1, n t = a(l,j) if ( l /= k ) then a(l,j) = a(k,j) a(k,j) = t end if call saxpy ( n-k, t, a(k+1,k), 1, a(k+1,j), 1 ) end do end do ipvt(n) = n if ( a(n,n) == 0.0E+00 ) then info = n end if return end subroutine sgesl ( a, lda, n, ipvt, b, job ) ! !******************************************************************************* ! !! SGESL solves a real general linear system A * X = B. ! ! ! Discussion: ! ! SGESL can solve either of the systems A * X = B or A' * X = B. ! ! The system matrix must have been factored by SGECO or SGEFA. ! ! A division by zero will occur if the input factor contains a ! zero on the diagonal. Technically this indicates singularity ! but it is often caused by improper arguments or improper ! setting of LDA. It will not occur if the subroutines are ! called correctly and if SGECO has set RCOND > 0.0E+00 ! or SGEFA has set INFO == 0. ! ! Reference: ! ! Dongarra, Moler, Bunch and Stewart, ! LINPACK User's Guide, ! SIAM, (Society for Industrial and Applied Mathematics), ! 3600 University City Science Center, ! Philadelphia, PA, 19104-2688. ! ISBN 0-89871-172-X ! ! Modified: ! ! 17 April 2002 ! ! Parameters: ! ! Input, real A(LDA,N), the output from SGECO or SGEFA. ! ! Input, integer LDA, the leading dimension of A. ! ! Input, integer N, the order of the matrix A. ! ! Input, integer IPVT(N), the pivot vector from SGECO or SGEFA. ! ! Input/output, real B(N). ! On input, the right hand side vector. ! On output, the solution vector. ! ! Input, integer JOB. ! 0, solve A * X = B; ! nonzero, solve A' * X = B. ! implicit none ! integer lda integer n ! real a(lda,n) real b(n) integer ipvt(n) integer job integer k integer l real t ! ! Solve A * X = B. ! if ( job == 0 ) then do k = 1, n-1 l = ipvt(k) t = b(l) if ( l /= k ) then b(l) = b(k) b(k) = t end if call saxpy ( n-k, t, a(k+1,k), 1, b(k+1), 1 ) end do do k = n, 1, -1 b(k) = b(k) / a(k,k) t = -b(k) call saxpy ( k-1, t, a(1,k), 1, b(1), 1 ) end do else ! ! Solve A' * X = B. ! do k = 1, n t = dot_product ( a(1:k-1,k), b(1:k-1) ) b(k) = ( b(k) - t ) / a(k,k) end do do k = n-1, 1, -1 b(k) = b(k) + dot_product ( a(k+1:n,k), b(k+1:n) ) l = ipvt(k) if ( l /= k ) then call r_swap ( b(l), b(k) ) end if end do end if return end subroutine sgtsl ( n, c, d, e, b, info ) ! !******************************************************************************* ! !! SGTSL solves a general tridiagonal linear system. ! ! ! Reference: ! ! Dongarra, Moler, Bunch and Stewart, ! LINPACK User's Guide, ! SIAM, (Society for Industrial and Applied Mathematics), ! 3600 University City Science Center, ! Philadelphia, PA, 19104-2688. ! ISBN 0-89871-172-X ! ! Modified: ! ! 31 October 2001 ! ! Parameters: ! ! Input, integer N, the order of the tridiagonal matrix. ! ! Input/output, real C(N), contains the subdiagonal of the tridiagonal ! matrix in entries C(2:N). On output, C is destroyed. ! ! Input/output, real D(N). On input, the diagonal of the matrix. ! On output, D is destroyed. ! ! Input/output, real E(N), contains the superdiagonal of the tridiagonal ! matrix in entries E(1:N-1). On output E is destroyed. ! ! Input/output, real B(N). On input, the right hand side. On output, ! the solution. ! ! Output, integer INFO, error flag. ! 0, normal value. ! K, the K-th element of the diagonal becomes exactly zero. The ! subroutine returns if this error condition is detected. ! implicit none ! integer n ! real b(n) real c(n) real d(n) real e(n) integer info integer k real t ! info = 0 c(1) = d(1) if ( n >= 2 ) then d(1) = e(1) e(1) = 0.0E+00 e(n) = 0.0E+00 do k = 1, n - 1 ! ! Find the larger of the two rows, and interchange if necessary. ! if ( abs ( c(k+1) ) >= abs ( c(k) ) ) then call r_swap ( c(k), c(k+1) ) call r_swap ( d(k), d(k+1) ) call r_swap ( e(k), e(k+1) ) call r_swap ( b(k), b(k+1) ) end if ! ! Fail if no nonzero pivot could be found. ! if ( c(k) == 0.0E+00 ) then info = k return end if ! ! Zero elements. ! t = -c(k+1) / c(k) c(k+1) = d(k+1) + t * d(k) d(k+1) = e(k+1) + t * e(k) e(k+1) = 0.0E+00 b(k+1) = b(k+1) + t * b(k) end do end if if ( c(n) == 0.0E+00 ) then info = n return end if ! ! Back solve. ! b(n) = b(n) / c(n) if ( n > 1 ) then b(n-1) = ( b(n-1) - d(n-1) * b(n) ) / c(n-1) do k = n-2, 1, -1 b(k) = ( b(k) - d(k) * b(k+1) - e(k) * b(k+2) ) / c(k) end do end if return end function snrm2 ( n, x, incx ) ! !******************************************************************************* ! !! SNRM2 computes the Euclidean norm of a vector. ! ! ! Discussion: ! ! The original SNRM2 algorithm is accurate but written in an ! obscure and obsolete format. This version goes for clarity. ! ! Modified: ! ! 01 June 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer N, the number of entries in the vector. ! ! Input, real X(*), the vector whose norm is to be computed. ! ! Input, integer INCX, the increment between successive entries of X. ! ! Output, real SNRM2, the Euclidean norm of X. ! implicit none ! integer i integer incx integer ix integer n real samax real snrm2 real stemp real x(*) real xmax ! if ( n <= 0 ) then snrm2 = 0.0E+00 else xmax = samax ( n, x, incx ) if ( xmax == 0.0E+00 ) then snrm2 = 0.0E+00 else if ( incx >= 0 ) then ix = 1 else ix = ( - n + 1 ) * incx + 1 end if stemp = 0.0E+00 do i = 1, n stemp = stemp + ( x(ix) / xmax )**2 ix = ix + incx end do snrm2 = xmax * sqrt ( stemp ) end if end if return end subroutine spbco ( abd, lda, n, m, rcond, z, info ) ! !******************************************************************************* ! !! SPBCO factors a real symmetric positive definite banded matrix. ! ! ! Discussion: ! ! SPBCO also estimates the condition of the matrix. ! ! If RCOND is not needed, SPBFA is slightly faster. ! ! To solve A*X = B, follow SPBCO by SPBSL. ! ! To compute inverse(A)*C, follow SPBCO by SPBSL. ! ! To compute determinant(A), follow SPBCO by SPBDI. ! ! Band storage: ! ! If A is a symmetric positive definite band matrix, the following ! program segment will set up the input. ! ! m = (band width above diagonal) ! do j = 1, n ! i1 = max (1, j-m) ! do i = i1, j ! k = i-j+m+1 ! abd(k,j) = a(i,j) ! end do ! end do ! ! This uses M + 1 rows of A, except for the M by M upper left triangle, ! which is ignored. ! ! For example, if the original matrix is ! ! 11 12 13 0 0 0 ! 12 22 23 24 0 0 ! 13 23 33 34 35 0 ! 0 24 34 44 45 46 ! 0 0 35 45 55 56 ! 0 0 0 46 56 66 ! ! then N = 6, M = 2 and ABD should contain ! ! * * 13 24 35 46 ! * 12 23 34 45 56 ! 11 22 33 44 55 66 ! ! Reference: ! ! Dongarra, Moler, Bunch and Stewart, ! LINPACK User's Guide, ! SIAM, (Society for Industrial and Applied Mathematics), ! 3600 University City Science Center, ! Philadelphia, PA, 19104-2688. ! ISBN 0-89871-172-X ! ! Parameters: ! ! Input/output, real ABD(LDA,N). On input, the matrix to be factored. ! The columns of the upper triangle are stored in the columns of ABD and the ! diagonals of the upper triangle are stored in the rows of ABD. ! On output, an upper triangular matrix R, stored in band form, so that ! A = R'*R. If INFO /= 0, the factorization is not complete. ! ! Input, integer LDA, the leading dimension of the array ABD. ! LDA must be >= M+1. ! ! Input, integer N, the order of the matrix. ! ! Input, integer M, the number of diagonals above the main diagonal. ! ! Output, real RCOND, an estimate of the reciprocal condition of A. ! For the system A*X = B, relative perturbations in A and B of size ! EPSILON may cause relative perturbations in X of size EPSILON/RCOND. ! If RCOND is so small that the logical expression ! 1.0E+00 + RCOND == 1.0E+00 ! is true, then A may be singular to working precision. In particular, ! RCOND is zero if exact singularity is detected or the estimate underflows. ! ! Output, real Z(N), a work vector whose contents are usually unimportant. ! If A is singular to working precision, then Z is an approximate null ! vector in the sense that ! norm(A*Z) = RCOND * norm(A) * norm(Z). ! If INFO /= 0, Z is unchanged. ! ! Output, integer INFO, error flag. ! 0, for normal return. ! K, signals an error condition. The leading minor of order K is not ! positive definite. ! implicit none ! integer lda integer n ! real abd(lda,n) real anorm real ek integer i integer info integer j integer j2 integer k integer kp1 integer l integer la integer lb integer lm integer m integer mu real rcond real s real sdot real sm real t real wk real wkm real ynorm real z(n) ! ! Find the norm of A. ! do j = 1, n l = min ( j, m+1 ) mu = max ( m+2-j, 1 ) z(j) = sum ( abs ( abd(mu:mu+l-1,j) ) ) k = j - l do i = mu, m k = k + 1 z(k) = z(k) + abs ( abd(i,j) ) end do end do anorm = maxval ( z(1:n) ) ! ! Factor. ! call spbfa ( abd, lda, n, m, info ) if ( info /= 0 ) then return end if ! ! RCOND = 1/(norm(A)*(estimate of norm(inverse(A)))). ! ! Estimate = norm(Z)/norm(Y) where A*Z = Y and A*Y = E. ! ! The components of E are chosen to cause maximum local ! growth in the elements of W where R'*W = E. ! ! The vectors are frequently rescaled to avoid overflow. ! ! Solve R'*W = E. ! ek = 1.0E+00 z(1:n) = 0.0E+00 do k = 1, n if ( z(k) /= 0.0E+00 ) then ek = sign ( ek, -z(k) ) end if if ( abs ( ek - z(k) ) > abd(m+1,k) ) then s = abd(m+1,k) / abs ( ek - z(k) ) z(1:n) = s * z(1:n) ek = s * ek end if wk = ek - z(k) wkm = -ek - z(k) s = abs ( wk ) sm = abs ( wkm ) wk = wk / abd(m+1,k) wkm = wkm / abd(m+1,k) kp1 = k + 1 j2 = min ( k+m, n ) i = m + 1 if ( k+1 <= j2 ) then do j = k+1, j2 i = i - 1 sm = sm + abs ( z(j) + wkm * abd(i,j) ) z(j) = z(j) + wk * abd(i,j) s = s + abs ( z(j) ) end do if ( s < sm ) then t = wkm - wk wk = wkm i = m + 1 do j = k+1, j2 i = i - 1 z(j) = z(j) + t * abd(i,j) end do end if end if z(k) = wk end do z(1:n) = z(1:n) / sum ( abs ( z(1:n) ) ) ! ! Solve R*Y = W. ! do k = n, 1, -1 if ( abs ( z(k) ) > abd(m+1,k) ) then s = abd(m+1,k) / abs ( z(k) ) z(1:n) = s * z(1:n) end if z(k) = z(k) / abd(m+1,k) lm = min ( k-1, m ) la = m + 1 - lm lb = k - lm t = -z(k) call saxpy ( lm, t, abd(la,k), 1, z(lb), 1 ) end do z(1:n) = z(1:n) / sum ( abs ( z(1:n) ) ) ynorm = 1.0E+00 ! ! Solve R'*V = Y. ! do k = 1, n lm = min ( k-1, m ) la = m + 1 - lm lb = k - lm z(k) = z(k) - dot_product ( abd(la:la+lm-1,k), z(lb:lb+lm-1) ) if ( abs ( z(k) ) > abd(m+1,k) ) then s = abd(m+1,k) / abs ( z(k) ) z(1:n) = s * z(1:n) ynorm = s * ynorm end if z(k) = z(k) / abd(m+1,k) end do s = 1.0E+00 / sum ( abs ( z(1:n) ) ) z(1:n) = s * z(1:n) ynorm = s * ynorm ! ! Solve R*Z = W. ! do k = n, 1, -1 if ( abs ( z(k) ) > abd(m+1,k) ) then s = abd(m+1,k) / abs ( z(k) ) z(1:n) = s * z(1:n) ynorm = s * ynorm end if z(k) = z(k) / abd(m+1,k) lm = min ( k-1, m ) la = m + 1 - lm lb = k - lm t = -z(k) call saxpy ( lm, t, abd(la,k), 1, z(lb), 1 ) end do ! ! Make ZNORM = 1.0. ! s = 1.0E+00 / sum ( abs ( z(1:n) ) ) z(1:n) = s * z(1:n) ynorm = s * ynorm if ( anorm /= 0.0E+00 ) then rcond = ynorm / anorm else rcond = 0.0E+00 end if return end subroutine spbdi ( abd, lda, n, m, det ) ! !******************************************************************************* ! !! SPBDI computes the determinant of a matrix factored by SPBCO or SPBFA. ! ! ! Discussion: ! ! If the inverse is needed, use SPBSL N times. ! ! Reference: ! ! Dongarra, Moler, Bunch and Stewart, ! LINPACK User's Guide, ! SIAM, (Society for Industrial and Applied Mathematics), ! 3600 University City Science Center, ! Philadelphia, PA, 19104-2688. ! ISBN 0-89871-172-X ! ! Parameters: ! ! Input, real ABD(LDA,N), the output from SPBCO or SPBFA. ! ! Input, integer LDA, the leading dimension of the array ABD. ! ! Input, integer N, the order of the matrix. ! ! Input, integer M, the number of diagonals above the main diagonal. ! ! Output, real DET(2), the determinant of the original matrix in the form ! determinant = DET(1) * 10.0**DET(2) ! with 1.0E+00 <= DET(1) < 10.0E+00 or DET(1) == 0.0E+00. ! implicit none ! integer lda integer n ! real abd(lda,n) real det(2) integer i integer m real s ! ! Compute the determinant. ! det(1) = 1.0E+00 det(2) = 0.0E+00 s = 10.0E+00 do i = 1, n det(1) = abd(m+1,i)**2 * det(1) if ( det(1) == 0.0E+00 ) then return end if do while ( det(1) < 1.0E+00 ) det(1) = s * det(1) det(2) = det(2) - 1.0E+00 end do do while ( det(1) >= s ) det(1) = det(1) / s det(2) = det(2) + 1.0E+00 end do end do return end subroutine spbfa ( abd, lda, n, m, info ) ! !******************************************************************************* ! !! SPBFA factors a real symmetric positive definite matrix stored in band form. ! ! ! Discussion: ! ! SPBFA is usually called by SPBCO, but it can be called ! directly with a saving in time if RCOND is not needed. ! ! Band storage: ! ! If A is a symmetric positive definite band matrix, ! the following program segment will set up the input. ! ! m = (band width above diagonal) ! do j = 1, n ! i1 = max (1, j-m) ! do i = i1, j ! k = i-j+m+1 ! abd(k,j) = a(i,j) ! end do ! end do ! ! Reference: ! ! Dongarra, Moler, Bunch and Stewart, ! LINPACK User's Guide, ! SIAM, (Society for Industrial and Applied Mathematics), ! 3600 University City Science Center, ! Philadelphia, PA, 19104-2688. ! ISBN 0-89871-172-X ! ! Parameters: ! ! Input/output, real ABD(LDA,N). On input, the matrix to be factored. ! The columns of the upper triangle are stored in the columns of ABD ! and the diagonals of the upper triangle are stored in the ! rows of ABD. On output, an upper triangular matrix R, stored in band ! form, so that A = R' * R. ! ! Input, integer LDA, the leading dimension of the array ABD. ! LDA must be >= M+1. ! ! Input, integer N, the order of the matrix. ! ! Input, integer M, the number of diagonals above the main diagonal. ! ! Output, integer INFO, error indicator. ! 0, for normal return. ! K, if the leading minor of order K is not positive definite. ! implicit none ! integer lda integer n ! real abd(lda,n) integer ik integer info integer j integer jk integer k integer m integer mu real sdot real s real t ! do j = 1, n s = 0.0E+00 ik = m + 1 jk = max (j-m,1) mu = max (m+2-j,1) do k = mu, m t = abd(k,j) - sdot ( k-mu, abd(ik,jk), 1, abd(mu,j), 1 ) t = t / abd(m+1,jk) abd(k,j) = t s = s + t * t ik = ik - 1 jk = jk + 1 end do s = abd(m+1,j) - s if ( s <= 0.0E+00 ) then info = j return end if abd(m+1,j) = sqrt ( s ) end do info = 0 return end subroutine spbsl ( abd, lda, n, m, b ) ! !******************************************************************************* ! !! SPBSL solves a real symmetric positive definite band factored by SPBCO or SPBFA. ! ! ! Discussion: ! ! To compute inverse(A) * C where C is a matrix with P columns: ! ! call spbco ( abd, lda, n, rcond, z, info ) ! ! if ( rcond is too small .or. info /= 0) go to ... ! ! do j = 1, p ! call spbsl ( abd, lda, n, c(1,j) ) ! end do ! ! A division by zero will occur if the input factor contains ! a zero on the diagonal. Technically this indicates ! singularity but it is usually caused by improper subroutine ! arguments. It will not occur if the subroutines are called ! correctly and INFO == 0. ! ! Reference: ! ! Dongarra, Moler, Bunch and Stewart, ! LINPACK User's Guide, ! SIAM, (Society for Industrial and Applied Mathematics), ! 3600 University City Science Center, ! Philadelphia, PA, 19104-2688. ! ISBN 0-89871-172-X ! ! Parameters: ! ! Input, real ABD(LDA,N), the output from SPBCO or SPBFA. ! ! Input, integer LDA, the leading dimension of the array ABD. ! ! Input, integer N, the order of the matrix. ! ! Input, integer M, the number of diagonals above the main diagonal. ! ! Input/output, real B(N). On input, the right hand side. ! On output, the solution. ! implicit none ! integer lda integer n ! real abd(lda,n) real b(n) integer k integer la integer lb integer lm integer m real t ! ! Solve R'*Y = B. ! do k = 1, n lm = min ( k-1, m ) t = dot_product ( abd(m+1-lm:m,k), b(k-lm:k-1) ) b(k) = ( b(k) - t ) / abd(m+1,k) end do ! ! Solve R*X = Y. ! do k = n, 1, -1 lm = min ( k-1, m ) la = m + 1 - lm lb = k - lm b(k) = b(k) / abd(m+1,k) t = -b(k) call saxpy ( lm, t, abd(la,k), 1, b(lb), 1 ) end do return end subroutine spoco ( a, lda, n, rcond, z, info ) ! !******************************************************************************* ! !! SPOCO factors a real symmetric positive definite matrix and estimates its condition. ! ! ! Discussion: ! ! If RCOND is not needed, SPOFA is slightly faster. ! ! To solve A*X = B, follow SPOCO by SPOSL. ! ! To compute inverse(A)*C, follow SPOCO by SPOSL. ! ! To compute determinant(A), follow SPOCO by SPODI. ! ! To compute inverse(A), follow SPOCO by SPODI. ! ! Reference: ! ! Dongarra, Moler, Bunch and Stewart, ! LINPACK User's Guide, ! SIAM, (Society for Industrial and Applied Mathematics), ! 3600 University City Science Center, ! Philadelphia, PA, 19104-2688. ! ISBN 0-89871-172-X ! ! Parameters: ! ! Input/output, real A(LDA,N). On input, the symmetric matrix to ! be factored. Only the diagonal and upper triangle are used. ! On output, an upper triangular matrix R so that A = R'*R where R' ! is the transpose. The strict lower triangle is unaltered. ! If INFO /= 0, the factorization is not complete. ! ! Input, integer LDA, the leading dimension of the array A. ! ! Input, integer N, the order of the matrix. ! ! Output, real RCOND, an estimate of the reciprocal condition of A. ! For the system A*X = B, relative perturbations in A and B of size ! EPSILON may cause relative perturbations in X of size EPSILON/RCOND. ! If RCOND is so small that the logical expression ! 1.0E+00 + RCOND == 1.0E+00 ! is true, then A may be singular to working precision. In particular, ! RCOND is zero if exact singularity is detected or the estimate underflows. ! ! Output, real Z(N), a work vector whose contents are usually unimportant. ! If A is close to a singular matrix, then Z is an approximate null vector ! in the sense that ! norm(A*Z) = RCOND * norm(A) * norm(Z). ! If INFO /= 0, Z is unchanged. ! ! Output, integer INFO, error flag. ! 0, for normal return. ! K, signals an error condition. The leading minor of order K is not ! positive definite. ! implicit none ! integer lda integer n real a(lda,n) real anorm real ek integer i integer info integer j integer k integer kp1 real rcond real s real sm real t real wk real wkm real ynorm real z(n) ! ! Find norm of A using only upper half. ! do j = 1, n z(j) = sum ( abs ( a(1:j,j) ) ) do i = 1, j-1 z(i) = z(i) + abs ( a(i,j) ) end do end do anorm = maxval ( z(1:n) ) ! ! Factor. ! call spofa ( a, lda, n, info ) if ( info /= 0 ) then return end if ! ! RCOND = 1/(norm(A)*(estimate of norm(inverse(A)))). ! ! Estimate = norm(Z)/norm(Y) where A*Z = Y and A*Y = E. ! ! The components of E are chosen to cause maximum local ! growth in the elements of W where R'*W = E. ! ! The vectors are frequently rescaled to avoid overflow. ! ! Solve R'*W = E. ! ek = 1.0E+00 z(1:n) = 0.0E+00 do k = 1, n if ( z(k) /= 0.0E+00 ) then ek = sign ( ek, -z(k) ) end if if ( abs ( ek - z(k) ) > a(k,k) ) then s = a(k,k) / abs ( ek - z(k) ) z(1:n) = s * z(1:n) ek = s * ek end if wk = ek - z(k) wkm = -ek - z(k) s = abs ( wk ) sm = abs ( wkm ) wk = wk / a(k,k) wkm = wkm / a(k,k) kp1 = k + 1 if ( k+1 <= n ) then do j = k+1, n sm = sm + abs ( z(j) + wkm * a(k,j) ) z(j) = z(j) + wk * a(k,j) s = s + abs ( z(j) ) end do if ( s < sm ) then t = wkm - wk wk = wkm do j = k+1, n z(j) = z(j) + t * a(k,j) end do end if end if z(k) = wk end do z(1:n) = z(1:n) / sum ( abs ( z(1:n) ) ) ! ! Solve R*Y = W. ! do k = n, 1, -1 if ( abs ( z(k) ) > a(k,k) ) then s = a(k,k) / abs ( z(k) ) z(1:n) = s * z(1:n) end if z(k) = z(k)/a(k,k) t = -z(k) call saxpy ( k-1, t, a(1,k), 1, z(1), 1 ) end do z(1:n) = z(1:n) / sum ( abs ( z(1:n) ) ) ynorm = 1.0E+00 ! ! Solve R'*V = Y. ! do k = 1, n z(k) = z(k) - dot_product ( a(1:k-1,k), z(1:k-1) ) if ( abs ( z(k) ) > a(k,k) ) then s = a(k,k) / abs ( z(k) ) z(1:n) = s * z(1:n) ynorm = s * ynorm end if z(k) = z(k) / a(k,k) end do s = 1.0E+00 / sum ( abs ( z(1:n) ) ) z(1:n) = s * z(1:n) ynorm = s * ynorm ! ! Solve R*Z = V. ! do k = n, 1, -1 if ( abs ( z(k) ) > a(k,k) ) then s = a(k,k) / abs ( z(k) ) z(1:n) = s * z(1:n) ynorm = s * ynorm end if z(k) = z(k)/a(k,k) t = -z(k) call saxpy ( k-1, t, a(1,k), 1, z(1), 1 ) end do ! ! Make ZNORM = 1.0. ! s = 1.0E+00 / sum ( abs ( z(1:n) ) ) z(1:n) = s * z(1:n) ynorm = s * ynorm if ( anorm /= 0.0E+00 ) then rcond = ynorm / anorm else rcond = 0.0E+00 end if return end subroutine spodi ( a, lda, n, det, job ) ! !******************************************************************************* ! !! SPODI computes the determinant and inverse of a certain matrix. ! ! ! Discussion: ! ! The matrix is real symmetric positive definite. ! SPODI uses the factors computed by SPOCO, SPOFA or SQRDC. ! ! A division by zero will occur if the input factor contains ! a zero on the diagonal and the inverse is requested. ! It will not occur if the subroutines are called correctly ! and if SPOCO or SPOFA has set INFO == 0. ! ! Reference: ! ! Dongarra, Moler, Bunch and Stewart, ! LINPACK User's Guide, ! SIAM, (Society for Industrial and Applied Mathematics), ! 3600 University City Science Center, ! Philadelphia, PA, 19104-2688. ! ISBN 0-89871-172-X ! ! Parameters: ! ! Input/output, real A(LDA,N). On input, the output A from SPOCO ! or SPOFA, or the output X from SQRDC. On output, if SPOCO or ! SPOFA was used to factor A then SPODI produces the upper half of ! inverse(A). If SQRDC was used to decompose X then SPODI produces ! the upper half of inverse(X'*X) where X' is the transpose. ! Elements of A below the diagonal are unchanged. If the units digit ! of JOB is zero, A is unchanged. ! ! Input, integer LDA, the leading dimension of the array A. ! ! Input, integer N, the order of the matrix A. ! ! Input, integer JOB, specifies the task. ! 11, both determinant and inverse. ! 01, inverse only. ! 10, determinant only. ! ! Output, real DET(2), the determinant of A or of X'*X if requested. ! determinant = DET(1) * 10.0**DET(2) ! with 1.0E+00 <= DET(1) < 10.0E+00 or DET(1) == 0.0E+00. ! implicit none ! integer lda integer n ! real a(lda,n) real det(2) integer i integer j integer job integer k integer kp1 real s real t ! ! Compute the determinant. ! if ( job / 10 /= 0 ) then det(1) = 1.0E+00 det(2) = 0.0E+00 s = 10.0E+00 do i = 1, n det(1) = a(i,i)**2 * det(1) if ( det(1) == 0.0E+00 ) then exit end if do while ( det(1) < 1.0E+00 ) det(1) = s * det(1) det(2) = det(2) - 1.0E+00 end do do while ( det(1) >= s ) det(1) = det(1) / s det(2) = det(2) + 1.0E+00 end do end do end if ! ! Compute inverse(R). ! if ( mod ( job, 10 ) /= 0 ) then do k = 1, n a(k,k) = 1.0E+00 / a(k,k) t = -a(k,k) call sscal ( k-1, t, a(1,k), 1 ) kp1 = k + 1 do j = k+1, n t = a(k,j) a(k,j) = 0.0E+00 call saxpy ( k, t, a(1,k), 1, a(1,j), 1 ) end do end do ! ! Form inverse(R) * transpose(inverse(R)). ! do j = 1, n do k = 1, j-1 t = a(k,j) call saxpy ( k, t, a(1,j), 1, a(1,k), 1 ) end do t = a(j,j) call sscal ( j, t, a(1,j), 1 ) end do end if return end subroutine spofa ( a, lda, n, info ) ! !******************************************************************************* ! !! SPOFA factors a real symmetric positive definite matrix. ! ! ! Discussion: ! ! SPOFA is usually called by SPOCO, but it can be called ! directly with a saving in time if RCOND is not needed. ! ! Modified: ! ! 17 April 2002 ! ! Reference: ! ! Dongarra, Moler, Bunch and Stewart, ! LINPACK User's Guide, ! SIAM, (Society for Industrial and Applied Mathematics), ! 3600 University City Science Center, ! Philadelphia, PA, 19104-2688. ! ISBN 0-89871-172-X ! ! Parameters: ! ! Input/output, real A(LDA,N). On input, the symmetric matrix to be ! factored. Only the diagonal and upper triangle are used. ! On output, an upper triangular matrix R so that A = R'*R ! where R' is the transpose. The strict lower triangle is unaltered. ! If INFO /= 0, the factorization is not complete. ! ! Input, integer LDA, the leading dimension of the array A. ! ! Input, integer N, the order of the matrix. ! ! Output, integer INFO, error flag. ! 0, for normal return. ! K, signals an error condition. The leading minor of order K is not ! positive definite. ! implicit none ! integer lda integer n ! real a(lda,n) integer info integer j integer k real s real t ! do j = 1, n s = 0.0E+00 do k = 1, j-1 t = a(k,j) - dot_product ( a(1:k-1,k), a(1:k-1,j) ) t = t / a(k,k) a(k,j) = t s = s + t * t end do s = a(j,j) - s if ( s <= 0.0E+00 ) then info = j return end if a(j,j) = sqrt ( s ) end do info = 0 return end subroutine sposl ( a, lda, n, b ) ! !******************************************************************************* ! !! SPOSL solves a linear system factored by SPOCO or SPOFA. ! ! ! Discussion: ! ! To compute inverse(A) * C where C is a matrix with P columns: ! ! call spoco ( a, lda, n, rcond, z, info ) ! ! if ( rcond is not too small .and. info == 0 ) then ! do j = 1, p ! call sposl ( a, lda, n, c(1,j) ) ! end do ! end if ! ! A division by zero will occur if the input factor contains ! a zero on the diagonal. Technically this indicates ! singularity but it is usually caused by improper subroutine ! arguments. It will not occur if the subroutines are called ! correctly and INFO == 0. ! ! Reference: ! ! Dongarra, Moler, Bunch and Stewart, ! LINPACK User's Guide, ! SIAM, (Society for Industrial and Applied Mathematics), ! 3600 University City Science Center, ! Philadelphia, PA, 19104-2688. ! ISBN 0-89871-172-X ! ! Parameters: ! ! Input, real A(LDA,N), the output from SPOCO or SPOFA. ! ! Input, integer LDA, the leading dimension of the array A. ! ! Input, integer N, the order of the matrix. ! ! Input/output, real B(N). On input, the right hand side. ! On output, the solution. ! implicit none ! integer lda integer n ! real a(lda,n) real b(n) integer k real t ! ! Solve R'*Y = B. ! do k = 1, n t = dot_product ( a(1:k-1,k), b(1:k-1) ) b(k) = ( b(k) - t ) / a(k,k) end do ! ! Solve R*X = Y. ! do k = n, 1, -1 b(k) = b(k)/a(k,k) t = -b(k) call saxpy ( k-1, t, a(1,k), 1, b(1), 1 ) end do return end subroutine sppco ( ap, n, rcond, z, info ) ! !******************************************************************************* ! !! SPPCO factors a real symmetric positive definite matrix in packed form. ! ! ! Discussion: ! ! SPPCO also estimates the condition of the matrix. ! ! If RCOND is not needed, SPPFA is slightly faster. ! ! To solve A*X = B, follow SPPCO by SPPSL. ! ! To compute inverse(A)*C, follow SPPCO by SPPSL. ! ! To compute determinant(A), follow SPPCO by SPPDI. ! ! To compute inverse(A), follow SPPCO by SPPDI. ! ! Packed storage: ! ! The following program segment will pack the upper triangle of ! a symmetric matrix. ! ! k = 0 ! do j = 1, n ! do i = 1, j ! k = k + 1 ! ap(k) = a(i,j) ! end do ! end do ! ! Reference: ! ! Dongarra, Moler, Bunch and Stewart, ! LINPACK User's Guide, ! SIAM, (Society for Industrial and Applied Mathematics), ! 3600 University City Science Center, ! Philadelphia, PA, 19104-2688. ! ISBN 0-89871-172-X ! ! Parameters: ! ! Input/output, real AP(N*(N+1)/2). On input, the packed form of a ! symmetric matrix A. The columns of the upper triangle are stored ! sequentially in a one-dimensional array. On output, an upper ! triangular matrix R, stored in packed form, so that A = R'*R. ! If INFO /= 0, the factorization is not complete. ! ! Input, integer N, the order of the matrix. ! ! Output, real RCOND, an estimate of the reciprocal condition of A. ! For the system A*X = B, relative perturbations in A and B of size ! EPSILON may cause relative perturbations in X of size EPSILON/RCOND. ! If RCOND is so small that the logical expression ! 1.0E+00 + RCOND == 1.0E+00 ! is true, then A may be singular to working precision. In particular, ! RCOND is zero if exact singularity is detected or the estimate underflows. ! ! Output, real Z(N), a work vector whose contents are usually unimportant. ! If A is singular to working precision, then Z is an approximate null ! vector in the sense that ! norm(A*Z) = RCOND * norm(A) * norm(Z). ! If INFO /= 0, Z is unchanged. ! ! Output, integer INFO, error flag. ! 0, for normal return. ! K, signals an error condition. The leading minor of order K is ! not positive definite. ! implicit none ! integer n ! real anorm real ap((n*(n+1))/2) real ek integer i integer ij integer info integer j integer j1 integer k integer kj integer kk integer kp1 real rcond real s real sdot real sm real t real wk real wkm real ynorm real z(n) ! ! Find the norm of A. ! j1 = 1 do j = 1, n z(j) = sum ( abs ( ap(j1:j1+j-1) ) ) ij = j1 j1 = j1 + j do i = 1, j-1 z(i) = z(i) + abs ( ap(ij) ) ij = ij + 1 end do end do anorm = maxval ( z(1:n) ) ! ! Factor. ! call sppfa ( ap, n, info ) if ( info /= 0 ) then return end if ! ! RCOND = 1/(norm(A)*(estimate of norm(inverse(A)))). ! ! Estimate = norm(Z)/norm(Y) where A*Z = Y and A*Y = E. ! ! The components of E are chosen to cause maximum local ! growth in the elements of W where R'*W = E. ! ! The vectors are frequently rescaled to avoid overflow. ! ! Solve R'*W = E. ! ek = 1.0E+00 z(1:n) = 0.0E+00 kk = 0 do k = 1, n kk = kk + k if ( z(k) /= 0.0E+00 ) then ek = sign ( ek, -z(k) ) end if if ( abs ( ek - z(k) ) > ap(kk) ) then s = ap(kk) / abs ( ek - z(k) ) z(1:n) = s * z(1:n) ek = s * ek end if wk = ek - z(k) wkm = -ek - z(k) s = abs ( wk ) sm = abs ( wkm ) wk = wk / ap(kk) wkm = wkm / ap(kk) kp1 = k + 1 kj = kk + k if ( k+1 <= n ) then do j = k+1, n sm = sm + abs ( z(j) + wkm * ap(kj) ) z(j) = z(j) + wk * ap(kj) s = s + abs ( z(j) ) kj = kj + j end do if ( s < sm ) then t = wkm - wk wk = wkm kj = kk + k do j = k+1, n z(j) = z(j) + t * ap(kj) kj = kj + j end do end if end if z(k) = wk end do z(1:n) = z(1:n) / sum ( abs ( z(1:n) ) ) ! ! Solve R*Y = W. ! do k = n, 1, -1 if ( abs ( z(k) ) > ap(kk) ) then s = ap(kk) / abs ( z(k) ) z(1:n) = s * z(1:n) end if z(k) = z(k) / ap(kk) kk = kk - k t = -z(k) call saxpy ( k-1, t, ap(kk+1), 1, z(1), 1 ) end do z(1:n) = z(1:n) / sum ( abs ( z(1:n) ) ) ynorm = 1.0E+00 ! ! Solve R'*V = Y. ! do k = 1, n z(k) = z(k) - sdot ( k-1, ap(kk+1), 1, z(1), 1 ) kk = kk + k if ( abs ( z(k) ) > ap(kk) ) then s = ap(kk) / abs ( z(k) ) z(1:n) = s * z(1:n) ynorm = s * ynorm end if z(k) = z(k) / ap(kk) end do s = 1.0E+00 / sum ( abs ( z(1:n) ) ) z(1:n) = s * z(1:n) ynorm = s * ynorm ! ! Solve R*Z = V. ! do k = n, 1, -1 if ( abs ( z(k) ) > ap(kk) ) then s = ap(kk) / abs ( z(k) ) z(1:n) = s * z(1:n) ynorm = s * ynorm end if z(k) = z(k) / ap(kk) kk = kk - k t = -z(k) call saxpy ( k-1, t, ap(kk+1), 1, z(1), 1 ) end do ! ! Make ZNORM = 1.0. ! s = 1.0E+00 / sum ( abs ( z(1:n) ) ) z(1:n) = s * z(1:n) ynorm = s * ynorm if ( anorm /= 0.0E+00 ) then rcond = ynorm / anorm else rcond = 0.0E+00 end if return end subroutine sppdi ( ap, n, det, job ) ! !******************************************************************************* ! !! SPPDI computes the determinant and inverse of a matrix factored by SPPCO or SPPFA. ! ! ! Discussion: ! ! A division by zero will occur if the input factor contains ! a zero on the diagonal and the inverse is requested. ! It will not occur if the subroutines are called correctly ! and if SPOCO or SPOFA has set INFO == 0. ! ! Reference: ! ! Dongarra, Moler, Bunch and Stewart, ! LINPACK User's Guide, ! SIAM, (Society for Industrial and Applied Mathematics), ! 3600 University City Science Center, ! Philadelphia, PA, 19104-2688. ! ISBN 0-89871-172-X ! ! Parameters: ! ! Input/output, real AP(N*(N+1)/2). On input, the output from ! SPPCO or SPPFA. On output, the upper triangular half of the ! inverse, if requested. ! ! Input, integer N, the order of the matrix. ! ! Output, real DET(2), the determinant of the original matrix if requested. ! determinant = DET(1) * 10.0**DET(2) ! with 1.0E+00 <= DET(1) < 10.0E+00 or DET(1) == 0.0E+00. ! ! Input, integer JOB, job request. ! 11, both determinant and inverse. ! 01, inverse only. ! 10, determinant only. ! implicit none ! integer n ! real ap((n*(n+1))/2) real det(2) integer i integer ii integer j integer j1 integer jj integer job integer k integer k1 integer kj integer kk integer kp1 real s real t ! ! Compute the determinant. ! if ( job/10 /= 0 ) then det(1) = 1.0E+00 det(2) = 0.0E+00 s = 10.0E+00 ii = 0 do i = 1, n ii = ii + i det(1) = ap(ii)**2 * det(1) if ( det(1) == 0.0E+00 ) then exit end if do while ( det(1) < 1.0E+00 ) det(1) = s * det(1) det(2) = det(2) - 1.0E+00 end do do while ( det(1) >= s ) det(1) = det(1) / s det(2) = det(2) + 1.0E+00 end do end do end if ! ! Compute inverse(R). ! if ( mod ( job, 10 ) /= 0 ) then kk = 0 do k = 1, n k1 = kk + 1 kk = kk + k ap(kk) = 1.0E+00 / ap(kk) t = - ap(kk) call sscal ( k-1, t, ap(k1), 1 ) kp1 = k + 1 j1 = kk + 1 kj = kk + k do j = k+1, n t = ap(kj) ap(kj) = 0.0E+00 call saxpy ( k, t, ap(k1), 1, ap(j1), 1 ) j1 = j1 + j kj = kj + j end do end do ! ! Form inverse(R) * transpose(inverse(R)). ! jj = 0 do j = 1, n j1 = jj + 1 jj = jj + j k1 = 1 kj = j1 do k = 1, j-1 t = ap(kj) call saxpy ( k, t, ap(j1), 1, ap(k1), 1 ) k1 = k1 + k kj = kj + 1 end do t = ap(jj) call sscal ( j, t, ap(j1), 1 ) end do end if return end subroutine sppfa ( ap, n, info ) ! !******************************************************************************* ! !! SPPFA factors a real symmetric positive definite matrix in packed form. ! ! ! Discussion: ! ! SPPFA is usually called by SPPCO, but it can be called ! directly with a saving in time if RCOND is not needed. ! ! Packed storage: ! ! The following program segment will pack the upper ! triangle of a symmetric matrix. ! ! k = 0 ! do j = 1, n ! do i = 1, j ! k = k + 1 ! ap(k) = a(i,j) ! end do ! end do ! ! Reference: ! ! Dongarra, Moler, Bunch and Stewart, ! LINPACK User's Guide, ! SIAM, (Society for Industrial and Applied Mathematics), ! 3600 University City Science Center, ! Philadelphia, PA, 19104-2688. ! ISBN 0-89871-172-X ! ! Parameters: ! ! Input/output, real AP(N*(N+1)/2). On input, the packed form of a ! symmetric matrix A. The columns of the upper triangle are stored ! sequentially in a one-dimensional array. On output, an upper ! triangular matrix R, stored in packed form, so that A = R'*R. ! ! Input, integer N, the order of the matrix. ! ! Output, integer INFO, error flag. ! 0, for normal return. ! K, if the leading minor of order K is not positive definite. ! implicit none ! integer n ! real ap((n*(n+1))/2) integer info integer j integer jj integer k integer kj integer kk real s real sdot real t ! info = 0 jj = 0 do j = 1, n s = 0.0E+00 kj = jj kk = 0 do k = 1, j-1 kj = kj + 1 t = ap(kj) - sdot ( k-1, ap(kk+1), 1, ap(jj+1), 1 ) kk = kk + k t = t / ap(kk) ap(kj) = t s = s + t * t end do jj = jj + j s = ap(jj) - s if ( s <= 0.0E+00 ) then info = j return end if ap(jj) = sqrt ( s ) end do return end subroutine sppsl ( ap, n, b ) ! !******************************************************************************* ! !! SPPSL solves a real symmetric positive definite system factored by SPPCO or SPPFA. ! ! ! Discussion: ! ! To compute inverse(A) * C where C is a matrix with P columns ! ! call sppco ( ap, n, rcond, z, info ) ! ! if ( rcond is too small .or. info /= 0 ) then ! exit ! end if ! ! do j = 1, p ! call sppsl ( ap, n, c(1,j) ) ! end do ! ! A division by zero will occur if the input factor contains ! a zero on the diagonal. Technically this indicates ! singularity but it is usually caused by improper subroutine ! arguments. It will not occur if the subroutines are called ! correctly and INFO == 0. ! ! Reference: ! ! Dongarra, Moler, Bunch and Stewart, ! LINPACK User's Guide, ! SIAM, (Society for Industrial and Applied Mathematics), ! 3600 University City Science Center, ! Philadelphia, PA, 19104-2688. ! ISBN 0-89871-172-X ! ! Parameters: ! ! Input, real AP(N*(N+1)/2), the output from SPPCO or SPPFA. ! ! Input, integer N, the order of the matrix. ! ! Input/output, real B(N). On input, the right hand side. ! On output, the solution. ! implicit none ! integer n ! real ap((n*(n+1))/2) real b(n) integer k integer kk real sdot real t ! kk = 0 do k = 1, n t = sdot ( k-1, ap(kk+1), 1, b(1), 1 ) kk = kk + k b(k) = (b(k) - t)/ap(kk) end do do k = n, 1, -1 b(k) = b(k) / ap(kk) kk = kk - k t = - b(k) call saxpy ( k-1, t, ap(kk+1), 1, b(1), 1 ) end do return end subroutine sptsl ( n, d, e, b ) ! !******************************************************************************* ! !! SPTSL solves a positive definite tridiagonal linear system. ! ! ! Reference: ! ! Dongarra, Moler, Bunch and Stewart, ! LINPACK User's Guide, ! SIAM, (Society for Industrial and Applied Mathematics), ! 3600 University City Science Center, ! Philadelphia, PA, 19104-2688. ! ISBN 0-89871-172-X ! ! Parameters: ! ! Input, integer N, the order of the matrix. ! ! Input/output, real D(N), on input the diagonal of the tridiagonal matrix. ! On output, D is destroyed. ! ! Input, real E(N), the offdiagonal of the tridiagonal matrix in ! entries E(1:N-1). ! ! Input/output, real B(N). On input, the right hand side. ! On output, the solution. ! implicit none ! integer n ! real b(n) real d(n) real e(n) integer k integer kbm1 integer ke integer kf integer kp1 integer nm1d2 real t1 real t2 ! ! Check for 1 x 1 case. ! if ( n == 1 ) then b(1) = b(1) / d(1) return end if nm1d2 = ( n - 1 ) / 2 if ( n > 2 ) then kbm1 = n - 1 ! ! Zero top half of subdiagonal and bottom half of superdiagonal. ! do k = 1, nm1d2 t1 = e(k) / d(k) d(k+1) = d(k+1) - t1 * e(k) b(k+1) = b(k+1) - t1 * b(k) t2 = e(kbm1) / d(kbm1+1) d(kbm1) = d(kbm1) - t2 * e(kbm1) b(kbm1) = b(kbm1) - t2 * b(kbm1+1) kbm1 = kbm1 - 1 end do end if kp1 = nm1d2 + 1 ! ! Clean up for possible 2 x 2 block at center. ! if ( mod ( n, 2 ) == 0 ) then t1 = e(kp1) / d(kp1) d(kp1+1) = d(kp1+1) - t1 * e(kp1) b(kp1+1) = b(kp1+1) - t1 * b(kp1) kp1 = kp1 + 1 end if ! ! Back solve starting at the center, going towards the top and bottom. ! b(kp1) = b(kp1)/d(kp1) if ( n > 2 ) then k = kp1 - 1 ke = kp1 + nm1d2 - 1 do kf = kp1, ke b(k) = ( b(k) - e(k) * b(k+1) ) / d(k) b(kf+1) = ( b(kf+1) - e(kf) * b(kf) ) / d(kf+1) k = k - 1 end do end if if ( mod ( n, 2 ) == 0 ) then b(1) = ( b(1) - e(1) * b(2) ) / d(1) end if return end subroutine sqrdc ( a, lda, n, p, qraux, jpvt, work, job ) ! !******************************************************************************* ! !! SQRDC computes the QR factorization of a real rectangular matrix. ! ! ! Discussion: ! ! SQRDC uses Householder transformations. ! ! Column pivoting based on the 2-norms of the reduced columns may be ! performed at the user's option. ! ! Reference: ! ! Dongarra, Moler, Bunch and Stewart, ! LINPACK User's Guide, ! SIAM, (Society for Industrial and Applied Mathematics), ! 3600 University City Science Center, ! Philadelphia, PA, 19104-2688. ! ISBN 0-89871-172-X ! ! Parameters: ! ! Input/output, real A(LDA,P). On input, the N by P matrix ! whose decomposition is to be computed. On output, A contains in ! its upper triangle the upper triangular matrix R of the QR ! factorization. Below its diagonal A contains information from ! which the orthogonal part of the decomposition can be recovered. ! Note that if pivoting has been requested, the decomposition is not that ! of the original matrix A but that of A with its columns permuted ! as described by JPVT. ! ! Input, integer LDA, the leading dimension of the array A. LDA must ! be at least N. ! ! Input, integer N, the number of rows of the matrix A. ! ! Input, integer P, the number of columns of the matrix A. ! ! Output, real QRAUX(P), contains further information required to recover ! the orthogonal part of the decomposition. ! ! Input/output, integer JPVT(P). On input, JPVT contains integers that ! control the selection of the pivot columns. The K-th column A(*,K) of A ! is placed in one of three classes according to the value of JPVT(K). ! > 0, then A(K) is an initial column. ! = 0, then A(K) is a free column. ! < 0, then A(K) is a final column. ! Before the decomposition is computed, initial columns are moved to ! the beginning of the array A and final columns to the end. Both ! initial and final columns are frozen in place during the computation ! and only free columns are moved. At the K-th stage of the ! reduction, if A(*,K) is occupied by a free column it is interchanged ! with the free column of largest reduced norm. JPVT is not referenced ! if JOB == 0. On output, JPVT(K) contains the index of the column of the ! original matrix that has been interchanged into the K-th column, if ! pivoting was requested. ! ! Workspace, real WORK(P). WORK is not referenced if JOB == 0. ! ! Input, integer JOB, initiates column pivoting. ! 0, no pivoting is done. ! nonzero, pivoting is done. ! implicit none ! integer lda integer n integer p ! real a(lda,p) integer jpvt(p) real qraux(p) real work(p) integer j integer jj integer job integer jp integer l integer lup integer maxj real maxnrm real nrmxl integer pl integer pu real sdot real snrm2 logical swapj real t real tt ! pl = 1 pu = 0 ! ! If pivoting is requested, rearrange the columns. ! if ( job /= 0 ) then do j = 1, p swapj = jpvt(j) > 0 if ( jpvt(j) < 0 ) then jpvt(j) = - j else jpvt(j) = j end if if ( swapj ) then if ( j /= pl ) then call sswap ( n, a(1,pl), 1, a(1,j), 1 ) end if jpvt(j) = jpvt(pl) jpvt(pl) = j pl = pl + 1 end if end do pu = p do j = p, 1, -1 if ( jpvt(j) < 0 ) then jpvt(j) = - jpvt(j) if ( j /= pu ) then call sswap ( n, a(1,pu), 1, a(1,j), 1 ) call i_swap ( jpvt(pu), jpvt(j) ) end if pu = pu - 1 end if end do end if ! ! Compute the norms of the free columns. ! do j = pl, pu qraux(j) = snrm2 ( n, a(1,j), 1 ) end do work(pl:pu) = qraux(pl:pu) ! ! Perform the Householder reduction of A. ! lup = min ( n, p ) do l = 1, lup ! ! Bring the column of largest norm into the pivot position. ! if ( l >= pl .and. l < pu ) then maxnrm = 0.0E+00 maxj = l do j = l, pu if ( qraux(j) > maxnrm ) then maxnrm = qraux(j) maxj = j end if end do if ( maxj /= l ) then call sswap ( n, a(1,l), 1, a(1,maxj), 1 ) qraux(maxj) = qraux(l) work(maxj) = work(l) call i_swap ( jpvt(maxj), jpvt(l) ) end if end if ! ! Compute the Householder transformation for column L. ! qraux(l) = 0.0E+00 if ( l /= n ) then nrmxl = snrm2 ( n-l+1, a(l,l), 1 ) if ( nrmxl /= 0.0E+00 ) then if ( a(l,l) /= 0.0E+00 ) then nrmxl = sign ( nrmxl, a(l,l) ) end if call sscal ( n-l+1, 1.0E+00 / nrmxl, a(l,l), 1 ) a(l,l) = 1.0E+00 + a(l,l) ! ! Apply the transformation to the remaining columns, updating the norms. ! do j = l + 1, p t = - sdot ( n-l+1, a(l,l), 1, a(l,j), 1 ) / a(l,l) call saxpy ( n-l+1, t, a(l,l), 1, a(l,j), 1 ) if ( j >= pl .and. j <= pu ) then if ( qraux(j) /= 0.0E+00 ) then tt = 1.0E+00 - ( abs ( a(l,j) ) / qraux(j) )**2 tt = max ( tt, 0.0E+00 ) t = tt tt = 1.0E+00 + 0.05E+00 * tt * ( qraux(j) / work(j) )**2 if ( tt /= 1.0E+00 ) then qraux(j) = qraux(j) * sqrt ( t ) else qraux(j) = snrm2 ( n-l, a(l+1,j), 1 ) work(j) = qraux(j) end if end if end if end do ! ! Save the transformation. ! qraux(l) = a(l,l) a(l,l) = - nrmxl end if end if end do return end subroutine sqrsl ( a, lda, n, k, qraux, y, qy, qty, b, rsd, ab, job, info ) ! !******************************************************************************* ! !! SQRSL computes transformations, projections, and least squares solutions. ! ! ! Discussion: ! ! SQRSL requires the output of SQRDC. ! ! For K <= min(N,P), let AK be the matrix ! ! AK = ( A(JPVT(1)), A(JPVT(2)), ..., A(JPVT(K)) ) ! ! formed from columns JPVT(1), ..., JPVT(K) of the original ! N by P matrix A that was input to SQRDC. If no pivoting was ! done, AK consists of the first K columns of A in their ! original order. SQRDC produces a factored orthogonal matrix Q ! and an upper triangular matrix R such that ! ! AK = Q * (R) ! (0) ! ! This information is contained in coded form in the arrays ! A and QRAUX. ! ! The parameters QY, QTY, B, RSD, and AB are not referenced ! if their computation is not requested and in this case ! can be replaced by dummy variables in the calling program. ! To save storage, the user may in some cases use the same ! array for different parameters in the calling sequence. A ! frequently occuring example is when one wishes to compute ! any of B, RSD, or AB and does not need Y or QTY. In this ! case one may identify Y, QTY, and one of B, RSD, or AB, while ! providing separate arrays for anything else that is to be ! computed. ! ! Thus the calling sequence ! ! call sqrsl ( a, lda, n, k, qraux, y, dum, y, b, y, dum, 110, info ) ! ! will result in the computation of B and RSD, with RSD ! overwriting Y. More generally, each item in the following ! list contains groups of permissible identifications for ! a single calling sequence. ! ! 1. (Y,QTY,B) (RSD) (AB) (QY) ! ! 2. (Y,QTY,RSD) (B) (AB) (QY) ! ! 3. (Y,QTY,AB) (B) (RSD) (QY) ! ! 4. (Y,QY) (QTY,B) (RSD) (AB) ! ! 5. (Y,QY) (QTY,RSD) (B) (AB) ! ! 6. (Y,QY) (QTY,AB) (B) (RSD) ! ! In any group the value returned in the array allocated to ! the group corresponds to the last member of the group. ! ! Reference: ! ! Dongarra, Moler, Bunch and Stewart, ! LINPACK User's Guide, ! SIAM, (Society for Industrial and Applied Mathematics), ! 3600 University City Science Center, ! Philadelphia, PA, 19104-2688. ! ISBN 0-89871-172-X ! ! Parameters: ! ! Input, real A(LDA,P), contains the output of SQRDC. ! ! Input, integer LDA, the leading dimension of the array A. ! ! Input, integer N, the number of rows of the matrix AK. It must ! have the same value as N in SQRDC. ! ! Input, integer K, the number of columns of the matrix AK. K ! must not be greater than min(N,P), where P is the same as in the ! calling sequence to SQRDC. ! ! Input, real QRAUX(P), the auxiliary output from SQRDC. ! ! Input, real Y(N), a vector to be manipulated by SQRSL. ! ! Output, real QY(N), contains Q * Y, if requested. ! ! Output, real QTY(N), contains Q' * Y, if requested. ! ! Output, real B(K), the solution of the least squares problem ! minimize norm2 ( Y - AK * B), ! if its computation has been requested. Note that if pivoting was ! requested in SQRDC, the J-th component of B will be associated with ! column JPVT(J) of the original matrix A that was input into SQRDC. ! ! Output, real RSD(N), the least squares residual Y - AK * B, ! if its computation has been requested. RSD is also the orthogonal ! projection of Y onto the orthogonal complement of the column space ! of AK. ! ! Output, real AB(N), the least squares approximation Ak * B, if its ! computation has been requested. AB is also the orthogonal projection ! of Y onto the column space of A. ! ! Input, integer JOB, specifies what is to be computed. JOB has ! the decimal expansion ABCDE, with the following meaning: ! ! if A /= 0, compute QY. ! if B /= 0, compute QTY. ! if C /= 0, compute QTY and B. ! if D /= 0, compute QTY and RSD. ! if E /= 0, compute QTY and AB. ! ! Note that a request to compute B, RSD, or AB automatically triggers ! the computation of QTY, for which an array must be provided in the ! calling sequence. ! ! Output, integer INFO, is zero unless the computation of B has ! been requested and R is exactly singular. In this case, INFO is the ! index of the first zero diagonal element of R, and B is left unaltered. ! implicit none ! integer k integer lda integer n ! real a(lda,*) real ab(n) real b(k) logical cab logical cb logical cqty logical cqy logical cr logical cxb integer i integer info integer j integer jj integer job integer ju integer kp1 real qraux(*) real qty(n) real qy(n) real rsd(n) real sdot real t real temp real y(n) ! ! set info flag. ! info = 0 ! ! Determine what is to be computed. ! cqy = job / 10000 /= 0 cqty = mod ( job, 10000 ) /= 0 cb = mod ( job, 1000 ) / 100 /= 0 cr = mod ( job, 100 ) / 10 /= 0 cab = mod ( job, 10 ) /= 0 ju = min ( k, n-1 ) ! ! Special action when N = 1. ! if ( ju == 0 ) then if ( cqy ) then qy(1) = y(1) end if if ( cqty ) then qty(1) = y(1) end if if ( cab ) then ab(1) = y(1) end if if ( cb ) then if ( a(1,1) == 0.0E+00 ) then info = 1 else b(1) = y(1) / a(1,1) end if end if if ( cr ) then rsd(1) = 0.0E+00 end if return end if ! ! Set up to compute QY or QTY. ! if ( cqy ) then qy(1:n) = y(1:n) end if if ( cqty ) then qty(1:n) = y(1:n) end if ! ! Compute QY. ! if ( cqy ) then do jj = 1, ju j = ju - jj + 1 if ( qraux(j) /= 0.0E+00 ) then temp = a(j,j) a(j,j) = qraux(j) t = - sdot ( n-j+1, a(j,j), 1, qy(j), 1 ) / a(j,j) call saxpy ( n-j+1, t, a(j,j), 1, qy(j), 1 ) a(j,j) = temp end if end do end if ! ! Compute Q'*Y. ! if ( cqty ) then do j = 1, ju if ( qraux(j) /= 0.0E+00 ) then temp = a(j,j) a(j,j) = qraux(j) t = - sdot ( n-j+1, a(j,j), 1, qty(j), 1 ) / a(j,j) call saxpy ( n-j+1, t, a(j,j), 1, qty(j), 1 ) a(j,j) = temp end if end do end if ! ! Set up to compute B, RSD, or AB. ! if ( cb ) then b(1:k) = qty(1:k) end if kp1 = k + 1 if ( cab ) then ab(1:k) = qty(1:k) end if if ( cr .and. k < n ) then rsd(k+1:n) = qty(k+1:n) end if if ( cab .and. k+1 <= n ) then ab(k+1:n) = 0.0E+00 end if if ( cr ) then rsd(1:k) = 0.0E+00 end if ! ! Compute B. ! if ( cb ) then do jj = 1, k j = k - jj + 1 if ( a(j,j) == 0.0E+00 ) then info = j exit end if b(j) = b(j)/a(j,j) if ( j /= 1 ) then t = -b(j) call saxpy ( j-1, t, a(1,j), 1, b, 1 ) end if end do end if if ( cr .or. cab ) then ! ! Compute RSD or AB as required. ! do jj = 1, ju j = ju - jj + 1 if ( qraux(j) /= 0.0E+00 ) then temp = a(j,j) a(j,j) = qraux(j) if ( cr ) then t = - sdot ( n-j+1, a(j,j), 1, rsd(j), 1 ) / a(j,j) call saxpy ( n-j+1, t, a(j,j), 1, rsd(j), 1 ) end if if ( cab ) then t = - sdot ( n-j+1, a(j,j), 1, ab(j), 1 ) / a(j,j) call saxpy ( n-j+1, t, a(j,j), 1, ab(j), 1 ) end if a(j,j) = temp end if end do end if return end subroutine srot ( n, x, incx, y, incy, c, s ) ! !******************************************************************************* ! !! SROT applies a plane rotation. ! ! ! Reference: ! ! Lawson, Hanson, Kincaid, Krogh, ! Basic Linear Algebra Subprograms for Fortran Usage, ! Algorithm 539, ! ACM Transactions on Mathematical Software, ! Volume 5, Number 3, September 1979, pages 308-323. ! ! Modified: ! ! 08 April 1999 ! ! Parameters: ! ! Input, integer N, the number of entries in the vectors. ! ! Input/output, real X(*), one of the vectors to be rotated. ! ! Input, integer INCX, the increment between successive entries of X. ! ! Input/output, real Y(*), one of the vectors to be rotated. ! ! Input, integer INCY, the increment between successive elements of Y. ! ! Input, real C, S, parameters (presumably the cosine and sine of ! some angle) that define a plane rotation. ! implicit none ! real c integer i integer incx integer incy integer ix integer iy integer n real s real stemp real x(*) real y(*) ! if ( n <= 0 ) then else if ( incx == 1 .and. incy == 1 ) then do i = 1, n stemp = c * x(i) + s * y(i) y(i) = c * y(i) - s * x(i) x(i) = stemp end do else if ( incx >= 0 ) then ix = 1 else ix = ( - n + 1 ) * incx + 1 end if if ( incy >= 0 ) then iy = 1 else iy = ( - n + 1 ) * incy + 1 end if do i = 1, n stemp = c * x(ix) + s * y(iy) y(iy) = c * y(iy) - s * x(ix) x(ix) = stemp ix = ix + incx iy = iy + incy end do end if return end subroutine srotg ( sa, sb, c, s ) ! !******************************************************************************* ! !! SROTG constructs a Givens plane rotation. ! ! ! Reference: ! ! Lawson, Hanson, Kincaid, Krogh, ! Basic Linear Algebra Subprograms for Fortran Usage, ! Algorithm 539, ! ACM Transactions on Mathematical Software, ! Volume 5, Number 3, September 1979, pages 308-323. ! ! Modified: ! ! 08 April 1999 ! ! Parameters: ! ! Input/output, real SA, SB, ... ! ! Output, real C, S, ... ! implicit none ! real c real r real roe real s real sa real sb real scale real z ! if ( abs ( sa ) > abs ( sb ) ) then roe = sa else roe = sb end if scale = abs ( sa ) + abs ( sb ) if ( scale == 0.0E+00 ) then c = 1.0E+00 s = 0.0E+00 r = 0.0E+00 else r = scale * sqrt ( ( sa / scale )**2 + ( sb / scale )**2 ) r = sign ( 1.0E+00, roe ) * r c = sa / r s = sb / r end if if ( abs ( c ) > 0.0E+00 .and. abs ( c ) <= s ) then z = 1.0E+00 / c else z = s end if sa = r sb = z return end subroutine sscal ( n, sa, x, incx ) ! !******************************************************************************* ! !! SSCAL scales a vector by a constant. ! ! ! Reference: ! ! Lawson, Hanson, Kincaid, Krogh, ! Basic Linear Algebra Subprograms for Fortran Usage, ! Algorithm 539, ! ACM Transactions on Mathematical Software, ! Volume 5, Number 3, September 1979, pages 308-323. ! ! Modified: ! ! 08 April 1999 ! ! Parameters: ! ! Input, integer N, the number of entries in the vector. ! ! Input, real SA, the multiplier. ! ! Input/output, real X(*), the vector to be scaled. ! ! Input, integer INCX, the increment between successive entries of X. ! implicit none ! integer i integer incx integer ix integer m integer n real sa real x(*) ! if ( n <= 0 ) then else if ( incx == 1 ) then m = mod ( n, 5 ) x(1:m) = sa * x(1:m) do i = m+1, n, 5 x(i) = sa * x(i) x(i+1) = sa * x(i+1) x(i+2) = sa * x(i+2) x(i+3) = sa * x(i+3) x(i+4) = sa * x(i+4) end do else if ( incx >= 0 ) then ix = 1 else ix = ( - n + 1 ) * incx + 1 end if do i = 1, n x(ix) = sa * x(ix) ix = ix + incx end do end if return end subroutine ssico ( a, lda, n, kpvt, rcond, z ) ! !******************************************************************************* ! !! SSICO factors a real symmetric matrix and estimates its condition. ! ! ! Discussion: ! ! If RCOND is not needed, SSIFA is slightly faster. ! ! To solve A*X = B, follow SSICO by SSISL. ! ! To compute inverse(A)*C, follow SSICO by SSISL. ! ! To compute inverse(A), follow SSICO by SSIDI. ! ! To compute determinant(A), follow SSICO by SSIDI. ! ! To compute inertia(A), follow SSICO by SSIDI. ! ! Reference: ! ! Dongarra, Moler, Bunch and Stewart, ! LINPACK User's Guide, ! SIAM, (Society for Industrial and Applied Mathematics), ! 3600 University City Science Center, ! Philadelphia, PA, 19104-2688. ! ISBN 0-89871-172-X ! ! Parameters: ! ! Input/output, real A(LDA,N). On input, the symmetric matrix to be ! factored. Only the diagonal and upper triangle are used. ! On output, a block diagonal matrix and the multipliers which ! were used to obtain it. The factorization can be written A = U*D*U' ! where U is a product of permutation and unit upper triangular ! matrices, U' is the transpose of U, and D is block diagonal ! with 1 by 1 and 2 by 2 blocks. ! ! Input, integer LDA, the leading dimension of the array A. ! ! Input, integer N, the order of the matrix. ! ! Output, integer KPVT(N), pivot indices. ! ! Output, real RCOND, an estimate of the reciprocal condition of A. ! For the system A*X = B, relative perturbations in A and B of size ! EPSILON may cause relative perturbations in X of size EPSILON/RCOND. ! If RCOND is so small that the logical expression ! 1.0E+00 + RCOND == 1.0E+00 ! is true, then A may be singular to working precision. In particular, ! RCOND is zero if exact singularity is detected or the estimate underflows. ! ! Output, real Z(N), a work vector whose contents are usually unimportant. ! If A is close to a singular matrix, then Z is an approximate null vector ! in the sense that ! norm(A*Z) = RCOND * norm(A) * norm(Z). ! implicit none ! integer lda integer n ! real a(lda,n) real ak real akm1 real anorm real bk real bkm1 real denom real ek integer i integer info integer j integer k integer kp integer kps integer kpvt(n) integer ks real rcond real s real sdot real t real ynorm real z(n) ! ! Find the norm of A, using only entries in the upper half of the matrix. ! do j = 1, n z(j) = sum ( abs ( a(1:j-1,j) ) ) do i = 1, j-1 z(i) = z(i) + abs ( a(i,j) ) end do end do anorm = maxval ( z(1:n) ) ! ! Factor. ! call ssifa ( a, lda, n, kpvt, info ) ! ! RCOND = 1/(norm(A)*(estimate of norm(inverse(A)))). ! ! Estimate = norm(Z)/norm(Y) where A*Z = Y and A*Y = E. ! ! The components of E are chosen to cause maximum local ! growth in the elements of W where U*D*W = E. ! ! The vectors are frequently rescaled to avoid overflow. ! ! Solve U*D*W = E. ! ek = 1.0E+00 z(1:n) = 0.0E+00 k = n do while ( k /= 0 ) if ( kpvt(k) < 0 ) then ks = 2 else ks = 1 end if kp = abs ( kpvt(k) ) kps = k + 1 - ks if ( kp /= kps ) then call r_swap ( z(kps), z(kp) ) end if if ( z(k) /= 0.0E+00 ) then ek = sign ( ek, z(k) ) end if z(k) = z(k) + ek call saxpy ( k-ks, z(k), a(1,k), 1, z(1), 1 ) if ( ks /= 1 ) then if ( z(k-1) /= 0.0E+00 ) then ek = sign ( ek, z(k-1) ) end if z(k-1) = z(k-1) + ek call saxpy ( k-ks, z(k-1), a(1,k-1), 1, z(1), 1 ) end if if ( ks /= 2 ) then if ( abs ( z(k) ) > abs ( a(k,k) ) ) then s = abs ( a(k,k) ) / abs ( z(k) ) z(1:n) = s * z(1:n) ek = s * ek end if if ( a(k,k) /= 0.0E+00 ) then z(k) = z(k) / a(k,k) else z(k) = 1.0E+00 end if else ak = a(k,k) / a(k-1,k) akm1 = a(k-1,k-1) / a(k-1,k) bk = z(k) / a(k-1,k) bkm1 = z(k-1) / a(k-1,k) denom = ak * akm1 - 1.0E+00 z(k) = ( akm1 * bk - bkm1 ) / denom z(k-1) = ( ak * bkm1 - bk ) / denom end if k = k - ks end do z(1:n) = z(1:n) / sum ( abs ( z(1:n) ) ) ! ! Solve U'*Y = W. ! k = 1 do while ( k <= n ) ks = 1 if ( kpvt(k) < 0 ) ks = 2 if ( k /= 1 ) then z(k) = z(k) + sdot ( k-1, a(1,k), 1, z(1), 1 ) if ( ks == 2 ) then z(k+1) = z(k+1) + sdot ( k-1, a(1,k+1), 1, z(1), 1 ) end if kp = abs ( kpvt(k) ) if ( kp /= k ) then call r_swap ( z(k), z(kp) ) end if end if k = k + ks end do z(1:n) = z(1:n) / sum ( abs ( z(1:n) ) ) ynorm = 1.0E+00 ! ! Solve U*D*V = Y. ! k = n do while ( k /= 0 ) if ( kpvt(k) < 0 ) then ks = 2 else ks = 1 end if if ( k /= ks ) then kp = abs ( kpvt(k) ) kps = k + 1 - ks if ( kp /= kps ) then call r_swap ( z(kps), z(kp) ) end if call saxpy ( k-ks, z(k), a(1,k), 1, z(1), 1 ) if ( ks == 2 ) then call saxpy ( k-ks, z(k-1), a(1,k-1), 1, z(1), 1 ) end if end if if ( ks /= 2 ) then if ( abs ( z(k) ) > abs ( a(k,k) ) ) then s = abs ( a(k,k) ) / abs ( z(k) ) z(1:n) = s * z(1:n) ynorm = s * ynorm end if if ( a(k,k) /= 0.0E+00 ) then z(k) = z(k) / a(k,k) else z(k) = 1.0E+00 end if else ak = a(k,k) / a(k-1,k) akm1 = a(k-1,k-1) / a(k-1,k) bk = z(k) / a(k-1,k) bkm1 = z(k-1) / a(k-1,k) denom = ak * akm1 - 1.0E+00 z(k) = ( akm1 * bk - bkm1 ) / denom z(k-1) = ( ak * bkm1 - bk ) / denom end if k = k - ks end do s = 1.0E+00 / sum ( abs ( z(1:n) ) ) z(1:n) = s * z(1:n) ynorm = s * ynorm ! ! Solve U'*Z = V. ! k = 1 do while ( k <= n ) if ( kpvt(k) < 0 ) then ks = 2 else ks = 1 end if if ( k /= 1 ) then z(k) = z(k) + sdot ( k-1, a(1,k), 1, z(1), 1 ) if ( ks == 2 ) then z(k+1) = z(k+1) + sdot ( k-1, a(1,k+1), 1, z(1), 1 ) end if kp = abs ( kpvt(k) ) if ( kp /= k ) then call r_swap ( z(k), z(kp) ) end if end if k = k + ks end do ! ! Make ZNORM = 1.0. ! s = 1.0E+00 / sum ( abs ( z(1:n) ) ) z(1:n) = s * z(1:n) ynorm = s * ynorm if ( anorm /= 0.0E+00 ) then rcond = ynorm / anorm else rcond = 0.0E+00 end if return end subroutine ssidi ( a, lda, n, kpvt, det, inert, work, job ) ! !******************************************************************************* ! !! SSIDI computes the determinant, inertia and inverse of a real symmetric matrix. ! ! ! Discussion: ! ! SSIDI uses the factors from SSIFA. ! ! A division by zero may occur if the inverse is requested ! and SSICO has set RCOND == 0.0E+00 or SSIFA has set INFO /= 0. ! ! Variables not requested by JOB are not used. ! ! Reference: ! ! Dongarra, Moler, Bunch and Stewart, ! LINPACK User's Guide, ! SIAM, (Society for Industrial and Applied Mathematics), ! 3600 University City Science Center, ! Philadelphia, PA, 19104-2688. ! ISBN 0-89871-172-X ! ! Parameters: ! ! Input/output, real A(LDA,N). On input, the output from SSIFA. ! On output, the upper triangle of the inverse of the original matrix. ! The strict lower triangle is never referenced. ! ! Input, integer LDA, the leading dimension of the array A. ! ! Input, integer N, the order of the matrix. ! ! Input, integer KPVT(N), the pivot vector from SSIFA. ! ! Output, real DET(2), the determinant of the original matrix. ! determinant = DET(1) * 10.0**DET(2) ! with 1.0E+00 <= abs ( DET(1) ) < 10.0E+00 or DET(1) = 0.0. ! ! Output, integer INERT(3), the inertia of the original matrix. ! INERT(1) = number of positive eigenvalues. ! INERT(2) = number of negative eigenvalues. ! INERT(3) = number of zero eigenvalues. ! ! Workspace, real WORK(N). ! ! Input, integer JOB, specifies the tasks. ! JOB has the decimal expansion ABC where ! If C /= 0, the inverse is computed, ! If B /= 0, the determinant is computed, ! If A /= 0, the inertia is computed. ! For example, JOB = 111 gives all three. ! implicit none ! integer lda integer n ! real a(lda,n) real ak real akkp1 real akp1 real d real det(2) logical dodet logical doert logical doinv integer inert(3) integer j integer jb integer job integer k integer kpvt(n) integer ks integer kstep real sdot real t real temp real, parameter :: ten = 10.0E+00 real work(n) ! doinv = mod ( job, 10 ) /= 0 dodet = mod ( job, 100 ) / 10 /= 0 doert = mod ( job, 1000 ) / 100 /= 0 if ( dodet .or. doert ) then if ( doert ) then inert(1) = 0 inert(2) = 0 inert(3) = 0 end if if ( dodet ) then det(1) = 1.0E+00 det(2) = 0.0E+00 end if t = 0.0E+00 do k = 1, n d = a(k,k) ! ! 2 by 2 block. ! ! use det (d s) = (d/t * c - t) * t, t = abs ( s ) ! (s c) ! to avoid underflow/overflow troubles. ! ! Take two passes through scaling. Use T for flag. ! if ( kpvt(k) <= 0 ) then if ( t == 0.0E+00 ) then t = abs ( a(k,k+1) ) d = ( d / t ) * a(k+1,k+1) - t else d = t t = 0.0E+00 end if end if if ( doert ) then if ( d > 0.0E+00 ) then inert(1) = inert(1) + 1 else if ( d < 0.0E+00 ) then inert(2) = inert(2) + 1 else if ( d == 0.0E+00 ) then inert(3) = inert(3) + 1 end if end if if ( dodet ) then det(1) = d * det(1) if ( det(1) /= 0.0E+00 ) then do while ( abs ( det(1) ) < 1.0E+00 ) det(1) = ten * det(1) det(2) = det(2) - 1.0E+00 end do do while ( abs ( det(1) ) >= ten ) det(1) = det(1) / ten det(2) = det(2) + 1.0E+00 end do end if end if end do end if ! ! Compute inverse(A). ! if ( doinv ) then k = 1 do while ( k <= n ) if ( kpvt(k) >= 0 ) then ! ! 1 by 1. ! a(k,k) = 1.0E+00 / a(k,k) if ( k >= 2 ) then work(1:k-1) = a(1:k-1,k) do j = 1, k-1 a(j,k) = sdot ( j, a(1,j), 1, work, 1 ) call saxpy ( j-1, work(j), a(1,j), 1, a(1,k), 1 ) end do a(k,k) = a(k,k) + sdot ( k-1, work, 1, a(1,k), 1 ) end if kstep = 1 ! ! 2 by 2. ! else t = abs ( a(k,k+1) ) ak = a(k,k) / t akp1 = a(k+1,k+1) / t akkp1 = a(k,k+1) / t d = t * ( ak * akp1 - 1.0E+00 ) a(k,k) = akp1 / d a(k+1,k+1) = ak / d a(k,k+1) = - akkp1 / d if ( k >= 2 ) then call scopy ( k-1, a(1,k+1), 1, work, 1 ) do j = 1, k-1 a(j,k+1) = sdot ( j, a(1,j), 1, work, 1 ) call saxpy ( j-1, work(j), a(1,j), 1, a(1,k+1), 1 ) end do a(k+1,k+1) = a(k+1,k+1) + sdot ( k-1, work, 1, a(1,k+1), 1 ) a(k,k+1) = a(k,k+1) + sdot ( k-1, a(1,k), 1, a(1,k+1), 1 ) call scopy ( k-1, a(1,k), 1, work, 1 ) do j = 1, k-1 a(j,k) = sdot ( j, a(1,j), 1, work, 1 ) call saxpy ( j-1, work(j), a(1,j), 1, a(1,k), 1 ) end do a(k,k) = a(k,k) + sdot ( k-1, work, 1, a(1,k), 1 ) end if kstep = 2 end if ! ! Swap. ! ks = abs ( kpvt(k) ) if ( ks /= k ) then call sswap ( ks, a(1,ks), 1, a(1,k), 1 ) do jb = ks, k j = k + ks - jb call r_swap ( a(j,k), a(ks,j) ) end do if ( kstep /= 1 ) then call r_swap ( a(ks,k+1), a(k,k+1) ) end if end if k = k + kstep end do end if return end subroutine ssifa ( a, lda, n, kpvt, info ) ! !******************************************************************************* ! !! SSIFA factors a real symmetric matrix. ! ! ! Discussion: ! ! To solve A*X = B, follow SSIFA by SSISL. ! ! To compute inverse(A)*C, follow SSIFA by SSISL. ! ! To compute determinant(A), follow SSIFA by SSIDI. ! ! To compute inertia(A), follow SSIFA by SSIDI. ! ! To compute inverse(A), follow SSIFA by SSIDI. ! ! Reference: ! ! Dongarra, Moler, Bunch and Stewart, ! LINPACK User's Guide, ! SIAM, (Society for Industrial and Applied Mathematics), ! 3600 University City Science Center, ! Philadelphia, PA, 19104-2688. ! ISBN 0-89871-172-X ! ! Parameters: ! ! Input/output, real A(LDA,N). On input, the symmetric matrix to be ! factored. Only the diagonal and upper triangle are used. ! On output, a block diagonal matrix and the multipliers which ! were used to obtain it. The factorization can be written A = U*D*U' ! where U is a product of permutation and unit upper triangular ! matrices, U' is the transpose of U, and D is block diagonal ! with 1 by 1 and 2 by 2 blocks. ! ! Input, integer LDA, the leading dimension of the array A. ! ! Input, integer N, the order of the matrix. ! ! Output, integer KPVT(N), the pivot indices. ! ! Output, integer INFO, error flag. ! 0, normal value. ! K, if the K-th pivot block is singular. This is not an error ! condition for this subroutine, but it does indicate that SSISL ! or SSIDI may divide by zero if called. ! implicit none ! integer lda integer n ! real a(lda,n) real absakk real ak real akm1 real alpha real bk real bkm1 real colmax real denom integer imax integer imaxp1 integer info integer isamax integer j integer jj integer jmax integer k integer km1 integer kpvt(n) integer kstep real mulk real mulkm1 real rowmax logical swap real t ! ! ALPHA is used in choosing pivot block size. ! alpha = ( 1.0E+00 + sqrt ( 17.0E+00 ) ) / 8.0E+00 info = 0 ! ! Main loop on K, which goes from N to 1. ! k = n do while ( k > 0 ) if ( k == 1 ) then kpvt(1) = 1 if ( a(1,1) == 0.0E+00 ) then info = 1 end if return end if ! ! This section of code determines the kind of ! elimination to be performed. When it is completed, ! KSTEP will be set to the size of the pivot block, and ! SWAP will be set to .true. if an interchange is required. ! km1 = k - 1 absakk = abs ( a(k,k) ) ! ! Determine the largest off-diagonal element in column K. ! imax = isamax ( k-1, a(1,k), 1 ) colmax = abs ( a(imax,k) ) if ( absakk >= alpha * colmax ) then kstep = 1 swap = .false. else ! ! Determine the largest off-diagonal element in row IMAX. ! rowmax = 0.0E+00 imaxp1 = imax + 1 do j = imaxp1, k rowmax = max ( rowmax, abs ( a(imax,j) ) ) end do if ( imax /= 1 ) then jmax = isamax ( imax-1, a(1,imax), 1 ) rowmax = max ( rowmax, abs ( a(jmax,imax) ) ) end if if ( abs ( a(imax,imax) ) >= alpha * rowmax ) then kstep = 1 swap = .true. else if ( absakk >= alpha * colmax * ( colmax / rowmax ) ) then kstep = 1 swap = .false. else kstep = 2 swap = ( imax /= k-1 ) end if end if ! ! Column K is zero. ! Set INFO and iterate the loop. ! if ( max ( absakk, colmax ) == 0.0E+00 ) then kpvt(k) = k info = k ! ! 1 x 1 pivot block. ! ! Perform an interchange. ! else if ( kstep /= 2 ) then if ( swap ) then call sswap ( imax, a(1,imax), 1, a(1,k), 1 ) do jj = imax, k j = k + imax - jj call r_swap ( a(j,k), a(imax,j) ) end do end if ! ! Perform the elimination. ! do jj = 1, k-1 j = k - jj mulk = - a(j,k) / a(k,k) t = mulk call saxpy ( j, t, a(1,k), 1, a(1,j), 1 ) a(j,k) = mulk end do ! ! Set the pivot array. ! if ( swap ) then kpvt(k) = imax else kpvt(k) = k end if ! ! 2 x 2 pivot block. ! ! Perform an interchange. ! else if ( swap ) then call sswap ( imax, a(1,imax), 1, a(1,k-1), 1 ) do jj = imax, k-1 j = k-1 + imax - jj call r_swap ( a(j,k-1), a(imax,j) ) end do t = a(k-1,k) a(k-1,k) = a(imax,k) a(imax,k) = t end if ! ! Perform the elimination. ! if ( k-2 /= 0 ) then ak = a(k,k) / a(k-1,k) akm1 = a(k-1,k-1) / a(k-1,k) denom = 1.0E+00 - ak * akm1 do jj = 1, k-2 j = k-1 - jj bk = a(j,k) / a(k-1,k) bkm1 = a(j,k-1) / a(k-1,k) mulk = ( akm1 * bk - bkm1 ) / denom mulkm1 = ( ak * bkm1 - bk ) / denom t = mulk call saxpy ( j, t, a(1,k), 1, a(1,j), 1 ) t = mulkm1 call saxpy ( j, t, a(1,k-1), 1, a(1,j), 1 ) a(j,k) = mulk a(j,k-1) = mulkm1 end do end if ! ! Set the pivot array. ! if ( swap ) then kpvt(k) = -imax else kpvt(k) = 1 - k end if kpvt(k-1) = kpvt(k) end if k = k - kstep end do return end subroutine ssisl ( a, lda, n, kpvt, b ) ! !******************************************************************************* ! !! SSISL solves a real symmetric system factored by SSIFA. ! ! ! Discussion: ! ! To compute inverse(A) * C where C is a matrix with P columns ! ! call ssifa ( a, lda, n, kpvt, info ) ! ! if ( info == 0 ) then ! do j = 1, p ! call ssisl ( a, lda, n, kpvt, c(1,j) ) ! end do ! end if ! ! A division by zero may occur if the inverse is requested ! and SSICO has set RCOND == 0.0E+00 or SSIFA has set INFO /= 0. ! ! Reference: ! ! Dongarra, Moler, Bunch and Stewart, ! LINPACK User's Guide, ! SIAM, (Society for Industrial and Applied Mathematics), ! 3600 University City Science Center, ! Philadelphia, PA, 19104-2688. ! ISBN 0-89871-172-X ! ! Parameters: ! ! Input, real A(LDA,N), the output from SSIFA. ! ! Input, integer LDA, the leading dimension of the array A. ! ! Input, integer N, the order of the matrix. ! ! Input, integer KPVT(N), the pivot vector from SSIFA. ! ! Input/output, real B(N). On input, the right hand side. ! On output, the solution. ! implicit none ! integer lda integer n ! real a(lda,n) real ak real akm1 real b(n) real bk real bkm1 real denom integer k integer kp integer kpvt(n) real sdot real temp ! ! Loop backward applying the transformations and D inverse to B. ! k = n do while ( k > 0 ) if ( kpvt(k) >= 0 ) then ! ! 1 x 1 pivot block. ! if ( k /= 1 ) then kp = kpvt(k) ! ! Interchange. ! if ( kp /= k ) then call r_swap ( b(k), b(kp) ) end if ! ! Apply the transformation. ! call saxpy ( k-1, b(k), a(1,k), 1, b(1), 1 ) end if ! ! Apply D inverse. ! b(k) = b(k)/a(k,k) k = k - 1 else ! ! 2 x 2 pivot block. ! if ( k /= 2 ) then kp = abs (kpvt(k)) ! ! Interchange. ! if ( kp /= k-1 ) then call r_swap ( b(k-1), b(kp) ) end if ! ! Apply the transformation. ! call saxpy ( k-2, b(k), a(1,k), 1, b(1), 1 ) call saxpy ( k-2, b(k-1), a(1,k-1), 1, b(1), 1 ) end if ! ! Apply D inverse. ! ak = a(k,k) / a(k-1,k) akm1 = a(k-1,k-1) / a(k-1,k) bk = b(k) / a(k-1,k) bkm1 = b(k-1) / a(k-1,k) denom = ak * akm1 - 1.0E+00 b(k) = ( akm1 * bk - bkm1 ) / denom b(k-1) = ( ak * bkm1 - bk ) / denom k = k - 2 end if end do ! ! Loop forward applying the transformations. ! k = 1 do while ( k <= n ) if ( kpvt(k) >= 0 ) then ! ! 1 x 1 pivot block. ! if ( k /= 1 ) then ! ! Apply the transformation. ! b(k) = b(k) + sdot ( k-1, a(1,k), 1, b(1), 1 ) kp = kpvt(k) ! ! Interchange. ! if ( kp /= k ) then call r_swap ( b(k), b(kp) ) end if end if k = k + 1 else ! ! 2 x 2 pivot block. ! if ( k /= 1 ) then ! ! Apply the transformation. ! b(k) = b(k) + sdot ( k-1, a(1,k), 1, b(1), 1 ) b(k+1) = b(k+1) + sdot ( k-1, a(1,k+1), 1, b(1), 1 ) kp = abs (kpvt(k)) ! ! Interchange. ! if ( kp /= k ) then call r_swap ( b(k), b(kp) ) end if end if k = k + 2 end if end do return end subroutine sspco ( ap, n, kpvt, rcond, z ) ! !******************************************************************************* ! !! SSPCO factors a real symmetric matrix stored in packed form. ! ! ! Discussion: ! ! SSPCO uses elimination with symmetric pivoting and estimates ! the condition of the matrix. ! ! If RCOND is not needed, SSPFA is slightly faster. ! ! To solve A*X = B, follow SSPCO by SSPSL. ! ! To compute inverse(A)*C, follow SSPCO by SSPSL. ! ! To compute inverse(A), follow SSPCO by SSPDI. ! ! To compute determinant(A), follow SSPCO by SSPDI. ! ! To compute inertia(A), follow SSPCO by SSPDI. ! ! Packed storage: ! ! The following program segment will pack the upper triangle of a ! symmetric matrix. ! ! k = 0 ! do j = 1, n ! do i = 1, j ! k = k + 1 ! ap(k) = a(i,j) ! end do ! end do ! ! Reference: ! ! Dongarra, Moler, Bunch and Stewart, ! LINPACK User's Guide, ! SIAM, (Society for Industrial and Applied Mathematics), ! 3600 University City Science Center, ! Philadelphia, PA, 19104-2688. ! ISBN 0-89871-172-X ! ! Parameters: ! ! Input/output, real AP(N*(N+1)/2). On input, the packed form of a ! symmetric matrix A. The columns of the upper triangle are stored ! sequentially in a one-dimensional array. On output, a block diagonal ! matrix and the multipliers which were used to obtain it, stored in ! packed form. The factorization can be written A = U*D*U' ! where U is a product of permutation and unit upper triangular ! matrices, U' is the transpose of U, and D is block diagonal ! with 1 by 1 and 2 by 2 blocks. ! ! Input, integer N, the order of the matrix. ! ! Output, integer KPVT(N), the pivot indices. ! ! Output, real RCOND, an estimate of the reciprocal condition of A. ! For the system A*X = B, relative perturbations in A and B of size ! EPSILON may cause relative perturbations in X of size EPSILON/RCOND. ! If RCOND is so small that the logical expression ! 1.0E+00 + RCOND == 1.0E+00 ! is true, then A may be singular to working precision. In particular, ! RCOND is zero if exact singularity is detected or the estimate underflows. ! ! Output, real Z(N) a work vector whose contents are usually unimportant. ! If A is close to a singular matrix, then Z is an approximate null ! vector in the sense that ! norm(A*Z) = RCOND * norm(A) * norm(Z). ! implicit none ! integer n ! real ak real akm1 real anorm real ap((n*(n+1))/2) real bk real bkm1 real denom real ek integer i integer ij integer ik integer ikm1 integer ikp1 integer info integer j integer j1 integer k integer kk integer km1k integer km1km1 integer kp integer kps integer kpvt(n) integer ks real rcond real s real sdot real t real ynorm real z(n) ! ! Find norm of A using only upper half. ! j1 = 1 do j = 1, n z(j) = sum ( abs ( ap(j1:j1+j-1) ) ) ij = j1 j1 = j1 + j do i = 1, j-1 z(i) = z(i) + abs ( ap(ij) ) ij = ij + 1 end do end do anorm = maxval ( z(1:n) ) ! ! Factor. ! call sspfa ( ap, n, kpvt, info ) ! ! RCOND = 1/(norm(A)*(estimate of norm(inverse(A)))). ! ! Estimate = norm(Z)/norm(Y) where A*Z = Y and A*Y = E. ! ! The components of E are chosen to cause maximum local ! growth in the elements of W where U*D*W = E. ! ! The vectors are frequently rescaled to avoid overflow. ! ! Solve U*D*W = E. ! ek = 1.0E+00 z(1:n) = 0.0E+00 k = n ik = (n*(n - 1))/2 do while ( k /= 0 ) kk = ik + k ikm1 = ik - (k - 1) ks = 1 if ( kpvt(k) < 0 ) ks = 2 kp = abs (kpvt(k)) kps = k + 1 - ks if ( kp /= kps ) then t = z(kps) z(kps) = z(kp) z(kp) = t end if if ( z(k) /= 0.0E+00 ) then ek = sign ( ek, z(k) ) end if z(k) = z(k) + ek call saxpy ( k-ks, z(k), ap(ik+1), 1, z(1), 1 ) if ( ks /= 1 ) then if ( z(k-1) /= 0.0E+00 ) then ek = sign ( ek, z(k-1) ) end if z(k-1) = z(k-1) + ek call saxpy ( k-ks, z(k-1), ap(ikm1+1), 1, z(1), 1 ) end if if ( ks /= 2 ) then if ( abs ( z(k) ) > abs ( ap(kk) ) ) then s = abs ( ap(kk) ) / abs ( z(k) ) z(1:n) = s * z(1:n) ek = s * ek end if if ( ap(kk) /= 0.0E+00 ) then z(k) = z(k) / ap(kk) else z(k) = 1.0E+00 end if else km1k = ik + k - 1 km1km1 = ikm1 + k - 1 ak = ap(kk) / ap(km1k) akm1 = ap(km1km1) / ap(km1k) bk = z(k) / ap(km1k) bkm1 = z(k-1) / ap(km1k) denom = ak * akm1 - 1.0E+00 z(k) = ( akm1 * bk - bkm1 ) / denom z(k-1) = ( ak * bkm1 - bk ) / denom end if k = k - ks ik = ik - k if ( ks == 2 ) then ik = ik - (k + 1) end if end do z(1:n) = z(1:n) / sum ( abs ( z(1:n) ) ) ! ! Solve U'*Y = W. ! k = 1 ik = 0 do while ( k <= n ) if ( kpvt(k) < 0 ) then ks = 2 else ks = 1 end if if ( k /= 1 ) then z(k) = z(k) + sdot ( k-1, ap(ik+1), 1, z(1), 1 ) ikp1 = ik + k if ( ks == 2 ) then z(k+1) = z(k+1) + sdot ( k-1, ap(ikp1+1), 1, z(1), 1 ) end if kp = abs (kpvt(k)) if ( kp /= k ) then t = z(k) z(k) = z(kp) z(kp) = t end if end if ik = ik + k if ( ks == 2 ) ik = ik + (k + 1) k = k + ks end do s = 1.0E+00 / sum ( abs ( z(1:n) ) ) z(1:n) = s * z(1:n) ynorm = 1.0E+00 ! ! Solve U*D*V = Y. ! k = n ik = n * ( n - 1 ) / 2 do while ( k > 0 ) kk = ik + k ikm1 = ik - (k - 1) if ( kpvt(k) < 0 ) then ks = 2 else ks = 1 end if if ( k /= ks ) then kp = abs ( kpvt(k) ) kps = k + 1 - ks if ( kp /= kps ) then t = z(kps) z(kps) = z(kp) z(kp) = t end if call saxpy ( k-ks, z(k), ap(ik+1), 1, z(1), 1 ) if ( ks == 2 ) then call saxpy ( k-ks, z(k-1), ap(ikm1+1), 1, z(1), 1 ) end if end if if ( ks /= 2 ) then if ( abs ( z(k) ) > abs ( ap(kk) ) ) then s = abs ( ap(kk) ) / abs ( z(k) ) z(1:n) = s * z(1:n) ynorm = s * ynorm end if if ( ap(kk) /= 0.0E+00 ) then z(k) = z(k) / ap(kk) else z(k) = 1.0E+00 end if else km1k = ik + k - 1 km1km1 = ikm1 + k - 1 ak = ap(kk) / ap(km1k) akm1 = ap(km1km1) / ap(km1k) bk = z(k) / ap(km1k) bkm1 = z(k-1) / ap(km1k) denom = ak * akm1 - 1.0E+00 z(k) = ( akm1 * bk - bkm1 ) / denom z(k-1) = ( ak * bkm1 - bk ) / denom end if k = k - ks ik = ik - k if ( ks == 2 ) ik = ik - (k + 1) end do s = 1.0E+00 / sum ( abs ( z(1:n) ) ) z(1:n) = s * z(1:n) ynorm = s * ynorm ! ! Solve U'*Z = V. ! k = 1 ik = 0 do while ( k <= n ) ks = 1 if ( kpvt(k) < 0 ) ks = 2 if ( k /= 1 ) then z(k) = z(k) + sdot ( k-1, ap(ik+1), 1, z(1), 1 ) ikp1 = ik + k if ( ks == 2 ) then z(k+1) = z(k+1) + sdot ( k-1, ap(ikp1+1), 1, z(1), 1 ) end if kp = abs ( kpvt(k) ) if ( kp /= k ) then t = z(k) z(k) = z(kp) z(kp) = t end if end if ik = ik + k if ( ks == 2 ) ik = ik + (k + 1) k = k + ks end do ! ! Make ZNORM = 1.0. ! s = 1.0E+00 / sum ( abs ( z(1:n) ) ) z(1:n) = s * z(1:n) ynorm = s * ynorm if ( anorm /= 0.0E+00 ) then rcond = ynorm / anorm else rcond = 0.0E+00 end if return end subroutine sspdi ( ap, n, kpvt, det, inert, work, job ) ! !******************************************************************************* ! !! SSPDI computes the determinant, inertia and inverse of a real symmetric matrix. ! ! ! Discussion: ! ! SSPDI uses the factors from SSPFA, where the matrix is stored in ! packed form. ! ! A division by zero will occur if the inverse is requested ! and SSPCO has set RCOND == 0.0E+00 or SSPFA has set INFO /= 0. ! ! Variables not requested by JOB are not used. ! ! Reference: ! ! Dongarra, Moler, Bunch and Stewart, ! LINPACK User's Guide, ! SIAM, (Society for Industrial and Applied Mathematics), ! 3600 University City Science Center, ! Philadelphia, PA, 19104-2688. ! ISBN 0-89871-172-X ! ! Parameters: ! ! Input/output, real AP(N*(N+1)/2). On input, the output from SSPFA. ! On output, the upper triangle of the inverse of the original matrix, ! stored in packed form. The columns of the upper triangle are stored ! sequentially in a one-dimensional array. ! ! Input, integer N, the order of the matrix. ! ! Input, integer KPVT(N), the pivot vector from SSPFA. ! ! Output, real DET(2), the determinant of the original matrix. ! determinant = DET(1) * 10.0**DET(2) ! with 1.0E+00 <= abs ( DET(1) ) < 10.0E+00 or DET(1) = 0.0. ! ! Output, integer INERT(3), the inertia of the original matrix. ! INERT(1) = number of positive eigenvalues. ! INERT(2) = number of negative eigenvalues. ! INERT(3) = number of zero eigenvalues. ! ! Workspace, real WORK(N). ! ! Input, integer JOB, has the decimal expansion ABC where: ! if A /= 0, the inertia is computed, ! if B /= 0, the determinant is computed, ! if C /= 0, the inverse is computed. ! For example, JOB = 111 gives all three. ! implicit none ! integer n ! real ak real akkp1 real akp1 real ap((n*(n+1))/2) real d real det(2) logical dodet logical doert logical doinv integer ij integer ik integer ikp1 integer iks integer inert(3) integer j integer jb integer jk integer jkp1 integer job integer k integer kk integer kkp1 integer km1 integer kpvt(n) integer ks integer ksj integer kskp1 integer kstep real sdot real t real temp real, parameter :: ten = 10.0E+00 real work(n) ! doinv = mod ( job, 10 ) /= 0 dodet = mod ( job, 100 ) / 10 /= 0 doert = mod ( job, 1000 ) / 100 /= 0 if ( dodet .or. doert ) then if ( doert ) then inert(1) = 0 inert(2) = 0 inert(3) = 0 end if if ( dodet ) then det(1) = 1.0E+00 det(2) = 0.0E+00 end if t = 0.0E+00 ik = 0 do k = 1, n kk = ik + k d = ap(kk) ! ! 2 by 2 block ! use det (d s) = (d/t * c - t) * t, t = abs ( s ) ! (s c) ! to avoid underflow/overflow troubles. ! ! Take two passes through scaling. Use T for flag. ! if ( kpvt(k) <= 0 ) then if ( t == 0.0E+00 ) then ikp1 = ik + k kkp1 = ikp1 + k t = abs ( ap(kkp1) ) d = ( d / t ) * ap(kkp1+1) - t else d = t t = 0.0E+00 end if end if if ( doert ) then if ( d > 0.0E+00 ) then inert(1) = inert(1) + 1 else if ( d < 0.0E+00 ) then inert(2) = inert(2) + 1 else if ( d == 0.0E+00 ) then inert(3) = inert(3) + 1 end if end if if ( dodet ) then det(1) = d * det(1) if ( det(1) /= 0.0E+00 ) then do while ( abs ( det(1) ) < 1.0E+00 ) det(1) = ten * det(1) det(2) = det(2) - 1.0E+00 end do do while ( abs ( det(1) ) >= ten ) det(1) = det(1) / ten det(2) = det(2) + 1.0E+00 end do end if end if ik = ik + k end do end if ! ! Compute inverse(A). ! if ( doinv ) then k = 1 ik = 0 do while ( k <= n ) km1 = k - 1 kk = ik + k ikp1 = ik + k kkp1 = ikp1 + k if ( kpvt(k) >= 0) then ! ! 1 by 1. ! ap(kk) = 1.0E+00 / ap(kk) if ( k-1 >= 1 ) then call scopy ( k-1, ap(ik+1), 1, work, 1 ) ij = 0 do j = 1, k-1 jk = ik + j ap(jk) = sdot ( j, ap(ij+1), 1, work, 1 ) call saxpy ( j-1, work(j), ap(ij+1), 1, ap(ik+1), 1 ) ij = ij + j end do ap(kk) = ap(kk) + sdot ( k-1, work, 1, ap(ik+1), 1 ) end if kstep = 1 else ! ! 2 by 2. ! t = abs ( ap(kkp1) ) ak = ap(kk) / t akp1 = ap(kkp1+1) / t akkp1 = ap(kkp1) / t d = t * ( ak * akp1 - 1.0E+00 ) ap(kk) = akp1 / d ap(kkp1+1) = ak / d ap(kkp1) = -akkp1 / d if ( km1 >= 1 ) then call scopy ( km1, ap(ikp1+1), 1, work, 1 ) ij = 0 do j = 1, km1 jkp1 = ikp1 + j ap(jkp1) = sdot ( j, ap(ij+1), 1, work, 1 ) call saxpy ( j-1, work(j), ap(ij+1), 1, ap(ikp1+1), 1 ) ij = ij + j end do ap(kkp1+1) = ap(kkp1+1) + sdot ( km1, work, 1, ap(ikp1+1), 1 ) ap(kkp1) = ap(kkp1) + sdot ( km1, ap(ik+1), 1, ap(ikp1+1), 1 ) call scopy ( km1, ap(ik+1), 1, work, 1 ) ij = 0 do j = 1, km1 jk = ik + j ap(jk) = sdot ( j, ap(ij+1), 1, work, 1 ) call saxpy ( j-1, work(j), ap(ij+1), 1, ap(ik+1), 1 ) ij = ij + j end do ap(kk) = ap(kk) + sdot ( km1, work, 1, ap(ik+1), 1 ) end if kstep = 2 end if ! ! Swap. ! ks = abs ( kpvt(k) ) if ( ks /= k ) then iks = ( ks * ( ks - 1 ) ) / 2 call sswap ( ks, ap(iks+1), 1, ap(ik+1), 1 ) ksj = ik + ks do jb = ks, k j = k + ks - jb jk = ik + j call r_swap ( ap(jk), ap(ksj) ) ksj = ksj - (j - 1) end do if ( kstep /= 1 ) then kskp1 = ikp1 + ks call r_swap ( ap(kskp1), ap(kkp1) ) end if end if ik = ik + k if (kstep == 2) ik = ik + k + 1 k = k + kstep end do end if return end subroutine sspfa ( ap, n, kpvt, info ) ! !******************************************************************************* ! !! SSPFA factors a real symmetric matrix stored in packed form. ! ! ! Discussion: ! ! To solve A*X = B, follow SSPFA by SSPSL. ! ! To compute inverse(A)*C, follow SSPFA by SSPSL. ! ! To compute determinant(A), follow SSPFA by SSPDI. ! ! To compute inertia(A), follow SSPFA by SSPDI. ! ! To compute inverse(A), follow SSPFA by SSPDI. ! ! Packed storage: ! ! The following program segment will pack the upper triangle of a ! symmetric matrix. ! ! k = 0 ! do j = 1, n ! do i = 1, j ! k = k + 1 ! ap(k) = a(i,j) ! end do ! end do ! ! Reference: ! ! Dongarra, Moler, Bunch and Stewart, ! LINPACK User's Guide, ! SIAM, (Society for Industrial and Applied Mathematics), ! 3600 University City Science Center, ! Philadelphia, PA, 19104-2688. ! ISBN 0-89871-172-X ! ! Parameters: ! ! Input, real AP(N*(N+1)/2). On input, the packed form of a symmetric ! matrix A. The columns of the upper triangle are stored sequentially ! in a one-dimensional array. On output, a block diagonal matrix and ! the multipliers which were used to obtain it stored in packed form. ! The factorization can be written A = U*D*U' where U is a product of ! permutation and unit upper triangular matrices, U' is the transpose ! of U, and D is block diagonal with 1 by 1 and 2 by 2 blocks. ! ! Input, integer N, the order of the matrix. ! ! Output, integer KPVT(N), the pivot indices. ! ! Output, integer INFO, error flag. ! 0, normal value. ! K, if the K-th pivot block is singular. This is not an error ! condition for this subroutine, but it does indicate that SSPSL or ! SSPDI may divide by zero if called. ! implicit none ! integer n ! real absakk real ak real akm1 real alpha real ap((n*(n+1))/2) real bk real bkm1 real colmax real denom integer ij integer ijj integer ik integer ikm1 integer im integer imax integer imaxp1 integer imim integer imj integer imk integer info integer isamax integer j integer jj integer jk integer jkm1 integer jmax integer jmim integer k integer kk integer km1 integer km1k integer km1km1 integer kpvt(n) integer kstep real mulk real mulkm1 real rowmax logical swap real t ! ! ALPHA is used in choosing pivot block size. ! alpha = ( 1.0E+00 + sqrt ( 17.0E+00 ) ) / 8.0E+00 info = 0 ! ! Main loop on K, which goes from N to 1. ! k = n ik = (n*(n - 1))/2 10 continue ! ! Leave the loop if K=0 or K=1. ! if ( k == 0 ) then return end if if ( k == 1 ) then kpvt(1) = 1 if ( ap(1) == 0.0E+00 ) then info = 1 end if return end if ! ! This section of code determines the kind of elimination to be performed. ! When it is completed, KSTEP will be set to the size of the pivot block, ! and SWAP will be set to .true. if an interchange is required. ! km1 = k - 1 kk = ik + k absakk = abs ( ap(kk) ) ! ! Determine the largest off-diagonal element in column K. ! imax = isamax ( k-1, ap(ik+1), 1 ) imk = ik + imax colmax = abs ( ap(imk) ) if ( absakk >= alpha * colmax ) then kstep = 1 swap = .false. ! ! Determine the largest off-diagonal element in row IMAX. ! else rowmax = 0.0E+00 imaxp1 = imax + 1 im = imax * ( imax - 1 ) / 2 imj = im + 2 * imax do j = imaxp1, k rowmax = max ( rowmax, abs ( ap(imj) ) ) imj = imj + j end do if ( imax /= 1 ) then jmax = isamax ( imax-1, ap(im+1), 1 ) jmim = jmax + im rowmax = max ( rowmax, abs ( ap(jmim) ) ) end if imim = imax + im if ( abs ( ap(imim) ) >= alpha * rowmax ) then kstep = 1 swap = .true. else if ( absakk >= alpha * colmax * ( colmax / rowmax ) ) then kstep = 1 swap = .false. else kstep = 2 swap = imax /= km1 end if end if ! ! Column K is zero. Set INFO and iterate the loop. ! if ( max ( absakk, colmax ) == 0.0E+00 ) then kpvt(k) = k info = k go to 190 end if if ( kstep /= 2 ) then ! ! 1 x 1 pivot block. ! if ( swap ) then ! ! Perform an interchange. ! call sswap ( imax, ap(im+1), 1, ap(ik+1), 1 ) imj = ik + imax do jj = imax, k j = k + imax - jj jk = ik + j t = ap(jk) ap(jk) = ap(imj) ap(imj) = t imj = imj - (j - 1) end do end if ! ! Perform the elimination. ! ij = ik - (k - 1) do jj = 1, km1 j = k - jj jk = ik + j mulk = -ap(jk)/ap(kk) t = mulk call saxpy ( j, t, ap(ik+1), 1, ap(ij+1), 1 ) ijj = ij + j ap(jk) = mulk ij = ij - (j - 1) end do ! ! Set the pivot array. ! if ( swap ) then kpvt(k) = imax else kpvt(k) = k end if else ! ! 2 x 2 pivot block. ! km1k = ik + k - 1 ikm1 = ik - (k - 1) ! ! Perform an interchange. ! if ( swap ) then call sswap ( imax, ap(im+1), 1, ap(ikm1+1), 1 ) imj = ikm1 + imax do jj = imax, km1 j = km1 + imax - jj jkm1 = ikm1 + j t = ap(jkm1) ap(jkm1) = ap(imj) ap(imj) = t imj = imj - ( j - 1 ) end do t = ap(km1k) ap(km1k) = ap(imk) ap(imk) = t end if ! ! Perform the elimination. ! if ( k-2 /= 0 ) then ak = ap(kk) / ap(km1k) km1km1 = ikm1 + k - 1 akm1 = ap(km1km1) / ap(km1k) denom = 1.0E+00 - ak * akm1 ij = ik - ( k - 1 ) - ( k - 2 ) do jj = 1, k-2 j = km1 - jj jk = ik + j bk = ap(jk) / ap(km1k) jkm1 = ikm1 + j bkm1 = ap(jkm1) / ap(km1k) mulk = ( akm1 * bk - bkm1 ) / denom mulkm1 = ( ak * bkm1 - bk ) / denom t = mulk call saxpy ( j, t, ap(ik+1), 1, ap(ij+1), 1 ) t = mulkm1 call saxpy ( j, t, ap(ikm1+1), 1, ap(ij+1), 1 ) ap(jk) = mulk ap(jkm1) = mulkm1 ijj = ij + j ij = ij - (j - 1) end do end if ! ! Set the pivot array. ! kpvt(k) = 1 - k if (swap) kpvt(k) = -imax kpvt(k-1) = kpvt(k) end if 190 continue ik = ik - (k - 1) if ( kstep == 2 ) then ik = ik - (k - 2) end if k = k - kstep go to 10 end subroutine sspsl ( ap, n, kpvt, b ) ! !******************************************************************************* ! !! SSPSL solves the real symmetric system factored by SSPFA. ! ! ! Discussion: ! ! To compute inverse(A) * C where C is a matrix with P columns: ! ! call sspfa ( ap, n, kpvt, info ) ! ! if ( info /= 0 ) go to ... ! ! do j = 1, p ! call sspsl ( ap, n, kpvt, c(1,j) ) ! end do ! ! A division by zero may occur if SSPCO has set RCOND == 0.0E+00 ! or SSPFA has set INFO /= 0. ! ! Reference: ! ! Dongarra, Moler, Bunch and Stewart, ! LINPACK User's Guide, ! SIAM, (Society for Industrial and Applied Mathematics), ! 3600 University City Science Center, ! Philadelphia, PA, 19104-2688. ! ISBN 0-89871-172-X ! ! Parameters: ! ! Input, real AP(N*(N+1)/2), the output from SSPFA. ! ! Input, integer N, the order of the matrix. ! ! Input, integer KPVT(N), the pivot vector from SSPFA. ! ! Input/output, real B(N). On input, the right hand side. ! On output, the solution. ! implicit none ! integer n ! real ak real akm1 real ap((n*(n+1))/2) real b(n) real bk real bkm1 real denom integer ik integer ikm1 integer ikp1 integer k integer kk integer km1k integer km1km1 integer kp integer kpvt(n) real sdot real temp ! ! Loop backward applying the transformations and D inverse to B. ! k = n ik = (n*(n - 1))/2 do while ( k > 0 ) kk = ik + k if ( kpvt(k) >= 0 ) then ! ! 1 x 1 pivot block. ! if (k /= 1) then kp = kpvt(k) ! ! Interchange. ! if ( kp /= k ) then call r_swap ( b(k), b(kp) ) end if ! ! Apply the transformation. ! call saxpy ( k-1, b(k), ap(ik+1), 1, b(1), 1 ) end if ! ! Apply D inverse. ! b(k) = b(k) / ap(kk) k = k - 1 ik = ik - k else ! ! 2 x 2 pivot block. ! ikm1 = ik - (k - 1) if ( k /= 2 ) then kp = abs ( kpvt(k) ) ! ! Interchange. ! if ( kp /= k-1 ) then call r_swap ( b(k-1), b(kp) ) end if ! ! Apply the transformation. ! call saxpy ( k-2, b(k), ap(ik+1), 1, b(1), 1 ) call saxpy ( k-2, b(k-1), ap(ikm1+1), 1, b(1), 1 ) end if ! ! Apply D inverse. ! km1k = ik + k - 1 kk = ik + k ak = ap(kk) / ap(km1k) km1km1 = ikm1 + k - 1 akm1 = ap(km1km1) / ap(km1k) bk = b(k) / ap(km1k) bkm1 = b(k-1) / ap(km1k) denom = ak * akm1 - 1.0E+00 b(k) = ( akm1 * bk - bkm1 ) / denom b(k-1) = ( ak * bkm1 - bk ) / denom k = k - 2 ik = ik - (k + 1) - k end if end do ! ! Loop forward applying the transformations. ! k = 1 ik = 0 do while ( k <= n ) if ( kpvt(k) >= 0 ) then ! ! 1 x 1 pivot block. ! if ( k /= 1 ) then ! ! Apply the transformation. ! b(k) = b(k) + sdot ( k-1, ap(ik+1), 1, b(1), 1 ) kp = kpvt(k) ! ! Interchange. ! if ( kp /= k ) then call r_swap ( b(k), b(kp) ) end if end if ik = ik + k k = k + 1 else ! ! 2 x 2 pivot block. ! if ( k /= 1 ) then ! ! Apply the transformation. ! b(k) = b(k) + sdot ( k-1, ap(ik+1), 1, b(1), 1 ) ikp1 = ik + k b(k+1) = b(k+1) + sdot ( k-1, ap(ikp1+1), 1, b(1), 1 ) kp = abs ( kpvt(k) ) ! ! Interchange. ! if ( kp /= k ) then call r_swap ( b(k), b(kp) ) end if end if ik = ik + k + k + 1 k = k + 2 end if end do return end subroutine ssvdc ( x, ldx, n, p, s, e, u, ldu, v, ldv, work, job, info ) ! !******************************************************************************* ! !! SSVDC computes the singular value decomposition of a real rectangular matrix. ! ! ! Discussion: ! ! SSVDC reduces a real N by P matrix X to diagonal form by orthogonal ! transformations U and V. The diagonal elements S(I) are the singular ! values of X. The columns of U are the corresponding left singular ! vectors, and the columns of V the right singular vectors. ! ! Reference: ! ! Dongarra, Moler, Bunch and Stewart, ! LINPACK User's Guide, ! SIAM, (Society for Industrial and Applied Mathematics), ! 3600 University City Science Center, ! Philadelphia, PA, 19104-2688. ! ISBN 0-89871-172-X ! ! Parameters: ! ! Input/output, real X(LDX,P). On input, the matrix whose singular value ! decomposition is to be computed. On output, the matrix has been ! destroyed. ! ! Input, integer LDX, the leading dimension of the array X. ! LDX must be at least N. ! ! Input, integer N, the number of rows of the matrix. ! ! Input, integer P, the number of columns of the matrix X. ! ! Output, real S(MM), where MM = min(N+1,P). The first min(N,P) entries ! of S contain the singular values of X arranged in descending ! order of magnitude. ! ! Output, real E(P), ordinarily contains zeros. However see the ! discussion of INFO for exceptions. ! ! Output, real U(LDU,K). If JOBA = 1 then K = N; if JOBA >= 2 then ! K = min(N,P). U contains the matrix of left singular vectors. ! U is not referenced if JOBA = 0. If N <= P or if JOBA = 2, then ! U may be identified with X in the subroutine call. ! ! Input, integer LDU, the leading dimension of the array U. ! LDU must be at least N. ! ! Output, real V(LDV,P), the matrix of right singular vectors. ! V is not referenced if JOB is 0. If P <= N, then V may be identified ! with X in the subroutine call. ! ! Input, integer LDV, the leading dimension of the array V. ! LDV must be at least P. ! ! Workspace, real WORK(N). ! ! Input, integer JOB, controls the computation of the singular ! vectors. It has the decimal expansion AB with the following meaning: ! ! A = 0, do not compute the left singular vectors. ! A = 1, return the N left singular vectors in U. ! A >= 2, return the first min(N,P) singular vectors in U. ! B = 0, do not compute the right singular vectors. ! B = 1, return the right singular vectors in V. ! ! Output, integer INFO, status indicator. ! The singular values (and their corresponding singular vectors) ! S(INFO+1), S(INFO+2),...,S(M) are correct (here M=min(N,P)). ! Thus if INFO is 0, all the singular values and their vectors are ! correct. In any event, the matrix B = U'*X*V is the bidiagonal matrix ! with the elements of S on its diagonal and the elements of E on ! its super-diagonal (U' is the transpose of U). Thus the singular ! values of X and B are the same. ! implicit none ! integer ldu integer ldv integer ldx integer n integer p ! real b real c real cs real e(p) real el real emm1 real f real g integer i integer info integer iter integer j integer job integer jobu integer k integer kase integer kk integer l integer ll integer lls integer lm1 integer ls integer lu integer m integer, parameter :: maxit = 30 integer mm integer mm1 integer mp1 integer nct integer nctp1 integer ncu integer nrt integer nrtp1 real s(*) real scale real sdot real shift real sl real sm real smm1 real sn real snrm2 real t real t1 real test real u(ldu,*) real v(ldv,p) logical wantu logical wantv real work(n) real x(ldx,p) real ztest ! ! Determine what is to be computed. ! wantu = .false. wantv = .false. jobu = mod ( job, 100 ) / 10 if (jobu > 1) then ncu = min ( n, p ) else ncu = n end if if ( jobu /= 0 ) then wantu = .true. end if if ( mod ( job, 10 ) /= 0 ) then wantv = .true. end if ! ! Reduce X to bidiagonal form, storing the diagonal elements ! in S and the super-diagonal elements in E. ! info = 0 nct = min ( n-1, p ) nrt = max ( 0, min ( p-2, n ) ) lu = max ( nct, nrt ) do l = 1, lu ! ! Compute the transformation for the L-th column and ! place the L-th diagonal in S(L). ! if ( l <= nct ) then s(l) = snrm2 ( n-l+1, x(l,l), 1 ) if ( s(l) /= 0.0E+00 ) then if (x(l,l) /= 0.0E+00 ) then s(l) = sign ( s(l), x(l,l) ) end if call sscal ( n-l+1, 1.0E+00 / s(l), x(l,l), 1 ) x(l,l) = 1.0E+00 + x(l,l) end if s(l) = -s(l) end if do j = l+1, p ! ! Apply the transformation. ! if ( l <= nct .and. s(l) /= 0.0E+00 ) then t = - sdot ( n-l+1, x(l,l), 1, x(l,j), 1 ) / x(l,l) call saxpy ( n-l+1, t, x(l,l), 1, x(l,j), 1 ) end if ! ! Place the L-th row of X into E for the ! subsequent calculation of the row transformation. ! e(j) = x(l,j) end do ! ! Place the transformation in U for subsequent back multiplication. ! if ( wantu .and. l <= nct ) then u(l:n,l) = x(l:n,l) end if if ( l <= nrt ) then ! ! Compute the L-th row transformation and place the ! L-th superdiagonal in E(L). ! e(l) = snrm2 ( p-l, e(l+1), 1 ) if ( e(l) /= 0.0E+00 ) then if ( e(l+1) /= 0.0E+00 ) then e(l) = sign ( e(l), e(l+1) ) end if call sscal ( p-l, 1.0E+00 / e(l), e(l+1), 1 ) e(l+1) = 1.0E+00 + e(l+1) end if e(l) = -e(l) ! ! Apply the transformation. ! if ( l+1 <= n .and. e(l) /= 0.0E+00 ) then work(l+1:n) = 0.0E+00 do j = l+1, p call saxpy ( n-l, e(j), x(l+1,j), 1, work(l+1), 1 ) end do do j = l+1, p call saxpy ( n-l, -e(j)/e(l+1), work(l+1), 1, x(l+1,j), 1 ) end do end if ! ! Place the transformation in V for subsequent back multiplication. ! if ( wantv ) then v(l+1:p,l) = e(l+1:p) end if end if end do ! ! Set up the final bidiagonal matrix of order M. ! m = min ( p, n+1 ) nctp1 = nct + 1 nrtp1 = nrt + 1 if ( nct < p ) then s(nctp1) = x(nctp1,nctp1) end if if ( n < m ) then s(m) = 0.0E+00 end if if ( nrtp1 < m ) then e(nrtp1) = x(nrtp1,m) end if e(m) = 0.0E+00 ! ! If required, generate U. ! if ( wantu ) then u(1:n,nctp1:ncu) = 0.0E+00 do j = nctp1, ncu u(j,j) = 1.0E+00 end do do ll = 1, nct l = nct - ll + 1 if ( s(l) /= 0.0E+00 ) then do j = l+1, ncu t = - sdot ( n-l+1, u(l,l), 1, u(l,j), 1 ) / u(l,l) call saxpy ( n-l+1, t, u(l,l), 1, u(l,j), 1 ) end do call sscal ( n-l+1, -1.0E+00, u(l,l), 1 ) u(l,l) = 1.0E+00 + u(l,l) u(1:l-1,l) = 0.0E+00 else u(1:n,l) = 0.0E+00 u(l,l) = 1.0E+00 end if end do end if ! ! If it is required, generate V. ! if ( wantv ) then do ll = 1, p l = p - ll + 1 if ( l <= nrt .and. e(l) /= 0.0E+00 ) then do j = l+1, p t = - sdot ( p-l, v(l+1,l), 1, v(l+1,j), 1 ) / v(l+1,l) call saxpy ( p-l, t, v(l+1,l), 1, v(l+1,j), 1 ) end do end if v(1:p,l) = 0.0E+00 v(l,l) = 1.0E+00 end do end if ! ! Main iteration loop for the singular values. ! mm = m iter = 0 do while ( m > 0 ) ! ! If too many iterations have been performed, set flag and return. ! if ( iter >= maxit ) then info = m return end if ! ! This section of the program inspects for ! negligible elements in the S and E arrays. ! ! On completion the variables KASE and L are set as follows: ! ! KASE = 1 if S(M) and E(L-1) are negligible and L= 0 ) then ix = 1 else ix = ( - n + 1 ) * incx + 1 end if if ( incy >= 0 ) then iy = 1 else iy = ( - n + 1 ) * incy + 1 end if do i = 1, n stemp = x(ix) x(ix) = y(iy) y(iy) = stemp ix = ix + incx iy = iy + incy end do end if return end subroutine strco ( t, ldt, n, rcond, z, job ) ! !******************************************************************************* ! !! STRCO estimates the condition of a real triangular matrix. ! ! ! Reference: ! ! Dongarra, Moler, Bunch and Stewart, ! LINPACK User's Guide, ! SIAM, (Society for Industrial and Applied Mathematics), ! 3600 University City Science Center, ! Philadelphia, PA, 19104-2688. ! ISBN 0-89871-172-X ! ! Parameters: ! ! Input, real T(LDT,N), the triangular matrix. The zero elements of ! the matrix are not referenced, and the corresponding elements of the ! array can be used to store other information. ! ! Input, integer LDT, the leading dimension of the array T. ! ! Input, integer N, the order of the matrix. ! ! Input, integer JOB, indicates the shape of T: ! 0, T is lower triangular. ! nonzero, T is upper triangular. ! ! Output, real RCOND, an estimate of the reciprocal condition of T. ! For the system T*X = B, relative perturbations in T and B of size ! EPSILON may cause relative perturbations in X of size EPSILON/RCOND. ! If RCOND is so small that the logical expression ! 1.0E+00 + RCOND == 1.0E+00 ! is true, then T may be singular to working precision. In particular, ! RCOND is zero if exact singularity is detected or the estimate underflows. ! ! Output, real Z(N) a work vector whose contents are usually unimportant. ! If T is close to a singular matrix, then Z is an approximate null vector ! in the sense that ! norm(A*Z) = RCOND * norm(A) * norm(Z). ! implicit none ! integer ldt integer n ! real ek integer i1 integer j integer j1 integer j2 integer job integer k integer kk integer l logical lower real rcond real s real sm real t(ldt,n) real tnorm real w real wk real wkm real ynorm real z(n) ! lower = job == 0 ! ! Compute the 1-norm of T. ! tnorm = 0.0E+00 do j = 1, n if ( lower ) then l = n + 1 - j i1 = j else l = j i1 = 1 end if tnorm = max ( tnorm, sum ( abs ( t(i1:i1+l-1,j) ) ) ) end do ! ! RCOND = 1/(norm(T)*(estimate of norm(inverse(T)))). ! ! Estimate = norm(Z)/norm(Y) where T*Z = Y and T'*Y = E. ! ! T' is the transpose of T. ! ! The components of E are chosen to cause maximum local ! growth in the elements of Y. ! ! The vectors are frequently rescaled to avoid overflow. ! ! Solve T'*Y = E. ! ek = 1.0E+00 z(1:n) = 0.0E+00 do kk = 1, n if ( lower ) then k = n + 1 - kk else k = kk end if if ( z(k) /= 0.0E+00 ) then ek = sign ( ek, -z(k) ) end if if ( abs ( ek-z(k) ) > abs ( t(k,k) ) ) then s = abs ( t(k,k) ) / abs ( ek - z(k) ) z(1:n) = s * z(1:n) ek = s * ek end if wk = ek - z(k) wkm = -ek - z(k) s = abs ( wk ) sm = abs ( wkm ) if ( t(k,k) /= 0.0E+00 ) then wk = wk / t(k,k) wkm = wkm / t(k,k) else wk = 1.0E+00 wkm = 1.0E+00 end if if ( kk /= n ) then if ( lower ) then j1 = 1 j2 = k - 1 else j1 = k + 1 j2 = n end if do j = j1, j2 sm = sm + abs ( z(j) + wkm * t(k,j) ) z(j) = z(j) + wk * t(k,j) s = s + abs ( z(j) ) end do if ( s < sm ) then w = wkm - wk wk = wkm z(j1:j2) = z(j1:j2) + w * t(k,j1:j2) end if end if z(k) = wk end do z(1:n) = z(1:n) / sum ( abs ( z(1:n) ) ) ynorm = 1.0E+00 ! ! Solve T*Z = Y. ! do kk = 1, n if ( lower ) then k = kk else k = n + 1 - kk end if if ( abs ( z(k) ) > abs ( t(k,k) ) ) then s = abs ( t(k,k) ) / abs ( z(k) ) z(1:n) = s * z(1:n) ynorm = s * ynorm end if if ( t(k,k) /= 0.0E+00 ) then z(k) = z(k) / t(k,k) else z(k) = 1.0E+00 end if if ( lower ) then i1 = k + 1 else i1 = 1 end if if ( kk < n ) then w = -z(k) call saxpy ( n-kk, w, t(i1,k), 1, z(i1), 1 ) end if end do ! ! Make ZNORM = 1.0. ! s = 1.0E+00 / sum ( abs ( z(1:n) ) ) z(1:n) = s * z(1:n) ynorm = s * ynorm if ( tnorm /= 0.0E+00 ) then rcond = ynorm / tnorm else rcond = 0.0E+00 end if return end subroutine strdi ( t, ldt, n, det, job, info ) ! !*********************************************************************** ! !! STRDI computes the determinant and inverse of a real triangular matrix. ! ! ! Modified: ! ! 07 March 2001 ! ! Reference: ! ! Dongarra, Moler, Bunch and Stewart, ! LINPACK User's Guide, ! SIAM, (Society for Industrial and Applied Mathematics), ! 3600 University City Science Center, ! Philadelphia, PA, 19104-2688. ! ISBN 0-89871-172-X ! ! Parameters: ! ! Input/output, real T(LDT,N). ! On input, T contains the triangular matrix. The zero elements of the ! matrix are not referenced, and the corresponding elements of the array ! can be used to store other information. ! On output, T contains the inverse matrix, if it was requested. ! ! Input, integer LDT, the leading dimension of T. ! ! Input, integer N, the order of the matrix. ! ! Output, real DET(2), the determinant of the matrix, if requested. ! The determinant = DET(1) * 10.0**DET(2), with 1.0 <= abs ( DET(1) ) < 10.0, ! or DET(1) == 0. ! ! Input, integer JOB, specifies the shape of T, and the task. ! 010, inverse of lower triangular matrix. ! 011, inverse of upper triangular matrix. ! 100, determinant only. ! 110, determinant and inverse of lower triangular. ! 111, determinant and inverse of upper triangular. ! ! Output, integer INFO. ! If the inverse was requested, then ! 0, if the system was nonsingular; ! nonzero, if the system was singular. ! implicit none ! integer ldt integer n ! real det(2) integer i integer info integer j integer job integer k real t(ldt,n) real temp real, parameter :: ten = 10.0E+00 ! ! Determinant. ! if ( job / 100 /= 0 ) then det(1) = 1.0E+00 det(2) = 0.0E+00 do i = 1, n det(1) = t(i,i) * det(1) if ( det(1) == 0.0E+00 ) then exit end if do while ( abs ( det(1) ) < 1.0E+00 ) det(1) = ten * det(1) det(2) = det(2) - 1.0E+00 end do do while ( abs ( det(1) ) >= ten ) det(1) = det(1) / ten det(2) = det(2) + 1.0E+00 end do end do end if if ( mod ( job / 10, 10 ) == 0 ) then return end if ! ! Inverse of an upper triangular matrix. ! if ( mod ( job,10 ) /= 0 ) then info = 0 do k = 1, n if ( t(k,k) == 0.0E+00 ) then info = k exit end if t(k,k) = 1.0E+00 / t(k,k) temp = -t(k,k) call sscal ( k-1, temp, t(1,k), 1 ) do j = k + 1, n temp = t(k,j) t(k,j) = 0.0E+00 call saxpy ( k, temp, t(1,k), 1, t(1,j), 1 ) end do end do ! ! Inverse of a lower triangular matrix. ! else info = 0 do k = n, 1, -1 if ( t(k,k) == 0.0E+00 ) then info = k exit end if t(k,k) = 1.0E+00 / t(k,k) temp = -t(k,k) if ( k /= n ) then call sscal ( n-k, temp, t(k+1,k), 1 ) end if do j = 1, k-1 temp = t(k,j) t(k,j) = 0.0E+00 call saxpy ( n-k+1, temp, t(k,k), 1, t(k,j), 1 ) end do end do end if return end subroutine strsl ( t, ldt, n, b, job, info ) ! !******************************************************************************* ! !! STRSL solves triangular linear systems. ! ! ! Discussion: ! ! STRSL can solve T * X = B or T' * X = B where T is a triangular ! matrix of order N. ! ! Here T' denotes the transpose of the matrix T. ! ! Reference: ! ! Dongarra, Moler, Bunch and Stewart, ! LINPACK User's Guide, ! SIAM, (Society for Industrial and Applied Mathematics), ! 3600 University City Science Center, ! Philadelphia, PA, 19104-2688. ! ISBN 0-89871-172-X ! ! Parameters: ! ! Input, real T(LDT,N), the matrix of the system. The zero elements of ! the matrix are not referenced, and the corresponding elements of the ! array can be used to store other information. ! ! Input, integer LDT, the leading dimension of the array T. ! ! Input, integer N, the order of the matrix. ! ! Input/output, real B(N). On input, the right hand side. ! On output, the solution. ! ! Input, integer JOB, specifies what kind of system is to be solved: ! ! 00, solve T * X = B, T lower triangular, ! 01, solve T * X = B, T upper triangular, ! 10, solve T'* X = B, T lower triangular, ! 11, solve T'* X = B, T upper triangular. ! ! Output, integer INFO, singularity indicator. ! 0, the system is nonsingular. ! nonzero, the index of the first zero diagonal element of T. ! implicit none ! integer ldt integer n ! real b(n) integer case integer info integer j integer jj integer job real sdot real t(ldt,n) real temp ! ! Check for zero diagonal elements. ! do j = 1, n if ( t(j,j) == 0.0E+00 ) then info = j return end if end do info = 0 ! ! Determine the task and go to it. ! if ( mod ( job, 10 ) == 0 ) then case = 1 else case = 2 end if if ( mod ( job, 100 ) / 10 /= 0 ) then case = case + 2 end if ! ! Solve T*X=B for T lower triangular. ! if ( case == 1 ) then b(1) = b(1) / t(1,1) do j = 2, n temp = -b(j-1) call saxpy ( n-j+1, temp, t(j,j-1), 1, b(j), 1 ) b(j) = b(j) / t(j,j) end do ! ! Solve T*X=B for T upper triangular. ! else if ( case == 2 ) then b(n) = b(n) / t(n,n) do jj = 2, n j = n - jj + 1 temp = -b(j+1) call saxpy ( j, temp, t(1,j+1), 1, b(1), 1 ) b(j) = b(j) / t(j,j) end do ! ! Solve T'*X=B for T lower triangular. ! else if ( case == 3 ) then b(n) = b(n) / t(n,n) do jj = 2, n j = n - jj + 1 b(j) = b(j) - sdot ( jj-1, t(j+1,j), 1, b(j+1), 1 ) b(j) = b(j) / t(j,j) end do ! ! Solve T'*X=B for T upper triangular. ! else if ( case == 4 ) then b(1) = b(1) / t(1,1) do j = 2, n b(j) = b(j) - sdot ( j-1, t(1,j), 1, b(1), 1 ) b(j) = b(j) / t(j,j) end do end if return end subroutine timestamp ( ) ! !******************************************************************************* ! !! TIMESTAMP prints the current YMDHMS date as a time stamp. ! ! ! Example: ! ! May 31 2001 9:45:54.872 AM ! ! Modified: ! ! 31 May 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! None ! character ( len = 8 ) ampm integer d character ( len = 8 ) date 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 character ( len = 10 ) time integer values(8) integer 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