www.pudn.com > sources_4610.rar > NUMcblas.c
/* NUMcblas.c
-- LAPACK helper routines -- Univ. of Tennessee, Univ. of
California Berkeley, NAG Ltd., Courant Institute, Argonne National Lab,
and Rice University October 31, 1999 -- translated by f2c (version
19990503)
Adapted by David Weenink 20021201.
*/
/*
djmw 20020813 GPL header
djmw 20030205 Latest modification
*/
/* #include "blaswrap.h" */
#include "melder.h"
#include "NUMcblas.h"
#include "NUMf2c.h"
#include "NUM2.h"
#define MAX(m,n) ((m) > (n) ? (m) : (n))
#define MIN(m,n) ((m) < (n) ? (m) : (n))
static int dlamc1_ (long *beta, long *t, long *rnd, long *ieee1);
static int dlamc2_ (long *beta, long *t, long *rnd, double *eps, long *emin, double *rmin, long *emax,
double *rmax);
static double dlamc3_ (double *, double *);
static int dlamc4_ (long *emin, double *start, long *base);
static int dlamc5_ (long *beta, long *p, long *emin, long *ieee, long *emax, double *rmax);
int NUMblas_daxpy (long *n, double *da, double *dx, long *incx, double *dy, long *incy)
{
/* System generated locals */
long i__1;
/* Local variables */
static long i__, m, ix, iy, mp1;
--dy;
--dx;
/* Function Body */
if (*n <= 0)
{
return 0;
}
if (*da == 0.)
{
return 0;
}
if (*incx == 1 && *incy == 1)
{
goto L20;
}
/* code for unequal increments or equal increments not equal to 1 */
ix = 1;
iy = 1;
if (*incx < 0)
{
ix = (-(*n) + 1) * *incx + 1;
}
if (*incy < 0)
{
iy = (-(*n) + 1) * *incy + 1;
}
i__1 = *n;
for (i__ = 1; i__ <= i__1; ++i__)
{
dy[iy] += *da * dx[ix];
ix += *incx;
iy += *incy;
/* L10: */
}
return 0;
/* code for both increments equal to 1 clean-up loop */
L20:
m = *n % 4;
if (m == 0)
{
goto L40;
}
i__1 = m;
for (i__ = 1; i__ <= i__1; ++i__)
{
dy[i__] += *da * dx[i__];
/* L30: */
}
if (*n < 4)
{
return 0;
}
L40:
mp1 = m + 1;
i__1 = *n;
for (i__ = mp1; i__ <= i__1; i__ += 4)
{
dy[i__] += *da * dx[i__];
dy[i__ + 1] += *da * dx[i__ + 1];
dy[i__ + 2] += *da * dx[i__ + 2];
dy[i__ + 3] += *da * dx[i__ + 3];
/* L50: */
}
return 0;
} /* NUMblas_daxpy */
int NUMblas_dcopy (long *n, double *dx, long *incx, double *dy, long *incy)
{
/* System generated locals */
long i__1;
/* Local variables */
static long i__, m, ix, iy, mp1;
--dy;
--dx;
/* Function Body */
if (*n <= 0)
{
return 0;
}
if (*incx == 1 && *incy == 1)
{
goto L20;
}
/* code for unequal increments or equal increments not equal to 1 */
ix = 1;
iy = 1;
if (*incx < 0)
{
ix = (-(*n) + 1) * *incx + 1;
}
if (*incy < 0)
{
iy = (-(*n) + 1) * *incy + 1;
}
i__1 = *n;
for (i__ = 1; i__ <= i__1; ++i__)
{
dy[iy] = dx[ix];
ix += *incx;
iy += *incy;
/* L10: */
}
return 0;
/* code for both increments equal to 1 clean-up loop */
L20:
m = *n % 7;
if (m == 0)
{
goto L40;
}
i__1 = m;
for (i__ = 1; i__ <= i__1; ++i__)
{
dy[i__] = dx[i__];
/* L30: */
}
if (*n < 7)
{
return 0;
}
L40:
mp1 = m + 1;
i__1 = *n;
for (i__ = mp1; i__ <= i__1; i__ += 7)
{
dy[i__] = dx[i__];
dy[i__ + 1] = dx[i__ + 1];
dy[i__ + 2] = dx[i__ + 2];
dy[i__ + 3] = dx[i__ + 3];
dy[i__ + 4] = dx[i__ + 4];
dy[i__ + 5] = dx[i__ + 5];
dy[i__ + 6] = dx[i__ + 6];
/* L50: */
}
return 0;
} /* NUMblas_dcopy */
double NUMblas_ddot (long *n, double *dx, long *incx, double *dy, long *incy)
{
/* System generated locals */
long i__1;
double ret_val;
/* Local variables */
static long i__, m;
static double dtemp;
static long ix, iy, mp1;
/* Parameter adjustments */
--dy;
--dx;
/* Function Body */
ret_val = 0.;
dtemp = 0.;
if (*n <= 0)
{
return ret_val;
}
if (*incx == 1 && *incy == 1)
{
goto L20;
}
/* code for unequal increments or equal increments not equal to 1 */
ix = 1;
iy = 1;
if (*incx < 0)
{
ix = (-(*n) + 1) * *incx + 1;
}
if (*incy < 0)
{
iy = (-(*n) + 1) * *incy + 1;
}
i__1 = *n;
for (i__ = 1; i__ <= i__1; ++i__)
{
dtemp += dx[ix] * dy[iy];
ix += *incx;
iy += *incy;
/* L10: */
}
ret_val = dtemp;
return ret_val;
/* code for both increments equal to 1 clean-up loop */
L20:
m = *n % 5;
if (m == 0)
{
goto L40;
}
i__1 = m;
for (i__ = 1; i__ <= i__1; ++i__)
{
dtemp += dx[i__] * dy[i__];
/* L30: */
}
if (*n < 5)
{
goto L60;
}
L40:
mp1 = m + 1;
i__1 = *n;
for (i__ = mp1; i__ <= i__1; i__ += 5)
{
dtemp =
dtemp + dx[i__] * dy[i__] + dx[i__ + 1] * dy[i__ + 1] + dx[i__ + 2] * dy[i__ + 2] + dx[i__ +
3] * dy[i__ + 3] + dx[i__ + 4] * dy[i__ + 4];
/* L50: */
}
L60:
ret_val = dtemp;
return ret_val;
} /* NUMblas_ddot */
int NUMblas_dgemm (char *transa, char *transb, long *m, long *n, long *k, double *alpha, double *a, long *lda,
double *b, long *ldb, double *beta, double *c__, long *ldc)
{
/* System generated locals */
long a_dim1, a_offset, b_dim1, b_offset, c_dim1, c_offset, i__1, i__2, i__3;
/* Local variables */
static long info;
static long nota, notb;
static double temp;
static long i__, j, l, ncola;
static long nrowa, nrowb;
#define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
#define b_ref(a_1,a_2) b[(a_2)*b_dim1 + a_1]
#define c___ref(a_1,a_2) c__[(a_2)*c_dim1 + a_1]
/*
Set NOTA and NOTB as true if A and B respectively are not transposed
and set NROWA, NCOLA and NROWB as the number of rows and columns of A
and the number of rows of B respectively. Parameter adjustments */
a_dim1 = *lda;
a_offset = 1 + a_dim1 * 1;
a -= a_offset;
b_dim1 = *ldb;
b_offset = 1 + b_dim1 * 1;
b -= b_offset;
c_dim1 = *ldc;
c_offset = 1 + c_dim1 * 1;
c__ -= c_offset;
/* Function Body */
nota = lsame_ (transa, "N");
notb = lsame_ (transb, "N");
if (nota)
{
nrowa = *m;
ncola = *k;
}
else
{
nrowa = *k;
ncola = *m;
}
if (notb)
{
nrowb = *k;
}
else
{
nrowb = *n;
}
/* Test the input parameters. */
info = 0;
if (!nota && !lsame_ (transa, "C") && !lsame_ (transa, "T"))
{
info = 1;
}
else if (!notb && !lsame_ (transb, "C") && !lsame_ (transb, "T"))
{
info = 2;
}
else if (*m < 0)
{
info = 3;
}
else if (*n < 0)
{
info = 4;
}
else if (*k < 0)
{
info = 5;
}
else if (*lda < MAX (1, nrowa))
{
info = 8;
}
else if (*ldb < MAX (1, nrowb))
{
info = 10;
}
else if (*ldc < MAX (1, *m))
{
info = 13;
}
if (info != 0)
{
xerbla_ ("DGEMM ", &info);
return 0;
}
/* Quick return if possible. */
if (*m == 0 || *n == 0 || (*alpha == 0. || *k == 0) && *beta == 1.)
{
return 0;
}
/* And if alpha.eq.zero. */
if (*alpha == 0.)
{
if (*beta == 0.)
{
i__1 = *n;
for (j = 1; j <= i__1; ++j)
{
i__2 = *m;
for (i__ = 1; i__ <= i__2; ++i__)
{
c___ref (i__, j) = 0.;
/* L10: */
}
/* L20: */
}
}
else
{
i__1 = *n;
for (j = 1; j <= i__1; ++j)
{
i__2 = *m;
for (i__ = 1; i__ <= i__2; ++i__)
{
c___ref (i__, j) = *beta * c___ref (i__, j);
/* L30: */
}
/* L40: */
}
}
return 0;
}
/* Start the operations. */
if (notb)
{
if (nota)
{
/* Form C := alpha*A*B + beta*C. */
i__1 = *n;
for (j = 1; j <= i__1; ++j)
{
if (*beta == 0.)
{
i__2 = *m;
for (i__ = 1; i__ <= i__2; ++i__)
{
c___ref (i__, j) = 0.;
/* L50: */
}
}
else if (*beta != 1.)
{
i__2 = *m;
for (i__ = 1; i__ <= i__2; ++i__)
{
c___ref (i__, j) = *beta * c___ref (i__, j);
/* L60: */
}
}
i__2 = *k;
for (l = 1; l <= i__2; ++l)
{
if (b_ref (l, j) != 0.)
{
temp = *alpha * b_ref (l, j);
i__3 = *m;
for (i__ = 1; i__ <= i__3; ++i__)
{
c___ref (i__, j) = c___ref (i__, j) + temp * a_ref (i__, l);
/* L70: */
}
}
/* L80: */
}
/* L90: */
}
}
else
{
/* Form C := alpha*A'*B + beta*C */
i__1 = *n;
for (j = 1; j <= i__1; ++j)
{
i__2 = *m;
for (i__ = 1; i__ <= i__2; ++i__)
{
temp = 0.;
i__3 = *k;
for (l = 1; l <= i__3; ++l)
{
temp += a_ref (l, i__) * b_ref (l, j);
/* L100: */
}
if (*beta == 0.)
{
c___ref (i__, j) = *alpha * temp;
}
else
{
c___ref (i__, j) = *alpha * temp + *beta * c___ref (i__, j);
}
/* L110: */
}
/* L120: */
}
}
}
else
{
if (nota)
{
/* Form C := alpha*A*B' + beta*C */
i__1 = *n;
for (j = 1; j <= i__1; ++j)
{
if (*beta == 0.)
{
i__2 = *m;
for (i__ = 1; i__ <= i__2; ++i__)
{
c___ref (i__, j) = 0.;
/* L130: */
}
}
else if (*beta != 1.)
{
i__2 = *m;
for (i__ = 1; i__ <= i__2; ++i__)
{
c___ref (i__, j) = *beta * c___ref (i__, j);
/* L140: */
}
}
i__2 = *k;
for (l = 1; l <= i__2; ++l)
{
if (b_ref (j, l) != 0.)
{
temp = *alpha * b_ref (j, l);
i__3 = *m;
for (i__ = 1; i__ <= i__3; ++i__)
{
c___ref (i__, j) = c___ref (i__, j) + temp * a_ref (i__, l);
/* L150: */
}
}
/* L160: */
}
/* L170: */
}
}
else
{
/* Form C := alpha*A'*B' + beta*C */
i__1 = *n;
for (j = 1; j <= i__1; ++j)
{
i__2 = *m;
for (i__ = 1; i__ <= i__2; ++i__)
{
temp = 0.;
i__3 = *k;
for (l = 1; l <= i__3; ++l)
{
temp += a_ref (l, i__) * b_ref (j, l);
/* L180: */
}
if (*beta == 0.)
{
c___ref (i__, j) = *alpha * temp;
}
else
{
c___ref (i__, j) = *alpha * temp + *beta * c___ref (i__, j);
}
/* L190: */
}
/* L200: */
}
}
}
return 0;
/* End of DGEMM . */
} /* NUMblas_dgemm */
#undef c___ref
#undef b_ref
#undef a_ref
int NUMblas_dger (long *m, long *n, double *alpha, double *x, long *incx, double *y, long *incy, double *a,
long *lda)
{
/* System generated locals */
long a_dim1, a_offset, i__1, i__2;
/* Local variables */
static long info;
static double temp;
static long i__, j, ix, jy, kx;
#define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
/* Test the input parameters. Parameter adjustments */
--x;
--y;
a_dim1 = *lda;
a_offset = 1 + a_dim1 * 1;
a -= a_offset;
/* Function Body */
info = 0;
if (*m < 0)
{
info = 1;
}
else if (*n < 0)
{
info = 2;
}
else if (*incx == 0)
{
info = 5;
}
else if (*incy == 0)
{
info = 7;
}
else if (*lda < MAX (1, *m))
{
info = 9;
}
if (info != 0)
{
xerbla_ ("DGER ", &info);
return 0;
}
/* Quick return if possible. */
if (*m == 0 || *n == 0 || *alpha == 0.)
{
return 0;
}
/* Start the operations. In this version the elements of A are accessed
sequentially with one pass through A. */
if (*incy > 0)
{
jy = 1;
}
else
{
jy = 1 - (*n - 1) * *incy;
}
if (*incx == 1)
{
i__1 = *n;
for (j = 1; j <= i__1; ++j)
{
if (y[jy] != 0.)
{
temp = *alpha * y[jy];
i__2 = *m;
for (i__ = 1; i__ <= i__2; ++i__)
{
a_ref (i__, j) = a_ref (i__, j) + x[i__] * temp;
/* L10: */
}
}
jy += *incy;
/* L20: */
}
}
else
{
if (*incx > 0)
{
kx = 1;
}
else
{
kx = 1 - (*m - 1) * *incx;
}
i__1 = *n;
for (j = 1; j <= i__1; ++j)
{
if (y[jy] != 0.)
{
temp = *alpha * y[jy];
ix = kx;
i__2 = *m;
for (i__ = 1; i__ <= i__2; ++i__)
{
a_ref (i__, j) = a_ref (i__, j) + x[ix] * temp;
ix += *incx;
/* L30: */
}
}
jy += *incy;
/* L40: */
}
}
return 0;
} /* NUMblas_dger */
#undef a_ref
int NUMblas_dgemv (char *trans, long *m, long *n, double *alpha, double *a, long *lda, double *x, long *incx,
double *beta, double *y, long *incy)
{
/* System generated locals */
long a_dim1, a_offset, i__1, i__2;
/* Local variables */
static long info;
static double temp;
static long lenx, leny, i__, j;
static long ix, iy, jx, jy, kx, ky;
#define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
/* Parameter adjustments */
a_dim1 = *lda;
a_offset = 1 + a_dim1 * 1;
a -= a_offset;
--x;
--y;
/* Function Body */
info = 0;
if (!lsame_ (trans, "N") && !lsame_ (trans, "T") && !lsame_ (trans, "C"))
{
info = 1;
}
else if (*m < 0)
{
info = 2;
}
else if (*n < 0)
{
info = 3;
}
else if (*lda < MAX (1, *m))
{
info = 6;
}
else if (*incx == 0)
{
info = 8;
}
else if (*incy == 0)
{
info = 11;
}
if (info != 0)
{
xerbla_ ("DGEMV ", &info);
return 0;
}
/* Quick return if possible. */
if (*m == 0 || *n == 0 || *alpha == 0. && *beta == 1.)
{
return 0;
}
/* Set LENX and LENY, the lengths of the vectors x and y, and set up the
start points in X and Y. */
if (lsame_ (trans, "N"))
{
lenx = *n;
leny = *m;
}
else
{
lenx = *m;
leny = *n;
}
if (*incx > 0)
{
kx = 1;
}
else
{
kx = 1 - (lenx - 1) * *incx;
}
if (*incy > 0)
{
ky = 1;
}
else
{
ky = 1 - (leny - 1) * *incy;
}
/* Start the operations. In this version the elements of A are accessed
sequentially with one pass through A. First form y := beta*y. */
if (*beta != 1.)
{
if (*incy == 1)
{
if (*beta == 0.)
{
i__1 = leny;
for (i__ = 1; i__ <= i__1; ++i__)
{
y[i__] = 0.;
/* L10: */
}
}
else
{
i__1 = leny;
for (i__ = 1; i__ <= i__1; ++i__)
{
y[i__] = *beta * y[i__];
/* L20: */
}
}
}
else
{
iy = ky;
if (*beta == 0.)
{
i__1 = leny;
for (i__ = 1; i__ <= i__1; ++i__)
{
y[iy] = 0.;
iy += *incy;
/* L30: */
}
}
else
{
i__1 = leny;
for (i__ = 1; i__ <= i__1; ++i__)
{
y[iy] = *beta * y[iy];
iy += *incy;
/* L40: */
}
}
}
}
if (*alpha == 0.)
{
return 0;
}
if (lsame_ (trans, "N"))
{
/* Form y := alpha*A*x + y. */
jx = kx;
if (*incy == 1)
{
i__1 = *n;
for (j = 1; j <= i__1; ++j)
{
if (x[jx] != 0.)
{
temp = *alpha * x[jx];
i__2 = *m;
for (i__ = 1; i__ <= i__2; ++i__)
{
y[i__] += temp * a_ref (i__, j);
/* L50: */
}
}
jx += *incx;
/* L60: */
}
}
else
{
i__1 = *n;
for (j = 1; j <= i__1; ++j)
{
if (x[jx] != 0.)
{
temp = *alpha * x[jx];
iy = ky;
i__2 = *m;
for (i__ = 1; i__ <= i__2; ++i__)
{
y[iy] += temp * a_ref (i__, j);
iy += *incy;
/* L70: */
}
}
jx += *incx;
/* L80: */
}
}
}
else
{
/* Form y := alpha*A'*x + y. */
jy = ky;
if (*incx == 1)
{
i__1 = *n;
for (j = 1; j <= i__1; ++j)
{
temp = 0.;
i__2 = *m;
for (i__ = 1; i__ <= i__2; ++i__)
{
temp += a_ref (i__, j) * x[i__];
/* L90: */
}
y[jy] += *alpha * temp;
jy += *incy;
/* L100: */
}
}
else
{
i__1 = *n;
for (j = 1; j <= i__1; ++j)
{
temp = 0.;
ix = kx;
i__2 = *m;
for (i__ = 1; i__ <= i__2; ++i__)
{
temp += a_ref (i__, j) * x[ix];
ix += *incx;
/* L110: */
}
y[jy] += *alpha * temp;
jy += *incy;
/* L120: */
}
}
}
return 0;
} /* NUMblas_dgemv */
#undef a_ref
double NUMblas_dlamch (char *cmach)
{
/* Initialized data */
static long first = TRUE;
/* System generated locals */
long i__1;
double ret_val;
/* Builtin functions */
/* Local variables */
static double base;
static long beta;
static double emin, prec, emax;
static long imin, imax;
static long lrnd;
static double rmin, rmax, t, rmach;
static double smal, sfmin;
static long it;
static double rnd, eps;
if (first)
{
first = FALSE;
dlamc2_ (&beta, &it, &lrnd, &eps, &imin, &rmin, &imax, &rmax);
base = (double) beta;
t = (double) it;
if (lrnd)
{
rnd = 1.;
i__1 = 1 - it;
eps = pow_di (&base, &i__1) / 2;
}
else
{
rnd = 0.;
i__1 = 1 - it;
eps = pow_di (&base, &i__1);
}
prec = eps * base;
emin = (double) imin;
emax = (double) imax;
sfmin = rmin;
smal = 1. / rmax;
if (smal >= sfmin)
{
/* Use smal plus a bit, to avoid the possibility of rounding
causing overflow when computing 1/sfmin. */
sfmin = smal * (eps + 1.);
}
}
if (lsame_ (cmach, "E"))
{
rmach = eps;
}
else if (lsame_ (cmach, "S"))
{
rmach = sfmin;
}
else if (lsame_ (cmach, "B"))
{
rmach = base;
}
else if (lsame_ (cmach, "P"))
{
rmach = prec;
}
else if (lsame_ (cmach, "N"))
{
rmach = t;
}
else if (lsame_ (cmach, "R"))
{
rmach = rnd;
}
else if (lsame_ (cmach, "M"))
{
rmach = emin;
}
else if (lsame_ (cmach, "U"))
{
rmach = rmin;
}
else if (lsame_ (cmach, "L"))
{
rmach = emax;
}
else if (lsame_ (cmach, "O"))
{
rmach = rmax;
}
ret_val = rmach;
return ret_val;
} /* NUMblas_dlamch */
static int dlamc1_ (long *beta, long *t, long *rnd, long *ieee1)
{
/* -- LAPACK auxiliary routine (version 3.0) -- Univ. of Tennessee, Univ.
of California Berkeley, NAG Ltd., Courant Institute, Argonne National
Lab, and Rice University October 31, 1992
Purpose =======
DLAMC1 determines the machine parameters given by BETA, T, RND, and
IEEE1.
Arguments =========
BETA (output) INTEGER The base of the machine.
T (output) INTEGER The number of ( BETA ) digits in the mantissa.
RND (output) LOGICAL Specifies whether proper rounding ( RND = .TRUE.
) or chopping ( RND = .FALSE. ) occurs in addition. This may not
be a reliable guide to the way in which the machine performs
its arithmetic.
IEEE1 (output) LOGICAL Specifies whether rounding appears to be done
in the IEEE 'round to nearest' style.
Further Details ===============
The routine is based on the routine ENVRON by Malcolm and incorporates
suggestions by Gentleman and Marovich. See
Malcolm M. A. (1972) Algorithms to reveal properties of floating-point
arithmetic. Comms. of the ACM, 15, 949-951.
Gentleman W. M. and Marovich S. B. (1974) More on algorithms that
reveal properties of floating point arithmetic units. Comms. of the
ACM, 17, 276-277.
===================================================================== */
/* Initialized data */
static long first = TRUE;
/* System generated locals */
double d__1, d__2;
/* Local variables */
static long lrnd;
static double a, b, c, f;
static long lbeta;
static double savec;
static long lieee1;
static double t1, t2;
static long lt;
static double one, qtr;
if (first)
{
first = FALSE;
one = 1.;
/* LBETA, LIEEE1, LT and LRND are the local values of BETA, IEEE1, T
and RND.
Throughout this routine we use the function DLAMC3 to ensure that
relevant values are stored and not held in registers, or are not
affected by optimizers. Compute a = 2.0**m with the smallest
positive integer m such that fl( a + 1.0 ) = a. */
a = 1.;
c = 1.;
/* + WHILE( C.EQ.ONE )LOOP */
L10:
if (c == one)
{
a *= 2;
c = dlamc3_ (&a, &one);
d__1 = -a;
c = dlamc3_ (&c, &d__1);
goto L10;
}
/* + END WHILE
Now compute b = 2.0**m with the smallest positive integer m such
that fl( a + b ) .gt. a. */
b = 1.;
c = dlamc3_ (&a, &b);
/* + WHILE( C.EQ.A )LOOP */
L20:
if (c == a)
{
b *= 2;
c = dlamc3_ (&a, &b);
goto L20;
}
/* + END WHILE
Now compute the base. a and c are neighbouring floating point
numbers in the interval ( beta**t, beta**( t + 1 ) ) and so their
difference is beta. Adding 0.25 to c is to ensure that it is
truncated to beta and not ( beta - 1 ). */
qtr = one / 4;
savec = c;
d__1 = -a;
c = dlamc3_ (&c, &d__1);
lbeta = (long) (c + qtr);
/* Now determine whether rounding or chopping occurs, by adding a bit
less than beta/2 and a bit more than beta/2 to a. */
b = (double) lbeta;
d__1 = b / 2;
d__2 = -b / 100;
f = dlamc3_ (&d__1, &d__2);
c = dlamc3_ (&f, &a);
if (c == a)
{
lrnd = TRUE;
}
else
{
lrnd = FALSE;
}
d__1 = b / 2;
d__2 = b / 100;
f = dlamc3_ (&d__1, &d__2);
c = dlamc3_ (&f, &a);
if (lrnd && c == a)
{
lrnd = FALSE;
}
/* Try and decide whether rounding is done in the IEEE 'round to
nearest' style. B/2 is half a unit in the last place of the two
numbers A and SAVEC. Furthermore, A is even, i.e. has last bit
zero, and SAVEC is odd. Thus adding B/2 to A should not change A,
but adding B/2 to SAVEC should change SAVEC. */
d__1 = b / 2;
t1 = dlamc3_ (&d__1, &a);
d__1 = b / 2;
t2 = dlamc3_ (&d__1, &savec);
lieee1 = t1 == a && t2 > savec && lrnd;
/* Now find the mantissa, t. It should be the integer part of log to
the base beta of a, however it is safer to determine t by
powering. So we find t as the smallest positive integer for which
fl( beta**t + 1.0 ) = 1.0. */
lt = 0;
a = 1.;
c = 1.;
/* + WHILE( C.EQ.ONE )LOOP */
L30:
if (c == one)
{
++lt;
a *= lbeta;
c = dlamc3_ (&a, &one);
d__1 = -a;
c = dlamc3_ (&c, &d__1);
goto L30;
}
/* + END WHILE */
}
*beta = lbeta;
*t = lt;
*rnd = lrnd;
*ieee1 = lieee1;
return 0;
} /* dlamc1_ */
static int dlamc2_ (long *beta, long *t, long *rnd, double *eps, long *emin, double *rmin, long *emax,
double *rmax)
{
/* -- LAPACK auxiliary routine (version 3.0) -- Univ. of Tennessee, Univ.
of California Berkeley, NAG Ltd., Courant Institute, Argonne National
Lab, and Rice University October 31, 1992
Purpose =======
DLAMC2 determines the machine parameters specified in its argument
list.
Arguments =========
BETA (output) INTEGER The base of the machine.
T (output) INTEGER The number of ( BETA ) digits in the mantissa.
RND (output) LOGICAL Specifies whether proper rounding ( RND = .TRUE.
) or chopping ( RND = .FALSE. ) occurs in addition. This may not
be a reliable guide to the way in which the machine performs
its arithmetic.
EPS (output) DOUBLE PRECISION The smallest positive number such that
fl( 1.0 - EPS ) .LT. 1.0, where fl denotes the computed value.
EMIN (output) INTEGER The minimum exponent before (gradual) underflow
occurs.
RMIN (output) DOUBLE PRECISION The smallest normalized number for the
machine, given by BASE**( EMIN - 1 ), where BASE is the floating point
value of BETA.
EMAX (output) INTEGER The maximum exponent before overflow occurs.
RMAX (output) DOUBLE PRECISION The largest positive number for the
machine, given by BASE**EMAX * ( 1 - EPS ), where BASE is the floating
point value of BETA.
Further Details ===============
The computation of EPS is based on a routine PARANOIA by W. Kahan of
the University of California at Berkeley.
===================================================================== */
/* Table of constant values */
/* Initialized data */
static long first = TRUE;
static long iwarn = FALSE;
/* System generated locals */
long i__1;
double d__1, d__2, d__3, d__4, d__5;
/* Builtin functions */
/* Local variables */
static long ieee;
static double half;
static long lrnd;
static double leps, zero, a, b, c;
static long i, lbeta;
static double rbase;
static long lemin, lemax, gnmin;
static double smal;
static long gpmin;
static double third, lrmin, lrmax, sixth;
static long lieee1;
static long lt, ngnmin, ngpmin;
static double one, two;
if (first)
{
first = FALSE;
zero = 0.;
one = 1.;
two = 2.;
/* LBETA, LT, LRND, LEPS, LEMIN and LRMIN are the local values of
BETA, T, RND, EPS, EMIN and RMIN. Throughout this routine we use
the function DLAMC3 to ensure that relevant values are stored and
not held in registers, or are not affected by optimizers. DLAMC1
returns the parameters LBETA, LT, LRND and LIEEE1. */
dlamc1_ (&lbeta, <, &lrnd, &lieee1);
/* Start to find EPS. */
b = (double) lbeta;
i__1 = -lt;
a = pow_di (&b, &i__1);
leps = a;
/* Try some tricks to see whether or not this is the correct EPS. */
b = two / 3;
half = one / 2;
d__1 = -half;
sixth = dlamc3_ (&b, &d__1);
third = dlamc3_ (&sixth, &sixth);
d__1 = -half;
b = dlamc3_ (&third, &d__1);
b = dlamc3_ (&b, &sixth);
b = fabs (b);
if (b < leps)
{
b = leps;
}
leps = 1.;
/* + WHILE( ( LEPS.GT.B ).AND.( B.GT.ZERO ) )LOOP */
L10:
if (leps > b && b > zero)
{
leps = b;
d__1 = half * leps;
/* Computing 5th power */
d__3 = two, d__4 = d__3, d__3 *= d__3;
/* Computing 2nd power */
d__5 = leps;
d__2 = d__4 * (d__3 * d__3) * (d__5 * d__5);
c = dlamc3_ (&d__1, &d__2);
d__1 = -c;
c = dlamc3_ (&half, &d__1);
b = dlamc3_ (&half, &c);
d__1 = -b;
c = dlamc3_ (&half, &d__1);
b = dlamc3_ (&half, &c);
goto L10;
}
/* + END WHILE */
if (a < leps)
{
leps = a;
}
/* Computation of EPS complete. Now find EMIN. Let A = + or - 1, and
+ or - (1 + BASE**(-3)). Keep dividing A by BETA until (gradual)
underflow occurs. This is detected when we cannot recover the
previous A. */
rbase = one / lbeta;
smal = one;
for (i = 1; i <= 3; ++i)
{
d__1 = smal * rbase;
smal = dlamc3_ (&d__1, &zero);
/* L20: */
}
a = dlamc3_ (&one, &smal);
dlamc4_ (&ngpmin, &one, &lbeta);
d__1 = -one;
dlamc4_ (&ngnmin, &d__1, &lbeta);
dlamc4_ (&gpmin, &a, &lbeta);
d__1 = -a;
dlamc4_ (&gnmin, &d__1, &lbeta);
ieee = FALSE;
if (ngpmin == ngnmin && gpmin == gnmin)
{
if (ngpmin == gpmin)
{
lemin = ngpmin;
/* ( Non twos-complement machines, no gradual underflow;
e.g., VAX ) */
}
else if (gpmin - ngpmin == 3)
{
lemin = ngpmin - 1 + lt;
ieee = TRUE;
/* ( Non twos-complement machines, with gradual underflow;
e.g., IEEE standard followers ) */
}
else
{
lemin = MIN (ngpmin, gpmin);
/* ( A guess; no known machine ) */
iwarn = TRUE;
}
}
else if (ngpmin == gpmin && ngnmin == gnmin)
{
if ((i__1 = ngpmin - ngnmin, abs (i__1)) == 1)
{
lemin = MAX (ngpmin, ngnmin);
/* ( Twos-complement machines, no gradual underflow; e.g.,
CYBER 205 ) */
}
else
{
lemin = MIN (ngpmin, ngnmin);
/* ( A guess; no known machine ) */
iwarn = TRUE;
}
}
else if ((i__1 = ngpmin - ngnmin, abs (i__1)) == 1 && gpmin == gnmin)
{
if (gpmin - MIN (ngpmin, ngnmin) == 3)
{
lemin = MAX (ngpmin, ngnmin) - 1 + lt;
/* ( Twos-complement machines with gradual underflow; no
known machine ) */
}
else
{
lemin = MIN (ngpmin, ngnmin);
/* ( A guess; no known machine ) */
iwarn = TRUE;
}
}
else
{
/* Computing MIN */
i__1 = MIN (ngpmin, ngnmin), i__1 = MIN (i__1, gpmin);
lemin = MIN (i__1, gnmin);
/* ( A guess; no known machine ) */
iwarn = TRUE;
}
/* Comment out this if block if EMIN is ok */
if (iwarn)
{
first = TRUE;
Melder_warning ("\n\n WARNING. The value EMIN may be incorrect:- " "EMIN = %8i\n"
"If, after inspection, the value EMIN looks acceptable"
"please comment out \n the IF block as marked within the"
"code of routine DLAMC2, \n otherwise supply EMIN" "explicitly.\n", lemin);
}
/* ** Assume IEEE arithmetic if we found denormalised numbers above,
or if arithmetic seems to round in the IEEE style, determined in
routine DLAMC1. A true IEEE machine should have both things true;
however, faulty machines may have one or the other. */
ieee = ieee || lieee1;
/* Compute RMIN by successive division by BETA. We could compute RMIN
as BASE**( EMIN - 1 ), but some machines underflow during this
computation. */
lrmin = 1.;
i__1 = 1 - lemin;
for (i = 1; i <= 1 - lemin; ++i)
{
d__1 = lrmin * rbase;
lrmin = dlamc3_ (&d__1, &zero);
/* L30: */
}
/* Finally, call DLAMC5 to compute EMAX and RMAX. */
dlamc5_ (&lbeta, <, &lemin, &ieee, &lemax, &lrmax);
}
*beta = lbeta;
*t = lt;
*rnd = lrnd;
*eps = leps;
*emin = lemin;
*rmin = lrmin;
*emax = lemax;
*rmax = lrmax;
return 0;
} /* dlamc2_ */
static double dlamc3_ (double *a, double *b)
/* Purpose =======
dlamc3_ is intended to force A and B to be stored prior to doing the
addition of A and B , for use in situations where optimizers might hold
one of these in a register.
Arguments =========
A, B (input) DOUBLE PRECISION The values A and B.
===================================================================== */
{
double ret_val;
ret_val = *a + *b;
return ret_val;
} /* dlamc3_ */
static int dlamc4_ (long *emin, double *start, long *base)
{
/* -- LAPACK auxiliary routine (version 2.0) -- Univ. of Tennessee, Univ.
of California Berkeley, NAG Ltd., Courant Institute, Argonne National
Lab, and Rice University October 31, 1992
Purpose =======
DLAMC4 is a service routine for DLAMC2.
Arguments =========
EMIN (output) EMIN The minimum exponent before (gradual) underflow,
computed by
setting A = START and dividing by BASE until the previous A can not be
recovered.
START (input) DOUBLE PRECISION The starting point for determining
EMIN.
BASE (input) INTEGER The base of the machine.
===================================================================== */
/* System generated locals */
long i__1;
double d__1;
/* Local variables */
static double zero, a;
static long i;
static double rbase, b1, b2, c1, c2, d1, d2;
static double one;
a = *start;
one = 1.;
rbase = one / *base;
zero = 0.;
*emin = 1;
d__1 = a * rbase;
b1 = dlamc3_ (&d__1, &zero);
c1 = a;
c2 = a;
d1 = a;
d2 = a;
/* + WHILE( ( C1.EQ.A ).AND.( C2.EQ.A ).AND. $ ( D1.EQ.A ).AND.( D2.EQ.A
) )LOOP */
L10:
if (c1 == a && c2 == a && d1 == a && d2 == a)
{
--(*emin);
a = b1;
d__1 = a / *base;
b1 = dlamc3_ (&d__1, &zero);
d__1 = b1 * *base;
c1 = dlamc3_ (&d__1, &zero);
d1 = zero;
i__1 = *base;
for (i = 1; i <= *base; ++i)
{
d1 += b1;
/* L20: */
}
d__1 = a * rbase;
b2 = dlamc3_ (&d__1, &zero);
d__1 = b2 / rbase;
c2 = dlamc3_ (&d__1, &zero);
d2 = zero;
i__1 = *base;
for (i = 1; i <= *base; ++i)
{
d2 += b2;
/* L30: */
}
goto L10;
}
/* + END WHILE */
return 0;
} /* dlamc4_ */
static int dlamc5_ (long *beta, long *p, long *emin, long *ieee, long *emax, double *rmax)
{
/*
First compute LEXP and UEXP, two powers of 2 that bound abs(EMIN). We
then assume that EMAX + abs(EMIN) will sum approximately to the bound
that is closest to abs(EMIN). (EMAX is the exponent of the required
number RMAX). */
/* Table of constant values */
static double c_b5 = 0.;
/* System generated locals */
long i__1;
double d__1;
/* Local variables */
static long lexp;
static double oldy;
static long uexp, i;
static double y, z;
static long nbits;
static double recbas;
static long exbits, expsum, try__;
lexp = 1;
exbits = 1;
L10:
try__ = lexp << 1;
if (try__ <= -(*emin))
{
lexp = try__;
++exbits;
goto L10;
}
if (lexp == -(*emin))
{
uexp = lexp;
}
else
{
uexp = try__;
++exbits;
}
/* Now -LEXP is less than or equal to EMIN, and -UEXP is greater than or
equal to EMIN. EXBITS is the number of bits needed to store the
exponent. */
if (uexp + *emin > -lexp - *emin)
{
expsum = lexp << 1;
}
else
{
expsum = uexp << 1;
}
/* EXPSUM is the exponent range, approximately equal to EMAX - EMIN + 1 .
*/
*emax = expsum + *emin - 1;
nbits = exbits + 1 + *p;
/* NBITS is the total number of bits needed to store a floating-point
number. */
if (nbits % 2 == 1 && *beta == 2)
{
/* Either there are an odd number of bits used to store a
floating-point number, which is unlikely, or some bits are not
used in the representation of numbers, which is possible, (e.g.
Cray machines) or the mantissa has an implicit bit, (e.g. IEEE
machines, Dec Vax machines), which is perhaps the most likely. We
have to assume the last alternative. If this is true, then we need
to reduce EMAX by one because there must be some way of
representing zero in an implicit-bit system. On machines like
Cray, we are reducing EMAX by one unnecessarily. */
--(*emax);
}
if (*ieee)
{
/* Assume we are on an IEEE machine which reserves one exponent for
infinity and NaN. */
--(*emax);
}
/* Now create RMAX, the largest machine number, which should be equal to
(1.0 - BETA**(-P)) * BETA**EMAX . First compute 1.0 - BETA**(-P),
being careful that the result is less than 1.0 . */
recbas = 1. / *beta;
z = *beta - 1.;
y = 0.;
i__1 = *p;
for (i = 1; i <= *p; ++i)
{
z *= recbas;
if (y < 1.)
{
oldy = y;
}
y = dlamc3_ (&y, &z);
/* L20: */
}
if (y >= 1.)
{
y = oldy;
}
/* Now multiply by BETA**EMAX to get RMAX. */
i__1 = *emax;
for (i = 1; i <= *emax; ++i)
{
d__1 = y * *beta;
y = dlamc3_ (&d__1, &c_b5);
/* L30: */
}
*rmax = y;
return 0;
} /* dlamc5_ */
double NUMblas_dnrm2 (long *n, double *x, long *incx)
{
/* The following loop is equivalent to this call to the LAPACK auxiliary
routine: CALL DLASSQ( N, X, INCX, SCALE, SSQ ) */
/* System generated locals */
long i__1, i__2;
double ret_val, d__1;
/* Local variables */
static double norm, scale, absxi;
static long ix;
static double ssq;
--x;
/* Function Body */
if (*n < 1 || *incx < 1)
{
norm = 0.;
}
else if (*n == 1)
{
norm = fabs (x[1]);
}
else
{
scale = 0.;
ssq = 1.;
i__1 = (*n - 1) * *incx + 1;
i__2 = *incx;
for (ix = 1; i__2 < 0 ? ix >= i__1 : ix <= i__1; ix += i__2)
{
if (x[ix] != 0.)
{
absxi = (d__1 = x[ix], fabs (d__1));
if (scale < absxi)
{
/* Computing 2nd power */
d__1 = scale / absxi;
ssq = ssq * (d__1 * d__1) + 1.;
scale = absxi;
}
else
{
/* Computing 2nd power */
d__1 = absxi / scale;
ssq += d__1 * d__1;
}
}
/* L10: */
}
norm = scale * sqrt (ssq);
}
ret_val = norm;
return ret_val;
} /* NUMblas_dnrm2 */
int NUMblas_drot (long *n, double *dx, long *incx, double *dy, long *incy, double *c__, double *s)
{
/* System generated locals */
long i__1;
/* Local variables */
static long i__;
static double dtemp;
static long ix, iy;
/* applies a plane rotation. jack dongarra, linpack, 3/11/78. modified
12/3/93, array(1) declarations changed to array(*) Parameter
adjustments */
--dy;
--dx;
/* Function Body */
if (*n <= 0)
{
return 0;
}
if (*incx == 1 && *incy == 1)
{
goto L20;
}
/* code for unequal increments or equal increments not equal to 1 */
ix = 1;
iy = 1;
if (*incx < 0)
{
ix = (-(*n) + 1) * *incx + 1;
}
if (*incy < 0)
{
iy = (-(*n) + 1) * *incy + 1;
}
i__1 = *n;
for (i__ = 1; i__ <= i__1; ++i__)
{
dtemp = *c__ * dx[ix] + *s * dy[iy];
dy[iy] = *c__ * dy[iy] - *s * dx[ix];
dx[ix] = dtemp;
ix += *incx;
iy += *incy;
/* L10: */
}
return 0;
/* code for both increments equal to 1 */
L20:
i__1 = *n;
for (i__ = 1; i__ <= i__1; ++i__)
{
dtemp = *c__ * dx[i__] + *s * dy[i__];
dy[i__] = *c__ * dy[i__] - *s * dx[i__];
dx[i__] = dtemp;
/* L30: */
}
return 0;
} /* NUMblas_drot */
int NUMblas_dscal (long *n, double *da, double *dx, long *incx)
{
/* System generated locals */
long i__1, i__2;
/* Local variables */
static long i__, m, nincx, mp1;
/* Parameter adjustments */
--dx;
/* Function Body */
if (*n <= 0 || *incx <= 0)
{
return 0;
}
if (*incx == 1)
{
goto L20;
}
/* code for increment not equal to 1 */
nincx = *n * *incx;
i__1 = nincx;
i__2 = *incx;
for (i__ = 1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ += i__2)
{
dx[i__] = *da * dx[i__];
/* L10: */
}
return 0;
/* code for increment equal to 1 clean-up loop */
L20:
m = *n % 5;
if (m == 0)
{
goto L40;
}
i__2 = m;
for (i__ = 1; i__ <= i__2; ++i__)
{
dx[i__] = *da * dx[i__];
/* L30: */
}
if (*n < 5)
{
return 0;
}
L40:
mp1 = m + 1;
i__2 = *n;
for (i__ = mp1; i__ <= i__2; i__ += 5)
{
dx[i__] = *da * dx[i__];
dx[i__ + 1] = *da * dx[i__ + 1];
dx[i__ + 2] = *da * dx[i__ + 2];
dx[i__ + 3] = *da * dx[i__ + 3];
dx[i__ + 4] = *da * dx[i__ + 4];
/* L50: */
}
return 0;
} /* dscal_ */
int NUMblas_dswap (long *n, double *dx, long *incx, double *dy, long *incy)
{
/* System generated locals */
long i__1;
/* Local variables */
static long i__, m;
static double dtemp;
static long ix, iy, mp1;
/* interchanges two vectors. uses unrolled loops for increments equal
one. jack dongarra, linpack, 3/11/78. modified 12/3/93, array(1)
declarations changed to array(*) Parameter adjustments */
--dy;
--dx;
/* Function Body */
if (*n <= 0)
{
return 0;
}
if (*incx == 1 && *incy == 1)
{
goto L20;
}
/* code for unequal increments or equal increments not equal to 1 */
ix = 1;
iy = 1;
if (*incx < 0)
{
ix = (-(*n) + 1) * *incx + 1;
}
if (*incy < 0)
{
iy = (-(*n) + 1) * *incy + 1;
}
i__1 = *n;
for (i__ = 1; i__ <= i__1; ++i__)
{
dtemp = dx[ix];
dx[ix] = dy[iy];
dy[iy] = dtemp;
ix += *incx;
iy += *incy;
/* L10: */
}
return 0;
/* code for both increments equal to 1 clean-up loop */
L20:
m = *n % 3;
if (m == 0)
{
goto L40;
}
i__1 = m;
for (i__ = 1; i__ <= i__1; ++i__)
{
dtemp = dx[i__];
dx[i__] = dy[i__];
dy[i__] = dtemp;
/* L30: */
}
if (*n < 3)
{
return 0;
}
L40:
mp1 = m + 1;
i__1 = *n;
for (i__ = mp1; i__ <= i__1; i__ += 3)
{
dtemp = dx[i__];
dx[i__] = dy[i__];
dy[i__] = dtemp;
dtemp = dx[i__ + 1];
dx[i__ + 1] = dy[i__ + 1];
dy[i__ + 1] = dtemp;
dtemp = dx[i__ + 2];
dx[i__ + 2] = dy[i__ + 2];
dy[i__ + 2] = dtemp;
/* L50: */
}
return 0;
} /* NUMblas_dswap */
int NUMblas_dsymv (char *uplo, long *n, double *alpha, double *a, long *lda, double *x, long *incx, double *beta,
double *y, long *incy)
{
/* System generated locals */
long a_dim1, a_offset, i__1, i__2;
/* Local variables */
static long info;
static double temp1, temp2;
static long i__, j;
static long ix, iy, jx, jy, kx, ky;
#define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
a_dim1 = *lda;
a_offset = 1 + a_dim1 * 1;
a -= a_offset;
--x;
--y;
/* Function Body */
info = 0;
if (!lsame_ (uplo, "U") && !lsame_ (uplo, "L"))
{
info = 1;
}
else if (*n < 0)
{
info = 2;
}
else if (*lda < MAX (1, *n))
{
info = 5;
}
else if (*incx == 0)
{
info = 7;
}
else if (*incy == 0)
{
info = 10;
}
if (info != 0)
{
xerbla_ ("DSYMV ", &info);
return 0;
}
/* Quick return if possible. */
if (*n == 0 || *alpha == 0. && *beta == 1.)
{
return 0;
}
/* Set up the start points in X and Y. */
if (*incx > 0)
{
kx = 1;
}
else
{
kx = 1 - (*n - 1) * *incx;
}
if (*incy > 0)
{
ky = 1;
}
else
{
ky = 1 - (*n - 1) * *incy;
}
/* Start the operations. In this version the elements of A are accessed
sequentially with one pass through the triangular part of A. First
form y := beta*y. */
if (*beta != 1.)
{
if (*incy == 1)
{
if (*beta == 0.)
{
i__1 = *n;
for (i__ = 1; i__ <= i__1; ++i__)
{
y[i__] = 0.;
/* L10: */
}
}
else
{
i__1 = *n;
for (i__ = 1; i__ <= i__1; ++i__)
{
y[i__] = *beta * y[i__];
/* L20: */
}
}
}
else
{
iy = ky;
if (*beta == 0.)
{
i__1 = *n;
for (i__ = 1; i__ <= i__1; ++i__)
{
y[iy] = 0.;
iy += *incy;
/* L30: */
}
}
else
{
i__1 = *n;
for (i__ = 1; i__ <= i__1; ++i__)
{
y[iy] = *beta * y[iy];
iy += *incy;
/* L40: */
}
}
}
}
if (*alpha == 0.)
{
return 0;
}
if (lsame_ (uplo, "U"))
{
/* Form y when A is stored in upper triangle. */
if (*incx == 1 && *incy == 1)
{
i__1 = *n;
for (j = 1; j <= i__1; ++j)
{
temp1 = *alpha * x[j];
temp2 = 0.;
i__2 = j - 1;
for (i__ = 1; i__ <= i__2; ++i__)
{
y[i__] += temp1 * a_ref (i__, j);
temp2 += a_ref (i__, j) * x[i__];
/* L50: */
}
y[j] = y[j] + temp1 * a_ref (j, j) + *alpha * temp2;
/* L60: */
}
}
else
{
jx = kx;
jy = ky;
i__1 = *n;
for (j = 1; j <= i__1; ++j)
{
temp1 = *alpha * x[jx];
temp2 = 0.;
ix = kx;
iy = ky;
i__2 = j - 1;
for (i__ = 1; i__ <= i__2; ++i__)
{
y[iy] += temp1 * a_ref (i__, j);
temp2 += a_ref (i__, j) * x[ix];
ix += *incx;
iy += *incy;
/* L70: */
}
y[jy] = y[jy] + temp1 * a_ref (j, j) + *alpha * temp2;
jx += *incx;
jy += *incy;
/* L80: */
}
}
}
else
{
/* Form y when A is stored in lower triangle. */
if (*incx == 1 && *incy == 1)
{
i__1 = *n;
for (j = 1; j <= i__1; ++j)
{
temp1 = *alpha * x[j];
temp2 = 0.;
y[j] += temp1 * a_ref (j, j);
i__2 = *n;
for (i__ = j + 1; i__ <= i__2; ++i__)
{
y[i__] += temp1 * a_ref (i__, j);
temp2 += a_ref (i__, j) * x[i__];
/* L90: */
}
y[j] += *alpha * temp2;
/* L100: */
}
}
else
{
jx = kx;
jy = ky;
i__1 = *n;
for (j = 1; j <= i__1; ++j)
{
temp1 = *alpha * x[jx];
temp2 = 0.;
y[jy] += temp1 * a_ref (j, j);
ix = jx;
iy = jy;
i__2 = *n;
for (i__ = j + 1; i__ <= i__2; ++i__)
{
ix += *incx;
iy += *incy;
y[iy] += temp1 * a_ref (i__, j);
temp2 += a_ref (i__, j) * x[ix];
/* L110: */
}
y[jy] += *alpha * temp2;
jx += *incx;
jy += *incy;
/* L120: */
}
}
}
return 0;
} /* NUMblas_dsymv */
#undef a_ref
int NUMblas_dsyr2 (char *uplo, long *n, double *alpha, double *x, long *incx, double *y, long *incy, double *a,
long *lda)
{
/* System generated locals */
long a_dim1, a_offset, i__1, i__2;
/* Local variables */
static long info;
static double temp1, temp2;
static long i__, j;
static long ix, iy, jx, jy, kx, ky;
#define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
--x;
--y;
a_dim1 = *lda;
a_offset = 1 + a_dim1 * 1;
a -= a_offset;
/* Function Body */
info = 0;
if (!lsame_ (uplo, "U") && !lsame_ (uplo, "L"))
{
info = 1;
}
else if (*n < 0)
{
info = 2;
}
else if (*incx == 0)
{
info = 5;
}
else if (*incy == 0)
{
info = 7;
}
else if (*lda < MAX (1, *n))
{
info = 9;
}
if (info != 0)
{
xerbla_ ("DSYR2 ", &info);
return 0;
}
/* Quick return if possible. */
if (*n == 0 || *alpha == 0.)
{
return 0;
}
/* Set up the start points in X and Y if the increments are not both
unity. */
if (*incx != 1 || *incy != 1)
{
if (*incx > 0)
{
kx = 1;
}
else
{
kx = 1 - (*n - 1) * *incx;
}
if (*incy > 0)
{
ky = 1;
}
else
{
ky = 1 - (*n - 1) * *incy;
}
jx = kx;
jy = ky;
}
/* Start the operations. In this version the elements of A are accessed
sequentially with one pass through the triangular part of A. */
if (lsame_ (uplo, "U"))
{
/* Form A when A is stored in the upper triangle. */
if (*incx == 1 && *incy == 1)
{
i__1 = *n;
for (j = 1; j <= i__1; ++j)
{
if (x[j] != 0. || y[j] != 0.)
{
temp1 = *alpha * y[j];
temp2 = *alpha * x[j];
i__2 = j;
for (i__ = 1; i__ <= i__2; ++i__)
{
a_ref (i__, j) = a_ref (i__, j) + x[i__] * temp1 + y[i__] * temp2;
/* L10: */
}
}
/* L20: */
}
}
else
{
i__1 = *n;
for (j = 1; j <= i__1; ++j)
{
if (x[jx] != 0. || y[jy] != 0.)
{
temp1 = *alpha * y[jy];
temp2 = *alpha * x[jx];
ix = kx;
iy = ky;
i__2 = j;
for (i__ = 1; i__ <= i__2; ++i__)
{
a_ref (i__, j) = a_ref (i__, j) + x[ix] * temp1 + y[iy] * temp2;
ix += *incx;
iy += *incy;
/* L30: */
}
}
jx += *incx;
jy += *incy;
/* L40: */
}
}
}
else
{
/* Form A when A is stored in the lower triangle. */
if (*incx == 1 && *incy == 1)
{
i__1 = *n;
for (j = 1; j <= i__1; ++j)
{
if (x[j] != 0. || y[j] != 0.)
{
temp1 = *alpha * y[j];
temp2 = *alpha * x[j];
i__2 = *n;
for (i__ = j; i__ <= i__2; ++i__)
{
a_ref (i__, j) = a_ref (i__, j) + x[i__] * temp1 + y[i__] * temp2;
/* L50: */
}
}
/* L60: */
}
}
else
{
i__1 = *n;
for (j = 1; j <= i__1; ++j)
{
if (x[jx] != 0. || y[jy] != 0.)
{
temp1 = *alpha * y[jy];
temp2 = *alpha * x[jx];
ix = jx;
iy = jy;
i__2 = *n;
for (i__ = j; i__ <= i__2; ++i__)
{
a_ref (i__, j) = a_ref (i__, j) + x[ix] * temp1 + y[iy] * temp2;
ix += *incx;
iy += *incy;
/* L70: */
}
}
jx += *incx;
jy += *incy;
/* L80: */
}
}
}
return 0;
} /* NUMblas_dsyr2 */
#undef a_ref
int NUMblas_dsyr2k (char *uplo, char *trans, long *n, long *k, double *alpha, double *a, long *lda, double *b,
long *ldb, double *beta, double *c__, long *ldc)
{
/* System generated locals */
long a_dim1, a_offset, b_dim1, b_offset, c_dim1, c_offset, i__1, i__2, i__3;
/* Local variables */
static long info;
static double temp1, temp2;
static long i__, j, l;
static long nrowa;
static long upper;
#define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
#define b_ref(a_1,a_2) b[(a_2)*b_dim1 + a_1]
#define c___ref(a_1,a_2) c__[(a_2)*c_dim1 + a_1]
a_dim1 = *lda;
a_offset = 1 + a_dim1 * 1;
a -= a_offset;
b_dim1 = *ldb;
b_offset = 1 + b_dim1 * 1;
b -= b_offset;
c_dim1 = *ldc;
c_offset = 1 + c_dim1 * 1;
c__ -= c_offset;
/* Function Body */
if (lsame_ (trans, "N"))
{
nrowa = *n;
}
else
{
nrowa = *k;
}
upper = lsame_ (uplo, "U");
info = 0;
if (!upper && !lsame_ (uplo, "L"))
{
info = 1;
}
else if (!lsame_ (trans, "N") && !lsame_ (trans, "T") && !lsame_ (trans, "C"))
{
info = 2;
}
else if (*n < 0)
{
info = 3;
}
else if (*k < 0)
{
info = 4;
}
else if (*lda < MAX (1, nrowa))
{
info = 7;
}
else if (*ldb < MAX (1, nrowa))
{
info = 9;
}
else if (*ldc < MAX (1, *n))
{
info = 12;
}
if (info != 0)
{
xerbla_ ("DSYR2K", &info);
return 0;
}
/* Quick return if possible. */
if (*n == 0 || (*alpha == 0. || *k == 0) && *beta == 1.)
{
return 0;
}
/* And when alpha.eq.zero. */
if (*alpha == 0.)
{
if (upper)
{
if (*beta == 0.)
{
i__1 = *n;
for (j = 1; j <= i__1; ++j)
{
i__2 = j;
for (i__ = 1; i__ <= i__2; ++i__)
{
c___ref (i__, j) = 0.;
/* L10: */
}
/* L20: */
}
}
else
{
i__1 = *n;
for (j = 1; j <= i__1; ++j)
{
i__2 = j;
for (i__ = 1; i__ <= i__2; ++i__)
{
c___ref (i__, j) = *beta * c___ref (i__, j);
/* L30: */
}
/* L40: */
}
}
}
else
{
if (*beta == 0.)
{
i__1 = *n;
for (j = 1; j <= i__1; ++j)
{
i__2 = *n;
for (i__ = j; i__ <= i__2; ++i__)
{
c___ref (i__, j) = 0.;
/* L50: */
}
/* L60: */
}
}
else
{
i__1 = *n;
for (j = 1; j <= i__1; ++j)
{
i__2 = *n;
for (i__ = j; i__ <= i__2; ++i__)
{
c___ref (i__, j) = *beta * c___ref (i__, j);
/* L70: */
}
/* L80: */
}
}
}
return 0;
}
/* Start the operations. */
if (lsame_ (trans, "N"))
{
/* Form C := alpha*A*B' + alpha*B*A' + C. */
if (upper)
{
i__1 = *n;
for (j = 1; j <= i__1; ++j)
{
if (*beta == 0.)
{
i__2 = j;
for (i__ = 1; i__ <= i__2; ++i__)
{
c___ref (i__, j) = 0.;
/* L90: */
}
}
else if (*beta != 1.)
{
i__2 = j;
for (i__ = 1; i__ <= i__2; ++i__)
{
c___ref (i__, j) = *beta * c___ref (i__, j);
/* L100: */
}
}
i__2 = *k;
for (l = 1; l <= i__2; ++l)
{
if (a_ref (j, l) != 0. || b_ref (j, l) != 0.)
{
temp1 = *alpha * b_ref (j, l);
temp2 = *alpha * a_ref (j, l);
i__3 = j;
for (i__ = 1; i__ <= i__3; ++i__)
{
c___ref (i__, j) =
c___ref (i__, j) + a_ref (i__, l) * temp1 + b_ref (i__, l) * temp2;
/* L110: */
}
}
/* L120: */
}
/* L130: */
}
}
else
{
i__1 = *n;
for (j = 1; j <= i__1; ++j)
{
if (*beta == 0.)
{
i__2 = *n;
for (i__ = j; i__ <= i__2; ++i__)
{
c___ref (i__, j) = 0.;
/* L140: */
}
}
else if (*beta != 1.)
{
i__2 = *n;
for (i__ = j; i__ <= i__2; ++i__)
{
c___ref (i__, j) = *beta * c___ref (i__, j);
/* L150: */
}
}
i__2 = *k;
for (l = 1; l <= i__2; ++l)
{
if (a_ref (j, l) != 0. || b_ref (j, l) != 0.)
{
temp1 = *alpha * b_ref (j, l);
temp2 = *alpha * a_ref (j, l);
i__3 = *n;
for (i__ = j; i__ <= i__3; ++i__)
{
c___ref (i__, j) =
c___ref (i__, j) + a_ref (i__, l) * temp1 + b_ref (i__, l) * temp2;
/* L160: */
}
}
/* L170: */
}
/* L180: */
}
}
}
else
{
/* Form C := alpha*A'*B + alpha*B'*A + C. */
if (upper)
{
i__1 = *n;
for (j = 1; j <= i__1; ++j)
{
i__2 = j;
for (i__ = 1; i__ <= i__2; ++i__)
{
temp1 = 0.;
temp2 = 0.;
i__3 = *k;
for (l = 1; l <= i__3; ++l)
{
temp1 += a_ref (l, i__) * b_ref (l, j);
temp2 += b_ref (l, i__) * a_ref (l, j);
/* L190: */
}
if (*beta == 0.)
{
c___ref (i__, j) = *alpha * temp1 + *alpha * temp2;
}
else
{
c___ref (i__, j) = *beta * c___ref (i__, j) + *alpha * temp1 + *alpha * temp2;
}
/* L200: */
}
/* L210: */
}
}
else
{
i__1 = *n;
for (j = 1; j <= i__1; ++j)
{
i__2 = *n;
for (i__ = j; i__ <= i__2; ++i__)
{
temp1 = 0.;
temp2 = 0.;
i__3 = *k;
for (l = 1; l <= i__3; ++l)
{
temp1 += a_ref (l, i__) * b_ref (l, j);
temp2 += b_ref (l, i__) * a_ref (l, j);
/* L220: */
}
if (*beta == 0.)
{
c___ref (i__, j) = *alpha * temp1 + *alpha * temp2;
}
else
{
c___ref (i__, j) = *beta * c___ref (i__, j) + *alpha * temp1 + *alpha * temp2;
}
/* L230: */
}
/* L240: */
}
}
}
return 0;
} /* NUMblas_dsyr2k */
#undef c___ref
#undef b_ref
#undef a_ref
int NUMblas_dtrmm (char *side, char *uplo, char *transa, char *diag, long *m, long *n, double *alpha, double *a,
long *lda, double *b, long *ldb)
{
/* System generated locals */
long a_dim1, a_offset, b_dim1, b_offset, i__1, i__2, i__3;
/* Local variables */
static long info;
static double temp;
static long i__, j, k;
static long lside;
static long nrowa;
static long upper;
static long nounit;
#define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
#define b_ref(a_1,a_2) b[(a_2)*b_dim1 + a_1]
/* Level 3 Blas routine. -- Written on 8-February-1989. Jack Dongarra,
Argonne National Laboratory. Iain Duff, AERE Harwell. Jeremy Du Croz,
Numerical Algorithms Group Ltd. Sven Hammarling, Numerical Algorithms
Group Ltd. Test the input parameters. Parameter adjustments */
a_dim1 = *lda;
a_offset = 1 + a_dim1 * 1;
a -= a_offset;
b_dim1 = *ldb;
b_offset = 1 + b_dim1 * 1;
b -= b_offset;
/* Function Body */
lside = lsame_ (side, "L");
if (lside)
{
nrowa = *m;
}
else
{
nrowa = *n;
}
nounit = lsame_ (diag, "N");
upper = lsame_ (uplo, "U");
info = 0;
if (!lside && !lsame_ (side, "R"))
{
info = 1;
}
else if (!upper && !lsame_ (uplo, "L"))
{
info = 2;
}
else if (!lsame_ (transa, "N") && !lsame_ (transa, "T") && !lsame_ (transa, "C"))
{
info = 3;
}
else if (!lsame_ (diag, "U") && !lsame_ (diag, "N"))
{
info = 4;
}
else if (*m < 0)
{
info = 5;
}
else if (*n < 0)
{
info = 6;
}
else if (*lda < MAX (1, nrowa))
{
info = 9;
}
else if (*ldb < MAX (1, *m))
{
info = 11;
}
if (info != 0)
{
xerbla_ ("DTRMM ", &info);
return 0;
}
/* Quick return if possible. */
if (*n == 0)
{
return 0;
}
/* And when alpha.eq.zero. */
if (*alpha == 0.)
{
i__1 = *n;
for (j = 1; j <= i__1; ++j)
{
i__2 = *m;
for (i__ = 1; i__ <= i__2; ++i__)
{
b_ref (i__, j) = 0.;
/* L10: */
}
/* L20: */
}
return 0;
}
/* Start the operations. */
if (lside)
{
if (lsame_ (transa, "N"))
{
/* Form B := alpha*A*B. */
if (upper)
{
i__1 = *n;
for (j = 1; j <= i__1; ++j)
{
i__2 = *m;
for (k = 1; k <= i__2; ++k)
{
if (b_ref (k, j) != 0.)
{
temp = *alpha * b_ref (k, j);
i__3 = k - 1;
for (i__ = 1; i__ <= i__3; ++i__)
{
b_ref (i__, j) = b_ref (i__, j) + temp * a_ref (i__, k);
/* L30: */
}
if (nounit)
{
temp *= a_ref (k, k);
}
b_ref (k, j) = temp;
}
/* L40: */
}
/* L50: */
}
}
else
{
i__1 = *n;
for (j = 1; j <= i__1; ++j)
{
for (k = *m; k >= 1; --k)
{
if (b_ref (k, j) != 0.)
{
temp = *alpha * b_ref (k, j);
b_ref (k, j) = temp;
if (nounit)
{
b_ref (k, j) = b_ref (k, j) * a_ref (k, k);
}
i__2 = *m;
for (i__ = k + 1; i__ <= i__2; ++i__)
{
b_ref (i__, j) = b_ref (i__, j) + temp * a_ref (i__, k);
/* L60: */
}
}
/* L70: */
}
/* L80: */
}
}
}
else
{
/* Form B := alpha*A'*B. */
if (upper)
{
i__1 = *n;
for (j = 1; j <= i__1; ++j)
{
for (i__ = *m; i__ >= 1; --i__)
{
temp = b_ref (i__, j);
if (nounit)
{
temp *= a_ref (i__, i__);
}
i__2 = i__ - 1;
for (k = 1; k <= i__2; ++k)
{
temp += a_ref (k, i__) * b_ref (k, j);
/* L90: */
}
b_ref (i__, j) = *alpha * temp;
/* L100: */
}
/* L110: */
}
}
else
{
i__1 = *n;
for (j = 1; j <= i__1; ++j)
{
i__2 = *m;
for (i__ = 1; i__ <= i__2; ++i__)
{
temp = b_ref (i__, j);
if (nounit)
{
temp *= a_ref (i__, i__);
}
i__3 = *m;
for (k = i__ + 1; k <= i__3; ++k)
{
temp += a_ref (k, i__) * b_ref (k, j);
/* L120: */
}
b_ref (i__, j) = *alpha * temp;
/* L130: */
}
/* L140: */
}
}
}
}
else
{
if (lsame_ (transa, "N"))
{
/* Form B := alpha*B*A. */
if (upper)
{
for (j = *n; j >= 1; --j)
{
temp = *alpha;
if (nounit)
{
temp *= a_ref (j, j);
}
i__1 = *m;
for (i__ = 1; i__ <= i__1; ++i__)
{
b_ref (i__, j) = temp * b_ref (i__, j);
/* L150: */
}
i__1 = j - 1;
for (k = 1; k <= i__1; ++k)
{
if (a_ref (k, j) != 0.)
{
temp = *alpha * a_ref (k, j);
i__2 = *m;
for (i__ = 1; i__ <= i__2; ++i__)
{
b_ref (i__, j) = b_ref (i__, j) + temp * b_ref (i__, k);
/* L160: */
}
}
/* L170: */
}
/* L180: */
}
}
else
{
i__1 = *n;
for (j = 1; j <= i__1; ++j)
{
temp = *alpha;
if (nounit)
{
temp *= a_ref (j, j);
}
i__2 = *m;
for (i__ = 1; i__ <= i__2; ++i__)
{
b_ref (i__, j) = temp * b_ref (i__, j);
/* L190: */
}
i__2 = *n;
for (k = j + 1; k <= i__2; ++k)
{
if (a_ref (k, j) != 0.)
{
temp = *alpha * a_ref (k, j);
i__3 = *m;
for (i__ = 1; i__ <= i__3; ++i__)
{
b_ref (i__, j) = b_ref (i__, j) + temp * b_ref (i__, k);
/* L200: */
}
}
/* L210: */
}
/* L220: */
}
}
}
else
{
/* Form B := alpha*B*A'. */
if (upper)
{
i__1 = *n;
for (k = 1; k <= i__1; ++k)
{
i__2 = k - 1;
for (j = 1; j <= i__2; ++j)
{
if (a_ref (j, k) != 0.)
{
temp = *alpha * a_ref (j, k);
i__3 = *m;
for (i__ = 1; i__ <= i__3; ++i__)
{
b_ref (i__, j) = b_ref (i__, j) + temp * b_ref (i__, k);
/* L230: */
}
}
/* L240: */
}
temp = *alpha;
if (nounit)
{
temp *= a_ref (k, k);
}
if (temp != 1.)
{
i__2 = *m;
for (i__ = 1; i__ <= i__2; ++i__)
{
b_ref (i__, k) = temp * b_ref (i__, k);
/* L250: */
}
}
/* L260: */
}
}
else
{
for (k = *n; k >= 1; --k)
{
i__1 = *n;
for (j = k + 1; j <= i__1; ++j)
{
if (a_ref (j, k) != 0.)
{
temp = *alpha * a_ref (j, k);
i__2 = *m;
for (i__ = 1; i__ <= i__2; ++i__)
{
b_ref (i__, j) = b_ref (i__, j) + temp * b_ref (i__, k);
/* L270: */
}
}
/* L280: */
}
temp = *alpha;
if (nounit)
{
temp *= a_ref (k, k);
}
if (temp != 1.)
{
i__1 = *m;
for (i__ = 1; i__ <= i__1; ++i__)
{
b_ref (i__, k) = temp * b_ref (i__, k);
/* L290: */
}
}
/* L300: */
}
}
}
}
return 0;
} /* NUMblas_dtrmm */
#undef b_ref
#undef a_ref
int NUMblas_dtrmv (char *uplo, char *trans, char *diag, long *n, double *a, long *lda, double *x, long *incx)
{
/* System generated locals */
long a_dim1, a_offset, i__1, i__2;
/* Local variables */
static long info;
static double temp;
static long i__, j;
static long ix, jx, kx;
static long nounit;
#define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
/* -- Written on 22-October-1986. Jack Dongarra, Argonne National Lab.
Jeremy Du Croz, Nag Central Office. Sven Hammarling, Nag Central
Office. Richard Hanson, Sandia National Labs. Test the input
parameters. Parameter adjustments */
a_dim1 = *lda;
a_offset = 1 + a_dim1 * 1;
a -= a_offset;
--x;
/* Function Body */
info = 0;
if (!lsame_ (uplo, "U") && !lsame_ (uplo, "L"))
{
info = 1;
}
else if (!lsame_ (trans, "N") && !lsame_ (trans, "T") && !lsame_ (trans, "C"))
{
info = 2;
}
else if (!lsame_ (diag, "U") && !lsame_ (diag, "N"))
{
info = 3;
}
else if (*n < 0)
{
info = 4;
}
else if (*lda < MAX (1, *n))
{
info = 6;
}
else if (*incx == 0)
{
info = 8;
}
if (info != 0)
{
xerbla_ ("DTRMV ", &info);
return 0;
}
/* Quick return if possible. */
if (*n == 0)
{
return 0;
}
nounit = lsame_ (diag, "N");
/* Set up the start point in X if the increment is not unity. This will
be ( N - 1 )*INCX too small for descending loops. */
if (*incx <= 0)
{
kx = 1 - (*n - 1) * *incx;
}
else if (*incx != 1)
{
kx = 1;
}
/* Start the operations. In this version the elements of A are accessed
sequentially with one pass through A. */
if (lsame_ (trans, "N"))
{
/* Form x := A*x. */
if (lsame_ (uplo, "U"))
{
if (*incx == 1)
{
i__1 = *n;
for (j = 1; j <= i__1; ++j)
{
if (x[j] != 0.)
{
temp = x[j];
i__2 = j - 1;
for (i__ = 1; i__ <= i__2; ++i__)
{
x[i__] += temp * a_ref (i__, j);
/* L10: */
}
if (nounit)
{
x[j] *= a_ref (j, j);
}
}
/* L20: */
}
}
else
{
jx = kx;
i__1 = *n;
for (j = 1; j <= i__1; ++j)
{
if (x[jx] != 0.)
{
temp = x[jx];
ix = kx;
i__2 = j - 1;
for (i__ = 1; i__ <= i__2; ++i__)
{
x[ix] += temp * a_ref (i__, j);
ix += *incx;
/* L30: */
}
if (nounit)
{
x[jx] *= a_ref (j, j);
}
}
jx += *incx;
/* L40: */
}
}
}
else
{
if (*incx == 1)
{
for (j = *n; j >= 1; --j)
{
if (x[j] != 0.)
{
temp = x[j];
i__1 = j + 1;
for (i__ = *n; i__ >= i__1; --i__)
{
x[i__] += temp * a_ref (i__, j);
/* L50: */
}
if (nounit)
{
x[j] *= a_ref (j, j);
}
}
/* L60: */
}
}
else
{
kx += (*n - 1) * *incx;
jx = kx;
for (j = *n; j >= 1; --j)
{
if (x[jx] != 0.)
{
temp = x[jx];
ix = kx;
i__1 = j + 1;
for (i__ = *n; i__ >= i__1; --i__)
{
x[ix] += temp * a_ref (i__, j);
ix -= *incx;
/* L70: */
}
if (nounit)
{
x[jx] *= a_ref (j, j);
}
}
jx -= *incx;
/* L80: */
}
}
}
}
else
{
/* Form x := A'*x. */
if (lsame_ (uplo, "U"))
{
if (*incx == 1)
{
for (j = *n; j >= 1; --j)
{
temp = x[j];
if (nounit)
{
temp *= a_ref (j, j);
}
for (i__ = j - 1; i__ >= 1; --i__)
{
temp += a_ref (i__, j) * x[i__];
/* L90: */
}
x[j] = temp;
/* L100: */
}
}
else
{
jx = kx + (*n - 1) * *incx;
for (j = *n; j >= 1; --j)
{
temp = x[jx];
ix = jx;
if (nounit)
{
temp *= a_ref (j, j);
}
for (i__ = j - 1; i__ >= 1; --i__)
{
ix -= *incx;
temp += a_ref (i__, j) * x[ix];
/* L110: */
}
x[jx] = temp;
jx -= *incx;
/* L120: */
}
}
}
else
{
if (*incx == 1)
{
i__1 = *n;
for (j = 1; j <= i__1; ++j)
{
temp = x[j];
if (nounit)
{
temp *= a_ref (j, j);
}
i__2 = *n;
for (i__ = j + 1; i__ <= i__2; ++i__)
{
temp += a_ref (i__, j) * x[i__];
/* L130: */
}
x[j] = temp;
/* L140: */
}
}
else
{
jx = kx;
i__1 = *n;
for (j = 1; j <= i__1; ++j)
{
temp = x[jx];
ix = jx;
if (nounit)
{
temp *= a_ref (j, j);
}
i__2 = *n;
for (i__ = j + 1; i__ <= i__2; ++i__)
{
ix += *incx;
temp += a_ref (i__, j) * x[ix];
/* L150: */
}
x[jx] = temp;
jx += *incx;
/* L160: */
}
}
}
}
return 0;
} /* NUMblas_dtrmv */
#undef a_ref
int NUMblas_dtrsm (char *side, char *uplo, char *transa, char *diag, long *m, long *n,
double *alpha, double *a, long *lda, double *b, long *ldb)
{
/* System generated locals */
long a_dim1, a_offset, b_dim1, b_offset, i__1, i__2, i__3;
/* Local variables */
static long info;
static double temp;
static long i__, j, k;
static long lside;
static long nrowa;
static long upper;
static long nounit;
#define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
#define b_ref(a_1,a_2) b[(a_2)*b_dim1 + a_1]
a_dim1 = *lda;
a_offset = 1 + a_dim1 * 1;
a -= a_offset;
b_dim1 = *ldb;
b_offset = 1 + b_dim1 * 1;
b -= b_offset;
/* Function Body */
lside = lsame_ (side, "L");
if (lside)
{
nrowa = *m;
}
else
{
nrowa = *n;
}
nounit = lsame_ (diag, "N");
upper = lsame_ (uplo, "U");
info = 0;
if (!lside && !lsame_ (side, "R"))
{
info = 1;
}
else if (!upper && !lsame_ (uplo, "L"))
{
info = 2;
}
else if (!lsame_ (transa, "N") && !lsame_ (transa, "T") && !lsame_ (transa, "C"))
{
info = 3;
}
else if (!lsame_ (diag, "U") && !lsame_ (diag, "N"))
{
info = 4;
}
else if (*m < 0)
{
info = 5;
}
else if (*n < 0)
{
info = 6;
}
else if (*lda < MAX (1, nrowa))
{
info = 9;
}
else if (*ldb < MAX (1, *m))
{
info = 11;
}
if (info != 0)
{
xerbla_ ("DTRSM ", &info);
return 0;
}
/* Quick return if possible. */
if (*n == 0)
{
return 0;
}
/* And when alpha.eq.zero. */
if (*alpha == 0.)
{
i__1 = *n;
for (j = 1; j <= i__1; ++j)
{
i__2 = *m;
for (i__ = 1; i__ <= i__2; ++i__)
{
b_ref (i__, j) = 0.;
/* L10: */
}
/* L20: */
}
return 0;
}
/* Start the operations. */
if (lside)
{
if (lsame_ (transa, "N"))
{
/* Form B := alpha*inv( A )*B. */
if (upper)
{
i__1 = *n;
for (j = 1; j <= i__1; ++j)
{
if (*alpha != 1.)
{
i__2 = *m;
for (i__ = 1; i__ <= i__2; ++i__)
{
b_ref (i__, j) = *alpha * b_ref (i__, j);
/* L30: */
}
}
for (k = *m; k >= 1; --k)
{
if (b_ref (k, j) != 0.)
{
if (nounit)
{
b_ref (k, j) = b_ref (k, j) / a_ref (k, k);
}
i__2 = k - 1;
for (i__ = 1; i__ <= i__2; ++i__)
{
b_ref (i__, j) = b_ref (i__, j) - b_ref (k, j) * a_ref (i__, k);
/* L40: */
}
}
/* L50: */
}
/* L60: */
}
}
else
{
i__1 = *n;
for (j = 1; j <= i__1; ++j)
{
if (*alpha != 1.)
{
i__2 = *m;
for (i__ = 1; i__ <= i__2; ++i__)
{
b_ref (i__, j) = *alpha * b_ref (i__, j);
/* L70: */
}
}
i__2 = *m;
for (k = 1; k <= i__2; ++k)
{
if (b_ref (k, j) != 0.)
{
if (nounit)
{
b_ref (k, j) = b_ref (k, j) / a_ref (k, k);
}
i__3 = *m;
for (i__ = k + 1; i__ <= i__3; ++i__)
{
b_ref (i__, j) = b_ref (i__, j) - b_ref (k, j) * a_ref (i__, k);
/* L80: */
}
}
/* L90: */
}
/* L100: */
}
}
}
else
{
/* Form B := alpha*inv( A' )*B. */
if (upper)
{
i__1 = *n;
for (j = 1; j <= i__1; ++j)
{
i__2 = *m;
for (i__ = 1; i__ <= i__2; ++i__)
{
temp = *alpha * b_ref (i__, j);
i__3 = i__ - 1;
for (k = 1; k <= i__3; ++k)
{
temp -= a_ref (k, i__) * b_ref (k, j);
/* L110: */
}
if (nounit)
{
temp /= a_ref (i__, i__);
}
b_ref (i__, j) = temp;
/* L120: */
}
/* L130: */
}
}
else
{
i__1 = *n;
for (j = 1; j <= i__1; ++j)
{
for (i__ = *m; i__ >= 1; --i__)
{
temp = *alpha * b_ref (i__, j);
i__2 = *m;
for (k = i__ + 1; k <= i__2; ++k)
{
temp -= a_ref (k, i__) * b_ref (k, j);
/* L140: */
}
if (nounit)
{
temp /= a_ref (i__, i__);
}
b_ref (i__, j) = temp;
/* L150: */
}
/* L160: */
}
}
}
}
else
{
if (lsame_ (transa, "N"))
{
/* Form B := alpha*B*inv( A ). */
if (upper)
{
i__1 = *n;
for (j = 1; j <= i__1; ++j)
{
if (*alpha != 1.)
{
i__2 = *m;
for (i__ = 1; i__ <= i__2; ++i__)
{
b_ref (i__, j) = *alpha * b_ref (i__, j);
/* L170: */
}
}
i__2 = j - 1;
for (k = 1; k <= i__2; ++k)
{
if (a_ref (k, j) != 0.)
{
i__3 = *m;
for (i__ = 1; i__ <= i__3; ++i__)
{
b_ref (i__, j) = b_ref (i__, j) - a_ref (k, j) * b_ref (i__, k);
/* L180: */
}
}
/* L190: */
}
if (nounit)
{
temp = 1. / a_ref (j, j);
i__2 = *m;
for (i__ = 1; i__ <= i__2; ++i__)
{
b_ref (i__, j) = temp * b_ref (i__, j);
/* L200: */
}
}
/* L210: */
}
}
else
{
for (j = *n; j >= 1; --j)
{
if (*alpha != 1.)
{
i__1 = *m;
for (i__ = 1; i__ <= i__1; ++i__)
{
b_ref (i__, j) = *alpha * b_ref (i__, j);
/* L220: */
}
}
i__1 = *n;
for (k = j + 1; k <= i__1; ++k)
{
if (a_ref (k, j) != 0.)
{
i__2 = *m;
for (i__ = 1; i__ <= i__2; ++i__)
{
b_ref (i__, j) = b_ref (i__, j) - a_ref (k, j) * b_ref (i__, k);
/* L230: */
}
}
/* L240: */
}
if (nounit)
{
temp = 1. / a_ref (j, j);
i__1 = *m;
for (i__ = 1; i__ <= i__1; ++i__)
{
b_ref (i__, j) = temp * b_ref (i__, j);
/* L250: */
}
}
/* L260: */
}
}
}
else
{
/* Form B := alpha*B*inv( A' ). */
if (upper)
{
for (k = *n; k >= 1; --k)
{
if (nounit)
{
temp = 1. / a_ref (k, k);
i__1 = *m;
for (i__ = 1; i__ <= i__1; ++i__)
{
b_ref (i__, k) = temp * b_ref (i__, k);
/* L270: */
}
}
i__1 = k - 1;
for (j = 1; j <= i__1; ++j)
{
if (a_ref (j, k) != 0.)
{
temp = a_ref (j, k);
i__2 = *m;
for (i__ = 1; i__ <= i__2; ++i__)
{
b_ref (i__, j) = b_ref (i__, j) - temp * b_ref (i__, k);
/* L280: */
}
}
/* L290: */
}
if (*alpha != 1.)
{
i__1 = *m;
for (i__ = 1; i__ <= i__1; ++i__)
{
b_ref (i__, k) = *alpha * b_ref (i__, k);
/* L300: */
}
}
/* L310: */
}
}
else
{
i__1 = *n;
for (k = 1; k <= i__1; ++k)
{
if (nounit)
{
temp = 1. / a_ref (k, k);
i__2 = *m;
for (i__ = 1; i__ <= i__2; ++i__)
{
b_ref (i__, k) = temp * b_ref (i__, k);
/* L320: */
}
}
i__2 = *n;
for (j = k + 1; j <= i__2; ++j)
{
if (a_ref (j, k) != 0.)
{
temp = a_ref (j, k);
i__3 = *m;
for (i__ = 1; i__ <= i__3; ++i__)
{
b_ref (i__, j) = b_ref (i__, j) - temp * b_ref (i__, k);
/* L330: */
}
}
/* L340: */
}
if (*alpha != 1.)
{
i__2 = *m;
for (i__ = 1; i__ <= i__2; ++i__)
{
b_ref (i__, k) = *alpha * b_ref (i__, k);
/* L350: */
}
}
/* L360: */
}
}
}
}
return 0;
} /* NUMblas_dtrsm */
#undef b_ref
#undef a_ref
long NUMblas_idamax (long *n, double *dx, long *incx)
{
/* System generated locals */
long ret_val, i__1;
double d__1;
/* Local variables */
static double dmax__;
static long i__, ix;
/* finds the index of element having max. absolute value. jack
dongarra, linpack, 3/11/78. modified 3/93 to return if incx .le. 0.
modified 12/3/93, array(1) declarations changed to array(*) Parameter
adjustments */
--dx;
/* Function Body */
ret_val = 0;
if (*n < 1 || *incx <= 0)
{
return ret_val;
}
ret_val = 1;
if (*n == 1)
{
return ret_val;
}
if (*incx == 1)
{
goto L20;
}
/* code for increment not equal to 1 */
ix = 1;
dmax__ = fabs (dx[1]);
ix += *incx;
i__1 = *n;
for (i__ = 2; i__ <= i__1; ++i__)
{
if ((d__1 = dx[ix], fabs (d__1)) <= dmax__)
{
goto L5;
}
ret_val = i__;
dmax__ = (d__1 = dx[ix], fabs (d__1));
L5:
ix += *incx;
/* L10: */
}
return ret_val;
/* code for increment equal to 1 */
L20:
dmax__ = fabs (dx[1]);
i__1 = *n;
for (i__ = 2; i__ <= i__1; ++i__)
{
if ((d__1 = dx[i__], fabs (d__1)) <= dmax__)
{
goto L30;
}
ret_val = i__;
dmax__ = (d__1 = dx[i__], fabs (d__1));
L30:
;
}
return ret_val;
} /* NUMblas_idamax */
#undef MAX
#undef MIN
/* End of file NUMcblas.c */