Annotation of rpl/lapack/lapack/zlalsa.f, revision 1.16
1.13 bertrand 1: *> \brief \b ZLALSA computes the SVD of the coefficient matrix in compact form. Used by sgelsd.
1.10 bertrand 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.
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: *
1.13 bertrand 254: *> \date September 2012
1.10 bertrand 255: *
256: *> \ingroup complex16OTHERcomputational
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: * =====================================================================
1.1 bertrand 266: SUBROUTINE ZLALSA( 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: *
1.13 bertrand 271: * -- LAPACK computational routine (version 3.4.2) --
1.1 bertrand 272: * -- LAPACK is a software package provided by Univ. of Tennessee, --
273: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.13 bertrand 274: * September 2012
1.1 bertrand 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: DOUBLE PRECISION C( * ), DIFL( LDU, * ), DIFR( LDU, * ),
284: $ GIVNUM( LDU, * ), POLES( LDU, * ), RWORK( * ),
285: $ S( * ), U( LDU, * ), VT( LDU, * ), Z( LDU, * )
286: COMPLEX*16 B( LDB, * ), BX( LDBX, * )
287: * ..
288: *
289: * =====================================================================
290: *
291: * .. Parameters ..
292: DOUBLE PRECISION ZERO, ONE
293: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
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 DGEMM, DLASDT, XERBLA, ZCOPY, ZLALS0
302: * ..
303: * .. Intrinsic Functions ..
304: INTRINSIC DBLE, DCMPLX, DIMAG
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( 'ZLALSA', -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 DLASDT( 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 DLASDQ. 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 DGEMM
372: * is performed in two steps (real and imaginary parts).
373: *
374: * CALL DGEMM( '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 ) = DBLE( B( JROW, JCOL ) )
382: 10 CONTINUE
383: 20 CONTINUE
384: CALL DGEMM( '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 ) = DIMAG( B( JROW, JCOL ) )
391: 30 CONTINUE
392: 40 CONTINUE
393: CALL DGEMM( '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 ) = DCMPLX( RWORK( JREAL ),
403: $ RWORK( JIMAG ) )
404: 50 CONTINUE
405: 60 CONTINUE
406: *
407: * Since B and BX are complex, the following call to DGEMM
408: * is performed in two steps (real and imaginary parts).
409: *
410: * CALL DGEMM( '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 ) = DBLE( B( JROW, JCOL ) )
418: 70 CONTINUE
419: 80 CONTINUE
420: CALL DGEMM( '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 ) = DIMAG( B( JROW, JCOL ) )
427: 90 CONTINUE
428: 100 CONTINUE
429: CALL DGEMM( '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 ) = DCMPLX( 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 ZCOPY( 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 ZLALS0( 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 ZLALS0( 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 DLASDQ. 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 DGEMM is
555: * performed in two steps (real and imaginary parts).
556: *
557: * CALL DGEMM( '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 ) = DBLE( B( JROW, JCOL ) )
565: 200 CONTINUE
566: 210 CONTINUE
567: CALL DGEMM( '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 ) = DIMAG( B( JROW, JCOL ) )
575: 220 CONTINUE
576: 230 CONTINUE
577: CALL DGEMM( '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 ) = DCMPLX( RWORK( JREAL ),
587: $ RWORK( JIMAG ) )
588: 240 CONTINUE
589: 250 CONTINUE
590: *
591: * Since B and BX are complex, the following call to DGEMM is
592: * performed in two steps (real and imaginary parts).
593: *
594: * CALL DGEMM( '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 ) = DBLE( B( JROW, JCOL ) )
602: 260 CONTINUE
603: 270 CONTINUE
604: CALL DGEMM( '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 ) = DIMAG( B( JROW, JCOL ) )
612: 280 CONTINUE
613: 290 CONTINUE
614: CALL DGEMM( '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 ) = DCMPLX( 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 ZLALSA
635: *
636: 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