summaryrefslogtreecommitdiff
path: root/src/lib/lapack/dspgst.f
diff options
context:
space:
mode:
authorjofret2009-04-28 07:17:00 +0000
committerjofret2009-04-28 07:17:00 +0000
commit8c8d2f518968ce7057eec6aa5cd5aec8faab861a (patch)
tree3dd1788b71d6a3ce2b73d2d475a3133580e17530 /src/lib/lapack/dspgst.f
parent9f652ffc16a310ac6641a9766c5b9e2671e0e9cb (diff)
downloadscilab2c-8c8d2f518968ce7057eec6aa5cd5aec8faab861a.tar.gz
scilab2c-8c8d2f518968ce7057eec6aa5cd5aec8faab861a.tar.bz2
scilab2c-8c8d2f518968ce7057eec6aa5cd5aec8faab861a.zip
Moving lapack to right place
Diffstat (limited to 'src/lib/lapack/dspgst.f')
-rw-r--r--src/lib/lapack/dspgst.f208
1 files changed, 0 insertions, 208 deletions
diff --git a/src/lib/lapack/dspgst.f b/src/lib/lapack/dspgst.f
deleted file mode 100644
index 8e121a94..00000000
--- a/src/lib/lapack/dspgst.f
+++ /dev/null
@@ -1,208 +0,0 @@
- SUBROUTINE DSPGST( ITYPE, UPLO, N, AP, BP, INFO )
-*
-* -- LAPACK routine (version 3.1) --
-* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
-* November 2006
-*
-* .. Scalar Arguments ..
- CHARACTER UPLO
- INTEGER INFO, ITYPE, N
-* ..
-* .. Array Arguments ..
- DOUBLE PRECISION AP( * ), BP( * )
-* ..
-*
-* Purpose
-* =======
-*
-* DSPGST reduces a real symmetric-definite generalized eigenproblem
-* to standard form, using packed storage.
-*
-* If ITYPE = 1, the problem is A*x = lambda*B*x,
-* and A is overwritten by inv(U**T)*A*inv(U) or inv(L)*A*inv(L**T)
-*
-* If ITYPE = 2 or 3, the problem is A*B*x = lambda*x or
-* B*A*x = lambda*x, and A is overwritten by U*A*U**T or L**T*A*L.
-*
-* B must have been previously factorized as U**T*U or L*L**T by DPPTRF.
-*
-* Arguments
-* =========
-*
-* ITYPE (input) INTEGER
-* = 1: compute inv(U**T)*A*inv(U) or inv(L)*A*inv(L**T);
-* = 2 or 3: compute U*A*U**T or L**T*A*L.
-*
-* UPLO (input) CHARACTER*1
-* = 'U': Upper triangle of A is stored and B is factored as
-* U**T*U;
-* = 'L': Lower triangle of A is stored and B is factored as
-* L*L**T.
-*
-* N (input) INTEGER
-* The order of the matrices A and B. N >= 0.
-*
-* AP (input/output) DOUBLE PRECISION array, dimension (N*(N+1)/2)
-* On entry, the upper or lower triangle of the symmetric matrix
-* A, packed columnwise in a linear array. The j-th column of A
-* is stored in the array AP as follows:
-* if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
-* if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.
-*
-* On exit, if INFO = 0, the transformed matrix, stored in the
-* same format as A.
-*
-* BP (input) DOUBLE PRECISION array, dimension (N*(N+1)/2)
-* The triangular factor from the Cholesky factorization of B,
-* stored in the same format as A, as returned by DPPTRF.
-*
-* INFO (output) INTEGER
-* = 0: successful exit
-* < 0: if INFO = -i, the i-th argument had an illegal value
-*
-* =====================================================================
-*
-* .. Parameters ..
- DOUBLE PRECISION ONE, HALF
- PARAMETER ( ONE = 1.0D0, HALF = 0.5D0 )
-* ..
-* .. Local Scalars ..
- LOGICAL UPPER
- INTEGER J, J1, J1J1, JJ, K, K1, K1K1, KK
- DOUBLE PRECISION AJJ, AKK, BJJ, BKK, CT
-* ..
-* .. External Subroutines ..
- EXTERNAL DAXPY, DSCAL, DSPMV, DSPR2, DTPMV, DTPSV,
- $ XERBLA
-* ..
-* .. External Functions ..
- LOGICAL LSAME
- DOUBLE PRECISION DDOT
- EXTERNAL LSAME, DDOT
-* ..
-* .. Executable Statements ..
-*
-* Test the input parameters.
-*
- INFO = 0
- UPPER = LSAME( UPLO, 'U' )
- IF( ITYPE.LT.1 .OR. ITYPE.GT.3 ) THEN
- INFO = -1
- ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
- INFO = -2
- ELSE IF( N.LT.0 ) THEN
- INFO = -3
- END IF
- IF( INFO.NE.0 ) THEN
- CALL XERBLA( 'DSPGST', -INFO )
- RETURN
- END IF
-*
- IF( ITYPE.EQ.1 ) THEN
- IF( UPPER ) THEN
-*
-* Compute inv(U')*A*inv(U)
-*
-* J1 and JJ are the indices of A(1,j) and A(j,j)
-*
- JJ = 0
- DO 10 J = 1, N
- J1 = JJ + 1
- JJ = JJ + J
-*
-* Compute the j-th column of the upper triangle of A
-*
- BJJ = BP( JJ )
- CALL DTPSV( UPLO, 'Transpose', 'Nonunit', J, BP,
- $ AP( J1 ), 1 )
- CALL DSPMV( UPLO, J-1, -ONE, AP, BP( J1 ), 1, ONE,
- $ AP( J1 ), 1 )
- CALL DSCAL( J-1, ONE / BJJ, AP( J1 ), 1 )
- AP( JJ ) = ( AP( JJ )-DDOT( J-1, AP( J1 ), 1, BP( J1 ),
- $ 1 ) ) / BJJ
- 10 CONTINUE
- ELSE
-*
-* Compute inv(L)*A*inv(L')
-*
-* KK and K1K1 are the indices of A(k,k) and A(k+1,k+1)
-*
- KK = 1
- DO 20 K = 1, N
- K1K1 = KK + N - K + 1
-*
-* Update the lower triangle of A(k:n,k:n)
-*
- AKK = AP( KK )
- BKK = BP( KK )
- AKK = AKK / BKK**2
- AP( KK ) = AKK
- IF( K.LT.N ) THEN
- CALL DSCAL( N-K, ONE / BKK, AP( KK+1 ), 1 )
- CT = -HALF*AKK
- CALL DAXPY( N-K, CT, BP( KK+1 ), 1, AP( KK+1 ), 1 )
- CALL DSPR2( UPLO, N-K, -ONE, AP( KK+1 ), 1,
- $ BP( KK+1 ), 1, AP( K1K1 ) )
- CALL DAXPY( N-K, CT, BP( KK+1 ), 1, AP( KK+1 ), 1 )
- CALL DTPSV( UPLO, 'No transpose', 'Non-unit', N-K,
- $ BP( K1K1 ), AP( KK+1 ), 1 )
- END IF
- KK = K1K1
- 20 CONTINUE
- END IF
- ELSE
- IF( UPPER ) THEN
-*
-* Compute U*A*U'
-*
-* K1 and KK are the indices of A(1,k) and A(k,k)
-*
- KK = 0
- DO 30 K = 1, N
- K1 = KK + 1
- KK = KK + K
-*
-* Update the upper triangle of A(1:k,1:k)
-*
- AKK = AP( KK )
- BKK = BP( KK )
- CALL DTPMV( UPLO, 'No transpose', 'Non-unit', K-1, BP,
- $ AP( K1 ), 1 )
- CT = HALF*AKK
- CALL DAXPY( K-1, CT, BP( K1 ), 1, AP( K1 ), 1 )
- CALL DSPR2( UPLO, K-1, ONE, AP( K1 ), 1, BP( K1 ), 1,
- $ AP )
- CALL DAXPY( K-1, CT, BP( K1 ), 1, AP( K1 ), 1 )
- CALL DSCAL( K-1, BKK, AP( K1 ), 1 )
- AP( KK ) = AKK*BKK**2
- 30 CONTINUE
- ELSE
-*
-* Compute L'*A*L
-*
-* JJ and J1J1 are the indices of A(j,j) and A(j+1,j+1)
-*
- JJ = 1
- DO 40 J = 1, N
- J1J1 = JJ + N - J + 1
-*
-* Compute the j-th column of the lower triangle of A
-*
- AJJ = AP( JJ )
- BJJ = BP( JJ )
- AP( JJ ) = AJJ*BJJ + DDOT( N-J, AP( JJ+1 ), 1,
- $ BP( JJ+1 ), 1 )
- CALL DSCAL( N-J, BJJ, AP( JJ+1 ), 1 )
- CALL DSPMV( UPLO, N-J, ONE, AP( J1J1 ), BP( JJ+1 ), 1,
- $ ONE, AP( JJ+1 ), 1 )
- CALL DTPMV( UPLO, 'Transpose', 'Non-unit', N-J+1,
- $ BP( JJ ), AP( JJ ), 1 )
- JJ = J1J1
- 40 CONTINUE
- END IF
- END IF
- RETURN
-*
-* End of DSPGST
-*
- END