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