LAPACK 3.11.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*> 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:111
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