Mercurial > hg > octave-lyh

comparison libcruft/lapack/dlatrz.f @ 6936:e92bc778c634

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2007-09-30 21:39:05 by dbateman]
author dbateman
date Sun, 30 Sep 2007 21:41:04 +0000
parents
children
comparison
equal deleted inserted replaced
6935:5cd272497aae 6936:e92bc778c634
1 SUBROUTINE DLATRZ( M, N, L, A, LDA, TAU, WORK )
2 *
3 * -- LAPACK routine (version 3.1) --
4 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
5 * November 2006
6 *
7 * .. Scalar Arguments ..
8 INTEGER L, LDA, M, N
9 * ..
10 * .. Array Arguments ..
11 DOUBLE PRECISION A( LDA, * ), TAU( * ), WORK( * )
12 * ..
13 *
14 * Purpose
15 * =======
16 *
17 * DLATRZ factors the M-by-(M+L) real upper trapezoidal matrix
18 * [ A1 A2 ] = [ A(1:M,1:M) A(1:M,N-L+1:N) ] as ( R 0 ) * Z, by means
19 * of orthogonal transformations. Z is an (M+L)-by-(M+L) orthogonal
20 * matrix and, R and A1 are M-by-M upper triangular matrices.
21 *
22 * Arguments
23 * =========
24 *
25 * M (input) INTEGER
26 * The number of rows of the matrix A. M >= 0.
27 *
28 * N (input) INTEGER
29 * The number of columns of the matrix A. N >= 0.
30 *
31 * L (input) INTEGER
32 * The number of columns of the matrix A containing the
33 * meaningful part of the Householder vectors. N-M >= L >= 0.
34 *
35 * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
36 * On entry, the leading M-by-N upper trapezoidal part of the
37 * array A must contain the matrix to be factorized.
38 * On exit, the leading M-by-M upper triangular part of A
39 * contains the upper triangular matrix R, and elements N-L+1 to
40 * N of the first M rows of A, with the array TAU, represent the
41 * orthogonal matrix Z as a product of M elementary reflectors.
42 *
43 * LDA (input) INTEGER
44 * The leading dimension of the array A. LDA >= max(1,M).
45 *
46 * TAU (output) DOUBLE PRECISION array, dimension (M)
47 * The scalar factors of the elementary reflectors.
48 *
49 * WORK (workspace) DOUBLE PRECISION array, dimension (M)
50 *
51 * Further Details
52 * ===============
53 *
54 * Based on contributions by
55 * A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA
56 *
57 * The factorization is obtained by Householder's method. The kth
58 * transformation matrix, Z( k ), which is used to introduce zeros into
59 * the ( m - k + 1 )th row of A, is given in the form
60 *
61 * Z( k ) = ( I 0 ),
62 * ( 0 T( k ) )
63 *
64 * where
65 *
66 * T( k ) = I - tau*u( k )*u( k )', u( k ) = ( 1 ),
67 * ( 0 )
68 * ( z( k ) )
69 *
70 * tau is a scalar and z( k ) is an l element vector. tau and z( k )
71 * are chosen to annihilate the elements of the kth row of A2.
72 *
73 * The scalar tau is returned in the kth element of TAU and the vector
74 * u( k ) in the kth row of A2, such that the elements of z( k ) are
75 * in a( k, l + 1 ), ..., a( k, n ). The elements of R are returned in
76 * the upper triangular part of A1.
77 *
78 * Z is given by
79 *
80 * Z = Z( 1 ) * Z( 2 ) * ... * Z( m ).
81 *
82 * =====================================================================
83 *
84 * .. Parameters ..
85 DOUBLE PRECISION ZERO
86 PARAMETER ( ZERO = 0.0D+0 )
87 * ..
88 * .. Local Scalars ..
89 INTEGER I
90 * ..
91 * .. External Subroutines ..
92 EXTERNAL DLARFG, DLARZ
93 * ..
94 * .. Executable Statements ..
95 *
96 * Test the input arguments
97 *
98 * Quick return if possible
99 *
100 IF( M.EQ.0 ) THEN
101 RETURN
102 ELSE IF( M.EQ.N ) THEN
103 DO 10 I = 1, N
104 TAU( I ) = ZERO
105 10 CONTINUE
106 RETURN
107 END IF
108 *
109 DO 20 I = M, 1, -1
110 *
111 * Generate elementary reflector H(i) to annihilate
112 * [ A(i,i) A(i,n-l+1:n) ]
113 *
114 CALL DLARFG( L+1, A( I, I ), A( I, N-L+1 ), LDA, TAU( I ) )
115 *
116 * Apply H(i) to A(1:i-1,i:n) from the right
117 *
118 CALL DLARZ( 'Right', I-1, N-I+1, L, A( I, N-L+1 ), LDA,
119 $ TAU( I ), A( 1, I ), LDA, WORK )
120 *
121 20 CONTINUE
122 *
123 RETURN
124 *
125 * End of DLATRZ
126 *
127 END