File:  [local] / rpl / lapack / lapack / zlags2.f
Revision 1.16: download - view: text, annotated - select for diffs - revision graph
Sat Jun 17 11:06:53 2017 UTC (6 years, 10 months ago) by bertrand
Branches: MAIN
CVS tags: rpl-4_1_27, rpl-4_1_26, HEAD
Cohérence.

    1: *> \brief \b ZLAGS2
    2: *
    3: *  =========== DOCUMENTATION ===========
    4: *
    5: * Online html documentation available at
    6: *            http://www.netlib.org/lapack/explore-html/
    7: *
    8: *> \htmlonly
    9: *> Download ZLAGS2 + dependencies
   10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlags2.f">
   11: *> [TGZ]</a>
   12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlags2.f">
   13: *> [ZIP]</a>
   14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlags2.f">
   15: *> [TXT]</a>
   16: *> \endhtmlonly
   17: *
   18: *  Definition:
   19: *  ===========
   20: *
   21: *       SUBROUTINE ZLAGS2( UPPER, A1, A2, A3, B1, B2, B3, CSU, SNU, CSV,
   22: *                          SNV, CSQ, SNQ )
   23: *
   24: *       .. Scalar Arguments ..
   25: *       LOGICAL            UPPER
   26: *       DOUBLE PRECISION   A1, A3, B1, B3, CSQ, CSU, CSV
   27: *       COMPLEX*16         A2, B2, SNQ, SNU, SNV
   28: *       ..
   29: *
   30: *
   31: *> \par Purpose:
   32: *  =============
   33: *>
   34: *> \verbatim
   35: *>
   36: *> ZLAGS2 computes 2-by-2 unitary matrices U, V and Q, such
   37: *> that if ( UPPER ) then
   38: *>
   39: *>           U**H *A*Q = U**H *( A1 A2 )*Q = ( x  0  )
   40: *>                             ( 0  A3 )     ( x  x  )
   41: *> and
   42: *>           V**H*B*Q = V**H *( B1 B2 )*Q = ( x  0  )
   43: *>                            ( 0  B3 )     ( x  x  )
   44: *>
   45: *> or if ( .NOT.UPPER ) then
   46: *>
   47: *>           U**H *A*Q = U**H *( A1 0  )*Q = ( x  x  )
   48: *>                             ( A2 A3 )     ( 0  x  )
   49: *> and
   50: *>           V**H *B*Q = V**H *( B1 0  )*Q = ( x  x  )
   51: *>                             ( B2 B3 )     ( 0  x  )
   52: *> where
   53: *>
   54: *>   U = (   CSU    SNU ), V = (  CSV    SNV ),
   55: *>       ( -SNU**H  CSU )      ( -SNV**H CSV )
   56: *>
   57: *>   Q = (   CSQ    SNQ )
   58: *>       ( -SNQ**H  CSQ )
   59: *>
   60: *> The rows of the transformed A and B are parallel. Moreover, if the
   61: *> input 2-by-2 matrix A is not zero, then the transformed (1,1) entry
   62: *> of A is not zero. If the input matrices A and B are both not zero,
   63: *> then the transformed (2,2) element of B is not zero, except when the
   64: *> first rows of input A and B are parallel and the second rows are
   65: *> zero.
   66: *> \endverbatim
   67: *
   68: *  Arguments:
   69: *  ==========
   70: *
   71: *> \param[in] UPPER
   72: *> \verbatim
   73: *>          UPPER is LOGICAL
   74: *>          = .TRUE.: the input matrices A and B are upper triangular.
   75: *>          = .FALSE.: the input matrices A and B are lower triangular.
   76: *> \endverbatim
   77: *>
   78: *> \param[in] A1
   79: *> \verbatim
   80: *>          A1 is DOUBLE PRECISION
   81: *> \endverbatim
   82: *>
   83: *> \param[in] A2
   84: *> \verbatim
   85: *>          A2 is COMPLEX*16
   86: *> \endverbatim
   87: *>
   88: *> \param[in] A3
   89: *> \verbatim
   90: *>          A3 is DOUBLE PRECISION
   91: *>          On entry, A1, A2 and A3 are elements of the input 2-by-2
   92: *>          upper (lower) triangular matrix A.
   93: *> \endverbatim
   94: *>
   95: *> \param[in] B1
   96: *> \verbatim
   97: *>          B1 is DOUBLE PRECISION
   98: *> \endverbatim
   99: *>
  100: *> \param[in] B2
  101: *> \verbatim
  102: *>          B2 is COMPLEX*16
  103: *> \endverbatim
  104: *>
  105: *> \param[in] B3
  106: *> \verbatim
  107: *>          B3 is DOUBLE PRECISION
  108: *>          On entry, B1, B2 and B3 are elements of the input 2-by-2
  109: *>          upper (lower) triangular matrix B.
  110: *> \endverbatim
  111: *>
  112: *> \param[out] CSU
  113: *> \verbatim
  114: *>          CSU is DOUBLE PRECISION
  115: *> \endverbatim
  116: *>
  117: *> \param[out] SNU
  118: *> \verbatim
  119: *>          SNU is COMPLEX*16
  120: *>          The desired unitary matrix U.
  121: *> \endverbatim
  122: *>
  123: *> \param[out] CSV
  124: *> \verbatim
  125: *>          CSV is DOUBLE PRECISION
  126: *> \endverbatim
  127: *>
  128: *> \param[out] SNV
  129: *> \verbatim
  130: *>          SNV is COMPLEX*16
  131: *>          The desired unitary matrix V.
  132: *> \endverbatim
  133: *>
  134: *> \param[out] CSQ
  135: *> \verbatim
  136: *>          CSQ is DOUBLE PRECISION
  137: *> \endverbatim
  138: *>
  139: *> \param[out] SNQ
  140: *> \verbatim
  141: *>          SNQ is COMPLEX*16
  142: *>          The desired unitary matrix Q.
  143: *> \endverbatim
  144: *
  145: *  Authors:
  146: *  ========
  147: *
  148: *> \author Univ. of Tennessee
  149: *> \author Univ. of California Berkeley
  150: *> \author Univ. of Colorado Denver
  151: *> \author NAG Ltd.
  152: *
  153: *> \date December 2016
  154: *
  155: *> \ingroup complex16OTHERauxiliary
  156: *
  157: *  =====================================================================
  158:       SUBROUTINE ZLAGS2( UPPER, A1, A2, A3, B1, B2, B3, CSU, SNU, CSV,
  159:      $                   SNV, CSQ, SNQ )
  160: *
  161: *  -- LAPACK auxiliary routine (version 3.7.0) --
  162: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  163: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  164: *     December 2016
  165: *
  166: *     .. Scalar Arguments ..
  167:       LOGICAL            UPPER
  168:       DOUBLE PRECISION   A1, A3, B1, B3, CSQ, CSU, CSV
  169:       COMPLEX*16         A2, B2, SNQ, SNU, SNV
  170: *     ..
  171: *
  172: *  =====================================================================
  173: *
  174: *     .. Parameters ..
  175:       DOUBLE PRECISION   ZERO, ONE
  176:       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
  177: *     ..
  178: *     .. Local Scalars ..
  179:       DOUBLE PRECISION   A, AUA11, AUA12, AUA21, AUA22, AVB12, AVB11,
  180:      $                   AVB21, AVB22, CSL, CSR, D, FB, FC, S1, S2,
  181:      $                   SNL, SNR, UA11R, UA22R, VB11R, VB22R
  182:       COMPLEX*16         B, C, D1, R, T, UA11, UA12, UA21, UA22, VB11,
  183:      $                   VB12, VB21, VB22
  184: *     ..
  185: *     .. External Subroutines ..
  186:       EXTERNAL           DLASV2, ZLARTG
  187: *     ..
  188: *     .. Intrinsic Functions ..
  189:       INTRINSIC          ABS, DBLE, DCMPLX, DCONJG, DIMAG
  190: *     ..
  191: *     .. Statement Functions ..
  192:       DOUBLE PRECISION   ABS1
  193: *     ..
  194: *     .. Statement Function definitions ..
  195:       ABS1( T ) = ABS( DBLE( T ) ) + ABS( DIMAG( T ) )
  196: *     ..
  197: *     .. Executable Statements ..
  198: *
  199:       IF( UPPER ) THEN
  200: *
  201: *        Input matrices A and B are upper triangular matrices
  202: *
  203: *        Form matrix C = A*adj(B) = ( a b )
  204: *                                   ( 0 d )
  205: *
  206:          A = A1*B3
  207:          D = A3*B1
  208:          B = A2*B1 - A1*B2
  209:          FB = ABS( B )
  210: *
  211: *        Transform complex 2-by-2 matrix C to real matrix by unitary
  212: *        diagonal matrix diag(1,D1).
  213: *
  214:          D1 = ONE
  215:          IF( FB.NE.ZERO )
  216:      $      D1 = B / FB
  217: *
  218: *        The SVD of real 2 by 2 triangular C
  219: *
  220: *         ( CSL -SNL )*( A B )*(  CSR  SNR ) = ( R 0 )
  221: *         ( SNL  CSL ) ( 0 D ) ( -SNR  CSR )   ( 0 T )
  222: *
  223:          CALL DLASV2( A, FB, D, S1, S2, SNR, CSR, SNL, CSL )
  224: *
  225:          IF( ABS( CSL ).GE.ABS( SNL ) .OR. ABS( CSR ).GE.ABS( SNR ) )
  226:      $        THEN
  227: *
  228: *           Compute the (1,1) and (1,2) elements of U**H *A and V**H *B,
  229: *           and (1,2) element of |U|**H *|A| and |V|**H *|B|.
  230: *
  231:             UA11R = CSL*A1
  232:             UA12 = CSL*A2 + D1*SNL*A3
  233: *
  234:             VB11R = CSR*B1
  235:             VB12 = CSR*B2 + D1*SNR*B3
  236: *
  237:             AUA12 = ABS( CSL )*ABS1( A2 ) + ABS( SNL )*ABS( A3 )
  238:             AVB12 = ABS( CSR )*ABS1( B2 ) + ABS( SNR )*ABS( B3 )
  239: *
  240: *           zero (1,2) elements of U**H *A and V**H *B
  241: *
  242:             IF( ( ABS( UA11R )+ABS1( UA12 ) ).EQ.ZERO ) THEN
  243:                CALL ZLARTG( -DCMPLX( VB11R ), DCONJG( VB12 ), CSQ, SNQ,
  244:      $                      R )
  245:             ELSE IF( ( ABS( VB11R )+ABS1( VB12 ) ).EQ.ZERO ) THEN
  246:                CALL ZLARTG( -DCMPLX( UA11R ), DCONJG( UA12 ), CSQ, SNQ,
  247:      $                      R )
  248:             ELSE IF( AUA12 / ( ABS( UA11R )+ABS1( UA12 ) ).LE.AVB12 /
  249:      $               ( ABS( VB11R )+ABS1( VB12 ) ) ) THEN
  250:                CALL ZLARTG( -DCMPLX( UA11R ), DCONJG( UA12 ), CSQ, SNQ,
  251:      $                      R )
  252:             ELSE
  253:                CALL ZLARTG( -DCMPLX( VB11R ), DCONJG( VB12 ), CSQ, SNQ,
  254:      $                      R )
  255:             END IF
  256: *
  257:             CSU = CSL
  258:             SNU = -D1*SNL
  259:             CSV = CSR
  260:             SNV = -D1*SNR
  261: *
  262:          ELSE
  263: *
  264: *           Compute the (2,1) and (2,2) elements of U**H *A and V**H *B,
  265: *           and (2,2) element of |U|**H *|A| and |V|**H *|B|.
  266: *
  267:             UA21 = -DCONJG( D1 )*SNL*A1
  268:             UA22 = -DCONJG( D1 )*SNL*A2 + CSL*A3
  269: *
  270:             VB21 = -DCONJG( D1 )*SNR*B1
  271:             VB22 = -DCONJG( D1 )*SNR*B2 + CSR*B3
  272: *
  273:             AUA22 = ABS( SNL )*ABS1( A2 ) + ABS( CSL )*ABS( A3 )
  274:             AVB22 = ABS( SNR )*ABS1( B2 ) + ABS( CSR )*ABS( B3 )
  275: *
  276: *           zero (2,2) elements of U**H *A and V**H *B, and then swap.
  277: *
  278:             IF( ( ABS1( UA21 )+ABS1( UA22 ) ).EQ.ZERO ) THEN
  279:                CALL ZLARTG( -DCONJG( VB21 ), DCONJG( VB22 ), CSQ, SNQ,
  280:      $                      R )
  281:             ELSE IF( ( ABS1( VB21 )+ABS( VB22 ) ).EQ.ZERO ) THEN
  282:                CALL ZLARTG( -DCONJG( UA21 ), DCONJG( UA22 ), CSQ, SNQ,
  283:      $                      R )
  284:             ELSE IF( AUA22 / ( ABS1( UA21 )+ABS1( UA22 ) ).LE.AVB22 /
  285:      $               ( ABS1( VB21 )+ABS1( VB22 ) ) ) THEN
  286:                CALL ZLARTG( -DCONJG( UA21 ), DCONJG( UA22 ), CSQ, SNQ,
  287:      $                      R )
  288:             ELSE
  289:                CALL ZLARTG( -DCONJG( VB21 ), DCONJG( VB22 ), CSQ, SNQ,
  290:      $                      R )
  291:             END IF
  292: *
  293:             CSU = SNL
  294:             SNU = D1*CSL
  295:             CSV = SNR
  296:             SNV = D1*CSR
  297: *
  298:          END IF
  299: *
  300:       ELSE
  301: *
  302: *        Input matrices A and B are lower triangular matrices
  303: *
  304: *        Form matrix C = A*adj(B) = ( a 0 )
  305: *                                   ( c d )
  306: *
  307:          A = A1*B3
  308:          D = A3*B1
  309:          C = A2*B3 - A3*B2
  310:          FC = ABS( C )
  311: *
  312: *        Transform complex 2-by-2 matrix C to real matrix by unitary
  313: *        diagonal matrix diag(d1,1).
  314: *
  315:          D1 = ONE
  316:          IF( FC.NE.ZERO )
  317:      $      D1 = C / FC
  318: *
  319: *        The SVD of real 2 by 2 triangular C
  320: *
  321: *         ( CSL -SNL )*( A 0 )*(  CSR  SNR ) = ( R 0 )
  322: *         ( SNL  CSL ) ( C D ) ( -SNR  CSR )   ( 0 T )
  323: *
  324:          CALL DLASV2( A, FC, D, S1, S2, SNR, CSR, SNL, CSL )
  325: *
  326:          IF( ABS( CSR ).GE.ABS( SNR ) .OR. ABS( CSL ).GE.ABS( SNL ) )
  327:      $        THEN
  328: *
  329: *           Compute the (2,1) and (2,2) elements of U**H *A and V**H *B,
  330: *           and (2,1) element of |U|**H *|A| and |V|**H *|B|.
  331: *
  332:             UA21 = -D1*SNR*A1 + CSR*A2
  333:             UA22R = CSR*A3
  334: *
  335:             VB21 = -D1*SNL*B1 + CSL*B2
  336:             VB22R = CSL*B3
  337: *
  338:             AUA21 = ABS( SNR )*ABS( A1 ) + ABS( CSR )*ABS1( A2 )
  339:             AVB21 = ABS( SNL )*ABS( B1 ) + ABS( CSL )*ABS1( B2 )
  340: *
  341: *           zero (2,1) elements of U**H *A and V**H *B.
  342: *
  343:             IF( ( ABS1( UA21 )+ABS( UA22R ) ).EQ.ZERO ) THEN
  344:                CALL ZLARTG( DCMPLX( VB22R ), VB21, CSQ, SNQ, R )
  345:             ELSE IF( ( ABS1( VB21 )+ABS( VB22R ) ).EQ.ZERO ) THEN
  346:                CALL ZLARTG( DCMPLX( UA22R ), UA21, CSQ, SNQ, R )
  347:             ELSE IF( AUA21 / ( ABS1( UA21 )+ABS( UA22R ) ).LE.AVB21 /
  348:      $               ( ABS1( VB21 )+ABS( VB22R ) ) ) THEN
  349:                CALL ZLARTG( DCMPLX( UA22R ), UA21, CSQ, SNQ, R )
  350:             ELSE
  351:                CALL ZLARTG( DCMPLX( VB22R ), VB21, CSQ, SNQ, R )
  352:             END IF
  353: *
  354:             CSU = CSR
  355:             SNU = -DCONJG( D1 )*SNR
  356:             CSV = CSL
  357:             SNV = -DCONJG( D1 )*SNL
  358: *
  359:          ELSE
  360: *
  361: *           Compute the (1,1) and (1,2) elements of U**H *A and V**H *B,
  362: *           and (1,1) element of |U|**H *|A| and |V|**H *|B|.
  363: *
  364:             UA11 = CSR*A1 + DCONJG( D1 )*SNR*A2
  365:             UA12 = DCONJG( D1 )*SNR*A3
  366: *
  367:             VB11 = CSL*B1 + DCONJG( D1 )*SNL*B2
  368:             VB12 = DCONJG( D1 )*SNL*B3
  369: *
  370:             AUA11 = ABS( CSR )*ABS( A1 ) + ABS( SNR )*ABS1( A2 )
  371:             AVB11 = ABS( CSL )*ABS( B1 ) + ABS( SNL )*ABS1( B2 )
  372: *
  373: *           zero (1,1) elements of U**H *A and V**H *B, and then swap.
  374: *
  375:             IF( ( ABS1( UA11 )+ABS1( UA12 ) ).EQ.ZERO ) THEN
  376:                CALL ZLARTG( VB12, VB11, CSQ, SNQ, R )
  377:             ELSE IF( ( ABS1( VB11 )+ABS1( VB12 ) ).EQ.ZERO ) THEN
  378:                CALL ZLARTG( UA12, UA11, CSQ, SNQ, R )
  379:             ELSE IF( AUA11 / ( ABS1( UA11 )+ABS1( UA12 ) ).LE.AVB11 /
  380:      $               ( ABS1( VB11 )+ABS1( VB12 ) ) ) THEN
  381:                CALL ZLARTG( UA12, UA11, CSQ, SNQ, R )
  382:             ELSE
  383:                CALL ZLARTG( VB12, VB11, CSQ, SNQ, R )
  384:             END IF
  385: *
  386:             CSU = SNR
  387:             SNU = DCONJG( D1 )*CSR
  388:             CSV = SNL
  389:             SNV = DCONJG( D1 )*CSL
  390: *
  391:          END IF
  392: *
  393:       END IF
  394: *
  395:       RETURN
  396: *
  397: *     End of ZLAGS2
  398: *
  399:       END

CVSweb interface <joel.bertrand@systella.fr>