* F08ZPF Example Program Text * Mark 17 Release. NAG Copyright 1995. * .. Parameters .. INTEGER NIN, NOUT PARAMETER (NIN=5,NOUT=6) INTEGER MMAX, NB, NMAX, PMAX PARAMETER (MMAX=10,NB=64,NMAX=10,PMAX=10) INTEGER LDA, LDB, LWORK PARAMETER (LDA=NMAX,LDB=NMAX,LWORK=MMAX+NMAX+NB*(NMAX+PMAX) + ) * .. Local Scalars .. DOUBLE PRECISION RNORM INTEGER I, INFO, J, M, N, P * .. Local Arrays .. COMPLEX *16 A(LDA,MMAX), B(LDB,PMAX), D(NMAX), WORK(LWORK), + X(MMAX), Y(PMAX) * .. External Functions .. DOUBLE PRECISION DZNRM2 EXTERNAL DZNRM2 * .. External Subroutines .. EXTERNAL ZGGGLM * .. Executable Statements .. WRITE (NOUT,*) 'F08ZPF Example Program Results' WRITE (NOUT,*) * Skip heading in data file READ (NIN,*) READ (NIN,*) N, M, P IF (M.LE.MMAX .AND. N.LE.NMAX .AND. P.LE.PMAX) THEN * * Read A, B and D from data file * READ (NIN,*) ((A(I,J),J=1,M),I=1,N) READ (NIN,*) ((B(I,J),J=1,P),I=1,N) READ (NIN,*) (D(I),I=1,N) * * Solve the weighted least-squares problem * * minimize ||inv(B)*(d - A*x)|| (in the 2-norm) * CALL ZGGGLM(N,M,P,A,LDA,B,LDB,D,X,Y,WORK,LWORK,INFO) * * Print least-squares solution * WRITE (NOUT,*) 'Weighted least-squares solution' WRITE (NOUT,99999) (X(I),I=1,M) * * Print residual vector y = inv(B)*(d - A*x) * WRITE (NOUT,*) WRITE (NOUT,*) 'Residual vector' WRITE (NOUT,99998) (Y(I),I=1,P) * * Compute and print the square root of the residual sum of * squares * RNORM = DZNRM2(P,Y,1) * WRITE (NOUT,*) WRITE (NOUT,*) 'Square root of the residual sum of squares' WRITE (NOUT,99997) RNORM ELSE WRITE (NOUT,*) + 'One or more of MMAX, NMAX and PMAX is too small' END IF STOP * 99999 FORMAT (3(' (',F9.4,',',F9.4,')',:)) 99998 FORMAT (3(' (',1P,E9.2,',',1P,E9.2,')',:)) 99997 FORMAT (1X,1P,E10.2) END