Original Version

1 files changed, 234 insertions, 0 deletions

diff --git a/2.3-1/src/fortran/lapack/zgeqpf.f b/2.3-1/src/fortran/lapack/zgeqpf.f
new file mode 100644
index 00000000..6d4f86f0
--- /dev/null
+++ b/2.3-1/src/fortran/lapack/zgeqpf.f
@@ -0,0 +1,234 @@
+      SUBROUTINE ZGEQPF( M, N, A, LDA, JPVT, TAU, WORK, RWORK, INFO )
+*
+*  -- LAPACK deprecated driver routine (version 3.1) --
+*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
+*     November 2006
+*
+*     .. Scalar Arguments ..
+      INTEGER            INFO, LDA, M, N
+*     ..
+*     .. Array Arguments ..
+      INTEGER            JPVT( * )
+      DOUBLE PRECISION   RWORK( * )
+      COMPLEX*16         A( LDA, * ), TAU( * ), WORK( * )
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  This routine is deprecated and has been replaced by routine ZGEQP3.
+*
+*  ZGEQPF computes a QR factorization with column pivoting of a
+*  complex M-by-N matrix A: A*P = Q*R.
+*
+*  Arguments
+*  =========
+*
+*  M       (input) INTEGER
+*          The number of rows of the matrix A. M >= 0.
+*
+*  N       (input) INTEGER
+*          The number of columns of the matrix A. N >= 0
+*
+*  A       (input/output) COMPLEX*16 array, dimension (LDA,N)
+*          On entry, the M-by-N matrix A.
+*          On exit, the upper triangle of the array contains the
+*          min(M,N)-by-N upper triangular matrix R; the elements
+*          below the diagonal, together with the array TAU,
+*          represent the unitary matrix Q as a product of
+*          min(m,n) elementary reflectors.
+*
+*  LDA     (input) INTEGER
+*          The leading dimension of the array A. LDA >= max(1,M).
+*
+*  JPVT    (input/output) INTEGER array, dimension (N)
+*          On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted
+*          to the front of A*P (a leading column); if JPVT(i) = 0,
+*          the i-th column of A is a free column.
+*          On exit, if JPVT(i) = k, then the i-th column of A*P
+*          was the k-th column of A.
+*
+*  TAU     (output) COMPLEX*16 array, dimension (min(M,N))
+*          The scalar factors of the elementary reflectors.
+*
+*  WORK    (workspace) COMPLEX*16 array, dimension (N)
+*
+*  RWORK   (workspace) DOUBLE PRECISION array, dimension (2*N)
+*
+*  INFO    (output) INTEGER
+*          = 0:  successful exit
+*          < 0:  if INFO = -i, the i-th argument had an illegal value
+*
+*  Further Details
+*  ===============
+*
+*  The matrix Q is represented as a product of elementary reflectors
+*
+*     Q = H(1) H(2) . . . H(n)
+*
+*  Each H(i) has the form
+*
+*     H = I - tau * v * v'
+*
+*  where tau is a complex scalar, and v is a complex vector with
+*  v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in A(i+1:m,i).
+*
+*  The matrix P is represented in jpvt as follows: If
+*     jpvt(j) = i
+*  then the jth column of P is the ith canonical unit vector.
+*
+*  Partial column norm updating strategy modified by
+*    Z. Drmac and Z. Bujanovic, Dept. of Mathematics,
+*    University of Zagreb, Croatia.
+*    June 2006.
+*  For more details see LAPACK Working Note 176.
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      DOUBLE PRECISION   ZERO, ONE
+      PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
+*     ..
+*     .. Local Scalars ..
+      INTEGER            I, ITEMP, J, MA, MN, PVT
+      DOUBLE PRECISION   TEMP, TEMP2, TOL3Z
+      COMPLEX*16         AII
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           XERBLA, ZGEQR2, ZLARF, ZLARFG, ZSWAP, ZUNM2R
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          ABS, DCMPLX, DCONJG, MAX, MIN, SQRT
+*     ..
+*     .. External Functions ..
+      INTEGER            IDAMAX
+      DOUBLE PRECISION   DLAMCH, DZNRM2
+      EXTERNAL           IDAMAX, DLAMCH, DZNRM2
+*     ..
+*     .. Executable Statements ..
+*
+*     Test the input arguments
+*
+      INFO = 0
+      IF( M.LT.0 ) THEN
+         INFO = -1
+      ELSE IF( N.LT.0 ) THEN
+         INFO = -2
+      ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
+         INFO = -4
+      END IF
+      IF( INFO.NE.0 ) THEN
+         CALL XERBLA( 'ZGEQPF', -INFO )
+         RETURN
+      END IF
+*
+      MN = MIN( M, N )
+      TOL3Z = SQRT(DLAMCH('Epsilon'))
+*
+*     Move initial columns up front
+*
+      ITEMP = 1
+      DO 10 I = 1, N
+         IF( JPVT( I ).NE.0 ) THEN
+            IF( I.NE.ITEMP ) THEN
+               CALL ZSWAP( M, A( 1, I ), 1, A( 1, ITEMP ), 1 )
+               JPVT( I ) = JPVT( ITEMP )
+               JPVT( ITEMP ) = I
+            ELSE
+               JPVT( I ) = I
+            END IF
+            ITEMP = ITEMP + 1
+         ELSE
+            JPVT( I ) = I
+         END IF
+   10 CONTINUE
+      ITEMP = ITEMP - 1
+*
+*     Compute the QR factorization and update remaining columns
+*
+      IF( ITEMP.GT.0 ) THEN
+         MA = MIN( ITEMP, M )
+         CALL ZGEQR2( M, MA, A, LDA, TAU, WORK, INFO )
+         IF( MA.LT.N ) THEN
+            CALL ZUNM2R( 'Left', 'Conjugate transpose', M, N-MA, MA, A,
+     $                   LDA, TAU, A( 1, MA+1 ), LDA, WORK, INFO )
+         END IF
+      END IF
+*
+      IF( ITEMP.LT.MN ) THEN
+*
+*        Initialize partial column norms. The first n elements of
+*        work store the exact column norms.
+*
+         DO 20 I = ITEMP + 1, N
+            RWORK( I ) = DZNRM2( M-ITEMP, A( ITEMP+1, I ), 1 )
+            RWORK( N+I ) = RWORK( I )
+   20    CONTINUE
+*
+*        Compute factorization
+*
+         DO 40 I = ITEMP + 1, MN
+*
+*           Determine ith pivot column and swap if necessary
+*
+            PVT = ( I-1 ) + IDAMAX( N-I+1, RWORK( I ), 1 )
+*
+            IF( PVT.NE.I ) THEN
+               CALL ZSWAP( M, A( 1, PVT ), 1, A( 1, I ), 1 )
+               ITEMP = JPVT( PVT )
+               JPVT( PVT ) = JPVT( I )
+               JPVT( I ) = ITEMP
+               RWORK( PVT ) = RWORK( I )
+               RWORK( N+PVT ) = RWORK( N+I )
+            END IF
+*
+*           Generate elementary reflector H(i)
+*
+            AII = A( I, I )
+            CALL ZLARFG( M-I+1, AII, A( MIN( I+1, M ), I ), 1,
+     $                   TAU( I ) )
+            A( I, I ) = AII
+*
+            IF( I.LT.N ) THEN
+*
+*              Apply H(i) to A(i:m,i+1:n) from the left
+*
+               AII = A( I, I )
+               A( I, I ) = DCMPLX( ONE )
+               CALL ZLARF( 'Left', M-I+1, N-I, A( I, I ), 1,
+     $                     DCONJG( TAU( I ) ), A( I, I+1 ), LDA, WORK )
+               A( I, I ) = AII
+            END IF
+*
+*           Update partial column norms
+*
+            DO 30 J = I + 1, N
+               IF( RWORK( J ).NE.ZERO ) THEN
+*
+*                 NOTE: The following 4 lines follow from the analysis in
+*                 Lapack Working Note 176.
+*                 
+                  TEMP = ABS( A( I, J ) ) / RWORK( J )
+                  TEMP = MAX( ZERO, ( ONE+TEMP )*( ONE-TEMP ) )
+                  TEMP2 = TEMP*( RWORK( J ) / RWORK( N+J ) )**2
+                  IF( TEMP2 .LE. TOL3Z ) THEN 
+                     IF( M-I.GT.0 ) THEN
+                        RWORK( J ) = DZNRM2( M-I, A( I+1, J ), 1 )
+                        RWORK( N+J ) = RWORK( J )
+                     ELSE
+                        RWORK( J ) = ZERO
+                        RWORK( N+J ) = ZERO
+                     END IF
+                  ELSE
+                     RWORK( J ) = RWORK( J )*SQRT( TEMP )
+                  END IF
+               END IF
+   30       CONTINUE
+*
+   40    CONTINUE
+      END IF
+      RETURN
+*
+*     End of ZGEQPF
+*
+      END