File:  [local] / rpl / lapack / lapack / zlalsd.f
Revision 1.20: download - view: text, annotated - select for diffs - revision graph
Tue May 29 07:18:26 2018 UTC (5 years, 11 months ago) by bertrand
Branches: MAIN
CVS tags: rpl-4_1_33, rpl-4_1_32, rpl-4_1_31, rpl-4_1_30, rpl-4_1_29, rpl-4_1_28, HEAD
Mise à jour de Lapack.

    1: *> \brief \b ZLALSD uses the singular value decomposition of A to solve the least squares problem.
    2: *
    3: *  =========== DOCUMENTATION ===========
    4: *
    5: * Online html documentation available at
    6: *            http://www.netlib.org/lapack/explore-html/
    7: *
    8: *> \htmlonly
    9: *> Download ZLALSD + dependencies
   10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlalsd.f">
   11: *> [TGZ]</a>
   12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlalsd.f">
   13: *> [ZIP]</a>
   14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlalsd.f">
   15: *> [TXT]</a>
   16: *> \endhtmlonly
   17: *
   18: *  Definition:
   19: *  ===========
   20: *
   21: *       SUBROUTINE ZLALSD( UPLO, SMLSIZ, N, NRHS, D, E, B, LDB, RCOND,
   22: *                          RANK, WORK, RWORK, IWORK, INFO )
   23: *
   24: *       .. Scalar Arguments ..
   25: *       CHARACTER          UPLO
   26: *       INTEGER            INFO, LDB, N, NRHS, RANK, SMLSIZ
   27: *       DOUBLE PRECISION   RCOND
   28: *       ..
   29: *       .. Array Arguments ..
   30: *       INTEGER            IWORK( * )
   31: *       DOUBLE PRECISION   D( * ), E( * ), RWORK( * )
   32: *       COMPLEX*16         B( LDB, * ), WORK( * )
   33: *       ..
   34: *
   35: *
   36: *> \par Purpose:
   37: *  =============
   38: *>
   39: *> \verbatim
   40: *>
   41: *> ZLALSD uses the singular value decomposition of A to solve the least
   42: *> squares problem of finding X to minimize the Euclidean norm of each
   43: *> column of A*X-B, where A is N-by-N upper bidiagonal, and X and B
   44: *> are N-by-NRHS. The solution X overwrites B.
   45: *>
   46: *> The singular values of A smaller than RCOND times the largest
   47: *> singular value are treated as zero in solving the least squares
   48: *> problem; in this case a minimum norm solution is returned.
   49: *> The actual singular values are returned in D in ascending order.
   50: *>
   51: *> This code makes very mild assumptions about floating point
   52: *> arithmetic. It will work on machines with a guard digit in
   53: *> add/subtract, or on those binary machines without guard digits
   54: *> which subtract like the Cray XMP, Cray YMP, Cray C 90, or Cray 2.
   55: *> It could conceivably fail on hexadecimal or decimal machines
   56: *> without guard digits, but we know of none.
   57: *> \endverbatim
   58: *
   59: *  Arguments:
   60: *  ==========
   61: *
   62: *> \param[in] UPLO
   63: *> \verbatim
   64: *>          UPLO is CHARACTER*1
   65: *>         = 'U': D and E define an upper bidiagonal matrix.
   66: *>         = 'L': D and E define a  lower bidiagonal 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 dimension of the  bidiagonal matrix.  N >= 0.
   80: *> \endverbatim
   81: *>
   82: *> \param[in] NRHS
   83: *> \verbatim
   84: *>          NRHS is INTEGER
   85: *>         The number of columns of B. NRHS must be at least 1.
   86: *> \endverbatim
   87: *>
   88: *> \param[in,out] D
   89: *> \verbatim
   90: *>          D is DOUBLE PRECISION array, dimension (N)
   91: *>         On entry D contains the main diagonal of the bidiagonal
   92: *>         matrix. On exit, if INFO = 0, D contains its singular values.
   93: *> \endverbatim
   94: *>
   95: *> \param[in,out] E
   96: *> \verbatim
   97: *>          E is DOUBLE PRECISION array, dimension (N-1)
   98: *>         Contains the super-diagonal entries of the bidiagonal matrix.
   99: *>         On exit, E has been destroyed.
  100: *> \endverbatim
  101: *>
  102: *> \param[in,out] B
  103: *> \verbatim
  104: *>          B is COMPLEX*16 array, dimension (LDB,NRHS)
  105: *>         On input, B contains the right hand sides of the least
  106: *>         squares problem. On output, B contains the solution X.
  107: *> \endverbatim
  108: *>
  109: *> \param[in] LDB
  110: *> \verbatim
  111: *>          LDB is INTEGER
  112: *>         The leading dimension of B in the calling subprogram.
  113: *>         LDB must be at least max(1,N).
  114: *> \endverbatim
  115: *>
  116: *> \param[in] RCOND
  117: *> \verbatim
  118: *>          RCOND is DOUBLE PRECISION
  119: *>         The singular values of A less than or equal to RCOND times
  120: *>         the largest singular value are treated as zero in solving
  121: *>         the least squares problem. If RCOND is negative,
  122: *>         machine precision is used instead.
  123: *>         For example, if diag(S)*X=B were the least squares problem,
  124: *>         where diag(S) is a diagonal matrix of singular values, the
  125: *>         solution would be X(i) = B(i) / S(i) if S(i) is greater than
  126: *>         RCOND*max(S), and X(i) = 0 if S(i) is less than or equal to
  127: *>         RCOND*max(S).
  128: *> \endverbatim
  129: *>
  130: *> \param[out] RANK
  131: *> \verbatim
  132: *>          RANK is INTEGER
  133: *>         The number of singular values of A greater than RCOND times
  134: *>         the largest singular value.
  135: *> \endverbatim
  136: *>
  137: *> \param[out] WORK
  138: *> \verbatim
  139: *>          WORK is COMPLEX*16 array, dimension (N * NRHS)
  140: *> \endverbatim
  141: *>
  142: *> \param[out] RWORK
  143: *> \verbatim
  144: *>          RWORK is DOUBLE PRECISION array, dimension at least
  145: *>         (9*N + 2*N*SMLSIZ + 8*N*NLVL + 3*SMLSIZ*NRHS +
  146: *>         MAX( (SMLSIZ+1)**2, N*(1+NRHS) + 2*NRHS ),
  147: *>         where
  148: *>         NLVL = MAX( 0, INT( LOG_2( MIN( M,N )/(SMLSIZ+1) ) ) + 1 )
  149: *> \endverbatim
  150: *>
  151: *> \param[out] IWORK
  152: *> \verbatim
  153: *>          IWORK is INTEGER array, dimension at least
  154: *>         (3*N*NLVL + 11*N).
  155: *> \endverbatim
  156: *>
  157: *> \param[out] INFO
  158: *> \verbatim
  159: *>          INFO is INTEGER
  160: *>         = 0:  successful exit.
  161: *>         < 0:  if INFO = -i, the i-th argument had an illegal value.
  162: *>         > 0:  The algorithm failed to compute a singular value while
  163: *>               working on the submatrix lying in rows and columns
  164: *>               INFO/(N+1) through MOD(INFO,N+1).
  165: *> \endverbatim
  166: *
  167: *  Authors:
  168: *  ========
  169: *
  170: *> \author Univ. of Tennessee
  171: *> \author Univ. of California Berkeley
  172: *> \author Univ. of Colorado Denver
  173: *> \author NAG Ltd.
  174: *
  175: *> \date June 2017
  176: *
  177: *> \ingroup complex16OTHERcomputational
  178: *
  179: *> \par Contributors:
  180: *  ==================
  181: *>
  182: *>     Ming Gu and Ren-Cang Li, Computer Science Division, University of
  183: *>       California at Berkeley, USA \n
  184: *>     Osni Marques, LBNL/NERSC, USA \n
  185: *
  186: *  =====================================================================
  187:       SUBROUTINE ZLALSD( UPLO, SMLSIZ, N, NRHS, D, E, B, LDB, RCOND,
  188:      $                   RANK, WORK, RWORK, IWORK, INFO )
  189: *
  190: *  -- LAPACK computational routine (version 3.7.1) --
  191: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  192: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  193: *     June 2017
  194: *
  195: *     .. Scalar Arguments ..
  196:       CHARACTER          UPLO
  197:       INTEGER            INFO, LDB, N, NRHS, RANK, SMLSIZ
  198:       DOUBLE PRECISION   RCOND
  199: *     ..
  200: *     .. Array Arguments ..
  201:       INTEGER            IWORK( * )
  202:       DOUBLE PRECISION   D( * ), E( * ), RWORK( * )
  203:       COMPLEX*16         B( LDB, * ), WORK( * )
  204: *     ..
  205: *
  206: *  =====================================================================
  207: *
  208: *     .. Parameters ..
  209:       DOUBLE PRECISION   ZERO, ONE, TWO
  210:       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0 )
  211:       COMPLEX*16         CZERO
  212:       PARAMETER          ( CZERO = ( 0.0D0, 0.0D0 ) )
  213: *     ..
  214: *     .. Local Scalars ..
  215:       INTEGER            BX, BXST, C, DIFL, DIFR, GIVCOL, GIVNUM,
  216:      $                   GIVPTR, I, ICMPQ1, ICMPQ2, IRWB, IRWIB, IRWRB,
  217:      $                   IRWU, IRWVT, IRWWRK, IWK, J, JCOL, JIMAG,
  218:      $                   JREAL, JROW, K, NLVL, NM1, NRWORK, NSIZE, NSUB,
  219:      $                   PERM, POLES, S, SIZEI, SMLSZP, SQRE, ST, ST1,
  220:      $                   U, VT, Z
  221:       DOUBLE PRECISION   CS, EPS, ORGNRM, RCND, R, SN, TOL
  222: *     ..
  223: *     .. External Functions ..
  224:       INTEGER            IDAMAX
  225:       DOUBLE PRECISION   DLAMCH, DLANST
  226:       EXTERNAL           IDAMAX, DLAMCH, DLANST
  227: *     ..
  228: *     .. External Subroutines ..
  229:       EXTERNAL           DGEMM, DLARTG, DLASCL, DLASDA, DLASDQ, DLASET,
  230:      $                   DLASRT, XERBLA, ZCOPY, ZDROT, ZLACPY, ZLALSA,
  231:      $                   ZLASCL, ZLASET
  232: *     ..
  233: *     .. Intrinsic Functions ..
  234:       INTRINSIC          ABS, DBLE, DCMPLX, DIMAG, INT, LOG, SIGN
  235: *     ..
  236: *     .. Executable Statements ..
  237: *
  238: *     Test the input parameters.
  239: *
  240:       INFO = 0
  241: *
  242:       IF( N.LT.0 ) THEN
  243:          INFO = -3
  244:       ELSE IF( NRHS.LT.1 ) THEN
  245:          INFO = -4
  246:       ELSE IF( ( LDB.LT.1 ) .OR. ( LDB.LT.N ) ) THEN
  247:          INFO = -8
  248:       END IF
  249:       IF( INFO.NE.0 ) THEN
  250:          CALL XERBLA( 'ZLALSD', -INFO )
  251:          RETURN
  252:       END IF
  253: *
  254:       EPS = DLAMCH( 'Epsilon' )
  255: *
  256: *     Set up the tolerance.
  257: *
  258:       IF( ( RCOND.LE.ZERO ) .OR. ( RCOND.GE.ONE ) ) THEN
  259:          RCND = EPS
  260:       ELSE
  261:          RCND = RCOND
  262:       END IF
  263: *
  264:       RANK = 0
  265: *
  266: *     Quick return if possible.
  267: *
  268:       IF( N.EQ.0 ) THEN
  269:          RETURN
  270:       ELSE IF( N.EQ.1 ) THEN
  271:          IF( D( 1 ).EQ.ZERO ) THEN
  272:             CALL ZLASET( 'A', 1, NRHS, CZERO, CZERO, B, LDB )
  273:          ELSE
  274:             RANK = 1
  275:             CALL ZLASCL( 'G', 0, 0, D( 1 ), ONE, 1, NRHS, B, LDB, INFO )
  276:             D( 1 ) = ABS( D( 1 ) )
  277:          END IF
  278:          RETURN
  279:       END IF
  280: *
  281: *     Rotate the matrix if it is lower bidiagonal.
  282: *
  283:       IF( UPLO.EQ.'L' ) THEN
  284:          DO 10 I = 1, N - 1
  285:             CALL DLARTG( D( I ), E( I ), CS, SN, R )
  286:             D( I ) = R
  287:             E( I ) = SN*D( I+1 )
  288:             D( I+1 ) = CS*D( I+1 )
  289:             IF( NRHS.EQ.1 ) THEN
  290:                CALL ZDROT( 1, B( I, 1 ), 1, B( I+1, 1 ), 1, CS, SN )
  291:             ELSE
  292:                RWORK( I*2-1 ) = CS
  293:                RWORK( I*2 ) = SN
  294:             END IF
  295:    10    CONTINUE
  296:          IF( NRHS.GT.1 ) THEN
  297:             DO 30 I = 1, NRHS
  298:                DO 20 J = 1, N - 1
  299:                   CS = RWORK( J*2-1 )
  300:                   SN = RWORK( J*2 )
  301:                   CALL ZDROT( 1, B( J, I ), 1, B( J+1, I ), 1, CS, SN )
  302:    20          CONTINUE
  303:    30       CONTINUE
  304:          END IF
  305:       END IF
  306: *
  307: *     Scale.
  308: *
  309:       NM1 = N - 1
  310:       ORGNRM = DLANST( 'M', N, D, E )
  311:       IF( ORGNRM.EQ.ZERO ) THEN
  312:          CALL ZLASET( 'A', N, NRHS, CZERO, CZERO, B, LDB )
  313:          RETURN
  314:       END IF
  315: *
  316:       CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, N, 1, D, N, INFO )
  317:       CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, NM1, 1, E, NM1, INFO )
  318: *
  319: *     If N is smaller than the minimum divide size SMLSIZ, then solve
  320: *     the problem with another solver.
  321: *
  322:       IF( N.LE.SMLSIZ ) THEN
  323:          IRWU = 1
  324:          IRWVT = IRWU + N*N
  325:          IRWWRK = IRWVT + N*N
  326:          IRWRB = IRWWRK
  327:          IRWIB = IRWRB + N*NRHS
  328:          IRWB = IRWIB + N*NRHS
  329:          CALL DLASET( 'A', N, N, ZERO, ONE, RWORK( IRWU ), N )
  330:          CALL DLASET( 'A', N, N, ZERO, ONE, RWORK( IRWVT ), N )
  331:          CALL DLASDQ( 'U', 0, N, N, N, 0, D, E, RWORK( IRWVT ), N,
  332:      $                RWORK( IRWU ), N, RWORK( IRWWRK ), 1,
  333:      $                RWORK( IRWWRK ), INFO )
  334:          IF( INFO.NE.0 ) THEN
  335:             RETURN
  336:          END IF
  337: *
  338: *        In the real version, B is passed to DLASDQ and multiplied
  339: *        internally by Q**H. Here B is complex and that product is
  340: *        computed below in two steps (real and imaginary parts).
  341: *
  342:          J = IRWB - 1
  343:          DO 50 JCOL = 1, NRHS
  344:             DO 40 JROW = 1, N
  345:                J = J + 1
  346:                RWORK( J ) = DBLE( B( JROW, JCOL ) )
  347:    40       CONTINUE
  348:    50    CONTINUE
  349:          CALL DGEMM( 'T', 'N', N, NRHS, N, ONE, RWORK( IRWU ), N,
  350:      $               RWORK( IRWB ), N, ZERO, RWORK( IRWRB ), N )
  351:          J = IRWB - 1
  352:          DO 70 JCOL = 1, NRHS
  353:             DO 60 JROW = 1, N
  354:                J = J + 1
  355:                RWORK( J ) = DIMAG( B( JROW, JCOL ) )
  356:    60       CONTINUE
  357:    70    CONTINUE
  358:          CALL DGEMM( 'T', 'N', N, NRHS, N, ONE, RWORK( IRWU ), N,
  359:      $               RWORK( IRWB ), N, ZERO, RWORK( IRWIB ), N )
  360:          JREAL = IRWRB - 1
  361:          JIMAG = IRWIB - 1
  362:          DO 90 JCOL = 1, NRHS
  363:             DO 80 JROW = 1, N
  364:                JREAL = JREAL + 1
  365:                JIMAG = JIMAG + 1
  366:                B( JROW, JCOL ) = DCMPLX( RWORK( JREAL ),
  367:      $                           RWORK( JIMAG ) )
  368:    80       CONTINUE
  369:    90    CONTINUE
  370: *
  371:          TOL = RCND*ABS( D( IDAMAX( N, D, 1 ) ) )
  372:          DO 100 I = 1, N
  373:             IF( D( I ).LE.TOL ) THEN
  374:                CALL ZLASET( 'A', 1, NRHS, CZERO, CZERO, B( I, 1 ), LDB )
  375:             ELSE
  376:                CALL ZLASCL( 'G', 0, 0, D( I ), ONE, 1, NRHS, B( I, 1 ),
  377:      $                      LDB, INFO )
  378:                RANK = RANK + 1
  379:             END IF
  380:   100    CONTINUE
  381: *
  382: *        Since B is complex, the following call to DGEMM is performed
  383: *        in two steps (real and imaginary parts). That is for V * B
  384: *        (in the real version of the code V**H is stored in WORK).
  385: *
  386: *        CALL DGEMM( 'T', 'N', N, NRHS, N, ONE, WORK, N, B, LDB, ZERO,
  387: *    $               WORK( NWORK ), N )
  388: *
  389:          J = IRWB - 1
  390:          DO 120 JCOL = 1, NRHS
  391:             DO 110 JROW = 1, N
  392:                J = J + 1
  393:                RWORK( J ) = DBLE( B( JROW, JCOL ) )
  394:   110       CONTINUE
  395:   120    CONTINUE
  396:          CALL DGEMM( 'T', 'N', N, NRHS, N, ONE, RWORK( IRWVT ), N,
  397:      $               RWORK( IRWB ), N, ZERO, RWORK( IRWRB ), N )
  398:          J = IRWB - 1
  399:          DO 140 JCOL = 1, NRHS
  400:             DO 130 JROW = 1, N
  401:                J = J + 1
  402:                RWORK( J ) = DIMAG( B( JROW, JCOL ) )
  403:   130       CONTINUE
  404:   140    CONTINUE
  405:          CALL DGEMM( 'T', 'N', N, NRHS, N, ONE, RWORK( IRWVT ), N,
  406:      $               RWORK( IRWB ), N, ZERO, RWORK( IRWIB ), N )
  407:          JREAL = IRWRB - 1
  408:          JIMAG = IRWIB - 1
  409:          DO 160 JCOL = 1, NRHS
  410:             DO 150 JROW = 1, N
  411:                JREAL = JREAL + 1
  412:                JIMAG = JIMAG + 1
  413:                B( JROW, JCOL ) = DCMPLX( RWORK( JREAL ),
  414:      $                           RWORK( JIMAG ) )
  415:   150       CONTINUE
  416:   160    CONTINUE
  417: *
  418: *        Unscale.
  419: *
  420:          CALL DLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, INFO )
  421:          CALL DLASRT( 'D', N, D, INFO )
  422:          CALL ZLASCL( 'G', 0, 0, ORGNRM, ONE, N, NRHS, B, LDB, INFO )
  423: *
  424:          RETURN
  425:       END IF
  426: *
  427: *     Book-keeping and setting up some constants.
  428: *
  429:       NLVL = INT( LOG( DBLE( N ) / DBLE( SMLSIZ+1 ) ) / LOG( TWO ) ) + 1
  430: *
  431:       SMLSZP = SMLSIZ + 1
  432: *
  433:       U = 1
  434:       VT = 1 + SMLSIZ*N
  435:       DIFL = VT + SMLSZP*N
  436:       DIFR = DIFL + NLVL*N
  437:       Z = DIFR + NLVL*N*2
  438:       C = Z + NLVL*N
  439:       S = C + N
  440:       POLES = S + N
  441:       GIVNUM = POLES + 2*NLVL*N
  442:       NRWORK = GIVNUM + 2*NLVL*N
  443:       BX = 1
  444: *
  445:       IRWRB = NRWORK
  446:       IRWIB = IRWRB + SMLSIZ*NRHS
  447:       IRWB = IRWIB + SMLSIZ*NRHS
  448: *
  449:       SIZEI = 1 + N
  450:       K = SIZEI + N
  451:       GIVPTR = K + N
  452:       PERM = GIVPTR + N
  453:       GIVCOL = PERM + NLVL*N
  454:       IWK = GIVCOL + NLVL*N*2
  455: *
  456:       ST = 1
  457:       SQRE = 0
  458:       ICMPQ1 = 1
  459:       ICMPQ2 = 0
  460:       NSUB = 0
  461: *
  462:       DO 170 I = 1, N
  463:          IF( ABS( D( I ) ).LT.EPS ) THEN
  464:             D( I ) = SIGN( EPS, D( I ) )
  465:          END IF
  466:   170 CONTINUE
  467: *
  468:       DO 240 I = 1, NM1
  469:          IF( ( ABS( E( I ) ).LT.EPS ) .OR. ( I.EQ.NM1 ) ) THEN
  470:             NSUB = NSUB + 1
  471:             IWORK( NSUB ) = ST
  472: *
  473: *           Subproblem found. First determine its size and then
  474: *           apply divide and conquer on it.
  475: *
  476:             IF( I.LT.NM1 ) THEN
  477: *
  478: *              A subproblem with E(I) small for I < NM1.
  479: *
  480:                NSIZE = I - ST + 1
  481:                IWORK( SIZEI+NSUB-1 ) = NSIZE
  482:             ELSE IF( ABS( E( I ) ).GE.EPS ) THEN
  483: *
  484: *              A subproblem with E(NM1) not too small but I = NM1.
  485: *
  486:                NSIZE = N - ST + 1
  487:                IWORK( SIZEI+NSUB-1 ) = NSIZE
  488:             ELSE
  489: *
  490: *              A subproblem with E(NM1) small. This implies an
  491: *              1-by-1 subproblem at D(N), which is not solved
  492: *              explicitly.
  493: *
  494:                NSIZE = I - ST + 1
  495:                IWORK( SIZEI+NSUB-1 ) = NSIZE
  496:                NSUB = NSUB + 1
  497:                IWORK( NSUB ) = N
  498:                IWORK( SIZEI+NSUB-1 ) = 1
  499:                CALL ZCOPY( NRHS, B( N, 1 ), LDB, WORK( BX+NM1 ), N )
  500:             END IF
  501:             ST1 = ST - 1
  502:             IF( NSIZE.EQ.1 ) THEN
  503: *
  504: *              This is a 1-by-1 subproblem and is not solved
  505: *              explicitly.
  506: *
  507:                CALL ZCOPY( NRHS, B( ST, 1 ), LDB, WORK( BX+ST1 ), N )
  508:             ELSE IF( NSIZE.LE.SMLSIZ ) THEN
  509: *
  510: *              This is a small subproblem and is solved by DLASDQ.
  511: *
  512:                CALL DLASET( 'A', NSIZE, NSIZE, ZERO, ONE,
  513:      $                      RWORK( VT+ST1 ), N )
  514:                CALL DLASET( 'A', NSIZE, NSIZE, ZERO, ONE,
  515:      $                      RWORK( U+ST1 ), N )
  516:                CALL DLASDQ( 'U', 0, NSIZE, NSIZE, NSIZE, 0, D( ST ),
  517:      $                      E( ST ), RWORK( VT+ST1 ), N, RWORK( U+ST1 ),
  518:      $                      N, RWORK( NRWORK ), 1, RWORK( NRWORK ),
  519:      $                      INFO )
  520:                IF( INFO.NE.0 ) THEN
  521:                   RETURN
  522:                END IF
  523: *
  524: *              In the real version, B is passed to DLASDQ and multiplied
  525: *              internally by Q**H. Here B is complex and that product is
  526: *              computed below in two steps (real and imaginary parts).
  527: *
  528:                J = IRWB - 1
  529:                DO 190 JCOL = 1, NRHS
  530:                   DO 180 JROW = ST, ST + NSIZE - 1
  531:                      J = J + 1
  532:                      RWORK( J ) = DBLE( B( JROW, JCOL ) )
  533:   180             CONTINUE
  534:   190          CONTINUE
  535:                CALL DGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE,
  536:      $                     RWORK( U+ST1 ), N, RWORK( IRWB ), NSIZE,
  537:      $                     ZERO, RWORK( IRWRB ), NSIZE )
  538:                J = IRWB - 1
  539:                DO 210 JCOL = 1, NRHS
  540:                   DO 200 JROW = ST, ST + NSIZE - 1
  541:                      J = J + 1
  542:                      RWORK( J ) = DIMAG( B( JROW, JCOL ) )
  543:   200             CONTINUE
  544:   210          CONTINUE
  545:                CALL DGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE,
  546:      $                     RWORK( U+ST1 ), N, RWORK( IRWB ), NSIZE,
  547:      $                     ZERO, RWORK( IRWIB ), NSIZE )
  548:                JREAL = IRWRB - 1
  549:                JIMAG = IRWIB - 1
  550:                DO 230 JCOL = 1, NRHS
  551:                   DO 220 JROW = ST, ST + NSIZE - 1
  552:                      JREAL = JREAL + 1
  553:                      JIMAG = JIMAG + 1
  554:                      B( JROW, JCOL ) = DCMPLX( RWORK( JREAL ),
  555:      $                                 RWORK( JIMAG ) )
  556:   220             CONTINUE
  557:   230          CONTINUE
  558: *
  559:                CALL ZLACPY( 'A', NSIZE, NRHS, B( ST, 1 ), LDB,
  560:      $                      WORK( BX+ST1 ), N )
  561:             ELSE
  562: *
  563: *              A large problem. Solve it using divide and conquer.
  564: *
  565:                CALL DLASDA( ICMPQ1, SMLSIZ, NSIZE, SQRE, D( ST ),
  566:      $                      E( ST ), RWORK( U+ST1 ), N, RWORK( VT+ST1 ),
  567:      $                      IWORK( K+ST1 ), RWORK( DIFL+ST1 ),
  568:      $                      RWORK( DIFR+ST1 ), RWORK( Z+ST1 ),
  569:      $                      RWORK( POLES+ST1 ), IWORK( GIVPTR+ST1 ),
  570:      $                      IWORK( GIVCOL+ST1 ), N, IWORK( PERM+ST1 ),
  571:      $                      RWORK( GIVNUM+ST1 ), RWORK( C+ST1 ),
  572:      $                      RWORK( S+ST1 ), RWORK( NRWORK ),
  573:      $                      IWORK( IWK ), INFO )
  574:                IF( INFO.NE.0 ) THEN
  575:                   RETURN
  576:                END IF
  577:                BXST = BX + ST1
  578:                CALL ZLALSA( ICMPQ2, SMLSIZ, NSIZE, NRHS, B( ST, 1 ),
  579:      $                      LDB, WORK( BXST ), N, RWORK( U+ST1 ), N,
  580:      $                      RWORK( VT+ST1 ), IWORK( K+ST1 ),
  581:      $                      RWORK( DIFL+ST1 ), RWORK( DIFR+ST1 ),
  582:      $                      RWORK( Z+ST1 ), RWORK( POLES+ST1 ),
  583:      $                      IWORK( GIVPTR+ST1 ), IWORK( GIVCOL+ST1 ), N,
  584:      $                      IWORK( PERM+ST1 ), RWORK( GIVNUM+ST1 ),
  585:      $                      RWORK( C+ST1 ), RWORK( S+ST1 ),
  586:      $                      RWORK( NRWORK ), IWORK( IWK ), INFO )
  587:                IF( INFO.NE.0 ) THEN
  588:                   RETURN
  589:                END IF
  590:             END IF
  591:             ST = I + 1
  592:          END IF
  593:   240 CONTINUE
  594: *
  595: *     Apply the singular values and treat the tiny ones as zero.
  596: *
  597:       TOL = RCND*ABS( D( IDAMAX( N, D, 1 ) ) )
  598: *
  599:       DO 250 I = 1, N
  600: *
  601: *        Some of the elements in D can be negative because 1-by-1
  602: *        subproblems were not solved explicitly.
  603: *
  604:          IF( ABS( D( I ) ).LE.TOL ) THEN
  605:             CALL ZLASET( 'A', 1, NRHS, CZERO, CZERO, WORK( BX+I-1 ), N )
  606:          ELSE
  607:             RANK = RANK + 1
  608:             CALL ZLASCL( 'G', 0, 0, D( I ), ONE, 1, NRHS,
  609:      $                   WORK( BX+I-1 ), N, INFO )
  610:          END IF
  611:          D( I ) = ABS( D( I ) )
  612:   250 CONTINUE
  613: *
  614: *     Now apply back the right singular vectors.
  615: *
  616:       ICMPQ2 = 1
  617:       DO 320 I = 1, NSUB
  618:          ST = IWORK( I )
  619:          ST1 = ST - 1
  620:          NSIZE = IWORK( SIZEI+I-1 )
  621:          BXST = BX + ST1
  622:          IF( NSIZE.EQ.1 ) THEN
  623:             CALL ZCOPY( NRHS, WORK( BXST ), N, B( ST, 1 ), LDB )
  624:          ELSE IF( NSIZE.LE.SMLSIZ ) THEN
  625: *
  626: *           Since B and BX are complex, the following call to DGEMM
  627: *           is performed in two steps (real and imaginary parts).
  628: *
  629: *           CALL DGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE,
  630: *    $                  RWORK( VT+ST1 ), N, RWORK( BXST ), N, ZERO,
  631: *    $                  B( ST, 1 ), LDB )
  632: *
  633:             J = BXST - N - 1
  634:             JREAL = IRWB - 1
  635:             DO 270 JCOL = 1, NRHS
  636:                J = J + N
  637:                DO 260 JROW = 1, NSIZE
  638:                   JREAL = JREAL + 1
  639:                   RWORK( JREAL ) = DBLE( WORK( J+JROW ) )
  640:   260          CONTINUE
  641:   270       CONTINUE
  642:             CALL DGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE,
  643:      $                  RWORK( VT+ST1 ), N, RWORK( IRWB ), NSIZE, ZERO,
  644:      $                  RWORK( IRWRB ), NSIZE )
  645:             J = BXST - N - 1
  646:             JIMAG = IRWB - 1
  647:             DO 290 JCOL = 1, NRHS
  648:                J = J + N
  649:                DO 280 JROW = 1, NSIZE
  650:                   JIMAG = JIMAG + 1
  651:                   RWORK( JIMAG ) = DIMAG( WORK( J+JROW ) )
  652:   280          CONTINUE
  653:   290       CONTINUE
  654:             CALL DGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE,
  655:      $                  RWORK( VT+ST1 ), N, RWORK( IRWB ), NSIZE, ZERO,
  656:      $                  RWORK( IRWIB ), NSIZE )
  657:             JREAL = IRWRB - 1
  658:             JIMAG = IRWIB - 1
  659:             DO 310 JCOL = 1, NRHS
  660:                DO 300 JROW = ST, ST + NSIZE - 1
  661:                   JREAL = JREAL + 1
  662:                   JIMAG = JIMAG + 1
  663:                   B( JROW, JCOL ) = DCMPLX( RWORK( JREAL ),
  664:      $                              RWORK( JIMAG ) )
  665:   300          CONTINUE
  666:   310       CONTINUE
  667:          ELSE
  668:             CALL ZLALSA( ICMPQ2, SMLSIZ, NSIZE, NRHS, WORK( BXST ), N,
  669:      $                   B( ST, 1 ), LDB, RWORK( U+ST1 ), N,
  670:      $                   RWORK( VT+ST1 ), IWORK( K+ST1 ),
  671:      $                   RWORK( DIFL+ST1 ), RWORK( DIFR+ST1 ),
  672:      $                   RWORK( Z+ST1 ), RWORK( POLES+ST1 ),
  673:      $                   IWORK( GIVPTR+ST1 ), IWORK( GIVCOL+ST1 ), N,
  674:      $                   IWORK( PERM+ST1 ), RWORK( GIVNUM+ST1 ),
  675:      $                   RWORK( C+ST1 ), RWORK( S+ST1 ),
  676:      $                   RWORK( NRWORK ), IWORK( IWK ), INFO )
  677:             IF( INFO.NE.0 ) THEN
  678:                RETURN
  679:             END IF
  680:          END IF
  681:   320 CONTINUE
  682: *
  683: *     Unscale and sort the singular values.
  684: *
  685:       CALL DLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, INFO )
  686:       CALL DLASRT( 'D', N, D, INFO )
  687:       CALL ZLASCL( 'G', 0, 0, ORGNRM, ONE, N, NRHS, B, LDB, INFO )
  688: *
  689:       RETURN
  690: *
  691: *     End of ZLALSD
  692: *
  693:       END

CVSweb interface <joel.bertrand@systella.fr>