blas/f2c/drotmg.c

/* drotmg.f -- translated by f2c (version 20100827).
   You must link the resulting object file with libf2c:
        on Microsoft Windows system, link with libf2c.lib;
        on Linux or Unix systems, link with .../path/to/libf2c.a -lm
        or, if you install libf2c.a in a standard place, with -lf2c -lm
        -- in that order, at the end of the command line, as in
                cc *.o -lf2c -lm
        Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,

                http://www.netlib.org/f2c/libf2c.zip
*/

#include "datatypes.h"

/* Subroutine */ void drotmg_(doublereal *dd1, doublereal *dd2, doublereal *dx1, doublereal *dy1, doublereal *dparam) {
  /* Initialized data */

  static doublereal zero = 0.;
  static doublereal one = 1.;
  static doublereal two = 2.;
  static doublereal gam = 4096.;
  static doublereal gamsq = 16777216.;
  static doublereal rgamsq = 5.9604645e-8;

  /* Format strings */
  static char fmt_120[] = "";
  static char fmt_150[] = "";
  static char fmt_180[] = "";
  static char fmt_210[] = "";

  /* System generated locals */
  doublereal d__1;

  /* Local variables */
  doublereal du, dp1, dp2, dq1, dq2, dh11, dh12, dh21, dh22;
  integer igo;
  doublereal dflag, dtemp;

  /* Assigned format variables */
  static char *igo_fmt;
  (void)igo_fmt;

  /*     .. Scalar Arguments .. */
  /*     .. */
  /*     .. Array Arguments .. */
  /*     .. */

  /*  Purpose */
  /*  ======= */

  /*     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. */

  /*  Arguments */
  /*  ========= */

  /*  DD1    (input/output) DOUBLE PRECISION */

  /*  DD2    (input/output) DOUBLE PRECISION */

  /*  DX1    (input/output) DOUBLE PRECISION */

  /*  DY1    (input) DOUBLE PRECISION */

  /*  DPARAM (input/output)  DOUBLE PRECISION array, dimension 5 */
  /*     DPARAM(1)=DFLAG */
  /*     DPARAM(2)=DH11 */
  /*     DPARAM(3)=DH21 */
  /*     DPARAM(4)=DH12 */
  /*     DPARAM(5)=DH22 */

  /*  ===================================================================== */

  /*     .. Local Scalars .. */
  /*     .. */
  /*     .. Intrinsic Functions .. */
  /*     .. */
  /*     .. Data statements .. */

  /* Parameter adjustments */
  --dparam;

  /* Function Body */
  /*     .. */
  if (!(*dd1 < zero)) {
    goto L10;
  }
  /*       GO ZERO-H-D-AND-DX1.. */
  goto L60;
L10:
  /*     CASE-DD1-NONNEGATIVE */
  dp2 = *dd2 * *dy1;
  if (!(dp2 == zero)) {
    goto L20;
  }
  dflag = -two;
  goto L260;
/*     REGULAR-CASE.. */
L20:
  dp1 = *dd1 * *dx1;
  dq2 = dp2 * *dy1;
  dq1 = dp1 * *dx1;

  if (!(abs(dq1) > abs(dq2))) {
    goto L40;
  }
  dh21 = -(*dy1) / *dx1;
  dh12 = dp2 / dp1;

  du = one - dh12 * dh21;

  if (!(du <= zero)) {
    goto L30;
  }
  /*         GO ZERO-H-D-AND-DX1.. */
  goto L60;
L30:
  dflag = zero;
  *dd1 /= du;
  *dd2 /= du;
  *dx1 *= du;
  /*         GO SCALE-CHECK.. */
  goto L100;
L40:
  if (!(dq2 < zero)) {
    goto L50;
  }
  /*         GO ZERO-H-D-AND-DX1.. */
  goto L60;
L50:
  dflag = one;
  dh11 = dp1 / dp2;
  dh22 = *dx1 / *dy1;
  du = one + dh11 * dh22;
  dtemp = *dd2 / du;
  *dd2 = *dd1 / du;
  *dd1 = dtemp;
  *dx1 = *dy1 * du;
  /*         GO SCALE-CHECK */
  goto L100;
/*     PROCEDURE..ZERO-H-D-AND-DX1.. */
L60:
  dflag = -one;
  dh11 = zero;
  dh12 = zero;
  dh21 = zero;
  dh22 = zero;

  *dd1 = zero;
  *dd2 = zero;
  *dx1 = zero;
  /*         RETURN.. */
  goto L220;
/*     PROCEDURE..FIX-H.. */
L70:
  if (!(dflag >= zero)) {
    goto L90;
  }

  if (!(dflag == zero)) {
    goto L80;
  }
  dh11 = one;
  dh22 = one;
  dflag = -one;
  goto L90;
L80:
  dh21 = -one;
  dh12 = one;
  dflag = -one;
L90:
  switch (igo) {
    case 0:
      goto L120;
    case 1:
      goto L150;
    case 2:
      goto L180;
    case 3:
      goto L210;
  }
/*     PROCEDURE..SCALE-CHECK */
L100:
L110:
  if (!(*dd1 <= rgamsq)) {
    goto L130;
  }
  if (*dd1 == zero) {
    goto L160;
  }
  igo = 0;
  igo_fmt = fmt_120;
  /*              FIX-H.. */
  goto L70;
L120:
  /* Computing 2nd power */
  d__1 = gam;
  *dd1 *= d__1 * d__1;
  *dx1 /= gam;
  dh11 /= gam;
  dh12 /= gam;
  goto L110;
L130:
L140:
  if (!(*dd1 >= gamsq)) {
    goto L160;
  }
  igo = 1;
  igo_fmt = fmt_150;
  /*              FIX-H.. */
  goto L70;
L150:
  /* Computing 2nd power */
  d__1 = gam;
  *dd1 /= d__1 * d__1;
  *dx1 *= gam;
  dh11 *= gam;
  dh12 *= gam;
  goto L140;
L160:
L170:
  if (!(abs(*dd2) <= rgamsq)) {
    goto L190;
  }
  if (*dd2 == zero) {
    goto L220;
  }
  igo = 2;
  igo_fmt = fmt_180;
  /*              FIX-H.. */
  goto L70;
L180:
  /* Computing 2nd power */
  d__1 = gam;
  *dd2 *= d__1 * d__1;
  dh21 /= gam;
  dh22 /= gam;
  goto L170;
L190:
L200:
  if (!(abs(*dd2) >= gamsq)) {
    goto L220;
  }
  igo = 3;
  igo_fmt = fmt_210;
  /*              FIX-H.. */
  goto L70;
L210:
  /* Computing 2nd power */
  d__1 = gam;
  *dd2 /= d__1 * d__1;
  dh21 *= gam;
  dh22 *= gam;
  goto L200;
L220:
  if (dflag < 0.) {
    goto L250;
  } else if (dflag == 0) {
    goto L230;
  } else {
    goto L240;
  }
L230:
  dparam[3] = dh21;
  dparam[4] = dh12;
  goto L260;
L240:
  dparam[2] = dh11;
  dparam[5] = dh22;
  goto L260;
L250:
  dparam[2] = dh11;
  dparam[3] = dh21;
  dparam[4] = dh12;
  dparam[5] = dh22;
L260:
  dparam[1] = dflag;
} /* drotmg_ */