LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
clalsa.f
Go to the documentation of this file.
1 *> \brief \b CLALSA computes the SVD of the coefficient matrix in compact form. Used by sgelsd.
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download CLALSA + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/clalsa.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/clalsa.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/clalsa.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE CLALSA( ICOMPQ, SMLSIZ, N, NRHS, B, LDB, BX, LDBX, U,
22 * LDU, VT, K, DIFL, DIFR, Z, POLES, GIVPTR,
23 * GIVCOL, LDGCOL, PERM, GIVNUM, C, S, RWORK,
24 * IWORK, INFO )
25 *
26 * .. Scalar Arguments ..
27 * INTEGER ICOMPQ, INFO, LDB, LDBX, LDGCOL, LDU, N, NRHS,
28 * $ SMLSIZ
29 * ..
30 * .. Array Arguments ..
31 * INTEGER GIVCOL( LDGCOL, * ), GIVPTR( * ), IWORK( * ),
32 * $ K( * ), PERM( LDGCOL, * )
33 * REAL C( * ), DIFL( LDU, * ), DIFR( LDU, * ),
34 * $ GIVNUM( LDU, * ), POLES( LDU, * ), RWORK( * ),
35 * $ S( * ), U( LDU, * ), VT( LDU, * ), Z( LDU, * )
36 * COMPLEX B( LDB, * ), BX( LDBX, * )
37 * ..
38 *
39 *
40 *> \par Purpose:
41 * =============
42 *>
43 *> \verbatim
44 *>
45 *> CLALSA is an itermediate step in solving the least squares problem
46 *> by computing the SVD of the coefficient matrix in compact form (The
47 *> singular vectors are computed as products of simple orthorgonal
48 *> matrices.).
49 *>
50 *> If ICOMPQ = 0, CLALSA applies the inverse of the left singular vector
51 *> matrix of an upper bidiagonal matrix to the right hand side; and if
52 *> ICOMPQ = 1, CLALSA applies the right singular vector matrix to the
53 *> right hand side. The singular vector matrices were generated in
54 *> compact form by CLALSA.
55 *> \endverbatim
56 *
57 * Arguments:
58 * ==========
59 *
60 *> \param[in] ICOMPQ
61 *> \verbatim
62 *> ICOMPQ is INTEGER
63 *> Specifies whether the left or the right singular vector
64 *> matrix is involved.
65 *> = 0: Left singular vector matrix
66 *> = 1: Right singular vector 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 row and column dimensions of the upper bidiagonal matrix.
80 *> \endverbatim
81 *>
82 *> \param[in] NRHS
83 *> \verbatim
84 *> NRHS is INTEGER
85 *> The number of columns of B and BX. NRHS must be at least 1.
86 *> \endverbatim
87 *>
88 *> \param[in,out] B
89 *> \verbatim
90 *> B is COMPLEX array, dimension ( LDB, NRHS )
91 *> On input, B contains the right hand sides of the least
92 *> squares problem in rows 1 through M.
93 *> On output, B contains the solution X in rows 1 through N.
94 *> \endverbatim
95 *>
96 *> \param[in] LDB
97 *> \verbatim
98 *> LDB is INTEGER
99 *> The leading dimension of B in the calling subprogram.
100 *> LDB must be at least max(1,MAX( M, N ) ).
101 *> \endverbatim
102 *>
103 *> \param[out] BX
104 *> \verbatim
105 *> BX is COMPLEX array, dimension ( LDBX, NRHS )
106 *> On exit, the result of applying the left or right singular
107 *> vector matrix to B.
108 *> \endverbatim
109 *>
110 *> \param[in] LDBX
111 *> \verbatim
112 *> LDBX is INTEGER
113 *> The leading dimension of BX.
114 *> \endverbatim
115 *>
116 *> \param[in] U
117 *> \verbatim
118 *> U is REAL array, dimension ( LDU, SMLSIZ ).
119 *> On entry, U contains the left singular vector matrices of all
120 *> subproblems at the bottom level.
121 *> \endverbatim
122 *>
123 *> \param[in] LDU
124 *> \verbatim
125 *> LDU is INTEGER, LDU = > N.
126 *> The leading dimension of arrays U, VT, DIFL, DIFR,
127 *> POLES, GIVNUM, and Z.
128 *> \endverbatim
129 *>
130 *> \param[in] VT
131 *> \verbatim
132 *> VT is REAL array, dimension ( LDU, SMLSIZ+1 ).
133 *> On entry, VT**H contains the right singular vector matrices of
134 *> all subproblems at the bottom level.
135 *> \endverbatim
136 *>
137 *> \param[in] K
138 *> \verbatim
139 *> K is INTEGER array, dimension ( N ).
140 *> \endverbatim
141 *>
142 *> \param[in] DIFL
143 *> \verbatim
144 *> DIFL is REAL array, dimension ( LDU, NLVL ).
145 *> where NLVL = INT(log_2 (N/(SMLSIZ+1))) + 1.
146 *> \endverbatim
147 *>
148 *> \param[in] DIFR
149 *> \verbatim
150 *> DIFR is REAL array, dimension ( LDU, 2 * NLVL ).
151 *> On entry, DIFL(*, I) and DIFR(*, 2 * I -1) record
152 *> distances between singular values on the I-th level and
153 *> singular values on the (I -1)-th level, and DIFR(*, 2 * I)
154 *> record the normalizing factors of the right singular vectors
155 *> matrices of subproblems on I-th level.
156 *> \endverbatim
157 *>
158 *> \param[in] Z
159 *> \verbatim
160 *> Z is REAL array, dimension ( LDU, NLVL ).
161 *> On entry, Z(1, I) contains the components of the deflation-
162 *> adjusted updating row vector for subproblems on the I-th
163 *> level.
164 *> \endverbatim
165 *>
166 *> \param[in] POLES
167 *> \verbatim
168 *> POLES is REAL array, dimension ( LDU, 2 * NLVL ).
169 *> On entry, POLES(*, 2 * I -1: 2 * I) contains the new and old
170 *> singular values involved in the secular equations on the I-th
171 *> level.
172 *> \endverbatim
173 *>
174 *> \param[in] GIVPTR
175 *> \verbatim
176 *> GIVPTR is INTEGER array, dimension ( N ).
177 *> On entry, GIVPTR( I ) records the number of Givens
178 *> rotations performed on the I-th problem on the computation
179 *> tree.
180 *> \endverbatim
181 *>
182 *> \param[in] GIVCOL
183 *> \verbatim
184 *> GIVCOL is INTEGER array, dimension ( LDGCOL, 2 * NLVL ).
185 *> On entry, for each I, GIVCOL(*, 2 * I - 1: 2 * I) records the
186 *> locations of Givens rotations performed on the I-th level on
187 *> the computation tree.
188 *> \endverbatim
189 *>
190 *> \param[in] LDGCOL
191 *> \verbatim
192 *> LDGCOL is INTEGER, LDGCOL = > N.
193 *> The leading dimension of arrays GIVCOL and PERM.
194 *> \endverbatim
195 *>
196 *> \param[in] PERM
197 *> \verbatim
198 *> PERM is INTEGER array, dimension ( LDGCOL, NLVL ).
199 *> On entry, PERM(*, I) records permutations done on the I-th
200 *> level of the computation tree.
201 *> \endverbatim
202 *>
203 *> \param[in] GIVNUM
204 *> \verbatim
205 *> GIVNUM is REAL array, dimension ( LDU, 2 * NLVL ).
206 *> On entry, GIVNUM(*, 2 *I -1 : 2 * I) records the C- and S-
207 *> values of Givens rotations performed on the I-th level on the
208 *> computation tree.
209 *> \endverbatim
210 *>
211 *> \param[in] C
212 *> \verbatim
213 *> C is REAL array, dimension ( N ).
214 *> On entry, if the I-th subproblem is not square,
215 *> C( I ) contains the C-value of a Givens rotation related to
216 *> the right null space of the I-th subproblem.
217 *> \endverbatim
218 *>
219 *> \param[in] S
220 *> \verbatim
221 *> S is REAL array, dimension ( N ).
222 *> On entry, if the I-th subproblem is not square,
223 *> S( I ) contains the S-value of a Givens rotation related to
224 *> the right null space of the I-th subproblem.
225 *> \endverbatim
226 *>
227 *> \param[out] RWORK
228 *> \verbatim
229 *> RWORK is REAL array, dimension at least
230 *> MAX( (SMLSZ+1)*NRHS*3, N*(1+NRHS) + 2*NRHS ).
231 *> \endverbatim
232 *>
233 *> \param[out] IWORK
234 *> \verbatim
235 *> IWORK is INTEGER array.
236 *> The dimension must be at least 3 * N
237 *> \endverbatim
238 *>
239 *> \param[out] INFO
240 *> \verbatim
241 *> INFO is INTEGER
242 *> = 0: successful exit.
243 *> < 0: if INFO = -i, the i-th argument had an illegal value.
244 *> \endverbatim
245 *
246 * Authors:
247 * ========
248 *
249 *> \author Univ. of Tennessee
250 *> \author Univ. of California Berkeley
251 *> \author Univ. of Colorado Denver
252 *> \author NAG Ltd.
253 *
254 *> \date September 2012
255 *
256 *> \ingroup complexOTHERcomputational
257 *
258 *> \par Contributors:
259 * ==================
260 *>
261 *> Ming Gu and Ren-Cang Li, Computer Science Division, University of
262 *> California at Berkeley, USA \n
263 *> Osni Marques, LBNL/NERSC, USA \n
264 *
265 * =====================================================================
266  SUBROUTINE clalsa( ICOMPQ, SMLSIZ, N, NRHS, B, LDB, BX, LDBX, U,
267  $ ldu, vt, k, difl, difr, z, poles, givptr,
268  $ givcol, ldgcol, perm, givnum, c, s, rwork,
269  $ iwork, info )
270 *
271 * -- LAPACK computational routine (version 3.4.2) --
272 * -- LAPACK is a software package provided by Univ. of Tennessee, --
273 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
274 * September 2012
275 *
276 * .. Scalar Arguments ..
277  INTEGER icompq, info, ldb, ldbx, ldgcol, ldu, n, nrhs,
278  $ smlsiz
279 * ..
280 * .. Array Arguments ..
281  INTEGER givcol( ldgcol, * ), givptr( * ), iwork( * ),
282  $ k( * ), perm( ldgcol, * )
283  REAL c( * ), difl( ldu, * ), difr( ldu, * ),
284  $ givnum( ldu, * ), poles( ldu, * ), rwork( * ),
285  $ s( * ), u( ldu, * ), vt( ldu, * ), z( ldu, * )
286  COMPLEX b( ldb, * ), bx( ldbx, * )
287 * ..
288 *
289 * =====================================================================
290 *
291 * .. Parameters ..
292  REAL zero, one
293  parameter( zero = 0.0e0, one = 1.0e0 )
294 * ..
295 * .. Local Scalars ..
296  INTEGER i, i1, ic, im1, inode, j, jcol, jimag, jreal,
297  $ jrow, lf, ll, lvl, lvl2, nd, ndb1, ndiml,
298  $ ndimr, nl, nlf, nlp1, nlvl, nr, nrf, nrp1, sqre
299 * ..
300 * .. External Subroutines ..
301  EXTERNAL ccopy, clals0, sgemm, slasdt, xerbla
302 * ..
303 * .. Intrinsic Functions ..
304  INTRINSIC aimag, cmplx, real
305 * ..
306 * .. Executable Statements ..
307 *
308 * Test the input parameters.
309 *
310  info = 0
311 *
312  IF( ( icompq.LT.0 ) .OR. ( icompq.GT.1 ) ) THEN
313  info = -1
314  ELSE IF( smlsiz.LT.3 ) THEN
315  info = -2
316  ELSE IF( n.LT.smlsiz ) THEN
317  info = -3
318  ELSE IF( nrhs.LT.1 ) THEN
319  info = -4
320  ELSE IF( ldb.LT.n ) THEN
321  info = -6
322  ELSE IF( ldbx.LT.n ) THEN
323  info = -8
324  ELSE IF( ldu.LT.n ) THEN
325  info = -10
326  ELSE IF( ldgcol.LT.n ) THEN
327  info = -19
328  END IF
329  IF( info.NE.0 ) THEN
330  CALL xerbla( 'CLALSA', -info )
331  return
332  END IF
333 *
334 * Book-keeping and setting up the computation tree.
335 *
336  inode = 1
337  ndiml = inode + n
338  ndimr = ndiml + n
339 *
340  CALL slasdt( n, nlvl, nd, iwork( inode ), iwork( ndiml ),
341  $ iwork( ndimr ), smlsiz )
342 *
343 * The following code applies back the left singular vector factors.
344 * For applying back the right singular vector factors, go to 170.
345 *
346  IF( icompq.EQ.1 ) THEN
347  go to 170
348  END IF
349 *
350 * The nodes on the bottom level of the tree were solved
351 * by SLASDQ. The corresponding left and right singular vector
352 * matrices are in explicit form. First apply back the left
353 * singular vector matrices.
354 *
355  ndb1 = ( nd+1 ) / 2
356  DO 130 i = ndb1, nd
357 *
358 * IC : center row of each node
359 * NL : number of rows of left subproblem
360 * NR : number of rows of right subproblem
361 * NLF: starting row of the left subproblem
362 * NRF: starting row of the right subproblem
363 *
364  i1 = i - 1
365  ic = iwork( inode+i1 )
366  nl = iwork( ndiml+i1 )
367  nr = iwork( ndimr+i1 )
368  nlf = ic - nl
369  nrf = ic + 1
370 *
371 * Since B and BX are complex, the following call to SGEMM
372 * is performed in two steps (real and imaginary parts).
373 *
374 * CALL SGEMM( 'T', 'N', NL, NRHS, NL, ONE, U( NLF, 1 ), LDU,
375 * $ B( NLF, 1 ), LDB, ZERO, BX( NLF, 1 ), LDBX )
376 *
377  j = nl*nrhs*2
378  DO 20 jcol = 1, nrhs
379  DO 10 jrow = nlf, nlf + nl - 1
380  j = j + 1
381  rwork( j ) = REAL( B( JROW, JCOL ) )
382  10 continue
383  20 continue
384  CALL sgemm( 'T', 'N', nl, nrhs, nl, one, u( nlf, 1 ), ldu,
385  $ rwork( 1+nl*nrhs*2 ), nl, zero, rwork( 1 ), nl )
386  j = nl*nrhs*2
387  DO 40 jcol = 1, nrhs
388  DO 30 jrow = nlf, nlf + nl - 1
389  j = j + 1
390  rwork( j ) = aimag( b( jrow, jcol ) )
391  30 continue
392  40 continue
393  CALL sgemm( 'T', 'N', nl, nrhs, nl, one, u( nlf, 1 ), ldu,
394  $ rwork( 1+nl*nrhs*2 ), nl, zero, rwork( 1+nl*nrhs ),
395  $ nl )
396  jreal = 0
397  jimag = nl*nrhs
398  DO 60 jcol = 1, nrhs
399  DO 50 jrow = nlf, nlf + nl - 1
400  jreal = jreal + 1
401  jimag = jimag + 1
402  bx( jrow, jcol ) = cmplx( rwork( jreal ),
403  $ rwork( jimag ) )
404  50 continue
405  60 continue
406 *
407 * Since B and BX are complex, the following call to SGEMM
408 * is performed in two steps (real and imaginary parts).
409 *
410 * CALL SGEMM( 'T', 'N', NR, NRHS, NR, ONE, U( NRF, 1 ), LDU,
411 * $ B( NRF, 1 ), LDB, ZERO, BX( NRF, 1 ), LDBX )
412 *
413  j = nr*nrhs*2
414  DO 80 jcol = 1, nrhs
415  DO 70 jrow = nrf, nrf + nr - 1
416  j = j + 1
417  rwork( j ) = REAL( B( JROW, JCOL ) )
418  70 continue
419  80 continue
420  CALL sgemm( 'T', 'N', nr, nrhs, nr, one, u( nrf, 1 ), ldu,
421  $ rwork( 1+nr*nrhs*2 ), nr, zero, rwork( 1 ), nr )
422  j = nr*nrhs*2
423  DO 100 jcol = 1, nrhs
424  DO 90 jrow = nrf, nrf + nr - 1
425  j = j + 1
426  rwork( j ) = aimag( b( jrow, jcol ) )
427  90 continue
428  100 continue
429  CALL sgemm( 'T', 'N', nr, nrhs, nr, one, u( nrf, 1 ), ldu,
430  $ rwork( 1+nr*nrhs*2 ), nr, zero, rwork( 1+nr*nrhs ),
431  $ nr )
432  jreal = 0
433  jimag = nr*nrhs
434  DO 120 jcol = 1, nrhs
435  DO 110 jrow = nrf, nrf + nr - 1
436  jreal = jreal + 1
437  jimag = jimag + 1
438  bx( jrow, jcol ) = cmplx( rwork( jreal ),
439  $ rwork( jimag ) )
440  110 continue
441  120 continue
442 *
443  130 continue
444 *
445 * Next copy the rows of B that correspond to unchanged rows
446 * in the bidiagonal matrix to BX.
447 *
448  DO 140 i = 1, nd
449  ic = iwork( inode+i-1 )
450  CALL ccopy( nrhs, b( ic, 1 ), ldb, bx( ic, 1 ), ldbx )
451  140 continue
452 *
453 * Finally go through the left singular vector matrices of all
454 * the other subproblems bottom-up on the tree.
455 *
456  j = 2**nlvl
457  sqre = 0
458 *
459  DO 160 lvl = nlvl, 1, -1
460  lvl2 = 2*lvl - 1
461 *
462 * find the first node LF and last node LL on
463 * the current level LVL
464 *
465  IF( lvl.EQ.1 ) THEN
466  lf = 1
467  ll = 1
468  ELSE
469  lf = 2**( lvl-1 )
470  ll = 2*lf - 1
471  END IF
472  DO 150 i = lf, ll
473  im1 = i - 1
474  ic = iwork( inode+im1 )
475  nl = iwork( ndiml+im1 )
476  nr = iwork( ndimr+im1 )
477  nlf = ic - nl
478  nrf = ic + 1
479  j = j - 1
480  CALL clals0( icompq, nl, nr, sqre, nrhs, bx( nlf, 1 ), ldbx,
481  $ b( nlf, 1 ), ldb, perm( nlf, lvl ),
482  $ givptr( j ), givcol( nlf, lvl2 ), ldgcol,
483  $ givnum( nlf, lvl2 ), ldu, poles( nlf, lvl2 ),
484  $ difl( nlf, lvl ), difr( nlf, lvl2 ),
485  $ z( nlf, lvl ), k( j ), c( j ), s( j ), rwork,
486  $ info )
487  150 continue
488  160 continue
489  go to 330
490 *
491 * ICOMPQ = 1: applying back the right singular vector factors.
492 *
493  170 continue
494 *
495 * First now go through the right singular vector matrices of all
496 * the tree nodes top-down.
497 *
498  j = 0
499  DO 190 lvl = 1, nlvl
500  lvl2 = 2*lvl - 1
501 *
502 * Find the first node LF and last node LL on
503 * the current level LVL.
504 *
505  IF( lvl.EQ.1 ) THEN
506  lf = 1
507  ll = 1
508  ELSE
509  lf = 2**( lvl-1 )
510  ll = 2*lf - 1
511  END IF
512  DO 180 i = ll, lf, -1
513  im1 = i - 1
514  ic = iwork( inode+im1 )
515  nl = iwork( ndiml+im1 )
516  nr = iwork( ndimr+im1 )
517  nlf = ic - nl
518  nrf = ic + 1
519  IF( i.EQ.ll ) THEN
520  sqre = 0
521  ELSE
522  sqre = 1
523  END IF
524  j = j + 1
525  CALL clals0( icompq, nl, nr, sqre, nrhs, b( nlf, 1 ), ldb,
526  $ bx( nlf, 1 ), ldbx, perm( nlf, lvl ),
527  $ givptr( j ), givcol( nlf, lvl2 ), ldgcol,
528  $ givnum( nlf, lvl2 ), ldu, poles( nlf, lvl2 ),
529  $ difl( nlf, lvl ), difr( nlf, lvl2 ),
530  $ z( nlf, lvl ), k( j ), c( j ), s( j ), rwork,
531  $ info )
532  180 continue
533  190 continue
534 *
535 * The nodes on the bottom level of the tree were solved
536 * by SLASDQ. The corresponding right singular vector
537 * matrices are in explicit form. Apply them back.
538 *
539  ndb1 = ( nd+1 ) / 2
540  DO 320 i = ndb1, nd
541  i1 = i - 1
542  ic = iwork( inode+i1 )
543  nl = iwork( ndiml+i1 )
544  nr = iwork( ndimr+i1 )
545  nlp1 = nl + 1
546  IF( i.EQ.nd ) THEN
547  nrp1 = nr
548  ELSE
549  nrp1 = nr + 1
550  END IF
551  nlf = ic - nl
552  nrf = ic + 1
553 *
554 * Since B and BX are complex, the following call to SGEMM is
555 * performed in two steps (real and imaginary parts).
556 *
557 * CALL SGEMM( 'T', 'N', NLP1, NRHS, NLP1, ONE, VT( NLF, 1 ), LDU,
558 * $ B( NLF, 1 ), LDB, ZERO, BX( NLF, 1 ), LDBX )
559 *
560  j = nlp1*nrhs*2
561  DO 210 jcol = 1, nrhs
562  DO 200 jrow = nlf, nlf + nlp1 - 1
563  j = j + 1
564  rwork( j ) = REAL( B( JROW, JCOL ) )
565  200 continue
566  210 continue
567  CALL sgemm( 'T', 'N', nlp1, nrhs, nlp1, one, vt( nlf, 1 ), ldu,
568  $ rwork( 1+nlp1*nrhs*2 ), nlp1, zero, rwork( 1 ),
569  $ nlp1 )
570  j = nlp1*nrhs*2
571  DO 230 jcol = 1, nrhs
572  DO 220 jrow = nlf, nlf + nlp1 - 1
573  j = j + 1
574  rwork( j ) = aimag( b( jrow, jcol ) )
575  220 continue
576  230 continue
577  CALL sgemm( 'T', 'N', nlp1, nrhs, nlp1, one, vt( nlf, 1 ), ldu,
578  $ rwork( 1+nlp1*nrhs*2 ), nlp1, zero,
579  $ rwork( 1+nlp1*nrhs ), nlp1 )
580  jreal = 0
581  jimag = nlp1*nrhs
582  DO 250 jcol = 1, nrhs
583  DO 240 jrow = nlf, nlf + nlp1 - 1
584  jreal = jreal + 1
585  jimag = jimag + 1
586  bx( jrow, jcol ) = cmplx( rwork( jreal ),
587  $ rwork( jimag ) )
588  240 continue
589  250 continue
590 *
591 * Since B and BX are complex, the following call to SGEMM is
592 * performed in two steps (real and imaginary parts).
593 *
594 * CALL SGEMM( 'T', 'N', NRP1, NRHS, NRP1, ONE, VT( NRF, 1 ), LDU,
595 * $ B( NRF, 1 ), LDB, ZERO, BX( NRF, 1 ), LDBX )
596 *
597  j = nrp1*nrhs*2
598  DO 270 jcol = 1, nrhs
599  DO 260 jrow = nrf, nrf + nrp1 - 1
600  j = j + 1
601  rwork( j ) = REAL( B( JROW, JCOL ) )
602  260 continue
603  270 continue
604  CALL sgemm( 'T', 'N', nrp1, nrhs, nrp1, one, vt( nrf, 1 ), ldu,
605  $ rwork( 1+nrp1*nrhs*2 ), nrp1, zero, rwork( 1 ),
606  $ nrp1 )
607  j = nrp1*nrhs*2
608  DO 290 jcol = 1, nrhs
609  DO 280 jrow = nrf, nrf + nrp1 - 1
610  j = j + 1
611  rwork( j ) = aimag( b( jrow, jcol ) )
612  280 continue
613  290 continue
614  CALL sgemm( 'T', 'N', nrp1, nrhs, nrp1, one, vt( nrf, 1 ), ldu,
615  $ rwork( 1+nrp1*nrhs*2 ), nrp1, zero,
616  $ rwork( 1+nrp1*nrhs ), nrp1 )
617  jreal = 0
618  jimag = nrp1*nrhs
619  DO 310 jcol = 1, nrhs
620  DO 300 jrow = nrf, nrf + nrp1 - 1
621  jreal = jreal + 1
622  jimag = jimag + 1
623  bx( jrow, jcol ) = cmplx( rwork( jreal ),
624  $ rwork( jimag ) )
625  300 continue
626  310 continue
627 *
628  320 continue
629 *
630  330 continue
631 *
632  return
633 *
634 * End of CLALSA
635 *
636  END