C ALGORITHM 747, COLLECTED ALGORITHMS FROM ACM. C THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE, C VOL. 21, NO. 3, September, 1995, P. 299-326. C C This file contains 36 files separated by lines of the form C C*** filename C C The filenames in this file are: C C blas.f ddemo.f ddemo.sh C dlog.ref dmevas.f dstair.f C eispk.f lapack.f makefile C readme sdemo.f sdemo.sh C slog.ref smevas.f sstair.f C test.doc test01.dat test02.dat C test03.dat test04.dat test05.dat C test06.dat test07.dat test08.dat C test09.dat test10.dat test11.dat C test12.dat test13.dat test14.dat C test15.dat test16.dat test17.dat C test18.dat test19.dat test20.dat C C C*** blas.f c c FILE: Blas.f c c********************************************************** c********************************************************** c 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 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 c c********************************************************** c********************************************************** c 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 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 c c********************************************************** c********************************************************** c subroutine drot (n,dx,incx,dy,incy,c,s) c c applies a plane rotation. c jack dongarra, linpack, 3/11/78. c double precision dx(*),dy(*),dtemp,c,s integer i,incx,incy,ix,iy,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 not equal c 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 = c*dx(ix) + s*dy(iy) dy(iy) = c*dy(iy) - s*dx(ix) dx(ix) = dtemp ix = ix + incx iy = iy + incy 10 continue return c c code for both increments equal to 1 c 20 do 30 i = 1,n dtemp = c*dx(i) + s*dy(i) dy(i) = c*dy(i) - s*dx(i) dx(i) = dtemp 30 continue return end c c********************************************************** c********************************************************** c subroutine drotg(da,db,c,s) c c construct givens plane rotation. c jack dongarra, linpack, 3/11/78. c modified 9/27/86. c double precision da,db,c,s,roe,scale,r,z c roe = db if( dabs(da) .gt. dabs(db) ) roe = da scale = dabs(da) + dabs(db) if( scale .ne. 0.0d0 ) go to 10 c = 1.0d0 s = 0.0d0 r = 0.0d0 go to 20 10 r = scale*dsqrt((da/scale)**2 + (db/scale)**2) r = dsign(1.0d0,roe)*r c = da/r s = db/r 20 z = s if( dabs(c) .gt. 0.0d0 .and. dabs(c) .le. s ) z = 1.0d0/c da = r db = z return end c c********************************************************** c********************************************************** c subroutine dswap (n,dx,incx,dy,incy) c c interchanges two vectors. c uses unrolled loops for increments equal one. c jack dongarra, linpack, 3/11/78. c double precision dx(*),dy(*),dtemp 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 not equal c 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 = dx(ix) dx(ix) = dy(iy) dy(iy) = dtemp 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,3) if( m .eq. 0 ) go to 40 do 30 i = 1,m dtemp = dx(i) dx(i) = dy(i) dy(i) = dtemp 30 continue if( n .lt. 3 ) return 40 mp1 = m + 1 do 50 i = mp1,n,3 dtemp = dx(i) dx(i) = dy(i) dy(i) = dtemp dtemp = dx(i + 1) dx(i + 1) = dy(i + 1) dy(i + 1) = dtemp dtemp = dx(i + 2) dx(i + 2) = dy(i + 2) dy(i + 2) = dtemp 50 continue return end c ccccc ****************************************************************** ccccc ****************************************************************** c 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 double precision da,dx(*) integer i,incx,m,mp1,n,nincx c if(n.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 c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c 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 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 c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c double precision function dnrm2 ( n, dx, incx) integer incx, n, next double precision dx(*), cutlo, cuthi, hitest, sum, xmax,zero,one data zero, one /0.0d0, 1.0d0/ c c euclidean norm of the n-vector stored in dx() with storage c increment incx . c if n .le. 0 return with result = 0. c if n .ge. 1 then incx must be .ge. 1 c c c.l.lawson, 1978 jan 08 c c four phase method using two built-in constants that are c hopefully applicable to all machines. c cutlo = maximum of dsqrt(u/eps) over all known machines. c cuthi = minimum of dsqrt(v) over all known machines. c where c eps = smallest no. such that eps + 1. .gt. 1. c u = smallest positive no. (underflow limit) c v = largest no. (overflow limit) c c brief outline of algorithm.. c c phase 1 scans zero components. c move to phase 2 when a component is nonzero and .le. cutlo c move to phase 3 when a component is .gt. cutlo c move to phase 4 when a component is .ge. cuthi/m c where m = n for x() real and m = 2*n for complex. c c values for cutlo and cuthi.. c from the environmental parameters listed in the imsl converter c document the limiting values are as follows.. c cutlo, s.p. u/eps = 2**(-102) for honeywell. close seconds are c univac and dec at 2**(-103) c thus cutlo = 2**(-51) = 4.44089e-16 c cuthi, s.p. v = 2**127 for univac, honeywell, and dec. c thus cuthi = 2**(63.5) = 1.30438e19 c cutlo, d.p. u/eps = 2**(-67) for honeywell and dec. c thus cutlo = 2**(-33.5) = 8.23181d-11 c cuthi, d.p. same as s.p. cuthi = 1.30438d19 c data cutlo, cuthi / 8.232d-11, 1.304d19 / c data cutlo, cuthi / 4.441e-16, 1.304e19 / data cutlo, cuthi / 8.232d-11, 1.304d19 / c if(n .gt. 0) go to 10 dnrm2 = zero go to 300 c 10 assign 30 to next sum = zero nn = n * incx c begin main loop i = 1 20 go to next,(30, 50, 70, 110) 30 if( dabs(dx(i)) .gt. cutlo) go to 85 assign 50 to next xmax = zero c c phase 1. sum is zero c 50 if( dx(i) .eq. zero) go to 200 if( dabs(dx(i)) .gt. cutlo) go to 85 c c prepare for phase 2. assign 70 to next go to 105 c c prepare for phase 4. c 100 i = j assign 110 to next sum = (sum / dx(i)) / dx(i) 105 xmax = dabs(dx(i)) go to 115 c c phase 2. sum is small. c scale to avoid destructive underflow. c 70 if( dabs(dx(i)) .gt. cutlo ) go to 75 c c common code for phases 2 and 4. c in phase 4 sum is large. scale to avoid overflow. c 110 if( dabs(dx(i)) .le. xmax ) go to 115 sum = one + sum * (xmax / dx(i))**2 xmax = dabs(dx(i)) go to 200 c 115 sum = sum + (dx(i)/xmax)**2 go to 200 c c c prepare for phase 3. c 75 sum = (sum * xmax) * xmax c c c for real or d.p. set hitest = cuthi/n c for complex set hitest = cuthi/(2*n) c 85 hitest = cuthi/float( n ) c c phase 3. sum is mid-range. no scaling. c do 95 j =i,nn,incx if(dabs(dx(j)) .ge. hitest) go to 100 95 sum = sum + dx(j)**2 dnrm2 = dsqrt( sum ) go to 300 c 200 continue i = i + incx if ( i .le. nn ) go to 20 c c end of main loop. c c compute square root and adjust for scaling. c dnrm2 = xmax * dsqrt(sum) 300 continue return end c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c double precision function dasum(n,dx,incx) c c takes the sum of the absolute values. c jack dongarra, linpack, 3/11/78. c double precision dx(*),dtemp integer i,incx,m,mp1,n,nincx c dasum = 0.0d0 dtemp = 0.0d0 if(n.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 dtemp = dtemp + dabs(dx(i)) 10 continue dasum = dtemp return c c code for increment equal to 1 c c c clean-up loop c 20 m = mod(n,6) if( m .eq. 0 ) go to 40 do 30 i = 1,m dtemp = dtemp + dabs(dx(i)) 30 continue if( n .lt. 6 ) go to 60 40 mp1 = m + 1 do 50 i = mp1,n,6 dtemp = dtemp + dabs(dx(i)) + dabs(dx(i + 1)) + dabs(dx(i + 2)) * + dabs(dx(i + 3)) + dabs(dx(i + 4)) + dabs(dx(i + 5)) 50 continue 60 dasum = dtemp return end c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c LOGICAL FUNCTION LSAME ( CA, CB ) * .. Scalar Arguments .. CHARACTER*1 CA, CB * .. * * Purpose * ======= * * LSAME tests if CA is the same letter as CB regardless of case. * CB is assumed to be an upper case letter. LSAME returns .TRUE. if * CA is either the same as CB or the equivalent lower case letter. * * N.B. This version of the routine is only correct for ASCII code. * Installers must modify the routine for other character-codes. * * For EBCDIC systems the constant IOFF must be changed to -64. * For CDC systems using 6-12 bit representations, the system- * specific code in comments must be activated. * * Parameters * ========== * * CA - CHARACTER*1 * CB - CHARACTER*1 * On entry, CA and CB specify characters to be compared. * Unchanged on exit. * * * Auxiliary routine for Level 2 Blas. * * -- Written on 20-July-1986 * Richard Hanson, Sandia National Labs. * Jeremy Du Croz, Nag Central Office. * * .. Parameters .. INTEGER IOFF PARAMETER ( IOFF=32 ) * .. Intrinsic Functions .. INTRINSIC ICHAR * .. Executable Statements .. * * Test if the characters are equal * LSAME = CA .EQ. CB * * Now test for equivalence * IF ( .NOT.LSAME ) THEN LSAME = ICHAR(CA) - IOFF .EQ. ICHAR(CB) END IF * RETURN * * The following comments contain code for CDC systems using 6-12 bit * representations. * * .. Parameters .. * INTEGER ICIRFX * PARAMETER ( ICIRFX=62 ) * .. Scalar Arguments .. * CHARACTER*1 CB * .. Array Arguments .. * CHARACTER*1 CA(*) * .. Local Scalars .. * INTEGER IVAL * .. Intrinsic Functions .. * INTRINSIC ICHAR, CHAR * .. Executable Statements .. * * See if the first character in string CA equals string CB. * * LSAME = CA(1) .EQ. CB .AND. CA(1) .NE. CHAR(ICIRFX) * * IF (LSAME) RETURN * * The characters are not identical. Now check them for equivalence. * Look for the 'escape' character, circumflex, followed by the * letter. * * IVAL = ICHAR(CA(2)) * IF (IVAL.GE.ICHAR('A') .AND. IVAL.LE.ICHAR('Z')) THEN * LSAME = CA(1) .EQ. CHAR(ICIRFX) .AND. CA(2) .EQ. CB * END IF * * RETURN * * End of LSAME. * END c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c SUBROUTINE XERBLA ( SRNAME, INFO ) * .. Scalar Arguments .. INTEGER INFO CHARACTER*6 SRNAME * .. * * Purpose * ======= * * XERBLA is an error handler for the Level 2 BLAS routines. * * It is called by the Level 2 BLAS routines if an input parameter is * invalid. * * Installers should consider modifying the STOP statement in order to * call system-specific exception-handling facilities. * * Parameters * ========== * * SRNAME - CHARACTER*6. * On entry, SRNAME specifies the name of the routine which * called XERBLA. * * INFO - INTEGER. * On entry, INFO specifies the position of the invalid * parameter in the parameter-list of the calling routine. * * * Auxiliary routine for Level 2 Blas. * * Written on 20-July-1986. * * .. Executable Statements .. * WRITE (*,99999) SRNAME, INFO * STOP * 99999 FORMAT ( ' ** On entry to ', A6, ' parameter number ', I2, $ ' had an illegal value' ) * * End of XERBLA. * END c ccccc ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc * ************************************************************************ * SUBROUTINE DTRSV ( UPLO, TRANS, DIAG, N, A, LDA, X, INCX ) * .. Scalar Arguments .. INTEGER INCX, LDA, N CHARACTER*1 DIAG, TRANS, UPLO * .. Array Arguments .. DOUBLE PRECISION A( LDA, * ), X( * ) * .. * * Purpose * ======= * * DTRSV solves one of the systems of equations * * A*x = b, or A'*x = b, * * where b and x are n element vectors and A is an n by n unit, or * non-unit, upper or lower triangular matrix. * * No test for singularity or near-singularity is included in this * routine. Such tests must be performed before calling this routine. * * Parameters * ========== * * UPLO - CHARACTER*1. * On entry, UPLO specifies whether the matrix is an upper or * lower triangular matrix as follows: * * UPLO = 'U' or 'u' A is an upper triangular matrix. * * UPLO = 'L' or 'l' A is a lower triangular matrix. * * Unchanged on exit. * * TRANS - CHARACTER*1. * On entry, TRANS specifies the equations to be solved as * follows: * * TRANS = 'N' or 'n' A*x = b. * * TRANS = 'T' or 't' A'*x = b. * * TRANS = 'C' or 'c' A'*x = b. * * Unchanged on exit. * * DIAG - CHARACTER*1. * On entry, DIAG specifies whether or not A is unit * triangular as follows: * * DIAG = 'U' or 'u' A is assumed to be unit triangular. * * DIAG = 'N' or 'n' A is not assumed to be unit * triangular. * * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the order of the matrix A. * N must be at least zero. * Unchanged on exit. * * A - DOUBLE PRECISION array of DIMENSION ( LDA, n ). * Before entry with UPLO = 'U' or 'u', the leading n by n * upper triangular part of the array A must contain the upper * triangular matrix and the strictly lower triangular part of * A is not referenced. * Before entry with UPLO = 'L' or 'l', the leading n by n * lower triangular part of the array A must contain the lower * triangular matrix and the strictly upper triangular part of * A is not referenced. * Note that when DIAG = 'U' or 'u', the diagonal elements of * A are not referenced either, but are assumed to be unity. * 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, n ). * Unchanged on exit. * * X - DOUBLE PRECISION array of dimension at least * ( 1 + ( n - 1 )*abs( INCX ) ). * Before entry, the incremented array X must contain the n * element right-hand side vector b. On exit, X is overwritten * with the solution vector x. * * INCX - INTEGER. * On entry, INCX specifies the increment for the elements of * X. INCX 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 ZERO PARAMETER ( ZERO = 0.0D+0 ) * .. Local Scalars .. DOUBLE PRECISION TEMP INTEGER I, INFO, IX, J, JX, KX LOGICAL NOUNIT * .. 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( UPLO , 'U' ).AND. $ .NOT.LSAME( UPLO , 'L' ) )THEN INFO = 1 ELSE IF( .NOT.LSAME( TRANS, 'N' ).AND. $ .NOT.LSAME( TRANS, 'T' ).AND. $ .NOT.LSAME( TRANS, 'C' ) )THEN INFO = 2 ELSE IF( .NOT.LSAME( DIAG , 'U' ).AND. $ .NOT.LSAME( DIAG , 'N' ) )THEN INFO = 3 ELSE IF( N.LT.0 )THEN INFO = 4 ELSE IF( LDA.LT.MAX( 1, N ) )THEN INFO = 6 ELSE IF( INCX.EQ.0 )THEN INFO = 8 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'DTRSV ', INFO ) RETURN END IF * * Quick return if possible. * IF( N.EQ.0 ) $ RETURN * NOUNIT = LSAME( DIAG, 'N' ) * * Set up the start point in X if the increment is not unity. This * will be ( N - 1 )*INCX too small for descending loops. * IF( INCX.LE.0 )THEN KX = 1 - ( N - 1 )*INCX ELSE IF( INCX.NE.1 )THEN KX = 1 END IF * * Start the operations. In this version the elements of A are * accessed sequentially with one pass through A. * IF( LSAME( TRANS, 'N' ) )THEN * * Form x := inv( A )*x. * IF( LSAME( UPLO, 'U' ) )THEN IF( INCX.EQ.1 )THEN DO 20, J = N, 1, -1 IF( X( J ).NE.ZERO )THEN IF( NOUNIT ) $ X( J ) = X( J )/A( J, J ) TEMP = X( J ) DO 10, I = J - 1, 1, -1 X( I ) = X( I ) - TEMP*A( I, J ) 10 CONTINUE END IF 20 CONTINUE ELSE JX = KX + ( N - 1 )*INCX DO 40, J = N, 1, -1 IF( X( JX ).NE.ZERO )THEN IF( NOUNIT ) $ X( JX ) = X( JX )/A( J, J ) TEMP = X( JX ) IX = JX DO 30, I = J - 1, 1, -1 IX = IX - INCX X( IX ) = X( IX ) - TEMP*A( I, J ) 30 CONTINUE END IF JX = JX - INCX 40 CONTINUE END IF ELSE IF( INCX.EQ.1 )THEN DO 60, J = 1, N IF( X( J ).NE.ZERO )THEN IF( NOUNIT ) $ X( J ) = X( J )/A( J, J ) TEMP = X( J ) DO 50, I = J + 1, N X( I ) = X( I ) - TEMP*A( I, J ) 50 CONTINUE END IF 60 CONTINUE ELSE JX = KX DO 80, J = 1, N IF( X( JX ).NE.ZERO )THEN IF( NOUNIT ) $ X( JX ) = X( JX )/A( J, J ) TEMP = X( JX ) IX = JX DO 70, I = J + 1, N IX = IX + INCX X( IX ) = X( IX ) - TEMP*A( I, J ) 70 CONTINUE END IF JX = JX + INCX 80 CONTINUE END IF END IF ELSE * * Form x := inv( A' )*x. * IF( LSAME( UPLO, 'U' ) )THEN IF( INCX.EQ.1 )THEN DO 100, J = 1, N TEMP = X( J ) DO 90, I = 1, J - 1 TEMP = TEMP - A( I, J )*X( I ) 90 CONTINUE IF( NOUNIT ) $ TEMP = TEMP/A( J, J ) X( J ) = TEMP 100 CONTINUE ELSE JX = KX DO 120, J = 1, N TEMP = X( JX ) IX = KX DO 110, I = 1, J - 1 TEMP = TEMP - A( I, J )*X( IX ) IX = IX + INCX 110 CONTINUE IF( NOUNIT ) $ TEMP = TEMP/A( J, J ) X( JX ) = TEMP JX = JX + INCX 120 CONTINUE END IF ELSE IF( INCX.EQ.1 )THEN DO 140, J = N, 1, -1 TEMP = X( J ) DO 130, I = N, J + 1, -1 TEMP = TEMP - A( I, J )*X( I ) 130 CONTINUE IF( NOUNIT ) $ TEMP = TEMP/A( J, J ) X( J ) = TEMP 140 CONTINUE ELSE KX = KX + ( N - 1 )*INCX JX = KX DO 160, J = N, 1, -1 TEMP = X( JX ) IX = KX DO 150, I = N, J + 1, -1 TEMP = TEMP - A( I, J )*X( IX ) IX = IX - INCX 150 CONTINUE IF( NOUNIT ) $ TEMP = TEMP/A( J, J ) X( JX ) = TEMP JX = JX - INCX 160 CONTINUE END IF END IF END IF * RETURN * * End of DTRSV . * END * ************************************************************************ * SUBROUTINE DTRSM ( SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, $ B, LDB ) * .. Scalar Arguments .. CHARACTER*1 SIDE, UPLO, TRANSA, DIAG INTEGER M, N, LDA, LDB DOUBLE PRECISION ALPHA * .. Array Arguments .. DOUBLE PRECISION A( LDA, * ), B( LDB, * ) * .. * * Purpose * ======= * * DTRSM solves one of the matrix equations * * op( A )*X = alpha*B, or X*op( A ) = alpha*B, * * where alpha is a scalar, X and B are m by n matrices, A is a unit, or * non-unit, upper or lower triangular matrix and op( A ) is one of * * op( A ) = A or op( A ) = A'. * * The matrix X is overwritten on B. * * Parameters * ========== * * SIDE - CHARACTER*1. * On entry, SIDE specifies whether op( A ) appears on the left * or right of X as follows: * * SIDE = 'L' or 'l' op( A )*X = alpha*B. * * SIDE = 'R' or 'r' X*op( A ) = alpha*B. * * Unchanged on exit. * * UPLO - CHARACTER*1. * On entry, UPLO specifies whether the matrix A is an upper or * lower triangular matrix as follows: * * UPLO = 'U' or 'u' A is an upper triangular matrix. * * UPLO = 'L' or 'l' A is a lower triangular matrix. * * Unchanged on exit. * * TRANSA - CHARACTER*1. * On entry, TRANSA specifies the form of op( A ) to be used in * the matrix multiplication as follows: * * TRANSA = 'N' or 'n' op( A ) = A. * * TRANSA = 'T' or 't' op( A ) = A'. * * TRANSA = 'C' or 'c' op( A ) = A'. * * Unchanged on exit. * * DIAG - CHARACTER*1. * On entry, DIAG specifies whether or not A is unit triangular * as follows: * * DIAG = 'U' or 'u' A is assumed to be unit triangular. * * DIAG = 'N' or 'n' A is not assumed to be unit * triangular. * * Unchanged on exit. * * M - INTEGER. * On entry, M specifies the number of rows of B. M must be at * least zero. * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the number of columns of B. N must be * at least zero. * Unchanged on exit. * * ALPHA - DOUBLE PRECISION. * On entry, ALPHA specifies the scalar alpha. When alpha is * zero then A is not referenced and B need not be set before * entry. * Unchanged on exit. * * A - DOUBLE PRECISION array of DIMENSION ( LDA, k ), where k is m * when SIDE = 'L' or 'l' and is n when SIDE = 'R' or 'r'. * Before entry with UPLO = 'U' or 'u', the leading k by k * upper triangular part of the array A must contain the upper * triangular matrix and the strictly lower triangular part of * A is not referenced. * Before entry with UPLO = 'L' or 'l', the leading k by k * lower triangular part of the array A must contain the lower * triangular matrix and the strictly upper triangular part of * A is not referenced. * Note that when DIAG = 'U' or 'u', the diagonal elements of * A are not referenced either, but are assumed to be unity. * Unchanged on exit. * * LDA - INTEGER. * On entry, LDA specifies the first dimension of A as declared * in the calling (sub) program. When SIDE = 'L' or 'l' then * LDA must be at least max( 1, m ), when SIDE = 'R' or 'r' * then LDA must be at least max( 1, n ). * Unchanged on exit. * * B - DOUBLE PRECISION array of DIMENSION ( LDB, n ). * Before entry, the leading m by n part of the array B must * contain the right-hand side matrix B, and on exit is * overwritten by the solution matrix X. * * LDB - INTEGER. * On entry, LDB specifies the first dimension of B as declared * in the calling (sub) program. LDB must be at least * max( 1, m ). * Unchanged on exit. * * * Level 3 Blas routine. * * * -- Written on 8-February-1989. * Jack Dongarra, Argonne National Laboratory. * Iain Duff, AERE Harwell. * Jeremy Du Croz, Numerical Algorithms Group Ltd. * Sven Hammarling, Numerical Algorithms Group Ltd. * * * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. External Subroutines .. EXTERNAL XERBLA * .. Intrinsic Functions .. INTRINSIC MAX * .. Local Scalars .. LOGICAL LSIDE, NOUNIT, UPPER INTEGER I, INFO, J, K, NROWA DOUBLE PRECISION TEMP * .. Parameters .. DOUBLE PRECISION ONE , ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) * .. * .. Executable Statements .. * * Test the input parameters. * LSIDE = LSAME( SIDE , 'L' ) IF( LSIDE )THEN NROWA = M ELSE NROWA = N END IF NOUNIT = LSAME( DIAG , 'N' ) UPPER = LSAME( UPLO , 'U' ) * INFO = 0 IF( ( .NOT.LSIDE ).AND. $ ( .NOT.LSAME( SIDE , 'R' ) ) )THEN INFO = 1 ELSE IF( ( .NOT.UPPER ).AND. $ ( .NOT.LSAME( UPLO , 'L' ) ) )THEN INFO = 2 ELSE IF( ( .NOT.LSAME( TRANSA, 'N' ) ).AND. $ ( .NOT.LSAME( TRANSA, 'T' ) ).AND. $ ( .NOT.LSAME( TRANSA, 'C' ) ) )THEN INFO = 3 ELSE IF( ( .NOT.LSAME( DIAG , 'U' ) ).AND. $ ( .NOT.LSAME( DIAG , 'N' ) ) )THEN INFO = 4 ELSE IF( M .LT.0 )THEN INFO = 5 ELSE IF( N .LT.0 )THEN INFO = 6 ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN INFO = 9 ELSE IF( LDB.LT.MAX( 1, M ) )THEN INFO = 11 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'DTRSM ', INFO ) RETURN END IF * * Quick return if possible. * IF( N.EQ.0 ) $ RETURN * * And when alpha.eq.zero. * IF( ALPHA.EQ.ZERO )THEN DO 20, J = 1, N DO 10, I = 1, M B( I, J ) = ZERO 10 CONTINUE 20 CONTINUE RETURN END IF * * Start the operations. * IF( LSIDE )THEN IF( LSAME( TRANSA, 'N' ) )THEN * * Form B := alpha*inv( A )*B. * IF( UPPER )THEN DO 60, J = 1, N IF( ALPHA.NE.ONE )THEN DO 30, I = 1, M B( I, J ) = ALPHA*B( I, J ) 30 CONTINUE END IF DO 50, K = M, 1, -1 IF( B( K, J ).NE.ZERO )THEN IF( NOUNIT ) $ B( K, J ) = B( K, J )/A( K, K ) DO 40, I = 1, K - 1 B( I, J ) = B( I, J ) - B( K, J )*A( I, K ) 40 CONTINUE END IF 50 CONTINUE 60 CONTINUE ELSE DO 100, J = 1, N IF( ALPHA.NE.ONE )THEN DO 70, I = 1, M B( I, J ) = ALPHA*B( I, J ) 70 CONTINUE END IF DO 90 K = 1, M IF( B( K, J ).NE.ZERO )THEN IF( NOUNIT ) $ B( K, J ) = B( K, J )/A( K, K ) DO 80, I = K + 1, M B( I, J ) = B( I, J ) - B( K, J )*A( I, K ) 80 CONTINUE END IF 90 CONTINUE 100 CONTINUE END IF ELSE * * Form B := alpha*inv( A' )*B. * IF( UPPER )THEN DO 130, J = 1, N DO 120, I = 1, M TEMP = ALPHA*B( I, J ) DO 110, K = 1, I - 1 TEMP = TEMP - A( K, I )*B( K, J ) 110 CONTINUE IF( NOUNIT ) $ TEMP = TEMP/A( I, I ) B( I, J ) = TEMP 120 CONTINUE 130 CONTINUE ELSE DO 160, J = 1, N DO 150, I = M, 1, -1 TEMP = ALPHA*B( I, J ) DO 140, K = I + 1, M TEMP = TEMP - A( K, I )*B( K, J ) 140 CONTINUE IF( NOUNIT ) $ TEMP = TEMP/A( I, I ) B( I, J ) = TEMP 150 CONTINUE 160 CONTINUE END IF END IF ELSE IF( LSAME( TRANSA, 'N' ) )THEN * * Form B := alpha*B*inv( A ). * IF( UPPER )THEN DO 210, J = 1, N IF( ALPHA.NE.ONE )THEN DO 170, I = 1, M B( I, J ) = ALPHA*B( I, J ) 170 CONTINUE END IF DO 190, K = 1, J - 1 IF( A( K, J ).NE.ZERO )THEN DO 180, I = 1, M B( I, J ) = B( I, J ) - A( K, J )*B( I, K ) 180 CONTINUE END IF 190 CONTINUE IF( NOUNIT )THEN TEMP = ONE/A( J, J ) DO 200, I = 1, M B( I, J ) = TEMP*B( I, J ) 200 CONTINUE END IF 210 CONTINUE ELSE DO 260, J = N, 1, -1 IF( ALPHA.NE.ONE )THEN DO 220, I = 1, M B( I, J ) = ALPHA*B( I, J ) 220 CONTINUE END IF DO 240, K = J + 1, N IF( A( K, J ).NE.ZERO )THEN DO 230, I = 1, M B( I, J ) = B( I, J ) - A( K, J )*B( I, K ) 230 CONTINUE END IF 240 CONTINUE IF( NOUNIT )THEN TEMP = ONE/A( J, J ) DO 250, I = 1, M B( I, J ) = TEMP*B( I, J ) 250 CONTINUE END IF 260 CONTINUE END IF ELSE * * Form B := alpha*B*inv( A' ). * IF( UPPER )THEN DO 310, K = N, 1, -1 IF( NOUNIT )THEN TEMP = ONE/A( K, K ) DO 270, I = 1, M B( I, K ) = TEMP*B( I, K ) 270 CONTINUE END IF DO 290, J = 1, K - 1 IF( A( J, K ).NE.ZERO )THEN TEMP = A( J, K ) DO 280, I = 1, M B( I, J ) = B( I, J ) - TEMP*B( I, K ) 280 CONTINUE END IF 290 CONTINUE IF( ALPHA.NE.ONE )THEN DO 300, I = 1, M B( I, K ) = ALPHA*B( I, K ) 300 CONTINUE END IF 310 CONTINUE ELSE DO 360, K = 1, N IF( NOUNIT )THEN TEMP = ONE/A( K, K ) DO 320, I = 1, M B( I, K ) = TEMP*B( I, K ) 320 CONTINUE END IF DO 340, J = K + 1, N IF( A( J, K ).NE.ZERO )THEN TEMP = A( J, K ) DO 330, I = 1, M B( I, J ) = B( I, J ) - TEMP*B( I, K ) 330 CONTINUE END IF 340 CONTINUE IF( ALPHA.NE.ONE )THEN DO 350, I = 1, M B( I, K ) = ALPHA*B( I, K ) 350 CONTINUE END IF 360 CONTINUE END IF END IF END IF * RETURN * * End of DTRSM . * END c ccccc cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c integer function idamax(n,dx,incx) c c finds the index of element having max. absolute value. c jack dongarra, linpack, 3/11/78. c double precision dx(*),dmax integer i,incx,ix,n c idamax = 0 if( n .lt. 1 ) return idamax = 1 if(n.eq.1)return if(incx.eq.1)go to 20 c c code for increment not equal to 1 c ix = 1 dmax = dabs(dx(1)) ix = ix + incx do 10 i = 2,n if(dabs(dx(ix)).le.dmax) go to 5 idamax = i dmax = dabs(dx(ix)) 5 ix = ix + incx 10 continue return c c code for increment equal to 1 c 20 dmax = dabs(dx(1)) do 30 i = 2,n if(dabs(dx(i)).le.dmax) go to 30 idamax = i dmax = dabs(dx(i)) 30 continue return end * ************************************************************************ * * File of the DOUBLE PRECISION Level-3 BLAS. * ========================================== * * SUBROUTINE DGEMM ( TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, * $ BETA, C, LDC ) * * SUBROUTINE DSYMM ( SIDE, UPLO, M, N, ALPHA, A, LDA, B, LDB, * $ BETA, C, LDC ) * * SUBROUTINE DSYRK ( UPLO, TRANS, N, K, ALPHA, A, LDA, * $ BETA, C, LDC ) * * SUBROUTINE DSYR2K( UPLO, TRANS, N, K, ALPHA, A, LDA, B, LDB, * $ BETA, C, LDC ) * * SUBROUTINE DTRMM ( SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, * $ B, LDB ) * * SUBROUTINE DTRSM ( SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, * $ B, LDB ) * * See: * * Dongarra J. J., Du Croz J. J., Duff I. and Hammarling S. * A set of Level 3 Basic Linear Algebra Subprograms. Technical * Memorandum No.88 (Revision 1), Mathematics and Computer Science * Division, Argonne National Laboratory, 9700 South Cass Avenue, * Argonne, Illinois 60439. * * ************************************************************************ * SUBROUTINE DGEMM ( TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, $ BETA, C, LDC ) * .. Scalar Arguments .. CHARACTER*1 TRANSA, TRANSB INTEGER M, N, K, LDA, LDB, LDC DOUBLE PRECISION ALPHA, BETA * .. Array Arguments .. DOUBLE PRECISION A( LDA, * ), B( LDB, * ), C( LDC, * ) * .. * * Purpose * ======= * * DGEMM performs one of the matrix-matrix operations * * C := alpha*op( A )*op( B ) + beta*C, * * where op( X ) is one of * * op( X ) = X or op( X ) = X', * * alpha and beta are scalars, and A, B and C are matrices, with op( A ) * an m by k matrix, op( B ) a k by n matrix and C an m by n matrix. * * Parameters * ========== * * TRANSA - CHARACTER*1. * On entry, TRANSA specifies the form of op( A ) to be used in * the matrix multiplication as follows: * * TRANSA = 'N' or 'n', op( A ) = A. * * TRANSA = 'T' or 't', op( A ) = A'. * * TRANSA = 'C' or 'c', op( A ) = A'. * * Unchanged on exit. * * TRANSB - CHARACTER*1. * On entry, TRANSB specifies the form of op( B ) to be used in * the matrix multiplication as follows: * * TRANSB = 'N' or 'n', op( B ) = B. * * TRANSB = 'T' or 't', op( B ) = B'. * * TRANSB = 'C' or 'c', op( B ) = B'. * * Unchanged on exit. * * M - INTEGER. * On entry, M specifies the number of rows of the matrix * op( A ) and of the matrix C. M must be at least zero. * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the number of columns of the matrix * op( B ) and the number of columns of the matrix C. N must be * at least zero. * Unchanged on exit. * * K - INTEGER. * On entry, K specifies the number of columns of the matrix * op( A ) and the number of rows of the matrix op( B ). K 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, ka ), where ka is * k when TRANSA = 'N' or 'n', and is m otherwise. * Before entry with TRANSA = 'N' or 'n', the leading m by k * part of the array A must contain the matrix A, otherwise * the leading k by m part of the array A must contain the * matrix A. * Unchanged on exit. * * LDA - INTEGER. * On entry, LDA specifies the first dimension of A as declared * in the calling (sub) program. When TRANSA = 'N' or 'n' then * LDA must be at least max( 1, m ), otherwise LDA must be at * least max( 1, k ). * Unchanged on exit. * * B - DOUBLE PRECISION array of DIMENSION ( LDB, kb ), where kb is * n when TRANSB = 'N' or 'n', and is k otherwise. * Before entry with TRANSB = 'N' or 'n', the leading k by n * part of the array B must contain the matrix B, otherwise * the leading n by k part of the array B must contain the * matrix B. * Unchanged on exit. * * LDB - INTEGER. * On entry, LDB specifies the first dimension of B as declared * in the calling (sub) program. When TRANSB = 'N' or 'n' then * LDB must be at least max( 1, k ), otherwise LDB must be at * least max( 1, n ). * Unchanged on exit. * * BETA - DOUBLE PRECISION. * On entry, BETA specifies the scalar beta. When BETA is * supplied as zero then C need not be set on input. * Unchanged on exit. * * C - DOUBLE PRECISION array of DIMENSION ( LDC, n ). * Before entry, the leading m by n part of the array C must * contain the matrix C, except when beta is zero, in which * case C need not be set on entry. * On exit, the array C is overwritten by the m by n matrix * ( alpha*op( A )*op( B ) + beta*C ). * * LDC - INTEGER. * On entry, LDC specifies the first dimension of C as declared * in the calling (sub) program. LDC must be at least * max( 1, m ). * Unchanged on exit. * * * Level 3 Blas routine. * * -- Written on 8-February-1989. * Jack Dongarra, Argonne National Laboratory. * Iain Duff, AERE Harwell. * Jeremy Du Croz, Numerical Algorithms Group Ltd. * Sven Hammarling, Numerical Algorithms Group Ltd. * * * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. External Subroutines .. EXTERNAL XERBLA * .. Intrinsic Functions .. INTRINSIC MAX * .. Local Scalars .. LOGICAL NOTA, NOTB INTEGER I, INFO, J, L, NCOLA, NROWA, NROWB DOUBLE PRECISION TEMP * .. Parameters .. DOUBLE PRECISION ONE , ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) * .. * .. Executable Statements .. * * Set NOTA and NOTB as true if A and B respectively are not * transposed and set NROWA, NCOLA and NROWB as the number of rows * and columns of A and the number of rows of B respectively. * NOTA = LSAME( TRANSA, 'N' ) NOTB = LSAME( TRANSB, 'N' ) IF( NOTA )THEN NROWA = M NCOLA = K ELSE NROWA = K NCOLA = M END IF IF( NOTB )THEN NROWB = K ELSE NROWB = N END IF * * Test the input parameters. * INFO = 0 IF( ( .NOT.NOTA ).AND. $ ( .NOT.LSAME( TRANSA, 'C' ) ).AND. $ ( .NOT.LSAME( TRANSA, 'T' ) ) )THEN INFO = 1 ELSE IF( ( .NOT.NOTB ).AND. $ ( .NOT.LSAME( TRANSB, 'C' ) ).AND. $ ( .NOT.LSAME( TRANSB, 'T' ) ) )THEN INFO = 2 ELSE IF( M .LT.0 )THEN INFO = 3 ELSE IF( N .LT.0 )THEN INFO = 4 ELSE IF( K .LT.0 )THEN INFO = 5 ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN INFO = 8 ELSE IF( LDB.LT.MAX( 1, NROWB ) )THEN INFO = 10 ELSE IF( LDC.LT.MAX( 1, M ) )THEN INFO = 13 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'DGEMM ', INFO ) RETURN END IF * * Quick return if possible. * IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR. $ ( ( ( ALPHA.EQ.ZERO ).OR.( K.EQ.0 ) ).AND.( BETA.EQ.ONE ) ) ) $ RETURN * * And if alpha.eq.zero. * IF( ALPHA.EQ.ZERO )THEN IF( BETA.EQ.ZERO )THEN DO 20, J = 1, N DO 10, I = 1, M C( I, J ) = ZERO 10 CONTINUE 20 CONTINUE ELSE DO 40, J = 1, N DO 30, I = 1, M C( I, J ) = BETA*C( I, J ) 30 CONTINUE 40 CONTINUE END IF RETURN END IF * * Start the operations. * IF( NOTB )THEN IF( NOTA )THEN * * Form C := alpha*A*B + beta*C. * DO 90, J = 1, N IF( BETA.EQ.ZERO )THEN DO 50, I = 1, M C( I, J ) = ZERO 50 CONTINUE ELSE IF( BETA.NE.ONE )THEN DO 60, I = 1, M C( I, J ) = BETA*C( I, J ) 60 CONTINUE END IF DO 80, L = 1, K IF( B( L, J ).NE.ZERO )THEN TEMP = ALPHA*B( L, J ) DO 70, I = 1, M C( I, J ) = C( I, J ) + TEMP*A( I, L ) 70 CONTINUE END IF 80 CONTINUE 90 CONTINUE ELSE * * Form C := alpha*A'*B + beta*C * DO 120, J = 1, N DO 110, I = 1, M TEMP = ZERO DO 100, L = 1, K TEMP = TEMP + A( L, I )*B( L, J ) 100 CONTINUE IF( BETA.EQ.ZERO )THEN C( I, J ) = ALPHA*TEMP ELSE C( I, J ) = ALPHA*TEMP + BETA*C( I, J ) END IF 110 CONTINUE 120 CONTINUE END IF ELSE IF( NOTA )THEN * * Form C := alpha*A*B' + beta*C * DO 170, J = 1, N IF( BETA.EQ.ZERO )THEN DO 130, I = 1, M C( I, J ) = ZERO 130 CONTINUE ELSE IF( BETA.NE.ONE )THEN DO 140, I = 1, M C( I, J ) = BETA*C( I, J ) 140 CONTINUE END IF DO 160, L = 1, K IF( B( J, L ).NE.ZERO )THEN TEMP = ALPHA*B( J, L ) DO 150, I = 1, M C( I, J ) = C( I, J ) + TEMP*A( I, L ) 150 CONTINUE END IF 160 CONTINUE 170 CONTINUE ELSE * * Form C := alpha*A'*B' + beta*C * DO 200, J = 1, N DO 190, I = 1, M TEMP = ZERO DO 180, L = 1, K TEMP = TEMP + A( L, I )*B( J, L ) 180 CONTINUE IF( BETA.EQ.ZERO )THEN C( I, J ) = ALPHA*TEMP ELSE C( I, J ) = ALPHA*TEMP + BETA*C( I, J ) END IF 190 CONTINUE 200 CONTINUE END IF END IF * RETURN * * End of DGEMM . * END c c c SINGLE PRECISION ROUTINES FOLLOW c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c real function sasum(n,sx,incx) c c takes the sum of the absolute values. c uses unrolled loops for increment equal to one. c jack dongarra, linpack, 3/11/78. c modified to correct problem with negative increments, 9/29/88. c real sx(*),stemp integer i,ix,incx,m,mp1,n c sasum = 0.0e0 stemp = 0.0e0 if(n.le.0)return if(incx.eq.1)go to 20 c c code for increment not equal to 1 c ix = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 do 10 i = 1,n stemp = stemp + abs(sx(ix)) ix = ix + incx 10 continue sasum = stemp return c c code for increment equal to 1 c c c clean-up loop c 20 m = mod(n,6) if( m .eq. 0 ) go to 40 do 30 i = 1,m stemp = stemp + abs(sx(i)) 30 continue if( n .lt. 6 ) go to 60 40 mp1 = m + 1 do 50 i = mp1,n,6 stemp = stemp + abs(sx(i)) + abs(sx(i + 1)) + abs(sx(i + 2)) * + abs(sx(i + 3)) + abs(sx(i + 4)) + abs(sx(i + 5)) 50 continue 60 sasum = stemp return end c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c subroutine saxpy(n,sa,sx,incx,sy,incy) c c constant times a vector plus a vector. c uses unrolled loop for increments equal to one. c jack dongarra, linpack, 3/11/78. c real sx(*),sy(*),sa integer i,incx,incy,ix,iy,m,mp1,n c if(n.le.0)return if (sa .eq. 0.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 sy(iy) = sy(iy) + sa*sx(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 sy(i) = sy(i) + sa*sx(i) 30 continue if( n .lt. 4 ) return 40 mp1 = m + 1 do 50 i = mp1,n,4 sy(i) = sy(i) + sa*sx(i) sy(i + 1) = sy(i + 1) + sa*sx(i + 1) sy(i + 2) = sy(i + 2) + sa*sx(i + 2) sy(i + 3) = sy(i + 3) + sa*sx(i + 3) 50 continue return end c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c subroutine scopy(n,sx,incx,sy,incy) c c copies a vector, x, to a vector, y. c uses unrolled loops for increments equal to 1. c jack dongarra, linpack, 3/11/78. c real sx(*),sy(*) 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 sy(iy) = sx(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 sy(i) = sx(i) 30 continue if( n .lt. 7 ) return 40 mp1 = m + 1 do 50 i = mp1,n,7 sy(i) = sx(i) sy(i + 1) = sx(i + 1) sy(i + 2) = sx(i + 2) sy(i + 3) = sx(i + 3) sy(i + 4) = sx(i + 4) sy(i + 5) = sx(i + 5) sy(i + 6) = sx(i + 6) 50 continue return end c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c real function sdot(n,sx,incx,sy,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 real sx(*),sy(*),stemp integer i,incx,incy,ix,iy,m,mp1,n c stemp = 0.0e0 sdot = 0.0e0 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 stemp = stemp + sx(ix)*sy(iy) ix = ix + incx iy = iy + incy 10 continue sdot = stemp 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 stemp = stemp + sx(i)*sy(i) 30 continue if( n .lt. 5 ) go to 60 40 mp1 = m + 1 do 50 i = mp1,n,5 stemp = stemp + sx(i)*sy(i) + sx(i + 1)*sy(i + 1) + * sx(i + 2)*sy(i + 2) + sx(i + 3)*sy(i + 3) + sx(i + 4)*sy(i + 4) 50 continue 60 sdot = stemp return end c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c real function snrm2 ( n, sx, incx) integer incx, n, next real sx(*), cutlo, cuthi, hitest, sum, xmax, zero, one data zero, one /0.0e0, 1.0e0/ c c euclidean norm of the n-vector stored in sx() with storage c increment incx . c if n .le. 0 return with result = 0. c if n .ge. 1 then incx must be .ge. 1 c c c.l.lawson, 1978 jan 08 c c four phase method using two built-in constants that are c hopefully applicable to all machines. c cutlo = maximum of sqrt(u/eps) over all known machines. c cuthi = minimum of sqrt(v) over all known machines. c where c eps = smallest no. such that eps + 1. .gt. 1. c u = smallest positive no. (underflow limit) c v = largest no. (overflow limit) c c brief outline of algorithm.. c c phase 1 scans zero components. c move to phase 2 when a component is nonzero and .le. cutlo c move to phase 3 when a component is .gt. cutlo c move to phase 4 when a component is .ge. cuthi/m c where m = n for x() real and m = 2*n for complex. c c values for cutlo and cuthi.. c from the environmental parameters listed in the imsl converter c document the limiting values are as follows.. c cutlo, s.p. u/eps = 2**(-102) for honeywell. close seconds are c univac and dec at 2**(-103) c thus cutlo = 2**(-51) = 4.44089e-16 c cuthi, s.p. v = 2**127 for univac, honeywell, and dec. c thus cuthi = 2**(63.5) = 1.30438e19 c cutlo, d.p. u/eps = 2**(-67) for honeywell and dec. c thus cutlo = 2**(-33.5) = 8.23181d-11 c cuthi, d.p. same as s.p. cuthi = 1.30438d19 c data cutlo, cuthi / 8.232d-11, 1.304d19 / c data cutlo, cuthi / 4.441e-16, 1.304e19 / data cutlo, cuthi / 4.441e-16, 1.304e19 / c if(n .gt. 0) go to 10 snrm2 = zero go to 300 c 10 assign 30 to next sum = zero nn = n * incx c begin main loop i = 1 20 go to next,(30, 50, 70, 110) 30 if( abs(sx(i)) .gt. cutlo) go to 85 assign 50 to next xmax = zero c c phase 1. sum is zero c 50 if( sx(i) .eq. zero) go to 200 if( abs(sx(i)) .gt. cutlo) go to 85 c c prepare for phase 2. assign 70 to next go to 105 c c prepare for phase 4. c 100 i = j assign 110 to next sum = (sum / sx(i)) / sx(i) 105 xmax = abs(sx(i)) go to 115 c c phase 2. sum is small. c scale to avoid destructive underflow. c 70 if( abs(sx(i)) .gt. cutlo ) go to 75 c c common code for phases 2 and 4. c in phase 4 sum is large. scale to avoid overflow. c 110 if( abs(sx(i)) .le. xmax ) go to 115 sum = one + sum * (xmax / sx(i))**2 xmax = abs(sx(i)) go to 200 c 115 sum = sum + (sx(i)/xmax)**2 go to 200 c c c prepare for phase 3. c 75 sum = (sum * xmax) * xmax c c c for real or d.p. set hitest = cuthi/n c for complex set hitest = cuthi/(2*n) c 85 hitest = cuthi/float( n ) c c phase 3. sum is mid-range. no scaling. c do 95 j =i,nn,incx if(abs(sx(j)) .ge. hitest) go to 100 95 sum = sum + sx(j)**2 snrm2 = sqrt( sum ) go to 300 c 200 continue i = i + incx if ( i .le. nn ) go to 20 c c end of main loop. c c compute square root and adjust for scaling. c snrm2 = xmax * sqrt(sum) 300 continue return end c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c subroutine srot (n,sx,incx,sy,incy,c,s) c c applies a plane rotation. c jack dongarra, linpack, 3/11/78. c real sx(*),sy(*),stemp,c,s integer i,incx,incy,ix,iy,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 not equal c 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 stemp = c*sx(ix) + s*sy(iy) sy(iy) = c*sy(iy) - s*sx(ix) sx(ix) = stemp ix = ix + incx iy = iy + incy 10 continue return c c code for both increments equal to 1 c 20 do 30 i = 1,n stemp = c*sx(i) + s*sy(i) sy(i) = c*sy(i) - s*sx(i) sx(i) = stemp 30 continue return end c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c subroutine srotg(sa,sb,c,s) c c construct givens plane rotation. c jack dongarra, linpack, 3/11/78. c modified 9/27/86. c real sa,sb,c,s,roe,scale,r,z c roe = sb if( abs(sa) .gt. abs(sb) ) roe = sa scale = abs(sa) + abs(sb) if( scale .ne. 0.0 ) go to 10 c = 1.0 s = 0.0 r = 0.0 go to 20 10 r = scale*sqrt((sa/scale)**2 + (sb/scale)**2) r = sign(1.0,roe)*r c = sa/r s = sb/r 20 z = s if( abs(c) .gt. 0.0 .and. abs(c) .le. s ) z = 1.0/c sa = r sb = z return end c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c subroutine sswap (n,sx,incx,sy,incy) c c interchanges two vectors. c uses unrolled loops for increments equal to 1. c jack dongarra, linpack, 3/11/78. c real sx(*),sy(*),stemp 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 not equal c 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 stemp = sx(ix) sx(ix) = sy(iy) sy(iy) = stemp 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,3) if( m .eq. 0 ) go to 40 do 30 i = 1,m stemp = sx(i) sx(i) = sy(i) sy(i) = stemp 30 continue if( n .lt. 3 ) return 40 mp1 = m + 1 do 50 i = mp1,n,3 stemp = sx(i) sx(i) = sy(i) sy(i) = stemp stemp = sx(i + 1) sx(i + 1) = sy(i + 1) sy(i + 1) = stemp stemp = sx(i + 2) sx(i + 2) = sy(i + 2) sy(i + 2) = stemp 50 continue return end * ************************************************************************ * * File of the REAL Level-3 BLAS. * ========================================== * * SUBROUTINE SGEMM ( TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, * $ BETA, C, LDC ) * * SUBROUTINE SSYMM ( SIDE, UPLO, M, N, ALPHA, A, LDA, B, LDB, * $ BETA, C, LDC ) * * SUBROUTINE SSYRK ( UPLO, TRANS, N, K, ALPHA, A, LDA, * $ BETA, C, LDC ) * * SUBROUTINE SSYR2K( UPLO, TRANS, N, K, ALPHA, A, LDA, B, LDB, * $ BETA, C, LDC ) * * SUBROUTINE STRMM ( SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, * $ B, LDB ) * * SUBROUTINE STRSM ( SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, * $ B, LDB ) * * See: * * Dongarra J. J., Du Croz J. J., Duff I. and Hammarling S. * A set of Level 3 Basic Linear Algebra Subprograms. Technical * Memorandum No.88 (Revision 1), Mathematics and Computer Science * Division, Argonne National Laboratory, 9700 South Cass Avenue, * Argonne, Illinois 60439. * * ************************************************************************ * SUBROUTINE SGEMM ( TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, $ BETA, C, LDC ) * .. Scalar Arguments .. CHARACTER*1 TRANSA, TRANSB INTEGER M, N, K, LDA, LDB, LDC REAL ALPHA, BETA * .. Array Arguments .. REAL A( LDA, * ), B( LDB, * ), C( LDC, * ) * .. * * Purpose * ======= * * SGEMM performs one of the matrix-matrix operations * * C := alpha*op( A )*op( B ) + beta*C, * * where op( X ) is one of * * op( X ) = X or op( X ) = X', * * alpha and beta are scalars, and A, B and C are matrices, with op( A ) * an m by k matrix, op( B ) a k by n matrix and C an m by n matrix. * * Parameters * ========== * * TRANSA - CHARACTER*1. * On entry, TRANSA specifies the form of op( A ) to be used in * the matrix multiplication as follows: * * TRANSA = 'N' or 'n', op( A ) = A. * * TRANSA = 'T' or 't', op( A ) = A'. * * TRANSA = 'C' or 'c', op( A ) = A'. * * Unchanged on exit. * * TRANSB - CHARACTER*1. * On entry, TRANSB specifies the form of op( B ) to be used in * the matrix multiplication as follows: * * TRANSB = 'N' or 'n', op( B ) = B. * * TRANSB = 'T' or 't', op( B ) = B'. * * TRANSB = 'C' or 'c', op( B ) = B'. * * Unchanged on exit. * * M - INTEGER. * On entry, M specifies the number of rows of the matrix * op( A ) and of the matrix C. M must be at least zero. * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the number of columns of the matrix * op( B ) and the number of columns of the matrix C. N must be * at least zero. * Unchanged on exit. * * K - INTEGER. * On entry, K specifies the number of columns of the matrix * op( A ) and the number of rows of the matrix op( B ). K must * be at least zero. * Unchanged on exit. * * ALPHA - REAL . * On entry, ALPHA specifies the scalar alpha. * Unchanged on exit. * * A - REAL array of DIMENSION ( LDA, ka ), where ka is * k when TRANSA = 'N' or 'n', and is m otherwise. * Before entry with TRANSA = 'N' or 'n', the leading m by k * part of the array A must contain the matrix A, otherwise * the leading k by m part of the array A must contain the * matrix A. * Unchanged on exit. * * LDA - INTEGER. * On entry, LDA specifies the first dimension of A as declared * in the calling (sub) program. When TRANSA = 'N' or 'n' then * LDA must be at least max( 1, m ), otherwise LDA must be at * least max( 1, k ). * Unchanged on exit. * * B - REAL array of DIMENSION ( LDB, kb ), where kb is * n when TRANSB = 'N' or 'n', and is k otherwise. * Before entry with TRANSB = 'N' or 'n', the leading k by n * part of the array B must contain the matrix B, otherwise * the leading n by k part of the array B must contain the * matrix B. * Unchanged on exit. * * LDB - INTEGER. * On entry, LDB specifies the first dimension of B as declared * in the calling (sub) program. When TRANSB = 'N' or 'n' then * LDB must be at least max( 1, k ), otherwise LDB must be at * least max( 1, n ). * Unchanged on exit. * * BETA - REAL . * On entry, BETA specifies the scalar beta. When BETA is * supplied as zero then C need not be set on input. * Unchanged on exit. * * C - REAL array of DIMENSION ( LDC, n ). * Before entry, the leading m by n part of the array C must * contain the matrix C, except when beta is zero, in which * case C need not be set on entry. * On exit, the array C is overwritten by the m by n matrix * ( alpha*op( A )*op( B ) + beta*C ). * * LDC - INTEGER. * On entry, LDC specifies the first dimension of C as declared * in the calling (sub) program. LDC must be at least * max( 1, m ). * Unchanged on exit. * * * Level 3 Blas routine. * * -- Written on 8-February-1989. * Jack Dongarra, Argonne National Laboratory. * Iain Duff, AERE Harwell. * Jeremy Du Croz, Numerical Algorithms Group Ltd. * Sven Hammarling, Numerical Algorithms Group Ltd. * * * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. External Subroutines .. EXTERNAL XERBLA * .. Intrinsic Functions .. INTRINSIC MAX * .. Local Scalars .. LOGICAL NOTA, NOTB INTEGER I, INFO, J, L, NCOLA, NROWA, NROWB REAL TEMP * .. Parameters .. REAL ONE , ZERO PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) * .. * .. Executable Statements .. * * Set NOTA and NOTB as true if A and B respectively are not * transposed and set NROWA, NCOLA and NROWB as the number of rows * and columns of A and the number of rows of B respectively. * NOTA = LSAME( TRANSA, 'N' ) NOTB = LSAME( TRANSB, 'N' ) IF( NOTA )THEN NROWA = M NCOLA = K ELSE NROWA = K NCOLA = M END IF IF( NOTB )THEN NROWB = K ELSE NROWB = N END IF * * Test the input parameters. * INFO = 0 IF( ( .NOT.NOTA ).AND. $ ( .NOT.LSAME( TRANSA, 'C' ) ).AND. $ ( .NOT.LSAME( TRANSA, 'T' ) ) )THEN INFO = 1 ELSE IF( ( .NOT.NOTB ).AND. $ ( .NOT.LSAME( TRANSB, 'C' ) ).AND. $ ( .NOT.LSAME( TRANSB, 'T' ) ) )THEN INFO = 2 ELSE IF( M .LT.0 )THEN INFO = 3 ELSE IF( N .LT.0 )THEN INFO = 4 ELSE IF( K .LT.0 )THEN INFO = 5 ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN INFO = 8 ELSE IF( LDB.LT.MAX( 1, NROWB ) )THEN INFO = 10 ELSE IF( LDC.LT.MAX( 1, M ) )THEN INFO = 13 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'SGEMM ', INFO ) RETURN END IF * * Quick return if possible. * IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR. $ ( ( ( ALPHA.EQ.ZERO ).OR.( K.EQ.0 ) ).AND.( BETA.EQ.ONE ) ) ) $ RETURN * * And if alpha.eq.zero. * IF( ALPHA.EQ.ZERO )THEN IF( BETA.EQ.ZERO )THEN DO 20, J = 1, N DO 10, I = 1, M C( I, J ) = ZERO 10 CONTINUE 20 CONTINUE ELSE DO 40, J = 1, N DO 30, I = 1, M C( I, J ) = BETA*C( I, J ) 30 CONTINUE 40 CONTINUE END IF RETURN END IF * * Start the operations. * IF( NOTB )THEN IF( NOTA )THEN * * Form C := alpha*A*B + beta*C. * DO 90, J = 1, N IF( BETA.EQ.ZERO )THEN DO 50, I = 1, M C( I, J ) = ZERO 50 CONTINUE ELSE IF( BETA.NE.ONE )THEN DO 60, I = 1, M C( I, J ) = BETA*C( I, J ) 60 CONTINUE END IF DO 80, L = 1, K IF( B( L, J ).NE.ZERO )THEN TEMP = ALPHA*B( L, J ) DO 70, I = 1, M C( I, J ) = C( I, J ) + TEMP*A( I, L ) 70 CONTINUE END IF 80 CONTINUE 90 CONTINUE ELSE * * Form C := alpha*A'*B + beta*C * DO 120, J = 1, N DO 110, I = 1, M TEMP = ZERO DO 100, L = 1, K TEMP = TEMP + A( L, I )*B( L, J ) 100 CONTINUE IF( BETA.EQ.ZERO )THEN C( I, J ) = ALPHA*TEMP ELSE C( I, J ) = ALPHA*TEMP + BETA*C( I, J ) END IF 110 CONTINUE 120 CONTINUE END IF ELSE IF( NOTA )THEN * * Form C := alpha*A*B' + beta*C * DO 170, J = 1, N IF( BETA.EQ.ZERO )THEN DO 130, I = 1, M C( I, J ) = ZERO 130 CONTINUE ELSE IF( BETA.NE.ONE )THEN DO 140, I = 1, M C( I, J ) = BETA*C( I, J ) 140 CONTINUE END IF DO 160, L = 1, K IF( B( J, L ).NE.ZERO )THEN TEMP = ALPHA*B( J, L ) DO 150, I = 1, M C( I, J ) = C( I, J ) + TEMP*A( I, L ) 150 CONTINUE END IF 160 CONTINUE 170 CONTINUE ELSE * * Form C := alpha*A'*B' + beta*C * DO 200, J = 1, N DO 190, I = 1, M TEMP = ZERO DO 180, L = 1, K TEMP = TEMP + A( L, I )*B( J, L ) 180 CONTINUE IF( BETA.EQ.ZERO )THEN C( I, J ) = ALPHA*TEMP ELSE C( I, J ) = ALPHA*TEMP + BETA*C( I, J ) END IF 190 CONTINUE 200 CONTINUE END IF END IF * RETURN * * End of SGEMM . * END c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c integer function isamax(n,sx,incx) c c finds the index of element having max. absolute value. c jack dongarra, linpack, 3/11/78. c modified to correct problem with negative increments, 9/29/88. c real sx(*),smax integer i,incx,ix,n c isamax = 0 if( n .lt. 1 ) return isamax = 1 if(n.eq.1)return if(incx.eq.1)go to 20 c c code for increment not equal to 1 c ix = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 smax = abs(sx(ix)) ix = ix + incx do 10 i = 2,n if(abs(sx(ix)).le.smax) go to 5 isamax = i smax = abs(sx(ix)) 5 ix = ix + incx 10 continue return c c code for increment equal to 1 c 20 smax = abs(sx(1)) do 30 i = 2,n if(abs(sx(i)).le.smax) go to 30 isamax = i smax = abs(sx(i)) 30 continue return end c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c subroutine sscal(n,sa,sx,incx) c c scales a vector by a constant. c uses unrolled loops for increment equal to 1. c jack dongarra, linpack, 3/11/78. c modified to correct problem with negative increments, 9/29/88. c real sa,sx(*) integer i,ix,incx,m,mp1,n c if(n.le.0)return if(incx.eq.1)go to 20 c c code for increment not equal to 1 c ix = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 do 10 i = 1,n sx(ix) = sa*sx(ix) ix = ix + incx 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 sx(i) = sa*sx(i) 30 continue if( n .lt. 5 ) return 40 mp1 = m + 1 do 50 i = mp1,n,5 sx(i) = sa*sx(i) sx(i + 1) = sa*sx(i + 1) sx(i + 2) = sa*sx(i + 2) sx(i + 3) = sa*sx(i + 3) sx(i + 4) = sa*sx(i + 4) 50 continue return end * ************************************************************************ * SUBROUTINE STRSM ( SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, $ B, LDB ) * .. Scalar Arguments .. CHARACTER*1 SIDE, UPLO, TRANSA, DIAG INTEGER M, N, LDA, LDB REAL ALPHA * .. Array Arguments .. REAL A( LDA, * ), B( LDB, * ) * .. * * Purpose * ======= * * STRSM solves one of the matrix equations * * op( A )*X = alpha*B, or X*op( A ) = alpha*B, * * where alpha is a scalar, X and B are m by n matrices, A is a unit, or * non-unit, upper or lower triangular matrix and op( A ) is one of * * op( A ) = A or op( A ) = A'. * * The matrix X is overwritten on B. * * Parameters * ========== * * SIDE - CHARACTER*1. * On entry, SIDE specifies whether op( A ) appears on the left * or right of X as follows: * * SIDE = 'L' or 'l' op( A )*X = alpha*B. * * SIDE = 'R' or 'r' X*op( A ) = alpha*B. * * Unchanged on exit. * * UPLO - CHARACTER*1. * On entry, UPLO specifies whether the matrix A is an upper or * lower triangular matrix as follows: * * UPLO = 'U' or 'u' A is an upper triangular matrix. * * UPLO = 'L' or 'l' A is a lower triangular matrix. * * Unchanged on exit. * * TRANSA - CHARACTER*1. * On entry, TRANSA specifies the form of op( A ) to be used in * the matrix multiplication as follows: * * TRANSA = 'N' or 'n' op( A ) = A. * * TRANSA = 'T' or 't' op( A ) = A'. * * TRANSA = 'C' or 'c' op( A ) = A'. * * Unchanged on exit. * * DIAG - CHARACTER*1. * On entry, DIAG specifies whether or not A is unit triangular * as follows: * * DIAG = 'U' or 'u' A is assumed to be unit triangular. * * DIAG = 'N' or 'n' A is not assumed to be unit * triangular. * * Unchanged on exit. * * M - INTEGER. * On entry, M specifies the number of rows of B. M must be at * least zero. * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the number of columns of B. N must be * at least zero. * Unchanged on exit. * * ALPHA - REAL . * On entry, ALPHA specifies the scalar alpha. When alpha is * zero then A is not referenced and B need not be set before * entry. * Unchanged on exit. * * A - REAL array of DIMENSION ( LDA, k ), where k is m * when SIDE = 'L' or 'l' and is n when SIDE = 'R' or 'r'. * Before entry with UPLO = 'U' or 'u', the leading k by k * upper triangular part of the array A must contain the upper * triangular matrix and the strictly lower triangular part of * A is not referenced. * Before entry with UPLO = 'L' or 'l', the leading k by k * lower triangular part of the array A must contain the lower * triangular matrix and the strictly upper triangular part of * A is not referenced. * Note that when DIAG = 'U' or 'u', the diagonal elements of * A are not referenced either, but are assumed to be unity. * Unchanged on exit. * * LDA - INTEGER. * On entry, LDA specifies the first dimension of A as declared * in the calling (sub) program. When SIDE = 'L' or 'l' then * LDA must be at least max( 1, m ), when SIDE = 'R' or 'r' * then LDA must be at least max( 1, n ). * Unchanged on exit. * * B - REAL array of DIMENSION ( LDB, n ). * Before entry, the leading m by n part of the array B must * contain the right-hand side matrix B, and on exit is * overwritten by the solution matrix X. * * LDB - INTEGER. * On entry, LDB specifies the first dimension of B as declared * in the calling (sub) program. LDB must be at least * max( 1, m ). * Unchanged on exit. * * * Level 3 Blas routine. * * * -- Written on 8-February-1989. * Jack Dongarra, Argonne National Laboratory. * Iain Duff, AERE Harwell. * Jeremy Du Croz, Numerical Algorithms Group Ltd. * Sven Hammarling, Numerical Algorithms Group Ltd. * * * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. External Subroutines .. EXTERNAL XERBLA * .. Intrinsic Functions .. INTRINSIC MAX * .. Local Scalars .. LOGICAL LSIDE, NOUNIT, UPPER INTEGER I, INFO, J, K, NROWA REAL TEMP * .. Parameters .. REAL ONE , ZERO PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) * .. * .. Executable Statements .. * * Test the input parameters. * LSIDE = LSAME( SIDE , 'L' ) IF( LSIDE )THEN NROWA = M ELSE NROWA = N END IF NOUNIT = LSAME( DIAG , 'N' ) UPPER = LSAME( UPLO , 'U' ) * INFO = 0 IF( ( .NOT.LSIDE ).AND. $ ( .NOT.LSAME( SIDE , 'R' ) ) )THEN INFO = 1 ELSE IF( ( .NOT.UPPER ).AND. $ ( .NOT.LSAME( UPLO , 'L' ) ) )THEN INFO = 2 ELSE IF( ( .NOT.LSAME( TRANSA, 'N' ) ).AND. $ ( .NOT.LSAME( TRANSA, 'T' ) ).AND. $ ( .NOT.LSAME( TRANSA, 'C' ) ) )THEN INFO = 3 ELSE IF( ( .NOT.LSAME( DIAG , 'U' ) ).AND. $ ( .NOT.LSAME( DIAG , 'N' ) ) )THEN INFO = 4 ELSE IF( M .LT.0 )THEN INFO = 5 ELSE IF( N .LT.0 )THEN INFO = 6 ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN INFO = 9 ELSE IF( LDB.LT.MAX( 1, M ) )THEN INFO = 11 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'STRSM ', INFO ) RETURN END IF * * Quick return if possible. * IF( N.EQ.0 ) $ RETURN * * And when alpha.eq.zero. * IF( ALPHA.EQ.ZERO )THEN DO 20, J = 1, N DO 10, I = 1, M B( I, J ) = ZERO 10 CONTINUE 20 CONTINUE RETURN END IF * * Start the operations. * IF( LSIDE )THEN IF( LSAME( TRANSA, 'N' ) )THEN * * Form B := alpha*inv( A )*B. * IF( UPPER )THEN DO 60, J = 1, N IF( ALPHA.NE.ONE )THEN DO 30, I = 1, M B( I, J ) = ALPHA*B( I, J ) 30 CONTINUE END IF DO 50, K = M, 1, -1 IF( B( K, J ).NE.ZERO )THEN IF( NOUNIT ) $ B( K, J ) = B( K, J )/A( K, K ) DO 40, I = 1, K - 1 B( I, J ) = B( I, J ) - B( K, J )*A( I, K ) 40 CONTINUE END IF 50 CONTINUE 60 CONTINUE ELSE DO 100, J = 1, N IF( ALPHA.NE.ONE )THEN DO 70, I = 1, M B( I, J ) = ALPHA*B( I, J ) 70 CONTINUE END IF DO 90 K = 1, M IF( B( K, J ).NE.ZERO )THEN IF( NOUNIT ) $ B( K, J ) = B( K, J )/A( K, K ) DO 80, I = K + 1, M B( I, J ) = B( I, J ) - B( K, J )*A( I, K ) 80 CONTINUE END IF 90 CONTINUE 100 CONTINUE END IF ELSE * * Form B := alpha*inv( A' )*B. * IF( UPPER )THEN DO 130, J = 1, N DO 120, I = 1, M TEMP = ALPHA*B( I, J ) DO 110, K = 1, I - 1 TEMP = TEMP - A( K, I )*B( K, J ) 110 CONTINUE IF( NOUNIT ) $ TEMP = TEMP/A( I, I ) B( I, J ) = TEMP 120 CONTINUE 130 CONTINUE ELSE DO 160, J = 1, N DO 150, I = M, 1, -1 TEMP = ALPHA*B( I, J ) DO 140, K = I + 1, M TEMP = TEMP - A( K, I )*B( K, J ) 140 CONTINUE IF( NOUNIT ) $ TEMP = TEMP/A( I, I ) B( I, J ) = TEMP 150 CONTINUE 160 CONTINUE END IF END IF ELSE IF( LSAME( TRANSA, 'N' ) )THEN * * Form B := alpha*B*inv( A ). * IF( UPPER )THEN DO 210, J = 1, N IF( ALPHA.NE.ONE )THEN DO 170, I = 1, M B( I, J ) = ALPHA*B( I, J ) 170 CONTINUE END IF DO 190, K = 1, J - 1 IF( A( K, J ).NE.ZERO )THEN DO 180, I = 1, M B( I, J ) = B( I, J ) - A( K, J )*B( I, K ) 180 CONTINUE END IF 190 CONTINUE IF( NOUNIT )THEN TEMP = ONE/A( J, J ) DO 200, I = 1, M B( I, J ) = TEMP*B( I, J ) 200 CONTINUE END IF 210 CONTINUE ELSE DO 260, J = N, 1, -1 IF( ALPHA.NE.ONE )THEN DO 220, I = 1, M B( I, J ) = ALPHA*B( I, J ) 220 CONTINUE END IF DO 240, K = J + 1, N IF( A( K, J ).NE.ZERO )THEN DO 230, I = 1, M B( I, J ) = B( I, J ) - A( K, J )*B( I, K ) 230 CONTINUE END IF 240 CONTINUE IF( NOUNIT )THEN TEMP = ONE/A( J, J ) DO 250, I = 1, M B( I, J ) = TEMP*B( I, J ) 250 CONTINUE END IF 260 CONTINUE END IF ELSE * * Form B := alpha*B*inv( A' ). * IF( UPPER )THEN DO 310, K = N, 1, -1 IF( NOUNIT )THEN TEMP = ONE/A( K, K ) DO 270, I = 1, M B( I, K ) = TEMP*B( I, K ) 270 CONTINUE END IF DO 290, J = 1, K - 1 IF( A( J, K ).NE.ZERO )THEN TEMP = A( J, K ) DO 280, I = 1, M B( I, J ) = B( I, J ) - TEMP*B( I, K ) 280 CONTINUE END IF 290 CONTINUE IF( ALPHA.NE.ONE )THEN DO 300, I = 1, M B( I, K ) = ALPHA*B( I, K ) 300 CONTINUE END IF 310 CONTINUE ELSE DO 360, K = 1, N IF( NOUNIT )THEN TEMP = ONE/A( K, K ) DO 320, I = 1, M B( I, K ) = TEMP*B( I, K ) 320 CONTINUE END IF DO 340, J = K + 1, N IF( A( J, K ).NE.ZERO )THEN TEMP = A( J, K ) DO 330, I = 1, M B( I, J ) = B( I, J ) - TEMP*B( I, K ) 330 CONTINUE END IF 340 CONTINUE IF( ALPHA.NE.ONE )THEN DO 350, I = 1, M B( I, K ) = ALPHA*B( I, K ) 350 CONTINUE END IF 360 CONTINUE END IF END IF END IF * RETURN * * End of STRSM . * END * ************************************************************************ * SUBROUTINE STRSV ( UPLO, TRANS, DIAG, N, A, LDA, X, INCX ) * .. Scalar Arguments .. INTEGER INCX, LDA, N CHARACTER*1 DIAG, TRANS, UPLO * .. Array Arguments .. REAL A( LDA, * ), X( * ) * .. * * Purpose * ======= * * STRSV solves one of the systems of equations * * A*x = b, or A'*x = b, * * where b and x are n element vectors and A is an n by n unit, or * non-unit, upper or lower triangular matrix. * * No test for singularity or near-singularity is included in this * routine. Such tests must be performed before calling this routine. * * Parameters * ========== * * UPLO - CHARACTER*1. * On entry, UPLO specifies whether the matrix is an upper or * lower triangular matrix as follows: * * UPLO = 'U' or 'u' A is an upper triangular matrix. * * UPLO = 'L' or 'l' A is a lower triangular matrix. * * Unchanged on exit. * * TRANS - CHARACTER*1. * On entry, TRANS specifies the equations to be solved as * follows: * * TRANS = 'N' or 'n' A*x = b. * * TRANS = 'T' or 't' A'*x = b. * * TRANS = 'C' or 'c' A'*x = b. * * Unchanged on exit. * * DIAG - CHARACTER*1. * On entry, DIAG specifies whether or not A is unit * triangular as follows: * * DIAG = 'U' or 'u' A is assumed to be unit triangular. * * DIAG = 'N' or 'n' A is not assumed to be unit * triangular. * * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the order of the matrix A. * N must be at least zero. * Unchanged on exit. * * A - REAL array of DIMENSION ( LDA, n ). * Before entry with UPLO = 'U' or 'u', the leading n by n * upper triangular part of the array A must contain the upper * triangular matrix and the strictly lower triangular part of * A is not referenced. * Before entry with UPLO = 'L' or 'l', the leading n by n * lower triangular part of the array A must contain the lower * triangular matrix and the strictly upper triangular part of * A is not referenced. * Note that when DIAG = 'U' or 'u', the diagonal elements of * A are not referenced either, but are assumed to be unity. * 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, n ). * Unchanged on exit. * * X - REAL array of dimension at least * ( 1 + ( n - 1 )*abs( INCX ) ). * Before entry, the incremented array X must contain the n * element right-hand side vector b. On exit, X is overwritten * with the solution vector x. * * INCX - INTEGER. * On entry, INCX specifies the increment for the elements of * X. INCX 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 .. REAL ZERO PARAMETER ( ZERO = 0.0E+0 ) * .. Local Scalars .. REAL TEMP INTEGER I, INFO, IX, J, JX, KX LOGICAL NOUNIT * .. 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( UPLO , 'U' ).AND. $ .NOT.LSAME( UPLO , 'L' ) )THEN INFO = 1 ELSE IF( .NOT.LSAME( TRANS, 'N' ).AND. $ .NOT.LSAME( TRANS, 'T' ).AND. $ .NOT.LSAME( TRANS, 'C' ) )THEN INFO = 2 ELSE IF( .NOT.LSAME( DIAG , 'U' ).AND. $ .NOT.LSAME( DIAG , 'N' ) )THEN INFO = 3 ELSE IF( N.LT.0 )THEN INFO = 4 ELSE IF( LDA.LT.MAX( 1, N ) )THEN INFO = 6 ELSE IF( INCX.EQ.0 )THEN INFO = 8 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'STRSV ', INFO ) RETURN END IF * * Quick return if possible. * IF( N.EQ.0 ) $ RETURN * NOUNIT = LSAME( DIAG, 'N' ) * * Set up the start point in X if the increment is not unity. This * will be ( N - 1 )*INCX too small for descending loops. * IF( INCX.LE.0 )THEN KX = 1 - ( N - 1 )*INCX ELSE IF( INCX.NE.1 )THEN KX = 1 END IF * * Start the operations. In this version the elements of A are * accessed sequentially with one pass through A. * IF( LSAME( TRANS, 'N' ) )THEN * * Form x := inv( A )*x. * IF( LSAME( UPLO, 'U' ) )THEN IF( INCX.EQ.1 )THEN DO 20, J = N, 1, -1 IF( X( J ).NE.ZERO )THEN IF( NOUNIT ) $ X( J ) = X( J )/A( J, J ) TEMP = X( J ) DO 10, I = J - 1, 1, -1 X( I ) = X( I ) - TEMP*A( I, J ) 10 CONTINUE END IF 20 CONTINUE ELSE JX = KX + ( N - 1 )*INCX DO 40, J = N, 1, -1 IF( X( JX ).NE.ZERO )THEN IF( NOUNIT ) $ X( JX ) = X( JX )/A( J, J ) TEMP = X( JX ) IX = JX DO 30, I = J - 1, 1, -1 IX = IX - INCX X( IX ) = X( IX ) - TEMP*A( I, J ) 30 CONTINUE END IF JX = JX - INCX 40 CONTINUE END IF ELSE IF( INCX.EQ.1 )THEN DO 60, J = 1, N IF( X( J ).NE.ZERO )THEN IF( NOUNIT ) $ X( J ) = X( J )/A( J, J ) TEMP = X( J ) DO 50, I = J + 1, N X( I ) = X( I ) - TEMP*A( I, J ) 50 CONTINUE END IF 60 CONTINUE ELSE JX = KX DO 80, J = 1, N IF( X( JX ).NE.ZERO )THEN IF( NOUNIT ) $ X( JX ) = X( JX )/A( J, J ) TEMP = X( JX ) IX = JX DO 70, I = J + 1, N IX = IX + INCX X( IX ) = X( IX ) - TEMP*A( I, J ) 70 CONTINUE END IF JX = JX + INCX 80 CONTINUE END IF END IF ELSE * * Form x := inv( A' )*x. * IF( LSAME( UPLO, 'U' ) )THEN IF( INCX.EQ.1 )THEN DO 100, J = 1, N TEMP = X( J ) DO 90, I = 1, J - 1 TEMP = TEMP - A( I, J )*X( I ) 90 CONTINUE IF( NOUNIT ) $ TEMP = TEMP/A( J, J ) X( J ) = TEMP 100 CONTINUE ELSE JX = KX DO 120, J = 1, N TEMP = X( JX ) IX = KX DO 110, I = 1, J - 1 TEMP = TEMP - A( I, J )*X( IX ) IX = IX + INCX 110 CONTINUE IF( NOUNIT ) $ TEMP = TEMP/A( J, J ) X( JX ) = TEMP JX = JX + INCX 120 CONTINUE END IF ELSE IF( INCX.EQ.1 )THEN DO 140, J = N, 1, -1 TEMP = X( J ) DO 130, I = N, J + 1, -1 TEMP = TEMP - A( I, J )*X( I ) 130 CONTINUE IF( NOUNIT ) $ TEMP = TEMP/A( J, J ) X( J ) = TEMP 140 CONTINUE ELSE KX = KX + ( N - 1 )*INCX JX = KX DO 160, J = N, 1, -1 TEMP = X( JX ) IX = KX DO 150, I = N, J + 1, -1 TEMP = TEMP - A( I, J )*X( IX ) IX = IX - INCX 150 CONTINUE IF( NOUNIT ) $ TEMP = TEMP/A( J, J ) X( JX ) = TEMP JX = JX - INCX 160 CONTINUE END IF END IF END IF * RETURN * * End of STRSV . * END C*** ddemo.f c c FILE: ddemo.f c c==== ============================================================== c program ddemo c c This program reads in data, calls staircase subroutine DSTAIR, c calls pole placement subroutine DMEVAS, and calls back c transformation routine DBKTRN. c Computes eigenvalues of closed loop, and writes results. c c .. Parameters .. c implicit none integer Nin, Nout parameter (Nin = 5, Nout = 6) integer Nmax, Mmax parameter (Nmax = 40, Mmax = 40) integer lda, ldb, ldf parameter (lda=Nmax, ldb=Nmax, ldf=Mmax) c .... for upper bounds on Givens and Householder transformations c with N in {1,..,Nmax} and M in {1,..,min(N,Mmax)} c the expressions for gmax and hmax yield c greatest gmax = 211 when N=40, M=20 c greatest hmax = 401 when N=40, M=1 integer gmax parameter (gmax = 211) integer hmax parameter (hmax = 401) c .... and for the work space bounds integer liwork c c c Parameter liwork should normaly be declared as c c parameter (liwork = max(4*Nmax, Nmax+Nmax/2+gmax+hmax)) c c Microsoft's FORTRAN 5.00 compiler however reports a parameter c error that seems to be coming from the use of max. We have therefore c replace the declaration with the one below which is fine so far c as the test*.dat are concerned. All UNIX based FORTRAN c compilers had no problem with the above declaration. c parameter (liwork = Nmax + Nmax/2 + gmax + hmax) c integer lrwork parameter (lrwork = 3*Nmax + 2*gmax + 3*hmax) c c .. Local Scalars .. double precision tol integer n, m, kmax, ncmplx, iwarn, ierr, i, j character*20 header c c .. Local Arrays .. double precision A(lda,Nmax), B(ldb,Mmax), F(ldf,Nmax) double precision eigs(Nmax), rwork(lrwork) integer kstair(Nmax+1), info(2), iwork(liwork) integer itrnsf(nmax*(mmax+1)/2 + mmax+2*nmax+3) double precision rtrnsf(nmax*(mmax+1)/2 + nmax*(nmax+1)/2) double precision AA(lda,Nmax), BB(ldb,Mmax) double precision reigs(Nmax), imeigs(Nmax) c c .. External Subroutines .. external dstair, dmevas, dbktrn, lpeigs c c .. Executable Statements .. c c .. read the two headings in the data file c .. echo the second heading read (Nin,FMT=99990) header read (Nin,FMT=99990) header write (Nout,FMT=99999) write (Nout,FMT=99990) header c c .. read the data .. read (Nin,FMT=*) n, m, tol if (n.le.0 .or. n.gt.Nmax) then write (Nout,FMT=99998) n else read (Nin,FMT=*) (( A(i,j), j=1,n), i=1,n) if (m.le.0 .or. m.gt.Mmax) then write (Nout,FMT=99997) m else read (Nin,FMT=*) (( B(i,j), j=1,m), i=1,n) read (Nin,FMT=*) ( eigs(i), i=1,n ) read (Nin,FMT=*) ncmplx c c .. make copies of A,B so we can compute eigenvalues of closed loop c .. copy A to AA .. do 100 j = 1, n call dcopy(n, A(1,j), 1, AA(1,j), 1) 100 continue c .. copy B to BB .. do 120 j = 1, m call dcopy(n, B(1,j), 1, BB(1,j), 1) 120 continue c c .. echo the eigenvalues to be allocated write(Nout,FMT=80058) do 150 i = 1, ncmplx, 2 write(Nout,FMT=80054) EIGS(i),EIGS(i+1) write(Nout,FMT=80055) EIGS(i),EIGS(i+1) 150 continue do 170 i = ncmplx+1, n write(Nout,FMT=80056) EIGS(I) 170 continue c c ..compute the staircase form and the ranks of the c staircase blocks.. call dstair(n,m,A,lda,B,ldb, kmax, kstair, itrnsf, & rtrnsf, iwork, rwork, tol, iwarn, ierr) c if(ierr .lt. 0) then write(Nout,FMT=80000) -ierr else if (iwarn .ne. 0) then write (Nout,FMT=80020) iwarn end if c c .. allocate the eigenvalues .. call dmevas (n,m, ncmplx, gmax, hmax, A, lda, & B,ldb, F,ldf, eigs, kmax, kstair, & info, iwork, rwork, tol, iwarn, ierr) c write (Nout,FMT='()') if (ierr .lt. 0) then write(Nout,FMT=80000) -ierr else c .. print results .. if (iwarn .ne. 0) then write (Nout,FMT=80020) iwarn end if if (ierr .ne. 0) then write(Nout,FMT=80010) ierr end if write (Nout,FMT=80030) tol if (info(2) .ne. n) then write (Nout,FMT=80040) write (Nout,FMT=80041) info(2) end if if (info(1) .ne. n) then write (Nout,FMT=80050) n, info(1) c .. print UNallocated eigenvalues .. write (Nout,FMT=80052) do 200 i=info(1)+1,info(1)+ncmplx,2 write(Nout,FMT=80054) EIGS(i),EIGS(i+1) write(Nout,FMT=80055) EIGS(i),EIGS(i+1) 200 continue do 220 i = info(1)+1+ncmplx, n write(Nout,FMT=80056) EIGS(I) 220 continue end if c c ..do the back transform on F1 call dbktrn(n,m,F,ldf,itrnsf,rtrnsf,rwork,ierr) c c .. before printing F compute and print eigenvalues c of the closed loop. (lpeigs will overwrite AA).. call lpeigs(n,m, AA,lda, BB, ldb, F,ldf, & reigs, imeigs, iwork, rwork) c c .. print computed eigenvalues of closed loop .. c .. imaginary parts with magnitude < tol are set to zero .. write (Nout,FMT=80060) DO 400 i=1,n if ( abs(imeigs(i)) .LE. tol ) then write(Nout,FMT=80056) reigs(i) else if ( imeigs(i) .GE. 0.0 ) then write(Nout,FMT=80054) reigs(i),imeigs(i) else write(Nout,FMT=80055) reigs(i),-imeigs(i) endif 400 continue c c .. print computed F .. write (Nout,FMT='()') write (Nout,FMT=80080) DO 500 i = 1, m write(Nout,FMT=88888) (F(i,j), j=1,n) 500 continue c end if end if end if end if c 80000 FORMAT (' ERROR: error on ENTRY with argument ', I2) 80010 FORMAT (' ERROR: on EXIT ierr = ', I2) 80020 FORMAT (' WARNING: on exit iwarn = ', I1) 80030 FORMAT (' tolerance used = ', E16.8) 80040 FORMAT (' eigenvalue stored at EIGS(N) on entry ') 80041 FORMAT (' now stored at EIGS(', I2, ')') 80050 FORMAT (' of', I3, ' eigenvalues, the number allocated = ', I2) 80052 FORMAT (' the following eigenvalues were NOT allocated') 80054 FORMAT (F8.4, ' + i*', F8.4) 80055 FORMAT (F8.4, ' - i*', F8.4) 80056 FORMAT (F8.4) 80058 FORMAT (' the eigenvalues to be allocated are:') 80060 FORMAT (' the eigenvalues of the closed loop are:') 80080 FORMAT (' computed gain matrix F:') 88888 FORMAT (20(1x,F9.4)) 99990 FORMAT (A20) 99996 FORMAT (' kmax is out of range: kmax = ', I2) 99997 FORMAT (' m is out of range: m = ', I2) 99998 FORMAT (' n is out of range: n = ', I2) 99999 FORMAT (' Demonstration Program Results') c stop end c=================================================================== c subroutine lpeigs(n, m, A,lda, B,ldb, F,ldf, reig, imeig, & iwork, rwork) c c Purpose c ======= c To call routines to compute A-BF and the eigenvalues of A-BF. c c Arguments c ========= c Arguments In c ------------ c N INTEGER. c Row and column dimension of matrix A, c row dimension of matrix B, c column dimension of matrix F. c c M INTEGER. c Column dimension of matrix B, c row dimension of matrix F. c c A DOUBLE PRECISION array of DIMENSION (LDA,N). c The leading N by N part of this array must contain the matrix A. c Note: this array is overwritten. c c LDA INTEGER. c Row dimension of array A, as declared in the calling program c LDA .ge. N c c B DOUBLE PRECISION array of DIMENSION (LDB,M). c The leading N by M part of this array must contain the matrix B. c c LDB INTEGER. c Row dimension of array B, as declared in the calling program c LDB .ge. N. c c F DOUBLE PRECISION array of DIMENSION (LDF,N). c The leading M by N part of this array must contain the matrix F. c c LDF INTEGER. c Row dimension of array F, as declared in the calling program c LDB .ge. M. c c Arguments Out c ------------- c REIG DOUBLE PRECISION array of DIMENSION(N). c Contains the real parts of the computed eigenvalues. c c IMEIG DOUBLE PRECISION array of DIMENSION(N). c Contains the imaginary parts of the computed eigenvalues. c c Workspace c --------- c IWORK INTEGER array of DIMENSION(N). c c RWORK DOUBLE PRECISION array of DIMENSION(N). c c Tolerances c ---------- c None. c c Warning Indicator c ----------------- c None. c c Error Indicator c --------------- c None. c c Warnings and Errors Detected by the Routine c =========================================== c None c c Method c ====== c Uses BLAS routine DGEMM to compute B-AF. c Subsequent calls to EISPACK routines BALANC, ELMHES, HQR c balance the matrix, reduce it to upper hessenberg form, and c compute the eigenvalues via the QR algorithm. c c References c ========== c 1. Golub, G.H. and Van Loan, C.F., Matrix Computations, 2-nd ed., c Johns Hopkins University Press, Baltimore, 1989, Chapter 7. c c 2. Press, W.H. et al, Numerical Recipes, Cambridge University Press, c 1986, pp.365-376 c c c Revisions c ========= c 1994 Feb 03 c c arguments c implicit none integer n, m, lda, ldb, ldf, iwork(*) double precision A(lda,*), B(ldb,*), F(ldf,*) double precision reig(*), imeig(*), rwork(*) c c parameters character*1 Tran parameter(Tran='n') c c local variables integer low,igh,ierr c c ..compute closed loop A-B*F and store in A call dgemm(Tran,Tran, n,n,m, -1.0d0, B,ldb, F,ldf, 1.0d0, A,lda) c c ..compute eigenvalues of the closed loop (stored in A) call balanc( lda, n, A, low, igh, rwork) call elmhes( lda, n, low, igh, A, iwork) call hqr( lda, n, low, igh, A, reig, imeig, ierr) c return end C*** ddemo.sh # FILE: ddemo.sh # #!/bin/sh rm -f ddemo.log for i in test*.dat do echo $i ddemo.x < $i >> ddemo.log done C*** dlog.ref FILE: dlog.ref Demonstration Program Results test01 the eigenvalues to be allocated are: .9544 + i* .8513 .9544 - i* .8513 .2893 + i* .5374 .2893 - i* .5374 .5144 + i* .1034 .5144 - i* .1034 .4140 .5767 .8766 .4400 .7298 tolerance used = .98703268E-15 eigenvalue stored at EIGS(N) on entry now stored at EIGS( 3) the eigenvalues of the closed loop are: .9544 + i* .8513 .9544 - i* .8513 .2893 + i* .5374 .2893 - i* .5374 .7298 .8766 .5144 + i* .1034 .5144 - i* .1034 .5767 .4140 .4400 computed gain matrix F: 1.2366 -.0621 -1.1794 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .6030 -.1087 1.2359 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000 -.1042 -.8806 -.3792 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000 -.1902 .6348 -.5265 .1831 -.4240 .1296 -.3876 .2118 .1149 .0152 .0374 .5347 -.3346 .1998 -.0448 1.3815 -.1915 1.2044 -.1380 -.4250 .0377 -.0325 .1073 -.3782 .1560 .0393 .3933 -.8726 .6780 -.6574 -.6013 -.0181 -.1687 .1528 -.0455 -.4707 .0090 .5769 .3959 -.2292 .1545 .2531 -.2305 -.6368 Demonstration Program Results test02 the eigenvalues to be allocated are: .1312 + i* .8857 .1312 - i* .8857 .0922 + i* .1622 .0922 - i* .1622 .0711 + i* .3653 .0711 - i* .3653 .2531 .1351 .7832 .4553 .3495 tolerance used = .10610132E-14 eigenvalue stored at EIGS(N) on entry now stored at EIGS( 5) the eigenvalues of the closed loop are: .1312 + i* .8857 .1312 - i* .8857 .7832 .0711 + i* .3653 .0711 - i* .3653 .4553 .0922 + i* .1622 .0922 - i* .1622 .3495 .2531 .1351 computed gain matrix F: .1001 1.2141 -.0966 .0711 5.0344 .0000 .0000 .0000 .0000 .0000 .0000 -.4456 .3309 -.9465 -.4678 .0858 -1.3265 .5564 -.6209 -.0797 -.2924 .2009 -1.1914 -1.3564 2.8527 1.5547 1.6460 2.1313 -1.7744 1.2853 -.1253 .7265 -.4825 .6008 -.3938 -.2453 -1.5008 -.2922 -.2535 .2332 -.0025 .0174 .0085 -.2005 Demonstration Program Results test03 the eigenvalues to be allocated are: .7298 + i* .8693 .7298 - i* .8693 .7156 + i* .8007 .7156 - i* .8007 .7065 + i* .7417 .7065 - i* .7417 .0191 .8860 .5250 .4633 .0652 .7134 .4889 tolerance used = .11399781E-14 the eigenvalues of the closed loop are: .7298 + i* .8693 .7298 - i* .8693 .7156 + i* .8007 .7156 - i* .8007 .7065 + i* .7417 .7065 - i* .7417 .0191 .8860 .0652 .7134 .4633 .4889 .5250 computed gain matrix F: 3.8000 2.6996 1.8115 -12.1405 -1.2960 -12.3325 35.2843 -12.7917 -12.3098 6.7445 -46.2238 .0000 .0000 -.1299 2.2954 1.3242 -2.0195 -1.3181 -7.1153 12.0253 -6.3949 -5.2201 .1470 -22.7628 .0000 .0000 .1054 .7980 -.7706 .1107 .1728 .7289 -.1496 .0750 .0050 .5155 -.1578 .0200 -.1811 Demonstration Program Results test04 the eigenvalues to be allocated are: .1236 + i* .9734 .1236 - i* .9734 .0296 .0804 .4942 .7694 .9340 .2502 .3597 .7691 .5000 .7493 .6719 .6817 .7568 .0364 .2306 .2217 .5626 tolerance used = .19538777E-14 eigenvalue stored at EIGS(N) on entry now stored at EIGS( 3) the eigenvalues of the closed loop are: .1236 + i* .9733 .1236 - i* .9733 .9340 .0296 .0804 .3597 .4942 .5000 .6817 .7694 .6719 .7691 .5626 .7493 .2502 .0364 .7568 .2217 .2306 computed gain matrix F: .6844 -.7049 -.8840 .0815 .4088 .2177 -.9352 .6891 .3737 -.0426 .1916 .1751 -.5344 .3431 -.3536 -.0610 .1607 .0000 .0000 .5377 -1.2897 -.3693 -.0620 -1.1130 .6192 -.8361 -1.4656 1.1668 -1.0074 -.5434 .1141 -.8746 .4087 -1.5800 -.3427 -.2769 .0000 .0000 -.3414 -.1132 -1.1144 -1.0554 -1.0811 .8136 .9407 -.4220 1.7581 .6850 .9596 .4047 .0269 -1.2368 -1.1616 .8826 -.0732 .0000 .0000 -.2265 -.0043 -.2213 -1.1224 .3588 -.2094 .4005 .2942 .3240 .0002 -.2662 -.2661 .6789 .0902 .2643 .0241 -.1756 .0000 .0000 .2899 -.0107 .0301 .1390 .4738 -.7477 .6664 .2663 .5584 .8707 -.7539 -.1938 .2546 .1765 .4003 .2046 -.0310 .0000 .0000 .1126 -.0106 .1036 .4032 .6697 -.8356 .0606 .0181 -.2327 .3014 -.2265 .1679 -.1145 .0397 .0835 .0852 -.1335 -.0229 .1828 Demonstration Program Results test05 the eigenvalues to be allocated are: .4679 + i* .2872 .4679 - i* .2872 .1783 .1537 .5717 .8024 .0331 .5345 tolerance used = .69008470E-15 the eigenvalues of the closed loop are: .0331 .4679 + i* .2872 .4679 - i* .2872 .1537 .1783 .5344 .5717 .8024 computed gain matrix F: -.3808 -.3672 -.9907 -2.1210 3.0556 -1.1266 -.2946 -1.3675 .0657 1.5508 -.9710 -.4255 8.9037 -3.4898 .0508 .2359 .3084 .3142 .3307 1.2718 .4935 -1.1860 .2385 1.1073 -.8329 -.1707 -.6103 .2794 -.0432 .0798 -.4747 -.6511 Demonstration Program Results test06 the eigenvalues to be allocated are: .6216 + i* .8031 .6216 - i* .8031 .2478 .4764 .3893 .2033 .0284 .9017 .4265 tolerance used = .88104024E-15 eigenvalue stored at EIGS(N) on entry now stored at EIGS( 7) the eigenvalues of the closed loop are: .6216 + i* .8031 .6216 - i* .8031 .4764 .2478 .4265 .3893 .2033 .9017 .0284 computed gain matrix F: -2.5724 .6769 -1.2574 .3528 -.4997 -.9656 1.1845 .0000 .0000 -.2849 -.5452 -1.2494 .1274 -.5668 -.6412 -.3376 .0000 .0000 1.9985 .1808 1.2753 .5530 -.1198 -.8179 1.8712 .0000 .0000 -.6681 -1.0619 -.6613 .0059 -.1181 .3824 -.7930 .0000 .0000 -1.3325 -1.0313 -1.1417 -.1453 -.6545 -.2722 .0440 .0000 .0000 1.9865 .9305 -.2219 .2053 -.6571 -1.4549 1.0235 .0000 .0000 .2828 -.3858 .4358 .0333 -.4449 .3022 -.1848 .0000 .0000 .0465 -.0552 .7267 .2009 -.3151 .0733 .5488 -.0808 -.6580 Demonstration Program Results test07 the eigenvalues to be allocated are: .4679 + i* .2872 .4679 - i* .2872 .1783 + i* .1537 .1783 - i* .1537 .5717 .8024 .0331 .5345 tolerance used = .70267792E-15 the eigenvalues of the closed loop are: .4679 + i* .2872 .4679 - i* .2872 .0331 .1783 + i* .1537 .1783 - i* .1537 .8024 .5717 .5344 computed gain matrix F: .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .8879 .1303 .1784 -.5409 .0000 .0000 .0000 .0000 .6435 1.7507 -2.2184 .2100 1.3477 -6.4379 -2.4985 -23.3845 .0380 -.0819 2.0681 -.0602 -2.0816 7.7146 2.8935 26.6307 Demonstration Program Results test08 the eigenvalues to be allocated are: .5045 + i* .5163 .5045 - i* .5163 .3190 .9866 .4940 .2661 .0907 tolerance used = .49349413E-15 the eigenvalues of the closed loop are: .9866 .5045 + i* .5163 .5045 - i* .5163 .4940 .0907 .3190 .2661 computed gain matrix F: .2489 -2.3904 3.2020 -.0094 -.9292 1.9323 .0990 Demonstration Program Results test09 the eigenvalues to be allocated are: .3888 + i* .9522 .3888 - i* .9522 .9476 + i* .3898 .9476 - i* .3898 .2692 + i* .6922 .2692 - i* .6922 .2840 .7769 tolerance used = .78095530E-15 the eigenvalues of the closed loop are: .3888 + i* .9522 .3888 - i* .9522 .2692 + i* .6922 .2692 - i* .6922 .9475 + i* .3898 .9475 - i* .3898 .2840 .7769 computed gain matrix F: .3892 .0044 .5935 -1.0112 -.9976 -.1775 .3247 -1.2121 -1.4578 .0233 .2380 .5325 .8624 -.5668 -.1984 -.2923 -.6029 -.7545 -2.1290 .0062 1.1338 -.0397 -.7087 2.4381 .1572 -.1593 -.5716 1.6415 .6629 -.1088 -.1790 1.1188 .3124 -.2093 1.4273 -.8024 -.4689 .7710 .0748 -1.4526 .2782 .2652 1.8953 -.1284 -2.0392 .1586 -.0217 -1.6027 -.8854 .6067 -.7371 .3861 .3295 -.4005 .7691 .8716 -.7646 -.2133 -1.9246 .0416 .0268 .0591 -.3556 2.0463 Demonstration Program Results test10 the eigenvalues to be allocated are: .8287 + i* .0945 .8287 - i* .0945 .0817 + i* .7640 .0817 - i* .7640 .6296 + i* .2139 .6296 - i* .2139 .2136 .0811 tolerance used = .74874551E-15 the eigenvalues of the closed loop are: .0817 + i* .7640 .0817 - i* .7640 .8287 + i* .0945 .8287 - i* .0945 .6296 + i* .2139 .6296 - i* .2139 .0811 .2135 computed gain matrix F: .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000 -.0992 .7753 1.2865 1.2901 .0150 -.5475 .0000 .0000 .1659 -.4842 .6537 -.0132 -.1184 -.0933 .0000 .0000 .3834 -.3197 .2206 .2743 .1338 .1617 .0000 .0000 .2288 -.0223 -.0122 -.4309 .6747 .5071 .0000 .0000 .1007 -.2591 .0032 -.7122 -.4544 .2169 .0000 .0000 .2359 .2166 -.0051 -.2988 .1992 -.7912 .0000 .0000 -.3387 -.2734 -.0615 .3452 -.1993 -.2667 .3435 -.3311 Demonstration Program Results test11 the eigenvalues to be allocated are: .9017 + i* .4265 .9017 - i* .4265 .1420 + i* .9475 .1420 - i* .9475 .4103 + i* .1312 .4103 - i* .1312 .8857 .0922 WARNING: on exit iwarn = 1 tolerance used = .79822926E-15 the eigenvalues of the closed loop are: .1420 + i* .9475 .1420 - i* .9475 .9017 + i* .4265 .9017 - i* .4265 .8856 .0922 .4103 + i* .1312 .4103 - i* .1312 computed gain matrix F: .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .2648 -.9294 .1674 .2470 -.9065 -.4097 -.1178 .1701 .7043 .2260 -.1191 -.1199 -.2421 .2311 -.1223 .0333 -.5400 .5475 .6162 -.1538 .6763 -.6769 -.2359 .2566 .3207 .0781 -.7164 .1622 .4649 -.7413 .2581 -.0665 -.0691 .2863 -.5469 -.8528 .9486 1.2726 .3715 -.6028 -.1964 .4532 -.0715 -.3426 .2323 -.0024 .7386 .3585 .8049 .5425 -.1440 -.1333 -.2610 .1466 .6978 .1215 -1.0747 .0777 .0194 -.5196 .8482 -.3907 -.4436 -1.2912 Demonstration Program Results test12 the eigenvalues to be allocated are: .9017 + i* .4265 .9017 - i* .4265 .1420 .9475 .4103 .1312 .8857 .0922 WARNING: on exit iwarn = 1 tolerance used = .89167673E-15 the eigenvalues of the closed loop are: .9017 + i* .4265 .9017 - i* .4265 .9475 .1420 .8856 .4103 .1312 .0922 computed gain matrix F: .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .5939 -.4165 -.5444 -.7489 .0000 .0000 .0000 .0000 -.0214 1.1300 -.0313 .7734 .0000 .0000 .0000 .0000 -.1453 -.1448 1.2168 .4517 .0000 .0000 .0000 .0000 -.2979 -.0431 .3522 .7659 .0000 .0000 .0000 .0000 -.1032 -.0728 -.5010 -.2337 .5275 .3007 .9273 -.3155 -.3026 -.1632 -.3029 .2156 .1829 -.7744 .0493 .1539 Demonstration Program Results test13 the eigenvalues to be allocated are: .1234 + i* .4321 .1234 - i* .4321 .6789 + i* .9876 .6789 - i* .9876 .2468 + i* .8642 .2468 - i* .8642 WARNING: on exit iwarn = 1 WARNING: on exit iwarn = 2 ERROR: on EXIT ierr = 2 tolerance used = .37747583E-14 of 6 eigenvalues, the number allocated = 4 the following eigenvalues were NOT allocated .2468 + i* .8642 .2468 - i* .8642 the eigenvalues of the closed loop are: 1.0000 .1234 + i* .4321 .1234 - i* .4321 .6789 + i* .9876 .6789 - i* .9876 8.0000 computed gain matrix F: 1.3871 -.1798 -1.3884 -1.5106 -1.7360 .0000 -.0764 .3491 .5212 1.1550 1.2771 .0000 -2.1276 .2688 2.1960 1.8309 2.2450 .0000 Demonstration Program Results test14 the eigenvalues to be allocated are: .4645 + i* .9410 .4645 - i* .9410 .0501 .7615 .7702 .8278 .1254 WARNING: on exit iwarn = 1 WARNING: on exit iwarn = 2 ERROR: on EXIT ierr = 2 tolerance used = .53935489E-15 of 7 eigenvalues, the number allocated = 6 the following eigenvalues were NOT allocated .1254 the eigenvalues of the closed loop are: .4644 + i* .9410 .4644 - i* .9410 .0501 .8278 .7702 .7615 -.1122 computed gain matrix F: -1.0843 .0882 .2301 .2246 .0000 .0000 .2557 -.2033 1.2802 .0694 -.3091 -.5666 .0000 -.7177 -2.0619 -.1116 2.5711 -1.8260 2.1958 -9.5969 -.1839 Demonstration Program Results test15 the eigenvalues to be allocated are: .7219 + i* .4966 .7219 - i* .4966 .0537 .4416 .5192 .7719 .0654 .4428 WARNING: on exit iwarn = 1 WARNING: on exit iwarn = 2 ERROR: on EXIT ierr = 2 tolerance used = .64061423E-15 of 8 eigenvalues, the number allocated = 7 the following eigenvalues were NOT allocated .4428 the eigenvalues of the closed loop are: .7219 + i* .4966 .7219 - i* .4966 .0654 .0537 .7719 .5192 .4416 -.3030 computed gain matrix F: .3227 -.5004 -.4283 -.4781 .0902 -.3277 -.1166 .0000 -.0068 1.5612 -2.2274 -4.2036 .0341 -9.1466 -4.7626 .0000 .0732 -.3890 .1507 -1.1243 -.2800 2.2841 -.3359 .0000 -.3282 -.0839 .3472 .1730 -.0922 .2528 .0081 .0000 Demonstration Program Results test16 the eigenvalues to be allocated are: .7219 + i* .4966 .7219 - i* .4966 .0537 + i* .4416 .0537 - i* .4416 .5192 + i* .7719 .5192 - i* .7719 .0654 + i* .4428 .0654 - i* .4428 WARNING: on exit iwarn = 1 WARNING: on exit iwarn = 2 ERROR: on EXIT ierr = 2 tolerance used = .64061423E-15 of 8 eigenvalues, the number allocated = 6 the following eigenvalues were NOT allocated .0654 + i* .4428 .0654 - i* .4428 the eigenvalues of the closed loop are: 1.5847 .5192 + i* .7719 .5192 - i* .7719 .0537 + i* .4416 .0537 - i* .4416 .7219 + i* .4966 .7219 - i* .4966 -.3030 computed gain matrix F: -.2632 -.4941 -.3114 -.5365 .0647 -.5069 .0518 .0000 .2203 1.1585 -2.2792 -4.0142 -3.1412 -8.6576 -3.4432 .0000 -.3401 -.4836 -.9029 -.5833 -.3461 3.9367 -1.7233 .0000 .1311 -.0904 .2283 .2323 -.0662 .4350 -.1631 .0000 Demonstration Program Results test17 the eigenvalues to be allocated are: .1312 + i* .8857 .1312 - i* .8857 .0922 .1622 .0711 .3653 .2531 .1351 .7832 .4553 .3495 WARNING: on exit iwarn = 1 WARNING: on exit iwarn = 2 ERROR: on EXIT ierr = 2 tolerance used = .10610132E-14 of 11 eigenvalues, the number allocated = 8 the following eigenvalues were NOT allocated .7832 .4553 .3495 the eigenvalues of the closed loop are: .1312 + i* .8856 .1312 - i* .8856 .3653 .0711 .0922 .1351 .1622 .2531 .4872 .0551 + i* .3677 .0551 - i* .3677 computed gain matrix F: .2857 .9355 -2.4787 -1.0795 .2068 -1.9580 1.4284 -.8409 .1869 .0000 .0000 -.3546 .5065 .2770 .3595 .0899 -.1065 .0110 -.0425 .3711 .0000 .0000 -.7520 .5500 1.7342 1.5942 .7685 .5034 -.0882 -.2786 -.3631 .0000 .0000 .8251 -.4281 -.0957 -1.5130 -.1432 -.0325 -.0890 .0828 .2589 .0000 .0000 Demonstration Program Results test18 the eigenvalues to be allocated are: -.0219 + i* 1.9999 -.0219 - i* 1.9999 -.0865 + i* 1.9981 -.0865 - i* 1.9981 -.1910 + i* 1.9909 -.1910 - i* 1.9909 -.3309 + i* 1.9724 -.3309 - i* 1.9724 -.5000 + i* 1.9365 -.5000 - i* 1.9365 -.6910 + i* 1.8768 -.6910 - i* 1.8768 -.8955 + i* 1.7883 -.8955 - i* 1.7883 -1.1045 + i* 1.6673 -1.1045 - i* 1.6673 -1.3090 + i* 1.5121 -1.3090 - i* 1.5121 -1.5000 + i* 1.3229 -1.5000 - i* 1.3229 -1.6691 + i* 1.1018 -1.6691 - i* 1.1018 -1.8090 + i* .8529 -1.8090 - i* .8529 -1.9135 + i* .5817 -1.9135 - i* .5817 -1.9781 + i* .2948 -1.9781 - i* .2948 -2.0000 2.0000 tolerance used = .51625371E-13 the eigenvalues of the closed loop are: -.0219 + i* 1.9999 -.0219 - i* 1.9999 -.0865 + i* 1.9981 -.0865 - i* 1.9981 -.1910 + i* 1.9909 -.1910 - i* 1.9909 -.3309 + i* 1.9724 -.3309 - i* 1.9724 -.5000 + i* 1.9365 -.5000 - i* 1.9365 -.6910 + i* 1.8768 -.6910 - i* 1.8768 -.8955 + i* 1.7883 -.8955 - i* 1.7883 -1.1045 + i* 1.6673 -1.1045 - i* 1.6673 -1.3090 + i* 1.5121 -1.3090 - i* 1.5121 -1.5000 + i* 1.3229 -1.5000 - i* 1.3229 -1.6691 + i* 1.1018 -1.6691 - i* 1.1018 -1.8090 + i* .8529 -1.8090 - i* .8529 -1.9135 + i* .5817 -1.9135 - i* .5817 -1.9781 + i* .2948 -1.9781 - i* .2948 -2.0000 2.0000 computed gain matrix F: .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000 -.0001 -.0001 -.0001 .0000 .0002 .0004 .0007 .0001 .0015 .0053 .0192 .0647 .2241 .7693 Demonstration Program Results test19 the eigenvalues to be allocated are: -.0123 + i* 2.0000 -.0123 - i* 2.0000 -.0489 + i* 1.9994 -.0489 - i* 1.9994 -.1090 + i* 1.9970 -.1090 - i* 1.9970 -.1910 + i* 1.9909 -.1910 - i* 1.9909 -.2929 + i* 1.9784 -.2929 - i* 1.9784 -.4122 + i* 1.9571 -.4122 - i* 1.9571 -.5460 + i* 1.9240 -.5460 - i* 1.9240 -.6910 + i* 1.8768 -.6910 - i* 1.8768 -.8436 + i* 1.