C ALGORITHM 767, COLLECTED ALGORITHMS FROM ACM. C THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE, C VOL. 23, NO. 1, March, 1997, P. 111--129. C #! /bin/sh # This is a shell archive, meaning: # 1. Remove everything above the #! /bin/sh line. # 2. Save the resulting text in a file. # 3. Execute the file with /bin/sh (not csh) to create the files: # Drivers # Src # This archive created: Sat Apr 19 22:37:02 1997 export PATH; PATH=/bin:$PATH if test ! -d 'Drivers' then mkdir 'Drivers' fi cd 'Drivers' if test ! -d 'Dp' then mkdir 'Dp' fi cd 'Dp' if test -f 'data' then echo shar: will not over-write existing file "'data'" else cat << \SHAR_EOF > 'data' 11 COLRED EXAMPLE PROGRAM DATA (1) 1 1 1 0.0D-00 P0 1.0 P1 2.0 COLRED EXAMPLE PROGRAM DATA (2) 2 2 3 1.0D-10 P0 1.0D-00 1.0D-00 -2.0D-00 -1.0D-00 P1 2.0D-00 1.0D-02 -3.0D-00 0.0D-00 P2 1.0D-00 0.0D-00 2.0D-00 0.0D-00 P3 1.0D-00 0.0D-00 0.0D-00 0.0D-00 COLRED EXAMPLE PROGRAM DATA (3) 2 2 3 1.0D-12 P0 1.0D-00 1.0D-00 -2.0D-00 -1.0D-00 P1 2.0D-00 1.0D-02 -3.0D-00 0.0D-00 P2 1.0D-00 0.0D-00 2.0D-00 0.0D-00 P3 1.0D-00 0.0D-00 0.0D-00 0.0D-00 COLRED EXAMPLE PROGRAM DATA (4) 3 3 6 0.0D-00 P0 0.0D-00 0.0D-00 1.0D-00 0.0D-00 1.0D-00 0.0D-00 1.0D-00 0.0D-00 1.0D-00 P1 0.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 P2 0.0D-00 1.0D-00 0.0D-00 1.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 P3 0.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 P4 1.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 1.0D-00 0.0D-00 0.0D-00 0.0D-00 P5 0.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 P6 0.0D-00 0.0D-00 1.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 COLRED EXAMPLE PROGRAM DATA (5) 4 9 2 0.0D-00 P0 0.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 -1.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 -1.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 -1.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 -1.0D-00 P1 0.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 1.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 1.0D-00 1.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 1.0D-00 0.0D-00 0.0D-00 0.0D-00 1.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 1.0D-00 0.0D-00 0.0D-00 0.0D-00 1.0D-00 P2 1.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 COLRED EXAMPLE PROGRAM DATA (6) 2 2 1 0.0D-00 P0 1.0D-00 1.0D-00 -2.0D-00 -2.0D-00 P1 1.0D-00 1.0D-00 0.0D-00 0.0D-00 COLRED EXAMPLE PROGRAM DATA (7) 3 2 0 0.0D-00 P0 0.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 COLRED EXAMPLE PROGRAM DATA (8) 1 2 0 0.0D-00 P0 1.5D-00 2.7D-00 COLRED EXAMPLE PROGRAM DATA (9) 1 3 0 0.0D-00 P0 1.2D-00 1.56D-00 -0.361D-00 COLRED EXAMPLE PROGRAM DATA (10) 4 4 0 0.0D-00 P0 0.0D-00 0.0D-00 1.0D-00 1.0D-00 1.0D-00 1.0D-00 2.0D-00 2.0D-00 0.0D-00 0.0D-00 0.0D-00 0.0D-00 1.0D-00 1.0D-00 -3.0D-00 -3.0D-00 COLRED EXAMPLE PROGRAM DATA (11) 0 0 0 0.0D-00 P0 1.0D-00 2.0D-00 3 3 6 0.0D-00 SHAR_EOF fi # end of overwriting check if test -f 'driver.f' then echo shar: will not over-write existing file "'driver.f'" else cat << \SHAR_EOF > 'driver.f' C C COLRED TEST DRIVER TEXT C C Version dd February 24, 1995. C Revised May 21, 1996. C C IMPLICIT NONE C .. Parameters .. INTEGER NIN, NOUT PARAMETER (NIN = 5, NOUT = 6) DOUBLE PRECISION ONE PARAMETER (ONE = 1.0D0) C INTEGER DPMAX, MPMAX, NPMAX PARAMETER (DPMAX = 6, MPMAX = 8, NPMAX = 10) INTEGER MAMAX, NAMAX, LIWORK, LRWORK, DRMAX PARAMETER (MAMAX = (NPMAX * DPMAX + 1) * MPMAX, * NAMAX = MAMAX + NPMAX, * LIWORK = 2 * MAMAX + 2 * NPMAX + 1, * LRWORK = (2*MAMAX + NPMAX*DPMAX + 5) * NAMAX + * MAMAX + (2*MPMAX + 1) * NPMAX, * DRMAX = (NPMAX+1)*DPMAX + 1) INTEGER LDP1, LDP2, LDR1, LDR2, LDU1, LDU2 PARAMETER (LDP1 = MPMAX, LDP2 = NPMAX, LDR1 = MPMAX, * LDR2 = NPMAX, LDU1 = NPMAX, LDU2 = NPMAX) C .. Local Scalars .. INTEGER MP, NP, DP, DR, DU, I, J, K, L, NUM, IERR DOUBLE PRECISION TOL C .. Local Arrays .. INTEGER IWORK(LIWORK) DOUBLE PRECISION P(LDP1,LDP2,DPMAX+1), R(LDR1,LDR2,DRMAX), * U(LDU1,LDU2,NPMAX*DPMAX+1), RWORK(LRWORK) LOGICAL ZERCOL(NPMAX) C .. External Subroutines/Functions .. EXTERNAL COLRED, MULTPM, PRMAPO C C .. Executable Statements .. C C Read NUM, that is the number of testcases to be run. C READ(NIN, FMT = *) NUM C DO 30 L = 1, NUM C WRITE(NOUT, FMT = 99999) L READ(NIN, FMT = '()') C C Read MP, NP, DP, TOL and next P(k), k = 0,....,DP row after row. C READ(NIN, FMT = *) MP, NP, DP, TOL C DO 20 K = 1, DP + 1 READ(NIN, FMT = '()') DO 10 I = 1, MP READ(NIN, FMT = *) (P(I,J,K), J = 1, NP) 10 CONTINUE 20 CONTINUE C WRITE(NOUT, FMT = 99998) DP, MP, NP, TOL CALL PRMAPO(MP, NP, DP, 5, NOUT, P, LDP1, LDP2, 'P', IERR) C CALL COLRED(MP, NP, DP, P, LDP1, LDP2, DR, DU, R, LDR1, LDR2, * U, LDU1, LDU2, ZERCOL, IWORK, RWORK, TOL, IERR) C IF (IERR .EQ. 0) THEN WRITE (NOUT, FMT = 99997) CALL PRMAPO(NP, NP, DU, 5, NOUT, U, LDU1, LDU2, 'U', IERR) C WRITE (NOUT, FMT = 99996) CALL PRMAPO(MP, NP, DR, 5, NOUT, R, LDR1, LDR2, 'R', IERR) WRITE (NOUT, FMT = 99995) (ZERCOL(J), J = 1, NP) C CALL MULTPM(-ONE, MP, NP, NP, DP, DU, DR, P, LDP1, LDP2, * U, LDU1, LDU2, R, LDR1, LDR2, RWORK, IERR) IF (DR .GE. 0) THEN WRITE (NOUT, FMT = 99994) CALL PRMAPO(MP, NP, DR, 5, NOUT, R, LDR1, LDR2, * '(PU - R)', IERR) ELSE WRITE (NOUT, FMT = 99993) END IF ELSE WRITE (NOUT, FMT = 99992) IERR END IF 30 CONTINUE STOP C 99999 FORMAT (//' COLRED TEST DRIVER RESULTS, EXAMPLE(', I2, ')', /1X) 99998 FORMAT (' The input polynomial matrix:', //, * ' P(s) = P(0) + P(1) * s + . . . + P(dp-1) * s**(dp-1)', * ' + P(dp) * s**dp', //, ' with degree DP =', I2, * ' and size MP =', I2, ', NP =', I2, '.', //, * ' The tolerance is:', D10.3, /1X) 99997 FORMAT (' The unimodular polynomial matrix U(s):') 99996 FORMAT (' The column reduced polynomial matrix R(s):') 99995 FORMAT (' ZERCOL(j), j = 1, NP:', 10(L2)) 99994 FORMAT (' The residual matrix P(s) * U(s) - R(s):') 99993 FORMAT (' PU - R is the ZERO polynomial matrix.') 99992 FORMAT (' COLRED has failed: IERR =', I2) END C SUBROUTINE MULTPM(ALPHA, RP1, CP1, CP2, DP1, DP2, DP3, P1, LDP11, * LDP12, P2, LDP21, LDP22, P3, LDP31, LDP32, * RWORK, IERR) C C PURPOSE C C To compute the coefficients of the real polynomial matrix C C P(s) = P1(s) * P2(s) + alpha * P3(s), (1) C C where P1(s), P2(s) and P3(s) are real polynomial matrices and C alpha is a real scalar. C Each of the polynomial matrices may by the zero matrix. C C ARGUMENTS IN C C ALPHA - DOUBLE PRECISION. C The value of the scalar factor alpha of the problem. C RP1 - INTEGER. C The number of rows of the matrices P1(s) and P3(s). C RP1 >= 1. C CP1 - INTEGER. C The number of columns of matrix P1(s) and the number of rows C of matrix P2(s). C CP1 >= 1. C CP2 - INTEGER. C The number of columns of the matrices P2(s) and P3(s). C CP2 >= 1. C DP1 - INTEGER. C The degree of the polynomial matrix P1(s). C DP1 >= -1. C DP2 - INTEGER. C The degree of the polynomial matrix P2(s). C DP2 >= -1. C DP3 - INTEGER. C The degree of the polynomial matrix P3(s). C DP3 >= -1. C Note: DP3 is overwritten. C P1 - DOUBLE PRECISION array of (LDP11,LDP12,*). C If DP1 >= 0, then the leading RP1 by CP1 by (DP1+1) part of C this array must contain the coefficients of the polynomial C matrix P1(s). Specifically, P1(i,j,k) must contain the C coefficient of s**(k-1) of the polynomial which is the C (i,j)-th element of P1(s), where i = 1,2, ...,RP1, C j = 1,2,...,CP1 and k = 1,2,...,DP1+1. C If DP1 = -1, then P1(s) is taken to be the zero polynomial C matrix and the array P1 is not referenced. C LDP11 - INTEGER. C The leading dimension of array P1 as declared in the calling C program. C LDP11 >= RP1 if DP1 >= 0. C LDP11 >= 1 if DP1 = -1. C LDP12 - INTEGER. C The second dimension of array P1 as declared in the calling C program. C LDP12 >= CP1 if DP1 >= 0, C LDP12 >= 1 if DP1 = -1. C P2 - DOUBLE PRECISION array of (LDP21,LDP22,*). C If DP2 >= 0, then the leading CP1 by CP2 by (DP2+1) part of C this array must contain the coefficients of the polynomial C matrix P2(s). Specifically, P2(i,j,k) must contain the C coefficient of s**(k-1) of the polynomial which is the C (i,j)-th element of P2(s), where i = 1,2, ...,CP1, C j = 1,2,...,CP2 and k = 1,2,...,DP2+1. C If DP2 = -1, then P2(s) is taken to be the zero polynomial C matrix and the array P2 is not referenced. C LDP21 - INTEGER. C The leading dimension of array P2 as declared in the calling C program. C LDP21 >= CP1 if DP2 >= 0. C LDP21 >= 1 if DP2 = -1. C LDP22 - INTEGER. C The second dimension of array P2 as declared in the calling C program. C LDP22 >= CP2 if DP2 >= 0, C LDP22 >= 1 if DP2 = -1. C P3 - DOUBLE PRECISION array of (LDP31,LDP32,lenp3), C where lenp3 = MAX(DP1+DP2,DP3,0) + 1. C If DP3 >= 0, then the leading RP1 by CP2 by (DP3+1) part of C this array must contain the coefficients of the polynomial C matrix P3(s). Specifically, P3(i,j,k) must contain the C coefficient of s**(k-1) of the polynomial which is the C (i,j)-th element of P3(s), where i = 1,2, ...,RP1, C j = 1,2,...,CP2 and k = 1,2,...,DP3+1. C If DP3 = -1, then P3(s) is taken to be the zero polynomial C matrix. C Note: this array is overwritten. C LDP31 - INTEGER. C The leading dimension of array P3 as declared in the calling C program. C LDP31 >= RP1. C LDP32 - INTEGER. C The second dimension of array P3 as declared in the calling C program. C LDP32 >= CP2. C C ARGUMENTS OUT C C DP3 - INTEGER. C The degree of the resulting polynomial matrix P(s). C P3 - DOUBLE PRECISION array of DIMENSION (LDP31,LDP32,lenp3). C If DP3 >= 0 on exit, then the leading RP1 by CP2 by (DP3+1) C part of this array contains the coefficients of P(s). C Specifically, P3(i,j,k) contains the coefficient of s**(k-1) C of the polynomial which is the (i,j)-th element of P(s), C where i = 1,2, ...,RP1, j = 1,2,...,CP2 and C k = 1,2,...,DP3+1. C If DP3 = -1 on exit, then P(s) is the zero polynomial C matrix and the contents of the array P3 are undefined. C C WORK SPACE C C RWORK - DOUBLE PRECISION array of DIMENSION at least (CP1). C C ERROR INDICATOR C C IERR - INTEGER. C Unless the routine detects an error (see next section), C IERR contains 0 on exit. C C WARNINGS AND ERRORS DETECTED BY THE ROUTINE C C IERR = 1 : Invalid input parameter(s). C C METHOD C C Given the real polynomial matrices C DP1 i DP2 i C P1(s) = SUM (a(i+1) * s ), P2(s) = SUM (b(i+1) * s ), C i=0 i=0 C DP3 i C P3(s) = SUM (c(i+1) * s ). C i=0 C and a real scalar alpha, the routine computes the coefficients C d(1), d(2), ... of the polynomial matrix (1) from the formula C s C d(i+1) := SUM (a(k+1) * b(i-k+1)) + alpha * c(i+1), C k=r C where i = 0,1,...,DP1+DP2 and r and s depend on the value of i, C i.e. for r <= k <= s both a(k+1) and b(i-k+1) must exist. C C CONTRIBUTOR C C A.J. Geurts (Eindhoven University of Technology). C C REVISIONS C C 1992, October 28. C 1994, December 8. C C IMPLICIT NONE C .. Parameters .. DOUBLE PRECISION ZERO PARAMETER (ZERO = 0.0D0) C .. Scalar Arguments .. INTEGER RP1, CP1, CP2, DP1, DP2, DP3, LDP11, LDP12, * LDP21, LDP22, LDP31, LDP32, IERR DOUBLE PRECISION ALPHA C .. Array Arguments .. DOUBLE PRECISION P1(LDP11, LDP12, *), P2(LDP21, LDP22, *), * P3(LDP31, LDP32, *), RWORK(*) C .. Local Scalars .. INTEGER H, I, J, K, DPOL3, E DOUBLE PRECISION W LOGICAL CFZERO C .. External Subroutines/Functions .. EXTERNAL DCOPY, DDOT, DLAVEC, DSCAL DOUBLE PRECISION DDOT C C .. Executable Statements .. C C Check input parameters. C IF ((RP1.LT.1) .OR. (CP1.LT.1) .OR. (CP2.LT.1) * .OR. (DP1.LT.-1) .OR. (DP2.LT.-1) .OR. (DP3.LT.-1) * .OR. ((LDP11.LT.RP1) .AND. (DP1.GE.0)) * .OR. ((LDP11.LT.1) .AND. (DP1.EQ.-1)) * .OR. ((LDP12.LT.CP1) .AND. (DP1.GE.0)) * .OR. ((LDP12.LT.1) .AND. (DP1.EQ.-1)) * .OR. ((LDP21.LT.CP1) .AND. (DP2.GE.0)) * .OR. ((LDP21.LT.1) .AND. (DP2.EQ.-1)) * .OR. ((LDP22.LT.CP2) .AND. (DP2.GE.0)) * .OR. ((LDP22.LT.1) .AND. (DP2.EQ.-1)) * .OR. (LDP31.LT.RP1) .OR. (LDP32.LT.CP2)) THEN IERR = 1 RETURN END IF C IERR = 0 IF (ALPHA .EQ. ZERO) THEN DP3 = -1 END IF C IF (DP3 .GE. 0) THEN C C P3(s) := ALPHA * P3(s). C DO 20 K = 1, DP3 + 1 DO 10 J = 1, CP2 CALL DSCAL(RP1, ALPHA, P3(1,J,K), 1) 10 CONTINUE 20 CONTINUE END IF C IF ((DP1 .EQ. -1) .OR. (DP2 .EQ. -1)) RETURN C C Neither of P1(s) and P2(s) is the zero polynomial. C DPOL3 = DP1 + DP2 IF (DPOL3 .GT. DP3) THEN C C Initialize the additional part of P3(s) to zero. C DO 40 K = DP3 + 2, DPOL3 + 1 DO 30 J = 1, CP2 CALL DLAVEC(RP1, ZERO, P3(1,J,K), 1) 30 CONTINUE 40 CONTINUE DP3 = DPOL3 END IF C k-1 C The inner product of the j-th row of the coefficient of s of C i-1 C P1(s) and the h-th column of the coefficient of s of P2(s) C k+i-2 C contributes to the (j,h)-th element of the coefficient of s C of P3(s). C DO 80 K = 1, DP1 + 1 DO 70 J = 1, RP1 CALL DCOPY(CP1, P1(J,1,K), LDP11, RWORK, 1) DO 60 I = 1, DP2 + 1 E = K + I - 1 DO 50 H = 1, CP2 W = DDOT(CP1, RWORK, 1, P2(1,H,I), 1) P3(J,H,E) = W + P3(J,H,E) 50 CONTINUE 60 CONTINUE 70 CONTINUE 80 CONTINUE C C Computation of the exact degree of P3(s). C CFZERO = .TRUE. C WHILE (DP3 >= 0 and CFZERO) DO 90 IF ((DP3 .GE. 0) .AND. CFZERO) THEN DPOL3 = DP3 + 1 DO 110 I = 1, RP1 DO 100 J = 1, CP2 IF (P3(I,J,DPOL3) .NE. ZERO) CFZERO = .FALSE. 100 CONTINUE 110 CONTINUE IF (CFZERO) DP3 = DP3 - 1 GO TO 90 END IF C END WHILE 90 C RETURN C *** Last line of MULTPM *** END C SUBROUTINE PRMAPO(MP, NP, DP, L, NOUT, P, LDP1, LDP2, TEXT, IERR) C C PURPOSE C C To print the MP by NP coefficient matrices of a matrix polynomial C dp-1 dp C P(s) = P(0) + P(1) * s + . . . P(dp-1) * s + P(dp) * s . C C The elements of the matrices are output to 7 significant figures. C C ARGUMENTS IN C C MP - INTEGER. C The number of rows of the matrix polynomial P(s). C MP >= 1. C NP - INTEGER. C The number of columns of the matrix polynomial P(s). C NP >= 1. C DP - INTEGER. C The degree of the matrix polynomial P(s). C DP >= 0. C L - INTEGER. C The number of elements of the coefficient matrices to be C printed per line. C 1 <= L <= 5. C NOUT - INTEGER. C The output channel to which the results are sent. C NOUT >= 0. C P - DOUBLE PRECISION array of DIMENSION (LDP1,LDP2,DP+1). C The leading MP by NP by (DP+1) part of this array must C contain the coefficients of the matrix polynomial P(s). C Specifically, P(i,j,k) must contain the coefficient of C s**(k-1) of the polynomial which is the (i,j)-th element C of P(s), where i = 1,2,...,MP, j = 1,2,...,NP and C k = 1,2,...,DP+1. C LDP1 - INTEGER. C The leading dimension of array P as declared in the calling C program. C LDP1 >= MP. C LDP2 - INTEGER. C The second dimension of array P as declared in the calling C program. C LDP2 >= NP. C TEXT - CHARACTER*72. C Title caption of the coefficient matrices to be printed. C TEXT is followed by the degree of the coefficient matrix, C within brackets. If TEXT = ' ', then the coefficient C matrices are separated by an empty line. C C ERROR INDICATOR C C IERR - INTEGER. C Unless the routine detects an error (see next section), C IERR contains 0 on exit. C C WARNINGS AND ERRORS DETECTED BY THE ROUTINE C C IERR = 1 : Invalid input parameter(s). C C METHOD C C For i = 1, 2, ..., DP + 1 the routine first prints the contents of C TEXT followed by (i-1) as a title, followed by the elements of the C MP by NP coefficient matrix P(i) such that C (i) if NP < L, then the leading MP by NP part is printed; C (ii) if NP = k*L + p (where k, p > 0), then k MP by L blocks of C consecutive columns of P(i) are printed one after another C followed by one MP by p block containing the last p columns of C P(i). C Row numbers are printed on the left of each row and a column C number on top of each column. C C CONTRIBUTOR C C A.J. Geurts (Eindhoven University of Technology). C C REVISIONS C C 1992, October 28. C 1995, February 6. C C IMPLICIT NONE C .. Scalar Arguments .. INTEGER MP, NP, DP, L, NOUT, LDP1, LDP2, IERR CHARACTER*(*) TEXT C .. Array Arguments .. DOUBLE PRECISION P(LDP1,LDP2,*) C .. Local Scalars .. INTEGER I, J, J1, J2, JJ, K, LENTXT, LTEXT, N1 C .. Intrinsic Functions .. INTRINSIC LEN, MIN C C .. Executable Statements .. C C Check input parameters. C LENTXT = LEN(TEXT) LTEXT = MIN(72,LENTXT) C WHILE (TEXT(LTEXT:LTEXT) = ' ') DO 10 IF (TEXT(LTEXT:LTEXT) .EQ. ' ') THEN LTEXT = LTEXT - 1 GO TO 10 END IF C END WHILE 10 C IF (MP.LT.1 .OR. NP.LT.1 .OR. DP.LT.0 .OR. L.LT.1 .OR. L.GT.5 * .OR. NOUT.LT.0 .OR. LDP1.LT.MP .OR. LDP2.LT.NP) THEN IERR = 1 RETURN END IF C IERR = 0 C DO 50 K = 1, DP + 1 IF (LTEXT .EQ. 0) THEN WRITE (NOUT, FMT = 99999) ELSE WRITE (NOUT, FMT = 99998) TEXT(1:LTEXT), K - 1 END IF N1 = (NP - 1)/L J1 = 1 J2 = L DO 30 J = 1, N1 WRITE (NOUT, FMT = 99997) (JJ, JJ = J1, J2) DO 20 I = 1, MP WRITE (NOUT, FMT = 99996) I, (P(I,JJ,K), JJ = J1, J2) 20 CONTINUE J1 = J1 + L J2 = J2 + L 30 CONTINUE WRITE (NOUT, FMT = 99997) (J, J = J1, NP) DO 40 I = 1, MP WRITE (NOUT, FMT = 99996) I, (P(I,JJ,K), JJ = J1, NP) 40 CONTINUE 50 CONTINUE WRITE (NOUT, FMT = 99999) C RETURN 99999 FORMAT (' ') 99998 FORMAT (/, 1X, A, '(', I2, ')') 99997 FORMAT (5X, 5(6X, I2, 7X)) 99996 FORMAT (1X, I2, 2X, 5F15.7) C *** Last line of PRMAPO *** END SHAR_EOF fi # end of overwriting check if test -f 'RES' then echo shar: will not over-write existing file "'RES'" else cat << \SHAR_EOF > 'RES' COLRED TEST DRIVER RESULTS, EXAMPLE( 1) The input polynomial matrix: P(s) = P(0) + P(1) * s + . . . + P(dp-1) * s**(dp-1) + P(dp) * s**dp with degree DP = 1 and size MP = 1, NP = 1. The tolerance is: 0.000D+00 P( 0) 1 1 1.0000000 P( 1) 1 1 2.0000000 The unimodular polynomial matrix U(s): U( 0) 1 1 1.0000000 The column reduced polynomial matrix R(s): R( 0) 1 1 1.0000000 R( 1) 1 1 2.0000000 ZERCOL(j), j = 1, NP: F PU - R is the ZERO polynomial matrix. COLRED TEST DRIVER RESULTS, EXAMPLE( 2) The input polynomial matrix: P(s) = P(0) + P(1) * s + . . . + P(dp-1) * s**(dp-1) + P(dp) * s**dp with degree DP = 3 and size MP = 2, NP = 2. The tolerance is: 0.100D-09 P( 0) 1 2 1 1.0000000 1.0000000 2 -2.0000000 -1.0000000 P( 1) 1 2 1 2.0000000 0.0100000 2 -3.0000000 0.0000000 P( 2) 1 2 1 1.0000000 0.0000000 2 2.0000000 0.0000000 P( 3) 1 2 1 1.0000000 0.0000000 2 0.0000000 0.0000000 The unimodular polynomial matrix U(s): U( 0) 1 2 1 0.0000000 1.4066253 2 151.9532407 0.0000000 U( 1) 1 2 1 0.0000000 0.0000000 2 0.0000000 -475.9339805 U( 2) 1 2 1 0.0000000 0.0000000 2 0.0000000 -140.6625313 The column reduced polynomial matrix R(s): R( 0) 1 2 1 151.9532407 1.4066253 2 -151.9532407 -2.8132506 R( 1) 1 2 1 1.5195324 -473.1207299 2 0.0000000 471.7141046 R( 2) 1 2 1 0.0000000 -144.0152458 2 0.0000000 143.4757819 ZERCOL(j), j = 1, NP: F F The residual matrix P(s) * U(s) - R(s): (PU - R)( 0) 1 2 1 0.0000000 0.0000000 2 0.0000000 0.0000000 (PU - R)( 1) 1 2 1 0.0000000 0.0000000 2 0.0000000 0.0000000 (PU - R)( 2) 1 2 1 0.0000000 0.0000000 2 0.0000000 0.0000000 (PU - R)( 3) 1 2 1 0.0000000 0.0000000 2 0.0000000 0.0000000 (PU - R)( 4) 1 2 1 0.0000000 0.0000000 2 0.0000000 0.0000000 (PU - R)( 5) 1 2 1 0.0000000 0.0000000 2 0.0000000 0.0000000 COLRED TEST DRIVER RESULTS, EXAMPLE( 3) The input polynomial matrix: P(s) = P(0) + P(1) * s + . . . + P(dp-1) * s**(dp-1) + P(dp) * s**dp with degree DP = 3 and size MP = 2, NP = 2. The tolerance is: 0.100D-11 P( 0) 1 2 1 1.0000000 1.0000000 2 -2.0000000 -1.0000000 P( 1) 1 2 1 2.0000000 0.0100000 2 -3.0000000 0.0000000 P( 2) 1 2 1 1.0000000 0.0000000 2 2.0000000 0.0000000 P( 3) 1 2 1 1.0000000 0.0000000 2 0.0000000 0.0000000 The unimodular polynomial matrix U(s): U( 0) 1 2 1 0.0000000 1.4066253 2 151.9532407 0.0000000 U( 1) 1 2 1 0.0000000 0.0000000 2 0.0000000 -475.9339805 U( 2) 1 2 1 0.0000000 0.0000000 2 0.0000000 -140.6625313 The column reduced polynomial matrix R(s): R( 0) 1 2 1 151.9532407 1.4066253 2 -151.9532407 -2.8132506 R( 1) 1 2 1 1.5195324 -473.1207299 2 0.0000000 471.7141046 R( 2) 1 2 1 0.0000000 -144.0152458 2 0.0000000 143.4757819 ZERCOL(j), j = 1, NP: F F The residual matrix P(s) * U(s) - R(s): (PU - R)( 0) 1 2 1 0.0000000 0.0000000 2 0.0000000 0.0000000 (PU - R)( 1) 1 2 1 0.0000000 0.0000000 2 0.0000000 0.0000000 (PU - R)( 2) 1 2 1 0.0000000 0.0000000 2 0.0000000 0.0000000 (PU - R)( 3) 1 2 1 0.0000000 0.0000000 2 0.0000000 0.0000000 (PU - R)( 4) 1 2 1 0.0000000 0.0000000 2 0.0000000 0.0000000 (PU - R)( 5) 1 2 1 0.0000000 0.0000000 2 0.0000000 0.0000000 COLRED TEST DRIVER RESULTS, EXAMPLE( 4) The input polynomial matrix: P(s) = P(0) + P(1) * s + . . . + P(dp-1) * s**(dp-1) + P(dp) * s**dp with degree DP = 6 and size MP = 3, NP = 3. The tolerance is: 0.000D+00 P( 0) 1 2 3 1 0.0000000 0.0000000 1.0000000 2 0.0000000 1.0000000 0.0000000 3 1.0000000 0.0000000 1.0000000 P( 1) 1 2 3 1 0.0000000 0.0000000 0.0000000 2 0.0000000 0.0000000 0.0000000 3 0.0000000 0.0000000 0.0000000 P( 2) 1 2 3 1 0.0000000 1.0000000 0.0000000 2 1.0000000 0.0000000 0.0000000 3 0.0000000 0.0000000 0.0000000 P( 3) 1 2 3 1 0.0000000 0.0000000 0.0000000 2 0.0000000 0.0000000 0.0000000 3 0.0000000 0.0000000 0.0000000 P( 4) 1 2 3 1 1.0000000 0.0000000 0.0000000 2 0.0000000 0.0000000 1.0000000 3 0.0000000 0.0000000 0.0000000 P( 5) 1 2 3 1 0.0000000 0.0000000 0.0000000 2 0.0000000 0.0000000 0.0000000 3 0.0000000 0.0000000 0.0000000 P( 6) 1 2 3 1 0.0000000 0.0000000 1.0000000 2 0.0000000 0.0000000 0.0000000 3 0.0000000 0.0000000 0.0000000 The unimodular polynomial matrix U(s): U( 0) 1 2 3 1 1.0000000 -1.0000000 0.0000000 2 0.0000000 0.0000000 -1.7029386 3 0.0000000 1.0000000 0.0000000 U( 1) 1 2 3 1 0.0000000 0.0000000 0.0000000 2 0.0000000 0.0000000 0.0000000 3 0.0000000 0.0000000 0.0000000 U( 2) 1 2 3 1 0.0000000 0.0000000 -1.7029386 2 -1.0000000 1.0000000 0.0000000 3 0.0000000 0.0000000 1.7029386 U( 3) 1 2 3 1 0.0000000 0.0000000 0.0000000 2 0.0000000 0.0000000 0.0000000 3 0.0000000 0.0000000 0.0000000 U( 4) 1 2 3 1 0.0000000 0.0000000 0.0000000 2 0.0000000 -1.0000000 1.7029386 3 0.0000000 0.0000000 0.0000000 U( 5) 1 2 3 1 0.0000000 0.0000000 0.0000000 2 0.0000000 0.0000000 0.0000000 3 0.0000000 0.0000000 0.0000000 U( 6) 1 2 3 1 0.0000000 0.0000000 0.0000000 2 0.0000000 0.0000000 -1.7029386 3 0.0000000 0.0000000 0.0000000 The column reduced polynomial matrix R(s): R( 0) 1 2 3 1 0.0000000 1.0000000 0.0000000 2 0.0000000 0.0000000 -1.7029386 3 1.0000000 0.0000000 0.0000000 ZERCOL(j), j = 1, NP: F F F The residual matrix P(s) * U(s) - R(s): (PU - R)( 0) 1 2 3 1 0.0000000 0.0000000 0.0000000 2 0.0000000 0.0000000 0.0000000 3 0.0000000 0.0000000 0.0000000 (PU - R)( 1) 1 2 3 1 0.0000000 0.0000000 0.0000000 2 0.0000000 0.0000000 0.0000000 3 0.0000000 0.0000000 0.0000000 (PU - R)( 2) 1 2 3 1 0.0000000 0.0000000 0.0000000 2 0.0000000 0.0000000 0.0000000 3 0.0000000 0.0000000 0.0000000 (PU - R)( 3) 1 2 3 1 0.0000000 0.0000000 0.0000000 2 0.0000000 0.0000000 0.0000000 3 0.0000000 0.0000000 0.0000000 (PU - R)( 4) 1 2 3 1 0.0000000 0.0000000 0.0000000 2 0.0000000 0.0000000 0.0000000 3 0.0000000 0.0000000 0.0000000 (PU - R)( 5) 1 2 3 1 0.0000000 0.0000000 0.0000000 2 0.0000000 0.0000000 0.0000000 3 0.0000000 0.0000000 0.0000000 (PU - R)( 6) 1 2 3 1 0.0000000 0.0000000 0.0000000 2 0.0000000 0.0000000 0.0000000 3 0.0000000 0.0000000 0.0000000 (PU - R)( 7) 1 2 3 1 0.0000000 0.0000000 0.0000000 2 0.0000000 0.0000000 0.0000000 3 0.0000000 0.0000000 0.0000000 (PU - R)( 8) 1 2 3 1 0.0000000 0.0000000 0.0000000 2 0.0000000 0.0000000 0.0000000 3 0.0000000 0.0000000 0.0000000 (PU - R)( 9) 1 2 3 1 0.0000000 0.0000000 0.0000000 2 0.0000000 0.0000000 0.0000000 3 0.0000000 0.0000000 0.0000000 (PU - R)(10) 1 2 3 1 0.0000000 0.0000000 0.0000000 2 0.0000000 0.0000000 0.0000000 3 0.0000000 0.0000000 0.0000000 (PU - R)(11) 1 2 3 1 0.0000000 0.0000000 0.0000000 2 0.0000000 0.0000000 0.0000000 3 0.0000000 0.0000000 0.0000000 (PU - R)(12) 1 2 3 1 0.0000000 0.0000000 0.0000000 2 0.0000000 0.0000000 0.0000000 3 0.0000000 0.0000000 0.0000000 COLRED TEST DRIVER RESULTS, EXAMPLE( 5) The input polynomial matrix: P(s) = P(0) + P(1) * s + . . . + P(dp-1) * s**(dp-1) + P(dp) * s**dp with degree DP = 2 and size MP = 4, NP = 9. The tolerance is: 0.000D+00 P( 0) 1 2 3 4 5 1 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 2 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 3 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 4 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 6 7 8 9 1 -1.0000000 0.0000000 0.0000000 0.0000000 2 0.0000000 -1.0000000 0.0000000 0.0000000 3 0.0000000 0.0000000 -1.0000000 0.0000000 4 0.0000000 0.0000000 0.0000000 -1.0000000 P( 1) 1 2 3 4 5 1 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 2 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 3 0.0000000 0.0000000 0.0000000 1.0000000 0.0000000 4 0.0000000 0.0000000 0.0000000 0.0000000 1.0000000 6 7 8 9 1 1.0000000 0.0000000 0.0000000 0.0000000 2 0.0000000 1.0000000 1.0000000 0.0000000 3 0.0000000 0.0000000 1.0000000 0.0000000 4 0.0000000 0.0000000 0.0000000 1.0000000 P( 2) 1 2 3 4 5 1 1.0000000 0.0000000 0.0000000 0.0000000 0.0000000 2 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 3 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 4 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 6 7 8 9 1 0.0000000 0.0000000 0.0000000 0.0000000 2 0.0000000 0.0000000 0.0000000 0.0000000 3 0.0000000 0.0000000 0.0000000 0.0000000 4 0.0000000 0.0000000 0.0000000 0.0000000 The unimodular polynomial matrix U(s): U( 0) 1 2 3 4 5 1 0.0000000 0.0000000 0.0000000 1.4142136 0.0000000 2 -1.0000000 0.0000000 0.0000000 0.0000000 0.0000000 3 0.0000000 -1.0000000 0.0000000 0.0000000 0.0000000 4 0.0000000 0.0000000 0.0000000 0.0000000 1.7320508 5 0.0000000 0.0000000 1.0000000 0.0000000 0.0000000 6 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 7 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 8 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 9 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 6 7 8 9 1 1.0000000 0.0000000 0.0000000 0.0000000 2 0.0000000 0.0000000 0.0000000 0.0000000 3 0.0000000 0.0000000 0.0000000 0.0000000 4 0.0000000 -1.0000000 2.0000000 0.0000000 5 0.0000000 0.0000000 0.0000000 1.0000000 6 -1.0000000 0.0000000 0.0000000 0.0000000 7 0.0000000 1.0000000 0.0000000 0.0000000 8 0.0000000 0.0000000 -1.0000000 0.0000000 9 0.0000000 0.0000000 0.0000000 -1.0000000 U( 1) 1 2 3 4 5 1 0.0000000 0.0000000 0.0000000 -1.4142136 0.0000000 2 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 3 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 4 0.0000000 0.0000000 0.0000000 0.0000000 -3.4641016 5 0.0000000 0.0000000 -1.0000000 0.0000000 0.0000000 6 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 7 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 8 0.0000000 0.0000000 0.0000000 0.0000000 1.7320508 9 0.0000000 0.0000000 1.0000000 0.0000000 0.0000000 6 7 8 9 1 0.0000000 0.0000000 0.0000000 0.0000000 2 0.0000000 0.0000000 0.0000000 0.0000000 3 0.0000000 0.0000000 0.0000000 0.0000000 4 0.0000000 1.0000000 -1.0000000 0.0000000 5 0.0000000 0.0000000 0.0000000 0.0000000 6 -1.0000000 0.0000000 0.0000000 0.0000000 7 0.0000000 1.0000000 -1.0000000 0.0000000 8 0.0000000 -1.0000000 1.0000000 0.0000000 9 0.0000000 0.0000000 0.0000000 0.0000000 U( 2) 1 2 3 4 5 1 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 2 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 3 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 4 0.0000000 0.0000000 0.0000000 0.0000000 1.7320508 5 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 6 0.0000000 0.0000000 0.0000000 1.4142136 0.0000000 7 0.0000000 0.0000000 0.0000000 0.0000000 1.7320508 8 0.0000000 0.0000000 0.0000000 0.0000000 -1.7320508 9 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 6 7 8 9 1 0.0000000 0.0000000 0.0000000 0.0000000 2 0.0000000 0.0000000 0.0000000 0.0000000 3 0.0000000 0.0000000 0.0000000 0.0000000 4 0.0000000 0.0000000 0.0000000 0.0000000 5 0.0000000 0.0000000 0.0000000 0.0000000 6 0.0000000 0.0000000 0.0000000 0.0000000 7 0.0000000 0.0000000 0.0000000 0.0000000 8 0.0000000 0.0000000 0.0000000 0.0000000 9 0.0000000 0.0000000 0.0000000 0.0000000 The column reduced polynomial matrix R(s): R( 0) 1 2 3 4 5 1 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 2 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 3 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 4 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 6 7 8 9 1 1.0000000 0.0000000 0.0000000 0.0000000 2 0.0000000 -1.0000000 0.0000000 0.0000000 3 0.0000000 0.0000000 1.0000000 0.0000000 4 0.0000000 0.0000000 0.0000000 1.0000000 ZERCOL(j), j = 1, NP: T T T T T F F F F The residual matrix P(s) * U(s) - R(s): (PU - R)( 0) 1 2 3 4 5 1 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 2 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 3 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 4 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 6 7 8 9 1 0.0000000 0.0000000 0.0000000 0.0000000 2 0.0000000 0.0000000 0.0000000 0.0000000 3 0.0000000 0.0000000 0.0000000 0.0000000 4 0.0000000 0.0000000 0.0000000 0.0000000 (PU - R)( 1) 1 2 3 4 5 1 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 2 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 3 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 4 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 6 7 8 9 1 0.0000000 0.0000000 0.0000000 0.0000000 2 0.0000000 0.0000000 0.0000000 0.0000000 3 0.0000000 0.0000000 0.0000000 0.0000000 4 0.0000000 0.0000000 0.0000000 0.0000000 (PU - R)( 2) 1 2 3 4 5 1 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 2 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 3 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 4 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 6 7 8 9 1 0.0000000 0.0000000 0.0000000 0.0000000 2 0.0000000 0.0000000 0.0000000 0.0000000 3 0.0000000 0.0000000 0.0000000 0.0000000 4 0.0000000 0.0000000 0.0000000 0.0000000 (PU - R)( 3) 1 2 3 4 5 1 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 2 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 3 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 4 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 6 7 8 9 1 0.0000000 0.0000000 0.0000000 0.0000000 2 0.0000000 0.0000000 0.0000000 0.0000000 3 0.0000000 0.0000000 0.0000000 0.0000000 4 0.0000000 0.0000000 0.0000000 0.0000000 COLRED TEST DRIVER RESULTS, EXAMPLE( 6) The input polynomial matrix: P(s) = P(0) + P(1) * s + . . . + P(dp-1) * s**(dp-1) + P(dp) * s**dp with degree DP = 1 and size MP = 2, NP = 2. The tolerance is: 0.000D+00 P( 0) 1 2 1 1.0000000 1.0000000 2 -2.0000000 -2.0000000 P( 1) 1 2 1 1.0000000 1.0000000 2 0.0000000 0.0000000 The unimodular polynomial matrix U(s): U( 0) 1 2 1 1.0000000 -1.0954451 2 -1.0000000 0.0000000 The column reduced polynomial matrix R(s): R( 0) 1 2 1 0.0000000 -1.0954451 2 0.0000000 2.1908902 R( 1) 1 2 1 0.0000000 -1.0954451 2 0.0000000 0.0000000 ZERCOL(j), j = 1, NP: T F The residual matrix P(s) * U(s) - R(s): (PU - R)( 0) 1 2 1 0.0000000 0.0000000 2 0.0000000 0.0000000 (PU - R)( 1) 1 2 1 0.0000000 0.0000000 2 0.0000000 0.0000000 COLRED TEST DRIVER RESULTS, EXAMPLE( 7) The input polynomial matrix: P(s) = P(0) + P(1) * s + . . . + P(dp-1) * s**(dp-1) + P(dp) * s**dp with degree DP = 0 and size MP = 3, NP = 2. The tolerance is: 0.000D+00 P( 0) 1 2 1 0.0000000 0.0000000 2 0.0000000 0.0000000 3 0.0000000 0.0000000 The unimodular polynomial matrix U(s): U( 0) 1 2 1 1.0000000 0.0000000 2 0.0000000 1.0000000 The column reduced polynomial matrix R(s): R( 0) 1 2 1 0.0000000 0.0000000 2 0.0000000 0.0000000 3 0.0000000 0.0000000 ZERCOL(j), j = 1, NP: F F PU - R is the ZERO polynomial matrix. COLRED TEST DRIVER RESULTS, EXAMPLE( 8) The input polynomial matrix: P(s) = P(0) + P(1) * s + . . . + P(dp-1) * s**(dp-1) + P(dp) * s**dp with degree DP = 0 and size MP = 1, NP = 2. The tolerance is: 0.000D+00 P( 0) 1 2 1 1.5000000 2.7000000 The unimodular polynomial matrix U(s): U( 0) 1 2 1 1.8000000 -0.6666667 2 -1.0000000 0.0000000 The column reduced polynomial matrix R(s): R( 0) 1 2 1 0.0000000 -1.0000000 ZERCOL(j), j = 1, NP: T F PU - R is the ZERO polynomial matrix. COLRED TEST DRIVER RESULTS, EXAMPLE( 9) The input polynomial matrix: P(s) = P(0) + P(1) * s + . . . + P(dp-1) * s**(dp-1) + P(dp) * s**dp with degree DP = 0 and size MP = 1, NP = 3. The tolerance is: 0.000D+00 P( 0) 1 2 3 1 1.2000000 1.5600000 -0.3610000 The unimodular polynomial matrix U(s): U( 0) 1 2 3 1 1.3000000 -0.3008333 -0.8333333 2 -1.0000000 0.0000000 0.0000000 3 0.0000000 -1.0000000 0.0000000 The column reduced polynomial matrix R(s): R( 0) 1 2 3 1 0.0000000 0.0000000 -1.0000000 ZERCOL(j), j = 1, NP: T T F PU - R is the ZERO polynomial matrix. COLRED TEST DRIVER RESULTS, EXAMPLE(10) The input polynomial matrix: P(s) = P(0) + P(1) * s + . . . + P(dp-1) * s**(dp-1) + P(dp) * s**dp with degree DP = 0 and size MP = 4, NP = 4. The tolerance is: 0.000D+00 P( 0) 1 2 3 4 1 0.0000000 0.0000000 1.0000000 1.0000000 2 1.0000000 1.0000000 2.0000000 2.0000000 3 0.0000000 0.0000000 0.0000000 0.0000000 4 1.0000000 1.0000000 -3.0000000 -3.0000000 The unimodular polynomial matrix U(s): U( 0) 1 2 3 4 1 1.0000000 0.0000000 0.7071068 0.1360828 2 -1.0000000 0.0000000 0.0000000 0.0000000 3 0.0000000 1.0000000 0.0000000 0.2721655 4 0.0000000 -1.0000000 0.0000000 0.0000000 The column reduced polynomial matrix R(s): R( 0) 1 2 3 4 1 0.0000000 0.0000000 0.0000000 0.2721655 2 0.0000000 0.0000000 0.7071068 0.6804138 3 0.0000000 0.0000000 0.0000000 0.0000000 4 0.0000000 0.0000000 0.7071068 -0.6804138 ZERCOL(j), j = 1, NP: T T F F The residual matrix P(s) * U(s) - R(s): (PU - R)( 0) 1 2 3 4 1 0.0000000 0.0000000 0.0000000 0.0000000 2 0.0000000 0.0000000 0.0000000 0.0000000 3 0.0000000 0.0000000 0.0000000 0.0000000 4 0.0000000 0.0000000 0.0000000 0.0000000 COLRED TEST DRIVER RESULTS, EXAMPLE(11) The input polynomial matrix: P(s) = P(0) + P(1) * s + . . . + P(dp-1) * s**(dp-1) + P(dp) * s**dp with degree DP = 0 and size MP = 0, NP = 0. The tolerance is: 0.000D+00 COLRED has failed: IERR = 1 SHAR_EOF fi # end of overwriting check cd .. cd .. if test ! -d 'Src' then mkdir 'Src' fi cd 'Src' if test ! -d 'Dp' then mkdir 'Dp' fi cd 'Dp' if test -f 'blas1.f' then echo shar: will not over-write existing file "'blas1.f'" else cat << \SHAR_EOF > 'blas1.f' subroutine daxpy(n,da,dx,incx,dy,incy) c c constant times a vector plus a vector. c uses unrolled loops for increments equal to one. c jack dongarra, linpack, 3/11/78. c modified 12/3/93, array(1) declarations changed to array(*) c double precision dx(*),dy(*),da integer i,incx,incy,ix,iy,m,mp1,n c if(n.le.0)return if (da .eq. 0.0d0) return if(incx.eq.1.and.incy.eq.1)go to 20 c c code for unequal increments or equal increments c not equal to 1 c ix = 1 iy = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 if(incy.lt.0)iy = (-n+1)*incy + 1 do 10 i = 1,n dy(iy) = dy(iy) + da*dx(ix) ix = ix + incx iy = iy + incy 10 continue return c c code for both increments equal to 1 c c c clean-up loop c 20 m = mod(n,4) if( m .eq. 0 ) go to 40 do 30 i = 1,m dy(i) = dy(i) + da*dx(i) 30 continue if( n .lt. 4 ) return 40 mp1 = m + 1 do 50 i = mp1,n,4 dy(i) = dy(i) + da*dx(i) dy(i + 1) = dy(i + 1) + da*dx(i + 1) dy(i + 2) = dy(i + 2) + da*dx(i + 2) dy(i + 3) = dy(i + 3) + da*dx(i + 3) 50 continue return end subroutine dcopy(n,dx,incx,dy,incy) c c copies a vector, x, to a vector, y. c uses unrolled loops for increments equal to one. c jack dongarra, linpack, 3/11/78. c modified 12/3/93, array(1) declarations changed to array(*) c double precision dx(*),dy(*) integer i,incx,incy,ix,iy,m,mp1,n c if(n.le.0)return if(incx.eq.1.and.incy.eq.1)go to 20 c c code for unequal increments or equal increments c not equal to 1 c ix = 1 iy = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 if(incy.lt.0)iy = (-n+1)*incy + 1 do 10 i = 1,n dy(iy) = dx(ix) ix = ix + incx iy = iy + incy 10 continue return c c code for both increments equal to 1 c c c clean-up loop c 20 m = mod(n,7) if( m .eq. 0 ) go to 40 do 30 i = 1,m dy(i) = dx(i) 30 continue if( n .lt. 7 ) return 40 mp1 = m + 1 do 50 i = mp1,n,7 dy(i) = dx(i) dy(i + 1) = dx(i + 1) dy(i + 2) = dx(i + 2) dy(i + 3) = dx(i + 3) dy(i + 4) = dx(i + 4) dy(i + 5) = dx(i + 5) dy(i + 6) = dx(i + 6) 50 continue return end double precision function ddot(n,dx,incx,dy,incy) c c forms the dot product of two vectors. c uses unrolled loops for increments equal to one. c jack dongarra, linpack, 3/11/78. c modified 12/3/93, array(1) declarations changed to array(*) c double precision dx(*),dy(*),dtemp integer i,incx,incy,ix,iy,m,mp1,n c ddot = 0.0d0 dtemp = 0.0d0 if(n.le.0)return if(incx.eq.1.and.incy.eq.1)go to 20 c c code for unequal increments or equal increments c not equal to 1 c ix = 1 iy = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 if(incy.lt.0)iy = (-n+1)*incy + 1 do 10 i = 1,n dtemp = dtemp + dx(ix)*dy(iy) ix = ix + incx iy = iy + incy 10 continue ddot = dtemp return c c code for both increments equal to 1 c c c clean-up loop c 20 m = mod(n,5) if( m .eq. 0 ) go to 40 do 30 i = 1,m dtemp = dtemp + dx(i)*dy(i) 30 continue if( n .lt. 5 ) go to 60 40 mp1 = m + 1 do 50 i = mp1,n,5 dtemp = dtemp + dx(i)*dy(i) + dx(i + 1)*dy(i + 1) + * dx(i + 2)*dy(i + 2) + dx(i + 3)*dy(i + 3) + dx(i + 4)*dy(i + 4) 50 continue 60 ddot = dtemp return end DOUBLE PRECISION FUNCTION DNRM2 ( N, X, INCX ) * .. Scalar Arguments .. INTEGER INCX, N * .. Array Arguments .. DOUBLE PRECISION X( * ) * .. * * DNRM2 returns the euclidean norm of a vector via the function * name, so that * * DNRM2 := sqrt( x'*x ) * * * * -- This version written on 25-October-1982. * Modified on 14-October-1993 to inline the call to DLASSQ. * Sven Hammarling, Nag Ltd. * * * .. Parameters .. DOUBLE PRECISION ONE , ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) * .. Local Scalars .. INTEGER IX DOUBLE PRECISION ABSXI, NORM, SCALE, SSQ * .. Intrinsic Functions .. INTRINSIC ABS, SQRT * .. * .. Executable Statements .. IF( N.LT.1 .OR. INCX.LT.1 )THEN NORM = ZERO ELSE IF( N.EQ.1 )THEN NORM = ABS( X( 1 ) ) ELSE SCALE = ZERO SSQ = ONE * The following loop is equivalent to this call to the LAPACK * auxiliary routine: * CALL DLASSQ( N, X, INCX, SCALE, SSQ ) * DO 10, IX = 1, 1 + ( N - 1 )*INCX, INCX IF( X( IX ).NE.ZERO )THEN ABSXI = ABS( X( IX ) ) IF( SCALE.LT.ABSXI )THEN SSQ = ONE + SSQ*( SCALE/ABSXI )**2 SCALE = ABSXI ELSE SSQ = SSQ + ( ABSXI/SCALE )**2 END IF END IF 10 CONTINUE NORM = SCALE * SQRT( SSQ ) END IF * DNRM2 = NORM RETURN * * End of DNRM2. * END subroutine dscal(n,da,dx,incx) c c scales a vector by a constant. c uses unrolled loops for increment equal to one. c jack dongarra, linpack, 3/11/78. c modified 3/93 to return if incx .le. 0. c modified 12/3/93, array(1) declarations changed to array(*) c double precision da,dx(*) integer i,incx,m,mp1,n,nincx c if( n.le.0 .or. incx.le.0 )return if(incx.eq.1)go to 20 c c code for increment not equal to 1 c nincx = n*incx do 10 i = 1,nincx,incx dx(i) = da*dx(i) 10 continue return c c code for increment equal to 1 c c c clean-up loop c 20 m = mod(n,5) if( m .eq. 0 ) go to 40 do 30 i = 1,m dx(i) = da*dx(i) 30 continue if( n .lt. 5 ) return 40 mp1 = m + 1 do 50 i = mp1,n,5 dx(i) = da*dx(i) dx(i + 1) = da*dx(i + 1) dx(i + 2) = da*dx(i + 2) dx(i + 3) = da*dx(i + 3) dx(i + 4) = da*dx(i + 4) 50 continue return end SHAR_EOF fi # end of overwriting check if test -f 'blas2.f' then echo shar: will not over-write existing file "'blas2.f'" else cat << \SHAR_EOF > 'blas2.f' SUBROUTINE DGEMV ( TRANS, M, N, ALPHA, A, LDA, X, INCX, $ BETA, Y, INCY ) * .. Scalar Arguments .. DOUBLE PRECISION ALPHA, BETA INTEGER INCX, INCY, LDA, M, N CHARACTER*1 TRANS * .. Array Arguments .. DOUBLE PRECISION A( LDA, * ), X( * ), Y( * ) * .. * * Purpose * ======= * * DGEMV performs one of the matrix-vector operations * * y := alpha*A*x + beta*y, or y := alpha*A'*x + beta*y, * * where alpha and beta are scalars, x and y are vectors and A is an * m by n matrix. * * Parameters * ========== * * TRANS - CHARACTER*1. * On entry, TRANS specifies the operation to be performed as * follows: * * TRANS = 'N' or 'n' y := alpha*A*x + beta*y. * * TRANS = 'T' or 't' y := alpha*A'*x + beta*y. * * TRANS = 'C' or 'c' y := alpha*A'*x + beta*y. * * Unchanged on exit. * * M - INTEGER. * On entry, M specifies the number of rows of the matrix A. * M must be at least zero. * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the number of columns of the matrix A. * N must be at least zero. * Unchanged on exit. * * ALPHA - DOUBLE PRECISION. * On entry, ALPHA specifies the scalar alpha. * Unchanged on exit. * * A - DOUBLE PRECISION array of DIMENSION ( LDA, n ). * Before entry, the leading m by n part of the array A must * contain the matrix of coefficients. * Unchanged on exit. * * LDA - INTEGER. * On entry, LDA specifies the first dimension of A as declared * in the calling (sub) program. LDA must be at least * max( 1, m ). * Unchanged on exit. * * X - DOUBLE PRECISION array of DIMENSION at least * ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n' * and at least * ( 1 + ( m - 1 )*abs( INCX ) ) otherwise. * Before entry, the incremented array X must contain the * vector x. * Unchanged on exit. * * INCX - INTEGER. * On entry, INCX specifies the increment for the elements of * X. INCX must not be zero. * Unchanged on exit. * * BETA - DOUBLE PRECISION. * On entry, BETA specifies the scalar beta. When BETA is * supplied as zero then Y need not be set on input. * Unchanged on exit. * * Y - DOUBLE PRECISION array of DIMENSION at least * ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n' * and at least * ( 1 + ( n - 1 )*abs( INCY ) ) otherwise. * Before entry with BETA non-zero, the incremented array Y * must contain the vector y. On exit, Y is overwritten by the * updated vector y. * * INCY - INTEGER. * On entry, INCY specifies the increment for the elements of * Y. INCY must not be zero. * Unchanged on exit. * * * Level 2 Blas routine. * * -- Written on 22-October-1986. * Jack Dongarra, Argonne National Lab. * Jeremy Du Croz, Nag Central Office. * Sven Hammarling, Nag Central Office. * Richard Hanson, Sandia National Labs. * * * .. Parameters .. DOUBLE PRECISION ONE , ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) * .. Local Scalars .. DOUBLE PRECISION TEMP INTEGER I, INFO, IX, IY, J, JX, JY, KX, KY, LENX, LENY * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. External Subroutines .. EXTERNAL XERBLA * .. Intrinsic Functions .. INTRINSIC MAX * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 IF ( .NOT.LSAME( TRANS, 'N' ).AND. $ .NOT.LSAME( TRANS, 'T' ).AND. $ .NOT.LSAME( TRANS, 'C' ) )THEN INFO = 1 ELSE IF( M.LT.0 )THEN INFO = 2 ELSE IF( N.LT.0 )THEN INFO = 3 ELSE IF( LDA.LT.MAX( 1, M ) )THEN INFO = 6 ELSE IF( INCX.EQ.0 )THEN INFO = 8 ELSE IF( INCY.EQ.0 )THEN INFO = 11 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'DGEMV ', INFO ) RETURN END IF * * Quick return if possible. * IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR. $ ( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) ) $ RETURN * * Set LENX and LENY, the lengths of the vectors x and y, and set * up the start points in X and Y. * IF( LSAME( TRANS, 'N' ) )THEN LENX = N LENY = M ELSE LENX = M LENY = N END IF IF( INCX.GT.0 )THEN KX = 1 ELSE KX = 1 - ( LENX - 1 )*INCX END IF IF( INCY.GT.0 )THEN KY = 1 ELSE KY = 1 - ( LENY - 1 )*INCY END IF * * Start the operations. In this version the elements of A are * accessed sequentially with one pass through A. * * First form y := beta*y. * IF( BETA.NE.ONE )THEN IF( INCY.EQ.1 )THEN IF( BETA.EQ.ZERO )THEN DO 10, I = 1, LENY Y( I ) = ZERO 10 CONTINUE ELSE DO 20, I = 1, LENY Y( I ) = BETA*Y( I ) 20 CONTINUE END IF ELSE IY = KY IF( BETA.EQ.ZERO )THEN DO 30, I = 1, LENY Y( IY ) = ZERO IY = IY + INCY 30 CONTINUE ELSE DO 40, I = 1, LENY Y( IY ) = BETA*Y( IY ) IY = IY + INCY 40 CONTINUE END IF END IF END IF IF( ALPHA.EQ.ZERO ) $ RETURN IF( LSAME( TRANS, 'N' ) )THEN * * Form y := alpha*A*x + y. * JX = KX IF( INCY.EQ.1 )THEN DO 60, J = 1, N IF( X( JX ).NE.ZERO )THEN TEMP = ALPHA*X( JX ) DO 50, I = 1, M Y( I ) = Y( I ) + TEMP*A( I, J ) 50 CONTINUE END IF JX = JX + INCX 60 CONTINUE ELSE DO 80, J = 1, N IF( X( JX ).NE.ZERO )THEN TEMP = ALPHA*X( JX ) IY = KY DO 70, I = 1, M Y( IY ) = Y( IY ) + TEMP*A( I, J ) IY = IY + INCY 70 CONTINUE END IF JX = JX + INCX 80 CONTINUE END IF ELSE * * Form y := alpha*A'*x + y. * JY = KY IF( INCX.EQ.1 )THEN DO 100, J = 1, N TEMP = ZERO DO 90, I = 1, M TEMP = TEMP + A( I, J )*X( I ) 90 CONTINUE Y( JY ) = Y( JY ) + ALPHA*TEMP JY = JY + INCY 100 CONTINUE ELSE DO 120, J = 1, N TEMP = ZERO IX = KX DO 110, I = 1, M TEMP = TEMP + A( I, J )*X( IX ) IX = IX + INCX 110 CONTINUE Y( JY ) = Y( JY ) + ALPHA*TEMP JY = JY + INCY 120 CONTINUE END IF END IF * RETURN * * End of DGEMV . * END SUBROUTINE DGER ( M, N, ALPHA, X, INCX, Y, INCY, A, LDA ) * .. Scalar Arguments .. DOUBLE PRECISION ALPHA INTEGER INCX, INCY, LDA, M, N * .. Array Arguments .. DOUBLE PRECISION A( LDA, * ), X( * ), Y( * ) * .. * * Purpose * ======= * * DGER performs the rank 1 operation * * A := alpha*x*y' + A, * * where alpha is a scalar, x is an m element vector, y is an n element * vector and A is an m by n matrix. * * Parameters * ========== * * M - INTEGER. * On entry, M specifies the number of rows of the matrix A. * M must be at least zero. * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the number of columns of the matrix A. * N must be at least zero. * Unchanged on exit. * * ALPHA - DOUBLE PRECISION. * On entry, ALPHA specifies the scalar alpha. * Unchanged on exit. * * X - DOUBLE PRECISION array of dimension at least * ( 1 + ( m - 1 )*abs( INCX ) ). * Before entry, the incremented array X must contain the m * element vector x. * Unchanged on exit. * * INCX - INTEGER. * On entry, INCX specifies the increment for the elements of * X. INCX must not be zero. * Unchanged on exit. * * Y - DOUBLE PRECISION array of dimension at least * ( 1 + ( n - 1 )*abs( INCY ) ). * Before entry, the incremented array Y must contain the n * element vector y. * Unchanged on exit. * * INCY - INTEGER. * On entry, INCY specifies the increment for the elements of * Y. INCY must not be zero. * Unchanged on exit. * * A - DOUBLE PRECISION array of DIMENSION ( LDA, n ). * Before entry, the leading m by n part of the array A must * contain the matrix of coefficients. On exit, A is * overwritten by the updated matrix. * * LDA - INTEGER. * On entry, LDA specifies the first dimension of A as declared * in the calling (sub) program. LDA must be at least * max( 1, m ). * Unchanged on exit. * * * Level 2 Blas routine. * * -- Written on 22-October-1986. * Jack Dongarra, Argonne National Lab. * Jeremy Du Croz, Nag Central Office. * Sven Hammarling, Nag Central Office. * Richard Hanson, Sandia National Labs. * * * .. Parameters .. DOUBLE PRECISION ZERO PARAMETER ( ZERO = 0.0D+0 ) * .. Local Scalars .. DOUBLE PRECISION TEMP INTEGER I, INFO, IX, J, JY, KX * .. External Subroutines .. EXTERNAL XERBLA * .. Intrinsic Functions .. INTRINSIC MAX * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 IF ( M.LT.0 )THEN INFO = 1 ELSE IF( N.LT.0 )THEN INFO = 2 ELSE IF( INCX.EQ.0 )THEN INFO = 5 ELSE IF( INCY.EQ.0 )THEN INFO = 7 ELSE IF( LDA.LT.MAX( 1, M ) )THEN INFO = 9 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'DGER ', INFO ) RETURN END IF * * Quick return if possible. * IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR.( ALPHA.EQ.ZERO ) ) $ RETURN * * Start the operations. In this version the elements of A are * accessed sequentially with one pass through A. * IF( INCY.GT.0 )THEN JY = 1 ELSE JY = 1 - ( N - 1 )*INCY END IF IF( INCX.EQ.1 )THEN DO 20, J = 1, N IF( Y( JY ).NE.ZERO )THEN TEMP = ALPHA*Y( JY ) DO 10, I = 1, M A( I, J ) = A( I, J ) + X( I )*TEMP 10 CONTINUE END IF JY = JY + INCY 20 CONTINUE ELSE IF( INCX.GT.0 )THEN KX = 1 ELSE KX = 1 - ( M - 1 )*INCX END IF DO 40, J = 1, N IF( Y( JY ).NE.ZERO )THEN TEMP = ALPHA*Y( JY ) IX = KX DO 30, I = 1, M A( I, J ) = A( I, J ) + X( IX )*TEMP IX = IX + INCX 30 CONTINUE END IF JY = JY + INCY 40 CONTINUE END IF * RETURN * * End of DGER . * END SHAR_EOF fi # end of overwriting check if test -f 'lapack.f' then echo shar: will not over-write existing file "'lapack.f'" else cat << \SHAR_EOF > 'lapack.f' SUBROUTINE DLACPY( UPLO, M, N, A, LDA, B, LDB ) * * -- LAPACK auxiliary routine (version 2.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * February 29, 1992 * * .. Scalar Arguments .. CHARACTER UPLO INTEGER LDA, LDB, M, N * .. * .. Array Arguments .. DOUBLE PRECISION A( LDA, * ), B( LDB, * ) * .. * * Purpose * ======= * * DLACPY copies all or part of a two-dimensional matrix A to another * matrix B. * * Arguments * ========= * * UPLO (input) CHARACTER*1 * Specifies the part of the matrix A to be copied to B. * = 'U': Upper triangular part * = 'L': Lower triangular part * Otherwise: All of the matrix A * * M (input) INTEGER * The number of rows of the matrix A. M >= 0. * * N (input) INTEGER * The number of columns of the matrix A. N >= 0. * * A (input) DOUBLE PRECISION array, dimension (LDA,N) * The m by n matrix A. If UPLO = 'U', only the upper triangle * or trapezoid is accessed; if UPLO = 'L', only the lower * triangle or trapezoid is accessed. * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,M). * * B (output) DOUBLE PRECISION array, dimension (LDB,N) * On exit, B = A in the locations specified by UPLO. * * LDB (input) INTEGER * The leading dimension of the array B. LDB >= max(1,M). * * ===================================================================== * * .. Local Scalars .. INTEGER I, J * .. * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. * .. Intrinsic Functions .. INTRINSIC MIN * .. * .. Executable Statements .. * IF( LSAME( UPLO, 'U' ) ) THEN DO 20 J = 1, N DO 10 I = 1, MIN( J, M ) B( I, J ) = A( I, J ) 10 CONTINUE 20 CONTINUE ELSE IF( LSAME( UPLO, 'L' ) ) THEN DO 40 J = 1, N DO 30 I = J, M B( I, J ) = A( I, J ) 30 CONTINUE 40 CONTINUE ELSE DO 60 J = 1, N DO 50 I = 1, M B( I, J ) = A( I, J ) 50 CONTINUE 60 CONTINUE END IF RETURN * * End of DLACPY * END DOUBLE PRECISION FUNCTION DLAMCH( CMACH ) * * -- LAPACK auxiliary routine (version 2.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * October 31, 1992 * * .. Scalar Arguments .. CHARACTER CMACH * .. * * Purpose * ======= * * DLAMCH determines double precision machine parameters. * * Arguments * ========= * * CMACH (input) CHARACTER*1 * Specifies the value to be returned by DLAMCH: * = 'E' or 'e', DLAMCH := eps * = 'S' or 's , DLAMCH := sfmin * = 'B' or 'b', DLAMCH := base * = 'P' or 'p', DLAMCH := eps*base * = 'N' or 'n', DLAMCH := t * = 'R' or 'r', DLAMCH := rnd * = 'M' or 'm', DLAMCH := emin * = 'U' or 'u', DLAMCH := rmin * = 'L' or 'l', DLAMCH := emax * = 'O' or 'o', DLAMCH := rmax * * where * * eps = relative machine precision * sfmin = safe minimum, such that 1/sfmin does not overflow * base = base of the machine * prec = eps*base * t = number of (base) digits in the mantissa * rnd = 1.0 when rounding occurs in addition, 0.0 otherwise * emin = minimum exponent before (gradual) underflow * rmin = underflow threshold - base**(emin-1) * emax = largest exponent before overflow * rmax = overflow threshold - (base**emax)*(1-eps) * * ===================================================================== * * .. Parameters .. DOUBLE PRECISION ONE, ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) * .. * .. Local Scalars .. LOGICAL FIRST, LRND INTEGER BETA, IMAX, IMIN, IT DOUBLE PRECISION BASE, EMAX, EMIN, EPS, PREC, RMACH, RMAX, RMIN, $ RND, SFMIN, SMALL, T * .. * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. * .. External Subroutines .. EXTERNAL DLAMC2 * .. * .. Save statement .. SAVE FIRST, EPS, SFMIN, BASE, T, RND, EMIN, RMIN, $ EMAX, RMAX, PREC * .. * .. Data statements .. DATA FIRST / .TRUE. / * .. * .. Executable Statements .. * IF( FIRST ) THEN FIRST = .FALSE. CALL DLAMC2( BETA, IT, LRND, EPS, IMIN, RMIN, IMAX, RMAX ) BASE = BETA T = IT IF( LRND ) THEN RND = ONE EPS = ( BASE**( 1-IT ) ) / 2 ELSE RND = ZERO EPS = BASE**( 1-IT ) END IF PREC = EPS*BASE EMIN = IMIN EMAX = IMAX SFMIN = RMIN SMALL = ONE / RMAX IF( SMALL.GE.SFMIN ) THEN * * Use SMALL plus a bit, to avoid the possibility of rounding * causing overflow when computing 1/sfmin. * SFMIN = SMALL*( ONE+EPS ) END IF END IF * IF( LSAME( CMACH, 'E' ) ) THEN RMACH = EPS ELSE IF( LSAME( CMACH, 'S' ) ) THEN RMACH = SFMIN ELSE IF( LSAME( CMACH, 'B' ) ) THEN RMACH = BASE ELSE IF( LSAME( CMACH, 'P' ) ) THEN RMACH = PREC ELSE IF( LSAME( CMACH, 'N' ) ) THEN RMACH = T ELSE IF( LSAME( CMACH, 'R' ) ) THEN RMACH = RND ELSE IF( LSAME( CMACH, 'M' ) ) THEN RMACH = EMIN ELSE IF( LSAME( CMACH, 'U' ) ) THEN RMACH = RMIN ELSE IF( LSAME( CMACH, 'L' ) ) THEN RMACH = EMAX ELSE IF( LSAME( CMACH, 'O' ) ) THEN RMACH = RMAX END IF * DLAMCH = RMACH RETURN * * End of DLAMCH * END * ************************************************************************ * SUBROUTINE DLAMC1( BETA, T, RND, IEEE1 ) * * -- LAPACK auxiliary routine (version 2.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * October 31, 1992 * * .. Scalar Arguments .. LOGICAL IEEE1, RND INTEGER BETA, T * .. * * Purpose * ======= * * DLAMC1 determines the machine parameters given by BETA, T, RND, and * IEEE1. * * Arguments * ========= * * BETA (output) INTEGER * The base of the machine. * * T (output) INTEGER * The number of ( BETA ) digits in the mantissa. * * RND (output) LOGICAL * Specifies whether proper rounding ( RND = .TRUE. ) or * chopping ( RND = .FALSE. ) occurs in addition. This may not * be a reliable guide to the way in which the machine performs * its arithmetic. * * IEEE1 (output) LOGICAL * Specifies whether rounding appears to be done in the IEEE * 'round to nearest' style. * * Further Details * =============== * * The routine is based on the routine ENVRON by Malcolm and * incorporates suggestions by Gentleman and Marovich. See * * Malcolm M. A. (1972) Algorithms to reveal properties of * floating-point arithmetic. Comms. of the ACM, 15, 949-951. * * Gentleman W. M. and Marovich S. B. (1974) More on algorithms * that reveal properties of floating point arithmetic units. * Comms. of the ACM, 17, 276-277. * * ===================================================================== * * .. Local Scalars .. LOGICAL FIRST, LIEEE1, LRND INTEGER LBETA, LT DOUBLE PRECISION A, B, C, F, ONE, QTR, SAVEC, T1, T2 * .. * .. External Functions .. DOUBLE PRECISION DLAMC3 EXTERNAL DLAMC3 * .. * .. Save statement .. SAVE FIRST, LIEEE1, LBETA, LRND, LT * .. * .. Data statements .. DATA FIRST / .TRUE. / * .. * .. Executable Statements .. * IF( FIRST ) THEN FIRST = .FALSE. ONE = 1 * * LBETA, LIEEE1, LT and LRND are the local values of BETA, * IEEE1, T and RND. * * Throughout this routine we use the function DLAMC3 to ensure * that relevant values are stored and not held in registers, or * are not affected by optimizers. * * Compute a = 2.0**m with the smallest positive integer m such * that * * fl( a + 1.0 ) = a. * A = 1 C = 1 * *+ WHILE( C.EQ.ONE )LOOP 10 CONTINUE IF( C.EQ.ONE ) THEN A = 2*A C = DLAMC3( A, ONE ) C = DLAMC3( C, -A ) GO TO 10 END IF *+ END WHILE * * Now compute b = 2.0**m with the smallest positive integer m * such that * * fl( a + b ) .gt. a. * B = 1 C = DLAMC3( A, B ) * *+ WHILE( C.EQ.A )LOOP 20 CONTINUE IF( C.EQ.A ) THEN B = 2*B C = DLAMC3( A, B ) GO TO 20 END IF *+ END WHILE * * Now compute the base. a and c are neighbouring floating point * numbers in the interval ( beta**t, beta**( t + 1 ) ) and so * their difference is beta. Adding 0.25 to c is to ensure that it * is truncated to beta and not ( beta - 1 ). * QTR = ONE / 4 SAVEC = C C = DLAMC3( C, -A ) LBETA = C + QTR * * Now determine whether rounding or chopping occurs, by adding a * bit less than beta/2 and a bit more than beta/2 to a. * B = LBETA F = DLAMC3( B / 2, -B / 100 ) C = DLAMC3( F, A ) IF( C.EQ.A ) THEN LRND = .TRUE. ELSE LRND = .FALSE. END IF F = DLAMC3( B / 2, B / 100 ) C = DLAMC3( F, A ) IF( ( LRND ) .AND. ( C.EQ.A ) ) $ LRND = .FALSE. * * Try and decide whether rounding is done in the IEEE 'round to * nearest' style. B/2 is half a unit in the last place of the two * numbers A and SAVEC. Furthermore, A is even, i.e. has last bit * zero, and SAVEC is odd. Thus adding B/2 to A should not change * A, but adding B/2 to SAVEC should change SAVEC. * T1 = DLAMC3( B / 2, A ) T2 = DLAMC3( B / 2, SAVEC ) LIEEE1 = ( T1.EQ.A ) .AND. ( T2.GT.SAVEC ) .AND. LRND * * Now find the mantissa, t. It should be the integer part of * log to the base beta of a, however it is safer to determine t * by powering. So we find t as the smallest positive integer for * which * * fl( beta**t + 1.0 ) = 1.0. * LT = 0 A = 1 C = 1 * *+ WHILE( C.EQ.ONE )LOOP 30 CONTINUE IF( C.EQ.ONE ) THEN LT = LT + 1 A = A*LBETA C = DLAMC3( A, ONE ) C = DLAMC3( C, -A ) GO TO 30 END IF *+ END WHILE * END IF * BETA = LBETA T = LT RND = LRND IEEE1 = LIEEE1 RETURN * * End of DLAMC1 * END * ************************************************************************ * SUBROUTINE DLAMC2( BETA, T, RND, EPS, EMIN, RMIN, EMAX, RMAX ) * * -- LAPACK auxiliary routine (version 2.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * October 31, 1992 * * .. Scalar Arguments .. LOGICAL RND INTEGER BETA, EMAX, EMIN, T DOUBLE PRECISION EPS, RMAX, RMIN * .. * * Purpose * ======= * * DLAMC2 determines the machine parameters specified in its argument * list. * * Arguments * ========= * * BETA (output) INTEGER * The base of the machine. * * T (output) INTEGER * The number of ( BETA ) digits in the mantissa. * * RND (output) LOGICAL * Specifies whether proper rounding ( RND = .TRUE. ) or * chopping ( RND = .FALSE. ) occurs in addition. This may not * be a reliable guide to the way in which the machine performs * its arithmetic. * * EPS (output) DOUBLE PRECISION * The smallest positive number such that * * fl( 1.0 - EPS ) .LT. 1.0, * * where fl denotes the computed value. * * EMIN (output) INTEGER * The minimum exponent before (gradual) underflow occurs. * * RMIN (output) DOUBLE PRECISION * The smallest normalized number for the machine, given by * BASE**( EMIN - 1 ), where BASE is the floating point value * of BETA. * * EMAX (output) INTEGER * The maximum exponent before overflow occurs. * * RMAX (output) DOUBLE PRECISION * The largest positive number for the machine, given by * BASE**EMAX * ( 1 - EPS ), where BASE is the floating point * value of BETA. * * Further Details * =============== * * The computation of EPS is based on a routine PARANOIA by * W. Kahan of the University of California at Berkeley. * * ===================================================================== * * .. Local Scalars .. LOGICAL FIRST, IEEE, IWARN, LIEEE1, LRND INTEGER GNMIN, GPMIN, I, LBETA, LEMAX, LEMIN, LT, $ NGNMIN, NGPMIN DOUBLE PRECISION A, B, C, HALF, LEPS, LRMAX, LRMIN, ONE, RBASE, $ SIXTH, SMALL, THIRD, TWO, ZERO * .. * .. External Functions .. DOUBLE PRECISION DLAMC3 EXTERNAL DLAMC3 * .. * .. External Subroutines .. EXTERNAL DLAMC1, DLAMC4, DLAMC5 * .. * .. Intrinsic Functions .. INTRINSIC ABS, MAX, MIN * .. * .. Save statement .. SAVE FIRST, IWARN, LBETA, LEMAX, LEMIN, LEPS, LRMAX, $ LRMIN, LT * .. * .. Data statements .. DATA FIRST / .TRUE. / , IWARN / .FALSE. / * .. * .. Executable Statements .. * IF( FIRST ) THEN FIRST = .FALSE. ZERO = 0 ONE = 1 TWO = 2 * * LBETA, LT, LRND, LEPS, LEMIN and LRMIN are the local values of * BETA, T, RND, EPS, EMIN and RMIN. * * Throughout this routine we use the function DLAMC3 to ensure * that relevant values are stored and not held in registers, or * are not affected by optimizers. * * DLAMC1 returns the parameters LBETA, LT, LRND and LIEEE1. * CALL DLAMC1( LBETA, LT, LRND, LIEEE1 ) * * Start to find EPS. * B = LBETA A = B**( -LT ) LEPS = A * * Try some tricks to see whether or not this is the correct EPS. * B = TWO / 3 HALF = ONE / 2 SIXTH = DLAMC3( B, -HALF ) THIRD = DLAMC3( SIXTH, SIXTH ) B = DLAMC3( THIRD, -HALF ) B = DLAMC3( B, SIXTH ) B = ABS( B ) IF( B.LT.LEPS ) $ B = LEPS * LEPS = 1 * *+ WHILE( ( LEPS.GT.B ).AND.( B.GT.ZERO ) )LOOP 10 CONTINUE IF( ( LEPS.GT.B ) .AND. ( B.GT.ZERO ) ) THEN LEPS = B C = DLAMC3( HALF*LEPS, ( TWO**5 )*( LEPS**2 ) ) C = DLAMC3( HALF, -C ) B = DLAMC3( HALF, C ) C = DLAMC3( HALF, -B ) B = DLAMC3( HALF, C ) GO TO 10 END IF *+ END WHILE * IF( A.LT.LEPS ) $ LEPS = A * * Computation of EPS complete. * * Now find EMIN. Let A = + or - 1, and + or - (1 + BASE**(-3)). * Keep dividing A by BETA until (gradual) underflow occurs. This * is detected when we cannot recover the previous A. * RBASE = ONE / LBETA SMALL = ONE DO 20 I = 1, 3 SMALL = DLAMC3( SMALL*RBASE, ZERO ) 20 CONTINUE A = DLAMC3( ONE, SMALL ) CALL DLAMC4( NGPMIN, ONE, LBETA ) CALL DLAMC4( NGNMIN, -ONE, LBETA ) CALL DLAMC4( GPMIN, A, LBETA ) CALL DLAMC4( GNMIN, -A, LBETA ) IEEE = .FALSE. * IF( ( NGPMIN.EQ.NGNMIN ) .AND. ( GPMIN.EQ.GNMIN ) ) THEN IF( NGPMIN.EQ.GPMIN ) THEN LEMIN = NGPMIN * ( Non twos-complement machines, no gradual underflow; * e.g., VAX ) ELSE IF( ( GPMIN-NGPMIN ).EQ.3 ) THEN LEMIN = NGPMIN - 1 + LT IEEE = .TRUE. * ( Non twos-complement machines, with gradual underflow; * e.g., IEEE standard followers ) ELSE LEMIN = MIN( NGPMIN, GPMIN ) * ( A guess; no known machine ) IWARN = .TRUE. END IF * ELSE IF( ( NGPMIN.EQ.GPMIN ) .AND. ( NGNMIN.EQ.GNMIN ) ) THEN IF( ABS( NGPMIN-NGNMIN ).EQ.1 ) THEN LEMIN = MAX( NGPMIN, NGNMIN ) * ( Twos-complement machines, no gradual underflow; * e.g., CYBER 205 ) ELSE LEMIN = MIN( NGPMIN, NGNMIN ) * ( A guess; no known machine ) IWARN = .TRUE. END IF * ELSE IF( ( ABS( NGPMIN-NGNMIN ).EQ.1 ) .AND. $ ( GPMIN.EQ.GNMIN ) ) THEN IF( ( GPMIN-MIN( NGPMIN, NGNMIN ) ).EQ.3 ) THEN LEMIN = MAX( NGPMIN, NGNMIN ) - 1 + LT * ( Twos-complement machines with gradual underflow; * no known machine ) ELSE LEMIN = MIN( NGPMIN, NGNMIN ) * ( A guess; no known machine ) IWARN = .TRUE. END IF * ELSE LEMIN = MIN( NGPMIN, NGNMIN, GPMIN, GNMIN ) * ( A guess; no known machine ) IWARN = .TRUE. END IF *** * Comment out this if block if EMIN is ok IF( IWARN ) THEN FIRST = .TRUE. WRITE( 6, FMT = 9999 )LEMIN END IF *** * * Assume IEEE arithmetic if we found denormalised numbers above, * or if arithmetic seems to round in the IEEE style, determined * in routine DLAMC1. A true IEEE machine should have both things * true; however, faulty machines may have one or the other. * IEEE = IEEE .OR. LIEEE1 * * Compute RMIN by successive division by BETA. We could compute * RMIN as BASE**( EMIN - 1 ), but some machines underflow during * this computation. * LRMIN = 1 DO 30 I = 1, 1 - LEMIN LRMIN = DLAMC3( LRMIN*RBASE, ZERO ) 30 CONTINUE * * Finally, call DLAMC5 to compute EMAX and RMAX. * CALL DLAMC5( LBETA, LT, LEMIN, IEEE, LEMAX, LRMAX ) END IF * BETA = LBETA T = LT RND = LRND EPS = LEPS EMIN = LEMIN RMIN = LRMIN EMAX = LEMAX RMAX = LRMAX * RETURN * 9999 FORMAT( / / ' WARNING. The value EMIN may be incorrect:-', $ ' EMIN = ', I8, / $ ' If, after inspection, the value EMIN looks', $ ' acceptable please comment out ', $ / ' the IF block as marked within the code of routine', $ ' DLAMC2,', / ' otherwise supply EMIN explicitly.', / ) * * End of DLAMC2 * END * ************************************************************************ * DOUBLE PRECISION FUNCTION DLAMC3( A, B ) * * -- LAPACK auxiliary routine (version 2.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * October 31, 1992 * * .. Scalar Arguments .. DOUBLE PRECISION A, B * .. * * Purpose * ======= * * DLAMC3 is intended to force A and B to be stored prior to doing * the addition of A and B , for use in situations where optimizers * might hold one of these in a register. * * Arguments * ========= * * A, B (input) DOUBLE PRECISION * The values A and B. * * ===================================================================== * * .. Executable Statements .. * DLAMC3 = A + B * RETURN * * End of DLAMC3 * END * ************************************************************************ * SUBROUTINE DLAMC4( EMIN, START, BASE ) * * -- LAPACK auxiliary routine (version 2.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * October 31, 1992 * * .. Scalar Arguments .. INTEGER BASE, EMIN DOUBLE PRECISION START * .. * * Purpose * ======= * * DLAMC4 is a service routine for DLAMC2. * * Arguments * ========= * * EMIN (output) EMIN * The minimum exponent before (gradual) underflow, computed by * setting A = START and dividing by BASE until the previous A * can not be recovered. * * START (input) DOUBLE PRECISION * The starting point for determining EMIN. * * BASE (input) INTEGER * The base of the machine. * * ===================================================================== * * .. Local Scalars .. INTEGER I DOUBLE PRECISION A, B1, B2, C1, C2, D1, D2, ONE, RBASE, ZERO * .. * .. External Functions .. DOUBLE PRECISION DLAMC3 EXTERNAL DLAMC3 * .. * .. Executable Statements .. * A = START ONE = 1 RBASE = ONE / BASE ZERO = 0 EMIN = 1 B1 = DLAMC3( A*RBASE, ZERO ) C1 = A C2 = A D1 = A D2 = A *+ WHILE( ( C1.EQ.A ).AND.( C2.EQ.A ).AND. * $ ( D1.EQ.A ).AND.( D2.EQ.A ) )LOOP 10 CONTINUE IF( ( C1.EQ.A ) .AND. ( C2.EQ.A ) .AND. ( D1.EQ.A ) .AND. $ ( D2.EQ.A ) ) THEN EMIN = EMIN - 1 A = B1 B1 = DLAMC3( A / BASE, ZERO ) C1 = DLAMC3( B1*BASE, ZERO ) D1 = ZERO DO 20 I = 1, BASE D1 = D1 + B1 20 CONTINUE B2 = DLAMC3( A*RBASE, ZERO ) C2 = DLAMC3( B2 / RBASE, ZERO ) D2 = ZERO DO 30 I = 1, BASE D2 = D2 + B2 30 CONTINUE GO TO 10 END IF *+ END WHILE * RETURN * * End of DLAMC4 * END * ************************************************************************ * SUBROUTINE DLAMC5( BETA, P, EMIN, IEEE, EMAX, RMAX ) * * -- LAPACK auxiliary routine (version 2.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * October 31, 1992 * * .. Scalar Arguments .. LOGICAL IEEE INTEGER BETA, EMAX, EMIN, P DOUBLE PRECISION RMAX * .. * * Purpose * ======= * * DLAMC5 attempts to compute RMAX, the largest machine floating-point * number, without overflow. It assumes that EMAX + abs(EMIN) sum * approximately to a power of 2. It will fail on machines where this * assumption does not hold, for example, the Cyber 205 (EMIN = -28625, * EMAX = 28718). It will also fail if the value supplied for EMIN is * too large (i.e. too close to zero), probably with overflow. * * Arguments * ========= * * BETA (input) INTEGER * The base of floating-point arithmetic. * * P (input) INTEGER * The number of base BETA digits in the mantissa of a * floating-point value. * * EMIN (input) INTEGER * The minimum exponent before (gradual) underflow. * * IEEE (input) LOGICAL * A logical flag specifying whether or not the arithmetic * system is thought to comply with the IEEE standard. * * EMAX (output) INTEGER * The largest exponent before overflow * * RMAX (output) DOUBLE PRECISION * The largest machine floating-point number. * * ===================================================================== * * .. Parameters .. DOUBLE PRECISION ZERO, ONE PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) * .. * .. Local Scalars .. INTEGER EXBITS, EXPSUM, I, LEXP, NBITS, TRY, UEXP DOUBLE PRECISION OLDY, RECBAS, Y, Z * .. * .. External Functions .. DOUBLE PRECISION DLAMC3 EXTERNAL DLAMC3 * .. * .. Intrinsic Functions .. INTRINSIC MOD * .. * .. Executable Statements .. * * First compute LEXP and UEXP, two powers of 2 that bound * abs(EMIN). We then assume that EMAX + abs(EMIN) will sum * approximately to the bound that is closest to abs(EMIN). * (EMAX is the exponent of the required number RMAX). * LEXP = 1 EXBITS = 1 10 CONTINUE TRY = LEXP*2 IF( TRY.LE.( -EMIN ) ) THEN LEXP = TRY EXBITS = EXBITS + 1 GO TO 10 END IF IF( LEXP.EQ.-EMIN ) THEN UEXP = LEXP ELSE UEXP = TRY EXBITS = EXBITS + 1 END IF * * Now -LEXP is less than or equal to EMIN, and -UEXP is greater * than or equal to EMIN. EXBITS is the number of bits needed to * store the exponent. * IF( ( UEXP+EMIN ).GT.( -LEXP-EMIN ) ) THEN EXPSUM = 2*LEXP ELSE EXPSUM = 2*UEXP END IF * * EXPSUM is the exponent range, approximately equal to * EMAX - EMIN + 1 . * EMAX = EXPSUM + EMIN - 1 NBITS = 1 + EXBITS + P * * NBITS is the total number of bits needed to store a * floating-point number. * IF( ( MOD( NBITS, 2 ).EQ.1 ) .AND. ( BETA.EQ.2 ) ) THEN * * Either there are an odd number of bits used to store a * floating-point number, which is unlikely, or some bits are * not used in the representation of numbers, which is possible, * (e.g. Cray machines) or the mantissa has an implicit bit, * (e.g. IEEE machines, Dec Vax machines), which is perhaps the * most likely. We have to assume the last alternative. * If this is true, then we need to reduce EMAX by one because * there must be some way of representing zero in an implicit-bit * system. On machines like Cray, we are reducing EMAX by one * unnecessarily. * EMAX = EMAX - 1 END IF * IF( IEEE ) THEN * * Assume we are on an IEEE machine which reserves one exponent * for infinity and NaN. * EMAX = EMAX - 1 END IF * * Now create RMAX, the largest machine number, which should * be equal to (1.0 - BETA**(-P)) * BETA**EMAX . * * First compute 1.0 - BETA**(-P), being careful that the * result is less than 1.0 . * RECBAS = ONE / BETA Z = BETA - ONE Y = ZERO DO 20 I = 1, P Z = Z*RECBAS IF( Y.LT.ONE ) $ OLDY = Y Y = DLAMC3( Y, Z ) 20 CONTINUE IF( Y.GE.ONE ) $ Y = OLDY * * Now multiply by BETA**EMAX to get RMAX. * DO 30 I = 1, EMAX Y = DLAMC3( Y*BETA, ZERO ) 30 CONTINUE * RMAX = Y RETURN * * End of DLAMC5 * END DOUBLE PRECISION FUNCTION DLANGE( NORM, M, N, A, LDA, WORK ) * * -- LAPACK auxiliary routine (version 2.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * October 31, 1992 * * .. Scalar Arguments .. CHARACTER NORM INTEGER LDA, M, N * .. * .. Array Arguments .. DOUBLE PRECISION A( LDA, * ), WORK( * ) * .. * * Purpose * ======= * * DLANGE returns the value of the one norm, or the Frobenius norm, or * the infinity norm, or the element of largest absolute value of a * real matrix A. * * Description * =========== * * DLANGE returns the value * * DLANGE = ( max(abs(A(i,j))), NORM = 'M' or 'm' * ( * ( norm1(A), NORM = '1', 'O' or 'o' * ( * ( normI(A), NORM = 'I' or 'i' * ( * ( normF(A), NORM = 'F', 'f', 'E' or 'e' * * where norm1 denotes the one norm of a matrix (maximum column sum), * normI denotes the infinity norm of a matrix (maximum row sum) and * normF denotes the Frobenius norm of a matrix (square root of sum of * squares). Note that max(abs(A(i,j))) is not a matrix norm. * * Arguments * ========= * * NORM (input) CHARACTER*1 * Specifies the value to be returned in DLANGE as described * above. * * M (input) INTEGER * The number of rows of the matrix A. M >= 0. When M = 0, * DLANGE is set to zero. * * N (input) INTEGER * The number of columns of the matrix A. N >= 0. When N = 0, * DLANGE is set to zero. * * A (input) DOUBLE PRECISION array, dimension (LDA,N) * The m by n matrix A. * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(M,1). * * WORK (workspace) DOUBLE PRECISION array, dimension (LWORK), * where LWORK >= M when NORM = 'I'; otherwise, WORK is not * referenced. * * ===================================================================== * * .. Parameters .. DOUBLE PRECISION ONE, ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) * .. * .. Local Scalars .. INTEGER I, J DOUBLE PRECISION SCALE, SUM, VALUE * .. * .. External Subroutines .. EXTERNAL DLASSQ * .. * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. * .. Intrinsic Functions .. INTRINSIC ABS, MAX, MIN, SQRT * .. * .. Executable Statements .. * IF( MIN( M, N ).EQ.0 ) THEN VALUE = ZERO ELSE IF( LSAME( NORM, 'M' ) ) THEN * * Find max(abs(A(i,j))). * VALUE = ZERO DO 20 J = 1, N DO 10 I = 1, M VALUE = MAX( VALUE, ABS( A( I, J ) ) ) 10 CONTINUE 20 CONTINUE ELSE IF( ( LSAME( NORM, 'O' ) ) .OR. ( NORM.EQ.'1' ) ) THEN * * Find norm1(A). * VALUE = ZERO DO 40 J = 1, N SUM = ZERO DO 30 I = 1, M SUM = SUM + ABS( A( I, J ) ) 30 CONTINUE VALUE = MAX( VALUE, SUM ) 40 CONTINUE ELSE IF( LSAME( NORM, 'I' ) ) THEN * * Find normI(A). * DO 50 I = 1, M WORK( I ) = ZERO 50 CONTINUE DO 70 J = 1, N DO 60 I = 1, M WORK( I ) = WORK( I ) + ABS( A( I, J ) ) 60 CONTINUE 70 CONTINUE VALUE = ZERO DO 80 I = 1, M VALUE = MAX( VALUE, WORK( I ) ) 80 CONTINUE ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN * * Find normF(A). * SCALE = ZERO SUM = ONE DO 90 J = 1, N CALL DLASSQ( M, A( 1, J ), 1, SCALE, SUM ) 90 CONTINUE VALUE = SCALE*SQRT( SUM ) END IF * DLANGE = VALUE RETURN * * End of DLANGE * END DOUBLE PRECISION FUNCTION DLAPY2( X, Y ) * * -- LAPACK auxiliary routine (version 2.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * October 31, 1992 * * .. Scalar Arguments .. DOUBLE PRECISION X, Y * .. * * Purpose * ======= * * DLAPY2 returns sqrt(x**2+y**2), taking care not to cause unnecessary * overflow. * * Arguments * ========= * * X (input) DOUBLE PRECISION * Y (input) DOUBLE PRECISION * X and Y specify the values x and y. * * ===================================================================== * * .. Parameters .. DOUBLE PRECISION ZERO PARAMETER ( ZERO = 0.0D0 ) DOUBLE PRECISION ONE PARAMETER ( ONE = 1.0D0 ) * .. * .. Local Scalars .. DOUBLE PRECISION W, XABS, YABS, Z * .. * .. Intrinsic Functions .. INTRINSIC ABS, MAX, MIN, SQRT * .. * .. Executable Statements .. * XABS = ABS( X ) YABS = ABS( Y ) W = MAX( XABS, YABS ) Z = MIN( XABS, YABS ) IF( Z.EQ.ZERO ) THEN DLAPY2 = W ELSE DLAPY2 = W*SQRT( ONE+( Z / W )**2 ) END IF RETURN * * End of DLAPY2 * END SUBROUTINE DLARFG( N, ALPHA, X, INCX, TAU ) * * -- LAPACK auxiliary routine (version 2.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * September 30, 1994 * * .. Scalar Arguments .. INTEGER INCX, N DOUBLE PRECISION ALPHA, TAU * .. * .. Array Arguments .. DOUBLE PRECISION X( * ) * .. * * Purpose * ======= * * DLARFG generates a real elementary reflector H of order n, such * that * * H * ( alpha ) = ( beta ), H' * H = I. * ( x ) ( 0 ) * * where alpha and beta are scalars, and x is an (n-1)-element real * vector. H is represented in the form * * H = I - tau * ( 1 ) * ( 1 v' ) , * ( v ) * * where tau is a real scalar and v is a real (n-1)-element * vector. * * If the elements of x are all zero, then tau = 0 and H is taken to be * the unit matrix. * * Otherwise 1 <= tau <= 2. * * Arguments * ========= * * N (input) INTEGER * The order of the elementary reflector. * * ALPHA (input/output) DOUBLE PRECISION * On entry, the value alpha. * On exit, it is overwritten with the value beta. * * X (input/output) DOUBLE PRECISION array, dimension * (1+(N-2)*abs(INCX)) * On entry, the vector x. * On exit, it is overwritten with the vector v. * * INCX (input) INTEGER * The increment between elements of X. INCX > 0. * * TAU (output) DOUBLE PRECISION * The value tau. * * ===================================================================== * * .. Parameters .. DOUBLE PRECISION ONE, ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) * .. * .. Local Scalars .. INTEGER J, KNT DOUBLE PRECISION BETA, RSAFMN, SAFMIN, XNORM * .. * .. External Functions .. DOUBLE PRECISION DLAMCH, DLAPY2, DNRM2 EXTERNAL DLAMCH, DLAPY2, DNRM2 * .. * .. Intrinsic Functions .. INTRINSIC ABS, SIGN * .. * .. External Subroutines .. EXTERNAL DSCAL * .. * .. Executable Statements .. * IF( N.LE.1 ) THEN TAU = ZERO RETURN END IF * XNORM = DNRM2( N-1, X, INCX ) * IF( XNORM.EQ.ZERO ) THEN * * H = I * TAU = ZERO ELSE * * general case * BETA = -SIGN( DLAPY2( ALPHA, XNORM ), ALPHA ) SAFMIN = DLAMCH( 'S' ) / DLAMCH( 'E' ) IF( ABS( BETA ).LT.SAFMIN ) THEN * * XNORM, BETA may be inaccurate; scale X and recompute them * RSAFMN = ONE / SAFMIN KNT = 0 10 CONTINUE KNT = KNT + 1 CALL DSCAL( N-1, RSAFMN, X, INCX ) BETA = BETA*RSAFMN ALPHA = ALPHA*RSAFMN IF( ABS( BETA ).LT.SAFMIN ) $ GO TO 10 * * New BETA is at most 1, at least SAFMIN * XNORM = DNRM2( N-1, X, INCX ) BETA = -SIGN( DLAPY2( ALPHA, XNORM ), ALPHA ) TAU = ( BETA-ALPHA ) / BETA CALL DSCAL( N-1, ONE / ( ALPHA-BETA ), X, INCX ) * * If ALPHA is subnormal, it may lose relative accuracy * ALPHA = BETA DO 20 J = 1, KNT ALPHA = ALPHA*SAFMIN 20 CONTINUE ELSE TAU = ( BETA-ALPHA ) / BETA CALL DSCAL( N-1, ONE / ( ALPHA-BETA ), X, INCX ) ALPHA = BETA END IF END IF * RETURN * * End of DLARFG * END SUBROUTINE DLASET( UPLO, M, N, ALPHA, BETA, A, LDA ) * * -- LAPACK auxiliary routine (version 2.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * October 31, 1992 * * .. Scalar Arguments .. CHARACTER UPLO INTEGER LDA, M, N DOUBLE PRECISION ALPHA, BETA * .. * .. Array Arguments .. DOUBLE PRECISION A( LDA, * ) * .. * * Purpose * ======= * * DLASET initializes an m-by-n matrix A to BETA on the diagonal and * ALPHA on the offdiagonals. * * Arguments * ========= * * UPLO (input) CHARACTER*1 * Specifies the part of the matrix A to be set. * = 'U': Upper triangular part is set; the strictly lower * triangular part of A is not changed. * = 'L': Lower triangular part is set; the strictly upper * triangular part of A is not changed. * Otherwise: All of the matrix A is set. * * M (input) INTEGER * The number of rows of the matrix A. M >= 0. * * N (input) INTEGER * The number of columns of the matrix A. N >= 0. * * ALPHA (input) DOUBLE PRECISION * The constant to which the offdiagonal elements are to be set. * * BETA (input) DOUBLE PRECISION * The constant to which the diagonal elements are to be set. * * A (input/output) DOUBLE PRECISION array, dimension (LDA,N) * On exit, the leading m-by-n submatrix of A is set as follows: * * if UPLO = 'U', A(i,j) = ALPHA, 1<=i<=j-1, 1<=j<=n, * if UPLO = 'L', A(i,j) = ALPHA, j+1<=i<=m, 1<=j<=n, * otherwise, A(i,j) = ALPHA, 1<=i<=m, 1<=j<=n, i.ne.j, * * and, for all UPLO, A(i,i) = BETA, 1<=i<=min(m,n). * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,M). * * ===================================================================== * * .. Local Scalars .. INTEGER I, J * .. * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. * .. Intrinsic Functions .. INTRINSIC MIN * .. * .. Executable Statements .. * IF( LSAME( UPLO, 'U' ) ) THEN * * Set the strictly upper triangular or trapezoidal part of the * array to ALPHA. * DO 20 J = 2, N DO 10 I = 1, MIN( J-1, M ) A( I, J ) = ALPHA 10 CONTINUE 20 CONTINUE * ELSE IF( LSAME( UPLO, 'L' ) ) THEN * * Set the strictly lower triangular or trapezoidal part of the * array to ALPHA. * DO 40 J = 1, MIN( M, N ) DO 30 I = J + 1, M A( I, J ) = ALPHA 30 CONTINUE 40 CONTINUE * ELSE * * Set the leading m-by-n submatrix to ALPHA. * DO 60 J = 1, N DO 50 I = 1, M A( I, J ) = ALPHA 50 CONTINUE 60 CONTINUE END IF * * Set the first min(M,N) diagonal elements to BETA. * DO 70 I = 1, MIN( M, N ) A( I, I ) = BETA 70 CONTINUE * RETURN * * End of DLASET * END SUBROUTINE DLASSQ( N, X, INCX, SCALE, SUMSQ ) * * -- LAPACK auxiliary routine (version 2.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * October 31, 1992 * * .. Scalar Arguments .. INTEGER INCX, N DOUBLE PRECISION SCALE, SUMSQ * .. * .. Array Arguments .. DOUBLE PRECISION X( * ) * .. * * Purpose * ======= * * DLASSQ returns the values scl and smsq such that * * ( scl**2 )*smsq = x( 1 )**2 +...+ x( n )**2 + ( scale**2 )*sumsq, * * where x( i ) = X( 1 + ( i - 1 )*INCX ). The value of sumsq is * assumed to be non-negative and scl returns the value * * scl = max( scale, abs( x( i ) ) ). * * scale and sumsq must be supplied in SCALE and SUMSQ and * scl and smsq are overwritten on SCALE and SUMSQ respectively. * * The routine makes only one pass through the vector x. * * Arguments * ========= * * N (input) INTEGER * The number of elements to be used from the vector X. * * X (input) DOUBLE PRECISION * The vector for which a scaled sum of squares is computed. * x( i ) = X( 1 + ( i - 1 )*INCX ), 1 <= i <= n. * * INCX (input) INTEGER * The increment between successive values of the vector X. * INCX > 0. * * SCALE (input/output) DOUBLE PRECISION * On entry, the value scale in the equation above. * On exit, SCALE is overwritten with scl , the scaling factor * for the sum of squares. * * SUMSQ (input/output) DOUBLE PRECISION * On entry, the value sumsq in the equation above. * On exit, SUMSQ is overwritten with smsq , the basic sum of * squares from which scl has been factored out. * * ===================================================================== * * .. Parameters .. DOUBLE PRECISION ZERO PARAMETER ( ZERO = 0.0D+0 ) * .. * .. Local Scalars .. INTEGER IX DOUBLE PRECISION ABSXI * .. * .. Intrinsic Functions .. INTRINSIC ABS * .. * .. Executable Statements .. * IF( N.GT.0 ) THEN DO 10 IX = 1, 1 + ( N-1 )*INCX, INCX IF( X( IX ).NE.ZERO ) THEN ABSXI = ABS( X( IX ) ) IF( SCALE.LT.ABSXI ) THEN SUMSQ = 1 + SUMSQ*( SCALE / ABSXI )**2 SCALE = ABSXI ELSE SUMSQ = SUMSQ + ( ABSXI / SCALE )**2 END IF END IF 10 CONTINUE END IF RETURN * * End of DLASSQ * END LOGICAL FUNCTION LSAME( CA, CB ) * * -- LAPACK auxiliary routine (version 2.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * September 30, 1994 * * .. Scalar Arguments .. CHARACTER CA, CB * .. * * Purpose * ======= * * LSAME returns .TRUE. if CA is the same letter as CB regardless of * case. * * Arguments * ========= * * CA (input) CHARACTER*1 * CB (input) CHARACTER*1 * CA and CB specify the single characters to be compared. * * ===================================================================== * * .. Intrinsic Functions .. INTRINSIC ICHAR * .. * .. Local Scalars .. INTEGER INTA, INTB, ZCODE * .. * .. Executable Statements .. * * Test if the characters are equal * LSAME = CA.EQ.CB IF( LSAME ) $ RETURN * * Now test for equivalence if both characters are alphabetic. * ZCODE = ICHAR( 'Z' ) * * Use 'Z' rather than 'A' so that ASCII can be detected on Prime * machines, on which ICHAR returns a value with bit 8 set. * ICHAR('A') on Prime machines returns 193 which is the same as * ICHAR('A') on an EBCDIC machine. * INTA = ICHAR( CA ) INTB = ICHAR( CB ) * IF( ZCODE.EQ.90 .OR. ZCODE.EQ.122 ) THEN * * ASCII is assumed - ZCODE is the ASCII code of either lower or * upper case 'Z'. * IF( INTA.GE.97 .AND. INTA.LE.122 ) INTA = INTA - 32 IF( INTB.GE.97 .AND. INTB.LE.122 ) INTB = INTB - 32 * ELSE IF( ZCODE.EQ.233 .OR. ZCODE.EQ.169 ) THEN * * EBCDIC is assumed - ZCODE is the EBCDIC code of either lower or * upper case 'Z'. * IF( INTA.GE.129 .AND. INTA.LE.137 .OR. $ INTA.GE.145 .AND. INTA.LE.153 .OR. $ INTA.GE.162 .AND. INTA.LE.169 ) INTA = INTA + 64 IF( INTB.GE.129 .AND. INTB.LE.137 .OR. $ INTB.GE.145 .AND. INTB.LE.153 .OR. $ INTB.GE.162 .AND. INTB.LE.169 ) INTB = INTB + 64 * ELSE IF( ZCODE.EQ.218 .OR. ZCODE.EQ.250 ) THEN * * ASCII is assumed, on Prime machines - ZCODE is the ASCII code * plus 128 of either lower or upper case 'Z'. * IF( INTA.GE.225 .AND. INTA.LE.250 ) INTA = INTA - 32 IF( INTB.GE.225 .AND. INTB.LE.250 ) INTB = INTB - 32 END IF LSAME = INTA.EQ.INTB * * RETURN * * End of LSAME * END SUBROUTINE XERBLA( SRNAME, INFO ) * * -- LAPACK auxiliary routine (version 2.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * September 30, 1994 * * .. Scalar Arguments .. CHARACTER*6 SRNAME INTEGER INFO * .. * * Purpose * ======= * * XERBLA is an error handler for the LAPACK routines. * It is called by an LAPACK routine if an input parameter has an * invalid value. A message is printed and execution stops. * * Installers may consider modifying the STOP statement in order to * call system-specific exception-handling facilities. * * Arguments * ========= * * SRNAME (input) CHARACTER*6 * The name of the routine which called XERBLA. * * INFO (input) INTEGER * The position of the invalid parameter in the parameter list * of the calling routine. * * ===================================================================== * * .. Executable Statements .. * WRITE( *, FMT = 9999 )SRNAME, INFO * STOP * 9999 FORMAT( ' ** On entry to ', A6, ' parameter number ', I2, ' had ', $ 'an illegal value' ) * * End of XERBLA * END SHAR_EOF fi # end of overwriting check if test -f 'src.f' then echo shar: will not over-write existing file "'src.f'" else cat << \SHAR_EOF > 'src.f' LOGICAL FUNCTION CKCOLR(MP, NP, DP, P, LDP1, LDP2, ZERCOL, Q, LDQ, * W, TOL, IERR) C C PURPOSE C C To check whether the polynomial matrix C dp-1 dp C P(s) = P(0) + P(1) * s + . . . + P(dp-1) * s + P(dp) * s , C C is column reduced. C C ARGUMENTS IN C C MP - INTEGER. C The number of rows of the polynomial matrix P(s). C MP >= 1. C NP - INTEGER. C The number of columns of the polynomial matrix P(s). C NP >= 1. C DP - INTEGER. C The degree of the polynomial matrix P(s). C DP >= 0. C P - DOUBLE PRECISION array of DIMENSION (LDP1,LDP2,DP+1). C The leading MP by NP by (DP+1) part of this array must C contain the coefficients of the polynomial matrix P(s). C Specifically, P(i,j,k) must contain the coefficient of C s**(k-1) of the polynomial which is the (i,j)-th element C of P(s), where i = 1,2,...,MP, j = 1,2,...,NP and C k = 1,2,...,DP+1. C LDP1 - INTEGER. C The leading dimension of array P as declared in the calling C program. C LDP1 >= MP. C LDP2 - INTEGER. C The second dimension of array P as declared in the calling C program. C LDP2 >= NP. C C ARGUMENTS OUT C C ZERCOL - LOGICAL array of DIMENSION at least (NP). C If ZERCOL(j) = .TRUE. then the j-th column of P(s) is zero; C otherwise the j-th column belongs to P1(s) (see METHOD). C C WORKSPACE C C Q - DOUBLE PRECISION array of DIMENSION (LDQ,NP). C LDQ - INTEGER. C The leading dimension of array Q as declared by the calling C program. C LDQ >= MP. C W - DOUBLE PRECISION array of DIMENSION (MP). C C TOLERANCES C C TOL - DOUBLE PRECISION. C A tolerance below which matrix elements are considered to C be zero. If the user sets TOL to be less than C EPS * (((DP+1)*MP)**2 * MAX(P(i,j,k))), then the tolerance C is taken as EPS * (((DP+1)*MP)**2 * MAX(P(i,j,k))), C i = 1,...,MP, j = 1,...,NP, k = 1,...,DP+1, where EPS is C the machine precision. C C ERROR INDICATOR C C IERR - INTEGER. C Unless the routine detects an error (see next section), C IERR contains 0 on exit. C C WARNINGS AND ERRORS DETECTED BY THE ROUTINE C C IERR = 1 : Invalid input parameter(s). C C METHOD C C Let GAMC(P) be the contant matrix such that each of its columns C contains the coefficients of the highest power of s occurring C in the corresponding column of P(s), the so-called leading column C coefficient matrix. Then P(s) is called column reduced if there C exists a permutation matrix T such that P(s) = ( Z , P1(s) ) * T, C where Z is a zero matrix and GAMC(P1) has full column rank. C C The algorithm used, which is in fact the QR decomposition of the C leading column coefficient matrix, is as follows: C The columns of the leading column coefficient matrix of P1(s) are C determined one by one, where a column is considered zero if its C Euclidean norm is less than TOL. To each new column the C Householder transformations are applied that have transformed the C submatrix of the former columns in upper triangular form. If the C new column is independent of its predecessors, then a new C Householder transformation is generated and applied such that the C augmented matrix is upper triangular. C The routine terminates after the last column of P(s) has been C treated or when the new column is dependent of its predecessors. C C CONTRIBUTOR C C A.J. Geurts (Eindhoven University of Technology). C C. Praagman (University of Groningen). C C REVISIONS C C 1994, February 11. C 1996, May 28. C C IMPLICIT NONE C .. Parameters .. DOUBLE PRECISION ZERO PARAMETER (ZERO = 0.0D0) C .. Scalar Arguments .. INTEGER MP, NP, DP, LDP1, LDP2, LDQ, IERR DOUBLE PRECISION TOL C .. Array Arguments .. DOUBLE PRECISION P(LDP1,LDP2,*), Q(LDQ,*), W(*) LOGICAL ZERCOL(*) C .. Local Scalars .. INTEGER H, J, J1, K DOUBLE PRECISION EPS, TOLER, NORM, TAU LOGICAL FULLRK, NOTCJ C .. External Subroutines/Functions .. EXTERNAL DCOPY, DLAMCH, DLANGE, DLARFG, DLARFV, DLASET, DNRM2 DOUBLE PRECISION DLAMCH, DLANGE, DNRM2 C .. Intrinsic Functions .. INTRINSIC DBLE, MAX, MIN C C .. Executable Statements .. C C Check input parameters C IF (MP.LT.1 .OR. NP.LT.1 .OR. DP.LT.0 .OR. LDP1.LT.MP .OR. * LDP2.LT.NP) THEN IERR = 1 RETURN END IF C C Computation of the tolerance. EPS is the machine precision. C EPS = DLAMCH('E') TOLER = ZERO DO 10 K = 1, DP + 1 TOLER = MAX(TOLER, DLANGE('M', MP, NP, P(1,1,K), LDP1, W)) 10 CONTINUE TOLER = DBLE(((DP+1) * MP))**2 * TOLER * EPS IF (TOLER .LT. TOL) TOLER = TOL C CALL DLASET('G', MP, NP, ZERO, ZERO, Q, LDQ) J = 1 J1 = 1 FULLRK = .TRUE. C WHILE (FULLRK and J1 <= NP) DO 20 IF (FULLRK .AND. J1.LE.NP) THEN C C Find the j-th column of the leading column coefficient matrix C of P1(s) and put it in W. C K = DP + 1 NOTCJ = .TRUE. C WHILE (j-th column not found) DO 30 IF (NOTCJ .AND. K.GE.1) THEN NORM = DNRM2(MP, P(1,J1,K), 1) IF (NORM .GE. TOLER) THEN CALL DCOPY(MP, P(1,J1,K), 1, W, 1) NOTCJ = .FALSE. END IF K = K - 1 GO TO 30 END IF C END WHILE 30 C C Check whether the j-th column is linearly independent of the C preceeding columns. C IF (NOTCJ) THEN ZERCOL(J1) = .TRUE. J1 = J1 + 1 ELSE ZERCOL(J1) = .FALSE. C C Apply the Householder transformations Qh, C h = 1,...,min(mp,j) - 1, to W. C DO 40 H = 1, MIN(MP,J) - 1 CALL DLARFV(MP-H+1, Q(H+1,H), 1, Q(H,H), W(H), W(H+1), 1) 40 CONTINUE NORM = DNRM2(MP-J+1, W(J), 1) IF (NORM .LT. TOLER) THEN FULLRK = .FALSE. ELSE C C Generate the Householder transformation Qj. C IF (J .LT. MP) THEN CALL DLARFG(MP-J+1, W(J), W(J+1), 1, TAU) CALL DCOPY(MP-J, W(J+1), 1, Q(J+1,J), 1) Q(J,J) = TAU ELSE Q(J,J) = ZERO END IF END IF J = J + 1 J1 = J1 + 1 END IF GO TO 20 END IF C END WHILE 20 CKCOLR = FULLRK RETURN C *** Last line of CKCOLR *** END C LOGICAL FUNCTION CKGAMC(MP, INV, GAMC, LDG, Q, LDQ, RWORK, TOL) C C PURPOSE C C To check whether the leading coefficient matrix has still full C column rank after a new column has been appended. C REMARK: This auxiliary routine is intended to be called only from C the routine COLRED. C C ARGUMENTS IN C C MP - INTEGER. C The number of rows of the matrix GAMC. C MP >= 1. C INV - INTEGER. C The number of columns of the matrix GAMC. C INV >= 1. C GAMC - DOUBLE PRECISION array of DIMENSION (LDG,INV) C The leading MP by INV part of this array must contain the C leading coefficient matrix GAMC of which the first INV - 1 C columns have been transformed in upper triangular form. C Note: this array is overwritten C LDG - INTEGER. C The leading dimension of array GAMC as declared in the C calling program. C LDG >= MP. C Q - DOUBLE PRECISION array of DIMENSION (LDQ,INV) C The leading MP by INV - 1 part of this array must contain C the constants and the vectors of the elementary Householder C transformations, as supplied by DLARFG, by which the first C INV - 1 columns of GAMC have been transformed into an upper C triangular matrix. C LDQ - INTEGER. C The leading dimension of array Q as declared in the calling C program. C LDQ >= MP. C C ARGUMENTS OUT C C GAMC - DOUBLE PRECISION array of DIMENSION (LDG,INV) C The leading MP by INV part of this array contains the C leading coefficient matrix GAMC transformed in upper C triangular form. C Q - DOUBLE PRECISION array of DIMENSION (LDQ,INV) C The leading MP by INV part of this array contains the C constants and the vectors of the elementary Householder C transformations by which GAMC has been transformed into an C upper triangular matrix. C C WORKSPACE C C RWORK - DOUBLE PRECISION array of DIMENSION (MP) C C TOLERANCES C C TOL - DOUBLE PRECISION. C A tolerance below which matrix elements are considered to C be zero. C C METHOD C C Let the first INV - 1 columns of GAMC be linearly independent, C which has been checked by former calls of CKGAMC. A new column is C appended. To this column the Householder transformations are C applied that have transformed the matrix of the former columns in C upper triangular form. If the new column is independent of its C predecessors, then a new Householder transformation is generated C and applied such that the augmented matrix is upper triangular. C C CONTRIBUTOR C C A.J. Geurts. C C REVISIONS C C 1993, October 29. C 1994, December 8. C C IMPLICIT NONE C .. Parameters .. DOUBLE PRECISION ZERO PARAMETER (ZERO = 0.0D0) C .. Scalar Arguments .. INTEGER MP, INV, LDG, LDQ DOUBLE PRECISION TOL C .. Array Arguments .. DOUBLE PRECISION GAMC(LDG,*), Q(LDQ,*), RWORK(*) C .. Local Scalars .. INTEGER J DOUBLE PRECISION NORM, TAU C .. External Subroutines/Functions .. EXTERNAL DCOPY, DLARFG, DLARFV, DLAVEC, DNRM2 DOUBLE PRECISION DNRM2 C C .. Executable Statements .. C IF (INV .LE. MP) THEN CALL DLAVEC(MP, ZERO, Q(1,INV), 1) C C Check whether the INV-th column is linearly independent of the C preceeding columns by applying the Householder transformations C Q(h), h = 1,...,INV-1, to the INV-th column. C CALL DCOPY(MP, GAMC(1,INV), 1, RWORK, 1) DO 10 J = 1, INV - 1 CALL DLARFV(MP-J+1, Q(J+1,J), 1, Q(J,J), RWORK(J), * RWORK(J+1), 1) 10 CONTINUE NORM = DNRM2(MP-INV+1, RWORK(INV), 1) IF (NORM .LT. TOL) THEN CKGAMC = .FALSE. ELSE C C Generate the Householder transformation Q(INV) if INV < MP. C CALL DCOPY(MP, RWORK, 1, GAMC(1,INV), 1) IF (INV .LT. MP) THEN CALL DLARFG(MP-INV+1, RWORK(INV), RWORK(INV+1), 1, TAU) GAMC(INV,INV) = RWORK(INV) CALL DLAVEC(MP-INV, ZERO, GAMC(INV+1,INV), 1) CALL DCOPY(MP-INV, RWORK(INV+1), 1, Q(INV+1,INV), 1) Q(INV,INV) = TAU ELSE Q(INV,INV) = ZERO END IF CKGAMC = .TRUE. END IF ELSE CKGAMC = .FALSE. END IF RETURN C *** Last line of CKGAMC *** END C SUBROUTINE COLRD1(MP, NP, DP, P, LDP1, LDP2, DR, DU, R, LDR1, * LDR2, U, LDU1, LDU2, ZERCOL, MU, S, SK, A, LDA, * AB, LDAB, Q, LDQ, Y, LDY, YI, GAMC, LDG, W, * TOL, IERR) C C PURPOSE C C To compute for a given polynomial matrix C dp-1 dp C P(s) = P(0) + P(1) * s + . . . + P(dp-1) * s + P(dp) * s , C C which is not column reduced, a unimodular polynomial matrix U(s) C such that R(s) = P(s) * U(s) is column reduced. C REMARK: This auxiliary routine is intended to be called only from C the routine COLRED. C C ARGUMENTS IN C C MP - INTEGER. C The number of rows of the polynomial matrix P(s). C MP >= 1. C NP - INTEGER. C The number of columns of the polynomial matrix P(s). C NP >= 1. C DP - INTEGER. C The degree of the polynomial matrix P(s). C DP >= 0. C P - DOUBLE PRECISION array of DIMENSION (LDP1,LDP2,DP+1). C The leading MP by NP by (DP+1) part of this array must C contain the coefficients of the polynomial matrix P(s). C Specifically, P(i,j,k) must contain the coefficient of C s**(k-1) of the polynomial which is the (i,j)-th element C of P(s), where i = 1,2,...,MP, j = 1,2,...,NP and C k = 1,2,...,DP+1. C LDP1 - INTEGER. C The leading dimension of array P as declared in the calling C program. C LDP1 >= MP. C LDP2 - INTEGER. C The second dimension of array P as declared in the calling C program. C LDP2 >= NP. C C ARGUMENTS OUT C C DR - INTEGER. C The degree of the column reduced polynomial matrix R(s). C DU - INTEGER. C The degree of the unimodular polynomial matrix U(s). C R - DOUBLE PRECISION array of DIMENSION (LDR1,LDR2,DP+1). C The leading MP by NP by (DR+1) part of this array contains C the coefficients of the column reduced polynomial matrix C R(s). Specifically, R(i,j,k) contains the coefficient of C s**(k-1) of the polynomial which is the (i,j)-th element of C R(s), where i = 1,2,...,MP, j = 1,2,...,NP and C k = 1,2,...,DR+1. C LDR1 - INTEGER. C The leading dimension of array R as declared in the calling C program. C LDR1 >= MP. C LDR2 - INTEGER. C The second dimension of array R as declared in the calling C program. C LDR2 >= NP. C U - DOUBLE PRECISION array of DIMENSION (LDU1,LDU2,NP*DP+1). C The leading NP by NP by (DU+1) part of this array contains C the coefficients of the unimodular polynomial matrix U(s). C Specifically, U(i,j,k) contains the coefficient of s**(k-1) C of the polynomial which is the (i,j)-th element of U(s), C where i = 1,2,...,NP, j = 1,2,...,NP and k = 1,2,...,DU+1. C LDU1 - INTEGER. C The leading dimension of array U as declared in the calling C program. C LDU1 >= NP. C LDU2 - INTEGER. C The second dimension of array U as declared in the calling C program. C LDU2 >= NP. C ZERCOL - LOGICAL array of DIMENSION at least (NP). C If ZERCOL(j) = .TRUE. then the j-th column of R(s) is zero; C otherwise the j-th column belongs to R1(s) (see METHOD). C C WORKSPACE C C MU - INTEGER array of DIMENSION at least (mamax+1), C where mamax = (NP * DP + 1) * MP. C On exit, this array contains the row indices of the left C upper elements of the right invertible diagonal submatrices C of A'. C S - INTEGER array of DIMENSION at least (mamax). C On exit, this array contains the column indices of the C pivots of A', that is A transformed in upper staircase C form. C SK - INTEGER array of DIMENSION at least (2*NP). C A - DOUBLE PRECISION array of DIMENSION (LDA,namax), C where namax = mamax + NP. C On exit, the upper block diagonal part of this array C contains the transformed matrix A', the lower part contains C the vectors of the Householder transformations. C LDA - INTEGER. C The leading dimension of array A as declared in the calling C program. C LDA >= mamax + 1. C AB - DOUBLE PRECISION array of DIMENSION (LDAB,namax). C Array in which the transformed matrix Ab is saved. C LDAB - INTEGER. C The leading dimension of array AB as declared in the C calling program. C LDAB >= mamax + 1. C Q - DOUBLE PRECISION array of DIMENSION (LDQ,NP). C LDQ - INTEGER. C The leading dimension of array Q as declared in the calling C program. C LDQ >= MP. C Y - DOUBLE PRECISION array of DIMENSION (LDY,NP*DP+2). C Array in which a null vector of sA - E or sA' - E is C stored. C LDY - INTEGER. C The leading dimension of array Y as declared in the calling C program. C LDA >= namax. C YI - DOUBLE PRECISION array of DIMENSION at least (NP). C GAMC - DOUBLE PRECISION array of DIMENSION (LDG,NP) C On exit, this array contains the leading column coefficient C matrix of R(s), which has full column rank. C LDG - INTEGER. C The leading dimension of array GAMC as declared in the C calling program. C LDG >= MP. C W - DOUBLE PRECISION array of DIMENSION (mamax+namax). C C TOLERANCES C C TOL - DOUBLE PRECISION. C A tolerance below which matrix elements are considered to C be zero. C C ERROR INDICATOR C C IERR - INTEGER. C Unless the routine detects an error (see next section), C IERR contains 0 on exit. C C WARNINGS AND ERRORS DETECTED BY THE ROUTINE C C IERR = 2 : No column reduced R(s) has been found for the C maximum b, (see METHOD). C IERR = 3 : The computation of a null vector has failed, C because a diagonal block of A' is not right C invertible. C IERR = 4 : The computation of R(s) has failed, because a C computed null vector is not s**B times a polynomial C vector. C C METHOD C C Let GAMC(P) be the constant matrix such that each of its columns C contains the coefficients of the highest power of s occurring in C the corresponding column of P(s), the so-called leading column C coefficient matrix. Then P(s) is called column reduced if there C exists a permutation matrix T such that P(s) = ( Z , P1(s) ) * T, C where Z is a zero matrix and GAMC(P1) has full column rank. C C Let (U(s),Z(s))' be a minimal polynomial basis (MPB) for C b C Ker(s P(s), -I), for some b > 0. It has been proved, see [1], that C if b is greater than d'c(P), the sum of all but the smallest C -b C column degrees of P(s), then U(s) is unimodular and R(s) = s Z(s) C is column reduced and P(s) * U(s) = R(s). C The routine computes an MPB for b = 1,2,... and checks for each b C whether R(s) is column reduced, i.e. whether GAMC(R1) has full C column rank. The algorithm finishes with U(s) and R(s) as soon as C R(s) is column reduced. C C REFERENCES C C [1] Neven, W.H.L. and Praagman, C. C Column Reduction of Polynomial Matrices. C Linear Algebra and its Applications 188, 189, pp. 569-589, C 1993. C C CONTRIBUTORS C C A.J. Geurts (Eindhoven University of Technology). C C. Praagman (University of Groningen). C C REVISIONS C C 1993, November 8. C 1996, May 21. C C IMPLICIT NONE C .. Parameters .. DOUBLE PRECISION ZERO, ONE PARAMETER (ZERO = 0.0D0, ONE = 1.0D0) C .. Scalar Arguments .. INTEGER MP, NP, DP, LDP1, LDP2, DR, DU, LDR1, LDR2, LDU1, LDU2, * LDA, LDAB, LDQ, LDY, LDG, IERR DOUBLE PRECISION TOL C .. Array Arguments .. INTEGER MU(*), S(*), SK(*) DOUBLE PRECISION P(LDP1,LDP2,*), R(LDR1,LDR2,*), U(LDU1,LDU2,*), * A(LDA,*), AB(LDAB,*), Q(LDQ,*), Y(LDY,*), YI(*), * GAMC(LDG,*), W(*) LOGICAL ZERCOL(*) C .. Local Scalars .. INTEGER BMAX, B, I, IC, IR, IND, NNV, NNVB, NCR1, J, K, KK, * L, C1AI, DUB, DUR1, MA, NA, MAI, NAI, MUI DOUBLE PRECISION NORM, TAU LOGICAL COLRDC C .. External Subroutines/Functions .. EXTERNAL CKGAMC, COMPTV, COMPTY, COMPYI, DCOPY, DGEMV, DGER, * DLACPY, DLAVEC, DLASET, DNRM2, HHTRAN, MKPENC DOUBLE PRECISION DNRM2 LOGICAL CKGAMC C .. Intrinsic Functions .. INTRINSIC MAX C C .. Executable Statements .. C BMAX = (NP-1) * DP + 1 CALL MKPENC(MP, NP, DP, P, LDP1, LDP2, AB, LDAB, MA, NA) NNVB = 0 DUB = -1 DO 10 K = 1, NP * DP + 1 CALL DLASET('G', NP, NP, ZERO, ZERO, U(1,1,K), LDU1) 10 CONTINUE COLRDC = .FALSE. MU(1) = 1 B = 0 C WHILE (.NOT.COLRDC and B < BMAX) DO 20 IF (.NOT.COLRDC .AND. B.LT.BMAX) THEN B = B + 1 C C Initialization of A. C IF (B .EQ. 1) THEN CALL DLACPY('G', MA, NA, AB, LDAB, A, LDA) CALL DLAVEC(NA, ZERO, A(MA+1,1), LDA) ELSE CALL DLACPY('G', MA+1, NA, AB, LDAB, A, LDA) J = NP + MU(B-1) L = MA - MP - MU(B-1) + 1 CALL DLASET('G', MP, J-1, ZERO, ZERO, A(MA+2,1), LDA) CALL DLASET('G', MP+1, L, ZERO, ZERO, A(MA+1,J), LDA) CALL DLASET('G', MP+1, MP, ZERO, ONE, A(MA+1,L+J), LDA) MA = MA + MP CALL DLASET('G', MA+1, MP, ZERO, ZERO, A(1,NA+1), LDA) NA = NA + MP DO 30 K = MU(B-1), MU(B) - 1 TAU = A(K+1,S(K)) IF (TAU .NE. ZERO) THEN J = NP + K L = MA - MP - K + 1 W(MA+1) = ONE CALL DCOPY(L-1, A(K+2,S(K)), 1, W(MA+2), 1) CALL DGEMV('N', MP, L, ONE, A(MA-MP+1,J), LDA, * W(MA+1), 1, ZERO, W, 1) CALL DGER(MP, L, -TAU, W, 1, W(MA+1), 1, * A(MA-MP+1,J), LDA) END IF 30 CONTINUE END IF C I = B DU = DUB NNV = NNVB DR = -1 DO 40 K = 1, DU + 1 CALL DLASET('G', NP, NP-NNV, ZERO, ZERO, U(1,NNV+1,K), * LDU1) 40 CONTINUE DO 50 K = 1, DP + 1 CALL DLASET('G', MP, NP, ZERO, ZERO, R(1,1,K), LDR1) 50 CONTINUE NCR1 = 0 COLRDC = .TRUE. C C NNV is the number of already found null vectors. With these C null vectors correspond the MP by NNV column reduced matrix C R(s). NCR1 is the number of columns in R1(s). C C WHILE (NNV < NP, which means that not all null vectors of C sA - E have been found, and R(s) is column reduced) DO 60 IF ((NNV.LT.NP) .AND. COLRDC) THEN C C Determine A(i,i), the i(th) right invertible diagonal block. C The left upper element of A(i,i) is A(MUI,C1AI) and the row C and column dimensions are MAI and NAI, respectively. C MUI = MU(I) MAI = 0 IF (I .EQ. 1) THEN C1AI = 1 NAI = NP ELSE C1AI = MU(I-1) + NP NAI = MU(I) - MU(I-1) END IF C C Compute the null vectors of A(i,i) and the corresponding C nullvectors of sA - E one by one, append the appropriate C parts of a computed null vector to U(s) and R(s) and check C whether R(s) remains column reduced. C K = 1 C WHILE (K <= NAI and R(s) still column reduced) DO 70 IF (K.LE.NAI .AND. COLRDC) THEN IR = MUI + MAI IC = C1AI + K - 1 L = MA - IR + 1 IF (L .GT. 0) THEN NORM = DNRM2(L, A(IR,IC), 1) ELSE NORM = ZERO END IF IF (NORM .GE. TOL) THEN C C Generate and apply the Householder transformation for C the IC-th column of A. C CALL HHTRAN(MA, NA, IC, IR, A, LDA, W, W(MA+1)) CALL DCOPY(L, W, 1, A(IR+1,IC), 1) S(MUI+MAI) = IC MAI = MAI + 1 SK(MAI) = IC - C1AI + 1 ELSE C C Compute the null vector Y of sA - E corresponding C to the IC-th column of A. C IF (L .GT. 0) THEN CALL DLAVEC(L, ZERO, A(IR,IC), 1) END IF IF (MAI .GT. 0) THEN CALL DCOPY(MAI, A(MUI,IC), 1, W, 1) CALL COMPYI(MAI, K-1, A(MUI,C1AI), LDA, SK, W, * YI, IERR) IF (IERR .NE. 0) RETURN ELSE CALL DLAVEC(K-1, ZERO, YI, 1) END IF YI(K) = -ONE CALL COMPTY(NP, NA, I, K, YI, MU, A, LDA, S, Y, LDY, * SK(NP+1), W, IERR) IF (IERR .NE. 0) THEN IERR = 3 RETURN END IF CALL COMPTV(MA, NA, I, IR, A, LDA, S, Y, LDY, W) C C Append the first NP by I block of Y to the unimodular C U, and the last MP by I block to the column reduced R C and the last MP elements of the I-th column of Y to C the leading column coefficient matrix GAMC. C NNV = NNV + 1 DUR1 = 0 DO 80 KK = 1, I NORM = DNRM2(NP, Y(1,KK), 1) IF (NORM .GE. TOL) THEN DUR1 = KK - 1 CALL DCOPY(NP, Y(1,KK), 1, U(1,NNV,KK), 1) ELSE CALL DLAVEC(NP, ZERO, U(1,NNV,KK), 1) END IF 80 CONTINUE DU = MAX(DU, DUR1) IND = NA - MP + 1 DO 90 KK = 1, B NORM = DNRM2(MP, Y(IND,KK), 1) IF (NORM .NE. ZERO) THEN IERR = 4 RETURN END IF 90 CONTINUE DUR1 = -1 DO 100 KK = B + 1, I NORM = DNRM2(MP, Y(IND,KK), 1) IF (NORM .GE. TOL) THEN DUR1 = KK - B - 1 CALL DCOPY(MP, Y(IND,KK), 1, R(1,NNV,KK-B), 1) ELSE CALL DLAVEC(MP, ZERO, R(1,NNV,KK-B), 1) END IF 100 CONTINUE IF (DUR1 .NE. -1) THEN NCR1 = NCR1 + 1 ZERCOL(NNV) = .FALSE. DR = MAX(DR, DUR1) CALL DCOPY(MP, Y(IND,DUR1+B+1), 1, GAMC(1,NCR1), 1) C C Check whether R is still column reduced. C COLRDC = CKGAMC(MP, NCR1, GAMC, LDG, Q, LDQ, W, * TOL) ELSE ZERCOL(NNV) = .TRUE. END IF END IF K = K + 1 GO TO 70 END IF C END WHILE 70 C C Save the transformed A if I = B. C IF (I .EQ. B) THEN CALL DLACPY('G', MA+1, NA, A, LDA, AB, LDAB) NNVB = NNV DUB = DU END IF C IF (NNV.LE.NP .AND. COLRDC) THEN C C If i > b + dp + 1, then the degree of the next computed C column in R(s) will be greater than dp. This will occur C if a null vector is not found, due to rounding errors. C COLRDC = I .LE. (B + DP + 1) I = I + 1 MU(I) = MU(I-1) + MAI END IF GO TO 60 END IF C END WHILE 60 GO TO 20 END IF C END WHILE 20 C IF (.NOT. COLRDC) IERR = 2 RETURN C *** Last line of COLRD1 *** END C SUBROUTINE COLRED(MP, NP, DP, P, LDP1, LDP2, DR, DU, R, LDR1, * LDR2, U, LDU1, LDU2, ZERCOL, IWORK, RWORK, * TOL, IERR) C C PURPOSE C C To compute for a given polynomial matrix C dp-1 dp C P(s) = P(0) + P(1) * s + . . . + P(dp-1) * s + P(dp) * s , C C a unimodular polynomial matrix U(s) such that R(s) = P(s) * U(s) C is column reduced. C C ARGUMENTS IN C C MP - INTEGER. C The number of rows of the polynomial matrix P(s). C MP >= 1. C NP - INTEGER. C The number of columns of the polynomial matrix P(s). C NP >= 1. C DP - INTEGER. C The degree of the polynomial matrix P(s). C DP >= 0. C P - DOUBLE PRECISION array of DIMENSION (LDP1,LDP2,DP+1). C The leading MP by NP by (DP+1) part of this array must C contain the coefficients of the polynomial matrix P(s). C Specifically, P(i,j,k) must contain the coefficient of C s**(k-1) of the polynomial which is the (i,j)-th element C of P(s), where i = 1,2,...,MP, j = 1,2,...,NP and C k = 1,2,...,DP+1. C LDP1 - INTEGER. C The leading dimension of array P as declared in the calling C program. C LDP1 >= MP. C LDP2 - INTEGER. C The second dimension of array P as declared in the calling C program. C LDP2 >= NP. C C ARGUMENTS OUT C C DR - INTEGER. C The degree of the column reduced polynomial matrix R(s). C DU - INTEGER. C The degree of the unimodular polynomial matrix U(s). C R - DOUBLE PRECISION array of DIMENSION (LDR1,LDR2,DP+1). C The leading MP by NP by (DR+1) part of this array contains C the coefficients of the column reduced polynomial matrix C R(s). Specifically, R(i,j,k) contains the coefficient of C s**(k-1) of the polynomial which is the (i,j)-th element of C R(s), where i = 1,2,...,MP, j = 1,2,...,NP and C k = 1,2,...,DR+1. C LDR1 - INTEGER. C The leading dimension of array R as declared in the calling C program. C LDR1 >= MP. C LDR2 - INTEGER. C The second dimension of array R as declared in the calling C program. C LDR2 >= NP. C U - DOUBLE PRECISION array of DIMENSION (LDU1,LDU2,NP*DP+1). C The leading NP by NP by (DU+1) part of this array contains C the coefficients of the unimodular polynomial matrix U(s). C Specifically, U(i,j,k) contains the coefficient of s**(k-1) C of the polynomial which is the (i,j)-th element of U(s), C where i = 1,2,...,NP, j = 1,2,...,NP and k = 1,2,...,DU+1. C LDU1 - INTEGER. C The leading dimension of array U as declared in the calling C program. C LDU1 >= NP. C LDU2 - INTEGER. C The second dimension of array U as declared in the calling C program. C LDU2 >= NP. C ZERCOL - LOGICAL array of DIMENSION at least (NP). C If ZERCOL(j) = .TRUE. then the j-th column of R(s) is zero; C otherwise the j-th column belongs to R1(s) (see METHOD). C C WORKSPACE C C IWORK - INTEGER array of DIMENSION at least (liwork), C where liwork = 2*mamax + 2 * NP + 1, C and mamax = (NP*DP+1) * MP. C RWORK - DOUBLE PRECISION array of DIMENSION at least (lrwork), C where lrwork = (2*mamax+NP*DP+5) * namax + mamax + (2*MP+1) * NP, C and namax = mamax + NP. C . C TOLERANCES C C TOL - DOUBLE PRECISION. C A tolerance below which matrix elements are considered to C be zero. If the user sets TOL to be less than C EPS * (((DP+1)*MP)**2 * MAX(P(i,j,k))), then the tolerance C is taken as EPS * (((DP+1)*MP)**2 * MAX(P(i,j,k))), C i = 1,...,MP, j = 1,...,NP, k = 1,..., DP + 1, where EPS is C the machine precision. C C ERROR INDICATOR C C IERR - INTEGER. C Unless the routine detects an error (see next section), C IERR contains 0 on exit. C C WARNINGS AND ERRORS DETECTED BY THE ROUTINE C C IERR = 1 : On entry, MP < 1 or NP < 1 or DP < 0 or C LDP1 < MP or LDP2 < NP or C LDR1 < MP or LDR2 < NP or C LDU1 < NP or LDU2 < NP. C IERR = 2 : P(s) is not column reduced, and no column reduced C R(s) has been found. C C METHOD C C Let GAMC(P) be the constant matrix such that each of its columns C contains the coefficients of the highest power of s occurring in C the corresponding column of P(s), the so-called leading column C coefficient matrix. Then P(s) is called column reduced if there C exists a permutation matrix T such that P(s) = ( Z , P1(s) ) * T, C where Z is a zero matrix and GAMC(P1) has full column rank. C C Let (U(s),Z(s))' be a minimal polynomial basis (MPB) for C b C Ker(s P(s), -I), for some b > 0. It has been proved, see [1], that C if b is greater than d'c(P), the sum of all but the smallest C -b C column degrees of P(s), then U(s) is unimodular and R(s) = s Z(s) C is column reduced and P(s) * U(s) = R(s). C b C The routine uses a linearization of (s P(s), -I) to compute an MPB C for b = 1,2,... and checks for each b whether R(s) is column C reduced, i.e. whether GAMC(R1) has full column rank. The algorithm C finishes with U(s) and R(s) as soon as R(s) is column reduced. C C REFERENCES C C [1] Neven, W.H.L. and Praagman, C. C Column Reduction of Polynomial Matrices. C Linear Algebra and its Applications 188, 189, pp. 569-589, C 1993. C [2] Geurts, A.J. and Praagman, C. C A Fortran subroutine for column reduction of polynomial C matrices. EUT Report 94-WSK-01, Eindhoven, June 1994. C C NUMERICAL ASPECTS C C The algorithm used by the routine involves the construction of a C b C special staircase form of a linearization of (s P(s), -I) with C pivots considered to be non-zero when they are greater than or C equal to TOL. These pivots are then inverted in order to construct C b C the columns of ker(s P(s), -I). C The user is recommended to choose TOL of the order of the relative C error in the elements of P(s). If TOL is chosen to be too small, C then a very small element of insignificant value may be taken as C pivot. As a consequence, the correct null-vectors, and hence R(s), C may not be found. In the case that R(s) has not been found and in C the case that the elements of the computed U(s) and R(s) are large C relative to the elements of P(s) the user should consider trying C several values of TOL. C C CONTRIBUTORS C C A.J. Geurts (Eindhoven University of Technology). C C. Praagman (University of Groningen). C C REVISIONS C C 1994, February 11. C 1996, June 5. C C IMPLICIT NONE C .. Parameters .. DOUBLE PRECISION ZERO, ONE PARAMETER (ZERO = 0.0D0, ONE = 1.0D0) C .. Scalar Arguments .. INTEGER MP, NP, DP, LDP1, LDP2, DR, DU, LDR1, LDR2, LDU1, LDU2, * IERR DOUBLE PRECISION TOL C .. Array Arguments .. INTEGER IWORK(*) DOUBLE PRECISION P(LDP1,LDP2,*), R(LDR1,LDR2,*), U(LDU1,LDU2,*), * RWORK(*) LOGICAL ZERCOL(*) C .. Local Scalars .. INTEGER LDA, LDAB, LDQ, LDY, LDG, MU, S, SK, A, AB, Q, Y, YI, * GAMC, BMAX, DP1, K, MAMAX, NAMAX DOUBLE PRECISION EPS, TOLER LOGICAL COLRDC, PKZERO C .. External Subroutines/Functions .. EXTERNAL COLRD1, CKCOLR, DLACPY, DLAMCH, DLANGE, DLASET DOUBLE PRECISION DLAMCH, DLANGE LOGICAL CKCOLR C .. Intrinsic Functions .. INTRINSIC DBLE, MAX C C .. Executable Statements .. C C Check the input parameters. C IF (MP.LT.1 .OR. NP.LT.1 .OR. DP.LT.0 .OR. LDP1.LT.MP .OR. * LDP2.LT.NP .OR. LDR1.LT.MP .OR. LDR2.LT.NP .OR. * LDU1.LT.NP .OR. LDU2.LT.NP) THEN IERR = 1 RETURN END IF C IERR = 0 C C Computation of the tolerance. EPS is the machine precision. C EPS = DLAMCH('E') TOLER = ZERO DO 10 K = 1, DP + 1 TOLER = MAX(TOLER, DLANGE('M', MP, NP, P(1,1,K), LDP1, RWORK)) 10 CONTINUE C TOLER = DBLE(((DP+1) * MP)**2) * TOLER * EPS IF (TOLER .LT. TOL) TOLER = TOL C C Computation of the true degree of P(s). C K = DP + 2 PKZERO = .TRUE. C WHILE (P(k) is a zero matrix) DO 20 IF (PKZERO .AND. (K.GT.1)) THEN K = K - 1 PKZERO = (DLANGE('M', MP, NP, P(1,1,K), LDP1, RWORK) .EQ. ZERO) GO TO 20 END IF C END WHILE 20 DP1 = K - 1 C C Check whether P(s) is already column reduced. C Q = MP + 1 LDQ = MP COLRDC = CKCOLR(MP, NP, DP1, P, LDP1, LDP2, ZERCOL, RWORK(Q), LDQ, * RWORK, TOLER, IERR) IF (COLRDC) THEN DR = DP1 DO 30 K = 1, DR + 1 CALL DLACPY('G', MP, NP, P(1,1,K), LDP1, R(1,1,K), LDR1) 30 CONTINUE DU = 0 CALL DLASET('G', NP, NP, ZERO, ONE, U, LDU1) RETURN END IF C BMAX = (NP - 1) * DP1 + 1 MAMAX = (DP1 + BMAX) * MP NAMAX = MAMAX + NP LDA = MAMAX + 1 LDAB = LDA LDY = NAMAX LDG = MP MU = 1 S = MU + MAMAX + 1 SK = S + MAMAX A = 1 + MAMAX + NAMAX AB = A + (MAMAX + 1) * NAMAX Q = AB + (MAMAX + 1) * NAMAX Y = Q + MP * NP YI = Y + (BMAX + DP1 + 1) * NAMAX GAMC = YI + NP C CALL COLRD1(MP, NP, DP1, P, LDP1, LDP2, DR, DU, R, LDR1, LDR2, * U, LDU1, LDU2, ZERCOL, IWORK(MU), IWORK(S), IWORK(SK), * RWORK(A), LDA, RWORK(AB), LDAB, RWORK(Q), LDQ, * RWORK(Y), LDY, RWORK(YI), RWORK(GAMC), LDG, RWORK, * TOLER, IERR) C C Check whether the computed R(s) is column reduced. C IF (IERR .EQ. 0) THEN Q = MP + 1 COLRDC = CKCOLR(MP, NP, DR, R, LDR1, LDR2, ZERCOL, RWORK(Q), * LDQ, RWORK, TOL, IERR) IF (.NOT. COLRDC) IERR = 2 ELSE IERR = 2 END IF RETURN C *** Last line of COLRED *** END C SUBROUTINE COMPTV(MA, NA, I, J, A, LDA, S, Y, LDY, RWORK) C C PURPOSE C C To apply the Householder reflections, thus far applied to sA - E, C to the null vector Y(s) corresponding to a given column of the C transformed sA' - E, which transforms this null vector into a null C vector V(s) of the original pencil sA - E. C REMARK: This auxiliary routine is intended to be called only from C the routine COLRED. C C ARGUMENTS IN C C MA - INTEGER. C The number of rows of matrix A. C MA >= 1. C NA - INTEGER. C The number of columns of matrix A. C NA >= MA. C I - INTEGER. C The actual number of columns in Y, which is also the index C of the diagonal block in A corresponding to the null vector C Y(s). C I >= 1. C J - INTEGER. C The row index in A which corresponds to the null vector C Y(s). C J >= 1. C A - DOUBLE PRECISION array of DIMENSION (LDA,NA). C The leading (MA + 1) by NA part of this array must contain C the transformed matrix A' in the upper part and the C Householder transformation vectors in the lower part. C LDA - INTEGER. C The leading dimension of array A as declared in the calling C program. C LDA >= MA + 1. C S - INTEGER ARRAY of DIMENSION at least (J-1). C The leading J - 1 elements of this array must contain the C indices of the pivots of the right invertible diagonal C submatrices, i.e., the pivot of A(m,m) is A(m,S(m)), C M=1,...,J-1. S(m) is also the index of the column in array C A in which the m-th non-trivial Householder transformation C vector is stored. C Y - DOUBLE PRECISION array of DIMENSION (LDY,I). C The leading NA by I part of this array must contain the C polynomial null vector Y(s) of sA' - E to be transformed, C where the t-th column (t = 1,...,I) must contain the C coefficient of s**(t-1). C Note: this array is overwritten. C LDY - INTEGER. c The leading dimension of array Y as declared in the calling C program. C LDY >= NA. C C ARGUMENTS OUT C C Y - DOUBLE PRECISION array of DIMENSION (LDY,I). C The leading NA by I part of this array contains the C transformed null vector V(s) = Q * Y(s), where Q is the C product of the (J-1) Householder transformations Q(m). C C WORKSPACE C C RWORK - DOUBLE PRECISION array of DIMENSION (2*MA). C C METHOD C C Let C Q(m) = ( I 0 ) C ( 0 P(m) ) C be the elementary Householder transformation corresponding to the C pivot A(m,S(m)), augmented such that Q(m) is NA by NA, then C Q(1) Q(2) ... Q(J-1) Y ==> Y is computed. C C CONTRIBUTOR C C A.J. Geurts. C C REVISIONS C C 1992, October 27. C 1994, December 8. C C IMPLICIT NONE C .. Parameters .. DOUBLE PRECISION ZERO, ONE PARAMETER (ZERO = 0.0D0, ONE = 1.0D0) C .. Scalar Arguments .. INTEGER MA, NA, I, J, LDA, LDY C .. Array Arguments .. INTEGER S(*) DOUBLE PRECISION A(LDA,NA), Y(LDY,I), RWORK(2*NA) C .. Local Scalars .. INTEGER LEN, M, M1, MY DOUBLE PRECISION TAU C .. External Subroutines/Functions .. EXTERNAL DCOPY, DGEMV, DGER, DLAVEC C C .. Executable Statements .. C C Compute the matrix P(m) Y by w = Y'u and next Y = Y - tau * uw'. C for m = j-1, ... , 1. C DO 20 M = J - 1, 1, -1 IF (A(M+1,S(M)) .NE. ZERO) THEN LEN = MA - M + 1 MY = NA - MA + M M1 = MA + 1 CALL DCOPY(LEN, A(M+1,S(M)), 1, RWORK, 1) CALL DLAVEC(LEN, ZERO, RWORK(M1), 1) TAU = RWORK(1) RWORK(1) = ONE CALL DGEMV('T', LEN, I, ONE, Y(MY,1), LDY, RWORK, 1, * ZERO, RWORK(M1), 1) CALL DGER(LEN, I, -TAU, RWORK, 1, RWORK(M1), 1, Y(MY,1), * LDY) END IF 20 CONTINUE RETURN C *** Last line of COMPTV *** END C SUBROUTINE COMPTY(NP, NA, I, NYI, YI, MU, A, LDA, S, Y, LDY, SK, * RWORK, IERR) C C PURPOSE C C To compute a right null vector of the pencil sA' - E, where the C left part of A' is in staircase form. Actually, the computed C vector is the null vector of the corresponding left part of the C pencil. C REMARK: This auxiliary routine is intended to be called only from C the routine COLRED. C C ARGUMENTS IN C C NP - INTEGER. C The number of columns of the polynomial matrix. C NP >= 1. C NA - INTEGER. C The number of columns of the matrix A. C NA >= 1. C I - INTEGER. C The index of the current diagonal block A(i,i) of the C matrix A being transformed into staircase form. C I >= 1. C NYI - INTEGER. C The length of the righthandside vector YI. C 1 <= NYI <= NP. C YI - DOUBLE PRECISION array of DIMENSION at least (NP). C The right null vector of A(i,i). C MU - INTEGER array of DIMENSION at least (MA). C MU(k), k = 1, ..., i must contain the row index of the left C upper element of A(k,k). C A - DOUBLE PRECISION array of DIMENSION (LDA,NA). C The leading MU(i) - 1 by MU(i) + NP - 1 part of this array C must contain the part of the matrix A' which is in C staircase form. C LDA - INTEGER. C The leading dimension of array A as declared in the calling C program. C LDA >= MU(i) - 1. C S - INTEGER array of DIMENSION at least (MA). C The leading MU(i) - 1 elements of this array must contain C the column indices of the pivots of the right invertible C diagonal matrices A(k,k), k = 1, ..., i-1. C C ARGUMENTS OUT C C Y - DOUBLE PRECISION array of DIMENSION (LDY,NA). C The leading NA by i part of this array contains the C computed polynomial right null vector Y(s) of sA' - E, C where the j-th column contains the coefficient of s**(j-1). C The last NA - MU(i) - NP + 1 components of Y(s) are zero. C LDY - INTEGER. C The leading dimension of array Y as declared in the calling C program. C LDY >= NA. C C WORK SPACE C C SK - INTEGER array of DIMENSION at least (NP). C RWORK - DOUBLE PRECISION array of DIMENSION at least (NP). C C ERROR INDICATOR C C IERR - INTEGER. C Unless the routine detects an error (see next section), C IERR contains 0 on exit. C C WARNINGS AND ERRORS DETECTED BY THE ROUTINE C C IERR = k : A(k,k) is not right invertible. C C METHOD C C Let the pencil sA' - E, partially transformed up to the i-th C block, be C C ( sA(1,1) sA(1,2)-E(1) sA(1,3) . . . sA(1,i) . . . ) C ( sA(2,2) sA(2,3)-E(2) . . . sA(2,i) . . . ) C ( . . ) C ( . . ) C ( . . ) C ( sA(i,i) . . . ) C ( . . . ) C ( . . . ) C C where A(k,k), k = 1,..., i is right invertible, A(k,k+1) is square C and E(k) = I of appropriate size. C Let Y(i,i), the (constant) right null vector of A(i,i), be given. C Then the routine computes a right null vector of sA' - E of the C form C C ( Y(1,1) . Y(1,k) . . . Y(1,j) . . Y(1,i) ) C ( . C ( . C ( Y(k,k) . . . Y(k,j) . . Y(k,i) ) C ( . C ( . C ( . Y(i-1,i) ) C ( Y(i,i) ) C ( 0 0 ) C C by comparing coefficients of equal power of s in the formula C i j k C A(k,k) SUM Y(k,j) * s = Y(k+1,k+1) * s + C j=k C i-1 j j i i C SUM (Y(k+1,j+1) - SUM A(k,l) Y(l,j)) * s - SUM A(k,l) Y(l,i) * s C j=k+1 l=k+1 l=k+1 C C Y(k,j), j = k,..., i, is a vector of length MU(k) - MU(k-1). C C CONTRIBUTOR C C A.J. Geurts. C C REVISIONS C C 1992, October 27. C 1994, December 8. C C IMPLICIT NONE C .. Parameters .. DOUBLE PRECISION ZERO, ONE PARAMETER (ZERO = 0.0D0, ONE = 1.0D0) C .. Scalar Arguments .. INTEGER NP, NA, I, NYI, LDA, LDY, IERR C .. Array Arguments .. INTEGER MU(*), S(*), SK(*) DOUBLE PRECISION YI(*), A(LDA,*), Y(LDY,*), RWORK(*) C .. Local Scalars .. INTEGER J, K, M, INDEXK, MUK, MUK1, NUK, MUK1NP C .. External Subroutines .. EXTERNAL COMPYI, DCOPY, DGEMV, DLAVEC, DLASET C C .. Executable Statements .. C IERR = 0 C C Initialization of the polynomial null vector Y(s). C CALL DLASET('G', NA, I, ZERO, ZERO, Y, LDY) K = I IF (I .EQ. 1) THEN CALL DCOPY(NYI, YI, 1, Y(1,1), 1) ELSE INDEXK = NP + MU(I-1) CALL DCOPY(NYI, YI, 1, Y(INDEXK,I), 1) DO 30 K = I - 1, 1, -1 MUK = MU(K) INDEXK = NP + MUK NUK = MU(K+1) - MUK DO 20 J = K, I C C Compute the righthandside for Y(k,j) and store the C result in RWORK. C IF (J .LT. I) THEN CALL DCOPY(NUK, Y(INDEXK,J+1), 1, RWORK, 1) ELSE CALL DLAVEC(NUK, ZERO, RWORK, 1) END IF IF (J .GT. K) THEN CALL DGEMV('N', NUK, MU(J)-MUK, -ONE, A(MUK,INDEXK), * LDA, Y(INDEXK,J), 1, ONE, RWORK, 1) END IF C C Solve A(k,k) * Y(k,j) = RWORK for Y(k,j). C IF (K .EQ. 1) THEN CALL COMPYI(NUK, NP, A(1,1), LDA, S, RWORK, Y(1,J), * IERR) ELSE MUK1 = MU(K-1) MUK1NP = MUK1 + NP DO 10 M = 1, NUK SK(M) = S(MUK-1+M) - MUK1NP + 1 10 CONTINUE CALL COMPYI(NUK, MUK-MUK1, A(MUK,MUK1NP), LDA, SK, * RWORK, Y(MUK1NP,J), IERR) END IF IF (IERR .NE. 0) THEN IERR = K RETURN END IF 20 CONTINUE 30 CONTINUE END IF RETURN C *** Last line of COMPTY *** END C SUBROUTINE COMPYI(M, N, A, LDA, S, V, Y, IERR) C C PURPOSE C C To compute a null vector yi of the right invertible diagonal C submatrix A(i,i), by solving an appropriate M by N system of C linear equations A y = v, where A is in staircase form. C REMARK: This auxiliary routine is intended to be called only from C the routine COLRED. C C ARGUMENTS IN C C M - INTEGER. C The number of rows of matrix A. C M >= 1. C N - INTEGER. C The number of columns of matrix A. C N >= M. C A - DOUBLE PRECISION array of DIMENSION (LDA,N). C The leading M by N part of this array must contain the C matrix A. C LDA - INTEGER. C The leading dimension of array A as declared by the calling C program. C LDA >= M. C S - INTEGER array of DIMENSION at least (M). C S(i), i = 1,...,M, must contain the column index of the corner C in the i-th row of A. C V - DOUBLE PRECISION array of DIMENSION at least (M). C The righthand-side of the system of linear equations. C C ARGUMENTS OUT C C Y - DOUBLE PRECISION array of DIMENSION at least (N). C The computed solution of the system of linear equation. C C ERROR INDICATOR C C IERR - INTEGER. C Unless the routine detects an error (see next section), C IERR contains 0 on exit. C C WARNINGS AMD ERRORS DETECTED BY THE ROUTINE C C IERR = 3 : The matrix A is not right invertible. C C METHOD C C Let A * P = ( B | Z ) where P is a permutation matrix such that B C is nonsingular upper triangular. Z contains the remaining columns C of A. Then the system B x = v is solved and y = P ( x | 0 )'. C C CONTRIBUTOR C C A.J. Geurts. C C REVISIONS C C 1992, October 27. C 1994, December 8. C C IMPLICIT NONE C .. Parameters .. DOUBLE PRECISION ZERO PARAMETER (ZERO = 0.0D0) C .. Scalar Arguments .. INTEGER M, N, LDA, IERR C .. Array Arguments .. INTEGER S(M) DOUBLE PRECISION A(LDA,N), V(M), Y(N) C .. Local Scalars .. INTEGER J, K, SI DOUBLE PRECISION SUM LOGICAL FAIL C .. External Subroutines .. EXTERNAL DLARDV, DLAVEC DOUBLE PRECISION DLARDV C C .. Executable Statements .. C C Check input parameters. C IF (N .LT. M) THEN IERR = 3 RETURN END IF C IERR = 0 CALL DLAVEC(N, ZERO, Y, 1) SI = S(M) Y(SI) = DLARDV(V(M), A(M,SI), FAIL) IF (FAIL) THEN IERR = 3 RETURN END IF DO 20 K = M - 1, 1, -1 SUM = V(K) DO 10 J = K + 1, M SI = S(J) SUM = SUM - A(K,SI) * Y(SI) 10 CONTINUE SI = S(K) Y(SI) = DLARDV(SUM, A(K,SI), FAIL) IF (FAIL) THEN IERR = 3 RETURN END IF 20 CONTINUE RETURN C *** Last line of COMPYI *** END DOUBLE PRECISION FUNCTION DLARDV(A, B, FAIL) C C -- LAPACK like auxiliary routine C December 8, 1994 C C .. Scalar Arguments .. DOUBLE PRECISION A, B LOGICAL FAIL C .. C C Purpose C ======= C C DLARDV returns the value div given by C C div = ( a/b if a/b does not overflow, C ( C ( 0.0 if a .eq. 0.0, C ( C ( sign( a/b )*flmax if a .ne. 0.0 and a/b would overflow, C C where flmax is a large value, via the function name. In addition if C a/b would overflow then fail is returned as true, otherwise fail is C returned as false. C C Note that when a and b are both zero, fail is returned as true, but C div is returned as 0.0. In all other cases of overflow div is such C that abs( div ) = flmax. C C When b = 0 then sign( a/b ) is taken as sign( a ). C C NB. This routine is derived from the Nag Fortran 77 O( 1 ) basic C linear algebra routine routine F06BLF. C C C .. Parameters .. DOUBLE PRECISION ONE, ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) C .. Local Scalars .. DOUBLE PRECISION ABSB, DIV, FLMAX, FLMIN LOGICAL FIRST C .. External Functions .. DOUBLE PRECISION DLAMCH EXTERNAL DLAMCH C .. Intrinsic Functions .. INTRINSIC ABS, SIGN C .. Save statement .. SAVE FIRST, FLMIN, FLMAX C .. Data statements .. DATA FIRST/ .TRUE. / C .. C .. Executable Statements .. C IF( A.EQ.ZERO )THEN DIV = ZERO IF( B.EQ.ZERO )THEN FAIL = .TRUE. ELSE FAIL = .FALSE. END IF ELSE C IF( FIRST )THEN FIRST = .FALSE. FLMIN = DLAMCH('S') FLMAX = 1/FLMIN END IF C IF( B.EQ.ZERO )THEN DIV = SIGN( FLMAX, A ) FAIL = .TRUE. ELSE ABSB = ABS( B ) IF( ABSB.GE.ONE )THEN FAIL = .FALSE. IF( ABS( A ).GE.ABSB*FLMIN )THEN DIV = A/B ELSE DIV = ZERO END IF ELSE IF( ABS( A ).LE.ABSB*FLMAX )THEN FAIL = .FALSE. DIV = A/B ELSE FAIL = .TRUE. DIV = FLMAX IF( ( ( A.LT.ZERO ).AND.( B.GT.ZERO ) ).OR. $ ( ( A.GT.ZERO ).AND.( B.LT.ZERO ) ) ) $ DIV = -DIV END IF END IF END IF END IF C DLARDV = DIV RETURN C C End of DLARDV C END SUBROUTINE DLARFV( N, V, INCV, TAU, ALPHA, X, INCX ) C C -- LAPACK like auxiliary routine C December 8, 1994 C C .. Scalar Arguments .. DOUBLE PRECISION ALPHA, TAU INTEGER INCX, INCV, N C .. Array Arguments .. DOUBLE PRECISION X( * ), V( * ) C .. C C Purpose C ======= C C DLARFV performs a Householder reflection given by C C ( alpha ) = H * ( alpha ) , C ( x ) ( x ) C C where the orthogonal matrix H is given in the form C C H = I - tau * ( 1 ) * ( 1 v') , C ( v ) C C where tau is a real scalar and v is a real (n-1)-element. C If tau is zero then H is assumed to be the unit matrix and the C transformation is skipped, otherwise tau must be in the range C ( 1.0, 2.0 ). tau and v will usually be supplied by routine DLARFG. C C NB. This routine is derived from the Nag Fortran 77 O( n ) basic C linear agebra routine F06FUF. C C .. Parameters .. DOUBLE PRECISION ZERO PARAMETER ( ZERO = 0.0D+0 ) C .. Local Scalars .. DOUBLE PRECISION BETA C .. External Functions .. DOUBLE PRECISION DDOT EXTERNAL DDOT C .. External Subroutines .. EXTERNAL DAXPY C .. C .. Executable Statements .. IF( TAU.NE.ZERO ) THEN BETA = ALPHA + DDOT( N, V, INCV, X, INCX ) ALPHA = ALPHA - TAU * BETA CALL DAXPY( N-1, -TAU * BETA, V, INCV, X, INCX ) END IF C RETURN C C End of DLARFV C END SUBROUTINE DLAVEC( N, CONST, X, INCX ) C C -- LAPACK like auxiliary routine C August 31, 1994 C C .. Scalar Arguments .. DOUBLE PRECISION CONST INTEGER INCX, N C .. Array Arguments .. DOUBLE PRECISION X( * ) C .. C C Purpose C ======= C C DLAVEC performs the operation C C x = const*e, e' = ( 1 1 ... 1 ). C C NB. This routine is derived from the Nag Fortran 77 O( n ) basic C linear algebra routine F06FBF. C C .. Parameters .. DOUBLE PRECISION ZERO PARAMETER ( ZERO = 0.0D+0 ) C .. Local Scalars .. INTEGER IX C .. C C .. Executable Statements .. IF( N.GT.0 ) THEN IF( CONST.NE.ZERO ) THEN DO 10, IX = 1, 1 + ( N - 1 )*INCX, INCX X( IX ) = CONST 10 CONTINUE ELSE DO 20, IX = 1, 1 + ( N - 1 )*INCX, INCX X( IX ) = ZERO 20 CONTINUE END IF END IF C RETURN C C End of DLAVEC C END C SUBROUTINE HHTRAN(MA, NA, L, K, A, LDA, Q, RWORK) C C PURPOSE C C To compute the Householder reflection Q which transforms a C subcolumn of A into the first unit vector and to apply Q left and C right to appropriate submatrices of A. C REMARK: This auxiliary routine is intended to be called only from C the routine COLRED. C C ARGUMENTS IN C C MA - INTEGER. C The number of rows of matrix A. C MA >= 1. C NA - INTEGER. C The number of columns of matrix A. C NA >= MA. C L - INTEGER. C The index of the column of A to be transformed. C L >= 1. C K - INTEGER. C The row index from which the column is to be transformed. C K >= 1. C A - DOUBLE PRECISION array of DIMENSION (LDA,NA). C The leading MA by NA part of this array must contain the C matrix A. The left lower MA - K + 1 by L - 1 block of the C matrix A is understood to be zero made by former C transformations. C Note: this array is overwritten. C LDA - INTEGER. C The leading dimension of array A as declared in the calling C program. C LDA >= MA. C C ARGUMENTS OUT C C A - DOUBLE PRECISION array of DIMENSION (LDA,NA). C The leading MA by NA part of this array contains the C transformed matrix Q A Q , where Q is the Householder C 1 2 i C transformation appropriately augmented with an identity C matrix. C Q - DOUBLE PRECISION array of DIMENSION at least (MA-K+1). C The leading MA - K + 1 elements of this array contain the C constant tau and the vector v which define the Householder C reflection (see DLARGF, i.e., Q(1) contains the value tau C and Q(i), i=2,...MA-K+1 contain the vector v). C C WORKSPACE C C RWORK - DOUBLE PRECISION array of DIMENSION at least (MA-K+1). C C METHOD C C Let w be the subcolumn A(i,L), i=K,...,MA. Then the Householder C reflection P = I - tau * u * u', where u = (1, v')', such that C Pw = e1 (the first unit vector) is computed. C Let C Q = ( I 0 ) Q = ( I 0 ) C 1 ( 0 P ) 2 ( 0 P ) C be such that Q is MA by MA and Q is NA by NA. C 1 2 C Then Q A Q is computed and stored in A. C 1 2 C C CONTRIBUTOR C C A.J. Geurts. C C REVISIONS C C 1992, October 27. C 1994, December 8. C C IMPLICIT NONE C .. Parameters .. DOUBLE PRECISION ZERO, ONE PARAMETER (ZERO = 0.0D0, ONE = 1.0D0) C .. Scalar Arguments .. INTEGER MA, NA, K, L, LDA C .. Array Arguments .. DOUBLE PRECISION A(LDA,NA), Q(MA-K+1), RWORK(MA-K+1) C .. Local Scalars .. INTEGER LEN, NMK, NL1 DOUBLE PRECISION TAU C .. External Subroutines/Functions .. EXTERNAL DCOPY, DGEMV, DGER, DLARFG, DLAVEC C C .. Executable Statements .. C LEN = MA - K + 1 IF (LEN .GT. 1) THEN C C Generate the Householder transformation. C CALL DCOPY(LEN, A(K,L), 1, Q, 1) CALL DLARFG(LEN, Q(1), Q(2), 1, TAU) C IF (TAU .NE. ZERO) THEN C C Compute matrix Q A by w = A'u and next A = A - tau * uw'. C 1 Q(1) = ONE NL1 = NA - L + 1 CALL DLAVEC(MA, ZERO, RWORK, 1) CALL DGEMV('T', LEN, NL1, ONE, A(K,L), LDA, Q, 1, ZERO, * RWORK, 1) CALL DGER(LEN, NL1, -TAU, Q, 1, RWORK, 1, A(K,L), LDA) C C Compute matrix Q A Q by w = Au and next A = A - tau * wu'. C 1 2 NMK = NA - MA + K CALL DGEMV('N', MA, LEN, ONE, A(1,NMK), LDA, Q, 1, ZERO, * RWORK, 1) CALL DGER(MA, LEN, -TAU, RWORK, 1, Q, 1, A(1,NMK), LDA) Q(1) = TAU END IF ELSE Q(1) = ZERO END IF RETURN C *** Last line of HHTRAN *** END C SUBROUTINE MKPENC(MP, NP, DP, P, LDP1, LDP2, A, LDA, MA, NA) C C PURPOSE C C Given an MP x NP polynomial matrix of degree dp C dp-1 dp C P(s) = P(0) + ... + P(dp-1) * s + P(dp) * s (1) C C the subroutine MKPENC constructs the first degree part A of the C linearization of the polynomial matrix P(s), where C C | P(dp) 0 . . 0 | C | P(dp-1) I 0 . 0 | C A = | . 0 I . . 0 | (2) C | . . . . . | C | . I 0 0 | C | P(0) 0 I 0 | C C REMARK: This auxiliary routine is intended to be called only from C the routine COLRED. C C ARGUMENTS IN C C MP - INTEGER. C The number of rows of the polynomial matrix P(s). C MP >= 1. C NP - INTEGER. C The number of columns of the polynomial matrix P(s). C NP >= 1. C DP - INTEGER. C The degree of the polynomial matrix P(s). C DP >= 1. C P - DOUBLE PRECISION array of DIMENSION (LDP1,LP2,DP+1). C The leading MP by NP by (DP+1) part of this array must C contain the coefficients of the polynomial matrix P(s). C Specifically, P(i,j,k) must contain the coefficient of C s**(k-1) of the polynomial which is the (i,j)-th element C of P(s), where i = 1,2,...,MP, j = 1,2,...,NP and C k = 1,2,...,DP+1. C LDP1 - INTEGER. C The leading dimension of array P as declared in the calling C program. C LDP1 >= MP. C LDP2 - INTEGER. C The second dimension of array P as declared in the calling C program. C LDP2 >= NP. C C ARGUMENTS OUT C C A - DOUBLE PRECISION array of DIMENSION (LDA,(DP+1)*MP+NP). C The leading (DP+1)*MP by (DP+1)*MP+NP part of this array C contains the matrix A as described in (2). C LDA - INTEGER. C The leading dimension of array A as declared in the calling C program. C LDA >= (DP+1)*MP. C MA - INTEGER. C The number of rows of matrix A. C NA - INTEGER. C The number of columns of matrix A. C C CONTRIBUTOR C C A.J.Geurts. C C REVISIONS C C 1992, October 27. C 1994, December 8. C C IMPLICIT NONE C .. Parameters .. DOUBLE PRECISION ZERO, ONE PARAMETER (ZERO = 0.0D0, ONE = 1.0D0) C .. Scalar Arguments .. INTEGER MP, NP, DP, LDP1, LDP2, LDA, MA, NA C .. Array Arguments .. DOUBLE PRECISION P(LDP1,LDP2,*), A(LDA,*) C .. Local Scalars .. INTEGER J, M1, NJ C .. External Subroutines/Functions .. EXTERNAL DLACPY, DLASET C C .. Executable Statements .. C C Initialisation of the matrix A. C M1 = DP * MP MA = M1 + MP NA = MA + NP CALL DLASET('G', MP, MA, ZERO, ZERO, A(1,NP+1), LDA) CALL DLASET('G', M1, MA, ZERO, ONE, A(MP+1,NP+1), LDA) C C Insert the matrices P(0), P(1), ..., P(pd) in array A. C NJ = M1 + 1 DO 20 J = 1, DP + 1 CALL DLACPY('G', MP, NP, P(1,1,J), LDP1, A(NJ,1), LDA) NJ = NJ - MP 20 CONTINUE C RETURN C *** Last line of MKPENC *** END SHAR_EOF fi # end of overwriting check cd .. cd .. # End of shell archive exit 0