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