Statistiques
| Révision :

root / src / lapack / double / dorm2r.f @ 1

Historique | Voir | Annoter | Télécharger (5,17 ko)

1 1 equemene
      SUBROUTINE DORM2R( SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC,
2 1 equemene
     $                   WORK, INFO )
3 1 equemene
*
4 1 equemene
*  -- LAPACK routine (version 3.2) --
5 1 equemene
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
6 1 equemene
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7 1 equemene
*     November 2006
8 1 equemene
*
9 1 equemene
*     .. Scalar Arguments ..
10 1 equemene
      CHARACTER          SIDE, TRANS
11 1 equemene
      INTEGER            INFO, K, LDA, LDC, M, N
12 1 equemene
*     ..
13 1 equemene
*     .. Array Arguments ..
14 1 equemene
      DOUBLE PRECISION   A( LDA, * ), C( LDC, * ), TAU( * ), WORK( * )
15 1 equemene
*     ..
16 1 equemene
*
17 1 equemene
*  Purpose
18 1 equemene
*  =======
19 1 equemene
*
20 1 equemene
*  DORM2R overwrites the general real m by n matrix C with
21 1 equemene
*
22 1 equemene
*        Q * C  if SIDE = 'L' and TRANS = 'N', or
23 1 equemene
*
24 1 equemene
*        Q'* C  if SIDE = 'L' and TRANS = 'T', or
25 1 equemene
*
26 1 equemene
*        C * Q  if SIDE = 'R' and TRANS = 'N', or
27 1 equemene
*
28 1 equemene
*        C * Q' if SIDE = 'R' and TRANS = 'T',
29 1 equemene
*
30 1 equemene
*  where Q is a real orthogonal matrix defined as the product of k
31 1 equemene
*  elementary reflectors
32 1 equemene
*
33 1 equemene
*        Q = H(1) H(2) . . . H(k)
34 1 equemene
*
35 1 equemene
*  as returned by DGEQRF. Q is of order m if SIDE = 'L' and of order n
36 1 equemene
*  if SIDE = 'R'.
37 1 equemene
*
38 1 equemene
*  Arguments
39 1 equemene
*  =========
40 1 equemene
*
41 1 equemene
*  SIDE    (input) CHARACTER*1
42 1 equemene
*          = 'L': apply Q or Q' from the Left
43 1 equemene
*          = 'R': apply Q or Q' from the Right
44 1 equemene
*
45 1 equemene
*  TRANS   (input) CHARACTER*1
46 1 equemene
*          = 'N': apply Q  (No transpose)
47 1 equemene
*          = 'T': apply Q' (Transpose)
48 1 equemene
*
49 1 equemene
*  M       (input) INTEGER
50 1 equemene
*          The number of rows of the matrix C. M >= 0.
51 1 equemene
*
52 1 equemene
*  N       (input) INTEGER
53 1 equemene
*          The number of columns of the matrix C. N >= 0.
54 1 equemene
*
55 1 equemene
*  K       (input) INTEGER
56 1 equemene
*          The number of elementary reflectors whose product defines
57 1 equemene
*          the matrix Q.
58 1 equemene
*          If SIDE = 'L', M >= K >= 0;
59 1 equemene
*          if SIDE = 'R', N >= K >= 0.
60 1 equemene
*
61 1 equemene
*  A       (input) DOUBLE PRECISION array, dimension (LDA,K)
62 1 equemene
*          The i-th column must contain the vector which defines the
63 1 equemene
*          elementary reflector H(i), for i = 1,2,...,k, as returned by
64 1 equemene
*          DGEQRF in the first k columns of its array argument A.
65 1 equemene
*          A is modified by the routine but restored on exit.
66 1 equemene
*
67 1 equemene
*  LDA     (input) INTEGER
68 1 equemene
*          The leading dimension of the array A.
69 1 equemene
*          If SIDE = 'L', LDA >= max(1,M);
70 1 equemene
*          if SIDE = 'R', LDA >= max(1,N).
71 1 equemene
*
72 1 equemene
*  TAU     (input) DOUBLE PRECISION array, dimension (K)
73 1 equemene
*          TAU(i) must contain the scalar factor of the elementary
74 1 equemene
*          reflector H(i), as returned by DGEQRF.
75 1 equemene
*
76 1 equemene
*  C       (input/output) DOUBLE PRECISION array, dimension (LDC,N)
77 1 equemene
*          On entry, the m by n matrix C.
78 1 equemene
*          On exit, C is overwritten by Q*C or Q'*C or C*Q' or C*Q.
79 1 equemene
*
80 1 equemene
*  LDC     (input) INTEGER
81 1 equemene
*          The leading dimension of the array C. LDC >= max(1,M).
82 1 equemene
*
83 1 equemene
*  WORK    (workspace) DOUBLE PRECISION array, dimension
84 1 equemene
*                                   (N) if SIDE = 'L',
85 1 equemene
*                                   (M) if SIDE = 'R'
86 1 equemene
*
87 1 equemene
*  INFO    (output) INTEGER
88 1 equemene
*          = 0: successful exit
89 1 equemene
*          < 0: if INFO = -i, the i-th argument had an illegal value
90 1 equemene
*
91 1 equemene
*  =====================================================================
92 1 equemene
*
93 1 equemene
*     .. Parameters ..
94 1 equemene
      DOUBLE PRECISION   ONE
95 1 equemene
      PARAMETER          ( ONE = 1.0D+0 )
96 1 equemene
*     ..
97 1 equemene
*     .. Local Scalars ..
98 1 equemene
      LOGICAL            LEFT, NOTRAN
99 1 equemene
      INTEGER            I, I1, I2, I3, IC, JC, MI, NI, NQ
100 1 equemene
      DOUBLE PRECISION   AII
101 1 equemene
*     ..
102 1 equemene
*     .. External Functions ..
103 1 equemene
      LOGICAL            LSAME
104 1 equemene
      EXTERNAL           LSAME
105 1 equemene
*     ..
106 1 equemene
*     .. External Subroutines ..
107 1 equemene
      EXTERNAL           DLARF, XERBLA
108 1 equemene
*     ..
109 1 equemene
*     .. Intrinsic Functions ..
110 1 equemene
      INTRINSIC          MAX
111 1 equemene
*     ..
112 1 equemene
*     .. Executable Statements ..
113 1 equemene
*
114 1 equemene
*     Test the input arguments
115 1 equemene
*
116 1 equemene
      INFO = 0
117 1 equemene
      LEFT = LSAME( SIDE, 'L' )
118 1 equemene
      NOTRAN = LSAME( TRANS, 'N' )
119 1 equemene
*
120 1 equemene
*     NQ is the order of Q
121 1 equemene
*
122 1 equemene
      IF( LEFT ) THEN
123 1 equemene
         NQ = M
124 1 equemene
      ELSE
125 1 equemene
         NQ = N
126 1 equemene
      END IF
127 1 equemene
      IF( .NOT.LEFT .AND. .NOT.LSAME( SIDE, 'R' ) ) THEN
128 1 equemene
         INFO = -1
129 1 equemene
      ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) ) THEN
130 1 equemene
         INFO = -2
131 1 equemene
      ELSE IF( M.LT.0 ) THEN
132 1 equemene
         INFO = -3
133 1 equemene
      ELSE IF( N.LT.0 ) THEN
134 1 equemene
         INFO = -4
135 1 equemene
      ELSE IF( K.LT.0 .OR. K.GT.NQ ) THEN
136 1 equemene
         INFO = -5
137 1 equemene
      ELSE IF( LDA.LT.MAX( 1, NQ ) ) THEN
138 1 equemene
         INFO = -7
139 1 equemene
      ELSE IF( LDC.LT.MAX( 1, M ) ) THEN
140 1 equemene
         INFO = -10
141 1 equemene
      END IF
142 1 equemene
      IF( INFO.NE.0 ) THEN
143 1 equemene
         CALL XERBLA( 'DORM2R', -INFO )
144 1 equemene
         RETURN
145 1 equemene
      END IF
146 1 equemene
*
147 1 equemene
*     Quick return if possible
148 1 equemene
*
149 1 equemene
      IF( M.EQ.0 .OR. N.EQ.0 .OR. K.EQ.0 )
150 1 equemene
     $   RETURN
151 1 equemene
*
152 1 equemene
      IF( ( LEFT .AND. .NOT.NOTRAN ) .OR. ( .NOT.LEFT .AND. NOTRAN ) )
153 1 equemene
     $     THEN
154 1 equemene
         I1 = 1
155 1 equemene
         I2 = K
156 1 equemene
         I3 = 1
157 1 equemene
      ELSE
158 1 equemene
         I1 = K
159 1 equemene
         I2 = 1
160 1 equemene
         I3 = -1
161 1 equemene
      END IF
162 1 equemene
*
163 1 equemene
      IF( LEFT ) THEN
164 1 equemene
         NI = N
165 1 equemene
         JC = 1
166 1 equemene
      ELSE
167 1 equemene
         MI = M
168 1 equemene
         IC = 1
169 1 equemene
      END IF
170 1 equemene
*
171 1 equemene
      DO 10 I = I1, I2, I3
172 1 equemene
         IF( LEFT ) THEN
173 1 equemene
*
174 1 equemene
*           H(i) is applied to C(i:m,1:n)
175 1 equemene
*
176 1 equemene
            MI = M - I + 1
177 1 equemene
            IC = I
178 1 equemene
         ELSE
179 1 equemene
*
180 1 equemene
*           H(i) is applied to C(1:m,i:n)
181 1 equemene
*
182 1 equemene
            NI = N - I + 1
183 1 equemene
            JC = I
184 1 equemene
         END IF
185 1 equemene
*
186 1 equemene
*        Apply H(i)
187 1 equemene
*
188 1 equemene
         AII = A( I, I )
189 1 equemene
         A( I, I ) = ONE
190 1 equemene
         CALL DLARF( SIDE, MI, NI, A( I, I ), 1, TAU( I ), C( IC, JC ),
191 1 equemene
     $               LDC, WORK )
192 1 equemene
         A( I, I ) = AII
193 1 equemene
   10 CONTINUE
194 1 equemene
      RETURN
195 1 equemene
*
196 1 equemene
*     End of DORM2R
197 1 equemene
*
198 1 equemene
      END