Mercurial > octave
view libcruft/npsol/lschol.f @ 2329:30c606bec7a8
[project @ 1996-07-19 01:29:05 by jwe]
Initial revision
author | jwe |
---|---|
date | Fri, 19 Jul 1996 01:29:55 +0000 |
parents | |
children |
line wrap: on
line source
*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ SUBROUTINE LSCHOL( NROWH, N, NRANK, TOLRNK, KX, H, INFORM ) IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER KX(*) DOUBLE PRECISION H(NROWH,*) ************************************************************************ * LSCHOL forms the Cholesky factorization of the positive * semi-definite matrix H such that * PHP' = R'R * where P is a permutation matrix and R is upper triangular. * The permutation P is chosen to maximize the diagonal of R at each * stage. Only the diagonal and super-diagonal elements of H are * used. * * Output: * * INFORM = 0 the factorization was computed successfully, * with the Cholesky factor written in the upper * triangular part of H and P stored in KX. * 1 the matrix H was indefinite. * * Systems Optimization Laboratory, Stanford University. * Original version of LSCHOL dated 2-February-1981. * Level 2 Blas added 29-June-1986. * This version of LSCHOL dated 30-June-1986. ************************************************************************ COMMON /SOL1CM/ NOUT INTRINSIC ABS , MAX , SQRT EXTERNAL IDAMAX PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) INFORM = 0 NRANK = 0 * Main loop for computing rows of R. DO 200 J = 1, N * Find maximum available diagonal. KMAX = J - 1 + IDAMAX( N-J+1, H(J,J), NROWH+1 ) DMAX = H(KMAX,KMAX) IF (DMAX .LE. TOLRNK*ABS(H(1,1))) GO TO 300 * Perform a symmetric interchange if necessary. IF (KMAX .NE. J) THEN K = KX(KMAX) KX(KMAX) = KX(J) KX(J) = K CALL DSWAP ( J , H(1,J) , 1, H(1,KMAX), 1 ) CALL DSWAP ( KMAX-J+1, H(J,KMAX), 1, H(J,J) , NROWH ) CALL DSWAP ( N-KMAX+1, H(KMAX,KMAX), NROWH, $ H(J,KMAX) , NROWH ) H(KMAX,KMAX) = H(J,J) END IF * Set the diagonal of R. D = SQRT( DMAX ) H(J,J) = D NRANK = NRANK + 1 IF (J .LT. N) THEN * Set the super-diagonal elements of this row of R and update * the elements of the block that is yet to be factorized. CALL DSCAL ( N-J, (ONE/D), H(J ,J+1), NROWH ) CALL DSYR ( 'U', N-J, -ONE, H(J ,J+1), NROWH, $ H(J+1,J+1), NROWH ) END IF 200 CONTINUE * ------------------------------------------------------------------ * Check for the semi-definite case. * ------------------------------------------------------------------ 300 IF (NRANK .LT. N) THEN * Find the largest element in the unfactorized block. SUPMAX = ZERO DO 310 I = J, N-1 K = I + IDAMAX( N-I, H(I,I+1), NROWH ) SUPMAX = MAX( SUPMAX, ABS(H(I,K)) ) 310 CONTINUE IF (SUPMAX .GT. TOLRNK*ABS(H(1,1))) THEN WRITE (NOUT, 1000) DMAX, SUPMAX INFORM = 1 END IF END IF RETURN 1000 FORMAT(' XXX Hessian appears to be indefinite.' $ /' XXX Maximum diagonal and off-diagonal ignored', $ ' in the Cholesky factorization:', 1P2E22.14 ) * End of LSCHOL. END