Annotation of rpl/lapack/lapack/zlalsa.f, revision 1.10
1.10 ! bertrand 1: *> \brief \b ZLALSA
! 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: *
! 254: *> \date November 2011
! 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.10 ! bertrand 271: * -- LAPACK computational routine (version 3.4.0) --
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.10 ! bertrand 274: * November 2011
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