* F07GFF Example Program Text * Mark 21 Release. NAG Copyright 2004. * .. Parameters .. INTEGER NIN, NOUT PARAMETER (NIN=5,NOUT=6) INTEGER NMAX PARAMETER (NMAX=8) CHARACTER UPLO PARAMETER (UPLO='U') * .. Local Scalars .. DOUBLE PRECISION AMAX, BIG, SCOND, SJ, SMALL INTEGER I, IFAIL, INFO, J, JJ, N * .. Local Arrays .. DOUBLE PRECISION AP((NMAX*(NMAX+1))/2), S(NMAX) * .. External Functions .. DOUBLE PRECISION X02AJF, X02AMF INTEGER X02BHF EXTERNAL X02AJF, X02AMF, X02BHF * .. External Subroutines .. EXTERNAL DPPEQU, X04CCF * .. Executable Statements .. WRITE (NOUT,*) 'F07GFF Example Program Results' WRITE (NOUT,*) * Skip heading in data file READ (NIN,*) READ (NIN,*) N IF (N.LE.NMAX) THEN * * Read the upper or lower triangular part of the matrix A from * data file * IF (UPLO.EQ.'U') THEN READ (NIN,*) ((AP(I+(J*(J-1))/2),J=I,N),I=1,N) ELSE IF (UPLO.EQ.'L') THEN READ (NIN,*) ((AP(I+((2*N-J)*(J-1))/2),J=1,I),I=1,N) END IF * * Print the matrix A * IFAIL = 0 CALL X04CCF(UPLO,'Non-unit diagonal',N,AP,'Matrix A',IFAIL) WRITE (NOUT,*) * * Compute diagonal scaling factors * CALL DPPEQU(UPLO,N,AP,S,SCOND,AMAX,INFO) * IF (INFO.GT.0) THEN WRITE (NOUT,99999) 'Diagonal element', INFO, + ' of A is non positive' ELSE * * Print SCOND, AMAX and the scale factors * WRITE (NOUT,99998) 'SCOND = ', SCOND, ', AMAX = ', AMAX WRITE (NOUT,*) WRITE (NOUT,*) 'Diagonal scaling factors' WRITE (NOUT,99997) (S(I),I=1,N) WRITE (NOUT,*) * * Compute values close to underflow and overflow * SMALL = X02AMF()/(X02AJF()*X02BHF()) BIG = 1.0D0/SMALL IF ((SCOND.LT.0.1D0) .OR. (AMAX.LT.SMALL) .OR. (AMAX.GT.BIG) + ) THEN * * Scale A * IF (UPLO.EQ.'U') THEN DO 40 J = 1, N SJ = S(J) JJ = (J*(J-1))/2 DO 20 I = 1, J AP(I+JJ) = S(I)*AP(I+JJ)*SJ 20 CONTINUE 40 CONTINUE ELSE IF (UPLO.EQ.'L') THEN DO 80 J = 1, N SJ = S(J) JJ = ((2*N-J)*(J-1))/2 DO 60 I = J, N AP(I+JJ) = S(I)*AP(I+JJ)*SJ 60 CONTINUE 80 CONTINUE END IF * * Print the scaled matrix * IFAIL = 0 CALL X04CCF(UPLO,'Non-unit diagonal',N,AP, + 'Scaled matrix',IFAIL) * END IF END IF ELSE WRITE (NOUT,*) 'NMAX too small' END IF * 99999 FORMAT (1X,A,I4,A) 99998 FORMAT (1X,2(A,1P,E7.1)) 99997 FORMAT ((1X,1P,7E11.1)) END