summaryrefslogtreecommitdiff
path: root/src/fortran/lapack/zunmrz.f
diff options
context:
space:
mode:
Diffstat (limited to 'src/fortran/lapack/zunmrz.f')
-rw-r--r--src/fortran/lapack/zunmrz.f296
1 files changed, 296 insertions, 0 deletions
diff --git a/src/fortran/lapack/zunmrz.f b/src/fortran/lapack/zunmrz.f
new file mode 100644
index 0000000..c763705
--- /dev/null
+++ b/src/fortran/lapack/zunmrz.f
@@ -0,0 +1,296 @@
+ SUBROUTINE ZUNMRZ( SIDE, TRANS, M, N, K, L, A, LDA, TAU, C, LDC,
+ $ WORK, LWORK, INFO )
+*
+* -- LAPACK routine (version 3.1) --
+* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
+* November 2006
+*
+* .. Scalar Arguments ..
+ CHARACTER SIDE, TRANS
+ INTEGER INFO, K, L, LDA, LDC, LWORK, M, N
+* ..
+* .. Array Arguments ..
+ COMPLEX*16 A( LDA, * ), C( LDC, * ), TAU( * ), WORK( * )
+* ..
+*
+* Purpose
+* =======
+*
+* ZUNMRZ overwrites the general complex M-by-N matrix C with
+*
+* SIDE = 'L' SIDE = 'R'
+* TRANS = 'N': Q * C C * Q
+* TRANS = 'C': Q**H * C C * Q**H
+*
+* where Q is a complex unitary matrix defined as the product of k
+* elementary reflectors
+*
+* Q = H(1) H(2) . . . H(k)
+*
+* as returned by ZTZRZF. Q is of order M if SIDE = 'L' and of order N
+* if SIDE = 'R'.
+*
+* Arguments
+* =========
+*
+* SIDE (input) CHARACTER*1
+* = 'L': apply Q or Q**H from the Left;
+* = 'R': apply Q or Q**H from the Right.
+*
+* TRANS (input) CHARACTER*1
+* = 'N': No transpose, apply Q;
+* = 'C': Conjugate transpose, apply Q**H.
+*
+* M (input) INTEGER
+* The number of rows of the matrix C. M >= 0.
+*
+* N (input) INTEGER
+* The number of columns of the matrix C. N >= 0.
+*
+* K (input) INTEGER
+* The number of elementary reflectors whose product defines
+* the matrix Q.
+* If SIDE = 'L', M >= K >= 0;
+* if SIDE = 'R', N >= K >= 0.
+*
+* L (input) INTEGER
+* The number of columns of the matrix A containing
+* the meaningful part of the Householder reflectors.
+* If SIDE = 'L', M >= L >= 0, if SIDE = 'R', N >= L >= 0.
+*
+* A (input) COMPLEX*16 array, dimension
+* (LDA,M) if SIDE = 'L',
+* (LDA,N) if SIDE = 'R'
+* The i-th row must contain the vector which defines the
+* elementary reflector H(i), for i = 1,2,...,k, as returned by
+* ZTZRZF in the last k rows of its array argument A.
+* A is modified by the routine but restored on exit.
+*
+* LDA (input) INTEGER
+* The leading dimension of the array A. LDA >= max(1,K).
+*
+* TAU (input) COMPLEX*16 array, dimension (K)
+* TAU(i) must contain the scalar factor of the elementary
+* reflector H(i), as returned by ZTZRZF.
+*
+* C (input/output) COMPLEX*16 array, dimension (LDC,N)
+* On entry, the M-by-N matrix C.
+* On exit, C is overwritten by Q*C or Q**H*C or C*Q**H or C*Q.
+*
+* LDC (input) INTEGER
+* The leading dimension of the array C. LDC >= max(1,M).
+*
+* WORK (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
+* On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
+*
+* LWORK (input) INTEGER
+* The dimension of the array WORK.
+* If SIDE = 'L', LWORK >= max(1,N);
+* if SIDE = 'R', LWORK >= max(1,M).
+* For optimum performance LWORK >= N*NB if SIDE = 'L', and
+* LWORK >= M*NB if SIDE = 'R', where NB is the optimal
+* blocksize.
+*
+* If LWORK = -1, then a workspace query is assumed; the routine
+* only calculates the optimal size of the WORK array, returns
+* this value as the first entry of the WORK array, and no error
+* message related to LWORK is issued by XERBLA.
+*
+* INFO (output) INTEGER
+* = 0: successful exit
+* < 0: if INFO = -i, the i-th argument had an illegal value
+*
+* Further Details
+* ===============
+*
+* Based on contributions by
+* A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA
+*
+* =====================================================================
+*
+* .. Parameters ..
+ INTEGER NBMAX, LDT
+ PARAMETER ( NBMAX = 64, LDT = NBMAX+1 )
+* ..
+* .. Local Scalars ..
+ LOGICAL LEFT, LQUERY, NOTRAN
+ CHARACTER TRANST
+ INTEGER I, I1, I2, I3, IB, IC, IINFO, IWS, JA, JC,
+ $ LDWORK, LWKOPT, MI, NB, NBMIN, NI, NQ, NW
+* ..
+* .. Local Arrays ..
+ COMPLEX*16 T( LDT, NBMAX )
+* ..
+* .. External Functions ..
+ LOGICAL LSAME
+ INTEGER ILAENV
+ EXTERNAL LSAME, ILAENV
+* ..
+* .. External Subroutines ..
+ EXTERNAL XERBLA, ZLARZB, ZLARZT, ZUNMR3
+* ..
+* .. Intrinsic Functions ..
+ INTRINSIC MAX, MIN
+* ..
+* .. Executable Statements ..
+*
+* Test the input arguments
+*
+ INFO = 0
+ LEFT = LSAME( SIDE, 'L' )
+ NOTRAN = LSAME( TRANS, 'N' )
+ LQUERY = ( LWORK.EQ.-1 )
+*
+* NQ is the order of Q and NW is the minimum dimension of WORK
+*
+ IF( LEFT ) THEN
+ NQ = M
+ NW = MAX( 1, N )
+ ELSE
+ NQ = N
+ NW = MAX( 1, M )
+ END IF
+ IF( .NOT.LEFT .AND. .NOT.LSAME( SIDE, 'R' ) ) THEN
+ INFO = -1
+ ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'C' ) ) THEN
+ INFO = -2
+ ELSE IF( M.LT.0 ) THEN
+ INFO = -3
+ ELSE IF( N.LT.0 ) THEN
+ INFO = -4
+ ELSE IF( K.LT.0 .OR. K.GT.NQ ) THEN
+ INFO = -5
+ ELSE IF( L.LT.0 .OR. ( LEFT .AND. ( L.GT.M ) ) .OR.
+ $ ( .NOT.LEFT .AND. ( L.GT.N ) ) ) THEN
+ INFO = -6
+ ELSE IF( LDA.LT.MAX( 1, K ) ) THEN
+ INFO = -8
+ ELSE IF( LDC.LT.MAX( 1, M ) ) THEN
+ INFO = -11
+ END IF
+*
+ IF( INFO.EQ.0 ) THEN
+ IF( M.EQ.0 .OR. N.EQ.0 ) THEN
+ LWKOPT = 1
+*
+* Determine the block size. NB may be at most NBMAX, where
+* NBMAX is used to define the local array T.
+*
+ NB = MIN( NBMAX, ILAENV( 1, 'ZUNMRQ', SIDE // TRANS, M, N,
+ $ K, -1 ) )
+ LWKOPT = NW*NB
+ END IF
+ WORK( 1 ) = LWKOPT
+*
+ IF( LWORK.LT.MAX( 1, NW ) .AND. .NOT.LQUERY ) THEN
+ INFO = -13
+ END IF
+ END IF
+*
+ IF( INFO.NE.0 ) THEN
+ CALL XERBLA( 'ZUNMRZ', -INFO )
+ RETURN
+ ELSE IF( LQUERY ) THEN
+ RETURN
+ END IF
+*
+* Quick return if possible
+*
+ IF( M.EQ.0 .OR. N.EQ.0 ) THEN
+ RETURN
+ END IF
+*
+* Determine the block size. NB may be at most NBMAX, where NBMAX
+* is used to define the local array T.
+*
+ NB = MIN( NBMAX, ILAENV( 1, 'ZUNMRQ', SIDE // TRANS, M, N, K,
+ $ -1 ) )
+ NBMIN = 2
+ LDWORK = NW
+ IF( NB.GT.1 .AND. NB.LT.K ) THEN
+ IWS = NW*NB
+ IF( LWORK.LT.IWS ) THEN
+ NB = LWORK / LDWORK
+ NBMIN = MAX( 2, ILAENV( 2, 'ZUNMRQ', SIDE // TRANS, M, N, K,
+ $ -1 ) )
+ END IF
+ ELSE
+ IWS = NW
+ END IF
+*
+ IF( NB.LT.NBMIN .OR. NB.GE.K ) THEN
+*
+* Use unblocked code
+*
+ CALL ZUNMR3( SIDE, TRANS, M, N, K, L, A, LDA, TAU, C, LDC,
+ $ WORK, IINFO )
+ ELSE
+*
+* Use blocked code
+*
+ IF( ( LEFT .AND. .NOT.NOTRAN ) .OR.
+ $ ( .NOT.LEFT .AND. NOTRAN ) ) THEN
+ I1 = 1
+ I2 = K
+ I3 = NB
+ ELSE
+ I1 = ( ( K-1 ) / NB )*NB + 1
+ I2 = 1
+ I3 = -NB
+ END IF
+*
+ IF( LEFT ) THEN
+ NI = N
+ JC = 1
+ JA = M - L + 1
+ ELSE
+ MI = M
+ IC = 1
+ JA = N - L + 1
+ END IF
+*
+ IF( NOTRAN ) THEN
+ TRANST = 'C'
+ ELSE
+ TRANST = 'N'
+ END IF
+*
+ DO 10 I = I1, I2, I3
+ IB = MIN( NB, K-I+1 )
+*
+* Form the triangular factor of the block reflector
+* H = H(i+ib-1) . . . H(i+1) H(i)
+*
+ CALL ZLARZT( 'Backward', 'Rowwise', L, IB, A( I, JA ), LDA,
+ $ TAU( I ), T, LDT )
+*
+ IF( LEFT ) THEN
+*
+* H or H' is applied to C(i:m,1:n)
+*
+ MI = M - I + 1
+ IC = I
+ ELSE
+*
+* H or H' is applied to C(1:m,i:n)
+*
+ NI = N - I + 1
+ JC = I
+ END IF
+*
+* Apply H or H'
+*
+ CALL ZLARZB( SIDE, TRANST, 'Backward', 'Rowwise', MI, NI,
+ $ IB, L, A( I, JA ), LDA, T, LDT, C( IC, JC ),
+ $ LDC, WORK, LDWORK )
+ 10 CONTINUE
+*
+ END IF
+*
+ WORK( 1 ) = LWKOPT
+*
+ RETURN
+*
+* End of ZUNMRZ
+*
+ END