| double precision function dasum(n,dx,incx) |
| c |
| c takes the sum of the absolute values. |
| c jack dongarra, linpack, 3/11/78. |
| c modified 3/93 to return if incx .le. 0. |
| c modified 12/3/93, array(1) declarations changed to array(*) |
| c |
| double precision dx(*),dtemp |
| integer i,incx,m,mp1,n,nincx |
| c |
| dasum = 0.0d0 |
| dtemp = 0.0d0 |
| if( n.le.0 .or. incx.le.0 )return |
| if(incx.eq.1)go to 20 |
| c |
| c code for increment not equal to 1 |
| c |
| nincx = n*incx |
| do 10 i = 1,nincx,incx |
| 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 |
| subroutine daxpy(n,da,dx,incx,dy,incy) |
| c |
| c constant times a vector plus a vector. |
| c uses unrolled loops for increments equal to one. |
| c jack dongarra, linpack, 3/11/78. |
| c modified 12/3/93, array(1) declarations changed to array(*) |
| c |
| double precision dx(*),dy(*),da |
| integer i,incx,incy,ix,iy,m,mp1,n |
| c |
| if(n.le.0)return |
| if (da .eq. 0.0d0) return |
| if(incx.eq.1.and.incy.eq.1)go to 20 |
| c |
| c code for unequal increments or equal increments |
| c not equal to 1 |
| c |
| ix = 1 |
| iy = 1 |
| if(incx.lt.0)ix = (-n+1)*incx + 1 |
| if(incy.lt.0)iy = (-n+1)*incy + 1 |
| do 10 i = 1,n |
| dy(iy) = dy(iy) + da*dx(ix) |
| ix = ix + incx |
| iy = iy + incy |
| 10 continue |
| return |
| c |
| c code for both increments equal to 1 |
| c |
| c |
| c clean-up loop |
| c |
| 20 m = mod(n,4) |
| if( m .eq. 0 ) go to 40 |
| do 30 i = 1,m |
| dy(i) = dy(i) + da*dx(i) |
| 30 continue |
| if( n .lt. 4 ) return |
| 40 mp1 = m + 1 |
| do 50 i = mp1,n,4 |
| dy(i) = dy(i) + da*dx(i) |
| dy(i + 1) = dy(i + 1) + da*dx(i + 1) |
| dy(i + 2) = dy(i + 2) + da*dx(i + 2) |
| dy(i + 3) = dy(i + 3) + da*dx(i + 3) |
| 50 continue |
| return |
| end |
| subroutine dcopy(n,dx,incx,dy,incy) |
| c |
| c copies a vector, x, to a vector, y. |
| c uses unrolled loops for increments equal to one. |
| c jack dongarra, linpack, 3/11/78. |
| c modified 12/3/93, array(1) declarations changed to array(*) |
| c |
| double precision dx(*),dy(*) |
| integer i,incx,incy,ix,iy,m,mp1,n |
| c |
| if(n.le.0)return |
| if(incx.eq.1.and.incy.eq.1)go to 20 |
| c |
| c code for unequal increments or equal increments |
| c not equal to 1 |
| c |
| ix = 1 |
| iy = 1 |
| if(incx.lt.0)ix = (-n+1)*incx + 1 |
| if(incy.lt.0)iy = (-n+1)*incy + 1 |
| do 10 i = 1,n |
| dy(iy) = dx(ix) |
| ix = ix + incx |
| iy = iy + incy |
| 10 continue |
| return |
| c |
| c code for both increments equal to 1 |
| c |
| c |
| c clean-up loop |
| c |
| 20 m = mod(n,7) |
| if( m .eq. 0 ) go to 40 |
| do 30 i = 1,m |
| dy(i) = dx(i) |
| 30 continue |
| if( n .lt. 7 ) return |
| 40 mp1 = m + 1 |
| do 50 i = mp1,n,7 |
| dy(i) = dx(i) |
| dy(i + 1) = dx(i + 1) |
| dy(i + 2) = dx(i + 2) |
| dy(i + 3) = dx(i + 3) |
| dy(i + 4) = dx(i + 4) |
| dy(i + 5) = dx(i + 5) |
| dy(i + 6) = dx(i + 6) |
| 50 continue |
| return |
| end |
| double precision function ddot(n,dx,incx,dy,incy) |
| c |
| c forms the dot product of two vectors. |
| c uses unrolled loops for increments equal to one. |
| c jack dongarra, linpack, 3/11/78. |
| c modified 12/3/93, array(1) declarations changed to array(*) |
| c |
| double precision dx(*),dy(*),dtemp |
| integer i,incx,incy,ix,iy,m,mp1,n |
| c |
| ddot = 0.0d0 |
| dtemp = 0.0d0 |
| if(n.le.0)return |
| if(incx.eq.1.and.incy.eq.1)go to 20 |
| c |
| c code for unequal increments or equal increments |
| c not equal to 1 |
| c |
| ix = 1 |
| iy = 1 |
| if(incx.lt.0)ix = (-n+1)*incx + 1 |
| if(incy.lt.0)iy = (-n+1)*incy + 1 |
| do 10 i = 1,n |
| dtemp = dtemp + dx(ix)*dy(iy) |
| ix = ix + incx |
| iy = iy + incy |
| 10 continue |
| ddot = dtemp |
| return |
| c |
| c code for both increments equal to 1 |
| c |
| c |
| c clean-up loop |
| c |
| 20 m = mod(n,5) |
| if( m .eq. 0 ) go to 40 |
| do 30 i = 1,m |
| dtemp = dtemp + dx(i)*dy(i) |
| 30 continue |
| if( n .lt. 5 ) go to 60 |
| 40 mp1 = m + 1 |
| do 50 i = mp1,n,5 |
| dtemp = dtemp + dx(i)*dy(i) + dx(i + 1)*dy(i + 1) + |
| * dx(i + 2)*dy(i + 2) + dx(i + 3)*dy(i + 3) + dx(i + 4)*dy(i + 4) |
| 50 continue |
| 60 ddot = dtemp |
| return |
| end |
| SUBROUTINE DGBMV ( TRANS, M, N, KL, KU, ALPHA, A, LDA, X, INCX, |
| $ BETA, Y, INCY ) |
| * .. Scalar Arguments .. |
| DOUBLE PRECISION ALPHA, BETA |
| INTEGER INCX, INCY, KL, KU, LDA, M, N |
| CHARACTER TRANS |
| * .. Array Arguments .. |
| DOUBLE PRECISION A( LDA, * ), X( * ), Y( * ) |
| * .. |
| * |
| * Purpose |
| * ======= |
| * |
| * DGBMV performs one of the matrix-vector operations |
| * |
| * y := alpha*A*x + beta*y, or y := alpha*A'*x + beta*y, |
| * |
| * where alpha and beta are scalars, x and y are vectors and A is an |
| * m by n band matrix, with kl sub-diagonals and ku super-diagonals. |
| * |
| * Parameters |
| * ========== |
| * |
| * TRANS - CHARACTER*1. |
| * On entry, TRANS specifies the operation to be performed as |
| * follows: |
| * |
| * TRANS = 'N' or 'n' y := alpha*A*x + beta*y. |
| * |
| * TRANS = 'T' or 't' y := alpha*A'*x + beta*y. |
| * |
| * TRANS = 'C' or 'c' y := alpha*A'*x + beta*y. |
| * |
| * Unchanged on exit. |
| * |
| * M - INTEGER. |
| * On entry, M specifies the number of rows of the matrix A. |
| * M must be at least zero. |
| * Unchanged on exit. |
| * |
| * N - INTEGER. |
| * On entry, N specifies the number of columns of the matrix A. |
| * N must be at least zero. |
| * Unchanged on exit. |
| * |
| * KL - INTEGER. |
| * On entry, KL specifies the number of sub-diagonals of the |
| * matrix A. KL must satisfy 0 .le. KL. |
| * Unchanged on exit. |
| * |
| * KU - INTEGER. |
| * On entry, KU specifies the number of super-diagonals of the |
| * matrix A. KU must satisfy 0 .le. KU. |
| * Unchanged on exit. |
| * |
| * ALPHA - DOUBLE PRECISION. |
| * On entry, ALPHA specifies the scalar alpha. |
| * Unchanged on exit. |
| * |
| * A - DOUBLE PRECISION array of DIMENSION ( LDA, n ). |
| * Before entry, the leading ( kl + ku + 1 ) by n part of the |
| * array A must contain the matrix of coefficients, supplied |
| * column by column, with the leading diagonal of the matrix in |
| * row ( ku + 1 ) of the array, the first super-diagonal |
| * starting at position 2 in row ku, the first sub-diagonal |
| * starting at position 1 in row ( ku + 2 ), and so on. |
| * Elements in the array A that do not correspond to elements |
| * in the band matrix (such as the top left ku by ku triangle) |
| * are not referenced. |
| * The following program segment will transfer a band matrix |
| * from conventional full matrix storage to band storage: |
| * |
| * DO 20, J = 1, N |
| * K = KU + 1 - J |
| * DO 10, I = MAX( 1, J - KU ), MIN( M, J + KL ) |
| * A( K + I, J ) = matrix( I, J ) |
| * 10 CONTINUE |
| * 20 CONTINUE |
| * |
| * 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 |
| * ( kl + ku + 1 ). |
| * Unchanged on exit. |
| * |
| * X - DOUBLE PRECISION array of DIMENSION at least |
| * ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n' |
| * and at least |
| * ( 1 + ( m - 1 )*abs( INCX ) ) otherwise. |
| * Before entry, the incremented array X must contain the |
| * vector x. |
| * Unchanged on exit. |
| * |
| * INCX - INTEGER. |
| * On entry, INCX specifies the increment for the elements of |
| * X. INCX must not be zero. |
| * Unchanged on exit. |
| * |
| * BETA - DOUBLE PRECISION. |
| * On entry, BETA specifies the scalar beta. When BETA is |
| * supplied as zero then Y need not be set on input. |
| * Unchanged on exit. |
| * |
| * Y - DOUBLE PRECISION array of DIMENSION at least |
| * ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n' |
| * and at least |
| * ( 1 + ( n - 1 )*abs( INCY ) ) otherwise. |
| * Before entry, the incremented array Y must contain the |
| * vector y. On exit, Y is overwritten by the updated vector y. |
| * |
| * INCY - INTEGER. |
| * On entry, INCY specifies the increment for the elements of |
| * Y. INCY must not be zero. |
| * Unchanged on exit. |
| * |
| * |
| * Level 2 Blas routine. |
| * |
| * -- Written on 22-October-1986. |
| * Jack Dongarra, Argonne National Lab. |
| * Jeremy Du Croz, Nag Central Office. |
| * Sven Hammarling, Nag Central Office. |
| * Richard Hanson, Sandia National Labs. |
| * |
| * .. Parameters .. |
| DOUBLE PRECISION ONE , ZERO |
| PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) |
| * .. Local Scalars .. |
| DOUBLE PRECISION TEMP |
| INTEGER I, INFO, IX, IY, J, JX, JY, K, KUP1, KX, KY, |
| $ LENX, LENY |
| * .. External Functions .. |
| LOGICAL LSAME |
| EXTERNAL LSAME |
| * .. External Subroutines .. |
| EXTERNAL XERBLA |
| * .. Intrinsic Functions .. |
| INTRINSIC MAX, MIN |
| * .. |
| * .. Executable Statements .. |
| * |
| * Test the input parameters. |
| * |
| INFO = 0 |
| IF ( .NOT.LSAME( TRANS, 'N' ).AND. |
| $ .NOT.LSAME( TRANS, 'T' ).AND. |
| $ .NOT.LSAME( TRANS, 'C' ) )THEN |
| INFO = 1 |
| ELSE IF( M.LT.0 )THEN |
| INFO = 2 |
| ELSE IF( N.LT.0 )THEN |
| INFO = 3 |
| ELSE IF( KL.LT.0 )THEN |
| INFO = 4 |
| ELSE IF( KU.LT.0 )THEN |
| INFO = 5 |
| ELSE IF( LDA.LT.( KL + KU + 1 ) )THEN |
| INFO = 8 |
| ELSE IF( INCX.EQ.0 )THEN |
| INFO = 10 |
| ELSE IF( INCY.EQ.0 )THEN |
| INFO = 13 |
| END IF |
| IF( INFO.NE.0 )THEN |
| CALL XERBLA( 'DGBMV ', INFO ) |
| RETURN |
| END IF |
| * |
| * Quick return if possible. |
| * |
| IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR. |
| $ ( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) ) |
| $ RETURN |
| * |
| * Set LENX and LENY, the lengths of the vectors x and y, and set |
| * up the start points in X and Y. |
| * |
| IF( LSAME( TRANS, 'N' ) )THEN |
| LENX = N |
| LENY = M |
| ELSE |
| LENX = M |
| LENY = N |
| END IF |
| IF( INCX.GT.0 )THEN |
| KX = 1 |
| ELSE |
| KX = 1 - ( LENX - 1 )*INCX |
| END IF |
| IF( INCY.GT.0 )THEN |
| KY = 1 |
| ELSE |
| KY = 1 - ( LENY - 1 )*INCY |
| END IF |
| * |
| * Start the operations. In this version the elements of A are |
| * accessed sequentially with one pass through the band part of A. |
| * |
| * First form y := beta*y. |
| * |
| IF( BETA.NE.ONE )THEN |
| IF( INCY.EQ.1 )THEN |
| IF( BETA.EQ.ZERO )THEN |
| DO 10, I = 1, LENY |
| Y( I ) = ZERO |
| 10 CONTINUE |
| ELSE |
| DO 20, I = 1, LENY |
| Y( I ) = BETA*Y( I ) |
| 20 CONTINUE |
| END IF |
| ELSE |
| IY = KY |
| IF( BETA.EQ.ZERO )THEN |
| DO 30, I = 1, LENY |
| Y( IY ) = ZERO |
| IY = IY + INCY |
| 30 CONTINUE |
| ELSE |
| DO 40, I = 1, LENY |
| Y( IY ) = BETA*Y( IY ) |
| IY = IY + INCY |
| 40 CONTINUE |
| END IF |
| END IF |
| END IF |
| IF( ALPHA.EQ.ZERO ) |
| $ RETURN |
| KUP1 = KU + 1 |
| IF( LSAME( TRANS, 'N' ) )THEN |
| * |
| * Form y := alpha*A*x + y. |
| * |
| JX = KX |
| IF( INCY.EQ.1 )THEN |
| DO 60, J = 1, N |
| c IF( X( JX ).NE.ZERO )THEN |
| TEMP = ALPHA*X( JX ) |
| K = KUP1 - J |
| DO 50, I = MAX( 1, J - KU ), MIN( M, J + KL ) |
| Y( I ) = Y( I ) + TEMP*A( K + I, J ) |
| 50 CONTINUE |
| c END IF |
| JX = JX + INCX |
| 60 CONTINUE |
| ELSE |
| DO 80, J = 1, N |
| c IF( X( JX ).NE.ZERO )THEN |
| TEMP = ALPHA*X( JX ) |
| IY = KY |
| K = KUP1 - J |
| DO 70, I = MAX( 1, J - KU ), MIN( M, J + KL ) |
| Y( IY ) = Y( IY ) + TEMP*A( K + I, J ) |
| IY = IY + INCY |
| 70 CONTINUE |
| c END IF |
| JX = JX + INCX |
| IF( J.GT.KU ) |
| $ KY = KY + INCY |
| 80 CONTINUE |
| END IF |
| ELSE |
| * |
| * Form y := alpha*A'*x + y. |
| * |
| JY = KY |
| IF( INCX.EQ.1 )THEN |
| DO 100, J = 1, N |
| TEMP = ZERO |
| K = KUP1 - J |
| DO 90, I = MAX( 1, J - KU ), MIN( M, J + KL ) |
| TEMP = TEMP + A( K + I, J )*X( I ) |
| 90 CONTINUE |
| Y( JY ) = Y( JY ) + ALPHA*TEMP |
| JY = JY + INCY |
| 100 CONTINUE |
| ELSE |
| DO 120, J = 1, N |
| TEMP = ZERO |
| IX = KX |
| K = KUP1 - J |
| DO 110, I = MAX( 1, J - KU ), MIN( M, J + KL ) |
| TEMP = TEMP + A( K + I, J )*X( IX ) |
| IX = IX + INCX |
| 110 CONTINUE |
| Y( JY ) = Y( JY ) + ALPHA*TEMP |
| JY = JY + INCY |
| IF( J.GT.KU ) |
| $ KX = KX + INCX |
| 120 CONTINUE |
| END IF |
| END IF |
| * |
| RETURN |
| * |
| * End of DGBMV . |
| * |
| END |
| SUBROUTINE DGEMM ( TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, |
| $ BETA, C, LDC ) |
| * .. Scalar Arguments .. |
| CHARACTER 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 |
| c 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 |
| c 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 |
| c 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 |
| c 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 |
| SUBROUTINE DGEMV ( TRANS, M, N, ALPHA, A, LDA, X, INCX, |
| $ BETA, Y, INCY ) |
| * .. Scalar Arguments .. |
| DOUBLE PRECISION ALPHA, BETA |
| INTEGER INCX, INCY, LDA, M, N |
| CHARACTER TRANS |
| * .. Array Arguments .. |
| DOUBLE PRECISION A( LDA, * ), X( * ), Y( * ) |
| * .. |
| * |
| * Purpose |
| * ======= |
| * |
| * DGEMV performs one of the matrix-vector operations |
| * |
| * y := alpha*A*x + beta*y, or y := alpha*A'*x + beta*y, |
| * |
| * where alpha and beta are scalars, x and y are vectors and A is an |
| * m by n matrix. |
| * |
| * Parameters |
| * ========== |
| * |
| * TRANS - CHARACTER*1. |
| * On entry, TRANS specifies the operation to be performed as |
| * follows: |
| * |
| * TRANS = 'N' or 'n' y := alpha*A*x + beta*y. |
| * |
| * TRANS = 'T' or 't' y := alpha*A'*x + beta*y. |
| * |
| * TRANS = 'C' or 'c' y := alpha*A'*x + beta*y. |
| * |
| * Unchanged on exit. |
| * |
| * M - INTEGER. |
| * On entry, M specifies the number of rows of the matrix A. |
| * M must be at least zero. |
| * Unchanged on exit. |
| * |
| * N - INTEGER. |
| * On entry, N specifies the number of columns of the matrix A. |
| * N must be at least zero. |
| * Unchanged on exit. |
| * |
| * ALPHA - DOUBLE PRECISION. |
| * On entry, ALPHA specifies the scalar alpha. |
| * Unchanged on exit. |
| * |
| * A - DOUBLE PRECISION array of DIMENSION ( LDA, n ). |
| * Before entry, the leading m by n part of the array A must |
| * contain the matrix of coefficients. |
| * Unchanged on exit. |
| * |
| * LDA - INTEGER. |
| * On entry, LDA specifies the first dimension of A as declared |
| * in the calling (sub) program. LDA must be at least |
| * max( 1, m ). |
| * Unchanged on exit. |
| * |
| * X - DOUBLE PRECISION array of DIMENSION at least |
| * ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n' |
| * and at least |
| * ( 1 + ( m - 1 )*abs( INCX ) ) otherwise. |
| * Before entry, the incremented array X must contain the |
| * vector x. |
| * Unchanged on exit. |
| * |
| * INCX - INTEGER. |
| * On entry, INCX specifies the increment for the elements of |
| * X. INCX must not be zero. |
| * Unchanged on exit. |
| * |
| * BETA - DOUBLE PRECISION. |
| * On entry, BETA specifies the scalar beta. When BETA is |
| * supplied as zero then Y need not be set on input. |
| * Unchanged on exit. |
| * |
| * Y - DOUBLE PRECISION array of DIMENSION at least |
| * ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n' |
| * and at least |
| * ( 1 + ( n - 1 )*abs( INCY ) ) otherwise. |
| * Before entry with BETA non-zero, the incremented array Y |
| * must contain the vector y. On exit, Y is overwritten by the |
| * updated vector y. |
| * |
| * INCY - INTEGER. |
| * On entry, INCY specifies the increment for the elements of |
| * Y. INCY must not be zero. |
| * Unchanged on exit. |
| * |
| * |
| * Level 2 Blas routine. |
| * |
| * -- Written on 22-October-1986. |
| * Jack Dongarra, Argonne National Lab. |
| * Jeremy Du Croz, Nag Central Office. |
| * Sven Hammarling, Nag Central Office. |
| * Richard Hanson, Sandia National Labs. |
| * |
| * |
| * .. Parameters .. |
| DOUBLE PRECISION ONE , ZERO |
| PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) |
| * .. Local Scalars .. |
| DOUBLE PRECISION TEMP |
| INTEGER I, INFO, IX, IY, J, JX, JY, KX, KY, LENX, LENY |
| * .. External Functions .. |
| LOGICAL LSAME |
| EXTERNAL LSAME |
| * .. External Subroutines .. |
| EXTERNAL XERBLA |
| * .. Intrinsic Functions .. |
| INTRINSIC MAX |
| * .. |
| * .. Executable Statements .. |
| * |
| * Test the input parameters. |
| * |
| INFO = 0 |
| IF ( .NOT.LSAME( TRANS, 'N' ).AND. |
| $ .NOT.LSAME( TRANS, 'T' ).AND. |
| $ .NOT.LSAME( TRANS, 'C' ) )THEN |
| INFO = 1 |
| ELSE IF( M.LT.0 )THEN |
| INFO = 2 |
| ELSE IF( N.LT.0 )THEN |
| INFO = 3 |
| ELSE IF( LDA.LT.MAX( 1, M ) )THEN |
| INFO = 6 |
| ELSE IF( INCX.EQ.0 )THEN |
| INFO = 8 |
| ELSE IF( INCY.EQ.0 )THEN |
| INFO = 11 |
| END IF |
| IF( INFO.NE.0 )THEN |
| CALL XERBLA( 'DGEMV ', INFO ) |
| RETURN |
| END IF |
| * |
| * Quick return if possible. |
| * |
| IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR. |
| $ ( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) ) |
| $ RETURN |
| * |
| * Set LENX and LENY, the lengths of the vectors x and y, and set |
| * up the start points in X and Y. |
| * |
| IF( LSAME( TRANS, 'N' ) )THEN |
| LENX = N |
| LENY = M |
| ELSE |
| LENX = M |
| LENY = N |
| END IF |
| IF( INCX.GT.0 )THEN |
| KX = 1 |
| ELSE |
| KX = 1 - ( LENX - 1 )*INCX |
| END IF |
| IF( INCY.GT.0 )THEN |
| KY = 1 |
| ELSE |
| KY = 1 - ( LENY - 1 )*INCY |
| END IF |
| * |
| * Start the operations. In this version the elements of A are |
| * accessed sequentially with one pass through A. |
| * |
| * First form y := beta*y. |
| * |
| IF( BETA.NE.ONE )THEN |
| IF( INCY.EQ.1 )THEN |
| IF( BETA.EQ.ZERO )THEN |
| DO 10, I = 1, LENY |
| Y( I ) = ZERO |
| 10 CONTINUE |
| ELSE |
| DO 20, I = 1, LENY |
| Y( I ) = BETA*Y( I ) |
| 20 CONTINUE |
| END IF |
| ELSE |
| IY = KY |
| IF( BETA.EQ.ZERO )THEN |
| DO 30, I = 1, LENY |
| Y( IY ) = ZERO |
| IY = IY + INCY |
| 30 CONTINUE |
| ELSE |
| DO 40, I = 1, LENY |
| Y( IY ) = BETA*Y( IY ) |
| IY = IY + INCY |
| 40 CONTINUE |
| END IF |
| END IF |
| END IF |
| IF( ALPHA.EQ.ZERO ) |
| $ RETURN |
| IF( LSAME( TRANS, 'N' ) )THEN |
| * |
| * Form y := alpha*A*x + y. |
| * |
| JX = KX |
| IF( INCY.EQ.1 )THEN |
| DO 60, J = 1, N |
| c IF( X( JX ).NE.ZERO )THEN |
| TEMP = ALPHA*X( JX ) |
| DO 50, I = 1, M |
| Y( I ) = Y( I ) + TEMP*A( I, J ) |
| 50 CONTINUE |
| c END IF |
| JX = JX + INCX |
| 60 CONTINUE |
| ELSE |
| DO 80, J = 1, N |
| c IF( X( JX ).NE.ZERO )THEN |
| TEMP = ALPHA*X( JX ) |
| IY = KY |
| DO 70, I = 1, M |
| Y( IY ) = Y( IY ) + TEMP*A( I, J ) |
| IY = IY + INCY |
| 70 CONTINUE |
| c END IF |
| JX = JX + INCX |
| 80 CONTINUE |
| END IF |
| ELSE |
| * |
| * Form y := alpha*A'*x + y. |
| * |
| JY = KY |
| IF( INCX.EQ.1 )THEN |
| DO 100, J = 1, N |
| TEMP = ZERO |
| DO 90, I = 1, M |
| TEMP = TEMP + A( I, J )*X( I ) |
| 90 CONTINUE |
| Y( JY ) = Y( JY ) + ALPHA*TEMP |
| JY = JY + INCY |
| 100 CONTINUE |
| ELSE |
| DO 120, J = 1, N |
| TEMP = ZERO |
| IX = KX |
| DO 110, I = 1, M |
| TEMP = TEMP + A( I, J )*X( IX ) |
| IX = IX + INCX |
| 110 CONTINUE |
| Y( JY ) = Y( JY ) + ALPHA*TEMP |
| JY = JY + INCY |
| 120 CONTINUE |
| END IF |
| END IF |
| * |
| RETURN |
| * |
| * End of DGEMV . |
| * |
| END |
| SUBROUTINE DGER ( M, N, ALPHA, X, INCX, Y, INCY, A, LDA ) |
| * .. Scalar Arguments .. |
| DOUBLE PRECISION ALPHA |
| INTEGER INCX, INCY, LDA, M, N |
| * .. Array Arguments .. |
| DOUBLE PRECISION A( LDA, * ), X( * ), Y( * ) |
| * .. |
| * |
| * Purpose |
| * ======= |
| * |
| * DGER performs the rank 1 operation |
| * |
| * A := alpha*x*y' + A, |
| * |
| * where alpha is a scalar, x is an m element vector, y is an n element |
| * vector and A is an m by n matrix. |
| * |
| * Parameters |
| * ========== |
| * |
| * M - INTEGER. |
| * On entry, M specifies the number of rows of the matrix A. |
| * M must be at least zero. |
| * Unchanged on exit. |
| * |
| * N - INTEGER. |
| * On entry, N specifies the number of columns of the matrix A. |
| * N must be at least zero. |
| * Unchanged on exit. |
| * |
| * ALPHA - DOUBLE PRECISION. |
| * On entry, ALPHA specifies the scalar alpha. |
| * Unchanged on exit. |
| * |
| * X - DOUBLE PRECISION array of dimension at least |
| * ( 1 + ( m - 1 )*abs( INCX ) ). |
| * Before entry, the incremented array X must contain the m |
| * element vector x. |
| * Unchanged on exit. |
| * |
| * INCX - INTEGER. |
| * On entry, INCX specifies the increment for the elements of |
| * X. INCX must not be zero. |
| * Unchanged on exit. |
| * |
| * Y - DOUBLE PRECISION array of dimension at least |
| * ( 1 + ( n - 1 )*abs( INCY ) ). |
| * Before entry, the incremented array Y must contain the n |
| * element vector y. |
| * Unchanged on exit. |
| * |
| * INCY - INTEGER. |
| * On entry, INCY specifies the increment for the elements of |
| * Y. INCY must not be zero. |
| * Unchanged on exit. |
| * |
| * A - DOUBLE PRECISION array of DIMENSION ( LDA, n ). |
| * Before entry, the leading m by n part of the array A must |
| * contain the matrix of coefficients. On exit, A is |
| * overwritten by the updated matrix. |
| * |
| * LDA - INTEGER. |
| * On entry, LDA specifies the first dimension of A as declared |
| * in the calling (sub) program. LDA must be at least |
| * max( 1, m ). |
| * Unchanged on exit. |
| * |
| * |
| * Level 2 Blas routine. |
| * |
| * -- Written on 22-October-1986. |
| * Jack Dongarra, Argonne National Lab. |
| * Jeremy Du Croz, Nag Central Office. |
| * Sven Hammarling, Nag Central Office. |
| * Richard Hanson, Sandia National Labs. |
| * |
| * |
| * .. Parameters .. |
| DOUBLE PRECISION ZERO |
| PARAMETER ( ZERO = 0.0D+0 ) |
| * .. Local Scalars .. |
| DOUBLE PRECISION TEMP |
| INTEGER I, INFO, IX, J, JY, KX |
| * .. External Subroutines .. |
| EXTERNAL XERBLA |
| * .. Intrinsic Functions .. |
| INTRINSIC MAX |
| * .. |
| * .. Executable Statements .. |
| * |
| * Test the input parameters. |
| * |
| INFO = 0 |
| IF ( M.LT.0 )THEN |
| INFO = 1 |
| ELSE IF( N.LT.0 )THEN |
| INFO = 2 |
| ELSE IF( INCX.EQ.0 )THEN |
| INFO = 5 |
| ELSE IF( INCY.EQ.0 )THEN |
| INFO = 7 |
| ELSE IF( LDA.LT.MAX( 1, M ) )THEN |
| INFO = 9 |
| END IF |
| IF( INFO.NE.0 )THEN |
| CALL XERBLA( 'DGER ', INFO ) |
| RETURN |
| END IF |
| * |
| * Quick return if possible. |
| * |
| IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR.( ALPHA.EQ.ZERO ) ) |
| $ RETURN |
| * |
| * Start the operations. In this version the elements of A are |
| * accessed sequentially with one pass through A. |
| * |
| IF( INCY.GT.0 )THEN |
| JY = 1 |
| ELSE |
| JY = 1 - ( N - 1 )*INCY |
| END IF |
| IF( INCX.EQ.1 )THEN |
| DO 20, J = 1, N |
| IF( Y( JY ).NE.ZERO )THEN |
| TEMP = ALPHA*Y( JY ) |
| DO 10, I = 1, M |
| A( I, J ) = A( I, J ) + X( I )*TEMP |
| 10 CONTINUE |
| END IF |
| JY = JY + INCY |
| 20 CONTINUE |
| ELSE |
| IF( INCX.GT.0 )THEN |
| KX = 1 |
| ELSE |
| KX = 1 - ( M - 1 )*INCX |
| END IF |
| DO 40, J = 1, N |
| IF( Y( JY ).NE.ZERO )THEN |
| TEMP = ALPHA*Y( JY ) |
| IX = KX |
| DO 30, I = 1, M |
| A( I, J ) = A( I, J ) + X( IX )*TEMP |
| IX = IX + INCX |
| 30 CONTINUE |
| END IF |
| JY = JY + INCY |
| 40 CONTINUE |
| END IF |
| * |
| RETURN |
| * |
| * End of DGER . |
| * |
| END |
| DOUBLE PRECISION FUNCTION DNRM2 ( N, X, INCX ) |
| * .. Scalar Arguments .. |
| INTEGER INCX, N |
| * .. Array Arguments .. |
| DOUBLE PRECISION X( * ) |
| * .. |
| * |
| * DNRM2 returns the euclidean norm of a vector via the function |
| * name, so that |
| * |
| * DNRM2 := sqrt( x'*x ) |
| * |
| * |
| * |
| * -- This version written on 25-October-1982. |
| * Modified on 14-October-1993 to inline the call to DLASSQ. |
| * Sven Hammarling, Nag Ltd. |
| * |
| * |
| * .. Parameters .. |
| DOUBLE PRECISION ONE , ZERO |
| PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) |
| * .. Local Scalars .. |
| INTEGER IX |
| DOUBLE PRECISION ABSXI, NORM, SCALE, SSQ |
| * .. Intrinsic Functions .. |
| INTRINSIC ABS, SQRT |
| * .. |
| * .. Executable Statements .. |
| IF( N.LT.1 .OR. INCX.LT.1 )THEN |
| NORM = ZERO |
| ELSE IF( N.EQ.1 )THEN |
| NORM = ABS( X( 1 ) ) |
| ELSE |
| SCALE = ZERO |
| SSQ = ONE |
| * The following loop is equivalent to this call to the LAPACK |
| * auxiliary routine: |
| * CALL DLASSQ( N, X, INCX, SCALE, SSQ ) |
| * |
| DO 10, IX = 1, 1 + ( N - 1 )*INCX, INCX |
| IF( X( IX ).NE.ZERO )THEN |
| ABSXI = ABS( X( IX ) ) |
| IF( SCALE.LT.ABSXI )THEN |
| SSQ = ONE + SSQ*( SCALE/ABSXI )**2 |
| SCALE = ABSXI |
| ELSE |
| SSQ = SSQ + ( ABSXI/SCALE )**2 |
| END IF |
| END IF |
| 10 CONTINUE |
| NORM = SCALE * SQRT( SSQ ) |
| END IF |
| * |
| DNRM2 = NORM |
| RETURN |
| * |
| * End of DNRM2. |
| * |
| END |
| subroutine drot (n,dx,incx,dy,incy,c,s) |
| c |
| c applies a plane rotation. |
| c jack dongarra, linpack, 3/11/78. |
| c modified 12/3/93, array(1) declarations changed to array(*) |
| 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 |
| subroutine drotg(da,db,c,s) |
| c |
| c construct givens plane rotation. |
| c jack dongarra, linpack, 3/11/78. |
| 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 |
| z = 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 |
| z = 1.0d0 |
| if( dabs(da) .gt. dabs(db) ) z = s |
| if( dabs(db) .ge. dabs(da) .and. c .ne. 0.0d0 ) z = 1.0d0/c |
| 20 da = r |
| db = z |
| return |
| end |
| |
| SUBROUTINE DROTM (N,DX,INCX,DY,INCY,DPARAM) |
| C |
| C APPLY THE MODIFIED GIVENS TRANSFORMATION, H, TO THE 2 BY N MATRIX |
| C |
| C (DX**T) , WHERE **T INDICATES TRANSPOSE. THE ELEMENTS OF DX ARE IN |
| C (DY**T) |
| C |
| C DX(LX+I*INCX), I = 0 TO N-1, WHERE LX = 1 IF INCX .GE. 0, ELSE |
| C LX = (-INCX)*N, AND SIMILARLY FOR SY USING LY AND INCY. |
| C WITH DPARAM(1)=DFLAG, H HAS ONE OF THE FOLLOWING FORMS.. |
| C |
| C DFLAG=-1.D0 DFLAG=0.D0 DFLAG=1.D0 DFLAG=-2.D0 |
| C |
| C (DH11 DH12) (1.D0 DH12) (DH11 1.D0) (1.D0 0.D0) |
| C H=( ) ( ) ( ) ( ) |
| C (DH21 DH22), (DH21 1.D0), (-1.D0 DH22), (0.D0 1.D0). |
| C SEE DROTMG FOR A DESCRIPTION OF DATA STORAGE IN DPARAM. |
| C |
| INTEGER N, INCX, INCY |
| INTEGER NSTEPS, I, KX,KY |
| DOUBLE PRECISION DFLAG,DH12,DH22,DX,TWO,Z,DH11,DH21, |
| 1 DPARAM,DY,W,ZERO |
| DIMENSION DX(*),DY(*),DPARAM(5) |
| DATA ZERO,TWO/0.D0,2.D0/ |
| C |
| DFLAG=DPARAM(1) |
| IF(N .LE. 0 .OR.(DFLAG+TWO.EQ.ZERO)) GO TO 140 |
| IF(.NOT.(INCX.EQ.INCY.AND. INCX .GT.0)) GO TO 70 |
| C |
| NSTEPS=N*INCX |
| C IF(DFLAG) 50,10,30 |
| IF(DFLAG .GT. 0.0) GOTO 50 |
| IF(DFLAG .LT. 0.0) GOTO 30 |
| DH12=DPARAM(4) |
| DH21=DPARAM(3) |
| DO 20 I=1,NSTEPS,INCX |
| W=DX(I) |
| Z=DY(I) |
| DX(I)=W+Z*DH12 |
| DY(I)=W*DH21+Z |
| 20 CONTINUE |
| GO TO 140 |
| 30 CONTINUE |
| DH11=DPARAM(2) |
| DH22=DPARAM(5) |
| DO 40 I=1,NSTEPS,INCX |
| W=DX(I) |
| Z=DY(I) |
| DX(I)=W*DH11+Z |
| DY(I)=-W+DH22*Z |
| 40 CONTINUE |
| GO TO 140 |
| 50 CONTINUE |
| DH11=DPARAM(2) |
| DH12=DPARAM(4) |
| DH21=DPARAM(3) |
| DH22=DPARAM(5) |
| DO 60 I=1,NSTEPS,INCX |
| W=DX(I) |
| Z=DY(I) |
| DX(I)=W*DH11+Z*DH12 |
| DY(I)=W*DH21+Z*DH22 |
| 60 CONTINUE |
| GO TO 140 |
| 70 CONTINUE |
| KX=1 |
| KY=1 |
| IF(INCX .LT. 0) KX=1+(1-N)*INCX |
| IF(INCY .LT. 0) KY=1+(1-N)*INCY |
| C |
| C IF(DFLAG)120,80,100 |
| IF(DFLAG. GT. ZERO) GOTO 120 |
| IF(DFLAG. LT. ZERO) GOTO 100 |
| DH12=DPARAM(4) |
| DH21=DPARAM(3) |
| DO 90 I=1,N |
| W=DX(KX) |
| Z=DY(KY) |
| DX(KX)=W+Z*DH12 |
| DY(KY)=W*DH21+Z |
| KX=KX+INCX |
| KY=KY+INCY |
| 90 CONTINUE |
| GO TO 140 |
| 100 CONTINUE |
| DH11=DPARAM(2) |
| DH22=DPARAM(5) |
| DO 110 I=1,N |
| W=DX(KX) |
| Z=DY(KY) |
| DX(KX)=W*DH11+Z |
| DY(KY)=-W+DH22*Z |
| KX=KX+INCX |
| KY=KY+INCY |
| 110 CONTINUE |
| GO TO 140 |
| 120 CONTINUE |
| DH11=DPARAM(2) |
| DH12=DPARAM(4) |
| DH21=DPARAM(3) |
| DH22=DPARAM(5) |
| DO 130 I=1,N |
| W=DX(KX) |
| Z=DY(KY) |
| DX(KX)=W*DH11+Z*DH12 |
| DY(KY)=W*DH21+Z*DH22 |
| KX=KX+INCX |
| KY=KY+INCY |
| 130 CONTINUE |
| 140 CONTINUE |
| RETURN |
| END |
| SUBROUTINE DROTMG (DD1,DD2,DX1,DY1,DPARAM) |
| C |
| C CONSTRUCT THE MODIFIED GIVENS TRANSFORMATION MATRIX H WHICH ZEROS |
| C THE SECOND COMPONENT OF THE 2-VECTOR (DSQRT(DD1)*DX1,DSQRT(DD2)* |
| C DY2)**T. |
| C WITH DPARAM(1)=DFLAG, H HAS ONE OF THE FOLLOWING FORMS.. |
| C |
| C DFLAG=-1.D0 DFLAG=0.D0 DFLAG=1.D0 DFLAG=-2.D0 |
| C |
| C (DH11 DH12) (1.D0 DH12) (DH11 1.D0) (1.D0 0.D0) |
| C H=( ) ( ) ( ) ( ) |
| C (DH21 DH22), (DH21 1.D0), (-1.D0 DH22), (0.D0 1.D0). |
| C LOCATIONS 2-4 OF DPARAM CONTAIN DH11, DH21, DH12, AND DH22 |
| C RESPECTIVELY. (VALUES OF 1.D0, -1.D0, OR 0.D0 IMPLIED BY THE |
| C VALUE OF DPARAM(1) ARE NOT STORED IN DPARAM.) |
| C |
| C THE VALUES OF GAMSQ AND RGAMSQ SET IN THE DATA STATEMENT MAY BE |
| C INEXACT. THIS IS OK AS THEY ARE ONLY USED FOR TESTING THE SIZE |
| C OF DD1 AND DD2. ALL ACTUAL SCALING OF DATA IS DONE USING GAM. |
| C |
| DOUBLE PRECISION GAM,ONE,RGAMSQ,DD2,DH11,DH21,DPARAM,DP2, |
| 1 DQ2,DU,DY1,ZERO,GAMSQ,DD1,DFLAG,DH12,DH22,DP1,DQ1, |
| 2 DTEMP,DX1,TWO |
| DIMENSION DPARAM(5) |
| C |
| INTEGER IGO |
| C |
| DATA ZERO,ONE,TWO /0.D0,1.D0,2.D0/ |
| DATA GAM,GAMSQ,RGAMSQ/4096.D0,16777216.D0,5.9604645D-8/ |
| C |
| |
| IF(.NOT. DD1 .LT. ZERO) GO TO 10 |
| C GO ZERO-H-D-AND-DX1.. |
| GO TO 60 |
| 10 CONTINUE |
| C CASE-DD1-NONNEGATIVE |
| DP2=DD2*DY1 |
| IF(.NOT. DP2 .EQ. ZERO) GO TO 20 |
| DFLAG=-TWO |
| GO TO 260 |
| C REGULAR-CASE.. |
| 20 CONTINUE |
| DP1=DD1*DX1 |
| DQ2=DP2*DY1 |
| DQ1=DP1*DX1 |
| C |
| IF(.NOT. DABS(DQ1) .GT. DABS(DQ2)) GO TO 40 |
| DH21=-DY1/DX1 |
| DH12=DP2/DP1 |
| C |
| DU=ONE-DH12*DH21 |
| C |
| IF(.NOT. DU .LE. ZERO) GO TO 30 |
| C GO ZERO-H-D-AND-DX1.. |
| GO TO 60 |
| 30 CONTINUE |
| DFLAG=ZERO |
| DD1=DD1/DU |
| DD2=DD2/DU |
| DX1=DX1*DU |
| C GO SCALE-CHECK.. |
| GO TO 100 |
| 40 CONTINUE |
| IF(.NOT. DQ2 .LT. ZERO) GO TO 50 |
| C GO ZERO-H-D-AND-DX1.. |
| GO TO 60 |
| 50 CONTINUE |
| DFLAG=ONE |
| DH11=DP1/DP2 |
| DH22=DX1/DY1 |
| DU=ONE+DH11*DH22 |
| DTEMP=DD2/DU |
| DD2=DD1/DU |
| DD1=DTEMP |
| DX1=DY1*DU |
| C GO SCALE-CHECK |
| GO TO 100 |
| C PROCEDURE..ZERO-H-D-AND-DX1.. |
| 60 CONTINUE |
| DFLAG=-ONE |
| DH11=ZERO |
| DH12=ZERO |
| DH21=ZERO |
| DH22=ZERO |
| C |
| DD1=ZERO |
| DD2=ZERO |
| DX1=ZERO |
| C RETURN.. |
| GO TO 220 |
| C PROCEDURE..FIX-H.. |
| 70 CONTINUE |
| IF(.NOT. DFLAG .GE. ZERO) GO TO 90 |
| C |
| IF(.NOT. DFLAG .EQ. ZERO) GO TO 80 |
| DH11=ONE |
| DH22=ONE |
| DFLAG=-ONE |
| GO TO 90 |
| 80 CONTINUE |
| DH21=-ONE |
| DH12=ONE |
| DFLAG=-ONE |
| 90 CONTINUE |
| IF (IGO .EQ. 120) GOTO 120 |
| IF (IGO .EQ. 150) GOTO 150 |
| IF (IGO .EQ. 180) GOTO 180 |
| IF (IGO .EQ. 210) GOTO 210 |
| C PROCEDURE..SCALE-CHECK |
| 100 CONTINUE |
| 110 CONTINUE |
| IF(.NOT. DD1 .LE. RGAMSQ) GO TO 130 |
| IF(DD1 .EQ. ZERO) GO TO 160 |
| IGO = 120 |
| C FIX-H.. |
| GO TO 70 |
| 120 CONTINUE |
| DD1=DD1*GAM**2 |
| DX1=DX1/GAM |
| DH11=DH11/GAM |
| DH12=DH12/GAM |
| GO TO 110 |
| 130 CONTINUE |
| 140 CONTINUE |
| IF(.NOT. DD1 .GE. GAMSQ) GO TO 160 |
| IGO = 150 |
| C FIX-H.. |
| GO TO 70 |
| 150 CONTINUE |
| DD1=DD1/GAM**2 |
| DX1=DX1*GAM |
| DH11=DH11*GAM |
| DH12=DH12*GAM |
| GO TO 140 |
| 160 CONTINUE |
| 170 CONTINUE |
| IF(.NOT. DABS(DD2) .LE. RGAMSQ) GO TO 190 |
| IF(DD2 .EQ. ZERO) GO TO 220 |
| IGO = 180 |
| C FIX-H.. |
| GO TO 70 |
| 180 CONTINUE |
| DD2=DD2*GAM**2 |
| DH21=DH21/GAM |
| DH22=DH22/GAM |
| GO TO 170 |
| 190 CONTINUE |
| 200 CONTINUE |
| IF(.NOT. DABS(DD2) .GE. GAMSQ) GO TO 220 |
| IGO = 210 |
| C FIX-H.. |
| GO TO 70 |
| 210 CONTINUE |
| DD2=DD2/GAM**2 |
| DH21=DH21*GAM |
| DH22=DH22*GAM |
| GO TO 200 |
| 220 CONTINUE |
| C IF(DFLAG)250,230,240 |
| IF(DFLAG .GT. ZERO) GOTO 250 |
| IF(DFLAG .LT. ZERO) GOTO 240 |
| DPARAM(3)=DH21 |
| DPARAM(4)=DH12 |
| GO TO 260 |
| 240 CONTINUE |
| DPARAM(2)=DH11 |
| DPARAM(5)=DH22 |
| GO TO 260 |
| 250 CONTINUE |
| DPARAM(2)=DH11 |
| DPARAM(3)=DH21 |
| DPARAM(4)=DH12 |
| DPARAM(5)=DH22 |
| 260 CONTINUE |
| DPARAM(1)=DFLAG |
| RETURN |
| END |
| SUBROUTINE DSBMV ( UPLO, N, K, ALPHA, A, LDA, X, INCX, |
| $ BETA, Y, INCY ) |
| * .. Scalar Arguments .. |
| DOUBLE PRECISION ALPHA, BETA |
| INTEGER INCX, INCY, K, LDA, N |
| CHARACTER UPLO |
| * .. Array Arguments .. |
| DOUBLE PRECISION A( LDA, * ), X( * ), Y( * ) |
| * .. |
| * |
| * Purpose |
| * ======= |
| * |
| * DSBMV performs the matrix-vector operation |
| * |
| * y := alpha*A*x + beta*y, |
| * |
| * where alpha and beta are scalars, x and y are n element vectors and |
| * A is an n by n symmetric band matrix, with k super-diagonals. |
| * |
| * Parameters |
| * ========== |
| * |
| * UPLO - CHARACTER*1. |
| * On entry, UPLO specifies whether the upper or lower |
| * triangular part of the band matrix A is being supplied as |
| * follows: |
| * |
| * UPLO = 'U' or 'u' The upper triangular part of A is |
| * being supplied. |
| * |
| * UPLO = 'L' or 'l' The lower triangular part of A is |
| * being supplied. |
| * |
| * Unchanged on exit. |
| * |
| * N - INTEGER. |
| * On entry, N specifies the order of the matrix A. |
| * N must be at least zero. |
| * Unchanged on exit. |
| * |
| * K - INTEGER. |
| * On entry, K specifies the number of super-diagonals of the |
| * matrix A. K must satisfy 0 .le. K. |
| * Unchanged on exit. |
| * |
| * ALPHA - DOUBLE PRECISION. |
| * On entry, ALPHA specifies the scalar alpha. |
| * Unchanged on exit. |
| * |
| * A - DOUBLE PRECISION array of DIMENSION ( LDA, n ). |
| * Before entry with UPLO = 'U' or 'u', the leading ( k + 1 ) |
| * by n part of the array A must contain the upper triangular |
| * band part of the symmetric matrix, supplied column by |
| * column, with the leading diagonal of the matrix in row |
| * ( k + 1 ) of the array, the first super-diagonal starting at |
| * position 2 in row k, and so on. The top left k by k triangle |
| * of the array A is not referenced. |
| * The following program segment will transfer the upper |
| * triangular part of a symmetric band matrix from conventional |
| * full matrix storage to band storage: |
| * |
| * DO 20, J = 1, N |
| * M = K + 1 - J |
| * DO 10, I = MAX( 1, J - K ), J |
| * A( M + I, J ) = matrix( I, J ) |
| * 10 CONTINUE |
| * 20 CONTINUE |
| * |
| * Before entry with UPLO = 'L' or 'l', the leading ( k + 1 ) |
| * by n part of the array A must contain the lower triangular |
| * band part of the symmetric matrix, supplied column by |
| * column, with the leading diagonal of the matrix in row 1 of |
| * the array, the first sub-diagonal starting at position 1 in |
| * row 2, and so on. The bottom right k by k triangle of the |
| * array A is not referenced. |
| * The following program segment will transfer the lower |
| * triangular part of a symmetric band matrix from conventional |
| * full matrix storage to band storage: |
| * |
| * DO 20, J = 1, N |
| * M = 1 - J |
| * DO 10, I = J, MIN( N, J + K ) |
| * A( M + I, J ) = matrix( I, J ) |
| * 10 CONTINUE |
| * 20 CONTINUE |
| * |
| * 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 |
| * ( k + 1 ). |
| * 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 |
| * vector x. |
| * Unchanged on exit. |
| * |
| * INCX - INTEGER. |
| * On entry, INCX specifies the increment for the elements of |
| * X. INCX must not be zero. |
| * Unchanged on exit. |
| * |
| * BETA - DOUBLE PRECISION. |
| * On entry, BETA specifies the scalar beta. |
| * Unchanged on exit. |
| * |
| * Y - DOUBLE PRECISION array of DIMENSION at least |
| * ( 1 + ( n - 1 )*abs( INCY ) ). |
| * Before entry, the incremented array Y must contain the |
| * vector y. On exit, Y is overwritten by the updated vector y. |
| * |
| * INCY - INTEGER. |
| * On entry, INCY specifies the increment for the elements of |
| * Y. INCY must not be zero. |
| * Unchanged on exit. |
| * |
| * |
| * Level 2 Blas routine. |
| * |
| * -- Written on 22-October-1986. |
| * Jack Dongarra, Argonne National Lab. |
| * Jeremy Du Croz, Nag Central Office. |
| * Sven Hammarling, Nag Central Office. |
| * Richard Hanson, Sandia National Labs. |
| * |
| * |
| * .. Parameters .. |
| DOUBLE PRECISION ONE , ZERO |
| PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) |
| * .. Local Scalars .. |
| DOUBLE PRECISION TEMP1, TEMP2 |
| INTEGER I, INFO, IX, IY, J, JX, JY, KPLUS1, KX, KY, L |
| * .. External Functions .. |
| LOGICAL LSAME |
| EXTERNAL LSAME |
| * .. External Subroutines .. |
| EXTERNAL XERBLA |
| * .. Intrinsic Functions .. |
| INTRINSIC MAX, MIN |
| * .. |
| * .. Executable Statements .. |
| * |
| * Test the input parameters. |
| * |
| INFO = 0 |
| IF ( .NOT.LSAME( UPLO, 'U' ).AND. |
| $ .NOT.LSAME( UPLO, 'L' ) )THEN |
| INFO = 1 |
| ELSE IF( N.LT.0 )THEN |
| INFO = 2 |
| ELSE IF( K.LT.0 )THEN |
| INFO = 3 |
| ELSE IF( LDA.LT.( K + 1 ) )THEN |
| INFO = 6 |
| ELSE IF( INCX.EQ.0 )THEN |
| INFO = 8 |
| ELSE IF( INCY.EQ.0 )THEN |
| INFO = 11 |
| END IF |
| IF( INFO.NE.0 )THEN |
| CALL XERBLA( 'DSBMV ', INFO ) |
| RETURN |
| END IF |
| * |
| * Quick return if possible. |
| * |
| IF( ( N.EQ.0 ).OR.( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) ) |
| $ RETURN |
| * |
| * Set up the start points in X and Y. |
| * |
| IF( INCX.GT.0 )THEN |
| KX = 1 |
| ELSE |
| KX = 1 - ( N - 1 )*INCX |
| END IF |
| IF( INCY.GT.0 )THEN |
| KY = 1 |
| ELSE |
| KY = 1 - ( N - 1 )*INCY |
| END IF |
| * |
| * Start the operations. In this version the elements of the array A |
| * are accessed sequentially with one pass through A. |
| * |
| * First form y := beta*y. |
| * |
| IF( BETA.NE.ONE )THEN |
| IF( INCY.EQ.1 )THEN |
| IF( BETA.EQ.ZERO )THEN |
| DO 10, I = 1, N |
| Y( I ) = ZERO |
| 10 CONTINUE |
| ELSE |
| DO 20, I = 1, N |
| Y( I ) = BETA*Y( I ) |
| 20 CONTINUE |
| END IF |
| ELSE |
| IY = KY |
| IF( BETA.EQ.ZERO )THEN |
| DO 30, I = 1, N |
| Y( IY ) = ZERO |
| IY = IY + INCY |
| 30 CONTINUE |
| ELSE |
| DO 40, I = 1, N |
| Y( IY ) = BETA*Y( IY ) |
| IY = IY + INCY |
| 40 CONTINUE |
| END IF |
| END IF |
| END IF |
| IF( ALPHA.EQ.ZERO ) |
| $ RETURN |
| IF( LSAME( UPLO, 'U' ) )THEN |
| * |
| * Form y when upper triangle of A is stored. |
| * |
| KPLUS1 = K + 1 |
| IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN |
| DO 60, J = 1, N |
| TEMP1 = ALPHA*X( J ) |
| TEMP2 = ZERO |
| L = KPLUS1 - J |
| DO 50, I = MAX( 1, J - K ), J - 1 |
| Y( I ) = Y( I ) + TEMP1*A( L + I, J ) |
| TEMP2 = TEMP2 + A( L + I, J )*X( I ) |
| 50 CONTINUE |
| Y( J ) = Y( J ) + TEMP1*A( KPLUS1, J ) + ALPHA*TEMP2 |
| 60 CONTINUE |
| ELSE |
| JX = KX |
| JY = KY |
| DO 80, J = 1, N |
| TEMP1 = ALPHA*X( JX ) |
| TEMP2 = ZERO |
| IX = KX |
| IY = KY |
| L = KPLUS1 - J |
| DO 70, I = MAX( 1, J - K ), J - 1 |
| Y( IY ) = Y( IY ) + TEMP1*A( L + I, J ) |
| TEMP2 = TEMP2 + A( L + I, J )*X( IX ) |
| IX = IX + INCX |
| IY = IY + INCY |
| 70 CONTINUE |
| Y( JY ) = Y( JY ) + TEMP1*A( KPLUS1, J ) + ALPHA*TEMP2 |
| JX = JX + INCX |
| JY = JY + INCY |
| IF( J.GT.K )THEN |
| KX = KX + INCX |
| KY = KY + INCY |
| END IF |
| 80 CONTINUE |
| END IF |
| ELSE |
| * |
| * Form y when lower triangle of A is stored. |
| * |
| IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN |
| DO 100, J = 1, N |
| TEMP1 = ALPHA*X( J ) |
| TEMP2 = ZERO |
| Y( J ) = Y( J ) + TEMP1*A( 1, J ) |
| L = 1 - J |
| DO 90, I = J + 1, MIN( N, J + K ) |
| Y( I ) = Y( I ) + TEMP1*A( L + I, J ) |
| TEMP2 = TEMP2 + A( L + I, J )*X( I ) |
| 90 CONTINUE |
| Y( J ) = Y( J ) + ALPHA*TEMP2 |
| 100 CONTINUE |
| ELSE |
| JX = KX |
| JY = KY |
| DO 120, J = 1, N |
| TEMP1 = ALPHA*X( JX ) |
| TEMP2 = ZERO |
| Y( JY ) = Y( JY ) + TEMP1*A( 1, J ) |
| L = 1 - J |
| IX = JX |
| IY = JY |
| DO 110, I = J + 1, MIN( N, J + K ) |
| IX = IX + INCX |
| IY = IY + INCY |
| Y( IY ) = Y( IY ) + TEMP1*A( L + I, J ) |
| TEMP2 = TEMP2 + A( L + I, J )*X( IX ) |
| 110 CONTINUE |
| Y( JY ) = Y( JY ) + ALPHA*TEMP2 |
| JX = JX + INCX |
| JY = JY + INCY |
| 120 CONTINUE |
| END IF |
| END IF |
| * |
| RETURN |
| * |
| * End of DSBMV . |
| * |
| END |
| subroutine dscal(n,da,dx,incx) |
| c |
| c scales a vector by a constant. |
| c uses unrolled loops for increment equal to one. |
| c jack dongarra, linpack, 3/11/78. |
| c modified 3/93 to return if incx .le. 0. |
| c modified 12/3/93, array(1) declarations changed to array(*) |
| c |
| double precision da,dx(*) |
| integer i,incx,m,mp1,n,nincx |
| c |
| if( n.le.0 .or. incx.le.0 )return |
| if(incx.eq.1)go to 20 |
| c |
| c code for increment not equal to 1 |
| c |
| nincx = n*incx |
| do 10 i = 1,nincx,incx |
| dx(i) = da*dx(i) |
| 10 continue |
| return |
| c |
| c code for increment equal to 1 |
| c |
| c |
| c clean-up loop |
| c |
| 20 m = mod(n,5) |
| if( m .eq. 0 ) go to 40 |
| do 30 i = 1,m |
| dx(i) = da*dx(i) |
| 30 continue |
| if( n .lt. 5 ) return |
| 40 mp1 = m + 1 |
| do 50 i = mp1,n,5 |
| dx(i) = da*dx(i) |
| dx(i + 1) = da*dx(i + 1) |
| dx(i + 2) = da*dx(i + 2) |
| dx(i + 3) = da*dx(i + 3) |
| dx(i + 4) = da*dx(i + 4) |
| 50 continue |
| return |
| end |
| *DECK DSDOT |
| DOUBLE PRECISION FUNCTION DSDOT (N, SX, INCX, SY, INCY) |
| C***BEGIN PROLOGUE DSDOT |
| C***PURPOSE Compute the inner product of two vectors with extended |
| C precision accumulation and result. |
| C***LIBRARY SLATEC (BLAS) |
| C***CATEGORY D1A4 |
| C***TYPE DOUBLE PRECISION (DSDOT-D, DCDOT-C) |
| C***KEYWORDS BLAS, COMPLEX VECTORS, DOT PRODUCT, INNER PRODUCT, |
| C LINEAR ALGEBRA, VECTOR |
| C***AUTHOR Lawson, C. L., (JPL) |
| C Hanson, R. J., (SNLA) |
| C Kincaid, D. R., (U. of Texas) |
| C Krogh, F. T., (JPL) |
| C***DESCRIPTION |
| C |
| C B L A S Subprogram |
| C Description of Parameters |
| C |
| C --Input-- |
| C N number of elements in input vector(s) |
| C SX single precision vector with N elements |
| C INCX storage spacing between elements of SX |
| C SY single precision vector with N elements |
| C INCY storage spacing between elements of SY |
| C |
| C --Output-- |
| C DSDOT double precision dot product (zero if N.LE.0) |
| C |
| C Returns D.P. dot product accumulated in D.P., for S.P. SX and SY |
| C DSDOT = sum for I = 0 to N-1 of SX(LX+I*INCX) * SY(LY+I*INCY), |
| C where LX = 1 if INCX .GE. 0, else LX = 1+(1-N)*INCX, and LY is |
| C defined in a similar way using INCY. |
| C |
| C***REFERENCES C. L. Lawson, R. J. Hanson, D. R. Kincaid and F. T. |
| C Krogh, Basic linear algebra subprograms for Fortran |
| C usage, Algorithm No. 539, Transactions on Mathematical |
| C Software 5, 3 (September 1979), pp. 308-323. |
| C***ROUTINES CALLED (NONE) |
| C***REVISION HISTORY (YYMMDD) |
| C 791001 DATE WRITTEN |
| C 890831 Modified array declarations. (WRB) |
| C 890831 REVISION DATE from Version 3.2 |
| C 891214 Prologue converted to Version 4.0 format. (BAB) |
| C 920310 Corrected definition of LX in DESCRIPTION. (WRB) |
| C 920501 Reformatted the REFERENCES section. (WRB) |
| C***END PROLOGUE DSDOT |
| REAL SX(*),SY(*) |
| INTEGER N, INCX, INCY |
| C |
| INTEGER KX,KY,I,NS |
| C |
| C***FIRST EXECUTABLE STATEMENT DSDOT |
| DSDOT = 0.0D0 |
| IF (N .LE. 0) RETURN |
| IF (INCX.EQ.INCY .AND. INCX.GT.0) GO TO 20 |
| C |
| C Code for unequal or nonpositive increments. |
| C |
| KX = 1 |
| KY = 1 |
| IF (INCX .LT. 0) KX = 1+(1-N)*INCX |
| IF (INCY .LT. 0) KY = 1+(1-N)*INCY |
| DO 10 I = 1,N |
| DSDOT = DSDOT + DBLE(SX(KX))*DBLE(SY(KY)) |
| KX = KX + INCX |
| KY = KY + INCY |
| 10 CONTINUE |
| RETURN |
| C |
| C Code for equal, positive, non-unit increments. |
| C |
| 20 NS = N*INCX |
| DO 30 I = 1,NS,INCX |
| DSDOT = DSDOT + DBLE(SX(I))*DBLE(SY(I)) |
| 30 CONTINUE |
| RETURN |
| END |
| SUBROUTINE DSPMV ( UPLO, N, ALPHA, AP, X, INCX, BETA, Y, INCY ) |
| * .. Scalar Arguments .. |
| DOUBLE PRECISION ALPHA, BETA |
| INTEGER INCX, INCY, N |
| CHARACTER UPLO |
| * .. Array Arguments .. |
| DOUBLE PRECISION AP( * ), X( * ), Y( * ) |
| * .. |
| * |
| * Purpose |
| * ======= |
| * |
| * DSPMV performs the matrix-vector operation |
| * |
| * y := alpha*A*x + beta*y, |
| * |
| * where alpha and beta are scalars, x and y are n element vectors and |
| * A is an n by n symmetric matrix, supplied in packed form. |
| * |
| * Parameters |
| * ========== |
| * |
| * UPLO - CHARACTER*1. |
| * On entry, UPLO specifies whether the upper or lower |
| * triangular part of the matrix A is supplied in the packed |
| * array AP as follows: |
| * |
| * UPLO = 'U' or 'u' The upper triangular part of A is |
| * supplied in AP. |
| * |
| * UPLO = 'L' or 'l' The lower triangular part of A is |
| * supplied in AP. |
| * |
| * Unchanged on exit. |
| * |
| * N - INTEGER. |
| * On entry, N specifies the order of the matrix A. |
| * N must be at least zero. |
| * Unchanged on exit. |
| * |
| * ALPHA - DOUBLE PRECISION. |
| * On entry, ALPHA specifies the scalar alpha. |
| * Unchanged on exit. |
| * |
| * AP - DOUBLE PRECISION array of DIMENSION at least |
| * ( ( n*( n + 1 ) )/2 ). |
| * Before entry with UPLO = 'U' or 'u', the array AP must |
| * contain the upper triangular part of the symmetric matrix |
| * packed sequentially, column by column, so that AP( 1 ) |
| * contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 ) |
| * and a( 2, 2 ) respectively, and so on. |
| * Before entry with UPLO = 'L' or 'l', the array AP must |
| * contain the lower triangular part of the symmetric matrix |
| * packed sequentially, column by column, so that AP( 1 ) |
| * contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 ) |
| * and a( 3, 1 ) respectively, and so on. |
| * 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 vector x. |
| * Unchanged on exit. |
| * |
| * INCX - INTEGER. |
| * On entry, INCX specifies the increment for the elements of |
| * X. INCX must not be zero. |
| * Unchanged on exit. |
| * |
| * BETA - DOUBLE PRECISION. |
| * On entry, BETA specifies the scalar beta. When BETA is |
| * supplied as zero then Y need not be set on input. |
| * Unchanged on exit. |
| * |
| * Y - DOUBLE PRECISION array of dimension at least |
| * ( 1 + ( n - 1 )*abs( INCY ) ). |
| * Before entry, the incremented array Y must contain the n |
| * element vector y. On exit, Y is overwritten by the updated |
| * vector y. |
| * |
| * INCY - INTEGER. |
| * On entry, INCY specifies the increment for the elements of |
| * Y. INCY must not be zero. |
| * Unchanged on exit. |
| * |
| * |
| * Level 2 Blas routine. |
| * |
| * -- Written on 22-October-1986. |
| * Jack Dongarra, Argonne National Lab. |
| * Jeremy Du Croz, Nag Central Office. |
| * Sven Hammarling, Nag Central Office. |
| * Richard Hanson, Sandia National Labs. |
| * |
| * |
| * .. Parameters .. |
| DOUBLE PRECISION ONE , ZERO |
| PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) |
| * .. Local Scalars .. |
| DOUBLE PRECISION TEMP1, TEMP2 |
| INTEGER I, INFO, IX, IY, J, JX, JY, K, KK, KX, KY |
| * .. External Functions .. |
| LOGICAL LSAME |
| EXTERNAL LSAME |
| * .. External Subroutines .. |
| EXTERNAL XERBLA |
| * .. |
| * .. Executable Statements .. |
| * |
| * Test the input parameters. |
| * |
| INFO = 0 |
| IF ( .NOT.LSAME( UPLO, 'U' ).AND. |
| $ .NOT.LSAME( UPLO, 'L' ) )THEN |
| INFO = 1 |
| ELSE IF( N.LT.0 )THEN |
| INFO = 2 |
| ELSE IF( INCX.EQ.0 )THEN |
| INFO = 6 |
| ELSE IF( INCY.EQ.0 )THEN |
| INFO = 9 |
| END IF |
| IF( INFO.NE.0 )THEN |
| CALL XERBLA( 'DSPMV ', INFO ) |
| RETURN |
| END IF |
| * |
| * Quick return if possible. |
| * |
| IF( ( N.EQ.0 ).OR.( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) ) |
| $ RETURN |
| * |
| * Set up the start points in X and Y. |
| * |
| IF( INCX.GT.0 )THEN |
| KX = 1 |
| ELSE |
| KX = 1 - ( N - 1 )*INCX |
| END IF |
| IF( INCY.GT.0 )THEN |
| KY = 1 |
| ELSE |
| KY = 1 - ( N - 1 )*INCY |
| END IF |
| * |
| * Start the operations. In this version the elements of the array AP |
| * are accessed sequentially with one pass through AP. |
| * |
| * First form y := beta*y. |
| * |
| IF( BETA.NE.ONE )THEN |
| IF( INCY.EQ.1 )THEN |
| IF( BETA.EQ.ZERO )THEN |
| DO 10, I = 1, N |
| Y( I ) = ZERO |
| 10 CONTINUE |
| ELSE |
| DO 20, I = 1, N |
| Y( I ) = BETA*Y( I ) |
| 20 CONTINUE |
| END IF |
| ELSE |
| IY = KY |
| IF( BETA.EQ.ZERO )THEN |
| DO 30, I = 1, N |
| Y( IY ) = ZERO |
| IY = IY + INCY |
| 30 CONTINUE |
| ELSE |
| DO 40, I = 1, N |
| Y( IY ) = BETA*Y( IY ) |
| IY = IY + INCY |
| 40 CONTINUE |
| END IF |
| END IF |
| END IF |
| IF( ALPHA.EQ.ZERO ) |
| $ RETURN |
| KK = 1 |
| IF( LSAME( UPLO, 'U' ) )THEN |
| * |
| * Form y when AP contains the upper triangle. |
| * |
| IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN |
| DO 60, J = 1, N |
| TEMP1 = ALPHA*X( J ) |
| TEMP2 = ZERO |
| K = KK |
| DO 50, I = 1, J - 1 |
| Y( I ) = Y( I ) + TEMP1*AP( K ) |
| TEMP2 = TEMP2 + AP( K )*X( I ) |
| K = K + 1 |
| 50 CONTINUE |
| Y( J ) = Y( J ) + TEMP1*AP( KK + J - 1 ) + ALPHA*TEMP2 |
| KK = KK + J |
| 60 CONTINUE |
| ELSE |
| JX = KX |
| JY = KY |
| DO 80, J = 1, N |
| TEMP1 = ALPHA*X( JX ) |
| TEMP2 = ZERO |
| IX = KX |
| IY = KY |
| DO 70, K = KK, KK + J - 2 |
| Y( IY ) = Y( IY ) + TEMP1*AP( K ) |
| TEMP2 = TEMP2 + AP( K )*X( IX ) |
| IX = IX + INCX |
| IY = IY + INCY |
| 70 CONTINUE |
| Y( JY ) = Y( JY ) + TEMP1*AP( KK + J - 1 ) + ALPHA*TEMP2 |
| JX = JX + INCX |
| JY = JY + INCY |
| KK = KK + J |
| 80 CONTINUE |
| END IF |
| ELSE |
| * |
| * Form y when AP contains the lower triangle. |
| * |
| IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN |
| DO 100, J = 1, N |
| TEMP1 = ALPHA*X( J ) |
| TEMP2 = ZERO |
| Y( J ) = Y( J ) + TEMP1*AP( KK ) |
| K = KK + 1 |
| DO 90, I = J + 1, N |
| Y( I ) = Y( I ) + TEMP1*AP( K ) |
| TEMP2 = TEMP2 + AP( K )*X( I ) |
| K = K + 1 |
| 90 CONTINUE |
| Y( J ) = Y( J ) + ALPHA*TEMP2 |
| KK = KK + ( N - J + 1 ) |
| 100 CONTINUE |
| ELSE |
| JX = KX |
| JY = KY |
| DO 120, J = 1, N |
| TEMP1 = ALPHA*X( JX ) |
| TEMP2 = ZERO |
| Y( JY ) = Y( JY ) + TEMP1*AP( KK ) |
| IX = JX |
| IY = JY |
| DO 110, K = KK + 1, KK + N - J |
| IX = IX + INCX |
| IY = IY + INCY |
| Y( IY ) = Y( IY ) + TEMP1*AP( K ) |
| TEMP2 = TEMP2 + AP( K )*X( IX ) |
| 110 CONTINUE |
| Y( JY ) = Y( JY ) + ALPHA*TEMP2 |
| JX = JX + INCX |
| JY = JY + INCY |
| KK = KK + ( N - J + 1 ) |
| 120 CONTINUE |
| END IF |
| END IF |
| * |
| RETURN |
| * |
| * End of DSPMV . |
| * |
| END |
| SUBROUTINE DSPR ( UPLO, N, ALPHA, X, INCX, AP ) |
| * .. Scalar Arguments .. |
| DOUBLE PRECISION ALPHA |
| INTEGER INCX, N |
| CHARACTER UPLO |
| * .. Array Arguments .. |
| DOUBLE PRECISION AP( * ), X( * ) |
| * .. |
| * |
| * Purpose |
| * ======= |
| * |
| * DSPR performs the symmetric rank 1 operation |
| * |
| * A := alpha*x*x' + A, |
| * |
| * where alpha is a real scalar, x is an n element vector and A is an |
| * n by n symmetric matrix, supplied in packed form. |
| * |
| * Parameters |
| * ========== |
| * |
| * UPLO - CHARACTER*1. |
| * On entry, UPLO specifies whether the upper or lower |
| * triangular part of the matrix A is supplied in the packed |
| * array AP as follows: |
| * |
| * UPLO = 'U' or 'u' The upper triangular part of A is |
| * supplied in AP. |
| * |
| * UPLO = 'L' or 'l' The lower triangular part of A is |
| * supplied in AP. |
| * |
| * Unchanged on exit. |
| * |
| * N - INTEGER. |
| * On entry, N specifies the order of the matrix A. |
| * N must be at least zero. |
| * Unchanged on exit. |
| * |
| * ALPHA - DOUBLE PRECISION. |
| * On entry, ALPHA specifies the scalar alpha. |
| * Unchanged on exit. |
| * |
| * X - DOUBLE PRECISION array of dimension at least |
| * ( 1 + ( n - 1 )*abs( INCX ) ). |
| * Before entry, the incremented array X must contain the n |
| * element vector x. |
| * Unchanged on exit. |
| * |
| * INCX - INTEGER. |
| * On entry, INCX specifies the increment for the elements of |
| * X. INCX must not be zero. |
| * Unchanged on exit. |
| * |
| * AP - DOUBLE PRECISION array of DIMENSION at least |
| * ( ( n*( n + 1 ) )/2 ). |
| * Before entry with UPLO = 'U' or 'u', the array AP must |
| * contain the upper triangular part of the symmetric matrix |
| * packed sequentially, column by column, so that AP( 1 ) |
| * contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 ) |
| * and a( 2, 2 ) respectively, and so on. On exit, the array |
| * AP is overwritten by the upper triangular part of the |
| * updated matrix. |
| * Before entry with UPLO = 'L' or 'l', the array AP must |
| * contain the lower triangular part of the symmetric matrix |
| * packed sequentially, column by column, so that AP( 1 ) |
| * contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 ) |
| * and a( 3, 1 ) respectively, and so on. On exit, the array |
| * AP is overwritten by the lower triangular part of the |
| * updated matrix. |
| * |
| * |
| * 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, K, KK, KX |
| * .. External Functions .. |
| LOGICAL LSAME |
| EXTERNAL LSAME |
| * .. External Subroutines .. |
| EXTERNAL XERBLA |
| * .. |
| * .. Executable Statements .. |
| * |
| * Test the input parameters. |
| * |
| INFO = 0 |
| IF ( .NOT.LSAME( UPLO, 'U' ).AND. |
| $ .NOT.LSAME( UPLO, 'L' ) )THEN |
| INFO = 1 |
| ELSE IF( N.LT.0 )THEN |
| INFO = 2 |
| ELSE IF( INCX.EQ.0 )THEN |
| INFO = 5 |
| END IF |
| IF( INFO.NE.0 )THEN |
| CALL XERBLA( 'DSPR ', INFO ) |
| RETURN |
| END IF |
| * |
| * Quick return if possible. |
| * |
| IF( ( N.EQ.0 ).OR.( ALPHA.EQ.ZERO ) ) |
| $ RETURN |
| * |
| * Set the start point in X if the increment is not unity. |
| * |
| 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 the array AP |
| * are accessed sequentially with one pass through AP. |
| * |
| KK = 1 |
| IF( LSAME( UPLO, 'U' ) )THEN |
| * |
| * Form A when upper triangle is stored in AP. |
| * |
| IF( INCX.EQ.1 )THEN |
| DO 20, J = 1, N |
| IF( X( J ).NE.ZERO )THEN |
| TEMP = ALPHA*X( J ) |
| K = KK |
| DO 10, I = 1, J |
| AP( K ) = AP( K ) + X( I )*TEMP |
| K = K + 1 |
| 10 CONTINUE |
| END IF |
| KK = KK + J |
| 20 CONTINUE |
| ELSE |
| JX = KX |
| DO 40, J = 1, N |
| IF( X( JX ).NE.ZERO )THEN |
| TEMP = ALPHA*X( JX ) |
| IX = KX |
| DO 30, K = KK, KK + J - 1 |
| AP( K ) = AP( K ) + X( IX )*TEMP |
| IX = IX + INCX |
| 30 CONTINUE |
| END IF |
| JX = JX + INCX |
| KK = KK + J |
| 40 CONTINUE |
| END IF |
| ELSE |
| * |
| * Form A when lower triangle is stored in AP. |
| * |
| IF( INCX.EQ.1 )THEN |
| DO 60, J = 1, N |
| IF( X( J ).NE.ZERO )THEN |
| TEMP = ALPHA*X( J ) |
| K = KK |
| DO 50, I = J, N |
| AP( K ) = AP( K ) + X( I )*TEMP |
| K = K + 1 |
| 50 CONTINUE |
| END IF |
| KK = KK + N - J + 1 |
| 60 CONTINUE |
| ELSE |
| JX = KX |
| DO 80, J = 1, N |
| IF( X( JX ).NE.ZERO )THEN |
| TEMP = ALPHA*X( JX ) |
| IX = JX |
| DO 70, K = KK, KK + N - J |
| AP( K ) = AP( K ) + X( IX )*TEMP |
| IX = IX + INCX |
| 70 CONTINUE |
| END IF |
| JX = JX + INCX |
| KK = KK + N - J + 1 |
| 80 CONTINUE |
| END IF |
| END IF |
| * |
| RETURN |
| * |
| * End of DSPR . |
| * |
| END |
| SUBROUTINE DSPR2 ( UPLO, N, ALPHA, X, INCX, Y, INCY, AP ) |
| * .. Scalar Arguments .. |
| DOUBLE PRECISION ALPHA |
| INTEGER INCX, INCY, N |
| CHARACTER UPLO |
| * .. Array Arguments .. |
| DOUBLE PRECISION AP( * ), X( * ), Y( * ) |
| * .. |
| * |
| * Purpose |
| * ======= |
| * |
| * DSPR2 performs the symmetric rank 2 operation |
| * |
| * A := alpha*x*y' + alpha*y*x' + A, |
| * |
| * where alpha is a scalar, x and y are n element vectors and A is an |
| * n by n symmetric matrix, supplied in packed form. |
| * |
| * Parameters |
| * ========== |
| * |
| * UPLO - CHARACTER*1. |
| * On entry, UPLO specifies whether the upper or lower |
| * triangular part of the matrix A is supplied in the packed |
| * array AP as follows: |
| * |
| * UPLO = 'U' or 'u' The upper triangular part of A is |
| * supplied in AP. |
| * |
| * UPLO = 'L' or 'l' The lower triangular part of A is |
| * supplied in AP. |
| * |
| * Unchanged on exit. |
| * |
| * N - INTEGER. |
| * On entry, N specifies the order of the matrix A. |
| * N must be at least zero. |
| * Unchanged on exit. |
| * |
| * ALPHA - DOUBLE PRECISION. |
| * On entry, ALPHA specifies the scalar alpha. |
| * Unchanged on exit. |
| * |
| * X - DOUBLE PRECISION array of dimension at least |
| * ( 1 + ( n - 1 )*abs( INCX ) ). |
| * Before entry, the incremented array X must contain the n |
| * element vector x. |
| * Unchanged on exit. |
| * |
| * INCX - INTEGER. |
| * On entry, INCX specifies the increment for the elements of |
| * X. INCX must not be zero. |
| * Unchanged on exit. |
| * |
| * Y - DOUBLE PRECISION array of dimension at least |
| * ( 1 + ( n - 1 )*abs( INCY ) ). |
| * Before entry, the incremented array Y must contain the n |
| * element vector y. |
| * Unchanged on exit. |
| * |
| * INCY - INTEGER. |
| * On entry, INCY specifies the increment for the elements of |
| * Y. INCY must not be zero. |
| * Unchanged on exit. |
| * |
| * AP - DOUBLE PRECISION array of DIMENSION at least |
| * ( ( n*( n + 1 ) )/2 ). |
| * Before entry with UPLO = 'U' or 'u', the array AP must |
| * contain the upper triangular part of the symmetric matrix |
| * packed sequentially, column by column, so that AP( 1 ) |
| * contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 ) |
| * and a( 2, 2 ) respectively, and so on. On exit, the array |
| * AP is overwritten by the upper triangular part of the |
| * updated matrix. |
| * Before entry with UPLO = 'L' or 'l', the array AP must |
| * contain the lower triangular part of the symmetric matrix |
| * packed sequentially, column by column, so that AP( 1 ) |
| * contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 ) |
| * and a( 3, 1 ) respectively, and so on. On exit, the array |
| * AP is overwritten by the lower triangular part of the |
| * updated matrix. |
| * |
| * |
| * 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 TEMP1, TEMP2 |
| INTEGER I, INFO, IX, IY, J, JX, JY, K, KK, KX, KY |
| * .. External Functions .. |
| LOGICAL LSAME |
| EXTERNAL LSAME |
| * .. External Subroutines .. |
| EXTERNAL XERBLA |
| * .. |
| * .. Executable Statements .. |
| * |
| * Test the input parameters. |
| * |
| INFO = 0 |
| IF ( .NOT.LSAME( UPLO, 'U' ).AND. |
| $ .NOT.LSAME( UPLO, 'L' ) )THEN |
| INFO = 1 |
| ELSE IF( N.LT.0 )THEN |
| INFO = 2 |
| ELSE IF( INCX.EQ.0 )THEN |
| INFO = 5 |
| ELSE IF( INCY.EQ.0 )THEN |
| INFO = 7 |
| END IF |
| IF( INFO.NE.0 )THEN |
| CALL XERBLA( 'DSPR2 ', INFO ) |
| RETURN |
| END IF |
| * |
| * Quick return if possible. |
| * |
| IF( ( N.EQ.0 ).OR.( ALPHA.EQ.ZERO ) ) |
| $ RETURN |
| * |
| * Set up the start points in X and Y if the increments are not both |
| * unity. |
| * |
| IF( ( INCX.NE.1 ).OR.( INCY.NE.1 ) )THEN |
| IF( INCX.GT.0 )THEN |
| KX = 1 |
| ELSE |
| KX = 1 - ( N - 1 )*INCX |
| END IF |
| IF( INCY.GT.0 )THEN |
| KY = 1 |
| ELSE |
| KY = 1 - ( N - 1 )*INCY |
| END IF |
| JX = KX |
| JY = KY |
| END IF |
| * |
| * Start the operations. In this version the elements of the array AP |
| * are accessed sequentially with one pass through AP. |
| * |
| KK = 1 |
| IF( LSAME( UPLO, 'U' ) )THEN |
| * |
| * Form A when upper triangle is stored in AP. |
| * |
| IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN |
| DO 20, J = 1, N |
| IF( ( X( J ).NE.ZERO ).OR.( Y( J ).NE.ZERO ) )THEN |
| TEMP1 = ALPHA*Y( J ) |
| TEMP2 = ALPHA*X( J ) |
| K = KK |
| DO 10, I = 1, J |
| AP( K ) = AP( K ) + X( I )*TEMP1 + Y( I )*TEMP2 |
| K = K + 1 |
| 10 CONTINUE |
| END IF |
| KK = KK + J |
| 20 CONTINUE |
| ELSE |
| DO 40, J = 1, N |
| IF( ( X( JX ).NE.ZERO ).OR.( Y( JY ).NE.ZERO ) )THEN |
| TEMP1 = ALPHA*Y( JY ) |
| TEMP2 = ALPHA*X( JX ) |
| IX = KX |
| IY = KY |
| DO 30, K = KK, KK + J - 1 |
| AP( K ) = AP( K ) + X( IX )*TEMP1 + Y( IY )*TEMP2 |
| IX = IX + INCX |
| IY = IY + INCY |
| 30 CONTINUE |
| END IF |
| JX = JX + INCX |
| JY = JY + INCY |
| KK = KK + J |
| 40 CONTINUE |
| END IF |
| ELSE |
| * |
| * Form A when lower triangle is stored in AP. |
| * |
| IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN |
| DO 60, J = 1, N |
| IF( ( X( J ).NE.ZERO ).OR.( Y( J ).NE.ZERO ) )THEN |
| TEMP1 = ALPHA*Y( J ) |
| TEMP2 = ALPHA*X( J ) |
| K = KK |
| DO 50, I = J, N |
| AP( K ) = AP( K ) + X( I )*TEMP1 + Y( I )*TEMP2 |
| K = K + 1 |
| 50 CONTINUE |
| END IF |
| KK = KK + N - J + 1 |
| 60 CONTINUE |
| ELSE |
| DO 80, J = 1, N |
| IF( ( X( JX ).NE.ZERO ).OR.( Y( JY ).NE.ZERO ) )THEN |
| TEMP1 = ALPHA*Y( JY ) |
| TEMP2 = ALPHA*X( JX ) |
| IX = JX |
| IY = JY |
| DO 70, K = KK, KK + N - J |
| AP( K ) = AP( K ) + X( IX )*TEMP1 + Y( IY )*TEMP2 |
| IX = IX + INCX |
| IY = IY + INCY |
| 70 CONTINUE |
| END IF |
| JX = JX + INCX |
| JY = JY + INCY |
| KK = KK + N - J + 1 |
| 80 CONTINUE |
| END IF |
| END IF |
| * |
| RETURN |
| * |
| * End of DSPR2 . |
| * |
| END |
| 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 modified 12/3/93, array(1) declarations changed to array(*) |
| c |
| double precision dx(*),dy(*),dtemp |
| integer i,incx,incy,ix,iy,m,mp1,n |
| c |
| 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 |
| SUBROUTINE DSYMM ( SIDE, UPLO, M, N, ALPHA, A, LDA, B, LDB, |
| $ BETA, C, LDC ) |
| * .. Scalar Arguments .. |
| CHARACTER SIDE, UPLO |
| INTEGER M, N, LDA, LDB, LDC |
| DOUBLE PRECISION ALPHA, BETA |
| * .. Array Arguments .. |
| DOUBLE PRECISION A( LDA, * ), B( LDB, * ), C( LDC, * ) |
| * .. |
| * |
| * Purpose |
| * ======= |
| * |
| * DSYMM performs one of the matrix-matrix operations |
| * |
| * C := alpha*A*B + beta*C, |
| * |
| * or |
| * |
| * C := alpha*B*A + beta*C, |
| * |
| * where alpha and beta are scalars, A is a symmetric matrix and B and |
| * C are m by n matrices. |
| * |
| * Parameters |
| * ========== |
| * |
| * SIDE - CHARACTER*1. |
| * On entry, SIDE specifies whether the symmetric matrix A |
| * appears on the left or right in the operation as follows: |
| * |
| * SIDE = 'L' or 'l' C := alpha*A*B + beta*C, |
| * |
| * SIDE = 'R' or 'r' C := alpha*B*A + beta*C, |
| * |
| * Unchanged on exit. |
| * |
| * UPLO - CHARACTER*1. |
| * On entry, UPLO specifies whether the upper or lower |
| * triangular part of the symmetric matrix A is to be |
| * referenced as follows: |
| * |
| * UPLO = 'U' or 'u' Only the upper triangular part of the |
| * symmetric matrix is to be referenced. |
| * |
| * UPLO = 'L' or 'l' Only the lower triangular part of the |
| * symmetric matrix is to be referenced. |
| * |
| * Unchanged on exit. |
| * |
| * M - INTEGER. |
| * On entry, M specifies the number of rows 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 C. |
| * N must be at least zero. |
| * Unchanged on exit. |
| * |
| * ALPHA - DOUBLE PRECISION. |
| * On entry, ALPHA specifies the scalar alpha. |
| * Unchanged on exit. |
| * |
| * A - DOUBLE PRECISION array of DIMENSION ( LDA, ka ), where ka is |
| * m when SIDE = 'L' or 'l' and is n otherwise. |
| * Before entry with SIDE = 'L' or 'l', the m by m part of |
| * the array A must contain the symmetric matrix, such that |
| * when UPLO = 'U' or 'u', the leading m by m upper triangular |
| * part of the array A must contain the upper triangular part |
| * of the symmetric matrix and the strictly lower triangular |
| * part of A is not referenced, and when UPLO = 'L' or 'l', |
| * the leading m by m lower triangular part of the array A |
| * must contain the lower triangular part of the symmetric |
| * matrix and the strictly upper triangular part of A is not |
| * referenced. |
| * Before entry with SIDE = 'R' or 'r', the n by n part of |
| * the array A must contain the symmetric matrix, such that |
| * when UPLO = 'U' or 'u', the leading n by n upper triangular |
| * part of the array A must contain the upper triangular part |
| * of the symmetric matrix and the strictly lower triangular |
| * part of A is not referenced, and when UPLO = 'L' or 'l', |
| * the leading n by n lower triangular part of the array A |
| * must contain the lower triangular part of the symmetric |
| * matrix and the strictly upper triangular part of A is not |
| * referenced. |
| * 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 ), otherwise 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 matrix B. |
| * Unchanged on exit. |
| * |
| * 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. |
| * |
| * 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 updated |
| * matrix. |
| * |
| * 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 UPPER |
| INTEGER I, INFO, J, K, NROWA |
| DOUBLE PRECISION TEMP1, TEMP2 |
| * .. Parameters .. |
| DOUBLE PRECISION ONE , ZERO |
| PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) |
| * .. |
| * .. Executable Statements .. |
| * |
| * Set NROWA as the number of rows of A. |
| * |
| IF( LSAME( SIDE, 'L' ) )THEN |
| NROWA = M |
| ELSE |
| NROWA = N |
| END IF |
| UPPER = LSAME( UPLO, 'U' ) |
| * |
| * Test the input parameters. |
| * |
| INFO = 0 |
| IF( ( .NOT.LSAME( SIDE, 'L' ) ).AND. |
| $ ( .NOT.LSAME( SIDE, 'R' ) ) )THEN |
| INFO = 1 |
| ELSE IF( ( .NOT.UPPER ).AND. |
| $ ( .NOT.LSAME( UPLO, 'L' ) ) )THEN |
| INFO = 2 |
| ELSE IF( M .LT.0 )THEN |
| INFO = 3 |
| ELSE IF( N .LT.0 )THEN |
| INFO = 4 |
| ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN |
| INFO = 7 |
| ELSE IF( LDB.LT.MAX( 1, M ) )THEN |
| INFO = 9 |
| ELSE IF( LDC.LT.MAX( 1, M ) )THEN |
| INFO = 12 |
| END IF |
| IF( INFO.NE.0 )THEN |
| CALL XERBLA( 'DSYMM ', INFO ) |
| RETURN |
| END IF |
| * |
| * Quick return if possible. |
| * |
| IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR. |
| $ ( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) ) |
| $ RETURN |
| * |
| * And when 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( LSAME( SIDE, 'L' ) )THEN |
| * |
| * Form C := alpha*A*B + beta*C. |
| * |
| IF( UPPER )THEN |
| DO 70, J = 1, N |
| DO 60, I = 1, M |
| TEMP1 = ALPHA*B( I, J ) |
| TEMP2 = ZERO |
| DO 50, K = 1, I - 1 |
| C( K, J ) = C( K, J ) + TEMP1 *A( K, I ) |
| TEMP2 = TEMP2 + B( K, J )*A( K, I ) |
| 50 CONTINUE |
| IF( BETA.EQ.ZERO )THEN |
| C( I, J ) = TEMP1*A( I, I ) + ALPHA*TEMP2 |
| ELSE |
| C( I, J ) = BETA *C( I, J ) + |
| $ TEMP1*A( I, I ) + ALPHA*TEMP2 |
| END IF |
| 60 CONTINUE |
| 70 CONTINUE |
| ELSE |
| DO 100, J = 1, N |
| DO 90, I = M, 1, -1 |
| TEMP1 = ALPHA*B( I, J ) |
| TEMP2 = ZERO |
| DO 80, K = I + 1, M |
| C( K, J ) = C( K, J ) + TEMP1 *A( K, I ) |
| TEMP2 = TEMP2 + B( K, J )*A( K, I ) |
| 80 CONTINUE |
| IF( BETA.EQ.ZERO )THEN |
| C( I, J ) = TEMP1*A( I, I ) + ALPHA*TEMP2 |
| ELSE |
| C( I, J ) = BETA *C( I, J ) + |
| $ TEMP1*A( I, I ) + ALPHA*TEMP2 |
| END IF |
| 90 CONTINUE |
| 100 CONTINUE |
| END IF |
| ELSE |
| * |
| * Form C := alpha*B*A + beta*C. |
| * |
| DO 170, J = 1, N |
| TEMP1 = ALPHA*A( J, J ) |
| IF( BETA.EQ.ZERO )THEN |
| DO 110, I = 1, M |
| C( I, J ) = TEMP1*B( I, J ) |
| 110 CONTINUE |
| ELSE |
| DO 120, I = 1, M |
| C( I, J ) = BETA*C( I, J ) + TEMP1*B( I, J ) |
| 120 CONTINUE |
| END IF |
| DO 140, K = 1, J - 1 |
| IF( UPPER )THEN |
| TEMP1 = ALPHA*A( K, J ) |
| ELSE |
| TEMP1 = ALPHA*A( J, K ) |
| END IF |
| DO 130, I = 1, M |
| C( I, J ) = C( I, J ) + TEMP1*B( I, K ) |
| 130 CONTINUE |
| 140 CONTINUE |
| DO 160, K = J + 1, N |
| IF( UPPER )THEN |
| TEMP1 = ALPHA*A( J, K ) |
| ELSE |
| TEMP1 = ALPHA*A( K, J ) |
| END IF |
| DO 150, I = 1, M |
| C( I, J ) = C( I, J ) + TEMP1*B( I, K ) |
| 150 CONTINUE |
| 160 CONTINUE |
| 170 CONTINUE |
| END IF |
| * |
| RETURN |
| * |
| * End of DSYMM . |
| * |
| END |
| SUBROUTINE DSYMV ( UPLO, N, ALPHA, A, LDA, X, INCX, |
| $ BETA, Y, INCY ) |
| * .. Scalar Arguments .. |
| DOUBLE PRECISION ALPHA, BETA |
| INTEGER INCX, INCY, LDA, N |
| CHARACTER UPLO |
| * .. Array Arguments .. |
| DOUBLE PRECISION A( LDA, * ), X( * ), Y( * ) |
| * .. |
| * |
| * Purpose |
| * ======= |
| * |
| * DSYMV performs the matrix-vector operation |
| * |
| * y := alpha*A*x + beta*y, |
| * |
| * where alpha and beta are scalars, x and y are n element vectors and |
| * A is an n by n symmetric matrix. |
| * |
| * Parameters |
| * ========== |
| * |
| * UPLO - CHARACTER*1. |
| * On entry, UPLO specifies whether the upper or lower |
| * triangular part of the array A is to be referenced as |
| * follows: |
| * |
| * UPLO = 'U' or 'u' Only the upper triangular part of A |
| * is to be referenced. |
| * |
| * UPLO = 'L' or 'l' Only the lower triangular part of A |
| * is to be referenced. |
| * |
| * Unchanged on exit. |
| * |
| * N - INTEGER. |
| * On entry, N specifies the order of the matrix A. |
| * N must be at least zero. |
| * Unchanged on exit. |
| * |
| * ALPHA - DOUBLE PRECISION. |
| * On entry, ALPHA specifies the scalar alpha. |
| * Unchanged on exit. |
| * |
| * A - DOUBLE PRECISION array of DIMENSION ( LDA, n ). |
| * Before entry with UPLO = 'U' or 'u', the leading n by n |
| * upper triangular part of the array A must contain the upper |
| * triangular part of the symmetric 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 part of the symmetric matrix and the strictly |
| * upper triangular part of A is not referenced. |
| * 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 vector x. |
| * Unchanged on exit. |
| * |
| * INCX - INTEGER. |
| * On entry, INCX specifies the increment for the elements of |
| * X. INCX must not be zero. |
| * Unchanged on exit. |
| * |
| * BETA - DOUBLE PRECISION. |
| * On entry, BETA specifies the scalar beta. When BETA is |
| * supplied as zero then Y need not be set on input. |
| * Unchanged on exit. |
| * |
| * Y - DOUBLE PRECISION array of dimension at least |
| * ( 1 + ( n - 1 )*abs( INCY ) ). |
| * Before entry, the incremented array Y must contain the n |
| * element vector y. On exit, Y is overwritten by the updated |
| * vector y. |
| * |
| * INCY - INTEGER. |
| * On entry, INCY specifies the increment for the elements of |
| * Y. INCY must not be zero. |
| * Unchanged on exit. |
| * |
| * |
| * Level 2 Blas routine. |
| * |
| * -- Written on 22-October-1986. |
| * Jack Dongarra, Argonne National Lab. |
| * Jeremy Du Croz, Nag Central Office. |
| * Sven Hammarling, Nag Central Office. |
| * Richard Hanson, Sandia National Labs. |
| * |
| * |
| * .. Parameters .. |
| DOUBLE PRECISION ONE , ZERO |
| PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) |
| * .. Local Scalars .. |
| DOUBLE PRECISION TEMP1, TEMP2 |
| INTEGER I, INFO, IX, IY, J, JX, JY, KX, KY |
| * .. 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( N.LT.0 )THEN |
| INFO = 2 |
| ELSE IF( LDA.LT.MAX( 1, N ) )THEN |
| INFO = 5 |
| ELSE IF( INCX.EQ.0 )THEN |
| INFO = 7 |
| ELSE IF( INCY.EQ.0 )THEN |
| INFO = 10 |
| END IF |
| IF( INFO.NE.0 )THEN |
| CALL XERBLA( 'DSYMV ', INFO ) |
| RETURN |
| END IF |
| * |
| * Quick return if possible. |
| * |
| IF( ( N.EQ.0 ).OR.( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) ) |
| $ RETURN |
| * |
| * Set up the start points in X and Y. |
| * |
| IF( INCX.GT.0 )THEN |
| KX = 1 |
| ELSE |
| KX = 1 - ( N - 1 )*INCX |
| END IF |
| IF( INCY.GT.0 )THEN |
| KY = 1 |
| ELSE |
| KY = 1 - ( N - 1 )*INCY |
| END IF |
| * |
| * Start the operations. In this version the elements of A are |
| * accessed sequentially with one pass through the triangular part |
| * of A. |
| * |
| * First form y := beta*y. |
| * |
| IF( BETA.NE.ONE )THEN |
| IF( INCY.EQ.1 )THEN |
| IF( BETA.EQ.ZERO )THEN |
| DO 10, I = 1, N |
| Y( I ) = ZERO |
| 10 CONTINUE |
| ELSE |
| DO 20, I = 1, N |
| Y( I ) = BETA*Y( I ) |
| 20 CONTINUE |
| END IF |
| ELSE |
| IY = KY |
| IF( BETA.EQ.ZERO )THEN |
| DO 30, I = 1, N |
| Y( IY ) = ZERO |
| IY = IY + INCY |
| 30 CONTINUE |
| ELSE |
| DO 40, I = 1, N |
| Y( IY ) = BETA*Y( IY ) |
| IY = IY + INCY |
| 40 CONTINUE |
| END IF |
| END IF |
| END IF |
| IF( ALPHA.EQ.ZERO ) |
| $ RETURN |
| IF( LSAME( UPLO, 'U' ) )THEN |
| * |
| * Form y when A is stored in upper triangle. |
| * |
| IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN |
| DO 60, J = 1, N |
| TEMP1 = ALPHA*X( J ) |
| TEMP2 = ZERO |
| DO 50, I = 1, J - 1 |
| Y( I ) = Y( I ) + TEMP1*A( I, J ) |
| TEMP2 = TEMP2 + A( I, J )*X( I ) |
| 50 CONTINUE |
| Y( J ) = Y( J ) + TEMP1*A( J, J ) + ALPHA*TEMP2 |
| 60 CONTINUE |
| ELSE |
| JX = KX |
| JY = KY |
| DO 80, J = 1, N |
| TEMP1 = ALPHA*X( JX ) |
| TEMP2 = ZERO |
| IX = KX |
| IY = KY |
| DO 70, I = 1, J - 1 |
| Y( IY ) = Y( IY ) + TEMP1*A( I, J ) |
| TEMP2 = TEMP2 + A( I, J )*X( IX ) |
| IX = IX + INCX |
| IY = IY + INCY |
| 70 CONTINUE |
| Y( JY ) = Y( JY ) + TEMP1*A( J, J ) + ALPHA*TEMP2 |
| JX = JX + INCX |
| JY = JY + INCY |
| 80 CONTINUE |
| END IF |
| ELSE |
| * |
| * Form y when A is stored in lower triangle. |
| * |
| IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN |
| DO 100, J = 1, N |
| TEMP1 = ALPHA*X( J ) |
| TEMP2 = ZERO |
| Y( J ) = Y( J ) + TEMP1*A( J, J ) |
| DO 90, I = J + 1, N |
| Y( I ) = Y( I ) + TEMP1*A( I, J ) |
| TEMP2 = TEMP2 + A( I, J )*X( I ) |
| 90 CONTINUE |
| Y( J ) = Y( J ) + ALPHA*TEMP2 |
| 100 CONTINUE |
| ELSE |
| JX = KX |
| JY = KY |
| DO 120, J = 1, N |
| TEMP1 = ALPHA*X( JX ) |
| TEMP2 = ZERO |
| Y( JY ) = Y( JY ) + TEMP1*A( J, J ) |
| IX = JX |
| IY = JY |
| DO 110, I = J + 1, N |
| IX = IX + INCX |
| IY = IY + INCY |
| Y( IY ) = Y( IY ) + TEMP1*A( I, J ) |
| TEMP2 = TEMP2 + A( I, J )*X( IX ) |
| 110 CONTINUE |
| Y( JY ) = Y( JY ) + ALPHA*TEMP2 |
| JX = JX + INCX |
| JY = JY + INCY |
| 120 CONTINUE |
| END IF |
| END IF |
| * |
| RETURN |
| * |
| * End of DSYMV . |
| * |
| END |
| SUBROUTINE DSYR ( UPLO, N, ALPHA, X, INCX, A, LDA ) |
| * .. Scalar Arguments .. |
| DOUBLE PRECISION ALPHA |
| INTEGER INCX, LDA, N |
| CHARACTER UPLO |
| * .. Array Arguments .. |
| DOUBLE PRECISION A( LDA, * ), X( * ) |
| * .. |
| * |
| * Purpose |
| * ======= |
| * |
| * DSYR performs the symmetric rank 1 operation |
| * |
| * A := alpha*x*x' + A, |
| * |
| * where alpha is a real scalar, x is an n element vector and A is an |
| * n by n symmetric matrix. |
| * |
| * Parameters |
| * ========== |
| * |
| * UPLO - CHARACTER*1. |
| * On entry, UPLO specifies whether the upper or lower |
| * triangular part of the array A is to be referenced as |
| * follows: |
| * |
| * UPLO = 'U' or 'u' Only the upper triangular part of A |
| * is to be referenced. |
| * |
| * UPLO = 'L' or 'l' Only the lower triangular part of A |
| * is to be referenced. |
| * |
| * Unchanged on exit. |
| * |
| * N - INTEGER. |
| * On entry, N specifies the order of the matrix A. |
| * N must be at least zero. |
| * Unchanged on exit. |
| * |
| * ALPHA - DOUBLE PRECISION. |
| * On entry, ALPHA specifies the scalar alpha. |
| * Unchanged on exit. |
| * |
| * X - DOUBLE PRECISION array of dimension at least |
| * ( 1 + ( n - 1 )*abs( INCX ) ). |
| * Before entry, the incremented array X must contain the n |
| * element vector x. |
| * Unchanged on exit. |
| * |
| * INCX - INTEGER. |
| * On entry, INCX specifies the increment for the elements of |
| * X. INCX must not be zero. |
| * Unchanged on exit. |
| * |
| * 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 part of the symmetric matrix and the strictly |
| * lower triangular part of A is not referenced. On exit, the |
| * upper triangular part of the array A is overwritten by the |
| * upper triangular part of the updated matrix. |
| * Before entry with UPLO = 'L' or 'l', the leading n by n |
| * lower triangular part of the array A must contain the lower |
| * triangular part of the symmetric matrix and the strictly |
| * upper triangular part of A is not referenced. On exit, the |
| * lower triangular part of the array A is overwritten by the |
| * lower triangular part of the updated matrix. |
| * |
| * LDA - INTEGER. |
| * On entry, LDA specifies the first dimension of A as declared |
| * in the calling (sub) program. LDA must be at least |
| * max( 1, n ). |
| * 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 |
| * .. 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( N.LT.0 )THEN |
| INFO = 2 |
| ELSE IF( INCX.EQ.0 )THEN |
| INFO = 5 |
| ELSE IF( LDA.LT.MAX( 1, N ) )THEN |
| INFO = 7 |
| END IF |
| IF( INFO.NE.0 )THEN |
| CALL XERBLA( 'DSYR ', INFO ) |
| RETURN |
| END IF |
| * |
| * Quick return if possible. |
| * |
| IF( ( N.EQ.0 ).OR.( ALPHA.EQ.ZERO ) ) |
| $ RETURN |
| * |
| * Set the start point in X if the increment is not unity. |
| * |
| 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 the triangular part |
| * of A. |
| * |
| IF( LSAME( UPLO, 'U' ) )THEN |
| * |
| * Form A when A is stored in upper triangle. |
| * |
| IF( INCX.EQ.1 )THEN |
| DO 20, J = 1, N |
| IF( X( J ).NE.ZERO )THEN |
| TEMP = ALPHA*X( J ) |
| DO 10, I = 1, J |
| A( I, J ) = A( I, J ) + X( I )*TEMP |
| 10 CONTINUE |
| END IF |
| 20 CONTINUE |
| ELSE |
| JX = KX |
| DO 40, J = 1, N |
| IF( X( JX ).NE.ZERO )THEN |
| TEMP = ALPHA*X( JX ) |
| IX = KX |
| DO 30, I = 1, J |
| A( I, J ) = A( I, J ) + X( IX )*TEMP |
| IX = IX + INCX |
| 30 CONTINUE |
| END IF |
| JX = JX + INCX |
| 40 CONTINUE |
| END IF |
| ELSE |
| * |
| * Form A when A is stored in lower triangle. |
| * |
| IF( INCX.EQ.1 )THEN |
| DO 60, J = 1, N |
| IF( X( J ).NE.ZERO )THEN |
| TEMP = ALPHA*X( J ) |
| DO 50, I = J, N |
| A( I, J ) = A( I, J ) + X( I )*TEMP |
| 50 CONTINUE |
| END IF |
| 60 CONTINUE |
| ELSE |
| JX = KX |
| DO 80, J = 1, N |
| IF( X( JX ).NE.ZERO )THEN |
| TEMP = ALPHA*X( JX ) |
| IX = JX |
| DO 70, I = J, N |
| A( I, J ) = A( I, J ) + X( IX )*TEMP |
| IX = IX + INCX |
| 70 CONTINUE |
| END IF |
| JX = JX + INCX |
| 80 CONTINUE |
| END IF |
| END IF |
| * |
| RETURN |
| * |
| * End of DSYR . |
| * |
| END |
| SUBROUTINE DSYR2 ( UPLO, N, ALPHA, X, INCX, Y, INCY, A, LDA ) |
| * .. Scalar Arguments .. |
| DOUBLE PRECISION ALPHA |
| INTEGER INCX, INCY, LDA, N |
| CHARACTER UPLO |
| * .. Array Arguments .. |
| DOUBLE PRECISION A( LDA, * ), X( * ), Y( * ) |
| * .. |
| * |
| * Purpose |
| * ======= |
| * |
| * DSYR2 performs the symmetric rank 2 operation |
| * |
| * A := alpha*x*y' + alpha*y*x' + A, |
| * |
| * where alpha is a scalar, x and y are n element vectors and A is an n |
| * by n symmetric matrix. |
| * |
| * Parameters |
| * ========== |
| * |
| * UPLO - CHARACTER*1. |
| * On entry, UPLO specifies whether the upper or lower |
| * triangular part of the array A is to be referenced as |
| * follows: |
| * |
| * UPLO = 'U' or 'u' Only the upper triangular part of A |
| * is to be referenced. |
| * |
| * UPLO = 'L' or 'l' Only the lower triangular part of A |
| * is to be referenced. |
| * |
| * Unchanged on exit. |
| * |
| * N - INTEGER. |
| * On entry, N specifies the order of the matrix A. |
| * N must be at least zero. |
| * Unchanged on exit. |
| * |
| * ALPHA - DOUBLE PRECISION. |
| * On entry, ALPHA specifies the scalar alpha. |
| * Unchanged on exit. |
| * |
| * X - DOUBLE PRECISION array of dimension at least |
| * ( 1 + ( n - 1 )*abs( INCX ) ). |
| * Before entry, the incremented array X must contain the n |
| * element vector x. |
| * Unchanged on exit. |
| * |
| * INCX - INTEGER. |
| * On entry, INCX specifies the increment for the elements of |
| * X. INCX must not be zero. |
| * Unchanged on exit. |
| * |
| * Y - DOUBLE PRECISION array of dimension at least |
| * ( 1 + ( n - 1 )*abs( INCY ) ). |
| * Before entry, the incremented array Y must contain the n |
| * element vector y. |
| * Unchanged on exit. |
| * |
| * INCY - INTEGER. |
| * On entry, INCY specifies the increment for the elements of |
| * Y. INCY must not be zero. |
| * Unchanged on exit. |
| * |
| * A - DOUBLE PRECISION array of DIMENSION ( LDA, n ). |
| * Before entry with UPLO = 'U' or 'u', the leading n by n |
| * upper triangular part of the array A must contain the upper |
| * triangular part of the symmetric matrix and the strictly |
| * lower triangular part of A is not referenced. On exit, the |
| * upper triangular part of the array A is overwritten by the |
| * upper triangular part of the updated matrix. |
| * Before entry with UPLO = 'L' or 'l', the leading n by n |
| * lower triangular part of the array A must contain the lower |
| * triangular part of the symmetric matrix and the strictly |
| * upper triangular part of A is not referenced. On exit, the |
| * lower triangular part of the array A is overwritten by the |
| * lower triangular part of the updated matrix. |
| * |
| * LDA - INTEGER. |
| * On entry, LDA specifies the first dimension of A as declared |
| * in the calling (sub) program. LDA must be at least |
| * max( 1, n ). |
| * 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 TEMP1, TEMP2 |
| INTEGER I, INFO, IX, IY, J, JX, JY, KX, KY |
| * .. 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( N.LT.0 )THEN |
| INFO = 2 |
| ELSE IF( INCX.EQ.0 )THEN |
| INFO = 5 |
| ELSE IF( INCY.EQ.0 )THEN |
| INFO = 7 |
| ELSE IF( LDA.LT.MAX( 1, N ) )THEN |
| INFO = 9 |
| END IF |
| IF( INFO.NE.0 )THEN |
| CALL XERBLA( 'DSYR2 ', INFO ) |
| RETURN |
| END IF |
| * |
| * Quick return if possible. |
| * |
| IF( ( N.EQ.0 ).OR.( ALPHA.EQ.ZERO ) ) |
| $ RETURN |
| * |
| * Set up the start points in X and Y if the increments are not both |
| * unity. |
| * |
| IF( ( INCX.NE.1 ).OR.( INCY.NE.1 ) )THEN |
| IF( INCX.GT.0 )THEN |
| KX = 1 |
| ELSE |
| KX = 1 - ( N - 1 )*INCX |
| END IF |
| IF( INCY.GT.0 )THEN |
| KY = 1 |
| ELSE |
| KY = 1 - ( N - 1 )*INCY |
| END IF |
| JX = KX |
| JY = KY |
| END IF |
| * |
| * Start the operations. In this version the elements of A are |
| * accessed sequentially with one pass through the triangular part |
| * of A. |
| * |
| IF( LSAME( UPLO, 'U' ) )THEN |
| * |
| * Form A when A is stored in the upper triangle. |
| * |
| IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN |
| DO 20, J = 1, N |
| IF( ( X( J ).NE.ZERO ).OR.( Y( J ).NE.ZERO ) )THEN |
| TEMP1 = ALPHA*Y( J ) |
| TEMP2 = ALPHA*X( J ) |
| DO 10, I = 1, J |
| A( I, J ) = A( I, J ) + X( I )*TEMP1 + Y( I )*TEMP2 |
| 10 CONTINUE |
| END IF |
| 20 CONTINUE |
| ELSE |
| DO 40, J = 1, N |
| IF( ( X( JX ).NE.ZERO ).OR.( Y( JY ).NE.ZERO ) )THEN |
| TEMP1 = ALPHA*Y( JY ) |
| TEMP2 = ALPHA*X( JX ) |
| IX = KX |
| IY = KY |
| DO 30, I = 1, J |
| A( I, J ) = A( I, J ) + X( IX )*TEMP1 |
| $ + Y( IY )*TEMP2 |
| IX = IX + INCX |
| IY = IY + INCY |
| 30 CONTINUE |
| END IF |
| JX = JX + INCX |
| JY = JY + INCY |
| 40 CONTINUE |
| END IF |
| ELSE |
| * |
| * Form A when A is stored in the lower triangle. |
| * |
| IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN |
| DO 60, J = 1, N |
| IF( ( X( J ).NE.ZERO ).OR.( Y( J ).NE.ZERO ) )THEN |
| TEMP1 = ALPHA*Y( J ) |
| TEMP2 = ALPHA*X( J ) |
| DO 50, I = J, N |
| A( I, J ) = A( I, J ) + X( I )*TEMP1 + Y( I )*TEMP2 |
| 50 CONTINUE |
| END IF |
| 60 CONTINUE |
| ELSE |
| DO 80, J = 1, N |
| IF( ( X( JX ).NE.ZERO ).OR.( Y( JY ).NE.ZERO ) )THEN |
| TEMP1 = ALPHA*Y( JY ) |
| TEMP2 = ALPHA*X( JX ) |
| IX = JX |
| IY = JY |
| DO 70, I = J, N |
| A( I, J ) = A( I, J ) + X( IX )*TEMP1 |
| $ + Y( IY )*TEMP2 |
| IX = IX + INCX |
| IY = IY + INCY |
| 70 CONTINUE |
| END IF |
| JX = JX + INCX |
| JY = JY + INCY |
| 80 CONTINUE |
| END IF |
| END IF |
| * |
| RETURN |
| * |
| * End of DSYR2 . |
| * |
| END |
| SUBROUTINE DSYR2K( UPLO, TRANS, N, K, ALPHA, A, LDA, B, LDB, |
| $ BETA, C, LDC ) |
| * .. Scalar Arguments .. |
| CHARACTER UPLO, TRANS |
| INTEGER N, K, LDA, LDB, LDC |
| DOUBLE PRECISION ALPHA, BETA |
| * .. Array Arguments .. |
| DOUBLE PRECISION A( LDA, * ), B( LDB, * ), C( LDC, * ) |
| * .. |
| * |
| * Purpose |
| * ======= |
| * |
| * DSYR2K performs one of the symmetric rank 2k operations |
| * |
| * C := alpha*A*B' + alpha*B*A' + beta*C, |
| * |
| * or |
| * |
| * C := alpha*A'*B + alpha*B'*A + beta*C, |
| * |
| * where alpha and beta are scalars, C is an n by n symmetric matrix |
| * and A and B are n by k matrices in the first case and k by n |
| * matrices in the second case. |
| * |
| * Parameters |
| * ========== |
| * |
| * UPLO - CHARACTER*1. |
| * On entry, UPLO specifies whether the upper or lower |
| * triangular part of the array C is to be referenced as |
| * follows: |
| * |
| * UPLO = 'U' or 'u' Only the upper triangular part of C |
| * is to be referenced. |
| * |
| * UPLO = 'L' or 'l' Only the lower triangular part of C |
| * is to be referenced. |
| * |
| * Unchanged on exit. |
| * |
| * TRANS - CHARACTER*1. |
| * On entry, TRANS specifies the operation to be performed as |
| * follows: |
| * |
| * TRANS = 'N' or 'n' C := alpha*A*B' + alpha*B*A' + |
| * beta*C. |
| * |
| * TRANS = 'T' or 't' C := alpha*A'*B + alpha*B'*A + |
| * beta*C. |
| * |
| * TRANS = 'C' or 'c' C := alpha*A'*B + alpha*B'*A + |
| * beta*C. |
| * |
| * Unchanged on exit. |
| * |
| * N - INTEGER. |
| * On entry, N specifies the order of the matrix C. N must be |
| * at least zero. |
| * Unchanged on exit. |
| * |
| * K - INTEGER. |
| * On entry with TRANS = 'N' or 'n', K specifies the number |
| * of columns of the matrices A and B, and on entry with |
| * TRANS = 'T' or 't' or 'C' or 'c', K specifies the number |
| * of rows of the matrices A and 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 TRANS = 'N' or 'n', and is n otherwise. |
| * Before entry with TRANS = 'N' or 'n', the leading n by k |
| * part of the array A must contain the matrix A, otherwise |
| * the leading k by n 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 TRANS = 'N' or 'n' |
| * then LDA must be at least max( 1, n ), otherwise LDA must |
| * be at least max( 1, k ). |
| * Unchanged on exit. |
| * |
| * B - DOUBLE PRECISION array of DIMENSION ( LDB, kb ), where kb is |
| * k when TRANS = 'N' or 'n', and is n otherwise. |
| * Before entry with TRANS = 'N' or 'n', the leading n by k |
| * part of the array B must contain the matrix B, otherwise |
| * the leading k by n 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 TRANS = 'N' or 'n' |
| * then LDB must be at least max( 1, n ), otherwise LDB must |
| * be at least max( 1, k ). |
| * Unchanged on exit. |
| * |
| * BETA - DOUBLE PRECISION. |
| * On entry, BETA specifies the scalar beta. |
| * Unchanged on exit. |
| * |
| * C - DOUBLE PRECISION array of DIMENSION ( LDC, n ). |
| * Before entry with UPLO = 'U' or 'u', the leading n by n |
| * upper triangular part of the array C must contain the upper |
| * triangular part of the symmetric matrix and the strictly |
| * lower triangular part of C is not referenced. On exit, the |
| * upper triangular part of the array C is overwritten by the |
| * upper triangular part of the updated matrix. |
| * Before entry with UPLO = 'L' or 'l', the leading n by n |
| * lower triangular part of the array C must contain the lower |
| * triangular part of the symmetric matrix and the strictly |
| * upper triangular part of C is not referenced. On exit, the |
| * lower triangular part of the array C is overwritten by the |
| * lower triangular part of the updated matrix. |
| * |
| * 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, n ). |
| * 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 UPPER |
| INTEGER I, INFO, J, L, NROWA |
| DOUBLE PRECISION TEMP1, TEMP2 |
| * .. Parameters .. |
| DOUBLE PRECISION ONE , ZERO |
| PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) |
| * .. |
| * .. Executable Statements .. |
| * |
| * Test the input parameters. |
| * |
| IF( LSAME( TRANS, 'N' ) )THEN |
| NROWA = N |
| ELSE |
| NROWA = K |
| END IF |
| UPPER = LSAME( UPLO, 'U' ) |
| * |
| INFO = 0 |
| IF( ( .NOT.UPPER ).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( N .LT.0 )THEN |
| INFO = 3 |
| ELSE IF( K .LT.0 )THEN |
| INFO = 4 |
| ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN |
| INFO = 7 |
| ELSE IF( LDB.LT.MAX( 1, NROWA ) )THEN |
| INFO = 9 |
| ELSE IF( LDC.LT.MAX( 1, N ) )THEN |
| INFO = 12 |
| END IF |
| IF( INFO.NE.0 )THEN |
| CALL XERBLA( 'DSYR2K', INFO ) |
| RETURN |
| END IF |
| * |
| * Quick return if possible. |
| * |
| IF( ( N.EQ.0 ).OR. |
| $ ( ( ( ALPHA.EQ.ZERO ).OR.( K.EQ.0 ) ).AND.( BETA.EQ.ONE ) ) ) |
| $ RETURN |
| * |
| * And when alpha.eq.zero. |
| * |
| IF( ALPHA.EQ.ZERO )THEN |
| IF( UPPER )THEN |
| IF( BETA.EQ.ZERO )THEN |
| DO 20, J = 1, N |
| DO 10, I = 1, J |
| C( I, J ) = ZERO |
| 10 CONTINUE |
| 20 CONTINUE |
| ELSE |
| DO 40, J = 1, N |
| DO 30, I = 1, J |
| C( I, J ) = BETA*C( I, J ) |
| 30 CONTINUE |
| 40 CONTINUE |
| END IF |
| ELSE |
| IF( BETA.EQ.ZERO )THEN |
| DO 60, J = 1, N |
| DO 50, I = J, N |
| C( I, J ) = ZERO |
| 50 CONTINUE |
| 60 CONTINUE |
| ELSE |
| DO 80, J = 1, N |
| DO 70, I = J, N |
| C( I, J ) = BETA*C( I, J ) |
| 70 CONTINUE |
| 80 CONTINUE |
| END IF |
| END IF |
| RETURN |
| END IF |
| * |
| * Start the operations. |
| * |
| IF( LSAME( TRANS, 'N' ) )THEN |
| * |
| * Form C := alpha*A*B' + alpha*B*A' + C. |
| * |
| IF( UPPER )THEN |
| DO 130, J = 1, N |
| IF( BETA.EQ.ZERO )THEN |
| DO 90, I = 1, J |
| C( I, J ) = ZERO |
| 90 CONTINUE |
| ELSE IF( BETA.NE.ONE )THEN |
| DO 100, I = 1, J |
| C( I, J ) = BETA*C( I, J ) |
| 100 CONTINUE |
| END IF |
| DO 120, L = 1, K |
| IF( ( A( J, L ).NE.ZERO ).OR. |
| $ ( B( J, L ).NE.ZERO ) )THEN |
| TEMP1 = ALPHA*B( J, L ) |
| TEMP2 = ALPHA*A( J, L ) |
| DO 110, I = 1, J |
| C( I, J ) = C( I, J ) + |
| $ A( I, L )*TEMP1 + B( I, L )*TEMP2 |
| 110 CONTINUE |
| END IF |
| 120 CONTINUE |
| 130 CONTINUE |
| ELSE |
| DO 180, J = 1, N |
| IF( BETA.EQ.ZERO )THEN |
| DO 140, I = J, N |
| C( I, J ) = ZERO |
| 140 CONTINUE |
| ELSE IF( BETA.NE.ONE )THEN |
| DO 150, I = J, N |
| C( I, J ) = BETA*C( I, J ) |
| 150 CONTINUE |
| END IF |
| DO 170, L = 1, K |
| IF( ( A( J, L ).NE.ZERO ).OR. |
| $ ( B( J, L ).NE.ZERO ) )THEN |
| TEMP1 = ALPHA*B( J, L ) |
| TEMP2 = ALPHA*A( J, L ) |
| DO 160, I = J, N |
| C( I, J ) = C( I, J ) + |
| $ A( I, L )*TEMP1 + B( I, L )*TEMP2 |
| 160 CONTINUE |
| END IF |
| 170 CONTINUE |
| 180 CONTINUE |
| END IF |
| ELSE |
| * |
| * Form C := alpha*A'*B + alpha*B'*A + C. |
| * |
| IF( UPPER )THEN |
| DO 210, J = 1, N |
| DO 200, I = 1, J |
| TEMP1 = ZERO |
| TEMP2 = ZERO |
| DO 190, L = 1, K |
| TEMP1 = TEMP1 + A( L, I )*B( L, J ) |
| TEMP2 = TEMP2 + B( L, I )*A( L, J ) |
| 190 CONTINUE |
| IF( BETA.EQ.ZERO )THEN |
| C( I, J ) = ALPHA*TEMP1 + ALPHA*TEMP2 |
| ELSE |
| C( I, J ) = BETA *C( I, J ) + |
| $ ALPHA*TEMP1 + ALPHA*TEMP2 |
| END IF |
| 200 CONTINUE |
| 210 CONTINUE |
| ELSE |
| DO 240, J = 1, N |
| DO 230, I = J, N |
| TEMP1 = ZERO |
| TEMP2 = ZERO |
| DO 220, L = 1, K |
| TEMP1 = TEMP1 + A( L, I )*B( L, J ) |
| TEMP2 = TEMP2 + B( L, I )*A( L, J ) |
| 220 CONTINUE |
| IF( BETA.EQ.ZERO )THEN |
| C( I, J ) = ALPHA*TEMP1 + ALPHA*TEMP2 |
| ELSE |
| C( I, J ) = BETA *C( I, J ) + |
| $ ALPHA*TEMP1 + ALPHA*TEMP2 |
| END IF |
| 230 CONTINUE |
| 240 CONTINUE |
| END IF |
| END IF |
| * |
| RETURN |
| * |
| * End of DSYR2K. |
| * |
| END |
| SUBROUTINE DSYRK ( UPLO, TRANS, N, K, ALPHA, A, LDA, |
| $ BETA, C, LDC ) |
| * .. Scalar Arguments .. |
| CHARACTER UPLO, TRANS |
| INTEGER N, K, LDA, LDC |
| DOUBLE PRECISION ALPHA, BETA |
| * .. Array Arguments .. |
| DOUBLE PRECISION A( LDA, * ), C( LDC, * ) |
| * .. |
| * |
| * Purpose |
| * ======= |
| * |
| * DSYRK performs one of the symmetric rank k operations |
| * |
| * C := alpha*A*A' + beta*C, |
| * |
| * or |
| * |
| * C := alpha*A'*A + beta*C, |
| * |
| * where alpha and beta are scalars, C is an n by n symmetric matrix |
| * and A is an n by k matrix in the first case and a k by n matrix |
| * in the second case. |
| * |
| * Parameters |
| * ========== |
| * |
| * UPLO - CHARACTER*1. |
| * On entry, UPLO specifies whether the upper or lower |
| * triangular part of the array C is to be referenced as |
| * follows: |
| * |
| * UPLO = 'U' or 'u' Only the upper triangular part of C |
| * is to be referenced. |
| * |
| * UPLO = 'L' or 'l' Only the lower triangular part of C |
| * is to be referenced. |
| * |
| * Unchanged on exit. |
| * |
| * TRANS - CHARACTER*1. |
| * On entry, TRANS specifies the operation to be performed as |
| * follows: |
| * |
| * TRANS = 'N' or 'n' C := alpha*A*A' + beta*C. |
| * |
| * TRANS = 'T' or 't' C := alpha*A'*A + beta*C. |
| * |
| * TRANS = 'C' or 'c' C := alpha*A'*A + beta*C. |
| * |
| * Unchanged on exit. |
| * |
| * N - INTEGER. |
| * On entry, N specifies the order of the matrix C. N must be |
| * at least zero. |
| * Unchanged on exit. |
| * |
| * K - INTEGER. |
| * On entry with TRANS = 'N' or 'n', K specifies the number |
| * of columns of the matrix A, and on entry with |
| * TRANS = 'T' or 't' or 'C' or 'c', K specifies the number |
| * of rows of the matrix A. 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 TRANS = 'N' or 'n', and is n otherwise. |
| * Before entry with TRANS = 'N' or 'n', the leading n by k |
| * part of the array A must contain the matrix A, otherwise |
| * the leading k by n 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 TRANS = 'N' or 'n' |
| * then LDA must be at least max( 1, n ), otherwise LDA must |
| * be at least max( 1, k ). |
| * Unchanged on exit. |
| * |
| * BETA - DOUBLE PRECISION. |
| * On entry, BETA specifies the scalar beta. |
| * Unchanged on exit. |
| * |
| * C - DOUBLE PRECISION array of DIMENSION ( LDC, n ). |
| * Before entry with UPLO = 'U' or 'u', the leading n by n |
| * upper triangular part of the array C must contain the upper |
| * triangular part of the symmetric matrix and the strictly |
| * lower triangular part of C is not referenced. On exit, the |
| * upper triangular part of the array C is overwritten by the |
| * upper triangular part of the updated matrix. |
| * Before entry with UPLO = 'L' or 'l', the leading n by n |
| * lower triangular part of the array C must contain the lower |
| * triangular part of the symmetric matrix and the strictly |
| * upper triangular part of C is not referenced. On exit, the |
| * lower triangular part of the array C is overwritten by the |
| * lower triangular part of the updated matrix. |
| * |
| * 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, n ). |
| * 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 UPPER |
| INTEGER I, INFO, J, L, NROWA |
| DOUBLE PRECISION TEMP |
| * .. Parameters .. |
| DOUBLE PRECISION ONE , ZERO |
| PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) |
| * .. |
| * .. Executable Statements .. |
| * |
| * Test the input parameters. |
| * |
| IF( LSAME( TRANS, 'N' ) )THEN |
| NROWA = N |
| ELSE |
| NROWA = K |
| END IF |
| UPPER = LSAME( UPLO, 'U' ) |
| * |
| INFO = 0 |
| IF( ( .NOT.UPPER ).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( N .LT.0 )THEN |
| INFO = 3 |
| ELSE IF( K .LT.0 )THEN |
| INFO = 4 |
| ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN |
| INFO = 7 |
| ELSE IF( LDC.LT.MAX( 1, N ) )THEN |
| INFO = 10 |
| END IF |
| IF( INFO.NE.0 )THEN |
| CALL XERBLA( 'DSYRK ', INFO ) |
| RETURN |
| END IF |
| * |
| * Quick return if possible. |
| * |
| IF( ( N.EQ.0 ).OR. |
| $ ( ( ( ALPHA.EQ.ZERO ).OR.( K.EQ.0 ) ).AND.( BETA.EQ.ONE ) ) ) |
| $ RETURN |
| * |
| * And when alpha.eq.zero. |
| * |
| IF( ALPHA.EQ.ZERO )THEN |
| IF( UPPER )THEN |
| IF( BETA.EQ.ZERO )THEN |
| DO 20, J = 1, N |
| DO 10, I = 1, J |
| C( I, J ) = ZERO |
| 10 CONTINUE |
| 20 CONTINUE |
| ELSE |
| DO 40, J = 1, N |
| DO 30, I = 1, J |
| C( I, J ) = BETA*C( I, J ) |
| 30 CONTINUE |
| 40 CONTINUE |
| END IF |
| ELSE |
| IF( BETA.EQ.ZERO )THEN |
| DO 60, J = 1, N |
| DO 50, I = J, N |
| C( I, J ) = ZERO |
| 50 CONTINUE |
| 60 CONTINUE |
| ELSE |
| DO 80, J = 1, N |
| DO 70, I = J, N |
| C( I, J ) = BETA*C( I, J ) |
| 70 CONTINUE |
| 80 CONTINUE |
| END IF |
| END IF |
| RETURN |
| END IF |
| * |
| * Start the operations. |
| * |
| IF( LSAME( TRANS, 'N' ) )THEN |
| * |
| * Form C := alpha*A*A' + beta*C. |
| * |
| IF( UPPER )THEN |
| DO 130, J = 1, N |
| IF( BETA.EQ.ZERO )THEN |
| DO 90, I = 1, J |
| C( I, J ) = ZERO |
| 90 CONTINUE |
| ELSE IF( BETA.NE.ONE )THEN |
| DO 100, I = 1, J |
| C( I, J ) = BETA*C( I, J ) |
| 100 CONTINUE |
| END IF |
| DO 120, L = 1, K |
| IF( A( J, L ).NE.ZERO )THEN |
| TEMP = ALPHA*A( J, L ) |
| DO 110, I = 1, J |
| C( I, J ) = C( I, J ) + TEMP*A( I, L ) |
| 110 CONTINUE |
| END IF |
| 120 CONTINUE |
| 130 CONTINUE |
| ELSE |
| DO 180, J = 1, N |
| IF( BETA.EQ.ZERO )THEN |
| DO 140, I = J, N |
| C( I, J ) = ZERO |
| 140 CONTINUE |
| ELSE IF( BETA.NE.ONE )THEN |
| DO 150, I = J, N |
| C( I, J ) = BETA*C( I, J ) |
| 150 CONTINUE |
| END IF |
| DO 170, L = 1, K |
| IF( A( J, L ).NE.ZERO )THEN |
| TEMP = ALPHA*A( J, L ) |
| DO 160, I = J, N |
| C( I, J ) = C( I, J ) + TEMP*A( I, L ) |
| 160 CONTINUE |
| END IF |
| 170 CONTINUE |
| 180 CONTINUE |
| END IF |
| ELSE |
| * |
| * Form C := alpha*A'*A + beta*C. |
| * |
| IF( UPPER )THEN |
| DO 210, J = 1, N |
| DO 200, I = 1, J |
| TEMP = ZERO |
| DO 190, L = 1, K |
| TEMP = TEMP + A( L, I )*A( L, J ) |
| 190 CONTINUE |
| IF( BETA.EQ.ZERO )THEN |
| C( I, J ) = ALPHA*TEMP |
| ELSE |
| C( I, J ) = ALPHA*TEMP + BETA*C( I, J ) |
| END IF |
| 200 CONTINUE |
| 210 CONTINUE |
| ELSE |
| DO 240, J = 1, N |
| DO 230, I = J, N |
| TEMP = ZERO |
| DO 220, L = 1, K |
| TEMP = TEMP + A( L, I )*A( L, J ) |
| 220 CONTINUE |
| IF( BETA.EQ.ZERO )THEN |
| C( I, J ) = ALPHA*TEMP |
| ELSE |
| C( I, J ) = ALPHA*TEMP + BETA*C( I, J ) |
| END IF |
| 230 CONTINUE |
| 240 CONTINUE |
| END IF |
| END IF |
| * |
| RETURN |
| * |
| * End of DSYRK . |
| * |
| END |
| SUBROUTINE DTBMV ( UPLO, TRANS, DIAG, N, K, A, LDA, X, INCX ) |
| * .. Scalar Arguments .. |
| INTEGER INCX, K, LDA, N |
| CHARACTER DIAG, TRANS, UPLO |
| * .. Array Arguments .. |
| DOUBLE PRECISION A( LDA, * ), X( * ) |
| * .. |
| * |
| * Purpose |
| * ======= |
| * |
| * DTBMV performs one of the matrix-vector operations |
| * |
| * x := A*x, or x := A'*x, |
| * |
| * where x is an n element vector and A is an n by n unit, or non-unit, |
| * upper or lower triangular band matrix, with ( k + 1 ) diagonals. |
| * |
| * 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 operation to be performed as |
| * follows: |
| * |
| * TRANS = 'N' or 'n' x := A*x. |
| * |
| * TRANS = 'T' or 't' x := A'*x. |
| * |
| * TRANS = 'C' or 'c' x := A'*x. |
| * |
| * 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. |
| * |
| * K - INTEGER. |
| * On entry with UPLO = 'U' or 'u', K specifies the number of |
| * super-diagonals of the matrix A. |
| * On entry with UPLO = 'L' or 'l', K specifies the number of |
| * sub-diagonals of the matrix A. |
| * K must satisfy 0 .le. K. |
| * Unchanged on exit. |
| * |
| * A - DOUBLE PRECISION array of DIMENSION ( LDA, n ). |
| * Before entry with UPLO = 'U' or 'u', the leading ( k + 1 ) |
| * by n part of the array A must contain the upper triangular |
| * band part of the matrix of coefficients, supplied column by |
| * column, with the leading diagonal of the matrix in row |
| * ( k + 1 ) of the array, the first super-diagonal starting at |
| * position 2 in row k, and so on. The top left k by k triangle |
| * of the array A is not referenced. |
| * The following program segment will transfer an upper |
| * triangular band matrix from conventional full matrix storage |
| * to band storage: |
| * |
| * DO 20, J = 1, N |
| * M = K + 1 - J |
| * DO 10, I = MAX( 1, J - K ), J |
| * A( M + I, J ) = matrix( I, J ) |
| * 10 CONTINUE |
| * 20 CONTINUE |
| * |
| * Before entry with UPLO = 'L' or 'l', the leading ( k + 1 ) |
| * by n part of the array A must contain the lower triangular |
| * band part of the matrix of coefficients, supplied column by |
| * column, with the leading diagonal of the matrix in row 1 of |
| * the array, the first sub-diagonal starting at position 1 in |
| * row 2, and so on. The bottom right k by k triangle of the |
| * array A is not referenced. |
| * The following program segment will transfer a lower |
| * triangular band matrix from conventional full matrix storage |
| * to band storage: |
| * |
| * DO 20, J = 1, N |
| * M = 1 - J |
| * DO 10, I = J, MIN( N, J + K ) |
| * A( M + I, J ) = matrix( I, J ) |
| * 10 CONTINUE |
| * 20 CONTINUE |
| * |
| * Note that when DIAG = 'U' or 'u' the elements of the array A |
| * corresponding to the diagonal elements of the matrix are not |
| * referenced, 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 |
| * ( k + 1 ). |
| * 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 vector x. On exit, X is overwritten with the |
| * tranformed 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, KPLUS1, KX, L |
| LOGICAL NOUNIT |
| * .. External Functions .. |
| LOGICAL LSAME |
| EXTERNAL LSAME |
| * .. External Subroutines .. |
| EXTERNAL XERBLA |
| * .. Intrinsic Functions .. |
| INTRINSIC MAX, MIN |
| * .. |
| * .. 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( K.LT.0 )THEN |
| INFO = 5 |
| ELSE IF( LDA.LT.( K + 1 ) )THEN |
| INFO = 7 |
| ELSE IF( INCX.EQ.0 )THEN |
| INFO = 9 |
| END IF |
| IF( INFO.NE.0 )THEN |
| CALL XERBLA( 'DTBMV ', 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 := A*x. |
| * |
| IF( LSAME( UPLO, 'U' ) )THEN |
| KPLUS1 = K + 1 |
| IF( INCX.EQ.1 )THEN |
| DO 20, J = 1, N |
| IF( X( J ).NE.ZERO )THEN |
| TEMP = X( J ) |
| L = KPLUS1 - J |
| DO 10, I = MAX( 1, J - K ), J - 1 |
| X( I ) = X( I ) + TEMP*A( L + I, J ) |
| 10 CONTINUE |
| IF( NOUNIT ) |
| $ X( J ) = X( J )*A( KPLUS1, J ) |
| END IF |
| 20 CONTINUE |
| ELSE |
| JX = KX |
| DO 40, J = 1, N |
| IF( X( JX ).NE.ZERO )THEN |
| TEMP = X( JX ) |
| IX = KX |
| L = KPLUS1 - J |
| DO 30, I = MAX( 1, J - K ), J - 1 |
| X( IX ) = X( IX ) + TEMP*A( L + I, J ) |
| IX = IX + INCX |
| 30 CONTINUE |
| IF( NOUNIT ) |
| $ X( JX ) = X( JX )*A( KPLUS1, J ) |
| END IF |
| JX = JX + INCX |
| IF( J.GT.K ) |
| $ KX = KX + INCX |
| 40 CONTINUE |
| END IF |
| ELSE |
| IF( INCX.EQ.1 )THEN |
| DO 60, J = N, 1, -1 |
| IF( X( J ).NE.ZERO )THEN |
| TEMP = X( J ) |
| L = 1 - J |
| DO 50, I = MIN( N, J + K ), J + 1, -1 |
| X( I ) = X( I ) + TEMP*A( L + I, J ) |
| 50 CONTINUE |
| IF( NOUNIT ) |
| $ X( J ) = X( J )*A( 1, J ) |
| END IF |
| 60 CONTINUE |
| ELSE |
| KX = KX + ( N - 1 )*INCX |
| JX = KX |
| DO 80, J = N, 1, -1 |
| IF( X( JX ).NE.ZERO )THEN |
| TEMP = X( JX ) |
| IX = KX |
| L = 1 - J |
| DO 70, I = MIN( N, J + K ), J + 1, -1 |
| X( IX ) = X( IX ) + TEMP*A( L + I, J ) |
| IX = IX - INCX |
| 70 CONTINUE |
| IF( NOUNIT ) |
| $ X( JX ) = X( JX )*A( 1, J ) |
| END IF |
| JX = JX - INCX |
| IF( ( N - J ).GE.K ) |
| $ KX = KX - INCX |
| 80 CONTINUE |
| END IF |
| END IF |
| ELSE |
| * |
| * Form x := A'*x. |
| * |
| IF( LSAME( UPLO, 'U' ) )THEN |
| KPLUS1 = K + 1 |
| IF( INCX.EQ.1 )THEN |
| DO 100, J = N, 1, -1 |
| TEMP = X( J ) |
| L = KPLUS1 - J |
| IF( NOUNIT ) |
| $ TEMP = TEMP*A( KPLUS1, J ) |
| DO 90, I = J - 1, MAX( 1, J - K ), -1 |
| TEMP = TEMP + A( L + I, J )*X( I ) |
| 90 CONTINUE |
| X( J ) = TEMP |
| 100 CONTINUE |
| ELSE |
| KX = KX + ( N - 1 )*INCX |
| JX = KX |
| DO 120, J = N, 1, -1 |
| TEMP = X( JX ) |
| KX = KX - INCX |
| IX = KX |
| L = KPLUS1 - J |
| IF( NOUNIT ) |
| $ TEMP = TEMP*A( KPLUS1, J ) |
| DO 110, I = J - 1, MAX( 1, J - K ), -1 |
| TEMP = TEMP + A( L + I, J )*X( IX ) |
| IX = IX - INCX |
| 110 CONTINUE |
| X( JX ) = TEMP |
| JX = JX - INCX |
| 120 CONTINUE |
| END IF |
| ELSE |
| IF( INCX.EQ.1 )THEN |
| DO 140, J = 1, N |
| TEMP = X( J ) |
| L = 1 - J |
| IF( NOUNIT ) |
| $ TEMP = TEMP*A( 1, J ) |
| DO 130, I = J + 1, MIN( N, J + K ) |
| TEMP = TEMP + A( L + I, J )*X( I ) |
| 130 CONTINUE |
| X( J ) = TEMP |
| 140 CONTINUE |
| ELSE |
| JX = KX |
| DO 160, J = 1, N |
| TEMP = X( JX ) |
| KX = KX + INCX |
| IX = KX |
| L = 1 - J |
| IF( NOUNIT ) |
| $ TEMP = TEMP*A( 1, J ) |
| DO 150, I = J + 1, MIN( N, J + K ) |
| TEMP = TEMP + A( L + I, J )*X( IX ) |
| IX = IX + INCX |
| 150 CONTINUE |
| X( JX ) = TEMP |
| JX = JX + INCX |
| 160 CONTINUE |
| END IF |
| END IF |
| END IF |
| * |
| RETURN |
| * |
| * End of DTBMV . |
| * |
| END |
| SUBROUTINE DTBSV ( UPLO, TRANS, DIAG, N, K, A, LDA, X, INCX ) |
| * .. Scalar Arguments .. |
| INTEGER INCX, K, LDA, N |
| CHARACTER DIAG, TRANS, UPLO |
| * .. Array Arguments .. |
| DOUBLE PRECISION A( LDA, * ), X( * ) |
| * .. |
| * |
| * Purpose |
| * ======= |
| * |
| * DTBSV 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 band matrix, with ( k + 1 ) |
| * diagonals. |
| * |
| * 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. |
| * |
| * K - INTEGER. |
| * On entry with UPLO = 'U' or 'u', K specifies the number of |
| * super-diagonals of the matrix A. |
| * On entry with UPLO = 'L' or 'l', K specifies the number of |
| * sub-diagonals of the matrix A. |
| * K must satisfy 0 .le. K. |
| * Unchanged on exit. |
| * |
| * A - DOUBLE PRECISION array of DIMENSION ( LDA, n ). |
| * Before entry with UPLO = 'U' or 'u', the leading ( k + 1 ) |
| * by n part of the array A must contain the upper triangular |
| * band part of the matrix of coefficients, supplied column by |
| * column, with the leading diagonal of the matrix in row |
| * ( k + 1 ) of the array, the first super-diagonal starting at |
| * position 2 in row k, and so on. The top left k by k triangle |
| * of the array A is not referenced. |
| * The following program segment will transfer an upper |
| * triangular band matrix from conventional full matrix storage |
| * to band storage: |
| * |
| * DO 20, J = 1, N |
| * M = K + 1 - J |
| * DO 10, I = MAX( 1, J - K ), J |
| * A( M + I, J ) = matrix( I, J ) |
| * 10 CONTINUE |
| * 20 CONTINUE |
| * |
| * Before entry with UPLO = 'L' or 'l', the leading ( k + 1 ) |
| * by n part of the array A must contain the lower triangular |
| * band part of the matrix of coefficients, supplied column by |
| * column, with the leading diagonal of the matrix in row 1 of |
| * the array, the first sub-diagonal starting at position 1 in |
| * row 2, and so on. The bottom right k by k triangle of the |
| * array A is not referenced. |
| * The following program segment will transfer a lower |
| * triangular band matrix from conventional full matrix storage |
| * to band storage: |
| * |
| * DO 20, J = 1, N |
| * M = 1 - J |
| * DO 10, I = J, MIN( N, J + K ) |
| * A( M + I, J ) = matrix( I, J ) |
| * 10 CONTINUE |
| * 20 CONTINUE |
| * |
| * Note that when DIAG = 'U' or 'u' the elements of the array A |
| * corresponding to the diagonal elements of the matrix are not |
| * referenced, 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 |
| * ( k + 1 ). |
| * 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, KPLUS1, KX, L |
| LOGICAL NOUNIT |
| * .. External Functions .. |
| LOGICAL LSAME |
| EXTERNAL LSAME |
| * .. External Subroutines .. |
| EXTERNAL XERBLA |
| * .. Intrinsic Functions .. |
| INTRINSIC MAX, MIN |
| * .. |
| * .. 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( K.LT.0 )THEN |
| INFO = 5 |
| ELSE IF( LDA.LT.( K + 1 ) )THEN |
| INFO = 7 |
| ELSE IF( INCX.EQ.0 )THEN |
| INFO = 9 |
| END IF |
| IF( INFO.NE.0 )THEN |
| CALL XERBLA( 'DTBSV ', 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 by sequentially with one pass through A. |
| * |
| IF( LSAME( TRANS, 'N' ) )THEN |
| * |
| * Form x := inv( A )*x. |
| * |
| IF( LSAME( UPLO, 'U' ) )THEN |
| KPLUS1 = K + 1 |
| IF( INCX.EQ.1 )THEN |
| DO 20, J = N, 1, -1 |
| IF( X( J ).NE.ZERO )THEN |
| L = KPLUS1 - J |
| IF( NOUNIT ) |
| $ X( J ) = X( J )/A( KPLUS1, J ) |
| TEMP = X( J ) |
| DO 10, I = J - 1, MAX( 1, J - K ), -1 |
| X( I ) = X( I ) - TEMP*A( L + I, J ) |
| 10 CONTINUE |
| END IF |
| 20 CONTINUE |
| ELSE |
| KX = KX + ( N - 1 )*INCX |
| JX = KX |
| DO 40, J = N, 1, -1 |
| KX = KX - INCX |
| IF( X( JX ).NE.ZERO )THEN |
| IX = KX |
| L = KPLUS1 - J |
| IF( NOUNIT ) |
| $ X( JX ) = X( JX )/A( KPLUS1, J ) |
| TEMP = X( JX ) |
| DO 30, I = J - 1, MAX( 1, J - K ), -1 |
| X( IX ) = X( IX ) - TEMP*A( L + I, J ) |
| IX = IX - INCX |
| 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 |
| L = 1 - J |
| IF( NOUNIT ) |
| $ X( J ) = X( J )/A( 1, J ) |
| TEMP = X( J ) |
| DO 50, I = J + 1, MIN( N, J + K ) |
| X( I ) = X( I ) - TEMP*A( L + I, J ) |
| 50 CONTINUE |
| END IF |
| 60 CONTINUE |
| ELSE |
| JX = KX |
| DO 80, J = 1, N |
| KX = KX + INCX |
| IF( X( JX ).NE.ZERO )THEN |
| IX = KX |
| L = 1 - J |
| IF( NOUNIT ) |
| $ X( JX ) = X( JX )/A( 1, J ) |
| TEMP = X( JX ) |
| DO 70, I = J + 1, MIN( N, J + K ) |
| X( IX ) = X( IX ) - TEMP*A( L + I, J ) |
| IX = IX + INCX |
| 70 CONTINUE |
| END IF |
| JX = JX + INCX |
| 80 CONTINUE |
| END IF |
| END IF |
| ELSE |
| * |
| * Form x := inv( A')*x. |
| * |
| IF( LSAME( UPLO, 'U' ) )THEN |
| KPLUS1 = K + 1 |
| IF( INCX.EQ.1 )THEN |
| DO 100, J = 1, N |
| TEMP = X( J ) |
| L = KPLUS1 - J |
| DO 90, I = MAX( 1, J - K ), J - 1 |
| TEMP = TEMP - A( L + I, J )*X( I ) |
| 90 CONTINUE |
| IF( NOUNIT ) |
| $ TEMP = TEMP/A( KPLUS1, J ) |
| X( J ) = TEMP |
| 100 CONTINUE |
| ELSE |
| JX = KX |
| DO 120, J = 1, N |
| TEMP = X( JX ) |
| IX = KX |
| L = KPLUS1 - J |
| DO 110, I = MAX( 1, J - K ), J - 1 |
| TEMP = TEMP - A( L + I, J )*X( IX ) |
| IX = IX + INCX |
| 110 CONTINUE |
| IF( NOUNIT ) |
| $ TEMP = TEMP/A( KPLUS1, J ) |
| X( JX ) = TEMP |
| JX = JX + INCX |
| IF( J.GT.K ) |
| $ KX = KX + INCX |
| 120 CONTINUE |
| END IF |
| ELSE |
| IF( INCX.EQ.1 )THEN |
| DO 140, J = N, 1, -1 |
| TEMP = X( J ) |
| L = 1 - J |
| DO 130, I = MIN( N, J + K ), J + 1, -1 |
| TEMP = TEMP - A( L + I, J )*X( I ) |
| 130 CONTINUE |
| IF( NOUNIT ) |
| $ TEMP = TEMP/A( 1, 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 |
| L = 1 - J |
| DO 150, I = MIN( N, J + K ), J + 1, -1 |
| TEMP = TEMP - A( L + I, J )*X( IX ) |
| IX = IX - INCX |
| 150 CONTINUE |
| IF( NOUNIT ) |
| $ TEMP = TEMP/A( 1, J ) |
| X( JX ) = TEMP |
| JX = JX - INCX |
| IF( ( N - J ).GE.K ) |
| $ KX = KX - INCX |
| 160 CONTINUE |
| END IF |
| END IF |
| END IF |
| * |
| RETURN |
| * |
| * End of DTBSV . |
| * |
| END |
| SUBROUTINE DTPMV ( UPLO, TRANS, DIAG, N, AP, X, INCX ) |
| * .. Scalar Arguments .. |
| INTEGER INCX, N |
| CHARACTER DIAG, TRANS, UPLO |
| * .. Array Arguments .. |
| DOUBLE PRECISION AP( * ), X( * ) |
| * .. |
| * |
| * Purpose |
| * ======= |
| * |
| * DTPMV performs one of the matrix-vector operations |
| * |
| * x := A*x, or x := A'*x, |
| * |
| * where x is an n element vector and A is an n by n unit, or non-unit, |
| * upper or lower triangular matrix, supplied in packed form. |
| * |
| * 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 operation to be performed as |
| * follows: |
| * |
| * TRANS = 'N' or 'n' x := A*x. |
| * |
| * TRANS = 'T' or 't' x := A'*x. |
| * |
| * TRANS = 'C' or 'c' x := A'*x. |
| * |
| * 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. |
| * |
| * AP - DOUBLE PRECISION array of DIMENSION at least |
| * ( ( n*( n + 1 ) )/2 ). |
| * Before entry with UPLO = 'U' or 'u', the array AP must |
| * contain the upper triangular matrix packed sequentially, |
| * column by column, so that AP( 1 ) contains a( 1, 1 ), |
| * AP( 2 ) and AP( 3 ) contain a( 1, 2 ) and a( 2, 2 ) |
| * respectively, and so on. |
| * Before entry with UPLO = 'L' or 'l', the array AP must |
| * contain the lower triangular matrix packed sequentially, |
| * column by column, so that AP( 1 ) contains a( 1, 1 ), |
| * AP( 2 ) and AP( 3 ) contain a( 2, 1 ) and a( 3, 1 ) |
| * respectively, and so on. |
| * Note that when DIAG = 'U' or 'u', the diagonal elements of |
| * A are not referenced, but are assumed to be unity. |
| * 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 vector x. On exit, X is overwritten with the |
| * tranformed 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, K, KK, KX |
| LOGICAL NOUNIT |
| * .. External Functions .. |
| LOGICAL LSAME |
| EXTERNAL LSAME |
| * .. External Subroutines .. |
| EXTERNAL XERBLA |
| * .. |
| * .. 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( INCX.EQ.0 )THEN |
| INFO = 7 |
| END IF |
| IF( INFO.NE.0 )THEN |
| CALL XERBLA( 'DTPMV ', 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 AP are |
| * accessed sequentially with one pass through AP. |
| * |
| IF( LSAME( TRANS, 'N' ) )THEN |
| * |
| * Form x:= A*x. |
| * |
| IF( LSAME( UPLO, 'U' ) )THEN |
| KK =1 |
| IF( INCX.EQ.1 )THEN |
| DO 20, J = 1, N |
| IF( X( J ).NE.ZERO )THEN |
| TEMP = X( J ) |
| K = KK |
| DO 10, I = 1, J - 1 |
| X( I ) = X( I ) + TEMP*AP( K ) |
| K = K + 1 |
| 10 CONTINUE |
| IF( NOUNIT ) |
| $ X( J ) = X( J )*AP( KK + J - 1 ) |
| END IF |
| KK = KK + J |
| 20 CONTINUE |
| ELSE |
| JX = KX |
| DO 40, J = 1, N |
| IF( X( JX ).NE.ZERO )THEN |
| TEMP = X( JX ) |
| IX = KX |
| DO 30, K = KK, KK + J - 2 |
| X( IX ) = X( IX ) + TEMP*AP( K ) |
| IX = IX + INCX |
| 30 CONTINUE |
| IF( NOUNIT ) |
| $ X( JX ) = X( JX )*AP( KK + J - 1 ) |
| END IF |
| JX = JX + INCX |
| KK = KK + J |
| 40 CONTINUE |
| END IF |
| ELSE |
| KK = ( N*( N + 1 ) )/2 |
| IF( INCX.EQ.1 )THEN |
| DO 60, J = N, 1, -1 |
| IF( X( J ).NE.ZERO )THEN |
| TEMP = X( J ) |
| K = KK |
| DO 50, I = N, J + 1, -1 |
| X( I ) = X( I ) + TEMP*AP( K ) |
| K = K - 1 |
| 50 CONTINUE |
| IF( NOUNIT ) |
| $ X( J ) = X( J )*AP( KK - N + J ) |
| END IF |
| KK = KK - ( N - J + 1 ) |
| 60 CONTINUE |
| ELSE |
| KX = KX + ( N - 1 )*INCX |
| JX = KX |
| DO 80, J = N, 1, -1 |
| IF( X( JX ).NE.ZERO )THEN |
| TEMP = X( JX ) |
| IX = KX |
| DO 70, K = KK, KK - ( N - ( J + 1 ) ), -1 |
| X( IX ) = X( IX ) + TEMP*AP( K ) |
| IX = IX - INCX |
| 70 CONTINUE |
| IF( NOUNIT ) |
| $ X( JX ) = X( JX )*AP( KK - N + J ) |
| END IF |
| JX = JX - INCX |
| KK = KK - ( N - J + 1 ) |
| 80 CONTINUE |
| END IF |
| END IF |
| ELSE |
| * |
| * Form x := A'*x. |
| * |
| IF( LSAME( UPLO, 'U' ) )THEN |
| KK = ( N*( N + 1 ) )/2 |
| IF( INCX.EQ.1 )THEN |
| DO 100, J = N, 1, -1 |
| TEMP = X( J ) |
| IF( NOUNIT ) |
| $ TEMP = TEMP*AP( KK ) |
| K = KK - 1 |
| DO 90, I = J - 1, 1, -1 |
| TEMP = TEMP + AP( K )*X( I ) |
| K = K - 1 |
| 90 CONTINUE |
| X( J ) = TEMP |
| KK = KK - J |
| 100 CONTINUE |
| ELSE |
| JX = KX + ( N - 1 )*INCX |
| DO 120, J = N, 1, -1 |
| TEMP = X( JX ) |
| IX = JX |
| IF( NOUNIT ) |
| $ TEMP = TEMP*AP( KK ) |
| DO 110, K = KK - 1, KK - J + 1, -1 |
| IX = IX - INCX |
| TEMP = TEMP + AP( K )*X( IX ) |
| 110 CONTINUE |
| X( JX ) = TEMP |
| JX = JX - INCX |
| KK = KK - J |
| 120 CONTINUE |
| END IF |
| ELSE |
| KK = 1 |
| IF( INCX.EQ.1 )THEN |
| DO 140, J = 1, N |
| TEMP = X( J ) |
| IF( NOUNIT ) |
| $ TEMP = TEMP*AP( KK ) |
| K = KK + 1 |
| DO 130, I = J + 1, N |
| TEMP = TEMP + AP( K )*X( I ) |
| K = K + 1 |
| 130 CONTINUE |
| X( J ) = TEMP |
| KK = KK + ( N - J + 1 ) |
| 140 CONTINUE |
| ELSE |
| JX = KX |
| DO 160, J = 1, N |
| TEMP = X( JX ) |
| IX = JX |
| IF( NOUNIT ) |
| $ TEMP = TEMP*AP( KK ) |
| DO 150, K = KK + 1, KK + N - J |
| IX = IX + INCX |
| TEMP = TEMP + AP( K )*X( IX ) |
| 150 CONTINUE |
| X( JX ) = TEMP |
| JX = JX + INCX |
| KK = KK + ( N - J + 1 ) |
| 160 CONTINUE |
| END IF |
| END IF |
| END IF |
| * |
| RETURN |
| * |
| * End of DTPMV . |
| * |
| END |
| SUBROUTINE DTPSV ( UPLO, TRANS, DIAG, N, AP, X, INCX ) |
| * .. Scalar Arguments .. |
| INTEGER INCX, N |
| CHARACTER DIAG, TRANS, UPLO |
| * .. Array Arguments .. |
| DOUBLE PRECISION AP( * ), X( * ) |
| * .. |
| * |
| * Purpose |
| * ======= |
| * |
| * DTPSV 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, supplied in packed form. |
| * |
| * 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. |
| * |
| * AP - DOUBLE PRECISION array of DIMENSION at least |
| * ( ( n*( n + 1 ) )/2 ). |
| * Before entry with UPLO = 'U' or 'u', the array AP must |
| * contain the upper triangular matrix packed sequentially, |
| * column by column, so that AP( 1 ) contains a( 1, 1 ), |
| * AP( 2 ) and AP( 3 ) contain a( 1, 2 ) and a( 2, 2 ) |
| * respectively, and so on. |
| * Before entry with UPLO = 'L' or 'l', the array AP must |
| * contain the lower triangular matrix packed sequentially, |
| * column by column, so that AP( 1 ) contains a( 1, 1 ), |
| * AP( 2 ) and AP( 3 ) contain a( 2, 1 ) and a( 3, 1 ) |
| * respectively, and so on. |
| * Note that when DIAG = 'U' or 'u', the diagonal elements of |
| * A are not referenced, but are assumed to be unity. |
| * 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, K, KK, KX |
| LOGICAL NOUNIT |
| * .. External Functions .. |
| LOGICAL LSAME |
| EXTERNAL LSAME |
| * .. External Subroutines .. |
| EXTERNAL XERBLA |
| * .. |
| * .. 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( INCX.EQ.0 )THEN |
| INFO = 7 |
| END IF |
| IF( INFO.NE.0 )THEN |
| CALL XERBLA( 'DTPSV ', 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 AP are |
| * accessed sequentially with one pass through AP. |
| * |
| IF( LSAME( TRANS, 'N' ) )THEN |
| * |
| * Form x := inv( A )*x. |
| * |
| IF( LSAME( UPLO, 'U' ) )THEN |
| KK = ( N*( N + 1 ) )/2 |
| IF( INCX.EQ.1 )THEN |
| DO 20, J = N, 1, -1 |
| IF( X( J ).NE.ZERO )THEN |
| IF( NOUNIT ) |
| $ X( J ) = X( J )/AP( KK ) |
| TEMP = X( J ) |
| K = KK - 1 |
| DO 10, I = J - 1, 1, -1 |
| X( I ) = X( I ) - TEMP*AP( K ) |
| K = K - 1 |
| 10 CONTINUE |
| END IF |
| KK = KK - J |
| 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 )/AP( KK ) |
| TEMP = X( JX ) |
| IX = JX |
| DO 30, K = KK - 1, KK - J + 1, -1 |
| IX = IX - INCX |
| X( IX ) = X( IX ) - TEMP*AP( K ) |
| 30 CONTINUE |
| END IF |
| JX = JX - INCX |
| KK = KK - J |
| 40 CONTINUE |
| END IF |
| ELSE |
| KK = 1 |
| IF( INCX.EQ.1 )THEN |
| DO 60, J = 1, N |
| IF( X( J ).NE.ZERO )THEN |
| IF( NOUNIT ) |
| $ X( J ) = X( J )/AP( KK ) |
| TEMP = X( J ) |
| K = KK + 1 |
| DO 50, I = J + 1, N |
| X( I ) = X( I ) - TEMP*AP( K ) |
| K = K + 1 |
| 50 CONTINUE |
| END IF |
| KK = KK + ( N - J + 1 ) |
| 60 CONTINUE |
| ELSE |
| JX = KX |
| DO 80, J = 1, N |
| IF( X( JX ).NE.ZERO )THEN |
| IF( NOUNIT ) |
| $ X( JX ) = X( JX )/AP( KK ) |
| TEMP = X( JX ) |
| IX = JX |
| DO 70, K = KK + 1, KK + N - J |
| IX = IX + INCX |
| X( IX ) = X( IX ) - TEMP*AP( K ) |
| 70 CONTINUE |
| END IF |
| JX = JX + INCX |
| KK = KK + ( N - J + 1 ) |
| 80 CONTINUE |
| END IF |
| END IF |
| ELSE |
| * |
| * Form x := inv( A' )*x. |
| * |
| IF( LSAME( UPLO, 'U' ) )THEN |
| KK = 1 |
| IF( INCX.EQ.1 )THEN |
| DO 100, J = 1, N |
| TEMP = X( J ) |
| K = KK |
| DO 90, I = 1, J - 1 |
| TEMP = TEMP - AP( K )*X( I ) |
| K = K + 1 |
| 90 CONTINUE |
| IF( NOUNIT ) |
| $ TEMP = TEMP/AP( KK + J - 1 ) |
| X( J ) = TEMP |
| KK = KK + J |
| 100 CONTINUE |
| ELSE |
| JX = KX |
| DO 120, J = 1, N |
| TEMP = X( JX ) |
| IX = KX |
| DO 110, K = KK, KK + J - 2 |
| TEMP = TEMP - AP( K )*X( IX ) |
| IX = IX + INCX |
| 110 CONTINUE |
| IF( NOUNIT ) |
| $ TEMP = TEMP/AP( KK + J - 1 ) |
| X( JX ) = TEMP |
| JX = JX + INCX |
| KK = KK + J |
| 120 CONTINUE |
| END IF |
| ELSE |
| KK = ( N*( N + 1 ) )/2 |
| IF( INCX.EQ.1 )THEN |
| DO 140, J = N, 1, -1 |
| TEMP = X( J ) |
| K = KK |
| DO 130, I = N, J + 1, -1 |
| TEMP = TEMP - AP( K )*X( I ) |
| K = K - 1 |
| 130 CONTINUE |
| IF( NOUNIT ) |
| $ TEMP = TEMP/AP( KK - N + J ) |
| X( J ) = TEMP |
| KK = KK - ( N - J + 1 ) |
| 140 CONTINUE |
| ELSE |
| KX = KX + ( N - 1 )*INCX |
| JX = KX |
| DO 160, J = N, 1, -1 |
| TEMP = X( JX ) |
| IX = KX |
| DO 150, K = KK, KK - ( N - ( J + 1 ) ), -1 |
| TEMP = TEMP - AP( K )*X( IX ) |
| IX = IX - INCX |
| 150 CONTINUE |
| IF( NOUNIT ) |
| $ TEMP = TEMP/AP( KK - N + J ) |
| X( JX ) = TEMP |
| JX = JX - INCX |
| KK = KK - (N - J + 1 ) |
| 160 CONTINUE |
| END IF |
| END IF |
| END IF |
| * |
| RETURN |
| * |
| * End of DTPSV . |
| * |
| END |
| SUBROUTINE DTRMM ( SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, |
| $ B, LDB ) |
| * .. Scalar Arguments .. |
| CHARACTER SIDE, UPLO, TRANSA, DIAG |
| INTEGER M, N, LDA, LDB |
| DOUBLE PRECISION ALPHA |
| * .. Array Arguments .. |
| DOUBLE PRECISION A( LDA, * ), B( LDB, * ) |
| * .. |
| * |
| * Purpose |
| * ======= |
| * |
| * DTRMM performs one of the matrix-matrix operations |
| * |
| * B := alpha*op( A )*B, or B := alpha*B*op( A ), |
| * |
| * where alpha is a scalar, B is an m by n matrix, 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'. |
| * |
| * Parameters |
| * ========== |
| * |
| * SIDE - CHARACTER*1. |
| * On entry, SIDE specifies whether op( A ) multiplies B from |
| * the left or right as follows: |
| * |
| * SIDE = 'L' or 'l' B := alpha*op( A )*B. |
| * |
| * SIDE = 'R' or 'r' B := alpha*B*op( A ). |
| * |
| * 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 matrix B, and on exit is overwritten by the |
| * transformed matrix. |
| * |
| * 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( 'DTRMM ', 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*A*B. |
| * |
| IF( UPPER )THEN |
| DO 50, J = 1, N |
| DO 40, K = 1, M |
| IF( B( K, J ).NE.ZERO )THEN |
| TEMP = ALPHA*B( K, J ) |
| DO 30, I = 1, K - 1 |
| B( I, J ) = B( I, J ) + TEMP*A( I, K ) |
| 30 CONTINUE |
| IF( NOUNIT ) |
| $ TEMP = TEMP*A( K, K ) |
| B( K, J ) = TEMP |
| END IF |
| 40 CONTINUE |
| 50 CONTINUE |
| ELSE |
| DO 80, J = 1, N |
| DO 70 K = M, 1, -1 |
| IF( B( K, J ).NE.ZERO )THEN |
| TEMP = ALPHA*B( K, J ) |
| B( K, J ) = TEMP |
| IF( NOUNIT ) |
| $ B( K, J ) = B( K, J )*A( K, K ) |
| DO 60, I = K + 1, M |
| B( I, J ) = B( I, J ) + TEMP*A( I, K ) |
| 60 CONTINUE |
| END IF |
| 70 CONTINUE |
| 80 CONTINUE |
| END IF |
| ELSE |
| * |
| * Form B := alpha*A'*B. |
| * |
| IF( UPPER )THEN |
| DO 110, J = 1, N |
| DO 100, I = M, 1, -1 |
| TEMP = B( I, J ) |
| IF( NOUNIT ) |
| $ TEMP = TEMP*A( I, I ) |
| DO 90, K = 1, I - 1 |
| TEMP = TEMP + A( K, I )*B( K, J ) |
| 90 CONTINUE |
| B( I, J ) = ALPHA*TEMP |
| 100 CONTINUE |
| 110 CONTINUE |
| ELSE |
| DO 140, J = 1, N |
| DO 130, I = 1, M |
| TEMP = B( I, J ) |
| IF( NOUNIT ) |
| $ TEMP = TEMP*A( I, I ) |
| DO 120, K = I + 1, M |
| TEMP = TEMP + A( K, I )*B( K, J ) |
| 120 CONTINUE |
| B( I, J ) = ALPHA*TEMP |
| 130 CONTINUE |
| 140 CONTINUE |
| END IF |
| END IF |
| ELSE |
| IF( LSAME( TRANSA, 'N' ) )THEN |
| * |
| * Form B := alpha*B*A. |
| * |
| IF( UPPER )THEN |
| DO 180, J = N, 1, -1 |
| TEMP = ALPHA |
| IF( NOUNIT ) |
| $ TEMP = TEMP*A( J, J ) |
| DO 150, I = 1, M |
| B( I, J ) = TEMP*B( I, J ) |
| 150 CONTINUE |
| DO 170, K = 1, J - 1 |
| IF( A( K, J ).NE.ZERO )THEN |
| TEMP = ALPHA*A( K, J ) |
| DO 160, I = 1, M |
| B( I, J ) = B( I, J ) + TEMP*B( I, K ) |
| 160 CONTINUE |
| END IF |
| 170 CONTINUE |
| 180 CONTINUE |
| ELSE |
| DO 220, J = 1, N |
| TEMP = ALPHA |
| IF( NOUNIT ) |
| $ TEMP = TEMP*A( J, J ) |
| DO 190, I = 1, M |
| B( I, J ) = TEMP*B( I, J ) |
| 190 CONTINUE |
| DO 210, K = J + 1, N |
| IF( A( K, J ).NE.ZERO )THEN |
| TEMP = ALPHA*A( K, J ) |
| DO 200, I = 1, M |
| B( I, J ) = B( I, J ) + TEMP*B( I, K ) |
| 200 CONTINUE |
| END IF |
| 210 CONTINUE |
| 220 CONTINUE |
| END IF |
| ELSE |
| * |
| * Form B := alpha*B*A'. |
| * |
| IF( UPPER )THEN |
| DO 260, K = 1, N |
| DO 240, J = 1, K - 1 |
| IF( A( J, K ).NE.ZERO )THEN |
| TEMP = ALPHA*A( J, K ) |
| DO 230, I = 1, M |
| B( I, J ) = B( I, J ) + TEMP*B( I, K ) |
| 230 CONTINUE |
| END IF |
| 240 CONTINUE |
| TEMP = ALPHA |
| IF( NOUNIT ) |
| $ TEMP = TEMP*A( K, K ) |
| IF( TEMP.NE.ONE )THEN |
| DO 250, I = 1, M |
| B( I, K ) = TEMP*B( I, K ) |
| 250 CONTINUE |
| END IF |
| 260 CONTINUE |
| ELSE |
| DO 300, K = N, 1, -1 |
| DO 280, J = K + 1, N |
| IF( A( J, K ).NE.ZERO )THEN |
| TEMP = ALPHA*A( J, K ) |
| DO 270, I = 1, M |
| B( I, J ) = B( I, J ) + TEMP*B( I, K ) |
| 270 CONTINUE |
| END IF |
| 280 CONTINUE |
| TEMP = ALPHA |
| IF( NOUNIT ) |
| $ TEMP = TEMP*A( K, K ) |
| IF( TEMP.NE.ONE )THEN |
| DO 290, I = 1, M |
| B( I, K ) = TEMP*B( I, K ) |
| 290 CONTINUE |
| END IF |
| 300 CONTINUE |
| END IF |
| END IF |
| END IF |
| * |
| RETURN |
| * |
| * End of DTRMM . |
| * |
| END |
| SUBROUTINE DTRMV ( UPLO, TRANS, DIAG, N, A, LDA, X, INCX ) |
| * .. Scalar Arguments .. |
| INTEGER INCX, LDA, N |
| CHARACTER DIAG, TRANS, UPLO |
| * .. Array Arguments .. |
| DOUBLE PRECISION A( LDA, * ), X( * ) |
| * .. |
| * |
| * Purpose |
| * ======= |
| * |
| * DTRMV performs one of the matrix-vector operations |
| * |
| * x := A*x, or x := A'*x, |
| * |
| * where x is an n element vector and A is an n by n unit, or non-unit, |
| * upper or lower triangular matrix. |
| * |
| * 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 operation to be performed as |
| * follows: |
| * |
| * TRANS = 'N' or 'n' x := A*x. |
| * |
| * TRANS = 'T' or 't' x := A'*x. |
| * |
| * TRANS = 'C' or 'c' x := A'*x. |
| * |
| * 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 vector x. On exit, X is overwritten with the |
| * tranformed 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( 'DTRMV ', 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 := A*x. |
| * |
| IF( LSAME( UPLO, 'U' ) )THEN |
| IF( INCX.EQ.1 )THEN |
| DO 20, J = 1, N |
| IF( X( J ).NE.ZERO )THEN |
| TEMP = X( J ) |
| DO 10, I = 1, J - 1 |
| X( I ) = X( I ) + TEMP*A( I, J ) |
| 10 CONTINUE |
| IF( NOUNIT ) |
| $ X( J ) = X( J )*A( J, J ) |
| END IF |
| 20 CONTINUE |
| ELSE |
| JX = KX |
| DO 40, J = 1, N |
| IF( X( JX ).NE.ZERO )THEN |
| TEMP = X( JX ) |
| IX = KX |
| DO 30, I = 1, J - 1 |
| X( IX ) = X( IX ) + TEMP*A( I, J ) |
| IX = IX + INCX |
| 30 CONTINUE |
| IF( NOUNIT ) |
| $ X( JX ) = X( JX )*A( J, J ) |
| END IF |
| JX = JX + INCX |
| 40 CONTINUE |
| END IF |
| ELSE |
| IF( INCX.EQ.1 )THEN |
| DO 60, J = N, 1, -1 |
| IF( X( J ).NE.ZERO )THEN |
| TEMP = X( J ) |
| DO 50, I = N, J + 1, -1 |
| X( I ) = X( I ) + TEMP*A( I, J ) |
| 50 CONTINUE |
| IF( NOUNIT ) |
| $ X( J ) = X( J )*A( J, J ) |
| END IF |
| 60 CONTINUE |
| ELSE |
| KX = KX + ( N - 1 )*INCX |
| JX = KX |
| DO 80, J = N, 1, -1 |
| IF( X( JX ).NE.ZERO )THEN |
| TEMP = X( JX ) |
| IX = KX |
| DO 70, I = N, J + 1, -1 |
| X( IX ) = X( IX ) + TEMP*A( I, J ) |
| IX = IX - INCX |
| 70 CONTINUE |
| IF( NOUNIT ) |
| $ X( JX ) = X( JX )*A( J, J ) |
| END IF |
| JX = JX - INCX |
| 80 CONTINUE |
| END IF |
| END IF |
| ELSE |
| * |
| * Form x := A'*x. |
| * |
| IF( LSAME( UPLO, 'U' ) )THEN |
| IF( INCX.EQ.1 )THEN |
| DO 100, J = N, 1, -1 |
| TEMP = X( J ) |
| IF( NOUNIT ) |
| $ TEMP = TEMP*A( J, J ) |
| DO 90, I = J - 1, 1, -1 |
| TEMP = TEMP + A( I, J )*X( I ) |
| 90 CONTINUE |
| X( J ) = TEMP |
| 100 CONTINUE |
| ELSE |
| JX = KX + ( N - 1 )*INCX |
| DO 120, J = N, 1, -1 |
| TEMP = X( JX ) |
| IX = JX |
| IF( NOUNIT ) |
| $ TEMP = TEMP*A( J, J ) |
| DO 110, I = J - 1, 1, -1 |
| IX = IX - INCX |
| TEMP = TEMP + A( I, J )*X( IX ) |
| 110 CONTINUE |
| X( JX ) = TEMP |
| JX = JX - INCX |
| 120 CONTINUE |
| END IF |
| ELSE |
| IF( INCX.EQ.1 )THEN |
| DO 140, J = 1, N |
| TEMP = X( J ) |
| IF( NOUNIT ) |
| $ TEMP = TEMP*A( J, J ) |
| DO 130, I = J + 1, N |
| TEMP = TEMP + A( I, J )*X( I ) |
| 130 CONTINUE |
| X( J ) = TEMP |
| 140 CONTINUE |
| ELSE |
| JX = KX |
| DO 160, J = 1, N |
| TEMP = X( JX ) |
| IX = JX |
| IF( NOUNIT ) |
| $ TEMP = TEMP*A( J, J ) |
| DO 150, I = J + 1, N |
| IX = IX + INCX |
| TEMP = TEMP + A( I, J )*X( IX ) |
| 150 CONTINUE |
| X( JX ) = TEMP |
| JX = JX + INCX |
| 160 CONTINUE |
| END IF |
| END IF |
| END IF |
| * |
| RETURN |
| * |
| * End of DTRMV . |
| * |
| END |
| SUBROUTINE DTRSM ( SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, |
| $ B, LDB ) |
| * .. Scalar Arguments .. |
| CHARACTER 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 |
| SUBROUTINE DTRSV ( UPLO, TRANS, DIAG, N, A, LDA, X, INCX ) |
| * .. Scalar Arguments .. |
| INTEGER INCX, LDA, N |
| CHARACTER 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 |
| 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 modified 3/93 to return if incx .le. 0. |
| c modified 12/3/93, array(1) declarations changed to array(*) |
| c |
| double precision dx(*),dmax |
| integer i,incx,ix,n |
| c |
| idamax = 0 |
| if( n.lt.1 .or. incx.le.0 ) 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 |
| LOGICAL FUNCTION LSAME( CA, CB ) |
| * |
| * -- LAPACK auxiliary routine (version 2.0) -- |
| * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., |
| * Courant Institute, Argonne National Lab, and Rice University |
| * January 31, 1994 |
| * |
| * .. Scalar Arguments .. |
| CHARACTER CA, CB |
| * .. |
| * |
| * Purpose |
| * ======= |
| * |
| * LSAME returns .TRUE. if CA is the same letter as CB regardless of |
| * case. |
| * |
| * Arguments |
| * ========= |
| * |
| * CA (input) CHARACTER*1 |
| * CB (input) CHARACTER*1 |
| * CA and CB specify the single characters to be compared. |
| * |
| * ===================================================================== |
| * |
| * .. Intrinsic Functions .. |
| INTRINSIC ICHAR |
| * .. |
| * .. Local Scalars .. |
| INTEGER INTA, INTB, ZCODE |
| * .. |
| * .. Executable Statements .. |
| * |
| * Test if the characters are equal |
| * |
| LSAME = CA.EQ.CB |
| IF( LSAME ) |
| $ RETURN |
| * |
| * Now test for equivalence if both characters are alphabetic. |
| * |
| ZCODE = ICHAR( 'Z' ) |
| * |
| * Use 'Z' rather than 'A' so that ASCII can be detected on Prime |
| * machines, on which ICHAR returns a value with bit 8 set. |
| * ICHAR('A') on Prime machines returns 193 which is the same as |
| * ICHAR('A') on an EBCDIC machine. |
| * |
| INTA = ICHAR( CA ) |
| INTB = ICHAR( CB ) |
| * |
| IF( ZCODE.EQ.90 .OR. ZCODE.EQ.122 ) THEN |
| * |
| * ASCII is assumed - ZCODE is the ASCII code of either lower or |
| * upper case 'Z'. |
| * |
| IF( INTA.GE.97 .AND. INTA.LE.122 ) INTA = INTA - 32 |
| IF( INTB.GE.97 .AND. INTB.LE.122 ) INTB = INTB - 32 |
| * |
| ELSE IF( ZCODE.EQ.233 .OR. ZCODE.EQ.169 ) THEN |
| * |
| * EBCDIC is assumed - ZCODE is the EBCDIC code of either lower or |
| * upper case 'Z'. |
| * |
| IF( INTA.GE.129 .AND. INTA.LE.137 .OR. |
| $ INTA.GE.145 .AND. INTA.LE.153 .OR. |
| $ INTA.GE.162 .AND. INTA.LE.169 ) INTA = INTA + 64 |
| IF( INTB.GE.129 .AND. INTB.LE.137 .OR. |
| $ INTB.GE.145 .AND. INTB.LE.153 .OR. |
| $ INTB.GE.162 .AND. INTB.LE.169 ) INTB = INTB + 64 |
| * |
| ELSE IF( ZCODE.EQ.218 .OR. ZCODE.EQ.250 ) THEN |
| * |
| * ASCII is assumed, on Prime machines - ZCODE is the ASCII code |
| * plus 128 of either lower or upper case 'Z'. |
| * |
| IF( INTA.GE.225 .AND. INTA.LE.250 ) INTA = INTA - 32 |
| IF( INTB.GE.225 .AND. INTB.LE.250 ) INTB = INTB - 32 |
| END IF |
| LSAME = INTA.EQ.INTB |
| * |
| * RETURN |
| * |
| * End of LSAME |
| * |
| END |
