Statistiques
| Révision :

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

Historique | Voir | Annoter | Télécharger (12,16 ko)

1 1 equemene
      SUBROUTINE DLALSA( ICOMPQ, SMLSIZ, N, NRHS, B, LDB, BX, LDBX, U,
2 1 equemene
     $                   LDU, VT, K, DIFL, DIFR, Z, POLES, GIVPTR,
3 1 equemene
     $                   GIVCOL, LDGCOL, PERM, GIVNUM, C, S, WORK,
4 1 equemene
     $                   IWORK, INFO )
5 1 equemene
*
6 1 equemene
*  -- LAPACK routine (version 3.2) --
7 1 equemene
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
8 1 equemene
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
9 1 equemene
*     November 2006
10 1 equemene
*
11 1 equemene
*     .. Scalar Arguments ..
12 1 equemene
      INTEGER            ICOMPQ, INFO, LDB, LDBX, LDGCOL, LDU, N, NRHS,
13 1 equemene
     $                   SMLSIZ
14 1 equemene
*     ..
15 1 equemene
*     .. Array Arguments ..
16 1 equemene
      INTEGER            GIVCOL( LDGCOL, * ), GIVPTR( * ), IWORK( * ),
17 1 equemene
     $                   K( * ), PERM( LDGCOL, * )
18 1 equemene
      DOUBLE PRECISION   B( LDB, * ), BX( LDBX, * ), C( * ),
19 1 equemene
     $                   DIFL( LDU, * ), DIFR( LDU, * ),
20 1 equemene
     $                   GIVNUM( LDU, * ), POLES( LDU, * ), S( * ),
21 1 equemene
     $                   U( LDU, * ), VT( LDU, * ), WORK( * ),
22 1 equemene
     $                   Z( LDU, * )
23 1 equemene
*     ..
24 1 equemene
*
25 1 equemene
*  Purpose
26 1 equemene
*  =======
27 1 equemene
*
28 1 equemene
*  DLALSA is an itermediate step in solving the least squares problem
29 1 equemene
*  by computing the SVD of the coefficient matrix in compact form (The
30 1 equemene
*  singular vectors are computed as products of simple orthorgonal
31 1 equemene
*  matrices.).
32 1 equemene
*
33 1 equemene
*  If ICOMPQ = 0, DLALSA applies the inverse of the left singular vector
34 1 equemene
*  matrix of an upper bidiagonal matrix to the right hand side; and if
35 1 equemene
*  ICOMPQ = 1, DLALSA applies the right singular vector matrix to the
36 1 equemene
*  right hand side. The singular vector matrices were generated in
37 1 equemene
*  compact form by DLALSA.
38 1 equemene
*
39 1 equemene
*  Arguments
40 1 equemene
*  =========
41 1 equemene
*
42 1 equemene
*
43 1 equemene
*  ICOMPQ (input) INTEGER
44 1 equemene
*         Specifies whether the left or the right singular vector
45 1 equemene
*         matrix is involved.
46 1 equemene
*         = 0: Left singular vector matrix
47 1 equemene
*         = 1: Right singular vector matrix
48 1 equemene
*
49 1 equemene
*  SMLSIZ (input) INTEGER
50 1 equemene
*         The maximum size of the subproblems at the bottom of the
51 1 equemene
*         computation tree.
52 1 equemene
*
53 1 equemene
*  N      (input) INTEGER
54 1 equemene
*         The row and column dimensions of the upper bidiagonal matrix.
55 1 equemene
*
56 1 equemene
*  NRHS   (input) INTEGER
57 1 equemene
*         The number of columns of B and BX. NRHS must be at least 1.
58 1 equemene
*
59 1 equemene
*  B      (input/output) DOUBLE PRECISION array, dimension ( LDB, NRHS )
60 1 equemene
*         On input, B contains the right hand sides of the least
61 1 equemene
*         squares problem in rows 1 through M.
62 1 equemene
*         On output, B contains the solution X in rows 1 through N.
63 1 equemene
*
64 1 equemene
*  LDB    (input) INTEGER
65 1 equemene
*         The leading dimension of B in the calling subprogram.
66 1 equemene
*         LDB must be at least max(1,MAX( M, N ) ).
67 1 equemene
*
68 1 equemene
*  BX     (output) DOUBLE PRECISION array, dimension ( LDBX, NRHS )
69 1 equemene
*         On exit, the result of applying the left or right singular
70 1 equemene
*         vector matrix to B.
71 1 equemene
*
72 1 equemene
*  LDBX   (input) INTEGER
73 1 equemene
*         The leading dimension of BX.
74 1 equemene
*
75 1 equemene
*  U      (input) DOUBLE PRECISION array, dimension ( LDU, SMLSIZ ).
76 1 equemene
*         On entry, U contains the left singular vector matrices of all
77 1 equemene
*         subproblems at the bottom level.
78 1 equemene
*
79 1 equemene
*  LDU    (input) INTEGER, LDU = > N.
80 1 equemene
*         The leading dimension of arrays U, VT, DIFL, DIFR,
81 1 equemene
*         POLES, GIVNUM, and Z.
82 1 equemene
*
83 1 equemene
*  VT     (input) DOUBLE PRECISION array, dimension ( LDU, SMLSIZ+1 ).
84 1 equemene
*         On entry, VT' contains the right singular vector matrices of
85 1 equemene
*         all subproblems at the bottom level.
86 1 equemene
*
87 1 equemene
*  K      (input) INTEGER array, dimension ( N ).
88 1 equemene
*
89 1 equemene
*  DIFL   (input) DOUBLE PRECISION array, dimension ( LDU, NLVL ).
90 1 equemene
*         where NLVL = INT(log_2 (N/(SMLSIZ+1))) + 1.
91 1 equemene
*
92 1 equemene
*  DIFR   (input) DOUBLE PRECISION array, dimension ( LDU, 2 * NLVL ).
93 1 equemene
*         On entry, DIFL(*, I) and DIFR(*, 2 * I -1) record
94 1 equemene
*         distances between singular values on the I-th level and
95 1 equemene
*         singular values on the (I -1)-th level, and DIFR(*, 2 * I)
96 1 equemene
*         record the normalizing factors of the right singular vectors
97 1 equemene
*         matrices of subproblems on I-th level.
98 1 equemene
*
99 1 equemene
*  Z      (input) DOUBLE PRECISION array, dimension ( LDU, NLVL ).
100 1 equemene
*         On entry, Z(1, I) contains the components of the deflation-
101 1 equemene
*         adjusted updating row vector for subproblems on the I-th
102 1 equemene
*         level.
103 1 equemene
*
104 1 equemene
*  POLES  (input) DOUBLE PRECISION array, dimension ( LDU, 2 * NLVL ).
105 1 equemene
*         On entry, POLES(*, 2 * I -1: 2 * I) contains the new and old
106 1 equemene
*         singular values involved in the secular equations on the I-th
107 1 equemene
*         level.
108 1 equemene
*
109 1 equemene
*  GIVPTR (input) INTEGER array, dimension ( N ).
110 1 equemene
*         On entry, GIVPTR( I ) records the number of Givens
111 1 equemene
*         rotations performed on the I-th problem on the computation
112 1 equemene
*         tree.
113 1 equemene
*
114 1 equemene
*  GIVCOL (input) INTEGER array, dimension ( LDGCOL, 2 * NLVL ).
115 1 equemene
*         On entry, for each I, GIVCOL(*, 2 * I - 1: 2 * I) records the
116 1 equemene
*         locations of Givens rotations performed on the I-th level on
117 1 equemene
*         the computation tree.
118 1 equemene
*
119 1 equemene
*  LDGCOL (input) INTEGER, LDGCOL = > N.
120 1 equemene
*         The leading dimension of arrays GIVCOL and PERM.
121 1 equemene
*
122 1 equemene
*  PERM   (input) INTEGER array, dimension ( LDGCOL, NLVL ).
123 1 equemene
*         On entry, PERM(*, I) records permutations done on the I-th
124 1 equemene
*         level of the computation tree.
125 1 equemene
*
126 1 equemene
*  GIVNUM (input) DOUBLE PRECISION array, dimension ( LDU, 2 * NLVL ).
127 1 equemene
*         On entry, GIVNUM(*, 2 *I -1 : 2 * I) records the C- and S-
128 1 equemene
*         values of Givens rotations performed on the I-th level on the
129 1 equemene
*         computation tree.
130 1 equemene
*
131 1 equemene
*  C      (input) DOUBLE PRECISION array, dimension ( N ).
132 1 equemene
*         On entry, if the I-th subproblem is not square,
133 1 equemene
*         C( I ) contains the C-value of a Givens rotation related to
134 1 equemene
*         the right null space of the I-th subproblem.
135 1 equemene
*
136 1 equemene
*  S      (input) DOUBLE PRECISION array, dimension ( N ).
137 1 equemene
*         On entry, if the I-th subproblem is not square,
138 1 equemene
*         S( I ) contains the S-value of a Givens rotation related to
139 1 equemene
*         the right null space of the I-th subproblem.
140 1 equemene
*
141 1 equemene
*  WORK   (workspace) DOUBLE PRECISION array.
142 1 equemene
*         The dimension must be at least N.
143 1 equemene
*
144 1 equemene
*  IWORK  (workspace) INTEGER array.
145 1 equemene
*         The dimension must be at least 3 * N
146 1 equemene
*
147 1 equemene
*  INFO   (output) INTEGER
148 1 equemene
*          = 0:  successful exit.
149 1 equemene
*          < 0:  if INFO = -i, the i-th argument had an illegal value.
150 1 equemene
*
151 1 equemene
*  Further Details
152 1 equemene
*  ===============
153 1 equemene
*
154 1 equemene
*  Based on contributions by
155 1 equemene
*     Ming Gu and Ren-Cang Li, Computer Science Division, University of
156 1 equemene
*       California at Berkeley, USA
157 1 equemene
*     Osni Marques, LBNL/NERSC, USA
158 1 equemene
*
159 1 equemene
*  =====================================================================
160 1 equemene
*
161 1 equemene
*     .. Parameters ..
162 1 equemene
      DOUBLE PRECISION   ZERO, ONE
163 1 equemene
      PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
164 1 equemene
*     ..
165 1 equemene
*     .. Local Scalars ..
166 1 equemene
      INTEGER            I, I1, IC, IM1, INODE, J, LF, LL, LVL, LVL2,
167 1 equemene
     $                   ND, NDB1, NDIML, NDIMR, NL, NLF, NLP1, NLVL,
168 1 equemene
     $                   NR, NRF, NRP1, SQRE
169 1 equemene
*     ..
170 1 equemene
*     .. External Subroutines ..
171 1 equemene
      EXTERNAL           DCOPY, DGEMM, DLALS0, DLASDT, XERBLA
172 1 equemene
*     ..
173 1 equemene
*     .. Executable Statements ..
174 1 equemene
*
175 1 equemene
*     Test the input parameters.
176 1 equemene
*
177 1 equemene
      INFO = 0
178 1 equemene
*
179 1 equemene
      IF( ( ICOMPQ.LT.0 ) .OR. ( ICOMPQ.GT.1 ) ) THEN
180 1 equemene
         INFO = -1
181 1 equemene
      ELSE IF( SMLSIZ.LT.3 ) THEN
182 1 equemene
         INFO = -2
183 1 equemene
      ELSE IF( N.LT.SMLSIZ ) THEN
184 1 equemene
         INFO = -3
185 1 equemene
      ELSE IF( NRHS.LT.1 ) THEN
186 1 equemene
         INFO = -4
187 1 equemene
      ELSE IF( LDB.LT.N ) THEN
188 1 equemene
         INFO = -6
189 1 equemene
      ELSE IF( LDBX.LT.N ) THEN
190 1 equemene
         INFO = -8
191 1 equemene
      ELSE IF( LDU.LT.N ) THEN
192 1 equemene
         INFO = -10
193 1 equemene
      ELSE IF( LDGCOL.LT.N ) THEN
194 1 equemene
         INFO = -19
195 1 equemene
      END IF
196 1 equemene
      IF( INFO.NE.0 ) THEN
197 1 equemene
         CALL XERBLA( 'DLALSA', -INFO )
198 1 equemene
         RETURN
199 1 equemene
      END IF
200 1 equemene
*
201 1 equemene
*     Book-keeping and  setting up the computation tree.
202 1 equemene
*
203 1 equemene
      INODE = 1
204 1 equemene
      NDIML = INODE + N
205 1 equemene
      NDIMR = NDIML + N
206 1 equemene
*
207 1 equemene
      CALL DLASDT( N, NLVL, ND, IWORK( INODE ), IWORK( NDIML ),
208 1 equemene
     $             IWORK( NDIMR ), SMLSIZ )
209 1 equemene
*
210 1 equemene
*     The following code applies back the left singular vector factors.
211 1 equemene
*     For applying back the right singular vector factors, go to 50.
212 1 equemene
*
213 1 equemene
      IF( ICOMPQ.EQ.1 ) THEN
214 1 equemene
         GO TO 50
215 1 equemene
      END IF
216 1 equemene
*
217 1 equemene
*     The nodes on the bottom level of the tree were solved
218 1 equemene
*     by DLASDQ. The corresponding left and right singular vector
219 1 equemene
*     matrices are in explicit form. First apply back the left
220 1 equemene
*     singular vector matrices.
221 1 equemene
*
222 1 equemene
      NDB1 = ( ND+1 ) / 2
223 1 equemene
      DO 10 I = NDB1, ND
224 1 equemene
*
225 1 equemene
*        IC : center row of each node
226 1 equemene
*        NL : number of rows of left  subproblem
227 1 equemene
*        NR : number of rows of right subproblem
228 1 equemene
*        NLF: starting row of the left   subproblem
229 1 equemene
*        NRF: starting row of the right  subproblem
230 1 equemene
*
231 1 equemene
         I1 = I - 1
232 1 equemene
         IC = IWORK( INODE+I1 )
233 1 equemene
         NL = IWORK( NDIML+I1 )
234 1 equemene
         NR = IWORK( NDIMR+I1 )
235 1 equemene
         NLF = IC - NL
236 1 equemene
         NRF = IC + 1
237 1 equemene
         CALL DGEMM( 'T', 'N', NL, NRHS, NL, ONE, U( NLF, 1 ), LDU,
238 1 equemene
     $               B( NLF, 1 ), LDB, ZERO, BX( NLF, 1 ), LDBX )
239 1 equemene
         CALL DGEMM( 'T', 'N', NR, NRHS, NR, ONE, U( NRF, 1 ), LDU,
240 1 equemene
     $               B( NRF, 1 ), LDB, ZERO, BX( NRF, 1 ), LDBX )
241 1 equemene
   10 CONTINUE
242 1 equemene
*
243 1 equemene
*     Next copy the rows of B that correspond to unchanged rows
244 1 equemene
*     in the bidiagonal matrix to BX.
245 1 equemene
*
246 1 equemene
      DO 20 I = 1, ND
247 1 equemene
         IC = IWORK( INODE+I-1 )
248 1 equemene
         CALL DCOPY( NRHS, B( IC, 1 ), LDB, BX( IC, 1 ), LDBX )
249 1 equemene
   20 CONTINUE
250 1 equemene
*
251 1 equemene
*     Finally go through the left singular vector matrices of all
252 1 equemene
*     the other subproblems bottom-up on the tree.
253 1 equemene
*
254 1 equemene
      J = 2**NLVL
255 1 equemene
      SQRE = 0
256 1 equemene
*
257 1 equemene
      DO 40 LVL = NLVL, 1, -1
258 1 equemene
         LVL2 = 2*LVL - 1
259 1 equemene
*
260 1 equemene
*        find the first node LF and last node LL on
261 1 equemene
*        the current level LVL
262 1 equemene
*
263 1 equemene
         IF( LVL.EQ.1 ) THEN
264 1 equemene
            LF = 1
265 1 equemene
            LL = 1
266 1 equemene
         ELSE
267 1 equemene
            LF = 2**( LVL-1 )
268 1 equemene
            LL = 2*LF - 1
269 1 equemene
         END IF
270 1 equemene
         DO 30 I = LF, LL
271 1 equemene
            IM1 = I - 1
272 1 equemene
            IC = IWORK( INODE+IM1 )
273 1 equemene
            NL = IWORK( NDIML+IM1 )
274 1 equemene
            NR = IWORK( NDIMR+IM1 )
275 1 equemene
            NLF = IC - NL
276 1 equemene
            NRF = IC + 1
277 1 equemene
            J = J - 1
278 1 equemene
            CALL DLALS0( ICOMPQ, NL, NR, SQRE, NRHS, BX( NLF, 1 ), LDBX,
279 1 equemene
     $                   B( NLF, 1 ), LDB, PERM( NLF, LVL ),
280 1 equemene
     $                   GIVPTR( J ), GIVCOL( NLF, LVL2 ), LDGCOL,
281 1 equemene
     $                   GIVNUM( NLF, LVL2 ), LDU, POLES( NLF, LVL2 ),
282 1 equemene
     $                   DIFL( NLF, LVL ), DIFR( NLF, LVL2 ),
283 1 equemene
     $                   Z( NLF, LVL ), K( J ), C( J ), S( J ), WORK,
284 1 equemene
     $                   INFO )
285 1 equemene
   30    CONTINUE
286 1 equemene
   40 CONTINUE
287 1 equemene
      GO TO 90
288 1 equemene
*
289 1 equemene
*     ICOMPQ = 1: applying back the right singular vector factors.
290 1 equemene
*
291 1 equemene
   50 CONTINUE
292 1 equemene
*
293 1 equemene
*     First now go through the right singular vector matrices of all
294 1 equemene
*     the tree nodes top-down.
295 1 equemene
*
296 1 equemene
      J = 0
297 1 equemene
      DO 70 LVL = 1, NLVL
298 1 equemene
         LVL2 = 2*LVL - 1
299 1 equemene
*
300 1 equemene
*        Find the first node LF and last node LL on
301 1 equemene
*        the current level LVL.
302 1 equemene
*
303 1 equemene
         IF( LVL.EQ.1 ) THEN
304 1 equemene
            LF = 1
305 1 equemene
            LL = 1
306 1 equemene
         ELSE
307 1 equemene
            LF = 2**( LVL-1 )
308 1 equemene
            LL = 2*LF - 1
309 1 equemene
         END IF
310 1 equemene
         DO 60 I = LL, LF, -1
311 1 equemene
            IM1 = I - 1
312 1 equemene
            IC = IWORK( INODE+IM1 )
313 1 equemene
            NL = IWORK( NDIML+IM1 )
314 1 equemene
            NR = IWORK( NDIMR+IM1 )
315 1 equemene
            NLF = IC - NL
316 1 equemene
            NRF = IC + 1
317 1 equemene
            IF( I.EQ.LL ) THEN
318 1 equemene
               SQRE = 0
319 1 equemene
            ELSE
320 1 equemene
               SQRE = 1
321 1 equemene
            END IF
322 1 equemene
            J = J + 1
323 1 equemene
            CALL DLALS0( ICOMPQ, NL, NR, SQRE, NRHS, B( NLF, 1 ), LDB,
324 1 equemene
     $                   BX( NLF, 1 ), LDBX, PERM( NLF, LVL ),
325 1 equemene
     $                   GIVPTR( J ), GIVCOL( NLF, LVL2 ), LDGCOL,
326 1 equemene
     $                   GIVNUM( NLF, LVL2 ), LDU, POLES( NLF, LVL2 ),
327 1 equemene
     $                   DIFL( NLF, LVL ), DIFR( NLF, LVL2 ),
328 1 equemene
     $                   Z( NLF, LVL ), K( J ), C( J ), S( J ), WORK,
329 1 equemene
     $                   INFO )
330 1 equemene
   60    CONTINUE
331 1 equemene
   70 CONTINUE
332 1 equemene
*
333 1 equemene
*     The nodes on the bottom level of the tree were solved
334 1 equemene
*     by DLASDQ. The corresponding right singular vector
335 1 equemene
*     matrices are in explicit form. Apply them back.
336 1 equemene
*
337 1 equemene
      NDB1 = ( ND+1 ) / 2
338 1 equemene
      DO 80 I = NDB1, ND
339 1 equemene
         I1 = I - 1
340 1 equemene
         IC = IWORK( INODE+I1 )
341 1 equemene
         NL = IWORK( NDIML+I1 )
342 1 equemene
         NR = IWORK( NDIMR+I1 )
343 1 equemene
         NLP1 = NL + 1
344 1 equemene
         IF( I.EQ.ND ) THEN
345 1 equemene
            NRP1 = NR
346 1 equemene
         ELSE
347 1 equemene
            NRP1 = NR + 1
348 1 equemene
         END IF
349 1 equemene
         NLF = IC - NL
350 1 equemene
         NRF = IC + 1
351 1 equemene
         CALL DGEMM( 'T', 'N', NLP1, NRHS, NLP1, ONE, VT( NLF, 1 ), LDU,
352 1 equemene
     $               B( NLF, 1 ), LDB, ZERO, BX( NLF, 1 ), LDBX )
353 1 equemene
         CALL DGEMM( 'T', 'N', NRP1, NRHS, NRP1, ONE, VT( NRF, 1 ), LDU,
354 1 equemene
     $               B( NRF, 1 ), LDB, ZERO, BX( NRF, 1 ), LDBX )
355 1 equemene
   80 CONTINUE
356 1 equemene
*
357 1 equemene
   90 CONTINUE
358 1 equemene
*
359 1 equemene
      RETURN
360 1 equemene
*
361 1 equemene
*     End of DLALSA
362 1 equemene
*
363 1 equemene
      END