version 1.1.1.1, 2010/01/26 15:22:45
|
version 1.17, 2023/08/07 08:38:43
|
Line 1
|
Line 1
|
SUBROUTINE DROTMG(DD1,DD2,DX1,DY1,DPARAM) |
*> \brief \b DROTMG |
* .. Scalar Arguments .. |
|
DOUBLE PRECISION DD1,DD2,DX1,DY1 |
|
* .. |
|
* .. Array Arguments .. |
|
DOUBLE PRECISION DPARAM(5) |
|
* .. |
|
* |
* |
* Purpose |
* =========== DOCUMENTATION =========== |
* ======= |
|
* |
* |
* CONSTRUCT THE MODIFIED GIVENS TRANSFORMATION MATRIX H WHICH ZEROS |
* Online html documentation available at |
* THE SECOND COMPONENT OF THE 2-VECTOR (DSQRT(DD1)*DX1,DSQRT(DD2)* |
* http://www.netlib.org/lapack/explore-html/ |
* DY2)**T. |
|
* WITH DPARAM(1)=DFLAG, H HAS ONE OF THE FOLLOWING FORMS.. |
|
* |
* |
* DFLAG=-1.D0 DFLAG=0.D0 DFLAG=1.D0 DFLAG=-2.D0 |
* Definition: |
|
* =========== |
|
* |
|
* SUBROUTINE DROTMG(DD1,DD2,DX1,DY1,DPARAM) |
|
* |
|
* .. Scalar Arguments .. |
|
* DOUBLE PRECISION DD1,DD2,DX1,DY1 |
|
* .. |
|
* .. Array Arguments .. |
|
* DOUBLE PRECISION DPARAM(5) |
|
* .. |
|
* |
|
* |
|
*> \par Purpose: |
|
* ============= |
|
*> |
|
*> \verbatim |
|
*> |
|
*> CONSTRUCT THE MODIFIED GIVENS TRANSFORMATION MATRIX H WHICH ZEROS |
|
*> THE SECOND COMPONENT OF THE 2-VECTOR (DSQRT(DD1)*DX1,DSQRT(DD2)*> DY2)**T. |
|
*> WITH DPARAM(1)=DFLAG, H HAS ONE OF THE FOLLOWING FORMS.. |
|
*> |
|
*> DFLAG=-1.D0 DFLAG=0.D0 DFLAG=1.D0 DFLAG=-2.D0 |
|
*> |
|
*> (DH11 DH12) (1.D0 DH12) (DH11 1.D0) (1.D0 0.D0) |
|
*> H=( ) ( ) ( ) ( ) |
|
*> (DH21 DH22), (DH21 1.D0), (-1.D0 DH22), (0.D0 1.D0). |
|
*> LOCATIONS 2-4 OF DPARAM CONTAIN DH11, DH21, DH12, AND DH22 |
|
*> RESPECTIVELY. (VALUES OF 1.D0, -1.D0, OR 0.D0 IMPLIED BY THE |
|
*> VALUE OF DPARAM(1) ARE NOT STORED IN DPARAM.) |
|
*> |
|
*> THE VALUES OF GAMSQ AND RGAMSQ SET IN THE DATA STATEMENT MAY BE |
|
*> INEXACT. THIS IS OK AS THEY ARE ONLY USED FOR TESTING THE SIZE |
|
*> OF DD1 AND DD2. ALL ACTUAL SCALING OF DATA IS DONE USING GAM. |
|
*> |
|
*> \endverbatim |
|
* |
|
* Arguments: |
|
* ========== |
|
* |
|
*> \param[in,out] DD1 |
|
*> \verbatim |
|
*> DD1 is DOUBLE PRECISION |
|
*> \endverbatim |
|
*> |
|
*> \param[in,out] DD2 |
|
*> \verbatim |
|
*> DD2 is DOUBLE PRECISION |
|
*> \endverbatim |
|
*> |
|
*> \param[in,out] DX1 |
|
*> \verbatim |
|
*> DX1 is DOUBLE PRECISION |
|
*> \endverbatim |
|
*> |
|
*> \param[in] DY1 |
|
*> \verbatim |
|
*> DY1 is DOUBLE PRECISION |
|
*> \endverbatim |
|
*> |
|
*> \param[out] DPARAM |
|
*> \verbatim |
|
*> DPARAM is DOUBLE PRECISION array, dimension (5) |
|
*> DPARAM(1)=DFLAG |
|
*> DPARAM(2)=DH11 |
|
*> DPARAM(3)=DH21 |
|
*> DPARAM(4)=DH12 |
|
*> DPARAM(5)=DH22 |
|
*> \endverbatim |
|
* |
|
* Authors: |
|
* ======== |
|
* |
|
*> \author Univ. of Tennessee |
|
*> \author Univ. of California Berkeley |
|
*> \author Univ. of Colorado Denver |
|
*> \author NAG Ltd. |
* |
* |
* (DH11 DH12) (1.D0 DH12) (DH11 1.D0) (1.D0 0.D0) |
*> \ingroup double_blas_level1 |
* H=( ) ( ) ( ) ( ) |
|
* (DH21 DH22), (DH21 1.D0), (-1.D0 DH22), (0.D0 1.D0). |
|
* LOCATIONS 2-4 OF DPARAM CONTAIN DH11, DH21, DH12, AND DH22 |
|
* RESPECTIVELY. (VALUES OF 1.D0, -1.D0, OR 0.D0 IMPLIED BY THE |
|
* VALUE OF DPARAM(1) ARE NOT STORED IN DPARAM.) |
|
* |
* |
* THE VALUES OF GAMSQ AND RGAMSQ SET IN THE DATA STATEMENT MAY BE |
* ===================================================================== |
* INEXACT. THIS IS OK AS THEY ARE ONLY USED FOR TESTING THE SIZE |
SUBROUTINE DROTMG(DD1,DD2,DX1,DY1,DPARAM) |
* OF DD1 AND DD2. ALL ACTUAL SCALING OF DATA IS DONE USING GAM. |
|
* |
|
* |
|
* Arguments |
|
* ========= |
|
* |
|
* DD1 (input/output) DOUBLE PRECISION |
|
* |
|
* DD2 (input/output) DOUBLE PRECISION |
|
* |
|
* DX1 (input/output) DOUBLE PRECISION |
|
* |
* |
* DY1 (input) DOUBLE PRECISION |
* -- Reference BLAS level1 routine -- |
|
* -- Reference BLAS is a software package provided by Univ. of Tennessee, -- |
|
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- |
* |
* |
* DPARAM (input/output) DOUBLE PRECISION array, dimension 5 |
* .. Scalar Arguments .. |
* DPARAM(1)=DFLAG |
DOUBLE PRECISION DD1,DD2,DX1,DY1 |
* DPARAM(2)=DH11 |
* .. |
* DPARAM(3)=DH21 |
* .. Array Arguments .. |
* DPARAM(4)=DH12 |
DOUBLE PRECISION DPARAM(5) |
* DPARAM(5)=DH22 |
* .. |
* |
* |
* ===================================================================== |
* ===================================================================== |
* |
* |
* .. Local Scalars .. |
* .. Local Scalars .. |
DOUBLE PRECISION DFLAG,DH11,DH12,DH21,DH22,DP1,DP2,DQ1,DQ2,DTEMP, |
DOUBLE PRECISION DFLAG,DH11,DH12,DH21,DH22,DP1,DP2,DQ1,DQ2,DTEMP, |
+ DU,GAM,GAMSQ,ONE,RGAMSQ,TWO,ZERO |
$ DU,GAM,GAMSQ,ONE,RGAMSQ,TWO,ZERO |
INTEGER IGO |
|
* .. |
* .. |
* .. Intrinsic Functions .. |
* .. Intrinsic Functions .. |
INTRINSIC DABS |
INTRINSIC DABS |
Line 62
|
Line 114
|
DATA GAM,GAMSQ,RGAMSQ/4096.D0,16777216.D0,5.9604645D-8/ |
DATA GAM,GAMSQ,RGAMSQ/4096.D0,16777216.D0,5.9604645D-8/ |
* .. |
* .. |
|
|
IF (.NOT.DD1.LT.ZERO) GO TO 10 |
IF (DD1.LT.ZERO) THEN |
* GO ZERO-H-D-AND-DX1.. |
* GO ZERO-H-D-AND-DX1.. |
GO TO 60 |
DFLAG = -ONE |
10 CONTINUE |
DH11 = ZERO |
* CASE-DD1-NONNEGATIVE |
DH12 = ZERO |
DP2 = DD2*DY1 |
DH21 = ZERO |
IF (.NOT.DP2.EQ.ZERO) GO TO 20 |
DH22 = ZERO |
DFLAG = -TWO |
* |
GO TO 260 |
DD1 = ZERO |
* REGULAR-CASE.. |
DD2 = ZERO |
20 CONTINUE |
DX1 = ZERO |
DP1 = DD1*DX1 |
ELSE |
DQ2 = DP2*DY1 |
* CASE-DD1-NONNEGATIVE |
DQ1 = DP1*DX1 |
DP2 = DD2*DY1 |
* |
IF (DP2.EQ.ZERO) THEN |
IF (.NOT.DABS(DQ1).GT.DABS(DQ2)) GO TO 40 |
DFLAG = -TWO |
DH21 = -DY1/DX1 |
DPARAM(1) = DFLAG |
DH12 = DP2/DP1 |
RETURN |
* |
END IF |
DU = ONE - DH12*DH21 |
* REGULAR-CASE.. |
* |
DP1 = DD1*DX1 |
IF (.NOT.DU.LE.ZERO) GO TO 30 |
DQ2 = DP2*DY1 |
* GO ZERO-H-D-AND-DX1.. |
DQ1 = DP1*DX1 |
GO TO 60 |
* |
30 CONTINUE |
IF (DABS(DQ1).GT.DABS(DQ2)) THEN |
DFLAG = ZERO |
DH21 = -DY1/DX1 |
DD1 = DD1/DU |
DH12 = DP2/DP1 |
DD2 = DD2/DU |
* |
DX1 = DX1*DU |
DU = ONE - DH12*DH21 |
* GO SCALE-CHECK.. |
* |
GO TO 100 |
IF (DU.GT.ZERO) THEN |
40 CONTINUE |
DFLAG = ZERO |
IF (.NOT.DQ2.LT.ZERO) GO TO 50 |
DD1 = DD1/DU |
* GO ZERO-H-D-AND-DX1.. |
DD2 = DD2/DU |
GO TO 60 |
DX1 = DX1*DU |
50 CONTINUE |
ELSE |
DFLAG = ONE |
* This code path if here for safety. We do not expect this |
DH11 = DP1/DP2 |
* condition to ever hold except in edge cases with rounding |
DH22 = DX1/DY1 |
* errors. See DOI: 10.1145/355841.355847 |
DU = ONE + DH11*DH22 |
DFLAG = -ONE |
DTEMP = DD2/DU |
DH11 = ZERO |
DD2 = DD1/DU |
DH12 = ZERO |
DD1 = DTEMP |
DH21 = ZERO |
DX1 = DY1*DU |
DH22 = ZERO |
* GO SCALE-CHECK |
* |
GO TO 100 |
DD1 = ZERO |
* PROCEDURE..ZERO-H-D-AND-DX1.. |
DD2 = ZERO |
60 CONTINUE |
DX1 = ZERO |
DFLAG = -ONE |
END IF |
DH11 = ZERO |
ELSE |
DH12 = ZERO |
|
DH21 = ZERO |
IF (DQ2.LT.ZERO) THEN |
DH22 = ZERO |
* GO ZERO-H-D-AND-DX1.. |
* |
DFLAG = -ONE |
DD1 = ZERO |
DH11 = ZERO |
DD2 = ZERO |
DH12 = ZERO |
DX1 = ZERO |
DH21 = ZERO |
* RETURN.. |
DH22 = ZERO |
GO TO 220 |
* |
* PROCEDURE..FIX-H.. |
DD1 = ZERO |
70 CONTINUE |
DD2 = ZERO |
IF (.NOT.DFLAG.GE.ZERO) GO TO 90 |
DX1 = ZERO |
* |
ELSE |
IF (.NOT.DFLAG.EQ.ZERO) GO TO 80 |
DFLAG = ONE |
DH11 = ONE |
DH11 = DP1/DP2 |
DH22 = ONE |
DH22 = DX1/DY1 |
DFLAG = -ONE |
DU = ONE + DH11*DH22 |
GO TO 90 |
DTEMP = DD2/DU |
80 CONTINUE |
DD2 = DD1/DU |
DH21 = -ONE |
DD1 = DTEMP |
DH12 = ONE |
DX1 = DY1*DU |
DFLAG = -ONE |
END IF |
90 CONTINUE |
END IF |
GO TO IGO(120,150,180,210) |
|
* PROCEDURE..SCALE-CHECK |
* PROCEDURE..SCALE-CHECK |
100 CONTINUE |
IF (DD1.NE.ZERO) THEN |
110 CONTINUE |
DO WHILE ((DD1.LE.RGAMSQ) .OR. (DD1.GE.GAMSQ)) |
IF (.NOT.DD1.LE.RGAMSQ) GO TO 130 |
IF (DFLAG.EQ.ZERO) THEN |
IF (DD1.EQ.ZERO) GO TO 160 |
DH11 = ONE |
ASSIGN 120 TO IGO |
DH22 = ONE |
* FIX-H.. |
DFLAG = -ONE |
GO TO 70 |
ELSE |
120 CONTINUE |
DH21 = -ONE |
DD1 = DD1*GAM**2 |
DH12 = ONE |
DX1 = DX1/GAM |
DFLAG = -ONE |
DH11 = DH11/GAM |
END IF |
DH12 = DH12/GAM |
IF (DD1.LE.RGAMSQ) THEN |
GO TO 110 |
DD1 = DD1*GAM**2 |
130 CONTINUE |
DX1 = DX1/GAM |
140 CONTINUE |
DH11 = DH11/GAM |
IF (.NOT.DD1.GE.GAMSQ) GO TO 160 |
DH12 = DH12/GAM |
ASSIGN 150 TO IGO |
ELSE |
* FIX-H.. |
DD1 = DD1/GAM**2 |
GO TO 70 |
DX1 = DX1*GAM |
150 CONTINUE |
DH11 = DH11*GAM |
DD1 = DD1/GAM**2 |
DH12 = DH12*GAM |
DX1 = DX1*GAM |
END IF |
DH11 = DH11*GAM |
ENDDO |
DH12 = DH12*GAM |
END IF |
GO TO 140 |
|
160 CONTINUE |
IF (DD2.NE.ZERO) THEN |
170 CONTINUE |
DO WHILE ( (DABS(DD2).LE.RGAMSQ) .OR. (DABS(DD2).GE.GAMSQ) ) |
IF (.NOT.DABS(DD2).LE.RGAMSQ) GO TO 190 |
IF (DFLAG.EQ.ZERO) THEN |
IF (DD2.EQ.ZERO) GO TO 220 |
DH11 = ONE |
ASSIGN 180 TO IGO |
DH22 = ONE |
* FIX-H.. |
DFLAG = -ONE |
GO TO 70 |
ELSE |
180 CONTINUE |
DH21 = -ONE |
DD2 = DD2*GAM**2 |
DH12 = ONE |
DH21 = DH21/GAM |
DFLAG = -ONE |
DH22 = DH22/GAM |
END IF |
GO TO 170 |
IF (DABS(DD2).LE.RGAMSQ) THEN |
190 CONTINUE |
DD2 = DD2*GAM**2 |
200 CONTINUE |
DH21 = DH21/GAM |
IF (.NOT.DABS(DD2).GE.GAMSQ) GO TO 220 |
DH22 = DH22/GAM |
ASSIGN 210 TO IGO |
ELSE |
* FIX-H.. |
DD2 = DD2/GAM**2 |
GO TO 70 |
DH21 = DH21*GAM |
210 CONTINUE |
DH22 = DH22*GAM |
DD2 = DD2/GAM**2 |
END IF |
DH21 = DH21*GAM |
END DO |
DH22 = DH22*GAM |
END IF |
GO TO 200 |
|
220 CONTINUE |
END IF |
IF (DFLAG) 250,230,240 |
|
230 CONTINUE |
IF (DFLAG.LT.ZERO) THEN |
DPARAM(3) = DH21 |
DPARAM(2) = DH11 |
DPARAM(4) = DH12 |
DPARAM(3) = DH21 |
GO TO 260 |
DPARAM(4) = DH12 |
240 CONTINUE |
DPARAM(5) = DH22 |
DPARAM(2) = DH11 |
ELSE IF (DFLAG.EQ.ZERO) THEN |
DPARAM(5) = DH22 |
DPARAM(3) = DH21 |
GO TO 260 |
DPARAM(4) = DH12 |
250 CONTINUE |
ELSE |
DPARAM(2) = DH11 |
DPARAM(2) = DH11 |
DPARAM(3) = DH21 |
DPARAM(5) = DH22 |
DPARAM(4) = DH12 |
END IF |
DPARAM(5) = DH22 |
|
260 CONTINUE |
|
DPARAM(1) = DFLAG |
DPARAM(1) = DFLAG |
RETURN |
RETURN |
|
* |
|
* End of DROTMG |
|
* |
END |
END |