LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
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*> Download CLALSD + dependencies
9*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/clalsd.f">
10*> [TGZ]</a>
11*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/clalsd.f">
12*> [ZIP]</a>
13*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/clalsd.f">
14*> [TXT]</a>
15*
16* Definition:
17* ===========
18*
19* SUBROUTINE CLALSD( UPLO, SMLSIZ, N, NRHS, D, E, B, LDB, RCOND,
20* RANK, WORK, RWORK, IWORK, INFO )
21*
22* .. Scalar Arguments ..
23* CHARACTER UPLO
24* INTEGER INFO, LDB, N, NRHS, RANK, SMLSIZ
25* REAL RCOND
26* ..
27* .. Array Arguments ..
28* INTEGER IWORK( * )
29* REAL D( * ), E( * ), RWORK( * )
30* COMPLEX B( LDB, * ), WORK( * )
31* ..
32*
33*
34*> \par Purpose:
35* =============
36*>
37*> \verbatim
38*>
39*> CLALSD uses the singular value decomposition of A to solve the least
40*> squares problem of finding X to minimize the Euclidean norm of each
41*> column of A*X-B, where A is N-by-N upper bidiagonal, and X and B
42*> are N-by-NRHS. The solution X overwrites B.
43*>
44*> The singular values of A smaller than RCOND times the largest
45*> singular value are treated as zero in solving the least squares
46*> problem; in this case a minimum norm solution is returned.
47*> The actual singular values are returned in D in ascending order.
48*>
49*> \endverbatim
50*
51* Arguments:
52* ==========
53*
54*> \param[in] UPLO
55*> \verbatim
56*> UPLO is CHARACTER*1
57*> = 'U': D and E define an upper bidiagonal matrix.
58*> = 'L': D and E define a lower bidiagonal matrix.
59*> \endverbatim
60*>
61*> \param[in] SMLSIZ
62*> \verbatim
63*> SMLSIZ is INTEGER
64*> The maximum size of the subproblems at the bottom of the
65*> computation tree.
66*> \endverbatim
67*>
68*> \param[in] N
69*> \verbatim
70*> N is INTEGER
71*> The dimension of the bidiagonal matrix. N >= 0.
72*> \endverbatim
73*>
74*> \param[in] NRHS
75*> \verbatim
76*> NRHS is INTEGER
77*> The number of columns of B. NRHS must be at least 1.
78*> \endverbatim
79*>
80*> \param[in,out] D
81*> \verbatim
82*> D is REAL array, dimension (N)
83*> On entry D contains the main diagonal of the bidiagonal
84*> matrix. On exit, if INFO = 0, D contains its singular values.
85*> \endverbatim
86*>
87*> \param[in,out] E
88*> \verbatim
89*> E is REAL array, dimension (N-1)
90*> Contains the super-diagonal entries of the bidiagonal matrix.
91*> On exit, E has been destroyed.
92*> \endverbatim
93*>
94*> \param[in,out] B
95*> \verbatim
96*> B is COMPLEX array, dimension (LDB,NRHS)
97*> On input, B contains the right hand sides of the least
98*> squares problem. On output, B contains the solution X.
99*> \endverbatim
100*>
101*> \param[in] LDB
102*> \verbatim
103*> LDB is INTEGER
104*> The leading dimension of B in the calling subprogram.
105*> LDB must be at least max(1,N).
106*> \endverbatim
107*>
108*> \param[in] RCOND
109*> \verbatim
110*> RCOND is REAL
111*> The singular values of A less than or equal to RCOND times
112*> the largest singular value are treated as zero in solving
113*> the least squares problem. If RCOND is negative,
114*> machine precision is used instead.
115*> For example, if diag(S)*X=B were the least squares problem,
116*> where diag(S) is a diagonal matrix of singular values, the
117*> solution would be X(i) = B(i) / S(i) if S(i) is greater than
118*> RCOND*max(S), and X(i) = 0 if S(i) is less than or equal to
119*> RCOND*max(S).
120*> \endverbatim
121*>
122*> \param[out] RANK
123*> \verbatim
124*> RANK is INTEGER
125*> The number of singular values of A greater than RCOND times
126*> the largest singular value.
127*> \endverbatim
128*>
129*> \param[out] WORK
130*> \verbatim
131*> WORK is COMPLEX array, dimension (N * NRHS).
132*> \endverbatim
133*>
134*> \param[out] RWORK
135*> \verbatim
136*> RWORK is REAL array, dimension at least
137*> (9*N + 2*N*SMLSIZ + 8*N*NLVL + 3*SMLSIZ*NRHS +
138*> MAX( (SMLSIZ+1)**2, N*(1+NRHS) + 2*NRHS ),
139*> where
140*> NLVL = MAX( 0, INT( LOG_2( MIN( M,N )/(SMLSIZ+1) ) ) + 1 )
141*> \endverbatim
142*>
143*> \param[out] IWORK
144*> \verbatim
145*> IWORK is INTEGER array, dimension (3*N*NLVL + 11*N).
146*> \endverbatim
147*>
148*> \param[out] INFO
149*> \verbatim
150*> INFO is INTEGER
151*> = 0: successful exit.
152*> < 0: if INFO = -i, the i-th argument had an illegal value.
153*> > 0: The algorithm failed to compute a singular value while
154*> working on the submatrix lying in rows and columns
155*> INFO/(N+1) through MOD(INFO,N+1).
156*> \endverbatim
157*
158* Authors:
159* ========
160*
161*> \author Univ. of Tennessee
162*> \author Univ. of California Berkeley
163*> \author Univ. of Colorado Denver
164*> \author NAG Ltd.
165*
166*> \ingroup lalsd
167*
168*> \par Contributors:
169* ==================
170*>
171*> Ming Gu and Ren-Cang Li, Computer Science Division, University of
172*> California at Berkeley, USA \n
173*> Osni Marques, LBNL/NERSC, USA \n
174*
175* =====================================================================
176 SUBROUTINE clalsd( UPLO, SMLSIZ, N, NRHS, D, E, B, LDB, RCOND,
177 $ RANK, WORK, RWORK, IWORK, INFO )
178*
179* -- LAPACK computational routine --
180* -- LAPACK is a software package provided by Univ. of Tennessee, --
181* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
182*
183* .. Scalar Arguments ..
184 CHARACTER UPLO
185 INTEGER INFO, LDB, N, NRHS, RANK, SMLSIZ
186 REAL RCOND
187* ..
188* .. Array Arguments ..
189 INTEGER IWORK( * )
190 REAL D( * ), E( * ), RWORK( * )
191 COMPLEX B( LDB, * ), WORK( * )
192* ..
193*
194* =====================================================================
195*
196* .. Parameters ..
197 REAL ZERO, ONE, TWO
198 parameter( zero = 0.0e0, one = 1.0e0, two = 2.0e0 )
199 COMPLEX CZERO
200 parameter( czero = ( 0.0e0, 0.0e0 ) )
201* ..
202* .. Local Scalars ..
203 INTEGER BX, BXST, C, DIFL, DIFR, GIVCOL, GIVNUM,
204 $ givptr, i, icmpq1, icmpq2, irwb, irwib, irwrb,
205 $ irwu, irwvt, irwwrk, iwk, j, jcol, jimag,
206 $ jreal, jrow, k, nlvl, nm1, nrwork, nsize, nsub,
207 $ perm, poles, s, sizei, smlszp, sqre, st, st1,
208 $ u, vt, z
209 REAL CS, EPS, ORGNRM, R, RCND, SN, TOL
210* ..
211* .. External Functions ..
212 INTEGER ISAMAX
213 REAL SLAMCH, SLANST
214 EXTERNAL isamax, slamch, slanst
215* ..
216* .. External Subroutines ..
217 EXTERNAL ccopy, clacpy, clalsa, clascl, claset,
218 $ csrot,
220 $ slasrt, xerbla
221* ..
222* .. Intrinsic Functions ..
223 INTRINSIC abs, aimag, cmplx, int, log, real, sign
224* ..
225* .. Executable Statements ..
226*
227* Test the input parameters.
228*
229 info = 0
230*
231 IF( n.LT.0 ) THEN
232 info = -3
233 ELSE IF( nrhs.LT.1 ) THEN
234 info = -4
235 ELSE IF( ( ldb.LT.1 ) .OR. ( ldb.LT.n ) ) THEN
236 info = -8
237 END IF
238 IF( info.NE.0 ) THEN
239 CALL xerbla( 'CLALSD', -info )
240 RETURN
241 END IF
242*
243 eps = slamch( 'Epsilon' )
244*
245* Set up the tolerance.
246*
247 IF( ( rcond.LE.zero ) .OR. ( rcond.GE.one ) ) THEN
248 rcnd = eps
249 ELSE
250 rcnd = rcond
251 END IF
252*
253 rank = 0
254*
255* Quick return if possible.
256*
257 IF( n.EQ.0 ) THEN
258 RETURN
259 ELSE IF( n.EQ.1 ) THEN
260 IF( d( 1 ).EQ.zero ) THEN
261 CALL claset( 'A', 1, nrhs, czero, czero, b, ldb )
262 ELSE
263 rank = 1
264 CALL clascl( 'G', 0, 0, d( 1 ), one, 1, nrhs, b, ldb,
265 $ info )
266 d( 1 ) = abs( d( 1 ) )
267 END IF
268 RETURN
269 END IF
270*
271* Rotate the matrix if it is lower bidiagonal.
272*
273 IF( uplo.EQ.'L' ) THEN
274 DO 10 i = 1, n - 1
275 CALL slartg( d( i ), e( i ), cs, sn, r )
276 d( i ) = r
277 e( i ) = sn*d( i+1 )
278 d( i+1 ) = cs*d( i+1 )
279 IF( nrhs.EQ.1 ) THEN
280 CALL csrot( 1, b( i, 1 ), 1, b( i+1, 1 ), 1, cs, sn )
281 ELSE
282 rwork( i*2-1 ) = cs
283 rwork( i*2 ) = sn
284 END IF
285 10 CONTINUE
286 IF( nrhs.GT.1 ) THEN
287 DO 30 i = 1, nrhs
288 DO 20 j = 1, n - 1
289 cs = rwork( j*2-1 )
290 sn = rwork( j*2 )
291 CALL csrot( 1, b( j, i ), 1, b( j+1, i ), 1, cs,
292 $ sn )
293 20 CONTINUE
294 30 CONTINUE
295 END IF
296 END IF
297*
298* Scale.
299*
300 nm1 = n - 1
301 orgnrm = slanst( 'M', n, d, e )
302 IF( orgnrm.EQ.zero ) THEN
303 CALL claset( 'A', n, nrhs, czero, czero, b, ldb )
304 RETURN
305 END IF
306*
307 CALL slascl( 'G', 0, 0, orgnrm, one, n, 1, d, n, info )
308 CALL slascl( '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 slaset( 'A', n, n, zero, one, rwork( irwu ), n )
321 CALL slaset( 'A', n, n, zero, one, rwork( irwvt ), n )
322 CALL slasdq( '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 SLASDQ 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 ) = real( b( jrow, jcol ) )
338 40 CONTINUE
339 50 CONTINUE
340 CALL sgemm( '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 ) = aimag( b( jrow, jcol ) )
347 60 CONTINUE
348 70 CONTINUE
349 CALL sgemm( '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 ) = cmplx( rwork( jreal ), rwork( jimag ) )
358 80 CONTINUE
359 90 CONTINUE
360*
361 tol = rcnd*abs( d( isamax( n, d, 1 ) ) )
362 DO 100 i = 1, n
363 IF( d( i ).LE.tol ) THEN
364 CALL claset( 'A', 1, nrhs, czero, czero, b( i, 1 ),
365 $ ldb )
366 ELSE
367 CALL clascl( 'G', 0, 0, d( i ), one, 1, nrhs, b( i,
368 $ 1 ),
369 $ ldb, info )
370 rank = rank + 1
371 END IF
372 100 CONTINUE
373*
374* Since B is complex, the following call to SGEMM is performed
375* in two steps (real and imaginary parts). That is for V * B
376* (in the real version of the code V**H is stored in WORK).
377*
378* CALL SGEMM( 'T', 'N', N, NRHS, N, ONE, WORK, N, B, LDB, ZERO,
379* $ WORK( NWORK ), N )
380*
381 j = irwb - 1
382 DO 120 jcol = 1, nrhs
383 DO 110 jrow = 1, n
384 j = j + 1
385 rwork( j ) = real( b( jrow, jcol ) )
386 110 CONTINUE
387 120 CONTINUE
388 CALL sgemm( 'T', 'N', n, nrhs, n, one, rwork( irwvt ), n,
389 $ rwork( irwb ), n, zero, rwork( irwrb ), n )
390 j = irwb - 1
391 DO 140 jcol = 1, nrhs
392 DO 130 jrow = 1, n
393 j = j + 1
394 rwork( j ) = aimag( b( jrow, jcol ) )
395 130 CONTINUE
396 140 CONTINUE
397 CALL sgemm( 'T', 'N', n, nrhs, n, one, rwork( irwvt ), n,
398 $ rwork( irwb ), n, zero, rwork( irwib ), n )
399 jreal = irwrb - 1
400 jimag = irwib - 1
401 DO 160 jcol = 1, nrhs
402 DO 150 jrow = 1, n
403 jreal = jreal + 1
404 jimag = jimag + 1
405 b( jrow, jcol ) = cmplx( rwork( jreal ), rwork( jimag ) )
406 150 CONTINUE
407 160 CONTINUE
408*
409* Unscale.
410*
411 CALL slascl( 'G', 0, 0, one, orgnrm, n, 1, d, n, info )
412 CALL slasrt( 'D', n, d, info )
413 CALL clascl( '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( real( n ) / real( 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 ccopy( 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 ccopy( 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 SLASDQ.
502*
503 CALL slaset( 'A', nsize, nsize, zero, one,
504 $ rwork( vt+st1 ), n )
505 CALL slaset( 'A', nsize, nsize, zero, one,
506 $ rwork( u+st1 ), n )
507 CALL slasdq( '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 SLASDQ 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 ) = real( b( jrow, jcol ) )
524 180 CONTINUE
525 190 CONTINUE
526 CALL sgemm( '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 ) = aimag( b( jrow, jcol ) )
534 200 CONTINUE
535 210 CONTINUE
536 CALL sgemm( '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 ) = cmplx( rwork( jreal ),
546 $ rwork( jimag ) )
547 220 CONTINUE
548 230 CONTINUE
549*
550 CALL clacpy( '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 slasda( 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 clalsa( 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( isamax( 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 claset( 'A', 1, nrhs, czero, czero, work( bx+i-1 ),
597 $ n )
598 ELSE
599 rank = rank + 1
600 CALL clascl( 'G', 0, 0, d( i ), one, 1, nrhs,
601 $ work( bx+i-1 ), n, info )
602 END IF
603 d( i ) = abs( d( i ) )
604 250 CONTINUE
605*
606* Now apply back the right singular vectors.
607*
608 icmpq2 = 1
609 DO 320 i = 1, nsub
610 st = iwork( i )
611 st1 = st - 1
612 nsize = iwork( sizei+i-1 )
613 bxst = bx + st1
614 IF( nsize.EQ.1 ) THEN
615 CALL ccopy( nrhs, work( bxst ), n, b( st, 1 ), ldb )
616 ELSE IF( nsize.LE.smlsiz ) THEN
617*
618* Since B and BX are complex, the following call to SGEMM
619* is performed in two steps (real and imaginary parts).
620*
621* CALL SGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE,
622* $ RWORK( VT+ST1 ), N, RWORK( BXST ), N, ZERO,
623* $ B( ST, 1 ), LDB )
624*
625 j = bxst - n - 1
626 jreal = irwb - 1
627 DO 270 jcol = 1, nrhs
628 j = j + n
629 DO 260 jrow = 1, nsize
630 jreal = jreal + 1
631 rwork( jreal ) = real( work( j+jrow ) )
632 260 CONTINUE
633 270 CONTINUE
634 CALL sgemm( 'T', 'N', nsize, nrhs, nsize, one,
635 $ rwork( vt+st1 ), n, rwork( irwb ), nsize, zero,
636 $ rwork( irwrb ), nsize )
637 j = bxst - n - 1
638 jimag = irwb - 1
639 DO 290 jcol = 1, nrhs
640 j = j + n
641 DO 280 jrow = 1, nsize
642 jimag = jimag + 1
643 rwork( jimag ) = aimag( work( j+jrow ) )
644 280 CONTINUE
645 290 CONTINUE
646 CALL sgemm( 'T', 'N', nsize, nrhs, nsize, one,
647 $ rwork( vt+st1 ), n, rwork( irwb ), nsize, zero,
648 $ rwork( irwib ), nsize )
649 jreal = irwrb - 1
650 jimag = irwib - 1
651 DO 310 jcol = 1, nrhs
652 DO 300 jrow = st, st + nsize - 1
653 jreal = jreal + 1
654 jimag = jimag + 1
655 b( jrow, jcol ) = cmplx( rwork( jreal ),
656 $ rwork( jimag ) )
657 300 CONTINUE
658 310 CONTINUE
659 ELSE
660 CALL clalsa( icmpq2, smlsiz, nsize, nrhs, work( bxst ),
661 $ n,
662 $ b( st, 1 ), ldb, rwork( u+st1 ), n,
663 $ rwork( vt+st1 ), iwork( k+st1 ),
664 $ rwork( difl+st1 ), rwork( difr+st1 ),
665 $ rwork( z+st1 ), rwork( poles+st1 ),
666 $ iwork( givptr+st1 ), iwork( givcol+st1 ), n,
667 $ iwork( perm+st1 ), rwork( givnum+st1 ),
668 $ rwork( c+st1 ), rwork( s+st1 ),
669 $ rwork( nrwork ), iwork( iwk ), info )
670 IF( info.NE.0 ) THEN
671 RETURN
672 END IF
673 END IF
674 320 CONTINUE
675*
676* Unscale and sort the singular values.
677*
678 CALL slascl( 'G', 0, 0, one, orgnrm, n, 1, d, n, info )
679 CALL slasrt( 'D', n, d, info )
680 CALL clascl( 'G', 0, 0, orgnrm, one, n, nrhs, b, ldb, info )
681*
682 RETURN
683*
684* End of CLALSD
685*
686 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine ccopy(n, cx, incx, cy, incy)
CCOPY
Definition ccopy.f:81
subroutine sgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
SGEMM
Definition sgemm.f:188
subroutine clacpy(uplo, m, n, a, lda, b, ldb)
CLACPY copies all or part of one two-dimensional array to another.
Definition clacpy.f:101
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:266
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:178
subroutine slartg(f, g, c, s, r)
SLARTG generates a plane rotation with real cosine and real sine.
Definition slartg.f90:111
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:142
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:142
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:272
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:210
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:108
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:104
subroutine slasrt(id, n, d, info)
SLASRT sorts numbers in increasing or decreasing order.
Definition slasrt.f:86
subroutine csrot(n, cx, incx, cy, incy, c, s)
CSROT
Definition csrot.f:98