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