1: *> \brief \b ZLALSA 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 ZLALSA + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlalsa.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlalsa.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlalsa.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE ZLALSA( 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: * DOUBLE PRECISION C( * ), DIFL( LDU, * ), DIFR( LDU, * ),
34: * $ GIVNUM( LDU, * ), POLES( LDU, * ), RWORK( * ),
35: * $ S( * ), U( LDU, * ), VT( LDU, * ), Z( LDU, * )
36: * COMPLEX*16 B( LDB, * ), BX( LDBX, * )
37: * ..
38: *
39: *
40: *> \par Purpose:
41: * =============
42: *>
43: *> \verbatim
44: *>
45: *> ZLALSA 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, ZLALSA 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, ZLALSA applies the right singular vector matrix to the
53: *> right hand side. The singular vector matrices were generated in
54: *> compact form by ZLALSA.
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*16 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*16 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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, dimension (3*N)
236: *> \endverbatim
237: *>
238: *> \param[out] INFO
239: *> \verbatim
240: *> INFO is INTEGER
241: *> = 0: successful exit.
242: *> < 0: if INFO = -i, the i-th argument had an illegal value.
243: *> \endverbatim
244: *
245: * Authors:
246: * ========
247: *
248: *> \author Univ. of Tennessee
249: *> \author Univ. of California Berkeley
250: *> \author Univ. of Colorado Denver
251: *> \author NAG Ltd.
252: *
253: *> \date June 2017
254: *
255: *> \ingroup complex16OTHERcomputational
256: *
257: *> \par Contributors:
258: * ==================
259: *>
260: *> Ming Gu and Ren-Cang Li, Computer Science Division, University of
261: *> California at Berkeley, USA \n
262: *> Osni Marques, LBNL/NERSC, USA \n
263: *
264: * =====================================================================
265: SUBROUTINE ZLALSA( ICOMPQ, SMLSIZ, N, NRHS, B, LDB, BX, LDBX, U,
266: $ LDU, VT, K, DIFL, DIFR, Z, POLES, GIVPTR,
267: $ GIVCOL, LDGCOL, PERM, GIVNUM, C, S, RWORK,
268: $ IWORK, INFO )
269: *
270: * -- LAPACK computational routine (version 3.7.1) --
271: * -- LAPACK is a software package provided by Univ. of Tennessee, --
272: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
273: * June 2017
274: *
275: * .. Scalar Arguments ..
276: INTEGER ICOMPQ, INFO, LDB, LDBX, LDGCOL, LDU, N, NRHS,
277: $ SMLSIZ
278: * ..
279: * .. Array Arguments ..
280: INTEGER GIVCOL( LDGCOL, * ), GIVPTR( * ), IWORK( * ),
281: $ K( * ), PERM( LDGCOL, * )
282: DOUBLE PRECISION C( * ), DIFL( LDU, * ), DIFR( LDU, * ),
283: $ GIVNUM( LDU, * ), POLES( LDU, * ), RWORK( * ),
284: $ S( * ), U( LDU, * ), VT( LDU, * ), Z( LDU, * )
285: COMPLEX*16 B( LDB, * ), BX( LDBX, * )
286: * ..
287: *
288: * =====================================================================
289: *
290: * .. Parameters ..
291: DOUBLE PRECISION ZERO, ONE
292: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
293: * ..
294: * .. Local Scalars ..
295: INTEGER I, I1, IC, IM1, INODE, J, JCOL, JIMAG, JREAL,
296: $ JROW, LF, LL, LVL, LVL2, ND, NDB1, NDIML,
297: $ NDIMR, NL, NLF, NLP1, NLVL, NR, NRF, NRP1, SQRE
298: * ..
299: * .. External Subroutines ..
300: EXTERNAL DGEMM, DLASDT, XERBLA, ZCOPY, ZLALS0
301: * ..
302: * .. Intrinsic Functions ..
303: INTRINSIC DBLE, DCMPLX, DIMAG
304: * ..
305: * .. Executable Statements ..
306: *
307: * Test the input parameters.
308: *
309: INFO = 0
310: *
311: IF( ( ICOMPQ.LT.0 ) .OR. ( ICOMPQ.GT.1 ) ) THEN
312: INFO = -1
313: ELSE IF( SMLSIZ.LT.3 ) THEN
314: INFO = -2
315: ELSE IF( N.LT.SMLSIZ ) THEN
316: INFO = -3
317: ELSE IF( NRHS.LT.1 ) THEN
318: INFO = -4
319: ELSE IF( LDB.LT.N ) THEN
320: INFO = -6
321: ELSE IF( LDBX.LT.N ) THEN
322: INFO = -8
323: ELSE IF( LDU.LT.N ) THEN
324: INFO = -10
325: ELSE IF( LDGCOL.LT.N ) THEN
326: INFO = -19
327: END IF
328: IF( INFO.NE.0 ) THEN
329: CALL XERBLA( 'ZLALSA', -INFO )
330: RETURN
331: END IF
332: *
333: * Book-keeping and setting up the computation tree.
334: *
335: INODE = 1
336: NDIML = INODE + N
337: NDIMR = NDIML + N
338: *
339: CALL DLASDT( N, NLVL, ND, IWORK( INODE ), IWORK( NDIML ),
340: $ IWORK( NDIMR ), SMLSIZ )
341: *
342: * The following code applies back the left singular vector factors.
343: * For applying back the right singular vector factors, go to 170.
344: *
345: IF( ICOMPQ.EQ.1 ) THEN
346: GO TO 170
347: END IF
348: *
349: * The nodes on the bottom level of the tree were solved
350: * by DLASDQ. The corresponding left and right singular vector
351: * matrices are in explicit form. First apply back the left
352: * singular vector matrices.
353: *
354: NDB1 = ( ND+1 ) / 2
355: DO 130 I = NDB1, ND
356: *
357: * IC : center row of each node
358: * NL : number of rows of left subproblem
359: * NR : number of rows of right subproblem
360: * NLF: starting row of the left subproblem
361: * NRF: starting row of the right subproblem
362: *
363: I1 = I - 1
364: IC = IWORK( INODE+I1 )
365: NL = IWORK( NDIML+I1 )
366: NR = IWORK( NDIMR+I1 )
367: NLF = IC - NL
368: NRF = IC + 1
369: *
370: * Since B and BX are complex, the following call to DGEMM
371: * is performed in two steps (real and imaginary parts).
372: *
373: * CALL DGEMM( 'T', 'N', NL, NRHS, NL, ONE, U( NLF, 1 ), LDU,
374: * $ B( NLF, 1 ), LDB, ZERO, BX( NLF, 1 ), LDBX )
375: *
376: J = NL*NRHS*2
377: DO 20 JCOL = 1, NRHS
378: DO 10 JROW = NLF, NLF + NL - 1
379: J = J + 1
380: RWORK( J ) = DBLE( B( JROW, JCOL ) )
381: 10 CONTINUE
382: 20 CONTINUE
383: CALL DGEMM( 'T', 'N', NL, NRHS, NL, ONE, U( NLF, 1 ), LDU,
384: $ RWORK( 1+NL*NRHS*2 ), NL, ZERO, RWORK( 1 ), NL )
385: J = NL*NRHS*2
386: DO 40 JCOL = 1, NRHS
387: DO 30 JROW = NLF, NLF + NL - 1
388: J = J + 1
389: RWORK( J ) = DIMAG( B( JROW, JCOL ) )
390: 30 CONTINUE
391: 40 CONTINUE
392: CALL DGEMM( 'T', 'N', NL, NRHS, NL, ONE, U( NLF, 1 ), LDU,
393: $ RWORK( 1+NL*NRHS*2 ), NL, ZERO, RWORK( 1+NL*NRHS ),
394: $ NL )
395: JREAL = 0
396: JIMAG = NL*NRHS
397: DO 60 JCOL = 1, NRHS
398: DO 50 JROW = NLF, NLF + NL - 1
399: JREAL = JREAL + 1
400: JIMAG = JIMAG + 1
401: BX( JROW, JCOL ) = DCMPLX( RWORK( JREAL ),
402: $ RWORK( JIMAG ) )
403: 50 CONTINUE
404: 60 CONTINUE
405: *
406: * Since B and BX are complex, the following call to DGEMM
407: * is performed in two steps (real and imaginary parts).
408: *
409: * CALL DGEMM( 'T', 'N', NR, NRHS, NR, ONE, U( NRF, 1 ), LDU,
410: * $ B( NRF, 1 ), LDB, ZERO, BX( NRF, 1 ), LDBX )
411: *
412: J = NR*NRHS*2
413: DO 80 JCOL = 1, NRHS
414: DO 70 JROW = NRF, NRF + NR - 1
415: J = J + 1
416: RWORK( J ) = DBLE( B( JROW, JCOL ) )
417: 70 CONTINUE
418: 80 CONTINUE
419: CALL DGEMM( 'T', 'N', NR, NRHS, NR, ONE, U( NRF, 1 ), LDU,
420: $ RWORK( 1+NR*NRHS*2 ), NR, ZERO, RWORK( 1 ), NR )
421: J = NR*NRHS*2
422: DO 100 JCOL = 1, NRHS
423: DO 90 JROW = NRF, NRF + NR - 1
424: J = J + 1
425: RWORK( J ) = DIMAG( B( JROW, JCOL ) )
426: 90 CONTINUE
427: 100 CONTINUE
428: CALL DGEMM( 'T', 'N', NR, NRHS, NR, ONE, U( NRF, 1 ), LDU,
429: $ RWORK( 1+NR*NRHS*2 ), NR, ZERO, RWORK( 1+NR*NRHS ),
430: $ NR )
431: JREAL = 0
432: JIMAG = NR*NRHS
433: DO 120 JCOL = 1, NRHS
434: DO 110 JROW = NRF, NRF + NR - 1
435: JREAL = JREAL + 1
436: JIMAG = JIMAG + 1
437: BX( JROW, JCOL ) = DCMPLX( RWORK( JREAL ),
438: $ RWORK( JIMAG ) )
439: 110 CONTINUE
440: 120 CONTINUE
441: *
442: 130 CONTINUE
443: *
444: * Next copy the rows of B that correspond to unchanged rows
445: * in the bidiagonal matrix to BX.
446: *
447: DO 140 I = 1, ND
448: IC = IWORK( INODE+I-1 )
449: CALL ZCOPY( NRHS, B( IC, 1 ), LDB, BX( IC, 1 ), LDBX )
450: 140 CONTINUE
451: *
452: * Finally go through the left singular vector matrices of all
453: * the other subproblems bottom-up on the tree.
454: *
455: J = 2**NLVL
456: SQRE = 0
457: *
458: DO 160 LVL = NLVL, 1, -1
459: LVL2 = 2*LVL - 1
460: *
461: * find the first node LF and last node LL on
462: * the current level LVL
463: *
464: IF( LVL.EQ.1 ) THEN
465: LF = 1
466: LL = 1
467: ELSE
468: LF = 2**( LVL-1 )
469: LL = 2*LF - 1
470: END IF
471: DO 150 I = LF, LL
472: IM1 = I - 1
473: IC = IWORK( INODE+IM1 )
474: NL = IWORK( NDIML+IM1 )
475: NR = IWORK( NDIMR+IM1 )
476: NLF = IC - NL
477: NRF = IC + 1
478: J = J - 1
479: CALL ZLALS0( ICOMPQ, NL, NR, SQRE, NRHS, BX( NLF, 1 ), LDBX,
480: $ B( NLF, 1 ), LDB, PERM( NLF, LVL ),
481: $ GIVPTR( J ), GIVCOL( NLF, LVL2 ), LDGCOL,
482: $ GIVNUM( NLF, LVL2 ), LDU, POLES( NLF, LVL2 ),
483: $ DIFL( NLF, LVL ), DIFR( NLF, LVL2 ),
484: $ Z( NLF, LVL ), K( J ), C( J ), S( J ), RWORK,
485: $ INFO )
486: 150 CONTINUE
487: 160 CONTINUE
488: GO TO 330
489: *
490: * ICOMPQ = 1: applying back the right singular vector factors.
491: *
492: 170 CONTINUE
493: *
494: * First now go through the right singular vector matrices of all
495: * the tree nodes top-down.
496: *
497: J = 0
498: DO 190 LVL = 1, NLVL
499: LVL2 = 2*LVL - 1
500: *
501: * Find the first node LF and last node LL on
502: * the current level LVL.
503: *
504: IF( LVL.EQ.1 ) THEN
505: LF = 1
506: LL = 1
507: ELSE
508: LF = 2**( LVL-1 )
509: LL = 2*LF - 1
510: END IF
511: DO 180 I = LL, LF, -1
512: IM1 = I - 1
513: IC = IWORK( INODE+IM1 )
514: NL = IWORK( NDIML+IM1 )
515: NR = IWORK( NDIMR+IM1 )
516: NLF = IC - NL
517: NRF = IC + 1
518: IF( I.EQ.LL ) THEN
519: SQRE = 0
520: ELSE
521: SQRE = 1
522: END IF
523: J = J + 1
524: CALL ZLALS0( ICOMPQ, NL, NR, SQRE, NRHS, B( NLF, 1 ), LDB,
525: $ BX( NLF, 1 ), LDBX, PERM( NLF, LVL ),
526: $ GIVPTR( J ), GIVCOL( NLF, LVL2 ), LDGCOL,
527: $ GIVNUM( NLF, LVL2 ), LDU, POLES( NLF, LVL2 ),
528: $ DIFL( NLF, LVL ), DIFR( NLF, LVL2 ),
529: $ Z( NLF, LVL ), K( J ), C( J ), S( J ), RWORK,
530: $ INFO )
531: 180 CONTINUE
532: 190 CONTINUE
533: *
534: * The nodes on the bottom level of the tree were solved
535: * by DLASDQ. The corresponding right singular vector
536: * matrices are in explicit form. Apply them back.
537: *
538: NDB1 = ( ND+1 ) / 2
539: DO 320 I = NDB1, ND
540: I1 = I - 1
541: IC = IWORK( INODE+I1 )
542: NL = IWORK( NDIML+I1 )
543: NR = IWORK( NDIMR+I1 )
544: NLP1 = NL + 1
545: IF( I.EQ.ND ) THEN
546: NRP1 = NR
547: ELSE
548: NRP1 = NR + 1
549: END IF
550: NLF = IC - NL
551: NRF = IC + 1
552: *
553: * Since B and BX are complex, the following call to DGEMM is
554: * performed in two steps (real and imaginary parts).
555: *
556: * CALL DGEMM( 'T', 'N', NLP1, NRHS, NLP1, ONE, VT( NLF, 1 ), LDU,
557: * $ B( NLF, 1 ), LDB, ZERO, BX( NLF, 1 ), LDBX )
558: *
559: J = NLP1*NRHS*2
560: DO 210 JCOL = 1, NRHS
561: DO 200 JROW = NLF, NLF + NLP1 - 1
562: J = J + 1
563: RWORK( J ) = DBLE( B( JROW, JCOL ) )
564: 200 CONTINUE
565: 210 CONTINUE
566: CALL DGEMM( 'T', 'N', NLP1, NRHS, NLP1, ONE, VT( NLF, 1 ), LDU,
567: $ RWORK( 1+NLP1*NRHS*2 ), NLP1, ZERO, RWORK( 1 ),
568: $ NLP1 )
569: J = NLP1*NRHS*2
570: DO 230 JCOL = 1, NRHS
571: DO 220 JROW = NLF, NLF + NLP1 - 1
572: J = J + 1
573: RWORK( J ) = DIMAG( B( JROW, JCOL ) )
574: 220 CONTINUE
575: 230 CONTINUE
576: CALL DGEMM( 'T', 'N', NLP1, NRHS, NLP1, ONE, VT( NLF, 1 ), LDU,
577: $ RWORK( 1+NLP1*NRHS*2 ), NLP1, ZERO,
578: $ RWORK( 1+NLP1*NRHS ), NLP1 )
579: JREAL = 0
580: JIMAG = NLP1*NRHS
581: DO 250 JCOL = 1, NRHS
582: DO 240 JROW = NLF, NLF + NLP1 - 1
583: JREAL = JREAL + 1
584: JIMAG = JIMAG + 1
585: BX( JROW, JCOL ) = DCMPLX( RWORK( JREAL ),
586: $ RWORK( JIMAG ) )
587: 240 CONTINUE
588: 250 CONTINUE
589: *
590: * Since B and BX are complex, the following call to DGEMM is
591: * performed in two steps (real and imaginary parts).
592: *
593: * CALL DGEMM( 'T', 'N', NRP1, NRHS, NRP1, ONE, VT( NRF, 1 ), LDU,
594: * $ B( NRF, 1 ), LDB, ZERO, BX( NRF, 1 ), LDBX )
595: *
596: J = NRP1*NRHS*2
597: DO 270 JCOL = 1, NRHS
598: DO 260 JROW = NRF, NRF + NRP1 - 1
599: J = J + 1
600: RWORK( J ) = DBLE( B( JROW, JCOL ) )
601: 260 CONTINUE
602: 270 CONTINUE
603: CALL DGEMM( 'T', 'N', NRP1, NRHS, NRP1, ONE, VT( NRF, 1 ), LDU,
604: $ RWORK( 1+NRP1*NRHS*2 ), NRP1, ZERO, RWORK( 1 ),
605: $ NRP1 )
606: J = NRP1*NRHS*2
607: DO 290 JCOL = 1, NRHS
608: DO 280 JROW = NRF, NRF + NRP1 - 1
609: J = J + 1
610: RWORK( J ) = DIMAG( B( JROW, JCOL ) )
611: 280 CONTINUE
612: 290 CONTINUE
613: CALL DGEMM( 'T', 'N', NRP1, NRHS, NRP1, ONE, VT( NRF, 1 ), LDU,
614: $ RWORK( 1+NRP1*NRHS*2 ), NRP1, ZERO,
615: $ RWORK( 1+NRP1*NRHS ), NRP1 )
616: JREAL = 0
617: JIMAG = NRP1*NRHS
618: DO 310 JCOL = 1, NRHS
619: DO 300 JROW = NRF, NRF + NRP1 - 1
620: JREAL = JREAL + 1
621: JIMAG = JIMAG + 1
622: BX( JROW, JCOL ) = DCMPLX( RWORK( JREAL ),
623: $ RWORK( JIMAG ) )
624: 300 CONTINUE
625: 310 CONTINUE
626: *
627: 320 CONTINUE
628: *
629: 330 CONTINUE
630: *
631: RETURN
632: *
633: * End of ZLALSA
634: *
635: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zlalsa.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack