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