LAPACK  3.9.0
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
9 *> Download ZLALSD + dependencies
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 (N * NRHS)
140 *> \endverbatim
141 *>
142 *> \param[out] RWORK
143 *> \verbatim
144 *> RWORK is DOUBLE PRECISION array, dimension at least
145 *> (9*N + 2*N*SMLSIZ + 8*N*NLVL + 3*SMLSIZ*NRHS +
146 *> MAX( (SMLSIZ+1)**2, N*(1+NRHS) + 2*NRHS ),
147 *> where
148 *> NLVL = MAX( 0, INT( LOG_2( MIN( M,N )/(SMLSIZ+1) ) ) + 1 )
149 *> \endverbatim
150 *>
151 *> \param[out] IWORK
152 *> \verbatim
153 *> IWORK is INTEGER array, dimension at least
154 *> (3*N*NLVL + 11*N).
155 *> \endverbatim
156 *>
157 *> \param[out] INFO
158 *> \verbatim
159 *> INFO is INTEGER
160 *> = 0: successful exit.
161 *> < 0: if INFO = -i, the i-th argument had an illegal value.
162 *> > 0: The algorithm failed to compute a singular value while
163 *> working on the submatrix lying in rows and columns
164 *> INFO/(N+1) through MOD(INFO,N+1).
165 *> \endverbatim
166 *
167 * Authors:
168 * ========
169 *
170 *> \author Univ. of Tennessee
171 *> \author Univ. of California Berkeley
172 *> \author Univ. of Colorado Denver
173 *> \author NAG Ltd.
174 *
175 *> \date June 2017
176 *
177 *> \ingroup complex16OTHERcomputational
178 *
179 *> \par Contributors:
180 * ==================
181 *>
182 *> Ming Gu and Ren-Cang Li, Computer Science Division, University of
183 *> California at Berkeley, USA \n
184 *> Osni Marques, LBNL/NERSC, USA \n
185 *
186 * =====================================================================
187  SUBROUTINE zlalsd( UPLO, SMLSIZ, N, NRHS, D, E, B, LDB, RCOND,
188  $ RANK, WORK, RWORK, IWORK, INFO )
189 *
190 * -- LAPACK computational routine (version 3.7.1) --
191 * -- LAPACK is a software package provided by Univ. of Tennessee, --
192 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
193 * June 2017
194 *
195 * .. Scalar Arguments ..
196  CHARACTER UPLO
197  INTEGER INFO, LDB, N, NRHS, RANK, SMLSIZ
198  DOUBLE PRECISION RCOND
199 * ..
200 * .. Array Arguments ..
201  INTEGER IWORK( * )
202  DOUBLE PRECISION D( * ), E( * ), RWORK( * )
203  COMPLEX*16 B( LDB, * ), WORK( * )
204 * ..
205 *
206 * =====================================================================
207 *
208 * .. Parameters ..
209  DOUBLE PRECISION ZERO, ONE, TWO
210  parameter( zero = 0.0d0, one = 1.0d0, two = 2.0d0 )
211  COMPLEX*16 CZERO
212  parameter( czero = ( 0.0d0, 0.0d0 ) )
213 * ..
214 * .. Local Scalars ..
215  INTEGER BX, BXST, C, DIFL, DIFR, GIVCOL, GIVNUM,
216  $ givptr, i, icmpq1, icmpq2, irwb, irwib, irwrb,
217  $ irwu, irwvt, irwwrk, iwk, j, jcol, jimag,
218  $ jreal, jrow, k, nlvl, nm1, nrwork, nsize, nsub,
219  $ perm, poles, s, sizei, smlszp, sqre, st, st1,
220  $ u, vt, z
221  DOUBLE PRECISION CS, EPS, ORGNRM, RCND, R, SN, TOL
222 * ..
223 * .. External Functions ..
224  INTEGER IDAMAX
225  DOUBLE PRECISION DLAMCH, DLANST
226  EXTERNAL idamax, dlamch, dlanst
227 * ..
228 * .. External Subroutines ..
229  EXTERNAL dgemm, dlartg, dlascl, dlasda, dlasdq, dlaset,
230  $ dlasrt, xerbla, zcopy, zdrot, zlacpy, zlalsa,
231  $ zlascl, zlaset
232 * ..
233 * .. Intrinsic Functions ..
234  INTRINSIC abs, dble, dcmplx, dimag, int, log, sign
235 * ..
236 * .. Executable Statements ..
237 *
238 * Test the input parameters.
239 *
240  info = 0
241 *
242  IF( n.LT.0 ) THEN
243  info = -3
244  ELSE IF( nrhs.LT.1 ) THEN
245  info = -4
246  ELSE IF( ( ldb.LT.1 ) .OR. ( ldb.LT.n ) ) THEN
247  info = -8
248  END IF
249  IF( info.NE.0 ) THEN
250  CALL xerbla( 'ZLALSD', -info )
251  RETURN
252  END IF
253 *
254  eps = dlamch( 'Epsilon' )
255 *
256 * Set up the tolerance.
257 *
258  IF( ( rcond.LE.zero ) .OR. ( rcond.GE.one ) ) THEN
259  rcnd = eps
260  ELSE
261  rcnd = rcond
262  END IF
263 *
264  rank = 0
265 *
266 * Quick return if possible.
267 *
268  IF( n.EQ.0 ) THEN
269  RETURN
270  ELSE IF( n.EQ.1 ) THEN
271  IF( d( 1 ).EQ.zero ) THEN
272  CALL zlaset( 'A', 1, nrhs, czero, czero, b, ldb )
273  ELSE
274  rank = 1
275  CALL zlascl( 'G', 0, 0, d( 1 ), one, 1, nrhs, b, ldb, info )
276  d( 1 ) = abs( d( 1 ) )
277  END IF
278  RETURN
279  END IF
280 *
281 * Rotate the matrix if it is lower bidiagonal.
282 *
283  IF( uplo.EQ.'L' ) THEN
284  DO 10 i = 1, n - 1
285  CALL dlartg( d( i ), e( i ), cs, sn, r )
286  d( i ) = r
287  e( i ) = sn*d( i+1 )
288  d( i+1 ) = cs*d( i+1 )
289  IF( nrhs.EQ.1 ) THEN
290  CALL zdrot( 1, b( i, 1 ), 1, b( i+1, 1 ), 1, cs, sn )
291  ELSE
292  rwork( i*2-1 ) = cs
293  rwork( i*2 ) = sn
294  END IF
295  10 CONTINUE
296  IF( nrhs.GT.1 ) THEN
297  DO 30 i = 1, nrhs
298  DO 20 j = 1, n - 1
299  cs = rwork( j*2-1 )
300  sn = rwork( j*2 )
301  CALL zdrot( 1, b( j, i ), 1, b( j+1, i ), 1, cs, sn )
302  20 CONTINUE
303  30 CONTINUE
304  END IF
305  END IF
306 *
307 * Scale.
308 *
309  nm1 = n - 1
310  orgnrm = dlanst( 'M', n, d, e )
311  IF( orgnrm.EQ.zero ) THEN
312  CALL zlaset( 'A', n, nrhs, czero, czero, b, ldb )
313  RETURN
314  END IF
315 *
316  CALL dlascl( 'G', 0, 0, orgnrm, one, n, 1, d, n, info )
317  CALL dlascl( 'G', 0, 0, orgnrm, one, nm1, 1, e, nm1, info )
318 *
319 * If N is smaller than the minimum divide size SMLSIZ, then solve
320 * the problem with another solver.
321 *
322  IF( n.LE.smlsiz ) THEN
323  irwu = 1
324  irwvt = irwu + n*n
325  irwwrk = irwvt + n*n
326  irwrb = irwwrk
327  irwib = irwrb + n*nrhs
328  irwb = irwib + n*nrhs
329  CALL dlaset( 'A', n, n, zero, one, rwork( irwu ), n )
330  CALL dlaset( 'A', n, n, zero, one, rwork( irwvt ), n )
331  CALL dlasdq( 'U', 0, n, n, n, 0, d, e, rwork( irwvt ), n,
332  $ rwork( irwu ), n, rwork( irwwrk ), 1,
333  $ rwork( irwwrk ), info )
334  IF( info.NE.0 ) THEN
335  RETURN
336  END IF
337 *
338 * In the real version, B is passed to DLASDQ and multiplied
339 * internally by Q**H. Here B is complex and that product is
340 * computed below in two steps (real and imaginary parts).
341 *
342  j = irwb - 1
343  DO 50 jcol = 1, nrhs
344  DO 40 jrow = 1, n
345  j = j + 1
346  rwork( j ) = dble( b( jrow, jcol ) )
347  40 CONTINUE
348  50 CONTINUE
349  CALL dgemm( 'T', 'N', n, nrhs, n, one, rwork( irwu ), n,
350  $ rwork( irwb ), n, zero, rwork( irwrb ), n )
351  j = irwb - 1
352  DO 70 jcol = 1, nrhs
353  DO 60 jrow = 1, n
354  j = j + 1
355  rwork( j ) = dimag( b( jrow, jcol ) )
356  60 CONTINUE
357  70 CONTINUE
358  CALL dgemm( 'T', 'N', n, nrhs, n, one, rwork( irwu ), n,
359  $ rwork( irwb ), n, zero, rwork( irwib ), n )
360  jreal = irwrb - 1
361  jimag = irwib - 1
362  DO 90 jcol = 1, nrhs
363  DO 80 jrow = 1, n
364  jreal = jreal + 1
365  jimag = jimag + 1
366  b( jrow, jcol ) = dcmplx( rwork( jreal ),
367  $ rwork( jimag ) )
368  80 CONTINUE
369  90 CONTINUE
370 *
371  tol = rcnd*abs( d( idamax( n, d, 1 ) ) )
372  DO 100 i = 1, n
373  IF( d( i ).LE.tol ) THEN
374  CALL zlaset( 'A', 1, nrhs, czero, czero, b( i, 1 ), ldb )
375  ELSE
376  CALL zlascl( 'G', 0, 0, d( i ), one, 1, nrhs, b( i, 1 ),
377  $ ldb, info )
378  rank = rank + 1
379  END IF
380  100 CONTINUE
381 *
382 * Since B is complex, the following call to DGEMM is performed
383 * in two steps (real and imaginary parts). That is for V * B
384 * (in the real version of the code V**H is stored in WORK).
385 *
386 * CALL DGEMM( 'T', 'N', N, NRHS, N, ONE, WORK, N, B, LDB, ZERO,
387 * $ WORK( NWORK ), N )
388 *
389  j = irwb - 1
390  DO 120 jcol = 1, nrhs
391  DO 110 jrow = 1, n
392  j = j + 1
393  rwork( j ) = dble( b( jrow, jcol ) )
394  110 CONTINUE
395  120 CONTINUE
396  CALL dgemm( 'T', 'N', n, nrhs, n, one, rwork( irwvt ), n,
397  $ rwork( irwb ), n, zero, rwork( irwrb ), n )
398  j = irwb - 1
399  DO 140 jcol = 1, nrhs
400  DO 130 jrow = 1, n
401  j = j + 1
402  rwork( j ) = dimag( b( jrow, jcol ) )
403  130 CONTINUE
404  140 CONTINUE
405  CALL dgemm( 'T', 'N', n, nrhs, n, one, rwork( irwvt ), n,
406  $ rwork( irwb ), n, zero, rwork( irwib ), n )
407  jreal = irwrb - 1
408  jimag = irwib - 1
409  DO 160 jcol = 1, nrhs
410  DO 150 jrow = 1, n
411  jreal = jreal + 1
412  jimag = jimag + 1
413  b( jrow, jcol ) = dcmplx( rwork( jreal ),
414  $ rwork( jimag ) )
415  150 CONTINUE
416  160 CONTINUE
417 *
418 * Unscale.
419 *
420  CALL dlascl( 'G', 0, 0, one, orgnrm, n, 1, d, n, info )
421  CALL dlasrt( 'D', n, d, info )
422  CALL zlascl( 'G', 0, 0, orgnrm, one, n, nrhs, b, ldb, info )
423 *
424  RETURN
425  END IF
426 *
427 * Book-keeping and setting up some constants.
428 *
429  nlvl = int( log( dble( n ) / dble( smlsiz+1 ) ) / log( two ) ) + 1
430 *
431  smlszp = smlsiz + 1
432 *
433  u = 1
434  vt = 1 + smlsiz*n
435  difl = vt + smlszp*n
436  difr = difl + nlvl*n
437  z = difr + nlvl*n*2
438  c = z + nlvl*n
439  s = c + n
440  poles = s + n
441  givnum = poles + 2*nlvl*n
442  nrwork = givnum + 2*nlvl*n
443  bx = 1
444 *
445  irwrb = nrwork
446  irwib = irwrb + smlsiz*nrhs
447  irwb = irwib + smlsiz*nrhs
448 *
449  sizei = 1 + n
450  k = sizei + n
451  givptr = k + n
452  perm = givptr + n
453  givcol = perm + nlvl*n
454  iwk = givcol + nlvl*n*2
455 *
456  st = 1
457  sqre = 0
458  icmpq1 = 1
459  icmpq2 = 0
460  nsub = 0
461 *
462  DO 170 i = 1, n
463  IF( abs( d( i ) ).LT.eps ) THEN
464  d( i ) = sign( eps, d( i ) )
465  END IF
466  170 CONTINUE
467 *
468  DO 240 i = 1, nm1
469  IF( ( abs( e( i ) ).LT.eps ) .OR. ( i.EQ.nm1 ) ) THEN
470  nsub = nsub + 1
471  iwork( nsub ) = st
472 *
473 * Subproblem found. First determine its size and then
474 * apply divide and conquer on it.
475 *
476  IF( i.LT.nm1 ) THEN
477 *
478 * A subproblem with E(I) small for I < NM1.
479 *
480  nsize = i - st + 1
481  iwork( sizei+nsub-1 ) = nsize
482  ELSE IF( abs( e( i ) ).GE.eps ) THEN
483 *
484 * A subproblem with E(NM1) not too small but I = NM1.
485 *
486  nsize = n - st + 1
487  iwork( sizei+nsub-1 ) = nsize
488  ELSE
489 *
490 * A subproblem with E(NM1) small. This implies an
491 * 1-by-1 subproblem at D(N), which is not solved
492 * explicitly.
493 *
494  nsize = i - st + 1
495  iwork( sizei+nsub-1 ) = nsize
496  nsub = nsub + 1
497  iwork( nsub ) = n
498  iwork( sizei+nsub-1 ) = 1
499  CALL zcopy( nrhs, b( n, 1 ), ldb, work( bx+nm1 ), n )
500  END IF
501  st1 = st - 1
502  IF( nsize.EQ.1 ) THEN
503 *
504 * This is a 1-by-1 subproblem and is not solved
505 * explicitly.
506 *
507  CALL zcopy( nrhs, b( st, 1 ), ldb, work( bx+st1 ), n )
508  ELSE IF( nsize.LE.smlsiz ) THEN
509 *
510 * This is a small subproblem and is solved by DLASDQ.
511 *
512  CALL dlaset( 'A', nsize, nsize, zero, one,
513  $ rwork( vt+st1 ), n )
514  CALL dlaset( 'A', nsize, nsize, zero, one,
515  $ rwork( u+st1 ), n )
516  CALL dlasdq( 'U', 0, nsize, nsize, nsize, 0, d( st ),
517  $ e( st ), rwork( vt+st1 ), n, rwork( u+st1 ),
518  $ n, rwork( nrwork ), 1, rwork( nrwork ),
519  $ info )
520  IF( info.NE.0 ) THEN
521  RETURN
522  END IF
523 *
524 * In the real version, B is passed to DLASDQ and multiplied
525 * internally by Q**H. Here B is complex and that product is
526 * computed below in two steps (real and imaginary parts).
527 *
528  j = irwb - 1
529  DO 190 jcol = 1, nrhs
530  DO 180 jrow = st, st + nsize - 1
531  j = j + 1
532  rwork( j ) = dble( b( jrow, jcol ) )
533  180 CONTINUE
534  190 CONTINUE
535  CALL dgemm( 'T', 'N', nsize, nrhs, nsize, one,
536  $ rwork( u+st1 ), n, rwork( irwb ), nsize,
537  $ zero, rwork( irwrb ), nsize )
538  j = irwb - 1
539  DO 210 jcol = 1, nrhs
540  DO 200 jrow = st, st + nsize - 1
541  j = j + 1
542  rwork( j ) = dimag( b( jrow, jcol ) )
543  200 CONTINUE
544  210 CONTINUE
545  CALL dgemm( 'T', 'N', nsize, nrhs, nsize, one,
546  $ rwork( u+st1 ), n, rwork( irwb ), nsize,
547  $ zero, rwork( irwib ), nsize )
548  jreal = irwrb - 1
549  jimag = irwib - 1
550  DO 230 jcol = 1, nrhs
551  DO 220 jrow = st, st + nsize - 1
552  jreal = jreal + 1
553  jimag = jimag + 1
554  b( jrow, jcol ) = dcmplx( rwork( jreal ),
555  $ rwork( jimag ) )
556  220 CONTINUE
557  230 CONTINUE
558 *
559  CALL zlacpy( 'A', nsize, nrhs, b( st, 1 ), ldb,
560  $ work( bx+st1 ), n )
561  ELSE
562 *
563 * A large problem. Solve it using divide and conquer.
564 *
565  CALL dlasda( icmpq1, smlsiz, nsize, sqre, d( st ),
566  $ e( st ), rwork( u+st1 ), n, rwork( vt+st1 ),
567  $ iwork( k+st1 ), rwork( difl+st1 ),
568  $ rwork( difr+st1 ), rwork( z+st1 ),
569  $ rwork( poles+st1 ), iwork( givptr+st1 ),
570  $ iwork( givcol+st1 ), n, iwork( perm+st1 ),
571  $ rwork( givnum+st1 ), rwork( c+st1 ),
572  $ rwork( s+st1 ), rwork( nrwork ),
573  $ iwork( iwk ), info )
574  IF( info.NE.0 ) THEN
575  RETURN
576  END IF
577  bxst = bx + st1
578  CALL zlalsa( icmpq2, smlsiz, nsize, nrhs, b( st, 1 ),
579  $ ldb, work( bxst ), n, rwork( u+st1 ), n,
580  $ rwork( vt+st1 ), iwork( k+st1 ),
581  $ rwork( difl+st1 ), rwork( difr+st1 ),
582  $ rwork( z+st1 ), rwork( poles+st1 ),
583  $ iwork( givptr+st1 ), iwork( givcol+st1 ), n,
584  $ iwork( perm+st1 ), rwork( givnum+st1 ),
585  $ rwork( c+st1 ), rwork( s+st1 ),
586  $ rwork( nrwork ), iwork( iwk ), info )
587  IF( info.NE.0 ) THEN
588  RETURN
589  END IF
590  END IF
591  st = i + 1
592  END IF
593  240 CONTINUE
594 *
595 * Apply the singular values and treat the tiny ones as zero.
596 *
597  tol = rcnd*abs( d( idamax( n, d, 1 ) ) )
598 *
599  DO 250 i = 1, n
600 *
601 * Some of the elements in D can be negative because 1-by-1
602 * subproblems were not solved explicitly.
603 *
604  IF( abs( d( i ) ).LE.tol ) THEN
605  CALL zlaset( 'A', 1, nrhs, czero, czero, work( bx+i-1 ), n )
606  ELSE
607  rank = rank + 1
608  CALL zlascl( 'G', 0, 0, d( i ), one, 1, nrhs,
609  $ work( bx+i-1 ), n, info )
610  END IF
611  d( i ) = abs( d( i ) )
612  250 CONTINUE
613 *
614 * Now apply back the right singular vectors.
615 *
616  icmpq2 = 1
617  DO 320 i = 1, nsub
618  st = iwork( i )
619  st1 = st - 1
620  nsize = iwork( sizei+i-1 )
621  bxst = bx + st1
622  IF( nsize.EQ.1 ) THEN
623  CALL zcopy( nrhs, work( bxst ), n, b( st, 1 ), ldb )
624  ELSE IF( nsize.LE.smlsiz ) THEN
625 *
626 * Since B and BX are complex, the following call to DGEMM
627 * is performed in two steps (real and imaginary parts).
628 *
629 * CALL DGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE,
630 * $ RWORK( VT+ST1 ), N, RWORK( BXST ), N, ZERO,
631 * $ B( ST, 1 ), LDB )
632 *
633  j = bxst - n - 1
634  jreal = irwb - 1
635  DO 270 jcol = 1, nrhs
636  j = j + n
637  DO 260 jrow = 1, nsize
638  jreal = jreal + 1
639  rwork( jreal ) = dble( work( j+jrow ) )
640  260 CONTINUE
641  270 CONTINUE
642  CALL dgemm( 'T', 'N', nsize, nrhs, nsize, one,
643  $ rwork( vt+st1 ), n, rwork( irwb ), nsize, zero,
644  $ rwork( irwrb ), nsize )
645  j = bxst - n - 1
646  jimag = irwb - 1
647  DO 290 jcol = 1, nrhs
648  j = j + n
649  DO 280 jrow = 1, nsize
650  jimag = jimag + 1
651  rwork( jimag ) = dimag( work( j+jrow ) )
652  280 CONTINUE
653  290 CONTINUE
654  CALL dgemm( 'T', 'N', nsize, nrhs, nsize, one,
655  $ rwork( vt+st1 ), n, rwork( irwb ), nsize, zero,
656  $ rwork( irwib ), nsize )
657  jreal = irwrb - 1
658  jimag = irwib - 1
659  DO 310 jcol = 1, nrhs
660  DO 300 jrow = st, st + nsize - 1
661  jreal = jreal + 1
662  jimag = jimag + 1
663  b( jrow, jcol ) = dcmplx( rwork( jreal ),
664  $ rwork( jimag ) )
665  300 CONTINUE
666  310 CONTINUE
667  ELSE
668  CALL zlalsa( icmpq2, smlsiz, nsize, nrhs, work( bxst ), n,
669  $ b( st, 1 ), ldb, rwork( u+st1 ), n,
670  $ rwork( vt+st1 ), iwork( k+st1 ),
671  $ rwork( difl+st1 ), rwork( difr+st1 ),
672  $ rwork( z+st1 ), rwork( poles+st1 ),
673  $ iwork( givptr+st1 ), iwork( givcol+st1 ), n,
674  $ iwork( perm+st1 ), rwork( givnum+st1 ),
675  $ rwork( c+st1 ), rwork( s+st1 ),
676  $ rwork( nrwork ), iwork( iwk ), info )
677  IF( info.NE.0 ) THEN
678  RETURN
679  END IF
680  END IF
681  320 CONTINUE
682 *
683 * Unscale and sort the singular values.
684 *
685  CALL dlascl( 'G', 0, 0, one, orgnrm, n, 1, d, n, info )
686  CALL dlasrt( 'D', n, d, info )
687  CALL zlascl( 'G', 0, 0, orgnrm, one, n, nrhs, b, ldb, info )
688 *
689  RETURN
690 *
691 * End of ZLALSD
692 *
693  END
dlascl
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
dlartg
subroutine dlartg(F, G, CS, SN, R)
DLARTG generates a plane rotation with real cosine and real sine.
Definition: dlartg.f:99
zlalsd
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:189
zcopy
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
Definition: zcopy.f:83
zlaset
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
zlascl
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
dlasdq
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
zlacpy
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
dgemm
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
Definition: dgemm.f:189
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
dlaset
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
dlasda
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:275
dlasrt
subroutine dlasrt(ID, N, D, INFO)
DLASRT sorts numbers in increasing or decreasing order.
Definition: dlasrt.f:90
zlalsa
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:269
zdrot
subroutine zdrot(N, CX, INCX, CY, INCY, C, S)
ZDROT
Definition: zdrot.f:100