File:  [local] / rpl / lapack / lapack / zlags2.f
Revision 1.18: download - view: text, annotated - select for diffs - revision graph
Mon Aug 7 08:39:29 2023 UTC (8 months, 3 weeks ago) by bertrand
Branches: MAIN
CVS tags: rpl-4_1_35, rpl-4_1_34, HEAD
Première mise à jour de lapack et blas.

    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: *> \ingroup complex16OTHERauxiliary
  154: *
  155: *  =====================================================================
  156:       SUBROUTINE ZLAGS2( UPPER, A1, A2, A3, B1, B2, B3, CSU, SNU, CSV,
  157:      $                   SNV, CSQ, SNQ )
  158: *
  159: *  -- LAPACK auxiliary routine --
  160: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  161: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  162: *
  163: *     .. Scalar Arguments ..
  164:       LOGICAL            UPPER
  165:       DOUBLE PRECISION   A1, A3, B1, B3, CSQ, CSU, CSV
  166:       COMPLEX*16         A2, B2, SNQ, SNU, SNV
  167: *     ..
  168: *
  169: *  =====================================================================
  170: *
  171: *     .. Parameters ..
  172:       DOUBLE PRECISION   ZERO, ONE
  173:       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
  174: *     ..
  175: *     .. Local Scalars ..
  176:       DOUBLE PRECISION   A, AUA11, AUA12, AUA21, AUA22, AVB12, AVB11,
  177:      $                   AVB21, AVB22, CSL, CSR, D, FB, FC, S1, S2,
  178:      $                   SNL, SNR, UA11R, UA22R, VB11R, VB22R
  179:       COMPLEX*16         B, C, D1, R, T, UA11, UA12, UA21, UA22, VB11,
  180:      $                   VB12, VB21, VB22
  181: *     ..
  182: *     .. External Subroutines ..
  183:       EXTERNAL           DLASV2, ZLARTG
  184: *     ..
  185: *     .. Intrinsic Functions ..
  186:       INTRINSIC          ABS, DBLE, DCMPLX, DCONJG, DIMAG
  187: *     ..
  188: *     .. Statement Functions ..
  189:       DOUBLE PRECISION   ABS1
  190: *     ..
  191: *     .. Statement Function definitions ..
  192:       ABS1( T ) = ABS( DBLE( T ) ) + ABS( DIMAG( T ) )
  193: *     ..
  194: *     .. Executable Statements ..
  195: *
  196:       IF( UPPER ) THEN
  197: *
  198: *        Input matrices A and B are upper triangular matrices
  199: *
  200: *        Form matrix C = A*adj(B) = ( a b )
  201: *                                   ( 0 d )
  202: *
  203:          A = A1*B3
  204:          D = A3*B1
  205:          B = A2*B1 - A1*B2
  206:          FB = ABS( B )
  207: *
  208: *        Transform complex 2-by-2 matrix C to real matrix by unitary
  209: *        diagonal matrix diag(1,D1).
  210: *
  211:          D1 = ONE
  212:          IF( FB.NE.ZERO )
  213:      $      D1 = B / FB
  214: *
  215: *        The SVD of real 2 by 2 triangular C
  216: *
  217: *         ( CSL -SNL )*( A B )*(  CSR  SNR ) = ( R 0 )
  218: *         ( SNL  CSL ) ( 0 D ) ( -SNR  CSR )   ( 0 T )
  219: *
  220:          CALL DLASV2( A, FB, D, S1, S2, SNR, CSR, SNL, CSL )
  221: *
  222:          IF( ABS( CSL ).GE.ABS( SNL ) .OR. ABS( CSR ).GE.ABS( SNR ) )
  223:      $        THEN
  224: *
  225: *           Compute the (1,1) and (1,2) elements of U**H *A and V**H *B,
  226: *           and (1,2) element of |U|**H *|A| and |V|**H *|B|.
  227: *
  228:             UA11R = CSL*A1
  229:             UA12 = CSL*A2 + D1*SNL*A3
  230: *
  231:             VB11R = CSR*B1
  232:             VB12 = CSR*B2 + D1*SNR*B3
  233: *
  234:             AUA12 = ABS( CSL )*ABS1( A2 ) + ABS( SNL )*ABS( A3 )
  235:             AVB12 = ABS( CSR )*ABS1( B2 ) + ABS( SNR )*ABS( B3 )
  236: *
  237: *           zero (1,2) elements of U**H *A and V**H *B
  238: *
  239:             IF( ( ABS( UA11R )+ABS1( UA12 ) ).EQ.ZERO ) THEN
  240:                CALL ZLARTG( -DCMPLX( VB11R ), DCONJG( VB12 ), CSQ, SNQ,
  241:      $                      R )
  242:             ELSE IF( ( ABS( VB11R )+ABS1( VB12 ) ).EQ.ZERO ) THEN
  243:                CALL ZLARTG( -DCMPLX( UA11R ), DCONJG( UA12 ), CSQ, SNQ,
  244:      $                      R )
  245:             ELSE IF( AUA12 / ( ABS( UA11R )+ABS1( UA12 ) ).LE.AVB12 /
  246:      $               ( ABS( VB11R )+ABS1( VB12 ) ) ) THEN
  247:                CALL ZLARTG( -DCMPLX( UA11R ), DCONJG( UA12 ), CSQ, SNQ,
  248:      $                      R )
  249:             ELSE
  250:                CALL ZLARTG( -DCMPLX( VB11R ), DCONJG( VB12 ), CSQ, SNQ,
  251:      $                      R )
  252:             END IF
  253: *
  254:             CSU = CSL
  255:             SNU = -D1*SNL
  256:             CSV = CSR
  257:             SNV = -D1*SNR
  258: *
  259:          ELSE
  260: *
  261: *           Compute the (2,1) and (2,2) elements of U**H *A and V**H *B,
  262: *           and (2,2) element of |U|**H *|A| and |V|**H *|B|.
  263: *
  264:             UA21 = -DCONJG( D1 )*SNL*A1
  265:             UA22 = -DCONJG( D1 )*SNL*A2 + CSL*A3
  266: *
  267:             VB21 = -DCONJG( D1 )*SNR*B1
  268:             VB22 = -DCONJG( D1 )*SNR*B2 + CSR*B3
  269: *
  270:             AUA22 = ABS( SNL )*ABS1( A2 ) + ABS( CSL )*ABS( A3 )
  271:             AVB22 = ABS( SNR )*ABS1( B2 ) + ABS( CSR )*ABS( B3 )
  272: *
  273: *           zero (2,2) elements of U**H *A and V**H *B, and then swap.
  274: *
  275:             IF( ( ABS1( UA21 )+ABS1( UA22 ) ).EQ.ZERO ) THEN
  276:                CALL ZLARTG( -DCONJG( VB21 ), DCONJG( VB22 ), CSQ, SNQ,
  277:      $                      R )
  278:             ELSE IF( ( ABS1( VB21 )+ABS( VB22 ) ).EQ.ZERO ) THEN
  279:                CALL ZLARTG( -DCONJG( UA21 ), DCONJG( UA22 ), CSQ, SNQ,
  280:      $                      R )
  281:             ELSE IF( AUA22 / ( ABS1( UA21 )+ABS1( UA22 ) ).LE.AVB22 /
  282:      $               ( ABS1( VB21 )+ABS1( VB22 ) ) ) THEN
  283:                CALL ZLARTG( -DCONJG( UA21 ), DCONJG( UA22 ), CSQ, SNQ,
  284:      $                      R )
  285:             ELSE
  286:                CALL ZLARTG( -DCONJG( VB21 ), DCONJG( VB22 ), CSQ, SNQ,
  287:      $                      R )
  288:             END IF
  289: *
  290:             CSU = SNL
  291:             SNU = D1*CSL
  292:             CSV = SNR
  293:             SNV = D1*CSR
  294: *
  295:          END IF
  296: *
  297:       ELSE
  298: *
  299: *        Input matrices A and B are lower triangular matrices
  300: *
  301: *        Form matrix C = A*adj(B) = ( a 0 )
  302: *                                   ( c d )
  303: *
  304:          A = A1*B3
  305:          D = A3*B1
  306:          C = A2*B3 - A3*B2
  307:          FC = ABS( C )
  308: *
  309: *        Transform complex 2-by-2 matrix C to real matrix by unitary
  310: *        diagonal matrix diag(d1,1).
  311: *
  312:          D1 = ONE
  313:          IF( FC.NE.ZERO )
  314:      $      D1 = C / FC
  315: *
  316: *        The SVD of real 2 by 2 triangular C
  317: *
  318: *         ( CSL -SNL )*( A 0 )*(  CSR  SNR ) = ( R 0 )
  319: *         ( SNL  CSL ) ( C D ) ( -SNR  CSR )   ( 0 T )
  320: *
  321:          CALL DLASV2( A, FC, D, S1, S2, SNR, CSR, SNL, CSL )
  322: *
  323:          IF( ABS( CSR ).GE.ABS( SNR ) .OR. ABS( CSL ).GE.ABS( SNL ) )
  324:      $        THEN
  325: *
  326: *           Compute the (2,1) and (2,2) elements of U**H *A and V**H *B,
  327: *           and (2,1) element of |U|**H *|A| and |V|**H *|B|.
  328: *
  329:             UA21 = -D1*SNR*A1 + CSR*A2
  330:             UA22R = CSR*A3
  331: *
  332:             VB21 = -D1*SNL*B1 + CSL*B2
  333:             VB22R = CSL*B3
  334: *
  335:             AUA21 = ABS( SNR )*ABS( A1 ) + ABS( CSR )*ABS1( A2 )
  336:             AVB21 = ABS( SNL )*ABS( B1 ) + ABS( CSL )*ABS1( B2 )
  337: *
  338: *           zero (2,1) elements of U**H *A and V**H *B.
  339: *
  340:             IF( ( ABS1( UA21 )+ABS( UA22R ) ).EQ.ZERO ) THEN
  341:                CALL ZLARTG( DCMPLX( VB22R ), VB21, CSQ, SNQ, R )
  342:             ELSE IF( ( ABS1( VB21 )+ABS( VB22R ) ).EQ.ZERO ) THEN
  343:                CALL ZLARTG( DCMPLX( UA22R ), UA21, CSQ, SNQ, R )
  344:             ELSE IF( AUA21 / ( ABS1( UA21 )+ABS( UA22R ) ).LE.AVB21 /
  345:      $               ( ABS1( VB21 )+ABS( VB22R ) ) ) THEN
  346:                CALL ZLARTG( DCMPLX( UA22R ), UA21, CSQ, SNQ, R )
  347:             ELSE
  348:                CALL ZLARTG( DCMPLX( VB22R ), VB21, CSQ, SNQ, R )
  349:             END IF
  350: *
  351:             CSU = CSR
  352:             SNU = -DCONJG( D1 )*SNR
  353:             CSV = CSL
  354:             SNV = -DCONJG( D1 )*SNL
  355: *
  356:          ELSE
  357: *
  358: *           Compute the (1,1) and (1,2) elements of U**H *A and V**H *B,
  359: *           and (1,1) element of |U|**H *|A| and |V|**H *|B|.
  360: *
  361:             UA11 = CSR*A1 + DCONJG( D1 )*SNR*A2
  362:             UA12 = DCONJG( D1 )*SNR*A3
  363: *
  364:             VB11 = CSL*B1 + DCONJG( D1 )*SNL*B2
  365:             VB12 = DCONJG( D1 )*SNL*B3
  366: *
  367:             AUA11 = ABS( CSR )*ABS( A1 ) + ABS( SNR )*ABS1( A2 )
  368:             AVB11 = ABS( CSL )*ABS( B1 ) + ABS( SNL )*ABS1( B2 )
  369: *
  370: *           zero (1,1) elements of U**H *A and V**H *B, and then swap.
  371: *
  372:             IF( ( ABS1( UA11 )+ABS1( UA12 ) ).EQ.ZERO ) THEN
  373:                CALL ZLARTG( VB12, VB11, CSQ, SNQ, R )
  374:             ELSE IF( ( ABS1( VB11 )+ABS1( VB12 ) ).EQ.ZERO ) THEN
  375:                CALL ZLARTG( UA12, UA11, CSQ, SNQ, R )
  376:             ELSE IF( AUA11 / ( ABS1( UA11 )+ABS1( UA12 ) ).LE.AVB11 /
  377:      $               ( ABS1( VB11 )+ABS1( VB12 ) ) ) THEN
  378:                CALL ZLARTG( UA12, UA11, CSQ, SNQ, R )
  379:             ELSE
  380:                CALL ZLARTG( VB12, VB11, CSQ, SNQ, R )
  381:             END IF
  382: *
  383:             CSU = SNR
  384:             SNU = DCONJG( D1 )*CSR
  385:             CSV = SNL
  386:             SNV = DCONJG( D1 )*CSL
  387: *
  388:          END IF
  389: *
  390:       END IF
  391: *
  392:       RETURN
  393: *
  394: *     End of ZLAGS2
  395: *
  396:       END

CVSweb interface <joel.bertrand@systella.fr>