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