← Back to branch summary

~ubuntu-branches/ubuntu/hoary/scilab/hoary

« back to all changes in this revision

Viewing changes to routines/lapack/zlatrz.f

  • Committer: Bazaar Package Importer
  • Author(s): Torsten Werner
  • Date: 2005-01-09 22:58:21 UTC
  • mfrom: (1.1.1 upstream)
  • Revision ID: james.westby@ubuntu.com-20050109225821-473xr8vhgugxxx5j
Tags: 3.0-12
http://bugs.debian.org/255869
changed configure.in to build scilab's own malloc.o, closes: #255869

Show diffs side-by-side

added added

removed removed

Lines of Context:
routines/lapack/zlatrz.f
 
1
      SUBROUTINE ZLATRZ( M, N, L, A, LDA, TAU, WORK )
 
2
*
 
3
*  -- LAPACK routine (version 3.0) --
 
4
*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
 
5
*     Courant Institute, Argonne National Lab, and Rice University
 
6
*     June 30, 1999
 
7
*
 
8
*     .. Scalar Arguments ..
 
9
      INTEGER            L, LDA, M, N
 
10
*     ..
 
11
*     .. Array Arguments ..
 
12
      COMPLEX*16         A( LDA, * ), TAU( * ), WORK( * )
 
13
*     ..
 
14
*
 
15
*  Purpose
 
16
*  =======
 
17
*
 
18
*  ZLATRZ factors the M-by-(M+L) complex upper trapezoidal matrix
 
19
*  [ A1 A2 ] = [ A(1:M,1:M) A(1:M,N-L+1:N) ] as ( R  0 ) * Z by means
 
20
*  of unitary transformations, where  Z is an (M+L)-by-(M+L) unitary
 
21
*  matrix and, R and A1 are M-by-M upper triangular matrices.
 
22
*
 
23
*  Arguments
 
24
*  =========
 
25
*
 
26
*  M       (input) INTEGER
 
27
*          The number of rows of the matrix A.  M >= 0.
 
28
*
 
29
*  N       (input) INTEGER
 
30
*          The number of columns of the matrix A.  N >= 0.
 
31
*
 
32
*  L       (input) INTEGER
 
33
*          The number of columns of the matrix A containing the
 
34
*          meaningful part of the Householder vectors. N-M >= L >= 0.
 
35
*
 
36
*  A       (input/output) COMPLEX*16 array, dimension (LDA,N)
 
37
*          On entry, the leading M-by-N upper trapezoidal part of the
 
38
*          array A must contain the matrix to be factorized.
 
39
*          On exit, the leading M-by-M upper triangular part of A
 
40
*          contains the upper triangular matrix R, and elements N-L+1 to
 
41
*          N of the first M rows of A, with the array TAU, represent the
 
42
*          unitary matrix Z as a product of M elementary reflectors.
 
43
*
 
44
*  LDA     (input) INTEGER
 
45
*          The leading dimension of the array A.  LDA >= max(1,M).
 
46
*
 
47
*  TAU     (output) COMPLEX*16 array, dimension (M)
 
48
*          The scalar factors of the elementary reflectors.
 
49
*
 
50
*  WORK    (workspace) COMPLEX*16 array, dimension (M)
 
51
*
 
52
*  Further Details
 
53
*  ===============
 
54
*
 
55
*  Based on contributions by
 
56
*    A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA
 
57
*
 
58
*  The factorization is obtained by Householder's method.  The kth
 
59
*  transformation matrix, Z( k ), which is used to introduce zeros into
 
60
*  the ( m - k + 1 )th row of A, is given in the form
 
61
*
 
62
*     Z( k ) = ( I     0   ),
 
63
*              ( 0  T( k ) )
 
64
*
 
65
*  where
 
66
*
 
67
*     T( k ) = I - tau*u( k )*u( k )',   u( k ) = (   1    ),
 
68
*                                                 (   0    )
 
69
*                                                 ( z( k ) )
 
70
*
 
71
*  tau is a scalar and z( k ) is an l element vector. tau and z( k )
 
72
*  are chosen to annihilate the elements of the kth row of A2.
 
73
*
 
74
*  The scalar tau is returned in the kth element of TAU and the vector
 
75
*  u( k ) in the kth row of A2, such that the elements of z( k ) are
 
76
*  in  a( k, l + 1 ), ..., a( k, n ). The elements of R are returned in
 
77
*  the upper triangular part of A1.
 
78
*
 
79
*  Z is given by
 
80
*
 
81
*     Z =  Z( 1 ) * Z( 2 ) * ... * Z( m ).
 
82
*
 
83
*  =====================================================================
 
84
*
 
85
*     .. Parameters ..
 
86
      COMPLEX*16         ZERO
 
87
      PARAMETER          ( ZERO = ( 0.0D+0, 0.0D+0 ) )
 
88
*     ..
 
89
*     .. Local Scalars ..
 
90
      INTEGER            I
 
91
      COMPLEX*16         ALPHA
 
92
*     ..
 
93
*     .. External Subroutines ..
 
94
      EXTERNAL           ZLACGV, ZLARFG, ZLARZ
 
95
*     ..
 
96
*     .. Intrinsic Functions ..
 
97
      INTRINSIC          DCONJG
 
98
*     ..
 
99
*     .. Executable Statements ..
 
100
*
 
101
*     Quick return if possible
 
102
*
 
103
      IF( M.EQ.0 ) THEN
 
104
         RETURN
 
105
      ELSE IF( M.EQ.N ) THEN
 
106
         DO 10 I = 1, N
 
107
            TAU( I ) = ZERO
 
108
   10    CONTINUE
 
109
         RETURN
 
110
      END IF
 
111
*
 
112
      DO 20 I = M, 1, -1
 
113
*
 
114
*        Generate elementary reflector H(i) to annihilate
 
115
*        [ A(i,i) A(i,n-l+1:n) ]
 
116
*
 
117
         CALL ZLACGV( L, A( I, N-L+1 ), LDA )
 
118
         ALPHA = DCONJG( A( I, I ) )
 
119
         CALL ZLARFG( L+1, ALPHA, A( I, N-L+1 ), LDA, TAU( I ) )
 
120
         TAU( I ) = DCONJG( TAU( I ) )
 
121
*
 
122
*        Apply H(i) to A(1:i-1,i:n) from the right
 
123
*
 
124
         CALL ZLARZ( 'Right', I-1, N-I+1, L, A( I, N-L+1 ), LDA,
 
125
     $               DCONJG( TAU( I ) ), A( 1, I ), LDA, WORK )
 
126
         A( I, I ) = DCONJG( ALPHA )
 
127
*
 
128
   20 CONTINUE
 
129
*
 
130
      RETURN
 
131
*
 
132
*     End of ZLATRZ
 
133
*
 
134
      END