LAPACK  3.6.1 LAPACK: Linear Algebra PACKage
zlalsd.f
Go to the documentation of this file.
1 *> \brief \b ZLALSD uses the singular value decomposition of A to solve the least squares problem.
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlalsd.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlalsd.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlalsd.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE ZLALSD( UPLO, SMLSIZ, N, NRHS, D, E, B, LDB, RCOND,
22 * RANK, WORK, RWORK, IWORK, INFO )
23 *
24 * .. Scalar Arguments ..
25 * CHARACTER UPLO
26 * INTEGER INFO, LDB, N, NRHS, RANK, SMLSIZ
27 * DOUBLE PRECISION RCOND
28 * ..
29 * .. Array Arguments ..
30 * INTEGER IWORK( * )
31 * DOUBLE PRECISION D( * ), E( * ), RWORK( * )
32 * COMPLEX*16 B( LDB, * ), WORK( * )
33 * ..
34 *
35 *
36 *> \par Purpose:
37 * =============
38 *>
39 *> \verbatim
40 *>
41 *> ZLALSD uses the singular value decomposition of A to solve the least
42 *> squares problem of finding X to minimize the Euclidean norm of each
43 *> column of A*X-B, where A is N-by-N upper bidiagonal, and X and B
44 *> are N-by-NRHS. The solution X overwrites B.
45 *>
46 *> The singular values of A smaller than RCOND times the largest
47 *> singular value are treated as zero in solving the least squares
48 *> problem; in this case a minimum norm solution is returned.
49 *> The actual singular values are returned in D in ascending order.
50 *>
51 *> This code makes very mild assumptions about floating point
52 *> arithmetic. It will work on machines with a guard digit in
53 *> add/subtract, or on those binary machines without guard digits
54 *> which subtract like the Cray XMP, Cray YMP, Cray C 90, or Cray 2.
55 *> It could conceivably fail on hexadecimal or decimal machines
56 *> without guard digits, but we know of none.
57 *> \endverbatim
58 *
59 * Arguments:
60 * ==========
61 *
62 *> \param[in] UPLO
63 *> \verbatim
64 *> UPLO is CHARACTER*1
65 *> = 'U': D and E define an upper bidiagonal matrix.
66 *> = 'L': D and E define a lower bidiagonal matrix.
67 *> \endverbatim
68 *>
69 *> \param[in] SMLSIZ
70 *> \verbatim
71 *> SMLSIZ is INTEGER
72 *> The maximum size of the subproblems at the bottom of the
73 *> computation tree.
74 *> \endverbatim
75 *>
76 *> \param[in] N
77 *> \verbatim
78 *> N is INTEGER
79 *> The dimension of the bidiagonal matrix. N >= 0.
80 *> \endverbatim
81 *>
82 *> \param[in] NRHS
83 *> \verbatim
84 *> NRHS is INTEGER
85 *> The number of columns of B. NRHS must be at least 1.
86 *> \endverbatim
87 *>
88 *> \param[in,out] D
89 *> \verbatim
90 *> D is DOUBLE PRECISION array, dimension (N)
91 *> On entry D contains the main diagonal of the bidiagonal
92 *> matrix. On exit, if INFO = 0, D contains its singular values.
93 *> \endverbatim
94 *>
95 *> \param[in,out] E
96 *> \verbatim
97 *> E is DOUBLE PRECISION array, dimension (N-1)
98 *> Contains the super-diagonal entries of the bidiagonal matrix.
99 *> On exit, E has been destroyed.
100 *> \endverbatim
101 *>
102 *> \param[in,out] B
103 *> \verbatim
104 *> B is COMPLEX*16 array, dimension (LDB,NRHS)
105 *> On input, B contains the right hand sides of the least
106 *> squares problem. On output, B contains the solution X.
107 *> \endverbatim
108 *>
109 *> \param[in] LDB
110 *> \verbatim
111 *> LDB is INTEGER
112 *> The leading dimension of B in the calling subprogram.
113 *> LDB must be at least max(1,N).
114 *> \endverbatim
115 *>
116 *> \param[in] RCOND
117 *> \verbatim
118 *> RCOND is DOUBLE PRECISION
119 *> The singular values of A less than or equal to RCOND times
120 *> the largest singular value are treated as zero in solving
121 *> the least squares problem. If RCOND is negative,
122 *> machine precision is used instead.
123 *> For example, if diag(S)*X=B were the least squares problem,
124 *> where diag(S) is a diagonal matrix of singular values, the
125 *> solution would be X(i) = B(i) / S(i) if S(i) is greater than
126 *> RCOND*max(S), and X(i) = 0 if S(i) is less than or equal to
127 *> RCOND*max(S).
128 *> \endverbatim
129 *>
130 *> \param[out] RANK
131 *> \verbatim
132 *> RANK is INTEGER
133 *> The number of singular values of A greater than RCOND times
134 *> the largest singular value.
135 *> \endverbatim
136 *>
137 *> \param[out] WORK
138 *> \verbatim
139 *> WORK is COMPLEX*16 array, dimension at least
140 *> (N * NRHS).
141 *> \endverbatim
142 *>
143 *> \param[out] RWORK
144 *> \verbatim
145 *> RWORK is DOUBLE PRECISION array, dimension at least
146 *> (9*N + 2*N*SMLSIZ + 8*N*NLVL + 3*SMLSIZ*NRHS +
147 *> MAX( (SMLSIZ+1)**2, N*(1+NRHS) + 2*NRHS ),
148 *> where
149 *> NLVL = MAX( 0, INT( LOG_2( MIN( M,N )/(SMLSIZ+1) ) ) + 1 )
150 *> \endverbatim
151 *>
152 *> \param[out] IWORK
153 *> \verbatim
154 *> IWORK is INTEGER array, dimension at least
155 *> (3*N*NLVL + 11*N).
156 *> \endverbatim
157 *>
158 *> \param[out] INFO
159 *> \verbatim
160 *> INFO is INTEGER
161 *> = 0: successful exit.
162 *> < 0: if INFO = -i, the i-th argument had an illegal value.
163 *> > 0: The algorithm failed to compute a singular value while
164 *> working on the submatrix lying in rows and columns
165 *> INFO/(N+1) through MOD(INFO,N+1).
166 *> \endverbatim
167 *
168 * Authors:
169 * ========
170 *
171 *> \author Univ. of Tennessee
172 *> \author Univ. of California Berkeley
173 *> \author Univ. of Colorado Denver
174 *> \author NAG Ltd.
175 *
176 *> \date September 2012
177 *
178 *> \ingroup complex16OTHERcomputational
179 *
180 *> \par Contributors:
181 * ==================
182 *>
183 *> Ming Gu and Ren-Cang Li, Computer Science Division, University of
184 *> California at Berkeley, USA \n
185 *> Osni Marques, LBNL/NERSC, USA \n
186 *
187 * =====================================================================
188  SUBROUTINE zlalsd( UPLO, SMLSIZ, N, NRHS, D, E, B, LDB, RCOND,
189  \$ rank, work, rwork, iwork, info )
190 *
191 * -- LAPACK computational routine (version 3.4.2) --
192 * -- LAPACK is a software package provided by Univ. of Tennessee, --
193 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
194 * September 2012
195 *
196 * .. Scalar Arguments ..
197  CHARACTER UPLO
198  INTEGER INFO, LDB, N, NRHS, RANK, SMLSIZ
199  DOUBLE PRECISION RCOND
200 * ..
201 * .. Array Arguments ..
202  INTEGER IWORK( * )
203  DOUBLE PRECISION D( * ), E( * ), RWORK( * )
204  COMPLEX*16 B( ldb, * ), WORK( * )
205 * ..
206 *
207 * =====================================================================
208 *
209 * .. Parameters ..
210  DOUBLE PRECISION ZERO, ONE, TWO
211  parameter ( zero = 0.0d0, one = 1.0d0, two = 2.0d0 )
212  COMPLEX*16 CZERO
213  parameter ( czero = ( 0.0d0, 0.0d0 ) )
214 * ..
215 * .. Local Scalars ..
216  INTEGER BX, BXST, C, DIFL, DIFR, GIVCOL, GIVNUM,
217  \$ givptr, i, icmpq1, icmpq2, irwb, irwib, irwrb,
218  \$ irwu, irwvt, irwwrk, iwk, j, jcol, jimag,
219  \$ jreal, jrow, k, nlvl, nm1, nrwork, nsize, nsub,
220  \$ perm, poles, s, sizei, smlszp, sqre, st, st1,
221  \$ u, vt, z
222  DOUBLE PRECISION CS, EPS, ORGNRM, RCND, R, SN, TOL
223 * ..
224 * .. External Functions ..
225  INTEGER IDAMAX
226  DOUBLE PRECISION DLAMCH, DLANST
227  EXTERNAL idamax, dlamch, dlanst
228 * ..
229 * .. External Subroutines ..
230  EXTERNAL dgemm, dlartg, dlascl, dlasda, dlasdq, dlaset,
232  \$ zlascl, zlaset
233 * ..
234 * .. Intrinsic Functions ..
235  INTRINSIC abs, dble, dcmplx, dimag, int, log, sign
236 * ..
237 * .. Executable Statements ..
238 *
239 * Test the input parameters.
240 *
241  info = 0
242 *
243  IF( n.LT.0 ) THEN
244  info = -3
245  ELSE IF( nrhs.LT.1 ) THEN
246  info = -4
247  ELSE IF( ( ldb.LT.1 ) .OR. ( ldb.LT.n ) ) THEN
248  info = -8
249  END IF
250  IF( info.NE.0 ) THEN
251  CALL xerbla( 'ZLALSD', -info )
252  RETURN
253  END IF
254 *
255  eps = dlamch( 'Epsilon' )
256 *
257 * Set up the tolerance.
258 *
259  IF( ( rcond.LE.zero ) .OR. ( rcond.GE.one ) ) THEN
260  rcnd = eps
261  ELSE
262  rcnd = rcond
263  END IF
264 *
265  rank = 0
266 *
267 * Quick return if possible.
268 *
269  IF( n.EQ.0 ) THEN
270  RETURN
271  ELSE IF( n.EQ.1 ) THEN
272  IF( d( 1 ).EQ.zero ) THEN
273  CALL zlaset( 'A', 1, nrhs, czero, czero, b, ldb )
274  ELSE
275  rank = 1
276  CALL zlascl( 'G', 0, 0, d( 1 ), one, 1, nrhs, b, ldb, info )
277  d( 1 ) = abs( d( 1 ) )
278  END IF
279  RETURN
280  END IF
281 *
282 * Rotate the matrix if it is lower bidiagonal.
283 *
284  IF( uplo.EQ.'L' ) THEN
285  DO 10 i = 1, n - 1
286  CALL dlartg( d( i ), e( i ), cs, sn, r )
287  d( i ) = r
288  e( i ) = sn*d( i+1 )
289  d( i+1 ) = cs*d( i+1 )
290  IF( nrhs.EQ.1 ) THEN
291  CALL zdrot( 1, b( i, 1 ), 1, b( i+1, 1 ), 1, cs, sn )
292  ELSE
293  rwork( i*2-1 ) = cs
294  rwork( i*2 ) = sn
295  END IF
296  10 CONTINUE
297  IF( nrhs.GT.1 ) THEN
298  DO 30 i = 1, nrhs
299  DO 20 j = 1, n - 1
300  cs = rwork( j*2-1 )
301  sn = rwork( j*2 )
302  CALL zdrot( 1, b( j, i ), 1, b( j+1, i ), 1, cs, sn )
303  20 CONTINUE
304  30 CONTINUE
305  END IF
306  END IF
307 *
308 * Scale.
309 *
310  nm1 = n - 1
311  orgnrm = dlanst( 'M', n, d, e )
312  IF( orgnrm.EQ.zero ) THEN
313  CALL zlaset( 'A', n, nrhs, czero, czero, b, ldb )
314  RETURN
315  END IF
316 *
317  CALL dlascl( 'G', 0, 0, orgnrm, one, n, 1, d, n, info )
318  CALL dlascl( 'G', 0, 0, orgnrm, one, nm1, 1, e, nm1, info )
319 *
320 * If N is smaller than the minimum divide size SMLSIZ, then solve
321 * the problem with another solver.
322 *
323  IF( n.LE.smlsiz ) THEN
324  irwu = 1
325  irwvt = irwu + n*n
326  irwwrk = irwvt + n*n
327  irwrb = irwwrk
328  irwib = irwrb + n*nrhs
329  irwb = irwib + n*nrhs
330  CALL dlaset( 'A', n, n, zero, one, rwork( irwu ), n )
331  CALL dlaset( 'A', n, n, zero, one, rwork( irwvt ), n )
332  CALL dlasdq( 'U', 0, n, n, n, 0, d, e, rwork( irwvt ), n,
333  \$ rwork( irwu ), n, rwork( irwwrk ), 1,
334  \$ rwork( irwwrk ), info )
335  IF( info.NE.0 ) THEN
336  RETURN
337  END IF
338 *
339 * In the real version, B is passed to DLASDQ and multiplied
340 * internally by Q**H. Here B is complex and that product is
341 * computed below in two steps (real and imaginary parts).
342 *
343  j = irwb - 1
344  DO 50 jcol = 1, nrhs
345  DO 40 jrow = 1, n
346  j = j + 1
347  rwork( j ) = dble( b( jrow, jcol ) )
348  40 CONTINUE
349  50 CONTINUE
350  CALL dgemm( 'T', 'N', n, nrhs, n, one, rwork( irwu ), n,
351  \$ rwork( irwb ), n, zero, rwork( irwrb ), n )
352  j = irwb - 1
353  DO 70 jcol = 1, nrhs
354  DO 60 jrow = 1, n
355  j = j + 1
356  rwork( j ) = dimag( b( jrow, jcol ) )
357  60 CONTINUE
358  70 CONTINUE
359  CALL dgemm( 'T', 'N', n, nrhs, n, one, rwork( irwu ), n,
360  \$ rwork( irwb ), n, zero, rwork( irwib ), n )
361  jreal = irwrb - 1
362  jimag = irwib - 1
363  DO 90 jcol = 1, nrhs
364  DO 80 jrow = 1, n
365  jreal = jreal + 1
366  jimag = jimag + 1
367  b( jrow, jcol ) = dcmplx( rwork( jreal ),
368  \$ rwork( jimag ) )
369  80 CONTINUE
370  90 CONTINUE
371 *
372  tol = rcnd*abs( d( idamax( n, d, 1 ) ) )
373  DO 100 i = 1, n
374  IF( d( i ).LE.tol ) THEN
375  CALL zlaset( 'A', 1, nrhs, czero, czero, b( i, 1 ), ldb )
376  ELSE
377  CALL zlascl( 'G', 0, 0, d( i ), one, 1, nrhs, b( i, 1 ),
378  \$ ldb, info )
379  rank = rank + 1
380  END IF
381  100 CONTINUE
382 *
383 * Since B is complex, the following call to DGEMM is performed
384 * in two steps (real and imaginary parts). That is for V * B
385 * (in the real version of the code V**H is stored in WORK).
386 *
387 * CALL DGEMM( 'T', 'N', N, NRHS, N, ONE, WORK, N, B, LDB, ZERO,
388 * \$ WORK( NWORK ), N )
389 *
390  j = irwb - 1
391  DO 120 jcol = 1, nrhs
392  DO 110 jrow = 1, n
393  j = j + 1
394  rwork( j ) = dble( b( jrow, jcol ) )
395  110 CONTINUE
396  120 CONTINUE
397  CALL dgemm( 'T', 'N', n, nrhs, n, one, rwork( irwvt ), n,
398  \$ rwork( irwb ), n, zero, rwork( irwrb ), n )
399  j = irwb - 1
400  DO 140 jcol = 1, nrhs
401  DO 130 jrow = 1, n
402  j = j + 1
403  rwork( j ) = dimag( b( jrow, jcol ) )
404  130 CONTINUE
405  140 CONTINUE
406  CALL dgemm( 'T', 'N', n, nrhs, n, one, rwork( irwvt ), n,
407  \$ rwork( irwb ), n, zero, rwork( irwib ), n )
408  jreal = irwrb - 1
409  jimag = irwib - 1
410  DO 160 jcol = 1, nrhs
411  DO 150 jrow = 1, n
412  jreal = jreal + 1
413  jimag = jimag + 1
414  b( jrow, jcol ) = dcmplx( rwork( jreal ),
415  \$ rwork( jimag ) )
416  150 CONTINUE
417  160 CONTINUE
418 *
419 * Unscale.
420 *
421  CALL dlascl( 'G', 0, 0, one, orgnrm, n, 1, d, n, info )
422  CALL dlasrt( 'D', n, d, info )
423  CALL zlascl( 'G', 0, 0, orgnrm, one, n, nrhs, b, ldb, info )
424 *
425  RETURN
426  END IF
427 *
428 * Book-keeping and setting up some constants.
429 *
430  nlvl = int( log( dble( n ) / dble( smlsiz+1 ) ) / log( two ) ) + 1
431 *
432  smlszp = smlsiz + 1
433 *
434  u = 1
435  vt = 1 + smlsiz*n
436  difl = vt + smlszp*n
437  difr = difl + nlvl*n
438  z = difr + nlvl*n*2
439  c = z + nlvl*n
440  s = c + n
441  poles = s + n
442  givnum = poles + 2*nlvl*n
443  nrwork = givnum + 2*nlvl*n
444  bx = 1
445 *
446  irwrb = nrwork
447  irwib = irwrb + smlsiz*nrhs
448  irwb = irwib + smlsiz*nrhs
449 *
450  sizei = 1 + n
451  k = sizei + n
452  givptr = k + n
453  perm = givptr + n
454  givcol = perm + nlvl*n
455  iwk = givcol + nlvl*n*2
456 *
457  st = 1
458  sqre = 0
459  icmpq1 = 1
460  icmpq2 = 0
461  nsub = 0
462 *
463  DO 170 i = 1, n
464  IF( abs( d( i ) ).LT.eps ) THEN
465  d( i ) = sign( eps, d( i ) )
466  END IF
467  170 CONTINUE
468 *
469  DO 240 i = 1, nm1
470  IF( ( abs( e( i ) ).LT.eps ) .OR. ( i.EQ.nm1 ) ) THEN
471  nsub = nsub + 1
472  iwork( nsub ) = st
473 *
474 * Subproblem found. First determine its size and then
475 * apply divide and conquer on it.
476 *
477  IF( i.LT.nm1 ) THEN
478 *
479 * A subproblem with E(I) small for I < NM1.
480 *
481  nsize = i - st + 1
482  iwork( sizei+nsub-1 ) = nsize
483  ELSE IF( abs( e( i ) ).GE.eps ) THEN
484 *
485 * A subproblem with E(NM1) not too small but I = NM1.
486 *
487  nsize = n - st + 1
488  iwork( sizei+nsub-1 ) = nsize
489  ELSE
490 *
491 * A subproblem with E(NM1) small. This implies an
492 * 1-by-1 subproblem at D(N), which is not solved
493 * explicitly.
494 *
495  nsize = i - st + 1
496  iwork( sizei+nsub-1 ) = nsize
497  nsub = nsub + 1
498  iwork( nsub ) = n
499  iwork( sizei+nsub-1 ) = 1
500  CALL zcopy( nrhs, b( n, 1 ), ldb, work( bx+nm1 ), n )
501  END IF
502  st1 = st - 1
503  IF( nsize.EQ.1 ) THEN
504 *
505 * This is a 1-by-1 subproblem and is not solved
506 * explicitly.
507 *
508  CALL zcopy( nrhs, b( st, 1 ), ldb, work( bx+st1 ), n )
509  ELSE IF( nsize.LE.smlsiz ) THEN
510 *
511 * This is a small subproblem and is solved by DLASDQ.
512 *
513  CALL dlaset( 'A', nsize, nsize, zero, one,
514  \$ rwork( vt+st1 ), n )
515  CALL dlaset( 'A', nsize, nsize, zero, one,
516  \$ rwork( u+st1 ), n )
517  CALL dlasdq( 'U', 0, nsize, nsize, nsize, 0, d( st ),
518  \$ e( st ), rwork( vt+st1 ), n, rwork( u+st1 ),
519  \$ n, rwork( nrwork ), 1, rwork( nrwork ),
520  \$ info )
521  IF( info.NE.0 ) THEN
522  RETURN
523  END IF
524 *
525 * In the real version, B is passed to DLASDQ and multiplied
526 * internally by Q**H. Here B is complex and that product is
527 * computed below in two steps (real and imaginary parts).
528 *
529  j = irwb - 1
530  DO 190 jcol = 1, nrhs
531  DO 180 jrow = st, st + nsize - 1
532  j = j + 1
533  rwork( j ) = dble( b( jrow, jcol ) )
534  180 CONTINUE
535  190 CONTINUE
536  CALL dgemm( 'T', 'N', nsize, nrhs, nsize, one,
537  \$ rwork( u+st1 ), n, rwork( irwb ), nsize,
538  \$ zero, rwork( irwrb ), nsize )
539  j = irwb - 1
540  DO 210 jcol = 1, nrhs
541  DO 200 jrow = st, st + nsize - 1
542  j = j + 1
543  rwork( j ) = dimag( b( jrow, jcol ) )
544  200 CONTINUE
545  210 CONTINUE
546  CALL dgemm( 'T', 'N', nsize, nrhs, nsize, one,
547  \$ rwork( u+st1 ), n, rwork( irwb ), nsize,
548  \$ zero, rwork( irwib ), nsize )
549  jreal = irwrb - 1
550  jimag = irwib - 1
551  DO 230 jcol = 1, nrhs
552  DO 220 jrow = st, st + nsize - 1
553  jreal = jreal + 1
554  jimag = jimag + 1
555  b( jrow, jcol ) = dcmplx( rwork( jreal ),
556  \$ rwork( jimag ) )
557  220 CONTINUE
558  230 CONTINUE
559 *
560  CALL zlacpy( 'A', nsize, nrhs, b( st, 1 ), ldb,
561  \$ work( bx+st1 ), n )
562  ELSE
563 *
564 * A large problem. Solve it using divide and conquer.
565 *
566  CALL dlasda( icmpq1, smlsiz, nsize, sqre, d( st ),
567  \$ e( st ), rwork( u+st1 ), n, rwork( vt+st1 ),
568  \$ iwork( k+st1 ), rwork( difl+st1 ),
569  \$ rwork( difr+st1 ), rwork( z+st1 ),
570  \$ rwork( poles+st1 ), iwork( givptr+st1 ),
571  \$ iwork( givcol+st1 ), n, iwork( perm+st1 ),
572  \$ rwork( givnum+st1 ), rwork( c+st1 ),
573  \$ rwork( s+st1 ), rwork( nrwork ),
574  \$ iwork( iwk ), info )
575  IF( info.NE.0 ) THEN
576  RETURN
577  END IF
578  bxst = bx + st1
579  CALL zlalsa( icmpq2, smlsiz, nsize, nrhs, b( st, 1 ),
580  \$ ldb, work( bxst ), n, rwork( u+st1 ), n,
581  \$ rwork( vt+st1 ), iwork( k+st1 ),
582  \$ rwork( difl+st1 ), rwork( difr+st1 ),
583  \$ rwork( z+st1 ), rwork( poles+st1 ),
584  \$ iwork( givptr+st1 ), iwork( givcol+st1 ), n,
585  \$ iwork( perm+st1 ), rwork( givnum+st1 ),
586  \$ rwork( c+st1 ), rwork( s+st1 ),
587  \$ rwork( nrwork ), iwork( iwk ), info )
588  IF( info.NE.0 ) THEN
589  RETURN
590  END IF
591  END IF
592  st = i + 1
593  END IF
594  240 CONTINUE
595 *
596 * Apply the singular values and treat the tiny ones as zero.
597 *
598  tol = rcnd*abs( d( idamax( n, d, 1 ) ) )
599 *
600  DO 250 i = 1, n
601 *
602 * Some of the elements in D can be negative because 1-by-1
603 * subproblems were not solved explicitly.
604 *
605  IF( abs( d( i ) ).LE.tol ) THEN
606  CALL zlaset( 'A', 1, nrhs, czero, czero, work( bx+i-1 ), n )
607  ELSE
608  rank = rank + 1
609  CALL zlascl( 'G', 0, 0, d( i ), one, 1, nrhs,
610  \$ work( bx+i-1 ), n, info )
611  END IF
612  d( i ) = abs( d( i ) )
613  250 CONTINUE
614 *
615 * Now apply back the right singular vectors.
616 *
617  icmpq2 = 1
618  DO 320 i = 1, nsub
619  st = iwork( i )
620  st1 = st - 1
621  nsize = iwork( sizei+i-1 )
622  bxst = bx + st1
623  IF( nsize.EQ.1 ) THEN
624  CALL zcopy( nrhs, work( bxst ), n, b( st, 1 ), ldb )
625  ELSE IF( nsize.LE.smlsiz ) THEN
626 *
627 * Since B and BX are complex, the following call to DGEMM
628 * is performed in two steps (real and imaginary parts).
629 *
630 * CALL DGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE,
631 * \$ RWORK( VT+ST1 ), N, RWORK( BXST ), N, ZERO,
632 * \$ B( ST, 1 ), LDB )
633 *
634  j = bxst - n - 1
635  jreal = irwb - 1
636  DO 270 jcol = 1, nrhs
637  j = j + n
638  DO 260 jrow = 1, nsize
639  jreal = jreal + 1
640  rwork( jreal ) = dble( work( j+jrow ) )
641  260 CONTINUE
642  270 CONTINUE
643  CALL dgemm( 'T', 'N', nsize, nrhs, nsize, one,
644  \$ rwork( vt+st1 ), n, rwork( irwb ), nsize, zero,
645  \$ rwork( irwrb ), nsize )
646  j = bxst - n - 1
647  jimag = irwb - 1
648  DO 290 jcol = 1, nrhs
649  j = j + n
650  DO 280 jrow = 1, nsize
651  jimag = jimag + 1
652  rwork( jimag ) = dimag( work( j+jrow ) )
653  280 CONTINUE
654  290 CONTINUE
655  CALL dgemm( 'T', 'N', nsize, nrhs, nsize, one,
656  \$ rwork( vt+st1 ), n, rwork( irwb ), nsize, zero,
657  \$ rwork( irwib ), nsize )
658  jreal = irwrb - 1
659  jimag = irwib - 1
660  DO 310 jcol = 1, nrhs
661  DO 300 jrow = st, st + nsize - 1
662  jreal = jreal + 1
663  jimag = jimag + 1
664  b( jrow, jcol ) = dcmplx( rwork( jreal ),
665  \$ rwork( jimag ) )
666  300 CONTINUE
667  310 CONTINUE
668  ELSE
669  CALL zlalsa( icmpq2, smlsiz, nsize, nrhs, work( bxst ), n,
670  \$ b( st, 1 ), ldb, rwork( u+st1 ), n,
671  \$ rwork( vt+st1 ), iwork( k+st1 ),
672  \$ rwork( difl+st1 ), rwork( difr+st1 ),
673  \$ rwork( z+st1 ), rwork( poles+st1 ),
674  \$ iwork( givptr+st1 ), iwork( givcol+st1 ), n,
675  \$ iwork( perm+st1 ), rwork( givnum+st1 ),
676  \$ rwork( c+st1 ), rwork( s+st1 ),
677  \$ rwork( nrwork ), iwork( iwk ), info )
678  IF( info.NE.0 ) THEN
679  RETURN
680  END IF
681  END IF
682  320 CONTINUE
683 *
684 * Unscale and sort the singular values.
685 *
686  CALL dlascl( 'G', 0, 0, one, orgnrm, n, 1, d, n, info )
687  CALL dlasrt( 'D', n, d, info )
688  CALL zlascl( 'G', 0, 0, orgnrm, one, n, nrhs, b, ldb, info )
689 *
690  RETURN
691 *
692 * End of ZLALSD
693 *
694  END
subroutine dlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
Definition: dlaset.f:112
subroutine dlasrt(ID, N, D, INFO)
DLASRT sorts numbers in increasing or decreasing order.
Definition: dlasrt.f:90
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
Definition: zlacpy.f:105
subroutine zlalsa(ICOMPQ, SMLSIZ, N, NRHS, B, LDB, BX, LDBX, U, LDU, VT, K, DIFL, DIFR, Z, POLES, GIVPTR, GIVCOL, LDGCOL, PERM, GIVNUM, C, S, RWORK, IWORK, INFO)
ZLALSA computes the SVD of the coefficient matrix in compact form. Used by sgelsd.
Definition: zlalsa.f:270
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
Definition: zcopy.f:52
subroutine dlasda(ICOMPQ, SMLSIZ, N, SQRE, D, E, U, LDU, VT, K, DIFL, DIFR, Z, POLES, GIVPTR, GIVCOL, LDGCOL, PERM, GIVNUM, C, S, WORK, IWORK, INFO)
DLASDA computes the singular value decomposition (SVD) of a real upper bidiagonal matrix with diagona...
Definition: dlasda.f:276
subroutine dlascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
DLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: dlascl.f:145
subroutine zlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
Definition: zlaset.f:108
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
Definition: dgemm.f:189
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine zdrot(N, CX, INCX, CY, INCY, C, S)
ZDROT
Definition: zdrot.f:100
subroutine dlasdq(UPLO, SQRE, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C, LDC, WORK, INFO)
DLASDQ computes the SVD of a real bidiagonal matrix with diagonal d and off-diagonal e...
Definition: dlasdq.f:213
subroutine zlascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
ZLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: zlascl.f:145
subroutine dlartg(F, G, CS, SN, R)
DLARTG generates a plane rotation with real cosine and real sine.
Definition: dlartg.f:99
subroutine zlalsd(UPLO, SMLSIZ, N, NRHS, D, E, B, LDB, RCOND, RANK, WORK, RWORK, IWORK, INFO)
ZLALSD uses the singular value decomposition of A to solve the least squares problem.
Definition: zlalsd.f:190