/* dtrmm.f -- translated by f2c (version of 23 April 1993 18:34:30).
You must link the resulting object file with the libraries:
-lf2c -lm (in that order)
*/
#include "f2c.h"
/* Subroutine */ int dtrmm_(side, uplo, transa, diag, m, n, alpha, a, lda, b,
ldb, side_len, uplo_len, transa_len, diag_len)
char *side, *uplo, *transa, *diag;
integer *m, *n;
doublereal *alpha, *a;
integer *lda;
doublereal *b;
integer *ldb;
ftnlen side_len;
ftnlen uplo_len;
ftnlen transa_len;
ftnlen diag_len;
{
/* System generated locals */
integer a_dim1, a_offset, b_dim1, b_offset, i__1, i__2, i__3;
/* Local variables */
static integer info;
static doublereal temp;
static integer i, j, k;
static logical lside;
extern logical lsame_();
static integer nrowa;
static logical upper;
extern /* Subroutine */ int xerbla_();
static logical nounit;
/* .. Scalar Arguments .. */
/* .. Array Arguments .. */
/* .. */
/* Purpose */
/* ======= */
/* DTRMM performs one of the matrix-matrix operations */
/* B := alpha*op( A )*B, or B := alpha*B*op( A ), */
/* where alpha is a scalar, B is an m by n matrix, A is a unit, or
*/
/* non-unit, upper or lower triangular matrix and op( A ) is one of
*/
/* op( A ) = A or op( A ) = A'. */
/* Parameters */
/* ========== */
/* SIDE - CHARACTER*1. */
/* On entry, SIDE specifies whether op( A ) multiplies B from
*/
/* the left or right as follows: */
/* SIDE = 'L' or 'l' B := alpha*op( A )*B. */
/* SIDE = 'R' or 'r' B := alpha*B*op( A ). */
/* Unchanged on exit. */
/* UPLO - CHARACTER*1. */
/* On entry, UPLO specifies whether the matrix A is an upper or
*/
/* lower triangular matrix as follows: */
/* UPLO = 'U' or 'u' A is an upper triangular matrix. */
/* UPLO = 'L' or 'l' A is a lower triangular matrix. */
/* Unchanged on exit. */
/* TRANSA - CHARACTER*1. */
/* On entry, TRANSA specifies the form of op( A ) to be used in
*/
/* the matrix multiplication as follows: */
/* TRANSA = 'N' or 'n' op( A ) = A. */
/* TRANSA = 'T' or 't' op( A ) = A'. */
/* TRANSA = 'C' or 'c' op( A ) = A'. */
/* Unchanged on exit. */
/* DIAG - CHARACTER*1. */
/* On entry, DIAG specifies whether or not A is unit triangular
*/
/* as follows: */
/* DIAG = 'U' or 'u' A is assumed to be unit triangular. */
/* DIAG = 'N' or 'n' A is not assumed to be unit */
/* triangular. */
/* Unchanged on exit. */
/* M - INTEGER. */
/* On entry, M specifies the number of rows of B. M must be at
*/
/* least zero. */
/* Unchanged on exit. */
/* N - INTEGER. */
/* On entry, N specifies the number of columns of B. N must be
*/
/* at least zero. */
/* Unchanged on exit. */
/* ALPHA - DOUBLE PRECISION. */
/* On entry, ALPHA specifies the scalar alpha. When alpha is
*/
/* zero then A is not referenced and B need not be set before
*/
/* entry. */
/* Unchanged on exit. */
/* A - DOUBLE PRECISION array of DIMENSION ( LDA, k ), where k is m
*/
/* when SIDE = 'L' or 'l' and is n when SIDE = 'R' or 'r'.
*/
/* Before entry with UPLO = 'U' or 'u', the leading k by k
*/
/* upper triangular part of the array A must contain the upper
*/
/* triangular matrix and the strictly lower triangular part of
*/
/* A is not referenced. */
/* Before entry with UPLO = 'L' or 'l', the leading k by k
*/
/* lower triangular part of the array A must contain the lower
*/
/* triangular matrix and the strictly upper triangular part of
*/
/* A is not referenced. */
/* Note that when DIAG = 'U' or 'u', the diagonal elements of
*/
/* A are not referenced either, but are assumed to be unity.
*/
/* Unchanged on exit. */
/* LDA - INTEGER. */
/* On entry, LDA specifies the first dimension of A as declared
*/
/* in the calling (sub) program. When SIDE = 'L' or 'l' then
*/
/* LDA must be at least max( 1, m ), when SIDE = 'R' or 'r'
*/
/* then LDA must be at least max( 1, n ). */
/* Unchanged on exit. */
/* B - DOUBLE PRECISION array of DIMENSION ( LDB, n ). */
/* Before entry, the leading m by n part of the array B must
*/
/* contain the matrix B, and on exit is overwritten by the
*/
/* transformed matrix. */
/* LDB - INTEGER. */
/* On entry, LDB specifies the first dimension of B as declared
*/
/* in the calling (sub) program. LDB must be at least
*/
/* max( 1, m ). */
/* Unchanged on exit. */
/* 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. */
/* .. External Functions .. */
/* .. External Subroutines .. */
/* .. Intrinsic Functions .. */
/* .. Local Scalars .. */
/* .. Parameters .. */
/* .. */
/* .. Executable Statements .. */
/* Test the input parameters. */
/* Parameter adjustments */
b_dim1 = *ldb;
b_offset = b_dim1 + 1;
b -= b_offset;
a_dim1 = *lda;
a_offset = a_dim1 + 1;
a -= a_offset;
/* Function Body */
lside = lsame_(side, "L", 1L, 1L);
if (lside) {
nrowa = *m;
} else {
nrowa = *n;
}
nounit = lsame_(diag, "N", 1L, 1L);
upper = lsame_(uplo, "U", 1L, 1L);
info = 0;
if (! lside && ! lsame_(side, "R", 1L, 1L)) {
info = 1;
} else if (! upper && ! lsame_(uplo, "L", 1L, 1L)) {
info = 2;
} else if (! lsame_(transa, "N", 1L, 1L) && ! lsame_(transa, "T", 1L, 1L)
&& ! lsame_(transa, "C", 1L, 1L)) {
info = 3;
} else if (! lsame_(diag, "U", 1L, 1L) && ! lsame_(diag, "N", 1L, 1L)) {
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, 6L);
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[i + j * b_dim1] = 0.;
/* L10: */
}
/* L20: */
}
return 0;
}
/* Start the operations. */
if (lside) {
if (lsame_(transa, "N", 1L, 1L)) {
/* 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[k + j * b_dim1] != 0.) {
temp = *alpha * b[k + j * b_dim1];
i__3 = k - 1;
for (i = 1; i <= i__3; ++i) {
b[i + j * b_dim1] += temp * a[i + k * a_dim1];
/* L30: */
}
if (nounit) {
temp *= a[k + k * a_dim1];
}
b[k + j * b_dim1] = temp;
}
/* L40: */
}
/* L50: */
}
} else {
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
for (k = *m; k >= 1; --k) {
if (b[k + j * b_dim1] != 0.) {
temp = *alpha * b[k + j * b_dim1];
b[k + j * b_dim1] = temp;
if (nounit) {
b[k + j * b_dim1] *= a[k + k * a_dim1];
}
i__2 = *m;
for (i = k + 1; i <= i__2; ++i) {
b[i + j * b_dim1] += temp * a[i + k * a_dim1];
/* L60: */
}
}
/* L70: */
}
/* L80: */
}
}
} else {
/* Form B := alpha*B*A'. */
if (upper) {
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
for (i = *m; i >= 1; --i) {
temp = b[i + j * b_dim1];
if (nounit) {
temp *= a[i + i * a_dim1];
}
i__2 = i - 1;
for (k = 1; k <= i__2; ++k) {
temp += a[k + i * a_dim1] * b[k + j * b_dim1];
/* L90: */
}
b[i + j * b_dim1] = *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[i + j * b_dim1];
if (nounit) {
temp *= a[i + i * a_dim1];
}
i__3 = *m;
for (k = i + 1; k <= i__3; ++k) {
temp += a[k + i * a_dim1] * b[k + j * b_dim1];
/* L120: */
}
b[i + j * b_dim1] = *alpha * temp;
/* L130: */
}
/* L140: */
}
}
}
} else {
if (lsame_(transa, "N", 1L, 1L)) {
/* Form B := alpha*B*A. */
if (upper) {
for (j = *n; j >= 1; --j) {
temp = *alpha;
if (nounit) {
temp *= a[j + j * a_dim1];
}
i__1 = *m;
for (i = 1; i <= i__1; ++i) {
b[i + j * b_dim1] = temp * b[i + j * b_dim1];
/* L150: */
}
i__1 = j - 1;
for (k = 1; k <= i__1; ++k) {
if (a[k + j * a_dim1] != 0.) {
temp = *alpha * a[k + j * a_dim1];
i__2 = *m;
for (i = 1; i <= i__2; ++i) {
b[i + j * b_dim1] += temp * b[i + k * b_dim1];
/* L160: */
}
}
/* L170: */
}
/* L180: */
}
} else {
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
temp = *alpha;
if (nounit) {
temp *= a[j + j * a_dim1];
}
i__2 = *m;
for (i = 1; i <= i__2; ++i) {
b[i + j * b_dim1] = temp * b[i + j * b_dim1];
/* L190: */
}
i__2 = *n;
for (k = j + 1; k <= i__2; ++k) {
if (a[k + j * a_dim1] != 0.) {
temp = *alpha * a[k + j * a_dim1];
i__3 = *m;
for (i = 1; i <= i__3; ++i) {
b[i + j * b_dim1] += temp * b[i + k * b_dim1];
/* 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[j + k * a_dim1] != 0.) {
temp = *alpha * a[j + k * a_dim1];
i__3 = *m;
for (i = 1; i <= i__3; ++i) {
b[i + j * b_dim1] += temp * b[i + k * b_dim1];
/* L230: */
}
}
/* L240: */
}
temp = *alpha;
if (nounit) {
temp *= a[k + k * a_dim1];
}
if (temp != 1.) {
i__2 = *m;
for (i = 1; i <= i__2; ++i) {
b[i + k * b_dim1] = temp * b[i + k * b_dim1];
/* L250: */
}
}
/* L260: */
}
} else {
for (k = *n; k >= 1; --k) {
i__1 = *n;
for (j = k + 1; j <= i__1; ++j) {
if (a[j + k * a_dim1] != 0.) {
temp = *alpha * a[j + k * a_dim1];
i__2 = *m;
for (i = 1; i <= i__2; ++i) {
b[i + j * b_dim1] += temp * b[i + k * b_dim1];
/* L270: */
}
}
/* L280: */
}
temp = *alpha;
if (nounit) {
temp *= a[k + k * a_dim1];
}
if (temp != 1.) {
i__1 = *m;
for (i = 1; i <= i__1; ++i) {
b[i + k * b_dim1] = temp * b[i + k * b_dim1];
/* L290: */
}
}
/* L300: */
}
}
}
}
return 0;
/* End of DTRMM . */
} /* dtrmm_ */
/* idamax.f -- translated by f2c (version of 23 April 1993 18:34:30).
You must link the resulting object file with the libraries:
-lf2c -lm (in that order)
*/
integer idamax_(n, dx, incx)
integer *n;
doublereal *dx;
integer *incx;
{
/* System generated locals */
integer ret_val, i__1;
doublereal d__1;
/* Local variables */
static doublereal dmax_;
static integer i, ix;
/* finds the index of element having max. absolute value. */
/* jack dongarra, linpack, 3/11/78. */
/* modified to correct problem with negative increment, 8/21/90. */
/* Parameter adjustments */
--dx;
/* Function Body */
ret_val = 0;
if (*n < 1) {
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;
if (*incx < 0) {
ix = (-(*n) + 1) * *incx + 1;
}
dmax_ = (d__1 = dx[ix], abs(d__1));
ix += *incx;
i__1 = *n;
for (i = 2; i <= i__1; ++i) {
if ((d__1 = dx[ix], abs(d__1)) <= dmax_) {
goto L5;
}
ret_val = i;
dmax_ = (d__1 = dx[ix], abs(d__1));
L5:
ix += *incx;
/* L10: */
}
return ret_val;
/* code for increment equal to 1 */
L20:
dmax_ = abs(dx[1]);
i__1 = *n;
for (i = 2; i <= i__1; ++i) {
if ((d__1 = dx[i], abs(d__1)) <= dmax_) {
goto L30;
}
ret_val = i;
dmax_ = (d__1 = dx[i], abs(d__1));
L30:
;
}
return ret_val;
} /* idamax_ */
/* drot.f -- translated by f2c (version of 23 April 1993 18:34:30).
You must link the resulting object file with the libraries:
-lf2c -lm (in that order)
*/
/* Subroutine */ int drot_(n, dx, incx, dy, incy, c, s)
integer *n;
doublereal *dx;
integer *incx;
doublereal *dy;
integer *incy;
doublereal *c, *s;
{
/* System generated locals */
integer i__1;
/* Local variables */
static integer i;
static doublereal dtemp;
static integer ix, iy;
/* applies a plane rotation. */
/* jack dongarra, linpack, 3/11/78. */
/* 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;
} /* drot_ */
/* dtrsm.f -- translated by f2c (version of 23 April 1993 18:34:30).
You must link the resulting object file with the libraries:
-lf2c -lm (in that order)
*/
/* Subroutine */ int dtrsm_(side, uplo, transa, diag, m, n, alpha, a, lda, b,
ldb, side_len, uplo_len, transa_len, diag_len)
char *side, *uplo, *transa, *diag;
integer *m, *n;
doublereal *alpha, *a;
integer *lda;
doublereal *b;
integer *ldb;
ftnlen side_len;
ftnlen uplo_len;
ftnlen transa_len;
ftnlen diag_len;
{
/* System generated locals */
integer a_dim1, a_offset, b_dim1, b_offset, i__1, i__2, i__3;
/* Local variables */
static integer info;
static doublereal temp;
static integer i, j, k;
static logical lside;
extern logical lsame_();
static integer nrowa;
static logical upper;
extern /* Subroutine */ int xerbla_();
static logical nounit;
/* .. Scalar Arguments .. */
/* .. Array Arguments .. */
/* .. */
/* Purpose */
/* ======= */
/* DTRSM solves one of the matrix equations */
/* op( A )*X = alpha*B, or X*op( A ) = alpha*B, */
/* where alpha is a scalar, X and B are m by n matrices, A is a unit, or
*/
/* non-unit, upper or lower triangular matrix and op( A ) is one of
*/
/* op( A ) = A or op( A ) = A'. */
/* The matrix X is overwritten on B. */
/* Parameters */
/* ========== */
/* SIDE - CHARACTER*1. */
/* On entry, SIDE specifies whether op( A ) appears on the left
*/
/* or right of X as follows: */
/* SIDE = 'L' or 'l' op( A )*X = alpha*B. */
/* SIDE = 'R' or 'r' X*op( A ) = alpha*B. */
/* Unchanged on exit. */
/* UPLO - CHARACTER*1. */
/* On entry, UPLO specifies whether the matrix A is an upper or
*/
/* lower triangular matrix as follows: */
/* UPLO = 'U' or 'u' A is an upper triangular matrix. */
/* UPLO = 'L' or 'l' A is a lower triangular matrix. */
/* Unchanged on exit. */
/* TRANSA - CHARACTER*1. */
/* On entry, TRANSA specifies the form of op( A ) to be used in
*/
/* the matrix multiplication as follows: */
/* TRANSA = 'N' or 'n' op( A ) = A. */
/* TRANSA = 'T' or 't' op( A ) = A'. */
/* TRANSA = 'C' or 'c' op( A ) = A'. */
/* Unchanged on exit. */
/* DIAG - CHARACTER*1. */
/* On entry, DIAG specifies whether or not A is unit triangular
*/
/* as follows: */
/* DIAG = 'U' or 'u' A is assumed to be unit triangular. */
/* DIAG = 'N' or 'n' A is not assumed to be unit */
/* triangular. */
/* Unchanged on exit. */
/* M - INTEGER. */
/* On entry, M specifies the number of rows of B. M must be at
*/
/* least zero. */
/* Unchanged on exit. */
/* N - INTEGER. */
/* On entry, N specifies the number of columns of B. N must be
*/
/* at least zero. */
/* Unchanged on exit. */
/* ALPHA - DOUBLE PRECISION. */
/* On entry, ALPHA specifies the scalar alpha. When alpha is
*/
/* zero then A is not referenced and B need not be set before
*/
/* entry. */
/* Unchanged on exit. */
/* A - DOUBLE PRECISION array of DIMENSION ( LDA, k ), where k is m
*/
/* when SIDE = 'L' or 'l' and is n when SIDE = 'R' or 'r'.
*/
/* Before entry with UPLO = 'U' or 'u', the leading k by k
*/
/* upper triangular part of the array A must contain the upper
*/
/* triangular matrix and the strictly lower triangular part of
*/
/* A is not referenced. */
/* Before entry with UPLO = 'L' or 'l', the leading k by k
*/
/* lower triangular part of the array A must contain the lower
*/
/* triangular matrix and the strictly upper triangular part of
*/
/* A is not referenced. */
/* Note that when DIAG = 'U' or 'u', the diagonal elements of
*/
/* A are not referenced either, but are assumed to be unity.
*/
/* Unchanged on exit. */
/* LDA - INTEGER. */
/* On entry, LDA specifies the first dimension of A as declared
*/
/* in the calling (sub) program. When SIDE = 'L' or 'l' then
*/
/* LDA must be at least max( 1, m ), when SIDE = 'R' or 'r'
*/
/* then LDA must be at least max( 1, n ). */
/* Unchanged on exit. */
/* B - DOUBLE PRECISION array of DIMENSION ( LDB, n ). */
/* Before entry, the leading m by n part of the array B must
*/
/* contain the right-hand side matrix B, and on exit is
*/
/* overwritten by the solution matrix X. */
/* LDB - INTEGER. */
/* On entry, LDB specifies the first dimension of B as declared
*/
/* in the calling (sub) program. LDB must be at least
*/
/* max( 1, m ). */
/* Unchanged on exit. */
/* 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. */
/* .. External Functions .. */
/* .. External Subroutines .. */
/* .. Intrinsic Functions .. */
/* .. Local Scalars .. */
/* .. Parameters .. */
/* .. */
/* .. Executable Statements .. */
/* Test the input parameters. */
/* Parameter adjustments */
b_dim1 = *ldb;
b_offset = b_dim1 + 1;
b -= b_offset;
a_dim1 = *lda;
a_offset = a_dim1 + 1;
a -= a_offset;
/* Function Body */
lside = lsame_(side, "L", 1L, 1L);
if (lside) {
nrowa = *m;
} else {
nrowa = *n;
}
nounit = lsame_(diag, "N", 1L, 1L);
upper = lsame_(uplo, "U", 1L, 1L);
info = 0;
if (! lside && ! lsame_(side, "R", 1L, 1L)) {
info = 1;
} else if (! upper && ! lsame_(uplo, "L", 1L, 1L)) {
info = 2;
} else if (! lsame_(transa, "N", 1L, 1L) && ! lsame_(transa, "T", 1L, 1L)
&& ! lsame_(transa, "C", 1L, 1L)) {
info = 3;
} else if (! lsame_(diag, "U", 1L, 1L) && ! lsame_(diag, "N", 1L, 1L)) {
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, 6L);
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[i + j * b_dim1] = 0.;
/* L10: */
}
/* L20: */
}
return 0;
}
/* Start the operations. */
if (lside) {
if (lsame_(transa, "N", 1L, 1L)) {
/* 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[i + j * b_dim1] = *alpha * b[i + j * b_dim1];
/* L30: */
}
}
for (k = *m; k >= 1; --k) {
if (b[k + j * b_dim1] != 0.) {
if (nounit) {
b[k + j * b_dim1] /= a[k + k * a_dim1];
}
i__2 = k - 1;
for (i = 1; i <= i__2; ++i) {
b[i + j * b_dim1] -= b[k + j * b_dim1] * a[i
+ k * a_dim1];
/* 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[i + j * b_dim1] = *alpha * b[i + j * b_dim1];
/* L70: */
}
}
i__2 = *m;
for (k = 1; k <= i__2; ++k) {
if (b[k + j * b_dim1] != 0.) {
if (nounit) {
b[k + j * b_dim1] /= a[k + k * a_dim1];
}
i__3 = *m;
for (i = k + 1; i <= i__3; ++i) {
b[i + j * b_dim1] -= b[k + j * b_dim1] * a[i
+ k * a_dim1];
/* 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[i + j * b_dim1];
i__3 = i - 1;
for (k = 1; k <= i__3; ++k) {
temp -= a[k + i * a_dim1] * b[k + j * b_dim1];
/* L110: */
}
if (nounit) {
temp /= a[i + i * a_dim1];
}
b[i + j * b_dim1] = temp;
/* L120: */
}
/* L130: */
}
} else {
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
for (i = *m; i >= 1; --i) {
temp = *alpha * b[i + j * b_dim1];
i__2 = *m;
for (k = i + 1; k <= i__2; ++k) {
temp -= a[k + i * a_dim1] * b[k + j * b_dim1];
/* L140: */
}
if (nounit) {
temp /= a[i + i * a_dim1];
}
b[i + j * b_dim1] = temp;
/* L150: */
}
/* L160: */
}
}
}
} else {
if (lsame_(transa, "N", 1L, 1L)) {
/* 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[i + j * b_dim1] = *alpha * b[i + j * b_dim1];
/* L170: */
}
}
i__2 = j - 1;
for (k = 1; k <= i__2; ++k) {
if (a[k + j * a_dim1] != 0.) {
i__3 = *m;
for (i = 1; i <= i__3; ++i) {
b[i + j * b_dim1] -= a[k + j * a_dim1] * b[i
+ k * b_dim1];
/* L180: */
}
}
/* L190: */
}
if (nounit) {
temp = 1. / a[j + j * a_dim1];
i__2 = *m;
for (i = 1; i <= i__2; ++i) {
b[i + j * b_dim1] = temp * b[i + j * b_dim1];
/* L200: */
}
}
/* L210: */
}
} else {
for (j = *n; j >= 1; --j) {
if (*alpha != 1.) {
i__1 = *m;
for (i = 1; i <= i__1; ++i) {
b[i + j * b_dim1] = *alpha * b[i + j * b_dim1];
/* L220: */
}
}
i__1 = *n;
for (k = j + 1; k <= i__1; ++k) {
if (a[k + j * a_dim1] != 0.) {
i__2 = *m;
for (i = 1; i <= i__2; ++i) {
b[i + j * b_dim1] -= a[k + j * a_dim1] * b[i
+ k * b_dim1];
/* L230: */
}
}
/* L240: */
}
if (nounit) {
temp = 1. / a[j + j * a_dim1];
i__1 = *m;
for (i = 1; i <= i__1; ++i) {
b[i + j * b_dim1] = temp * b[i + j * b_dim1];
/* L250: */
}
}
/* L260: */
}
}
} else {
/* Form B := alpha*B*inv( A' ). */
if (upper) {
for (k = *n; k >= 1; --k) {
if (nounit) {
temp = 1. / a[k + k * a_dim1];
i__1 = *m;
for (i = 1; i <= i__1; ++i) {
b[i + k * b_dim1] = temp * b[i + k * b_dim1];
/* L270: */
}
}
i__1 = k - 1;
for (j = 1; j <= i__1; ++j) {
if (a[j + k * a_dim1] != 0.) {
temp = a[j + k * a_dim1];
i__2 = *m;
for (i = 1; i <= i__2; ++i) {
b[i + j * b_dim1] -= temp * b[i + k * b_dim1];
/* L280: */
}
}
/* L290: */
}
if (*alpha != 1.) {
i__1 = *m;
for (i = 1; i <= i__1; ++i) {
b[i + k * b_dim1] = *alpha * b[i + k * b_dim1];
/* L300: */
}
}
/* L310: */
}
} else {
i__1 = *n;
for (k = 1; k <= i__1; ++k) {
if (nounit) {
temp = 1. / a[k + k * a_dim1];
i__2 = *m;
for (i = 1; i <= i__2; ++i) {
b[i + k * b_dim1] = temp * b[i + k * b_dim1];
/* L320: */
}
}
i__2 = *n;
for (j = k + 1; j <= i__2; ++j) {
if (a[j + k * a_dim1] != 0.) {
temp = a[j + k * a_dim1];
i__3 = *m;
for (i = 1; i <= i__3; ++i) {
b[i + j * b_dim1] -= temp * b[i + k * b_dim1];
/* L330: */
}
}
/* L340: */
}
if (*alpha != 1.) {
i__2 = *m;
for (i = 1; i <= i__2; ++i) {
b[i + k * b_dim1] = *alpha * b[i + k * b_dim1];
/* L350: */
}
}
/* L360: */
}
}
}
}
return 0;
/* End of DTRSM . */
} /* dtrsm_ */
/* dger.f -- translated by f2c (version of 23 April 1993 18:34:30).
You must link the resulting object file with the libraries:
-lf2c -lm (in that order)
*/
/* Subroutine */ int dger_(m, n, alpha, x, incx, y, incy, a, lda)
integer *m, *n;
doublereal *alpha, *x;
integer *incx;
doublereal *y;
integer *incy;
doublereal *a;
integer *lda;
{
/* System generated locals */
integer a_dim1, a_offset, i__1, i__2;
/* Local variables */
static integer info;
static doublereal temp;
static integer i, j, ix, jy, kx;
extern /* Subroutine */ int xerbla_();
/* .. Scalar Arguments .. */
/* .. Array Arguments .. */
/* .. */
/* Purpose */
/* ======= */
/* DGER performs the rank 1 operation */
/* A := alpha*x*y' + A, */
/* where alpha is a scalar, x is an m element vector, y is an n element
*/
/* vector and A is an m by n matrix. */
/* Parameters */
/* ========== */
/* M - INTEGER. */
/* On entry, M specifies the number of rows of the matrix A. */
/* M must be at least zero. */
/* Unchanged on exit. */
/* N - INTEGER. */
/* On entry, N specifies the number of columns of the matrix A.
*/
/* N must be at least zero. */
/* Unchanged on exit. */
/* ALPHA - DOUBLE PRECISION. */
/* On entry, ALPHA specifies the scalar alpha. */
/* Unchanged on exit. */
/* X - DOUBLE PRECISION array of dimension at least */
/* ( 1 + ( m - 1 )*abs( INCX ) ). */
/* Before entry, the incremented array X must contain the m */
/* element vector x. */
/* Unchanged on exit. */
/* INCX - INTEGER. */
/* On entry, INCX specifies the increment for the elements of */
/* X. INCX must not be zero. */
/* Unchanged on exit. */
/* Y - DOUBLE PRECISION array of dimension at least */
/* ( 1 + ( n - 1 )*abs( INCY ) ). */
/* Before entry, the incremented array Y must contain the n */
/* element vector y. */
/* Unchanged on exit. */
/* INCY - INTEGER. */
/* On entry, INCY specifies the increment for the elements of */
/* Y. INCY must not be zero. */
/* Unchanged on exit. */
/* A - DOUBLE PRECISION array of DIMENSION ( LDA, n ). */
/* Before entry, the leading m by n part of the array A must */
/* contain the matrix of coefficients. On exit, A is */
/* overwritten by the updated matrix. */
/* LDA - INTEGER. */
/* On entry, LDA specifies the first dimension of A as declared
*/
/* in the calling (sub) program. LDA must be at least */
/* max( 1, m ). */
/* Unchanged on exit. */
/* Level 2 Blas routine. */
/* -- 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. */
/* .. Parameters .. */
/* .. Local Scalars .. */
/* .. External Subroutines .. */
/* .. Intrinsic Functions .. */
/* .. */
/* .. Executable Statements .. */
/* Test the input parameters. */
/* Parameter adjustments */
a_dim1 = *lda;
a_offset = a_dim1 + 1;
a -= a_offset;
--y;
--x;
/* 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, 6L);
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[i + j * a_dim1] += 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[i + j * a_dim1] += x[ix] * temp;
ix += *incx;
/* L30: */
}
}
jy += *incy;
/* L40: */
}
}
return 0;
/* End of DGER . */
} /* dger_ */
/* dcopy.f -- translated by f2c (version of 23 April 1993 18:34:30).
You must link the resulting object file with the libraries:
-lf2c -lm (in that order)
*/
/* Subroutine */ int dcopy_(n, dx, incx, dy, incy)
integer *n;
doublereal *dx;
integer *incx;
doublereal *dy;
integer *incy;
{
/* System generated locals */
integer i__1;
/* Local variables */
static integer i, m, ix, iy, mp1;
/* copies a vector, x, to a vector, y. */
/* uses unrolled loops for increments equal to one. */
/* jack dongarra, linpack, 3/11/78. */
/* 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) {
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;
} /* dcopy_ */
/* zswap.f -- translated by f2c (version of 23 April 1993 18:34:30).
You must link the resulting object file with the libraries:
-lf2c -lm (in that order)
*/
/* Subroutine */ int zswap_(n, zx, incx, zy, incy)
integer *n;
doublecomplex *zx;
integer *incx;
doublecomplex *zy;
integer *incy;
{
/* System generated locals */
integer i__1, i__2, i__3;
/* Local variables */
static integer i;
static doublecomplex ztemp;
static integer ix, iy;
/* interchanges two vectors. */
/* jack dongarra, 3/11/78. */
/* Parameter adjustments */
--zy;
--zx;
/* 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) {
i__2 = ix;
ztemp.r = zx[i__2].r, ztemp.i = zx[i__2].i;
i__2 = ix;
i__3 = iy;
zx[i__2].r = zy[i__3].r, zx[i__2].i = zy[i__3].i;
i__2 = iy;
zy[i__2].r = ztemp.r, zy[i__2].i = ztemp.i;
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) {
i__2 = i;
ztemp.r = zx[i__2].r, ztemp.i = zx[i__2].i;
i__2 = i;
i__3 = i;
zx[i__2].r = zy[i__3].r, zx[i__2].i = zy[i__3].i;
i__2 = i;
zy[i__2].r = ztemp.r, zy[i__2].i = ztemp.i;
/* L30: */
}
return 0;
} /* zswap_ */
/* ddot.f -- translated by f2c (version of 23 April 1993 18:34:30).
You must link the resulting object file with the libraries:
-lf2c -lm (in that order)
*/
doublereal ddot_(n, dx, incx, dy, incy)
integer *n;
doublereal *dx;
integer *incx;
doublereal *dy;
integer *incy;
{
/* System generated locals */
integer i__1;
doublereal ret_val;
/* Local variables */
static integer i, m;
static doublereal dtemp;
static integer ix, iy, mp1;
/* forms the dot product of two vectors. */
/* uses unrolled loops for increments equal to one. */
/* jack dongarra, linpack, 3/11/78. */
/* 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;
} /* ddot_ */
/* zscal.f -- translated by f2c (version of 23 April 1993 18:34:30).
You must link the resulting object file with the libraries:
-lf2c -lm (in that order)
*/
/* Subroutine */ int zscal_(n, za, zx, incx)
integer *n;
doublecomplex *za, *zx;
integer *incx;
{
/* System generated locals */
integer i__1, i__2, i__3;
doublecomplex z__1;
/* Local variables */
static integer i, ix;
/* scales a vector by a constant. */
/* jack dongarra, 3/11/78. */
/* modified to correct problem with negative increment, 8/21/90. */
/* Parameter adjustments */
--zx;
/* Function Body */
if (*n <= 0) {
return 0;
}
if (*incx == 1) {
goto L20;
}
/* code for increment not equal to 1 */
ix = 1;
if (*incx < 0) {
ix = (-(*n) + 1) * *incx + 1;
}
i__1 = *n;
for (i = 1; i <= i__1; ++i) {
i__2 = ix;
i__3 = ix;
z__1.r = za->r * zx[i__3].r - za->i * zx[i__3].i, z__1.i = za->r * zx[
i__3].i + za->i * zx[i__3].r;
zx[i__2].r = z__1.r, zx[i__2].i = z__1.i;
ix += *incx;
/* L10: */
}
return 0;
/* code for increment equal to 1 */
L20:
i__1 = *n;
for (i = 1; i <= i__1; ++i) {
i__2 = i;
i__3 = i;
z__1.r = za->r * zx[i__3].r - za->i * zx[i__3].i, z__1.i = za->r * zx[
i__3].i + za->i * zx[i__3].r;
zx[i__2].r = z__1.r, zx[i__2].i = z__1.i;
/* L30: */
}
return 0;
} /* zscal_ */
/* dcabs1.f -- translated by f2c (version of 23 April 1993 18:34:30).
You must link the resulting object file with the libraries:
-lf2c -lm (in that order)
*/
doublereal dcabs1_(z)
doublecomplex *z;
{
/* System generated locals */
doublereal ret_val;
static doublecomplex equiv_0[1];
/* Local variables */
#define t ((doublereal *)equiv_0)
#define zz (equiv_0)
zz->r = z->r, zz->i = z->i;
ret_val = abs(t[0]) + abs(t[1]);
return ret_val;
} /* dcabs1_ */
#undef zz
#undef t
/* ztrsv.f -- translated by f2c (version of 23 April 1993 18:34:30).
You must link the resulting object file with the libraries:
-lf2c -lm (in that order)
*/
/* Subroutine */ int ztrsv_(uplo, trans, diag, n, a, lda, x, incx, uplo_len,
trans_len, diag_len)
char *uplo, *trans, *diag;
integer *n;
doublecomplex *a;
integer *lda;
doublecomplex *x;
integer *incx;
ftnlen uplo_len;
ftnlen trans_len;
ftnlen diag_len;
{
/* System generated locals */
integer a_dim1, a_offset, i__1, i__2, i__3, i__4, i__5;
doublecomplex z__1, z__2, z__3;
/* Builtin functions */
void z_div(), d_cnjg();
/* Local variables */
static integer info;
static doublecomplex temp;
static integer i, j;
extern logical lsame_();
static integer ix, jx, kx;
extern /* Subroutine */ int xerbla_();
static logical noconj, nounit;
/* .. Scalar Arguments .. */
/* .. Array Arguments .. */
/* .. */
/* Purpose */
/* ======= */
/* ZTRSV solves one of the systems of equations */
/* A*x = b, or A'*x = b, or conjg( A' )*x = b, */
/* where b and x are n element vectors and A is an n by n unit, or */
/* non-unit, upper or lower triangular matrix. */
/* No test for singularity or near-singularity is included in this */
/* routine. Such tests must be performed before calling this routine. */
/* Parameters */
/* ========== */
/* UPLO - CHARACTER*1. */
/* On entry, UPLO specifies whether the matrix is an upper or */
/* lower triangular matrix as follows: */
/* UPLO = 'U' or 'u' A is an upper triangular matrix. */
/* UPLO = 'L' or 'l' A is a lower triangular matrix. */
/* Unchanged on exit. */
/* TRANS - CHARACTER*1. */
/* On entry, TRANS specifies the equations to be solved as */
/* follows: */
/* TRANS = 'N' or 'n' A*x = b. */
/* TRANS = 'T' or 't' A'*x = b. */
/* TRANS = 'C' or 'c' conjg( A' )*x = b. */
/* Unchanged on exit. */
/* DIAG - CHARACTER*1. */
/* On entry, DIAG specifies whether or not A is unit */
/* triangular as follows: */
/* DIAG = 'U' or 'u' A is assumed to be unit triangular. */
/* DIAG = 'N' or 'n' A is not assumed to be unit */
/* triangular. */
/* Unchanged on exit. */
/* N - INTEGER. */
/* On entry, N specifies the order of the matrix A. */
/* N must be at least zero. */
/* Unchanged on exit. */
/* A - COMPLEX*16 array of DIMENSION ( LDA, n ). */
/* Before entry with UPLO = 'U' or 'u', the leading n by n */
/* upper triangular part of the array A must contain the upper
*/
/* triangular matrix and the strictly lower triangular part of
*/
/* A is not referenced. */
/* Before entry with UPLO = 'L' or 'l', the leading n by n */
/* lower triangular part of the array A must contain the lower
*/
/* triangular matrix and the strictly upper triangular part of
*/
/* A is not referenced. */
/* Note that when DIAG = 'U' or 'u', the diagonal elements of
*/
/* A are not referenced either, but are assumed to be unity. */
/* Unchanged on exit. */
/* LDA - INTEGER. */
/* On entry, LDA specifies the first dimension of A as declared
*/
/* in the calling (sub) program. LDA must be at least */
/* max( 1, n ). */
/* Unchanged on exit. */
/* X - COMPLEX*16 array of dimension at least */
/* ( 1 + ( n - 1 )*abs( INCX ) ). */
/* Before entry, the incremented array X must contain the n */
/* element right-hand side vector b. On exit, X is overwritten
*/
/* with the solution vector x. */
/* INCX - INTEGER. */
/* On entry, INCX specifies the increment for the elements of */
/* X. INCX must not be zero. */
/* Unchanged on exit. */
/* Level 2 Blas routine. */
/* -- 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. */
/* .. Parameters .. */
/* .. Local Scalars .. */
/* .. External Functions .. */
/* .. External Subroutines .. */
/* .. Intrinsic Functions .. */
/* .. */
/* .. Executable Statements .. */
/* Test the input parameters. */
/* Parameter adjustments */
--x;
a_dim1 = *lda;
a_offset = a_dim1 + 1;
a -= a_offset;
/* Function Body */
info = 0;
if (! lsame_(uplo, "U", 1L, 1L) && ! lsame_(uplo, "L", 1L, 1L)) {
info = 1;
} else if (! lsame_(trans, "N", 1L, 1L) && ! lsame_(trans, "T", 1L, 1L) &&
! lsame_(trans, "C", 1L, 1L)) {
info = 2;
} else if (! lsame_(diag, "U", 1L, 1L) && ! lsame_(diag, "N", 1L, 1L)) {
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_("ZTRSV ", &info, 6L);
return 0;
}
/* Quick return if possible. */
if (*n == 0) {
return 0;
}
noconj = lsame_(trans, "T", 1L, 1L);
nounit = lsame_(diag, "N", 1L, 1L);
/* 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", 1L, 1L)) {
/* Form x := inv( A )*x. */
if (lsame_(uplo, "U", 1L, 1L)) {
if (*incx == 1) {
for (j = *n; j >= 1; --j) {
i__1 = j;
if (x[i__1].r != 0. || x[i__1].i != 0.) {
if (nounit) {
i__1 = j;
z_div(&z__1, &x[j], &a[j + j * a_dim1]);
x[i__1].r = z__1.r, x[i__1].i = z__1.i;
}
i__1 = j;
temp.r = x[i__1].r, temp.i = x[i__1].i;
for (i = j - 1; i >= 1; --i) {
i__1 = i;
i__2 = i;
i__3 = i + j * a_dim1;
z__2.r = temp.r * a[i__3].r - temp.i * a[i__3].i,
z__2.i = temp.r * a[i__3].i + temp.i * a[
i__3].r;
z__1.r = x[i__2].r - z__2.r, z__1.i = x[i__2].i -
z__2.i;
x[i__1].r = z__1.r, x[i__1].i = z__1.i;
/* L10: */
}
}
/* L20: */
}
} else {
jx = kx + (*n - 1) * *incx;
for (j = *n; j >= 1; --j) {
i__1 = jx;
if (x[i__1].r != 0. || x[i__1].i != 0.) {
if (nounit) {
i__1 = jx;
z_div(&z__1, &x[jx], &a[j + j * a_dim1]);
x[i__1].r = z__1.r, x[i__1].i = z__1.i;
}
i__1 = jx;
temp.r = x[i__1].r, temp.i = x[i__1].i;
ix = jx;
for (i = j - 1; i >= 1; --i) {
ix -= *incx;
i__1 = ix;
i__2 = ix;
i__3 = i + j * a_dim1;
z__2.r = temp.r * a[i__3].r - temp.i * a[i__3].i,
z__2.i = temp.r * a[i__3].i + temp.i * a[
i__3].r;
z__1.r = x[i__2].r - z__2.r, z__1.i = x[i__2].i -
z__2.i;
x[i__1].r = z__1.r, x[i__1].i = z__1.i;
/* L30: */
}
}
jx -= *incx;
/* L40: */
}
}
} else {
if (*incx == 1) {
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
i__2 = j;
if (x[i__2].r != 0. || x[i__2].i != 0.) {
if (nounit) {
i__2 = j;
z_div(&z__1, &x[j], &a[j + j * a_dim1]);
x[i__2].r = z__1.r, x[i__2].i = z__1.i;
}
i__2 = j;
temp.r = x[i__2].r, temp.i = x[i__2].i;
i__2 = *n;
for (i = j + 1; i <= i__2; ++i) {
i__3 = i;
i__4 = i;
i__5 = i + j * a_dim1;
z__2.r = temp.r * a[i__5].r - temp.i * a[i__5].i,
z__2.i = temp.r * a[i__5].i + temp.i * a[
i__5].r;
z__1.r = x[i__4].r - z__2.r, z__1.i = x[i__4].i -
z__2.i;
x[i__3].r = z__1.r, x[i__3].i = z__1.i;
/* L50: */
}
}
/* L60: */
}
} else {
jx = kx;
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
i__2 = jx;
if (x[i__2].r != 0. || x[i__2].i != 0.) {
if (nounit) {
i__2 = jx;
z_div(&z__1, &x[jx], &a[j + j * a_dim1]);
x[i__2].r = z__1.r, x[i__2].i = z__1.i;
}
i__2 = jx;
temp.r = x[i__2].r, temp.i = x[i__2].i;
ix = jx;
i__2 = *n;
for (i = j + 1; i <= i__2; ++i) {
ix += *incx;
i__3 = ix;
i__4 = ix;
i__5 = i + j * a_dim1;
z__2.r = temp.r * a[i__5].r - temp.i * a[i__5].i,
z__2.i = temp.r * a[i__5].i + temp.i * a[
i__5].r;
z__1.r = x[i__4].r - z__2.r, z__1.i = x[i__4].i -
z__2.i;
x[i__3].r = z__1.r, x[i__3].i = z__1.i;
/* L70: */
}
}
jx += *incx;
/* L80: */
}
}
}
} else {
/* Form x := inv( A' )*x or x := inv( conjg( A' ) )*x. */
if (lsame_(uplo, "U", 1L, 1L)) {
if (*incx == 1) {
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
i__2 = j;
temp.r = x[i__2].r, temp.i = x[i__2].i;
if (noconj) {
i__2 = j - 1;
for (i = 1; i <= i__2; ++i) {
i__3 = i + j * a_dim1;
i__4 = i;
z__2.r = a[i__3].r * x[i__4].r - a[i__3].i * x[
i__4].i, z__2.i = a[i__3].r * x[i__4].i +
a[i__3].i * x[i__4].r;
z__1.r = temp.r - z__2.r, z__1.i = temp.i -
z__2.i;
temp.r = z__1.r, temp.i = z__1.i;
/* L90: */
}
if (nounit) {
z_div(&z__1, &temp, &a[j + j * a_dim1]);
temp.r = z__1.r, temp.i = z__1.i;
}
} else {
i__2 = j - 1;
for (i = 1; i <= i__2; ++i) {
d_cnjg(&z__3, &a[i + j * a_dim1]);
i__3 = i;
z__2.r = z__3.r * x[i__3].r - z__3.i * x[i__3].i,
z__2.i = z__3.r * x[i__3].i + z__3.i * x[
i__3].r;
z__1.r = temp.r - z__2.r, z__1.i = temp.i -
z__2.i;
temp.r = z__1.r, temp.i = z__1.i;
/* L100: */
}
if (nounit) {
d_cnjg(&z__2, &a[j + j * a_dim1]);
z_div(&z__1, &temp, &z__2);
temp.r = z__1.r, temp.i = z__1.i;
}
}
i__2 = j;
x[i__2].r = temp.r, x[i__2].i = temp.i;
/* L110: */
}
} else {
jx = kx;
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
ix = kx;
i__2 = jx;
temp.r = x[i__2].r, temp.i = x[i__2].i;
if (noconj) {
i__2 = j - 1;
for (i = 1; i <= i__2; ++i) {
i__3 = i + j * a_dim1;
i__4 = ix;
z__2.r = a[i__3].r * x[i__4].r - a[i__3].i * x[
i__4].i, z__2.i = a[i__3].r * x[i__4].i +
a[i__3].i * x[i__4].r;
z__1.r = temp.r - z__2.r, z__1.i = temp.i -
z__2.i;
temp.r = z__1.r, temp.i = z__1.i;
ix += *incx;
/* L120: */
}
if (nounit) {
z_div(&z__1, &temp, &a[j + j * a_dim1]);
temp.r = z__1.r, temp.i = z__1.i;
}
} else {
i__2 = j - 1;
for (i = 1; i <= i__2; ++i) {
d_cnjg(&z__3, &a[i + j * a_dim1]);
i__3 = ix;
z__2.r = z__3.r * x[i__3].r - z__3.i * x[i__3].i,
z__2.i = z__3.r * x[i__3].i + z__3.i * x[
i__3].r;
z__1.r = temp.r - z__2.r, z__1.i = temp.i -
z__2.i;
temp.r = z__1.r, temp.i = z__1.i;
ix += *incx;
/* L130: */
}
if (nounit) {
d_cnjg(&z__2, &a[j + j * a_dim1]);
z_div(&z__1, &temp, &z__2);
temp.r = z__1.r, temp.i = z__1.i;
}
}
i__2 = jx;
x[i__2].r = temp.r, x[i__2].i = temp.i;
jx += *incx;
/* L140: */
}
}
} else {
if (*incx == 1) {
for (j = *n; j >= 1; --j) {
i__1 = j;
temp.r = x[i__1].r, temp.i = x[i__1].i;
if (noconj) {
i__1 = j + 1;
for (i = *n; i >= i__1; --i) {
i__2 = i + j * a_dim1;
i__3 = i;
z__2.r = a[i__2].r * x[i__3].r - a[i__2].i * x[
i__3].i, z__2.i = a[i__2].r * x[i__3].i +
a[i__2].i * x[i__3].r;
z__1.r = temp.r - z__2.r, z__1.i = temp.i -
z__2.i;
temp.r = z__1.r, temp.i = z__1.i;
/* L150: */
}
if (nounit) {
z_div(&z__1, &temp, &a[j + j * a_dim1]);
temp.r = z__1.r, temp.i = z__1.i;
}
} else {
i__1 = j + 1;
for (i = *n; i >= i__1; --i) {
d_cnjg(&z__3, &a[i + j * a_dim1]);
i__2 = i;
z__2.r = z__3.r * x[i__2].r - z__3.i * x[i__2].i,
z__2.i = z__3.r * x[i__2].i + z__3.i * x[
i__2].r;
z__1.r = temp.r - z__2.r, z__1.i = temp.i -
z__2.i;
temp.r = z__1.r, temp.i = z__1.i;
/* L160: */
}
if (nounit) {
d_cnjg(&z__2, &a[j + j * a_dim1]);
z_div(&z__1, &temp, &z__2);
temp.r = z__1.r, temp.i = z__1.i;
}
}
i__1 = j;
x[i__1].r = temp.r, x[i__1].i = temp.i;
/* L170: */
}
} else {
kx += (*n - 1) * *incx;
jx = kx;
for (j = *n; j >= 1; --j) {
ix = kx;
i__1 = jx;
temp.r = x[i__1].r, temp.i = x[i__1].i;
if (noconj) {
i__1 = j + 1;
for (i = *n; i >= i__1; --i) {
i__2 = i + j * a_dim1;
i__3 = ix;
z__2.r = a[i__2].r * x[i__3].r - a[i__2].i * x[
i__3].i, z__2.i = a[i__2].r * x[i__3].i +
a[i__2].i * x[i__3].r;
z__1.r = temp.r - z__2.r, z__1.i = temp.i -
z__2.i;
temp.r = z__1.r, temp.i = z__1.i;
ix -= *incx;
/* L180: */
}
if (nounit) {
z_div(&z__1, &temp, &a[j + j * a_dim1]);
temp.r = z__1.r, temp.i = z__1.i;
}
} else {
i__1 = j + 1;
for (i = *n; i >= i__1; --i) {
d_cnjg(&z__3, &a[i + j * a_dim1]);
i__2 = ix;
z__2.r = z__3.r * x[i__2].r - z__3.i * x[i__2].i,
z__2.i = z__3.r * x[i__2].i + z__3.i * x[
i__2].r;
z__1.r = temp.r - z__2.r, z__1.i = temp.i -
z__2.i;
temp.r = z__1.r, temp.i = z__1.i;
ix -= *incx;
/* L190: */
}
if (nounit) {
d_cnjg(&z__2, &a[j + j * a_dim1]);
z_div(&z__1, &temp, &z__2);
temp.r = z__1.r, temp.i = z__1.i;
}
}
i__1 = jx;
x[i__1].r = temp.r, x[i__1].i = temp.i;
jx -= *incx;
/* L200: */
}
}
}
}
return 0;
/* End of ZTRSV . */
} /* ztrsv_ */
/* dznrm2.f -- translated by f2c (version of 23 April 1993 18:34:30).
You must link the resulting object file with the libraries:
-lf2c -lm (in that order)
*/
doublereal dznrm2_(n, zx, incx)
integer *n;
doublecomplex *zx;
integer *incx;
{
/* Initialized data */
static doublereal zero = 0.;
static doublereal one = 1.;
static doublereal cutlo = 8.232e-11;
static doublereal cuthi = 1.304e19;
/* Format strings */
static char fmt_30[] = "";
static char fmt_50[] = "";
static char fmt_70[] = "";
static char fmt_90[] = "";
static char fmt_110[] = "";
/* System generated locals */
integer i__1, i__2;
doublereal ret_val, d__1;
doublecomplex z__1;
/* Builtin functions */
double sqrt();
/* Local variables */
static logical imag;
static doublereal absx, xmax;
static integer next, i;
static logical scale;
static integer ix;
static doublereal hitest, sum;
/* Assigned format variables */
char *next_fmt;
/* Parameter adjustments */
--zx;
/* Function Body */
/* unitary norm of the complex n-vector stored in zx() with storage */
/* increment incx . */
/* if n .le. 0 return with result = 0. */
/* if n .ge. 1 then incx must be .ge. 1 */
/* c.l.lawson , 1978 jan 08 */
/* modified to correct problem with negative increment, 8/21/90. */
/* four phase method using two built-in constants that are */
/* hopefully applicable to all machines. */
/* cutlo = maximum of sqrt(u/eps) over all known machines. */
/* cuthi = minimum of sqrt(v) over all known machines. */
/* where */
/* eps = smallest no. such that eps + 1. .gt. 1. */
/* u = smallest positive no. (underflow limit) */
/* v = largest no. (overflow limit) */
/* brief outline of algorithm.. */
/* phase 1 scans zero components. */
/* move to phase 2 when a component is nonzero and .le. cutlo */
/* move to phase 3 when a component is .gt. cutlo */
/* move to phase 4 when a component is .ge. cuthi/m */
/* where m = n for x() real and m = 2*n for complex. */
/* values for cutlo and cuthi.. */
/* from the environmental parameters listed in the imsl converter */
/* document the limiting values are as follows.. */
/* cutlo, s.p. u/eps = 2**(-102) for honeywell. close seconds are
*/
/* univac and dec at 2**(-103) */
/* thus cutlo = 2**(-51) = 4.44089e-16 */
/* cuthi, s.p. v = 2**127 for univac, honeywell, and dec. */
/* thus cuthi = 2**(63.5) = 1.30438e19 */
/* cutlo, d.p. u/eps = 2**(-67) for honeywell and dec. */
/* thus cutlo = 2**(-33.5) = 8.23181d-11 */
/* cuthi, d.p. same as s.p. cuthi = 1.30438d19 */
/* data cutlo, cuthi / 8.232d-11, 1.304d19 / */
/* data cutlo, cuthi / 4.441e-16, 1.304e19 / */
if (*n > 0) {
goto L10;
}
ret_val = zero;
goto L300;
L10:
next = 0;
next_fmt = fmt_30;
sum = zero;
i = 1;
if (*incx < 0) {
i = (-(*n) + 1) * *incx + 1;
}
/* begin main loop */
i__1 = *n;
for (ix = 1; ix <= i__1; ++ix) {
i__2 = i;
absx = (d__1 = zx[i__2].r, abs(d__1));
imag = FALSE_;
switch ((int)next) {
case 0: goto L30;
case 1: goto L50;
case 2: goto L70;
case 3: goto L110;
case 4: goto L90;
}
L30:
if (absx > cutlo) {
goto L85;
}
next = 1;
next_fmt = fmt_50;
scale = FALSE_;
/* phase 1. sum is zero */
L50:
if (absx == zero) {
goto L200;
}
if (absx > cutlo) {
goto L85;
}
/* prepare for phase 2. */
next = 2;
next_fmt = fmt_70;
goto L105;
/* prepare for phase 4. */
L100:
next = 3;
next_fmt = fmt_110;
sum = sum / absx / absx;
L105:
scale = TRUE_;
xmax = absx;
goto L115;
/* phase 2. sum is small. */
/* scale to avoid destructive underflow.
*/
L70:
if (absx > cutlo) {
goto L75;
}
/* common code for phases 2 and 4. */
/* in phase 4 sum is large. scale to avoid overfl
ow. */
L110:
if (absx <= xmax) {
goto L115;
}
/* Computing 2nd power */
d__1 = xmax / absx;
sum = one + sum * (d__1 * d__1);
xmax = absx;
goto L200;
L115:
/* Computing 2nd power */
d__1 = absx / xmax;
sum += d__1 * d__1;
goto L200;
/* prepare for phase 3. */
L75:
sum = sum * xmax * xmax;
L85:
next = 4;
next_fmt = fmt_90;
scale = FALSE_;
/* for real or d.p. set hitest = cuthi/n */
/* for complex set hitest = cuthi/(2*n) */
hitest = cuthi / (doublereal) (*n << 1);
/* phase 3. sum is mid-range. no scaling. */
L90:
if (absx >= hitest) {
goto L100;
}
/* Computing 2nd power */
d__1 = absx;
sum += d__1 * d__1;
L200:
/* control selection of real and imaginary parts. */
if (imag) {
goto L210;
}
i__2 = i;
z__1.r = zx[i__2].r * 0. - zx[i__2].i * -1., z__1.i = zx[i__2].r *
-1. + zx[i__2].i * 0.;
absx = (d__1 = z__1.r, abs(d__1));
imag = TRUE_;
switch ((int)next) {
case 0: goto L30;
case 1: goto L50;
case 2: goto L70;
case 3: goto L110;
case 4: goto L90;
}
L210:
i += *incx;
/* L220: */
}
/* end of main loop. */
/* compute square root and adjust for scaling. */
ret_val = sqrt(sum);
if (scale) {
ret_val *= xmax;
}
L300:
return ret_val;
} /* dznrm2_ */
/* ztrmm.f -- translated by f2c (version of 23 April 1993 18:34:30).
You must link the resulting object file with the libraries:
-lf2c -lm (in that order)
*/
/* Subroutine */ int ztrmm_(side, uplo, transa, diag, m, n, alpha, a, lda, b,
ldb, side_len, uplo_len, transa_len, diag_len)
char *side, *uplo, *transa, *diag;
integer *m, *n;
doublecomplex *alpha, *a;
integer *lda;
doublecomplex *b;
integer *ldb;
ftnlen side_len;
ftnlen uplo_len;
ftnlen transa_len;
ftnlen diag_len;
{
/* System generated locals */
integer a_dim1, a_offset, b_dim1, b_offset, i__1, i__2, i__3, i__4, i__5,
i__6;
doublecomplex z__1, z__2, z__3;
/* Builtin functions */
void d_cnjg();
/* Local variables */
static integer info;
static doublecomplex temp;
static integer i, j, k;
static logical lside;
extern logical lsame_();
static integer nrowa;
static logical upper;
extern /* Subroutine */ int xerbla_();
static logical noconj, nounit;
/* .. Scalar Arguments .. */
/* .. Array Arguments .. */
/* .. */
/* Purpose */
/* ======= */
/* ZTRMM performs one of the matrix-matrix operations */
/* B := alpha*op( A )*B, or B := alpha*B*op( A ) */
/* where alpha is a scalar, B is an m by n matrix, A is a unit, or
*/
/* non-unit, upper or lower triangular matrix and op( A ) is one of
*/
/* op( A ) = A or op( A ) = A' or op( A ) = conjg( A' ). */
/* Parameters */
/* ========== */
/* SIDE - CHARACTER*1. */
/* On entry, SIDE specifies whether op( A ) multiplies B from
*/
/* the left or right as follows: */
/* SIDE = 'L' or 'l' B := alpha*op( A )*B. */
/* SIDE = 'R' or 'r' B := alpha*B*op( A ). */
/* Unchanged on exit. */
/* UPLO - CHARACTER*1. */
/* On entry, UPLO specifies whether the matrix A is an upper or
*/
/* lower triangular matrix as follows: */
/* UPLO = 'U' or 'u' A is an upper triangular matrix. */
/* UPLO = 'L' or 'l' A is a lower triangular matrix. */
/* Unchanged on exit. */
/* TRANSA - CHARACTER*1. */
/* On entry, TRANSA specifies the form of op( A ) to be used in
*/
/* the matrix multiplication as follows: */
/* TRANSA = 'N' or 'n' op( A ) = A. */
/* TRANSA = 'T' or 't' op( A ) = A'. */
/* TRANSA = 'C' or 'c' op( A ) = conjg( A' ). */
/* Unchanged on exit. */
/* DIAG - CHARACTER*1. */
/* On entry, DIAG specifies whether or not A is unit triangular
*/
/* as follows: */
/* DIAG = 'U' or 'u' A is assumed to be unit triangular. */
/* DIAG = 'N' or 'n' A is not assumed to be unit */
/* triangular. */
/* Unchanged on exit. */
/* M - INTEGER. */
/* On entry, M specifies the number of rows of B. M must be at
*/
/* least zero. */
/* Unchanged on exit. */
/* N - INTEGER. */
/* On entry, N specifies the number of columns of B. N must be
*/
/* at least zero. */
/* Unchanged on exit. */
/* ALPHA - COMPLEX*16 . */
/* On entry, ALPHA specifies the scalar alpha. When alpha is
*/
/* zero then A is not referenced and B need not be set before
*/
/* entry. */
/* Unchanged on exit. */
/* A - COMPLEX*16 array of DIMENSION ( LDA, k ), where k is m
*/
/* when SIDE = 'L' or 'l' and is n when SIDE = 'R' or 'r'.
*/
/* Before entry with UPLO = 'U' or 'u', the leading k by k
*/
/* upper triangular part of the array A must contain the upper
*/
/* triangular matrix and the strictly lower triangular part of
*/
/* A is not referenced. */
/* Before entry with UPLO = 'L' or 'l', the leading k by k
*/
/* lower triangular part of the array A must contain the lower
*/
/* triangular matrix and the strictly upper triangular part of
*/
/* A is not referenced. */
/* Note that when DIAG = 'U' or 'u', the diagonal elements of
*/
/* A are not referenced either, but are assumed to be unity.
*/
/* Unchanged on exit. */
/* LDA - INTEGER. */
/* On entry, LDA specifies the first dimension of A as declared
*/
/* in the calling (sub) program. When SIDE = 'L' or 'l' then
*/
/* LDA must be at least max( 1, m ), when SIDE = 'R' or 'r'
*/
/* then LDA must be at least max( 1, n ). */
/* Unchanged on exit. */
/* B - COMPLEX*16 array of DIMENSION ( LDB, n ). */
/* Before entry, the leading m by n part of the array B must
*/
/* contain the matrix B, and on exit is overwritten by the
*/
/* transformed matrix. */
/* LDB - INTEGER. */
/* On entry, LDB specifies the first dimension of B as declared
*/
/* in the calling (sub) program. LDB must be at least
*/
/* max( 1, m ). */
/* Unchanged on exit. */
/* 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. */
/* .. External Functions .. */
/* .. External Subroutines .. */
/* .. Intrinsic Functions .. */
/* .. Local Scalars .. */
/* .. Parameters .. */
/* .. */
/* .. Executable Statements .. */
/* Test the input parameters. */
/* Parameter adjustments */
b_dim1 = *ldb;
b_offset = b_dim1 + 1;
b -= b_offset;
a_dim1 = *lda;
a_offset = a_dim1 + 1;
a -= a_offset;
/* Function Body */
lside = lsame_(side, "L", 1L, 1L);
if (lside) {
nrowa = *m;
} else {
nrowa = *n;
}
noconj = lsame_(transa, "T", 1L, 1L);
nounit = lsame_(diag, "N", 1L, 1L);
upper = lsame_(uplo, "U", 1L, 1L);
info = 0;
if (! lside && ! lsame_(side, "R", 1L, 1L)) {
info = 1;
} else if (! upper && ! lsame_(uplo, "L", 1L, 1L)) {
info = 2;
} else if (! lsame_(transa, "N", 1L, 1L) && ! lsame_(transa, "T", 1L, 1L)
&& ! lsame_(transa, "C", 1L, 1L)) {
info = 3;
} else if (! lsame_(diag, "U", 1L, 1L) && ! lsame_(diag, "N", 1L, 1L)) {
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_("ZTRMM ", &info, 6L);
return 0;
}
/* Quick return if possible. */
if (*n == 0) {
return 0;
}
/* And when alpha.eq.zero. */
if (alpha->r == 0. && alpha->i == 0.) {
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
i__2 = *m;
for (i = 1; i <= i__2; ++i) {
i__3 = i + j * b_dim1;
b[i__3].r = 0., b[i__3].i = 0.;
/* L10: */
}
/* L20: */
}
return 0;
}
/* Start the operations. */
if (lside) {
if (lsame_(transa, "N", 1L, 1L)) {
/* 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) {
i__3 = k + j * b_dim1;
if (b[i__3].r != 0. || b[i__3].i != 0.) {
i__3 = k + j * b_dim1;
z__1.r = alpha->r * b[i__3].r - alpha->i * b[i__3]
.i, z__1.i = alpha->r * b[i__3].i +
alpha->i * b[i__3].r;
temp.r = z__1.r, temp.i = z__1.i;
i__3 = k - 1;
for (i = 1; i <= i__3; ++i) {
i__4 = i + j * b_dim1;
i__5 = i + j * b_dim1;
i__6 = i + k * a_dim1;
z__2.r = temp.r * a[i__6].r - temp.i * a[i__6]
.i, z__2.i = temp.r * a[i__6].i +
temp.i * a[i__6].r;
z__1.r = b[i__5].r + z__2.r, z__1.i = b[i__5]
.i + z__2.i;
b[i__4].r = z__1.r, b[i__4].i = z__1.i;
/* L30: */
}
if (nounit) {
i__3 = k + k * a_dim1;
z__1.r = temp.r * a[i__3].r - temp.i * a[i__3]
.i, z__1.i = temp.r * a[i__3].i +
temp.i * a[i__3].r;
temp.r = z__1.r, temp.i = z__1.i;
}
i__3 = k + j * b_dim1;
b[i__3].r = temp.r, b[i__3].i = temp.i;
}
/* L40: */
}
/* L50: */
}
} else {
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
for (k = *m; k >= 1; --k) {
i__2 = k + j * b_dim1;
if (b[i__2].r != 0. || b[i__2].i != 0.) {
i__2 = k + j * b_dim1;
z__1.r = alpha->r * b[i__2].r - alpha->i * b[i__2]
.i, z__1.i = alpha->r * b[i__2].i +
alpha->i * b[i__2].r;
temp.r = z__1.r, temp.i = z__1.i;
i__2 = k + j * b_dim1;
b[i__2].r = temp.r, b[i__2].i = temp.i;
if (nounit) {
i__2 = k + j * b_dim1;
i__3 = k + j * b_dim1;
i__4 = k + k * a_dim1;
z__1.r = b[i__3].r * a[i__4].r - b[i__3].i *
a[i__4].i, z__1.i = b[i__3].r * a[
i__4].i + b[i__3].i * a[i__4].r;
b[i__2].r = z__1.r, b[i__2].i = z__1.i;
}
i__2 = *m;
for (i = k + 1; i <= i__2; ++i) {
i__3 = i + j * b_dim1;
i__4 = i + j * b_dim1;
i__5 = i + k * a_dim1;
z__2.r = temp.r * a[i__5].r - temp.i * a[i__5]
.i, z__2.i = temp.r * a[i__5].i +
temp.i * a[i__5].r;
z__1.r = b[i__4].r + z__2.r, z__1.i = b[i__4]
.i + z__2.i;
b[i__3].r = z__1.r, b[i__3].i = z__1.i;
/* L60: */
}
}
/* L70: */
}
/* L80: */
}
}
} else {
/* Form B := alpha*B*A' or B := alpha*B*conjg( A' )
. */
if (upper) {
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
for (i = *m; i >= 1; --i) {
i__2 = i + j * b_dim1;
temp.r = b[i__2].r, temp.i = b[i__2].i;
if (noconj) {
if (nounit) {
i__2 = i + i * a_dim1;
z__1.r = temp.r * a[i__2].r - temp.i * a[i__2]
.i, z__1.i = temp.r * a[i__2].i +
temp.i * a[i__2].r;
temp.r = z__1.r, temp.i = z__1.i;
}
i__2 = i - 1;
for (k = 1; k <= i__2; ++k) {
i__3 = k + i * a_dim1;
i__4 = k + j * b_dim1;
z__2.r = a[i__3].r * b[i__4].r - a[i__3].i *
b[i__4].i, z__2.i = a[i__3].r * b[
i__4].i + a[i__3].i * b[i__4].r;
z__1.r = temp.r + z__2.r, z__1.i = temp.i +
z__2.i;
temp.r = z__1.r, temp.i = z__1.i;
/* L90: */
}
} else {
if (nounit) {
d_cnjg(&z__2, &a[i + i * a_dim1]);
z__1.r = temp.r * z__2.r - temp.i * z__2.i,
z__1.i = temp.r * z__2.i + temp.i *
z__2.r;
temp.r = z__1.r, temp.i = z__1.i;
}
i__2 = i - 1;
for (k = 1; k <= i__2; ++k) {
d_cnjg(&z__3, &a[k + i * a_dim1]);
i__3 = k + j * b_dim1;
z__2.r = z__3.r * b[i__3].r - z__3.i * b[i__3]
.i, z__2.i = z__3.r * b[i__3].i +
z__3.i * b[i__3].r;
z__1.r = temp.r + z__2.r, z__1.i = temp.i +
z__2.i;
temp.r = z__1.r, temp.i = z__1.i;
/* L100: */
}
}
i__2 = i + j * b_dim1;
z__1.r = alpha->r * temp.r - alpha->i * temp.i,
z__1.i = alpha->r * temp.i + alpha->i *
temp.r;
b[i__2].r = z__1.r, b[i__2].i = z__1.i;
/* L110: */
}
/* L120: */
}
} else {
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
i__2 = *m;
for (i = 1; i <= i__2; ++i) {
i__3 = i + j * b_dim1;
temp.r = b[i__3].r, temp.i = b[i__3].i;
if (noconj) {
if (nounit) {
i__3 = i + i * a_dim1;
z__1.r = temp.r * a[i__3].r - temp.i * a[i__3]
.i, z__1.i = temp.r * a[i__3].i +
temp.i * a[i__3].r;
temp.r = z__1.r, temp.i = z__1.i;
}
i__3 = *m;
for (k = i + 1; k <= i__3; ++k) {
i__4 = k + i * a_dim1;
i__5 = k + j * b_dim1;
z__2.r = a[i__4].r * b[i__5].r - a[i__4].i *
b[i__5].i, z__2.i = a[i__4].r * b[
i__5].i + a[i__4].i * b[i__5].r;
z__1.r = temp.r + z__2.r, z__1.i = temp.i +
z__2.i;
temp.r = z__1.r, temp.i = z__1.i;
/* L130: */
}
} else {
if (nounit) {
d_cnjg(&z__2, &a[i + i * a_dim1]);
z__1.r = temp.r * z__2.r - temp.i * z__2.i,
z__1.i = temp.r * z__2.i + temp.i *
z__2.r;
temp.r = z__1.r, temp.i = z__1.i;
}
i__3 = *m;
for (k = i + 1; k <= i__3; ++k) {
d_cnjg(&z__3, &a[k + i * a_dim1]);
i__4 = k + j * b_dim1;
z__2.r = z__3.r * b[i__4].r - z__3.i * b[i__4]
.i, z__2.i = z__3.r * b[i__4].i +
z__3.i * b[i__4].r;
z__1.r = temp.r + z__2.r, z__1.i = temp.i +
z__2.i;
temp.r = z__1.r, temp.i = z__1.i;
/* L140: */
}
}
i__3 = i + j * b_dim1;
z__1.r = alpha->r * temp.r - alpha->i * temp.i,
z__1.i = alpha->r * temp.i + alpha->i *
temp.r;
b[i__3].r = z__1.r, b[i__3].i = z__1.i;
/* L150: */
}
/* L160: */
}
}
}
} else {
if (lsame_(transa, "N", 1L, 1L)) {
/* Form B := alpha*B*A. */
if (upper) {
for (j = *n; j >= 1; --j) {
temp.r = alpha->r, temp.i = alpha->i;
if (nounit) {
i__1 = j + j * a_dim1;
z__1.r = temp.r * a[i__1].r - temp.i * a[i__1].i,
z__1.i = temp.r * a[i__1].i + temp.i * a[i__1]
.r;
temp.r = z__1.r, temp.i = z__1.i;
}
i__1 = *m;
for (i = 1; i <= i__1; ++i) {
i__2 = i + j * b_dim1;
i__3 = i + j * b_dim1;
z__1.r = temp.r * b[i__3].r - temp.i * b[i__3].i,
z__1.i = temp.r * b[i__3].i + temp.i * b[i__3]
.r;
b[i__2].r = z__1.r, b[i__2].i = z__1.i;
/* L170: */
}
i__1 = j - 1;
for (k = 1; k <= i__1; ++k) {
i__2 = k + j * a_dim1;
if (a[i__2].r != 0. || a[i__2].i != 0.) {
i__2 = k + j * a_dim1;
z__1.r = alpha->r * a[i__2].r - alpha->i * a[i__2]
.i, z__1.i = alpha->r * a[i__2].i +
alpha->i * a[i__2].r;
temp.r = z__1.r, temp.i = z__1.i;
i__2 = *m;
for (i = 1; i <= i__2; ++i) {
i__3 = i + j * b_dim1;
i__4 = i + j * b_dim1;
i__5 = i + k * b_dim1;
z__2.r = temp.r * b[i__5].r - temp.i * b[i__5]
.i, z__2.i = temp.r * b[i__5].i +
temp.i * b[i__5].r;
z__1.r = b[i__4].r + z__2.r, z__1.i = b[i__4]
.i + z__2.i;
b[i__3].r = z__1.r, b[i__3].i = z__1.i;
/* L180: */
}
}
/* L190: */
}
/* L200: */
}
} else {
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
temp.r = alpha->r, temp.i = alpha->i;
if (nounit) {
i__2 = j + j * a_dim1;
z__1.r = temp.r * a[i__2].r - temp.i * a[i__2].i,
z__1.i = temp.r * a[i__2].i + temp.i * a[i__2]
.r;
temp.r = z__1.r, temp.i = z__1.i;
}
i__2 = *m;
for (i = 1; i <= i__2; ++i) {
i__3 = i + j * b_dim1;
i__4 = i + j * b_dim1;
z__1.r = temp.r * b[i__4].r - temp.i * b[i__4].i,
z__1.i = temp.r * b[i__4].i + temp.i * b[i__4]
.r;
b[i__3].r = z__1.r, b[i__3].i = z__1.i;
/* L210: */
}
i__2 = *n;
for (k = j + 1; k <= i__2; ++k) {
i__3 = k + j * a_dim1;
if (a[i__3].r != 0. || a[i__3].i != 0.) {
i__3 = k + j * a_dim1;
z__1.r = alpha->r * a[i__3].r - alpha->i * a[i__3]
.i, z__1.i = alpha->r * a[i__3].i +
alpha->i * a[i__3].r;
temp.r = z__1.r, temp.i = z__1.i;
i__3 = *m;
for (i = 1; i <= i__3; ++i) {
i__4 = i + j * b_dim1;
i__5 = i + j * b_dim1;
i__6 = i + k * b_dim1;
z__2.r = temp.r * b[i__6].r - temp.i * b[i__6]
.i, z__2.i = temp.r * b[i__6].i +
temp.i * b[i__6].r;
z__1.r = b[i__5].r + z__2.r, z__1.i = b[i__5]
.i + z__2.i;
b[i__4].r = z__1.r, b[i__4].i = z__1.i;
/* L220: */
}
}
/* L230: */
}
/* L240: */
}
}
} else {
/* Form B := alpha*B*A' or B := alpha*B*conjg( A' )
. */
if (upper) {
i__1 = *n;
for (k = 1; k <= i__1; ++k) {
i__2 = k - 1;
for (j = 1; j <= i__2; ++j) {
i__3 = j + k * a_dim1;
if (a[i__3].r != 0. || a[i__3].i != 0.) {
if (noconj) {
i__3 = j + k * a_dim1;
z__1.r = alpha->r * a[i__3].r - alpha->i * a[
i__3].i, z__1.i = alpha->r * a[i__3]
.i + alpha->i * a[i__3].r;
temp.r = z__1.r, temp.i = z__1.i;
} else {
d_cnjg(&z__2, &a[j + k * a_dim1]);
z__1.r = alpha->r * z__2.r - alpha->i *
z__2.i, z__1.i = alpha->r * z__2.i +
alpha->i * z__2.r;
temp.r = z__1.r, temp.i = z__1.i;
}
i__3 = *m;
for (i = 1; i <= i__3; ++i) {
i__4 = i + j * b_dim1;
i__5 = i + j * b_dim1;
i__6 = i + k * b_dim1;
z__2.r = temp.r * b[i__6].r - temp.i * b[i__6]
.i, z__2.i = temp.r * b[i__6].i +
temp.i * b[i__6].r;
z__1.r = b[i__5].r + z__2.r, z__1.i = b[i__5]
.i + z__2.i;
b[i__4].r = z__1.r, b[i__4].i = z__1.i;
/* L250: */
}
}
/* L260: */
}
temp.r = alpha->r, temp.i = alpha->i;
if (nounit) {
if (noconj) {
i__2 = k + k * a_dim1;
z__1.r = temp.r * a[i__2].r - temp.i * a[i__2].i,
z__1.i = temp.r * a[i__2].i + temp.i * a[
i__2].r;
temp.r = z__1.r, temp.i = z__1.i;
} else {
d_cnjg(&z__2, &a[k + k * a_dim1]);
z__1.r = temp.r * z__2.r - temp.i * z__2.i,
z__1.i = temp.r * z__2.i + temp.i *
z__2.r;
temp.r = z__1.r, temp.i = z__1.i;
}
}
if (temp.r != 1. || temp.i != 0.) {
i__2 = *m;
for (i = 1; i <= i__2; ++i) {
i__3 = i + k * b_dim1;
i__4 = i + k * b_dim1;
z__1.r = temp.r * b[i__4].r - temp.i * b[i__4].i,
z__1.i = temp.r * b[i__4].i + temp.i * b[
i__4].r;
b[i__3].r = z__1.r, b[i__3].i = z__1.i;
/* L270: */
}
}
/* L280: */
}
} else {
for (k = *n; k >= 1; --k) {
i__1 = *n;
for (j = k + 1; j <= i__1; ++j) {
i__2 = j + k * a_dim1;
if (a[i__2].r != 0. || a[i__2].i != 0.) {
if (noconj) {
i__2 = j + k * a_dim1;
z__1.r = alpha->r * a[i__2].r - alpha->i * a[
i__2].i, z__1.i = alpha->r * a[i__2]
.i + alpha->i * a[i__2].r;
temp.r = z__1.r, temp.i = z__1.i;
} else {
d_cnjg(&z__2, &a[j + k * a_dim1]);
z__1.r = alpha->r * z__2.r - alpha->i *
z__2.i, z__1.i = alpha->r * z__2.i +
alpha->i * z__2.r;
temp.r = z__1.r, temp.i = z__1.i;
}
i__2 = *m;
for (i = 1; i <= i__2; ++i) {
i__3 = i + j * b_dim1;
i__4 = i + j * b_dim1;
i__5 = i + k * b_dim1;
z__2.r = temp.r * b[i__5].r - temp.i * b[i__5]
.i, z__2.i = temp.r * b[i__5].i +
temp.i * b[i__5].r;
z__1.r = b[i__4].r + z__2.r, z__1.i = b[i__4]
.i + z__2.i;
b[i__3].r = z__1.r, b[i__3].i = z__1.i;
/* L290: */
}
}
/* L300: */
}
temp.r = alpha->r, temp.i = alpha->i;
if (nounit) {
if (noconj) {
i__1 = k + k * a_dim1;
z__1.r = temp.r * a[i__1].r - temp.i * a[i__1].i,
z__1.i = temp.r * a[i__1].i + temp.i * a[
i__1].r;
temp.r = z__1.r, temp.i = z__1.i;
} else {
d_cnjg(&z__2, &a[k + k * a_dim1]);
z__1.r = temp.r * z__2.r - temp.i * z__2.i,
z__1.i = temp.r * z__2.i + temp.i *
z__2.r;
temp.r = z__1.r, temp.i = z__1.i;
}
}
if (temp.r != 1. || temp.i != 0.) {
i__1 = *m;
for (i = 1; i <= i__1; ++i) {
i__2 = i + k * b_dim1;
i__3 = i + k * b_dim1;
z__1.r = temp.r * b[i__3].r - temp.i * b[i__3].i,
z__1.i = temp.r * b[i__3].i + temp.i * b[
i__3].r;
b[i__2].r = z__1.r, b[i__2].i = z__1.i;
/* L310: */
}
}
/* L320: */
}
}
}
}
return 0;
/* End of ZTRMM . */
} /* ztrmm_ */
/* daxpy.f -- translated by f2c (version of 23 April 1993 18:34:30).
You must link the resulting object file with the libraries:
-lf2c -lm (in that order)
*/
/* Subroutine */ int daxpy_(n, da, dx, incx, dy, incy)
integer *n;
doublereal *da, *dx;
integer *incx;
doublereal *dy;
integer *incy;
{
/* System generated locals */
integer i__1;
/* Local variables */
static integer i, m, ix, iy, mp1;
/* constant times a vector plus a vector. */
/* uses unrolled loops for increments equal to one. */
/* jack dongarra, linpack, 3/11/78. */
/* Parameter adjustments */
--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;
} /* daxpy_ */
/* dnrm2.f -- translated by f2c (version of 23 April 1993 18:34:30).
You must link the resulting object file with the libraries:
-lf2c -lm (in that order)
*/
doublereal dnrm2_(n, dx, incx)
integer *n;
doublereal *dx;
integer *incx;
{
/* Initialized data */
static doublereal zero = 0.;
static doublereal one = 1.;
static doublereal cutlo = 8.232e-11;
static doublereal cuthi = 1.304e19;
/* Format strings */
static char fmt_30[] = "";
static char fmt_50[] = "";
static char fmt_70[] = "";
static char fmt_110[] = "";
/* System generated locals */
integer i__1;
doublereal ret_val, d__1;
/* Builtin functions */
double sqrt();
/* Local variables */
static doublereal xmax;
static integer next, i, j, ix;
static doublereal hitest, sum;
/* Assigned format variables */
char *next_fmt;
/* Parameter adjustments */
--dx;
/* Function Body */
/* euclidean norm of the n-vector stored in dx() with storage */
/* increment incx . */
/* if n .le. 0 return with result = 0. */
/* if n .ge. 1 then incx must be .ge. 1 */
/* c.l.lawson, 1978 jan 08 */
/* modified to correct problem with negative increment, 8/21/90. */
/* modified to correct failure to update ix, 1/25/92. */
/* four phase method using two built-in constants that are */
/* hopefully applicable to all machines. */
/* cutlo = maximum of dsqrt(u/eps) over all known machines. */
/* cuthi = minimum of dsqrt(v) over all known machines. */
/* where */
/* eps = smallest no. such that eps + 1. .gt. 1. */
/* u = smallest positive no. (underflow limit) */
/* v = largest no. (overflow limit) */
/* brief outline of algorithm.. */
/* phase 1 scans zero components. */
/* move to phase 2 when a component is nonzero and .le. cutlo */
/* move to phase 3 when a component is .gt. cutlo */
/* move to phase 4 when a component is .ge. cuthi/m */
/* where m = n for x() real and m = 2*n for complex. */
/* values for cutlo and cuthi.. */
/* from the environmental parameters listed in the imsl converter */
/* document the limiting values are as follows.. */
/* cutlo, s.p. u/eps = 2**(-102) for honeywell. close seconds are
*/
/* univac and dec at 2**(-103) */
/* thus cutlo = 2**(-51) = 4.44089e-16 */
/* cuthi, s.p. v = 2**127 for univac, honeywell, and dec. */
/* thus cuthi = 2**(63.5) = 1.30438e19 */
/* cutlo, d.p. u/eps = 2**(-67) for honeywell and dec. */
/* thus cutlo = 2**(-33.5) = 8.23181d-11 */
/* cuthi, d.p. same as s.p. cuthi = 1.30438d19 */
/* data cutlo, cuthi / 8.232d-11, 1.304d19 / */
/* data cutlo, cuthi / 4.441e-16, 1.304e19 / */
if (*n > 0) {
goto L10;
}
ret_val = zero;
goto L300;
L10:
next = 0;
next_fmt = fmt_30;
sum = zero;
i = 1;
if (*incx < 0) {
i = (-(*n) + 1) * *incx + 1;
}
ix = 1;
/* begin main loop */
L20:
switch ((int)next) {
case 0: goto L30;
case 1: goto L50;
case 2: goto L70;
case 3: goto L110;
}
L30:
if ((d__1 = dx[i], abs(d__1)) > cutlo) {
goto L85;
}
next = 1;
next_fmt = fmt_50;
xmax = zero;
/* phase 1. sum is zero */
L50:
if (dx[i] == zero) {
goto L200;
}
if ((d__1 = dx[i], abs(d__1)) > cutlo) {
goto L85;
}
/* prepare for phase 2. */
next = 2;
next_fmt = fmt_70;
goto L105;
/* prepare for phase 4. */
L100:
ix = j;
next = 3;
next_fmt = fmt_110;
sum = sum / dx[i] / dx[i];
L105:
xmax = (d__1 = dx[i], abs(d__1));
goto L115;
/* phase 2. sum is small. */
/* scale to avoid destructive underflow. */
L70:
if ((d__1 = dx[i], abs(d__1)) > cutlo) {
goto L75;
}
/* common code for phases 2 and 4. */
/* in phase 4 sum is large. scale to avoid overflow.
*/
L110:
if ((d__1 = dx[i], abs(d__1)) <= xmax) {
goto L115;
}
/* Computing 2nd power */
d__1 = xmax / dx[i];
sum = one + sum * (d__1 * d__1);
xmax = (d__1 = dx[i], abs(d__1));
goto L200;
L115:
/* Computing 2nd power */
d__1 = dx[i] / xmax;
sum += d__1 * d__1;
goto L200;
/* prepare for phase 3. */
L75:
sum = sum * xmax * xmax;
/* for real or d.p. set hitest = cuthi/n */
/* for complex set hitest = cuthi/(2*n) */
L85:
hitest = cuthi / (real) (*n);
/* phase 3. sum is mid-range. no scaling. */
i__1 = *n;
for (j = ix; j <= i__1; ++j) {
if ((d__1 = dx[i], abs(d__1)) >= hitest) {
goto L100;
}
/* Computing 2nd power */
d__1 = dx[i];
sum += d__1 * d__1;
i += *incx;
/* L95: */
}
ret_val = sqrt(sum);
goto L300;
L200:
++ix;
i += *incx;
if (ix <= *n) {
goto L20;
}
/* end of main loop. */
/* compute square root and adjust for scaling. */
ret_val = xmax * sqrt(sum);
L300:
return ret_val;
} /* dnrm2_ */
/* zdotc.f -- translated by f2c (version of 23 April 1993 18:34:30).
You must link the resulting object file with the libraries:
-lf2c -lm (in that order)
*/
/* Double Complex */ int zdotc_( ret_val, n, zx, incx, zy, incy)
doublecomplex * ret_val;
integer *n;
doublecomplex *zx;
integer *incx;
doublecomplex *zy;
integer *incy;
{
/* System generated locals */
integer i__1, i__2;
doublecomplex z__1, z__2, z__3;
/* Builtin functions */
void d_cnjg();
/* Local variables */
static integer i;
static doublecomplex ztemp;
static integer ix, iy;
/* forms the dot product of a vector. */
/* jack dongarra, 3/11/78. */
/* Parameter adjustments */
--zy;
--zx;
/* Function Body */
ztemp.r = 0., ztemp.i = 0.;
ret_val->r = 0., ret_val->i = 0.;
if (*n <= 0) {
return ;
}
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) {
d_cnjg(&z__3, &zx[ix]);
i__2 = iy;
z__2.r = z__3.r * zy[i__2].r - z__3.i * zy[i__2].i, z__2.i = z__3.r *
zy[i__2].i + z__3.i * zy[i__2].r;
z__1.r = ztemp.r + z__2.r, z__1.i = ztemp.i + z__2.i;
ztemp.r = z__1.r, ztemp.i = z__1.i;
ix += *incx;
iy += *incy;
/* L10: */
}
ret_val->r = ztemp.r, ret_val->i = ztemp.i;
return ;
/* code for both increments equal to 1 */
L20:
i__1 = *n;
for (i = 1; i <= i__1; ++i) {
d_cnjg(&z__3, &zx[i]);
i__2 = i;
z__2.r = z__3.r * zy[i__2].r - z__3.i * zy[i__2].i, z__2.i = z__3.r *
zy[i__2].i + z__3.i * zy[i__2].r;
z__1.r = ztemp.r + z__2.r, z__1.i = ztemp.i + z__2.i;
ztemp.r = z__1.r, ztemp.i = z__1.i;
/* L30: */
}
ret_val->r = ztemp.r, ret_val->i = ztemp.i;
return ;
} /* zdotc_ */
/* zgemv.f -- translated by f2c (version of 23 April 1993 18:34:30).
You must link the resulting object file with the libraries:
-lf2c -lm (in that order)
*/
/* Subroutine */ int zgemv_(trans, m, n, alpha, a, lda, x, incx, beta, y,
incy, trans_len)
char *trans;
integer *m, *n;
doublecomplex *alpha, *a;
integer *lda;
doublecomplex *x;
integer *incx;
doublecomplex *beta, *y;
integer *incy;
ftnlen trans_len;
{
/* System generated locals */
integer a_dim1, a_offset, i__1, i__2, i__3, i__4, i__5;
doublecomplex z__1, z__2, z__3;
/* Builtin functions */
void d_cnjg();
/* Local variables */
static integer info;
static doublecomplex temp;
static integer lenx, leny, i, j;
extern logical lsame_();
static integer ix, iy, jx, jy, kx, ky;
extern /* Subroutine */ int xerbla_();
static logical noconj;
/* .. Scalar Arguments .. */
/* .. Array Arguments .. */
/* .. */
/* Purpose */
/* ======= */
/* ZGEMV performs one of the matrix-vector operations */
/* y := alpha*A*x + beta*y, or y := alpha*A'*x + beta*y, or */
/* y := alpha*conjg( A' )*x + beta*y, */
/* where alpha and beta are scalars, x and y are vectors and A is an */
/* m by n matrix. */
/* Parameters */
/* ========== */
/* TRANS - CHARACTER*1. */
/* On entry, TRANS specifies the operation to be performed as */
/* follows: */
/* TRANS = 'N' or 'n' y := alpha*A*x + beta*y. */
/* TRANS = 'T' or 't' y := alpha*A'*x + beta*y. */
/* TRANS = 'C' or 'c' y := alpha*conjg( A' )*x + beta*y. */
/* Unchanged on exit. */
/* M - INTEGER. */
/* On entry, M specifies the number of rows of the matrix A. */
/* M must be at least zero. */
/* Unchanged on exit. */
/* N - INTEGER. */
/* On entry, N specifies the number of columns of the matrix A.
*/
/* N must be at least zero. */
/* Unchanged on exit. */
/* ALPHA - COMPLEX*16 . */
/* On entry, ALPHA specifies the scalar alpha. */
/* Unchanged on exit. */
/* A - COMPLEX*16 array of DIMENSION ( LDA, n ). */
/* Before entry, the leading m by n part of the array A must */
/* contain the matrix of coefficients. */
/* Unchanged on exit. */
/* LDA - INTEGER. */
/* On entry, LDA specifies the first dimension of A as declared
*/
/* in the calling (sub) program. LDA must be at least */
/* max( 1, m ). */
/* Unchanged on exit. */
/* X - COMPLEX*16 array of DIMENSION at least */
/* ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n' */
/* and at least */
/* ( 1 + ( m - 1 )*abs( INCX ) ) otherwise. */
/* Before entry, the incremented array X must contain the */
/* vector x. */
/* Unchanged on exit. */
/* INCX - INTEGER. */
/* On entry, INCX specifies the increment for the elements of */
/* X. INCX must not be zero. */
/* Unchanged on exit. */
/* BETA - COMPLEX*16 . */
/* On entry, BETA specifies the scalar beta. When BETA is */
/* supplied as zero then Y need not be set on input. */
/* Unchanged on exit. */
/* Y - COMPLEX*16 array of DIMENSION at least */
/* ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n' */
/* and at least */
/* ( 1 + ( n - 1 )*abs( INCY ) ) otherwise. */
/* Before entry with BETA non-zero, the incremented array Y */
/* must contain the vector y. On exit, Y is overwritten by the
*/
/* updated vector y. */
/* INCY - INTEGER. */
/* On entry, INCY specifies the increment for the elements of */
/* Y. INCY must not be zero. */
/* Unchanged on exit. */
/* Level 2 Blas routine. */
/* -- 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. */
/* .. Parameters .. */
/* .. Local Scalars .. */
/* .. External Functions .. */
/* .. External Subroutines .. */
/* .. Intrinsic Functions .. */
/* .. */
/* .. Executable Statements .. */
/* Test the input parameters. */
/* Parameter adjustments */
--y;
--x;
a_dim1 = *lda;
a_offset = a_dim1 + 1;
a -= a_offset;
/* Function Body */
info = 0;
if (! lsame_(trans, "N", 1L, 1L) && ! lsame_(trans, "T", 1L, 1L) && !
lsame_(trans, "C", 1L, 1L)) {
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_("ZGEMV ", &info, 6L);
return 0;
}
/* Quick return if possible. */
if (*m == 0 || *n == 0 || alpha->r == 0. && alpha->i == 0. && (beta->r ==
1. && beta->i == 0.)) {
return 0;
}
noconj = lsame_(trans, "T", 1L, 1L);
/* 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", 1L, 1L)) {
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->r != 1. || beta->i != 0.) {
if (*incy == 1) {
if (beta->r == 0. && beta->i == 0.) {
i__1 = leny;
for (i = 1; i <= i__1; ++i) {
i__2 = i;
y[i__2].r = 0., y[i__2].i = 0.;
/* L10: */
}
} else {
i__1 = leny;
for (i = 1; i <= i__1; ++i) {
i__2 = i;
i__3 = i;
z__1.r = beta->r * y[i__3].r - beta->i * y[i__3].i,
z__1.i = beta->r * y[i__3].i + beta->i * y[i__3]
.r;
y[i__2].r = z__1.r, y[i__2].i = z__1.i;
/* L20: */
}
}
} else {
iy = ky;
if (beta->r == 0. && beta->i == 0.) {
i__1 = leny;
for (i = 1; i <= i__1; ++i) {
i__2 = iy;
y[i__2].r = 0., y[i__2].i = 0.;
iy += *incy;
/* L30: */
}
} else {
i__1 = leny;
for (i = 1; i <= i__1; ++i) {
i__2 = iy;
i__3 = iy;
z__1.r = beta->r * y[i__3].r - beta->i * y[i__3].i,
z__1.i = beta->r * y[i__3].i + beta->i * y[i__3]
.r;
y[i__2].r = z__1.r, y[i__2].i = z__1.i;
iy += *incy;
/* L40: */
}
}
}
}
if (alpha->r == 0. && alpha->i == 0.) {
return 0;
}
if (lsame_(trans, "N", 1L, 1L)) {
/* Form y := alpha*A*x + y. */
jx = kx;
if (*incy == 1) {
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
i__2 = jx;
if (x[i__2].r != 0. || x[i__2].i != 0.) {
i__2 = jx;
z__1.r = alpha->r * x[i__2].r - alpha->i * x[i__2].i,
z__1.i = alpha->r * x[i__2].i + alpha->i * x[i__2]
.r;
temp.r = z__1.r, temp.i = z__1.i;
i__2 = *m;
for (i = 1; i <= i__2; ++i) {
i__3 = i;
i__4 = i;
i__5 = i + j * a_dim1;
z__2.r = temp.r * a[i__5].r - temp.i * a[i__5].i,
z__2.i = temp.r * a[i__5].i + temp.i * a[i__5]
.r;
z__1.r = y[i__4].r + z__2.r, z__1.i = y[i__4].i +
z__2.i;
y[i__3].r = z__1.r, y[i__3].i = z__1.i;
/* L50: */
}
}
jx += *incx;
/* L60: */
}
} else {
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
i__2 = jx;
if (x[i__2].r != 0. || x[i__2].i != 0.) {
i__2 = jx;
z__1.r = alpha->r * x[i__2].r - alpha->i * x[i__2].i,
z__1.i = alpha->r * x[i__2].i + alpha->i * x[i__2]
.r;
temp.r = z__1.r, temp.i = z__1.i;
iy = ky;
i__2 = *m;
for (i = 1; i <= i__2; ++i) {
i__3 = iy;
i__4 = iy;
i__5 = i + j * a_dim1;
z__2.r = temp.r * a[i__5].r - temp.i * a[i__5].i,
z__2.i = temp.r * a[i__5].i + temp.i * a[i__5]
.r;
z__1.r = y[i__4].r + z__2.r, z__1.i = y[i__4].i +
z__2.i;
y[i__3].r = z__1.r, y[i__3].i = z__1.i;
iy += *incy;
/* L70: */
}
}
jx += *incx;
/* L80: */
}
}
} else {
/* Form y := alpha*A'*x + y or y := alpha*conjg( A' )*x + y.
*/
jy = ky;
if (*incx == 1) {
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
temp.r = 0., temp.i = 0.;
if (noconj) {
i__2 = *m;
for (i = 1; i <= i__2; ++i) {
i__3 = i + j * a_dim1;
i__4 = i;
z__2.r = a[i__3].r * x[i__4].r - a[i__3].i * x[i__4]
.i, z__2.i = a[i__3].r * x[i__4].i + a[i__3]
.i * x[i__4].r;
z__1.r = temp.r + z__2.r, z__1.i = temp.i + z__2.i;
temp.r = z__1.r, temp.i = z__1.i;
/* L90: */
}
} else {
i__2 = *m;
for (i = 1; i <= i__2; ++i) {
d_cnjg(&z__3, &a[i + j * a_dim1]);
i__3 = i;
z__2.r = z__3.r * x[i__3].r - z__3.i * x[i__3].i,
z__2.i = z__3.r * x[i__3].i + z__3.i * x[i__3]
.r;
z__1.r = temp.r + z__2.r, z__1.i = temp.i + z__2.i;
temp.r = z__1.r, temp.i = z__1.i;
/* L100: */
}
}
i__2 = jy;
i__3 = jy;
z__2.r = alpha->r * temp.r - alpha->i * temp.i, z__2.i =
alpha->r * temp.i + alpha->i * temp.r;
z__1.r = y[i__3].r + z__2.r, z__1.i = y[i__3].i + z__2.i;
y[i__2].r = z__1.r, y[i__2].i = z__1.i;
jy += *incy;
/* L110: */
}
} else {
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
temp.r = 0., temp.i = 0.;
ix = kx;
if (noconj) {
i__2 = *m;
for (i = 1; i <= i__2; ++i) {
i__3 = i + j * a_dim1;
i__4 = ix;
z__2.r = a[i__3].r * x[i__4].r - a[i__3].i * x[i__4]
.i, z__2.i = a[i__3].r * x[i__4].i + a[i__3]
.i * x[i__4].r;
z__1.r = temp.r + z__2.r, z__1.i = temp.i + z__2.i;
temp.r = z__1.r, temp.i = z__1.i;
ix += *incx;
/* L120: */
}
} else {
i__2 = *m;
for (i = 1; i <= i__2; ++i) {
d_cnjg(&z__3, &a[i + j * a_dim1]);
i__3 = ix;
z__2.r = z__3.r * x[i__3].r - z__3.i * x[i__3].i,
z__2.i = z__3.r * x[i__3].i + z__3.i * x[i__3]
.r;
z__1.r = temp.r + z__2.r, z__1.i = temp.i + z__2.i;
temp.r = z__1.r, temp.i = z__1.i;
ix += *incx;
/* L130: */
}
}
i__2 = jy;
i__3 = jy;
z__2.r = alpha->r * temp.r - alpha->i * temp.i, z__2.i =
alpha->r * temp.i + alpha->i * temp.r;
z__1.r = y[i__3].r + z__2.r, z__1.i = y[i__3].i + z__2.i;
y[i__2].r = z__1.r, y[i__2].i = z__1.i;
jy += *incy;
/* L140: */
}
}
}
return 0;
/* End of ZGEMV . */
} /* zgemv_ */
/* zgeru.f -- translated by f2c (version of 23 April 1993 18:34:30).
You must link the resulting object file with the libraries:
-lf2c -lm (in that order)
*/
/* Subroutine */ int zgeru_(m, n, alpha, x, incx, y, incy, a, lda)
integer *m, *n;
doublecomplex *alpha, *x;
integer *incx;
doublecomplex *y;
integer *incy;
doublecomplex *a;
integer *lda;
{
/* System generated locals */
integer a_dim1, a_offset, i__1, i__2, i__3, i__4, i__5;
doublecomplex z__1, z__2;
/* Local variables */
static integer info;
static doublecomplex temp;
static integer i, j, ix, jy, kx;
extern /* Subroutine */ int xerbla_();
/* .. Scalar Arguments .. */
/* .. Array Arguments .. */
/* .. */
/* Purpose */
/* ======= */
/* ZGERU performs the rank 1 operation */
/* A := alpha*x*y' + A, */
/* where alpha is a scalar, x is an m element vector, y is an n element
*/
/* vector and A is an m by n matrix. */
/* Parameters */
/* ========== */
/* M - INTEGER. */
/* On entry, M specifies the number of rows of the matrix A. */
/* M must be at least zero. */
/* Unchanged on exit. */
/* N - INTEGER. */
/* On entry, N specifies the number of columns of the matrix A.
*/
/* N must be at least zero. */
/* Unchanged on exit. */
/* ALPHA - COMPLEX*16 . */
/* On entry, ALPHA specifies the scalar alpha. */
/* Unchanged on exit. */
/* X - COMPLEX*16 array of dimension at least */
/* ( 1 + ( m - 1 )*abs( INCX ) ). */
/* Before entry, the incremented array X must contain the m */
/* element vector x. */
/* Unchanged on exit. */
/* INCX - INTEGER. */
/* On entry, INCX specifies the increment for the elements of */
/* X. INCX must not be zero. */
/* Unchanged on exit. */
/* Y - COMPLEX*16 array of dimension at least */
/* ( 1 + ( n - 1 )*abs( INCY ) ). */
/* Before entry, the incremented array Y must contain the n */
/* element vector y. */
/* Unchanged on exit. */
/* INCY - INTEGER. */
/* On entry, INCY specifies the increment for the elements of */
/* Y. INCY must not be zero. */
/* Unchanged on exit. */
/* A - COMPLEX*16 array of DIMENSION ( LDA, n ). */
/* Before entry, the leading m by n part of the array A must */
/* contain the matrix of coefficients. On exit, A is */
/* overwritten by the updated matrix. */
/* LDA - INTEGER. */
/* On entry, LDA specifies the first dimension of A as declared
*/
/* in the calling (sub) program. LDA must be at least */
/* max( 1, m ). */
/* Unchanged on exit. */
/* Level 2 Blas routine. */
/* -- 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. */
/* .. Parameters .. */
/* .. Local Scalars .. */
/* .. External Subroutines .. */
/* .. Intrinsic Functions .. */
/* .. */
/* .. Executable Statements .. */
/* Test the input parameters. */
/* Parameter adjustments */
a_dim1 = *lda;
a_offset = a_dim1 + 1;
a -= a_offset;
--y;
--x;
/* 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_("ZGERU ", &info, 6L);
return 0;
}
/* Quick return if possible. */
if (*m == 0 || *n == 0 || alpha->r == 0. && alpha->i == 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) {
i__2 = jy;
if (y[i__2].r != 0. || y[i__2].i != 0.) {
i__2 = jy;
z__1.r = alpha->r * y[i__2].r - alpha->i * y[i__2].i, z__1.i =
alpha->r * y[i__2].i + alpha->i * y[i__2].r;
temp.r = z__1.r, temp.i = z__1.i;
i__2 = *m;
for (i = 1; i <= i__2; ++i) {
i__3 = i + j * a_dim1;
i__4 = i + j * a_dim1;
i__5 = i;
z__2.r = x[i__5].r * temp.r - x[i__5].i * temp.i, z__2.i =
x[i__5].r * temp.i + x[i__5].i * temp.r;
z__1.r = a[i__4].r + z__2.r, z__1.i = a[i__4].i + z__2.i;
a[i__3].r = z__1.r, a[i__3].i = z__1.i;
/* 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) {
i__2 = jy;
if (y[i__2].r != 0. || y[i__2].i != 0.) {
i__2 = jy;
z__1.r = alpha->r * y[i__2].r - alpha->i * y[i__2].i, z__1.i =
alpha->r * y[i__2].i + alpha->i * y[i__2].r;
temp.r = z__1.r, temp.i = z__1.i;
ix = kx;
i__2 = *m;
for (i = 1; i <= i__2; ++i) {
i__3 = i + j * a_dim1;
i__4 = i + j * a_dim1;
i__5 = ix;
z__2.r = x[i__5].r * temp.r - x[i__5].i * temp.i, z__2.i =
x[i__5].r * temp.i + x[i__5].i * temp.r;
z__1.r = a[i__4].r + z__2.r, z__1.i = a[i__4].i + z__2.i;
a[i__3].r = z__1.r, a[i__3].i = z__1.i;
ix += *incx;
/* L30: */
}
}
jy += *incy;
/* L40: */
}
}
return 0;
/* End of ZGERU . */
} /* zgeru_ */
/* zgemm.f -- translated by f2c (version of 23 April 1993 18:34:30).
You must link the resulting object file with the libraries:
-lf2c -lm (in that order)
*/
/* Subroutine */ int zgemm_(transa, transb, m, n, k, alpha, a, lda, b, ldb,
beta, c, ldc, transa_len, transb_len)
char *transa, *transb;
integer *m, *n, *k;
doublecomplex *alpha, *a;
integer *lda;
doublecomplex *b;
integer *ldb;
doublecomplex *beta, *c;
integer *ldc;
ftnlen transa_len;
ftnlen transb_len;
{
/* System generated locals */
integer a_dim1, a_offset, b_dim1, b_offset, c_dim1, c_offset, i__1, i__2,
i__3, i__4, i__5, i__6;
doublecomplex z__1, z__2, z__3, z__4;
/* Builtin functions */
void d_cnjg();
/* Local variables */
static integer info;
static logical nota, notb;
static doublecomplex temp;
static integer i, j, l;
static logical conja, conjb;
static integer ncola;
extern logical lsame_();
static integer nrowa, nrowb;
extern /* Subroutine */ int xerbla_();
/* .. Scalar Arguments .. */
/* .. Array Arguments .. */
/* .. */
/* Purpose */
/* ======= */
/* ZGEMM performs one of the matrix-matrix operations */
/* C := alpha*op( A )*op( B ) + beta*C, */
/* where op( X ) is one of */
/* op( X ) = X or op( X ) = X' or op( X ) = conjg( X' ), */
/* alpha and beta are scalars, and A, B and C are matrices, with op( A )
*/
/* an m by k matrix, op( B ) a k by n matrix and C an m by n matrix.
*/
/* Parameters */
/* ========== */
/* TRANSA - CHARACTER*1. */
/* On entry, TRANSA specifies the form of op( A ) to be used in
*/
/* the matrix multiplication as follows: */
/* TRANSA = 'N' or 'n', op( A ) = A. */
/* TRANSA = 'T' or 't', op( A ) = A'. */
/* TRANSA = 'C' or 'c', op( A ) = conjg( A' ). */
/* Unchanged on exit. */
/* TRANSB - CHARACTER*1. */
/* On entry, TRANSB specifies the form of op( B ) to be used in
*/
/* the matrix multiplication as follows: */
/* TRANSB = 'N' or 'n', op( B ) = B. */
/* TRANSB = 'T' or 't', op( B ) = B'. */
/* TRANSB = 'C' or 'c', op( B ) = conjg( B' ). */
/* Unchanged on exit. */
/* M - INTEGER. */
/* On entry, M specifies the number of rows of the matrix
*/
/* op( A ) and of the matrix C. M must be at least zero.
*/
/* Unchanged on exit. */
/* N - INTEGER. */
/* On entry, N specifies the number of columns of the matrix
*/
/* op( B ) and the number of columns of the matrix C. N must be
*/
/* at least zero. */
/* Unchanged on exit. */
/* K - INTEGER. */
/* On entry, K specifies the number of columns of the matrix
*/
/* op( A ) and the number of rows of the matrix op( B ). K must
*/
/* be at least zero. */
/* Unchanged on exit. */
/* ALPHA - COMPLEX*16 . */
/* On entry, ALPHA specifies the scalar alpha. */
/* Unchanged on exit. */
/* A - COMPLEX*16 array of DIMENSION ( LDA, ka ), where ka is
*/
/* k when TRANSA = 'N' or 'n', and is m otherwise. */
/* Before entry with TRANSA = 'N' or 'n', the leading m by k
*/
/* part of the array A must contain the matrix A, otherwise
*/
/* the leading k by m part of the array A must contain the
*/
/* matrix A. */
/* Unchanged on exit. */
/* LDA - INTEGER. */
/* On entry, LDA specifies the first dimension of A as declared
*/
/* in the calling (sub) program. When TRANSA = 'N' or 'n' then
*/
/* LDA must be at least max( 1, m ), otherwise LDA must be at
*/
/* least max( 1, k ). */
/* Unchanged on exit. */
/* B - COMPLEX*16 array of DIMENSION ( LDB, kb ), where kb is
*/
/* n when TRANSB = 'N' or 'n', and is k otherwise. */
/* Before entry with TRANSB = 'N' or 'n', the leading k by n
*/
/* part of the array B must contain the matrix B, otherwise
*/
/* the leading n by k part of the array B must contain the
*/
/* matrix B. */
/* Unchanged on exit. */
/* LDB - INTEGER. */
/* On entry, LDB specifies the first dimension of B as declared
*/
/* in the calling (sub) program. When TRANSB = 'N' or 'n' then
*/
/* LDB must be at least max( 1, k ), otherwise LDB must be at
*/
/* least max( 1, n ). */
/* Unchanged on exit. */
/* BETA - COMPLEX*16 . */
/* On entry, BETA specifies the scalar beta. When BETA is
*/
/* supplied as zero then C need not be set on input. */
/* Unchanged on exit. */
/* C - COMPLEX*16 array of DIMENSION ( LDC, n ). */
/* Before entry, the leading m by n part of the array C must
*/
/* contain the matrix C, except when beta is zero, in which
*/
/* case C need not be set on entry. */
/* On exit, the array C is overwritten by the m by n matrix
*/
/* ( alpha*op( A )*op( B ) + beta*C ). */
/* LDC - INTEGER. */
/* On entry, LDC specifies the first dimension of C as declared
*/
/* in the calling (sub) program. LDC must be at least
*/
/* max( 1, m ). */
/* Unchanged on exit. */
/* 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. */
/* .. External Functions .. */
/* .. External Subroutines .. */
/* .. Intrinsic Functions .. */
/* .. Local Scalars .. */
/* .. Parameters .. */
/* .. */
/* .. Executable Statements .. */
/* Set NOTA and NOTB as true if A and B respectively are not
*/
/* conjugated or transposed, set CONJA and CONJB as true if A and
*/
/* B respectively are to be transposed but not conjugated 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 */
c_dim1 = *ldc;
c_offset = c_dim1 + 1;
c -= c_offset;
b_dim1 = *ldb;
b_offset = b_dim1 + 1;
b -= b_offset;
a_dim1 = *lda;
a_offset = a_dim1 + 1;
a -= a_offset;
/* Function Body */
nota = lsame_(transa, "N", 1L, 1L);
notb = lsame_(transb, "N", 1L, 1L);
conja = lsame_(transa, "C", 1L, 1L);
conjb = lsame_(transb, "C", 1L, 1L);
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 && ! conja && ! lsame_(transa, "T", 1L, 1L)) {
info = 1;
} else if (! notb && ! conjb && ! lsame_(transb, "T", 1L, 1L)) {
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_("ZGEMM ", &info, 6L);
return 0;
}
/* Quick return if possible. */
if (*m == 0 || *n == 0 || (alpha->r == 0. && alpha->i == 0. || *k == 0) &&
(beta->r == 1. && beta->i == 0.)) {
return 0;
}
/* And when alpha.eq.zero. */
if (alpha->r == 0. && alpha->i == 0.) {
if (beta->r == 0. && beta->i == 0.) {
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
i__2 = *m;
for (i = 1; i <= i__2; ++i) {
i__3 = i + j * c_dim1;
c[i__3].r = 0., c[i__3].i = 0.;
/* L10: */
}
/* L20: */
}
} else {
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
i__2 = *m;
for (i = 1; i <= i__2; ++i) {
i__3 = i + j * c_dim1;
i__4 = i + j * c_dim1;
z__1.r = beta->r * c[i__4].r - beta->i * c[i__4].i,
z__1.i = beta->r * c[i__4].i + beta->i * c[i__4]
.r;
c[i__3].r = z__1.r, c[i__3].i = z__1.i;
/* 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->r == 0. && beta->i == 0.) {
i__2 = *m;
for (i = 1; i <= i__2; ++i) {
i__3 = i + j * c_dim1;
c[i__3].r = 0., c[i__3].i = 0.;
/* L50: */
}
} else if (beta->r != 1. || beta->i != 0.) {
i__2 = *m;
for (i = 1; i <= i__2; ++i) {
i__3 = i + j * c_dim1;
i__4 = i + j * c_dim1;
z__1.r = beta->r * c[i__4].r - beta->i * c[i__4].i,
z__1.i = beta->r * c[i__4].i + beta->i * c[
i__4].r;
c[i__3].r = z__1.r, c[i__3].i = z__1.i;
/* L60: */
}
}
i__2 = *k;
for (l = 1; l <= i__2; ++l) {
i__3 = l + j * b_dim1;
if (b[i__3].r != 0. || b[i__3].i != 0.) {
i__3 = l + j * b_dim1;
z__1.r = alpha->r * b[i__3].r - alpha->i * b[i__3].i,
z__1.i = alpha->r * b[i__3].i + alpha->i * b[
i__3].r;
temp.r = z__1.r, temp.i = z__1.i;
i__3 = *m;
for (i = 1; i <= i__3; ++i) {
i__4 = i + j * c_dim1;
i__5 = i + j * c_dim1;
i__6 = i + l * a_dim1;
z__2.r = temp.r * a[i__6].r - temp.i * a[i__6].i,
z__2.i = temp.r * a[i__6].i + temp.i * a[
i__6].r;
z__1.r = c[i__5].r + z__2.r, z__1.i = c[i__5].i +
z__2.i;
c[i__4].r = z__1.r, c[i__4].i = z__1.i;
/* L70: */
}
}
/* L80: */
}
/* L90: */
}
} else if (conja) {
/* Form C := alpha*conjg( 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.r = 0., temp.i = 0.;
i__3 = *k;
for (l = 1; l <= i__3; ++l) {
d_cnjg(&z__3, &a[l + i * a_dim1]);
i__4 = l + j * b_dim1;
z__2.r = z__3.r * b[i__4].r - z__3.i * b[i__4].i,
z__2.i = z__3.r * b[i__4].i + z__3.i * b[i__4]
.r;
z__1.r = temp.r + z__2.r, z__1.i = temp.i + z__2.i;
temp.r = z__1.r, temp.i = z__1.i;
/* L100: */
}
if (beta->r == 0. && beta->i == 0.) {
i__3 = i + j * c_dim1;
z__1.r = alpha->r * temp.r - alpha->i * temp.i,
z__1.i = alpha->r * temp.i + alpha->i *
temp.r;
c[i__3].r = z__1.r, c[i__3].i = z__1.i;
} else {
i__3 = i + j * c_dim1;
z__2.r = alpha->r * temp.r - alpha->i * temp.i,
z__2.i = alpha->r * temp.i + alpha->i *
temp.r;
i__4 = i + j * c_dim1;
z__3.r = beta->r * c[i__4].r - beta->i * c[i__4].i,
z__3.i = beta->r * c[i__4].i + beta->i * c[
i__4].r;
z__1.r = z__2.r + z__3.r, z__1.i = z__2.i + z__3.i;
c[i__3].r = z__1.r, c[i__3].i = z__1.i;
}
/* L110: */
}
/* L120: */
}
} 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.r = 0., temp.i = 0.;
i__3 = *k;
for (l = 1; l <= i__3; ++l) {
i__4 = l + i * a_dim1;
i__5 = l + j * b_dim1;
z__2.r = a[i__4].r * b[i__5].r - a[i__4].i * b[i__5]
.i, z__2.i = a[i__4].r * b[i__5].i + a[i__4]
.i * b[i__5].r;
z__1.r = temp.r + z__2.r, z__1.i = temp.i + z__2.i;
temp.r = z__1.r, temp.i = z__1.i;
/* L130: */
}
if (beta->r == 0. && beta->i == 0.) {
i__3 = i + j * c_dim1;
z__1.r = alpha->r * temp.r - alpha->i * temp.i,
z__1.i = alpha->r * temp.i + alpha->i *
temp.r;
c[i__3].r = z__1.r, c[i__3].i = z__1.i;
} else {
i__3 = i + j * c_dim1;
z__2.r = alpha->r * temp.r - alpha->i * temp.i,
z__2.i = alpha->r * temp.i + alpha->i *
temp.r;
i__4 = i + j * c_dim1;
z__3.r = beta->r * c[i__4].r - beta->i * c[i__4].i,
z__3.i = beta->r * c[i__4].i + beta->i * c[
i__4].r;
z__1.r = z__2.r + z__3.r, z__1.i = z__2.i + z__3.i;
c[i__3].r = z__1.r, c[i__3].i = z__1.i;
}
/* L140: */
}
/* L150: */
}
}
} else if (nota) {
if (conjb) {
/* Form C := alpha*A*conjg( B' ) + beta*C. */
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
if (beta->r == 0. && beta->i == 0.) {
i__2 = *m;
for (i = 1; i <= i__2; ++i) {
i__3 = i + j * c_dim1;
c[i__3].r = 0., c[i__3].i = 0.;
/* L160: */
}
} else if (beta->r != 1. || beta->i != 0.) {
i__2 = *m;
for (i = 1; i <= i__2; ++i) {
i__3 = i + j * c_dim1;
i__4 = i + j * c_dim1;
z__1.r = beta->r * c[i__4].r - beta->i * c[i__4].i,
z__1.i = beta->r * c[i__4].i + beta->i * c[
i__4].r;
c[i__3].r = z__1.r, c[i__3].i = z__1.i;
/* L170: */
}
}
i__2 = *k;
for (l = 1; l <= i__2; ++l) {
i__3 = j + l * b_dim1;
if (b[i__3].r != 0. || b[i__3].i != 0.) {
d_cnjg(&z__2, &b[j + l * b_dim1]);
z__1.r = alpha->r * z__2.r - alpha->i * z__2.i,
z__1.i = alpha->r * z__2.i + alpha->i *
z__2.r;
temp.r = z__1.r, temp.i = z__1.i;
i__3 = *m;
for (i = 1; i <= i__3; ++i) {
i__4 = i + j * c_dim1;
i__5 = i + j * c_dim1;
i__6 = i + l * a_dim1;
z__2.r = temp.r * a[i__6].r - temp.i * a[i__6].i,
z__2.i = temp.r * a[i__6].i + temp.i * a[
i__6].r;
z__1.r = c[i__5].r + z__2.r, z__1.i = c[i__5].i +
z__2.i;
c[i__4].r = z__1.r, c[i__4].i = z__1.i;
/* L180: */
}
}
/* L190: */
}
/* L200: */
}
} else {
/* Form C := alpha*A*B' + beta*C */
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
if (beta->r == 0. && beta->i == 0.) {
i__2 = *m;
for (i = 1; i <= i__2; ++i) {
i__3 = i + j * c_dim1;
c[i__3].r = 0., c[i__3].i = 0.;
/* L210: */
}
} else if (beta->r != 1. || beta->i != 0.) {
i__2 = *m;
for (i = 1; i <= i__2; ++i) {
i__3 = i + j * c_dim1;
i__4 = i + j * c_dim1;
z__1.r = beta->r * c[i__4].r - beta->i * c[i__4].i,
z__1.i = beta->r * c[i__4].i + beta->i * c[
i__4].r;
c[i__3].r = z__1.r, c[i__3].i = z__1.i;
/* L220: */
}
}
i__2 = *k;
for (l = 1; l <= i__2; ++l) {
i__3 = j + l * b_dim1;
if (b[i__3].r != 0. || b[i__3].i != 0.) {
i__3 = j + l * b_dim1;
z__1.r = alpha->r * b[i__3].r - alpha->i * b[i__3].i,
z__1.i = alpha->r * b[i__3].i + alpha->i * b[
i__3].r;
temp.r = z__1.r, temp.i = z__1.i;
i__3 = *m;
for (i = 1; i <= i__3; ++i) {
i__4 = i + j * c_dim1;
i__5 = i + j * c_dim1;
i__6 = i + l * a_dim1;
z__2.r = temp.r * a[i__6].r - temp.i * a[i__6].i,
z__2.i = temp.r * a[i__6].i + temp.i * a[
i__6].r;
z__1.r = c[i__5].r + z__2.r, z__1.i = c[i__5].i +
z__2.i;
c[i__4].r = z__1.r, c[i__4].i = z__1.i;
/* L230: */
}
}
/* L240: */
}
/* L250: */
}
}
} else if (conja) {
if (conjb) {
/* Form C := alpha*conjg( A' )*conjg( B' ) + beta*C. */
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
i__2 = *m;
for (i = 1; i <= i__2; ++i) {
temp.r = 0., temp.i = 0.;
i__3 = *k;
for (l = 1; l <= i__3; ++l) {
d_cnjg(&z__3, &a[l + i * a_dim1]);
d_cnjg(&z__4, &b[j + l * b_dim1]);
z__2.r = z__3.r * z__4.r - z__3.i * z__4.i, z__2.i =
z__3.r * z__4.i + z__3.i * z__4.r;
z__1.r = temp.r + z__2.r, z__1.i = temp.i + z__2.i;
temp.r = z__1.r, temp.i = z__1.i;
/* L260: */
}
if (beta->r == 0. && beta->i == 0.) {
i__3 = i + j * c_dim1;
z__1.r = alpha->r * temp.r - alpha->i * temp.i,
z__1.i = alpha->r * temp.i + alpha->i *
temp.r;
c[i__3].r = z__1.r, c[i__3].i = z__1.i;
} else {
i__3 = i + j * c_dim1;
z__2.r = alpha->r * temp.r - alpha->i * temp.i,
z__2.i = alpha->r * temp.i + alpha->i *
temp.r;
i__4 = i + j * c_dim1;
z__3.r = beta->r * c[i__4].r - beta->i * c[i__4].i,
z__3.i = beta->r * c[i__4].i + beta->i * c[
i__4].r;
z__1.r = z__2.r + z__3.r, z__1.i = z__2.i + z__3.i;
c[i__3].r = z__1.r, c[i__3].i = z__1.i;
}
/* L270: */
}
/* L280: */
}
} else {
/* Form C := alpha*conjg( 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.r = 0., temp.i = 0.;
i__3 = *k;
for (l = 1; l <= i__3; ++l) {
d_cnjg(&z__3, &a[l + i * a_dim1]);
i__4 = j + l * b_dim1;
z__2.r = z__3.r * b[i__4].r - z__3.i * b[i__4].i,
z__2.i = z__3.r * b[i__4].i + z__3.i * b[i__4]
.r;
z__1.r = temp.r + z__2.r, z__1.i = temp.i + z__2.i;
temp.r = z__1.r, temp.i = z__1.i;
/* L290: */
}
if (beta->r == 0. && beta->i == 0.) {
i__3 = i + j * c_dim1;
z__1.r = alpha->r * temp.r - alpha->i * temp.i,
z__1.i = alpha->r * temp.i + alpha->i *
temp.r;
c[i__3].r = z__1.r, c[i__3].i = z__1.i;
} else {
i__3 = i + j * c_dim1;
z__2.r = alpha->r * temp.r - alpha->i * temp.i,
z__2.i = alpha->r * temp.i + alpha->i *
temp.r;
i__4 = i + j * c_dim1;
z__3.r = beta->r * c[i__4].r - beta->i * c[i__4].i,
z__3.i = beta->r * c[i__4].i + beta->i * c[
i__4].r;
z__1.r = z__2.r + z__3.r, z__1.i = z__2.i + z__3.i;
c[i__3].r = z__1.r, c[i__3].i = z__1.i;
}
/* L300: */
}
/* L310: */
}
}
} else {
if (conjb) {
/* Form C := alpha*A'*conjg( B' ) + beta*C */
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
i__2 = *m;
for (i = 1; i <= i__2; ++i) {
temp.r = 0., temp.i = 0.;
i__3 = *k;
for (l = 1; l <= i__3; ++l) {
i__4 = l + i * a_dim1;
d_cnjg(&z__3, &b[j + l * b_dim1]);
z__2.r = a[i__4].r * z__3.r - a[i__4].i * z__3.i,
z__2.i = a[i__4].r * z__3.i + a[i__4].i *
z__3.r;
z__1.r = temp.r + z__2.r, z__1.i = temp.i + z__2.i;
temp.r = z__1.r, temp.i = z__1.i;
/* L320: */
}
if (beta->r == 0. && beta->i == 0.) {
i__3 = i + j * c_dim1;
z__1.r = alpha->r * temp.r - alpha->i * temp.i,
z__1.i = alpha->r * temp.i + alpha->i *
temp.r;
c[i__3].r = z__1.r, c[i__3].i = z__1.i;
} else {
i__3 = i + j * c_dim1;
z__2.r = alpha->r * temp.r - alpha->i * temp.i,
z__2.i = alpha->r * temp.i + alpha->i *
temp.r;
i__4 = i + j * c_dim1;
z__3.r = beta->r * c[i__4].r - beta->i * c[i__4].i,
z__3.i = beta->r * c[i__4].i + beta->i * c[
i__4].r;
z__1.r = z__2.r + z__3.r, z__1.i = z__2.i + z__3.i;
c[i__3].r = z__1.r, c[i__3].i = z__1.i;
}
/* L330: */
}
/* L340: */
}
} 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.r = 0., temp.i = 0.;
i__3 = *k;
for (l = 1; l <= i__3; ++l) {
i__4 = l + i * a_dim1;
i__5 = j + l * b_dim1;
z__2.r = a[i__4].r * b[i__5].r - a[i__4].i * b[i__5]
.i, z__2.i = a[i__4].r * b[i__5].i + a[i__4]
.i * b[i__5].r;
z__1.r = temp.r + z__2.r, z__1.i = temp.i + z__2.i;
temp.r = z__1.r, temp.i = z__1.i;
/* L350: */
}
if (beta->r == 0. && beta->i == 0.) {
i__3 = i + j * c_dim1;
z__1.r = alpha->r * temp.r - alpha->i * temp.i,
z__1.i = alpha->r * temp.i + alpha->i *
temp.r;
c[i__3].r = z__1.r, c[i__3].i = z__1.i;
} else {
i__3 = i + j * c_dim1;
z__2.r = alpha->r * temp.r - alpha->i * temp.i,
z__2.i = alpha->r * temp.i + alpha->i *
temp.r;
i__4 = i + j * c_dim1;
z__3.r = beta->r * c[i__4].r - beta->i * c[i__4].i,
z__3.i = beta->r * c[i__4].i + beta->i * c[
i__4].r;
z__1.r = z__2.r + z__3.r, z__1.i = z__2.i + z__3.i;
c[i__3].r = z__1.r, c[i__3].i = z__1.i;
}
/* L360: */
}
/* L370: */
}
}
}
return 0;
/* End of ZGEMM . */
} /* zgemm_ */
/* zcopy.f -- translated by f2c (version of 23 April 1993 18:34:30).
You must link the resulting object file with the libraries:
-lf2c -lm (in that order)
*/
/* Subroutine */ int zcopy_(n, zx, incx, zy, incy)
integer *n;
doublecomplex *zx;
integer *incx;
doublecomplex *zy;
integer *incy;
{
/* System generated locals */
integer i__1, i__2, i__3;
/* Local variables */
static integer i, ix, iy;
/* copies a vector, x, to a vector, y. */
/* jack dongarra, linpack, 4/11/78. */
/* Parameter adjustments */
--zy;
--zx;
/* 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) {
i__2 = iy;
i__3 = ix;
zy[i__2].r = zx[i__3].r, zy[i__2].i = zx[i__3].i;
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) {
i__2 = i;
i__3 = i;
zy[i__2].r = zx[i__3].r, zy[i__2].i = zx[i__3].i;
/* L30: */
}
return 0;
} /* zcopy_ */
/* xerbla.f -- translated by f2c (version of 23 April 1993 18:34:30).
You must link the resulting object file with the libraries:
-lf2c -lm (in that order)
*/
/* Subroutine */ int xerbla_(srname, info, srname_len)
char *srname;
integer *info;
ftnlen srname_len;
{
/* -- LAPACK auxiliary routine (preliminary version) -- */
/* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
/* Courant Institute, Argonne National Lab, and Rice University */
/* February 29, 1992 */
/* .. Scalar Arguments .. */
/* .. */
/* Purpose */
/* ======= */
/* XERBLA is an error handler for the LAPACK routines. */
/* It is called by an LAPACK routine if an input parameter has an */
/* invalid value. A message is printed and execution stops. */
/* Installers may consider modifying the STOP statement in order to */
/* call system-specific exception-handling facilities. */
/* Arguments */
/* ========= */
/* SRNAME (input) CHARACTER*6 */
/* The name of the routine which called XERBLA. */
/* INFO (input) INTEGER */
/* The position of the invalid parameter in the parameter list */
/* of the calling routine. */
/* WRITE( *, FMT = 9999 )SRNAME, INFO */
/* STOP */
/* 9999 FORMAT( ' ** On entry to ', A6, ' parameter number ', I2, ' had ',
*/
/* $ 'an illegal value' ) */
/* End of XERBLA */
return 0;
} /* xerbla_ */
/* ztrsm.f -- translated by f2c (version of 23 April 1993 18:34:30).
You must link the resulting object file with the libraries:
-lf2c -lm (in that order)
*/
/* Table of constant values */
static doublecomplex c_b20 = {1.,0.};
/* Subroutine */ int ztrsm_(side, uplo, transa, diag, m, n, alpha, a, lda, b,
ldb, side_len, uplo_len, transa_len, diag_len)
char *side, *uplo, *transa, *diag;
integer *m, *n;
doublecomplex *alpha, *a;
integer *lda;
doublecomplex *b;
integer *ldb;
ftnlen side_len;
ftnlen uplo_len;
ftnlen transa_len;
ftnlen diag_len;
{
/* System generated locals */
integer a_dim1, a_offset, b_dim1, b_offset, i__1, i__2, i__3, i__4, i__5,
i__6, i__7;
doublecomplex z__1, z__2, z__3;
/* Builtin functions */
void z_div(), d_cnjg();
/* Local variables */
static integer info;
static doublecomplex temp;
static integer i, j, k;
static logical lside;
extern logical lsame_();
static integer nrowa;
static logical upper;
extern /* Subroutine */ int xerbla_();
static logical noconj, nounit;
/* .. Scalar Arguments .. */
/* .. Array Arguments .. */
/* .. */
/* Purpose */
/* ======= */
/* ZTRSM solves one of the matrix equations */
/* op( A )*X = alpha*B, or X*op( A ) = alpha*B, */
/* where alpha is a scalar, X and B are m by n matrices, A is a unit, or
*/
/* non-unit, upper or lower triangular matrix and op( A ) is one of
*/
/* op( A ) = A or op( A ) = A' or op( A ) = conjg( A' ). */
/* The matrix X is overwritten on B. */
/* Parameters */
/* ========== */
/* SIDE - CHARACTER*1. */
/* On entry, SIDE specifies whether op( A ) appears on the left
*/
/* or right of X as follows: */
/* SIDE = 'L' or 'l' op( A )*X = alpha*B. */
/* SIDE = 'R' or 'r' X*op( A ) = alpha*B. */
/* Unchanged on exit. */
/* UPLO - CHARACTER*1. */
/* On entry, UPLO specifies whether the matrix A is an upper or
*/
/* lower triangular matrix as follows: */
/* UPLO = 'U' or 'u' A is an upper triangular matrix. */
/* UPLO = 'L' or 'l' A is a lower triangular matrix. */
/* Unchanged on exit. */
/* TRANSA - CHARACTER*1. */
/* On entry, TRANSA specifies the form of op( A ) to be used in
*/
/* the matrix multiplication as follows: */
/* TRANSA = 'N' or 'n' op( A ) = A. */
/* TRANSA = 'T' or 't' op( A ) = A'. */
/* TRANSA = 'C' or 'c' op( A ) = conjg( A' ). */
/* Unchanged on exit. */
/* DIAG - CHARACTER*1. */
/* On entry, DIAG specifies whether or not A is unit triangular
*/
/* as follows: */
/* DIAG = 'U' or 'u' A is assumed to be unit triangular. */
/* DIAG = 'N' or 'n' A is not assumed to be unit */
/* triangular. */
/* Unchanged on exit. */
/* M - INTEGER. */
/* On entry, M specifies the number of rows of B. M must be at
*/
/* least zero. */
/* Unchanged on exit. */
/* N - INTEGER. */
/* On entry, N specifies the number of columns of B. N must be
*/
/* at least zero. */
/* Unchanged on exit. */
/* ALPHA - COMPLEX*16 . */
/* On entry, ALPHA specifies the scalar alpha. When alpha is
*/
/* zero then A is not referenced and B need not be set before
*/
/* entry. */
/* Unchanged on exit. */
/* A - COMPLEX*16 array of DIMENSION ( LDA, k ), where k is m
*/
/* when SIDE = 'L' or 'l' and is n when SIDE = 'R' or 'r'.
*/
/* Before entry with UPLO = 'U' or 'u', the leading k by k
*/
/* upper triangular part of the array A must contain the upper
*/
/* triangular matrix and the strictly lower triangular part of
*/
/* A is not referenced. */
/* Before entry with UPLO = 'L' or 'l', the leading k by k
*/
/* lower triangular part of the array A must contain the lower
*/
/* triangular matrix and the strictly upper triangular part of
*/
/* A is not referenced. */
/* Note that when DIAG = 'U' or 'u', the diagonal elements of
*/
/* A are not referenced either, but are assumed to be unity.
*/
/* Unchanged on exit. */
/* LDA - INTEGER. */
/* On entry, LDA specifies the first dimension of A as declared
*/
/* in the calling (sub) program. When SIDE = 'L' or 'l' then
*/
/* LDA must be at least max( 1, m ), when SIDE = 'R' or 'r'
*/
/* then LDA must be at least max( 1, n ). */
/* Unchanged on exit. */
/* B - COMPLEX*16 array of DIMENSION ( LDB, n ). */
/* Before entry, the leading m by n part of the array B must
*/
/* contain the right-hand side matrix B, and on exit is
*/
/* overwritten by the solution matrix X. */
/* LDB - INTEGER. */
/* On entry, LDB specifies the first dimension of B as declared
*/
/* in the calling (sub) program. LDB must be at least
*/
/* max( 1, m ). */
/* Unchanged on exit. */
/* 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. */
/* .. External Functions .. */
/* .. External Subroutines .. */
/* .. Intrinsic Functions .. */
/* .. Local Scalars .. */
/* .. Parameters .. */
/* .. */
/* .. Executable Statements .. */
/* Test the input parameters. */
/* Parameter adjustments */
b_dim1 = *ldb;
b_offset = b_dim1 + 1;
b -= b_offset;
a_dim1 = *lda;
a_offset = a_dim1 + 1;
a -= a_offset;
/* Function Body */
lside = lsame_(side, "L", 1L, 1L);
if (lside) {
nrowa = *m;
} else {
nrowa = *n;
}
noconj = lsame_(transa, "T", 1L, 1L);
nounit = lsame_(diag, "N", 1L, 1L);
upper = lsame_(uplo, "U", 1L, 1L);
info = 0;
if (! lside && ! lsame_(side, "R", 1L, 1L)) {
info = 1;
} else if (! upper && ! lsame_(uplo, "L", 1L, 1L)) {
info = 2;
} else if (! lsame_(transa, "N", 1L, 1L) && ! lsame_(transa, "T", 1L, 1L)
&& ! lsame_(transa, "C", 1L, 1L)) {
info = 3;
} else if (! lsame_(diag, "U", 1L, 1L) && ! lsame_(diag, "N", 1L, 1L)) {
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_("ZTRSM ", &info, 6L);
return 0;
}
/* Quick return if possible. */
if (*n == 0) {
return 0;
}
/* And when alpha.eq.zero. */
if (alpha->r == 0. && alpha->i == 0.) {
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
i__2 = *m;
for (i = 1; i <= i__2; ++i) {
i__3 = i + j * b_dim1;
b[i__3].r = 0., b[i__3].i = 0.;
/* L10: */
}
/* L20: */
}
return 0;
}
/* Start the operations. */
if (lside) {
if (lsame_(transa, "N", 1L, 1L)) {
/* Form B := alpha*inv( A )*B. */
if (upper) {
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
if (alpha->r != 1. || alpha->i != 0.) {
i__2 = *m;
for (i = 1; i <= i__2; ++i) {
i__3 = i + j * b_dim1;
i__4 = i + j * b_dim1;
z__1.r = alpha->r * b[i__4].r - alpha->i * b[i__4]
.i, z__1.i = alpha->r * b[i__4].i +
alpha->i * b[i__4].r;
b[i__3].r = z__1.r, b[i__3].i = z__1.i;
/* L30: */
}
}
for (k = *m; k >= 1; --k) {
i__2 = k + j * b_dim1;
if (b[i__2].r != 0. || b[i__2].i != 0.) {
if (nounit) {
i__2 = k + j * b_dim1;
z_div(&z__1, &b[k + j * b_dim1], &a[k + k *
a_dim1]);
b[i__2].r = z__1.r, b[i__2].i = z__1.i;
}
i__2 = k - 1;
for (i = 1; i <= i__2; ++i) {
i__3 = i + j * b_dim1;
i__4 = i + j * b_dim1;
i__5 = k + j * b_dim1;
i__6 = i + k * a_dim1;
z__2.r = b[i__5].r * a[i__6].r - b[i__5].i *
a[i__6].i, z__2.i = b[i__5].r * a[
i__6].i + b[i__5].i * a[i__6].r;
z__1.r = b[i__4].r - z__2.r, z__1.i = b[i__4]
.i - z__2.i;
b[i__3].r = z__1.r, b[i__3].i = z__1.i;
/* L40: */
}
}
/* L50: */
}
/* L60: */
}
} else {
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
if (alpha->r != 1. || alpha->i != 0.) {
i__2 = *m;
for (i = 1; i <= i__2; ++i) {
i__3 = i + j * b_dim1;
i__4 = i + j * b_dim1;
z__1.r = alpha->r * b[i__4].r - alpha->i * b[i__4]
.i, z__1.i = alpha->r * b[i__4].i +
alpha->i * b[i__4].r;
b[i__3].r = z__1.r, b[i__3].i = z__1.i;
/* L70: */
}
}
i__2 = *m;
for (k = 1; k <= i__2; ++k) {
i__3 = k + j * b_dim1;
if (b[i__3].r != 0. || b[i__3].i != 0.) {
if (nounit) {
i__3 = k + j * b_dim1;
z_div(&z__1, &b[k + j * b_dim1], &a[k + k *
a_dim1]);
b[i__3].r = z__1.r, b[i__3].i = z__1.i;
}
i__3 = *m;
for (i = k + 1; i <= i__3; ++i) {
i__4 = i + j * b_dim1;
i__5 = i + j * b_dim1;
i__6 = k + j * b_dim1;
i__7 = i + k * a_dim1;
z__2.r = b[i__6].r * a[i__7].r - b[i__6].i *
a[i__7].i, z__2.i = b[i__6].r * a[
i__7].i + b[i__6].i * a[i__7].r;
z__1.r = b[i__5].r - z__2.r, z__1.i = b[i__5]
.i - z__2.i;
b[i__4].r = z__1.r, b[i__4].i = z__1.i;
/* L80: */
}
}
/* L90: */
}
/* L100: */
}
}
} else {
/* Form B := alpha*inv( A' )*B */
/* or B := alpha*inv( conjg( A' ) )*B. */
if (upper) {
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
i__2 = *m;
for (i = 1; i <= i__2; ++i) {
i__3 = i + j * b_dim1;
z__1.r = alpha->r * b[i__3].r - alpha->i * b[i__3].i,
z__1.i = alpha->r * b[i__3].i + alpha->i * b[
i__3].r;
temp.r = z__1.r, temp.i = z__1.i;
if (noconj) {
i__3 = i - 1;
for (k = 1; k <= i__3; ++k) {
i__4 = k + i * a_dim1;
i__5 = k + j * b_dim1;
z__2.r = a[i__4].r * b[i__5].r - a[i__4].i *
b[i__5].i, z__2.i = a[i__4].r * b[
i__5].i + a[i__4].i * b[i__5].r;
z__1.r = temp.r - z__2.r, z__1.i = temp.i -
z__2.i;
temp.r = z__1.r, temp.i = z__1.i;
/* L110: */
}
if (nounit) {
z_div(&z__1, &temp, &a[i + i * a_dim1]);
temp.r = z__1.r, temp.i = z__1.i;
}
} else {
i__3 = i - 1;
for (k = 1; k <= i__3; ++k) {
d_cnjg(&z__3, &a[k + i * a_dim1]);
i__4 = k + j * b_dim1;
z__2.r = z__3.r * b[i__4].r - z__3.i * b[i__4]
.i, z__2.i = z__3.r * b[i__4].i +
z__3.i * b[i__4].r;
z__1.r = temp.r - z__2.r, z__1.i = temp.i -
z__2.i;
temp.r = z__1.r, temp.i = z__1.i;
/* L120: */
}
if (nounit) {
d_cnjg(&z__2, &a[i + i * a_dim1]);
z_div(&z__1, &temp, &z__2);
temp.r = z__1.r, temp.i = z__1.i;
}
}
i__3 = i + j * b_dim1;
b[i__3].r = temp.r, b[i__3].i = temp.i;
/* L130: */
}
/* L140: */
}
} else {
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
for (i = *m; i >= 1; --i) {
i__2 = i + j * b_dim1;
z__1.r = alpha->r * b[i__2].r - alpha->i * b[i__2].i,
z__1.i = alpha->r * b[i__2].i + alpha->i * b[
i__2].r;
temp.r = z__1.r, temp.i = z__1.i;
if (noconj) {
i__2 = *m;
for (k = i + 1; k <= i__2; ++k) {
i__3 = k + i * a_dim1;
i__4 = k + j * b_dim1;
z__2.r = a[i__3].r * b[i__4].r - a[i__3].i *
b[i__4].i, z__2.i = a[i__3].r * b[
i__4].i + a[i__3].i * b[i__4].r;
z__1.r = temp.r - z__2.r, z__1.i = temp.i -
z__2.i;
temp.r = z__1.r, temp.i = z__1.i;
/* L150: */
}
if (nounit) {
z_div(&z__1, &temp, &a[i + i * a_dim1]);
temp.r = z__1.r, temp.i = z__1.i;
}
} else {
i__2 = *m;
for (k = i + 1; k <= i__2; ++k) {
d_cnjg(&z__3, &a[k + i * a_dim1]);
i__3 = k + j * b_dim1;
z__2.r = z__3.r * b[i__3].r - z__3.i * b[i__3]
.i, z__2.i = z__3.r * b[i__3].i +
z__3.i * b[i__3].r;
z__1.r = temp.r - z__2.r, z__1.i = temp.i -
z__2.i;
temp.r = z__1.r, temp.i = z__1.i;
/* L160: */
}
if (nounit) {
d_cnjg(&z__2, &a[i + i * a_dim1]);
z_div(&z__1, &temp, &z__2);
temp.r = z__1.r, temp.i = z__1.i;
}
}
i__2 = i + j * b_dim1;
b[i__2].r = temp.r, b[i__2].i = temp.i;
/* L170: */
}
/* L180: */
}
}
}
} else {
if (lsame_(transa, "N", 1L, 1L)) {
/* Form B := alpha*B*inv( A ). */
if (upper) {
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
if (alpha->r != 1. || alpha->i != 0.) {
i__2 = *m;
for (i = 1; i <= i__2; ++i) {
i__3 = i + j * b_dim1;
i__4 = i + j * b_dim1;
z__1.r = alpha->r * b[i__4].r - alpha->i * b[i__4]
.i, z__1.i = alpha->r * b[i__4].i +
alpha->i * b[i__4].r;
b[i__3].r = z__1.r, b[i__3].i = z__1.i;
/* L190: */
}
}
i__2 = j - 1;
for (k = 1; k <= i__2; ++k) {
i__3 = k + j * a_dim1;
if (a[i__3].r != 0. || a[i__3].i != 0.) {
i__3 = *m;
for (i = 1; i <= i__3; ++i) {
i__4 = i + j * b_dim1;
i__5 = i + j * b_dim1;
i__6 = k + j * a_dim1;
i__7 = i + k * b_dim1;
z__2.r = a[i__6].r * b[i__7].r - a[i__6].i *
b[i__7].i, z__2.i = a[i__6].r * b[
i__7].i + a[i__6].i * b[i__7].r;
z__1.r = b[i__5].r - z__2.r, z__1.i = b[i__5]
.i - z__2.i;
b[i__4].r = z__1.r, b[i__4].i = z__1.i;
/* L200: */
}
}
/* L210: */
}
if (nounit) {
z_div(&z__1, &c_b20, &a[j + j * a_dim1]);
temp.r = z__1.r, temp.i = z__1.i;
i__2 = *m;
for (i = 1; i <= i__2; ++i) {
i__3 = i + j * b_dim1;
i__4 = i + j * b_dim1;
z__1.r = temp.r * b[i__4].r - temp.i * b[i__4].i,
z__1.i = temp.r * b[i__4].i + temp.i * b[
i__4].r;
b[i__3].r = z__1.r, b[i__3].i = z__1.i;
/* L220: */
}
}
/* L230: */
}
} else {
for (j = *n; j >= 1; --j) {
if (alpha->r != 1. || alpha->i != 0.) {
i__1 = *m;
for (i = 1; i <= i__1; ++i) {
i__2 = i + j * b_dim1;
i__3 = i + j * b_dim1;
z__1.r = alpha->r * b[i__3].r - alpha->i * b[i__3]
.i, z__1.i = alpha->r * b[i__3].i +
alpha->i * b[i__3].r;
b[i__2].r = z__1.r, b[i__2].i = z__1.i;
/* L240: */
}
}
i__1 = *n;
for (k = j + 1; k <= i__1; ++k) {
i__2 = k + j * a_dim1;
if (a[i__2].r != 0. || a[i__2].i != 0.) {
i__2 = *m;
for (i = 1; i <= i__2; ++i) {
i__3 = i + j * b_dim1;
i__4 = i + j * b_dim1;
i__5 = k + j * a_dim1;
i__6 = i + k * b_dim1;
z__2.r = a[i__5].r * b[i__6].r - a[i__5].i *
b[i__6].i, z__2.i = a[i__5].r * b[
i__6].i + a[i__5].i * b[i__6].r;
z__1.r = b[i__4].r - z__2.r, z__1.i = b[i__4]
.i - z__2.i;
b[i__3].r = z__1.r, b[i__3].i = z__1.i;
/* L250: */
}
}
/* L260: */
}
if (nounit) {
z_div(&z__1, &c_b20, &a[j + j * a_dim1]);
temp.r = z__1.r, temp.i = z__1.i;
i__1 = *m;
for (i = 1; i <= i__1; ++i) {
i__2 = i + j * b_dim1;
i__3 = i + j * b_dim1;
z__1.r = temp.r * b[i__3].r - temp.i * b[i__3].i,
z__1.i = temp.r * b[i__3].i + temp.i * b[
i__3].r;
b[i__2].r = z__1.r, b[i__2].i = z__1.i;
/* L270: */
}
}
/* L280: */
}
}
} else {
/* Form B := alpha*B*inv( A' ) */
/* or B := alpha*B*inv( conjg( A' ) ). */
if (upper) {
for (k = *n; k >= 1; --k) {
if (nounit) {
if (noconj) {
z_div(&z__1, &c_b20, &a[k + k * a_dim1]);
temp.r = z__1.r, temp.i = z__1.i;
} else {
d_cnjg(&z__2, &a[k + k * a_dim1]);
z_div(&z__1, &c_b20, &z__2);
temp.r = z__1.r, temp.i = z__1.i;
}
i__1 = *m;
for (i = 1; i <= i__1; ++i) {
i__2 = i + k * b_dim1;
i__3 = i + k * b_dim1;
z__1.r = temp.r * b[i__3].r - temp.i * b[i__3].i,
z__1.i = temp.r * b[i__3].i + temp.i * b[
i__3].r;
b[i__2].r = z__1.r, b[i__2].i = z__1.i;
/* L290: */
}
}
i__1 = k - 1;
for (j = 1; j <= i__1; ++j) {
i__2 = j + k * a_dim1;
if (a[i__2].r != 0. || a[i__2].i != 0.) {
if (noconj) {
i__2 = j + k * a_dim1;
temp.r = a[i__2].r, temp.i = a[i__2].i;
} else {
d_cnjg(&z__1, &a[j + k * a_dim1]);
temp.r = z__1.r, temp.i = z__1.i;
}
i__2 = *m;
for (i = 1; i <= i__2; ++i) {
i__3 = i + j * b_dim1;
i__4 = i + j * b_dim1;
i__5 = i + k * b_dim1;
z__2.r = temp.r * b[i__5].r - temp.i * b[i__5]
.i, z__2.i = temp.r * b[i__5].i +
temp.i * b[i__5].r;
z__1.r = b[i__4].r - z__2.r, z__1.i = b[i__4]
.i - z__2.i;
b[i__3].r = z__1.r, b[i__3].i = z__1.i;
/* L300: */
}
}
/* L310: */
}
if (alpha->r != 1. || alpha->i != 0.) {
i__1 = *m;
for (i = 1; i <= i__1; ++i) {
i__2 = i + k * b_dim1;
i__3 = i + k * b_dim1;
z__1.r = alpha->r * b[i__3].r - alpha->i * b[i__3]
.i, z__1.i = alpha->r * b[i__3].i +
alpha->i * b[i__3].r;
b[i__2].r = z__1.r, b[i__2].i = z__1.i;
/* L320: */
}
}
/* L330: */
}
} else {
i__1 = *n;
for (k = 1; k <= i__1; ++k) {
if (nounit) {
if (noconj) {
z_div(&z__1, &c_b20, &a[k + k * a_dim1]);
temp.r = z__1.r, temp.i = z__1.i;
} else {
d_cnjg(&z__2, &a[k + k * a_dim1]);
z_div(&z__1, &c_b20, &z__2);
temp.r = z__1.r, temp.i = z__1.i;
}
i__2 = *m;
for (i = 1; i <= i__2; ++i) {
i__3 = i + k * b_dim1;
i__4 = i + k * b_dim1;
z__1.r = temp.r * b[i__4].r - temp.i * b[i__4].i,
z__1.i = temp.r * b[i__4].i + temp.i * b[
i__4].r;
b[i__3].r = z__1.r, b[i__3].i = z__1.i;
/* L340: */
}
}
i__2 = *n;
for (j = k + 1; j <= i__2; ++j) {
i__3 = j + k * a_dim1;
if (a[i__3].r != 0. || a[i__3].i != 0.) {
if (noconj) {
i__3 = j + k * a_dim1;
temp.r = a[i__3].r, temp.i = a[i__3].i;
} else {
d_cnjg(&z__1, &a[j + k * a_dim1]);
temp.r = z__1.r, temp.i = z__1.i;
}
i__3 = *m;
for (i = 1; i <= i__3; ++i) {
i__4 = i + j * b_dim1;
i__5 = i + j * b_dim1;
i__6 = i + k * b_dim1;
z__2.r = temp.r * b[i__6].r - temp.i * b[i__6]
.i, z__2.i = temp.r * b[i__6].i +
temp.i * b[i__6].r;
z__1.r = b[i__5].r - z__2.r, z__1.i = b[i__5]
.i - z__2.i;
b[i__4].r = z__1.r, b[i__4].i = z__1.i;
/* L350: */
}
}
/* L360: */
}
if (alpha->r != 1. || alpha->i != 0.) {
i__2 = *m;
for (i = 1; i <= i__2; ++i) {
i__3 = i + k * b_dim1;
i__4 = i + k * b_dim1;
z__1.r = alpha->r * b[i__4].r - alpha->i * b[i__4]
.i, z__1.i = alpha->r * b[i__4].i +
alpha->i * b[i__4].r;
b[i__3].r = z__1.r, b[i__3].i = z__1.i;
/* L370: */
}
}
/* L380: */
}
}
}
}
return 0;
/* End of ZTRSM . */
} /* ztrsm_ */
/* dzasum.f -- translated by f2c (version of 23 April 1993 18:34:30).
You must link the resulting object file with the libraries:
-lf2c -lm (in that order)
*/
doublereal dzasum_(n, zx, incx)
integer *n;
doublecomplex *zx;
integer *incx;
{
/* System generated locals */
integer i__1;
doublereal ret_val;
/* Local variables */
static integer i;
static doublereal stemp;
extern doublereal dcabs1_();
static integer ix;
/* takes the sum of the absolute values. */
/* jack dongarra, 3/11/78. */
/* modified to correct problem with negative increment, 8/21/90. */
/* Parameter adjustments */
--zx;
/* Function Body */
ret_val = 0.;
stemp = 0.;
if (*n <= 0) {
return ret_val;
}
if (*incx == 1) {
goto L20;
}
/* code for increment not equal to 1 */
ix = 1;
if (*incx < 0) {
ix = (-(*n) + 1) * *incx + 1;
}
i__1 = *n;
for (i = 1; i <= i__1; ++i) {
stemp += dcabs1_(&zx[ix]);
ix += *incx;
/* L10: */
}
ret_val = stemp;
return ret_val;
/* code for increment equal to 1 */
L20:
i__1 = *n;
for (i = 1; i <= i__1; ++i) {
stemp += dcabs1_(&zx[i]);
/* L30: */
}
ret_val = stemp;
return ret_val;
} /* dzasum_ */
/* dswap.f -- translated by f2c (version of 23 April 1993 18:34:30).
You must link the resulting object file with the libraries:
-lf2c -lm (in that order)
*/
/* Subroutine */ int dswap_(n, dx, incx, dy, incy)
integer *n;
doublereal *dx;
integer *incx;
doublereal *dy;
integer *incy;
{
/* System generated locals */
integer i__1;
/* Local variables */
static integer i, m;
static doublereal dtemp;
static integer ix, iy, mp1;
/* interchanges two vectors. */
/* uses unrolled loops for increments equal one. */
/* jack dongarra, linpack, 3/11/78. */
/* 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;
} /* dswap_ */
/* dscal.f -- translated by f2c (version of 23 April 1993 18:34:30).
You must link the resulting object file with the libraries:
-lf2c -lm (in that order)
*/
/* Subroutine */ int dscal_(n, da, dx, incx)
integer *n;
doublereal *da, *dx;
integer *incx;
{
/* System generated locals */
integer i__1;
/* Local variables */
static integer i, m, ix, mp1;
/* scales a vector by a constant. */
/* uses unrolled loops for increment equal to one. */
/* jack dongarra, linpack, 3/11/78. */
/* modified to correct problem with negative increment, 8/21/90. */
/* Parameter adjustments */
--dx;
/* Function Body */
if (*n <= 0) {
return 0;
}
if (*incx == 1) {
goto L20;
}
/* code for increment not equal to 1 */
ix = 1;
if (*incx < 0) {
ix = (-(*n) + 1) * *incx + 1;
}
i__1 = *n;
for (i = 1; i <= i__1; ++i) {
dx[ix] = *da * dx[ix];
ix += *incx;
/* L10: */
}
return 0;
/* code for increment 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) {
dx[i] = *da * dx[i];
/* L30: */
}
if (*n < 5) {
return 0;
}
L40:
mp1 = m + 1;
i__1 = *n;
for (i = mp1; i <= i__1; 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_ */
/* ztrmv.f -- translated by f2c (version of 23 April 1993 18:34:30).
You must link the resulting object file with the libraries:
-lf2c -lm (in that order)
*/
/* Subroutine */ int ztrmv_(uplo, trans, diag, n, a, lda, x, incx, uplo_len,
trans_len, diag_len)
char *uplo, *trans, *diag;
integer *n;
doublecomplex *a;
integer *lda;
doublecomplex *x;
integer *incx;
ftnlen uplo_len;
ftnlen trans_len;
ftnlen diag_len;
{
/* System generated locals */
integer a_dim1, a_offset, i__1, i__2, i__3, i__4, i__5;
doublecomplex z__1, z__2, z__3;
/* Builtin functions */
void d_cnjg();
/* Local variables */
static integer info;
static doublecomplex temp;
static integer i, j;
extern logical lsame_();
static integer ix, jx, kx;
extern /* Subroutine */ int xerbla_();
static logical noconj, nounit;
/* .. Scalar Arguments .. */
/* .. Array Arguments .. */
/* .. */
/* Purpose */
/* ======= */
/* ZTRMV performs one of the matrix-vector operations */
/* x := A*x, or x := A'*x, or x := conjg( A' )*x, */
/* where x is an n element vector and A is an n by n unit, or non-unit,
*/
/* upper or lower triangular matrix. */
/* Parameters */
/* ========== */
/* UPLO - CHARACTER*1. */
/* On entry, UPLO specifies whether the matrix is an upper or */
/* lower triangular matrix as follows: */
/* UPLO = 'U' or 'u' A is an upper triangular matrix. */
/* UPLO = 'L' or 'l' A is a lower triangular matrix. */
/* Unchanged on exit. */
/* TRANS - CHARACTER*1. */
/* On entry, TRANS specifies the operation to be performed as */
/* follows: */
/* TRANS = 'N' or 'n' x := A*x. */
/* TRANS = 'T' or 't' x := A'*x. */
/* TRANS = 'C' or 'c' x := conjg( A' )*x. */
/* Unchanged on exit. */
/* DIAG - CHARACTER*1. */
/* On entry, DIAG specifies whether or not A is unit */
/* triangular as follows: */
/* DIAG = 'U' or 'u' A is assumed to be unit triangular. */
/* DIAG = 'N' or 'n' A is not assumed to be unit */
/* triangular. */
/* Unchanged on exit. */
/* N - INTEGER. */
/* On entry, N specifies the order of the matrix A. */
/* N must be at least zero. */
/* Unchanged on exit. */
/* A - COMPLEX*16 array of DIMENSION ( LDA, n ). */
/* Before entry with UPLO = 'U' or 'u', the leading n by n */
/* upper triangular part of the array A must contain the upper
*/
/* triangular matrix and the strictly lower triangular part of
*/
/* A is not referenced. */
/* Before entry with UPLO = 'L' or 'l', the leading n by n */
/* lower triangular part of the array A must contain the lower
*/
/* triangular matrix and the strictly upper triangular part of
*/
/* A is not referenced. */
/* Note that when DIAG = 'U' or 'u', the diagonal elements of
*/
/* A are not referenced either, but are assumed to be unity. */
/* Unchanged on exit. */
/* LDA - INTEGER. */
/* On entry, LDA specifies the first dimension of A as declared
*/
/* in the calling (sub) program. LDA must be at least */
/* max( 1, n ). */
/* Unchanged on exit. */
/* X - COMPLEX*16 array of dimension at least */
/* ( 1 + ( n - 1 )*abs( INCX ) ). */
/* Before entry, the incremented array X must contain the n */
/* element vector x. On exit, X is overwritten with the */
/* tranformed vector x. */
/* INCX - INTEGER. */
/* On entry, INCX specifies the increment for the elements of */
/* X. INCX must not be zero. */
/* Unchanged on exit. */
/* Level 2 Blas routine. */
/* -- 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. */
/* .. Parameters .. */
/* .. Local Scalars .. */
/* .. External Functions .. */
/* .. External Subroutines .. */
/* .. Intrinsic Functions .. */
/* .. */
/* .. Executable Statements .. */
/* Test the input parameters. */
/* Parameter adjustments */
--x;
a_dim1 = *lda;
a_offset = a_dim1 + 1;
a -= a_offset;
/* Function Body */
info = 0;
if (! lsame_(uplo, "U", 1L, 1L) && ! lsame_(uplo, "L", 1L, 1L)) {
info = 1;
} else if (! lsame_(trans, "N", 1L, 1L) && ! lsame_(trans, "T", 1L, 1L) &&
! lsame_(trans, "C", 1L, 1L)) {
info = 2;
} else if (! lsame_(diag, "U", 1L, 1L) && ! lsame_(diag, "N", 1L, 1L)) {
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_("ZTRMV ", &info, 6L);
return 0;
}
/* Quick return if possible. */
if (*n == 0) {
return 0;
}
noconj = lsame_(trans, "T", 1L, 1L);
nounit = lsame_(diag, "N", 1L, 1L);
/* 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", 1L, 1L)) {
/* Form x := A*x. */
if (lsame_(uplo, "U", 1L, 1L)) {
if (*incx == 1) {
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
i__2 = j;
if (x[i__2].r != 0. || x[i__2].i != 0.) {
i__2 = j;
temp.r = x[i__2].r, temp.i = x[i__2].i;
i__2 = j - 1;
for (i = 1; i <= i__2; ++i) {
i__3 = i;
i__4 = i;
i__5 = i + j * a_dim1;
z__2.r = temp.r * a[i__5].r - temp.i * a[i__5].i,
z__2.i = temp.r * a[i__5].i + temp.i * a[
i__5].r;
z__1.r = x[i__4].r + z__2.r, z__1.i = x[i__4].i +
z__2.i;
x[i__3].r = z__1.r, x[i__3].i = z__1.i;
/* L10: */
}
if (nounit) {
i__2 = j;
i__3 = j;
i__4 = j + j * a_dim1;
z__1.r = x[i__3].r * a[i__4].r - x[i__3].i * a[
i__4].i, z__1.i = x[i__3].r * a[i__4].i +
x[i__3].i * a[i__4].r;
x[i__2].r = z__1.r, x[i__2].i = z__1.i;
}
}
/* L20: */
}
} else {
jx = kx;
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
i__2 = jx;
if (x[i__2].r != 0. || x[i__2].i != 0.) {
i__2 = jx;
temp.r = x[i__2].r, temp.i = x[i__2].i;
ix = kx;
i__2 = j - 1;
for (i = 1; i <= i__2; ++i) {
i__3 = ix;
i__4 = ix;
i__5 = i + j * a_dim1;
z__2.r = temp.r * a[i__5].r - temp.i * a[i__5].i,
z__2.i = temp.r * a[i__5].i + temp.i * a[
i__5].r;
z__1.r = x[i__4].r + z__2.r, z__1.i = x[i__4].i +
z__2.i;
x[i__3].r = z__1.r, x[i__3].i = z__1.i;
ix += *incx;
/* L30: */
}
if (nounit) {
i__2 = jx;
i__3 = jx;
i__4 = j + j * a_dim1;
z__1.r = x[i__3].r * a[i__4].r - x[i__3].i * a[
i__4].i, z__1.i = x[i__3].r * a[i__4].i +
x[i__3].i * a[i__4].r;
x[i__2].r = z__1.r, x[i__2].i = z__1.i;
}
}
jx += *incx;
/* L40: */
}
}
} else {
if (*incx == 1) {
for (j = *n; j >= 1; --j) {
i__1 = j;
if (x[i__1].r != 0. || x[i__1].i != 0.) {
i__1 = j;
temp.r = x[i__1].r, temp.i = x[i__1].i;
i__1 = j + 1;
for (i = *n; i >= i__1; --i) {
i__2 = i;
i__3 = i;
i__4 = i + j * a_dim1;
z__2.r = temp.r * a[i__4].r - temp.i * a[i__4].i,
z__2.i = temp.r * a[i__4].i + temp.i * a[
i__4].r;
z__1.r = x[i__3].r + z__2.r, z__1.i = x[i__3].i +
z__2.i;
x[i__2].r = z__1.r, x[i__2].i = z__1.i;
/* L50: */
}
if (nounit) {
i__1 = j;
i__2 = j;
i__3 = j + j * a_dim1;
z__1.r = x[i__2].r * a[i__3].r - x[i__2].i * a[
i__3].i, z__1.i = x[i__2].r * a[i__3].i +
x[i__2].i * a[i__3].r;
x[i__1].r = z__1.r, x[i__1].i = z__1.i;
}
}
/* L60: */
}
} else {
kx += (*n - 1) * *incx;
jx = kx;
for (j = *n; j >= 1; --j) {
i__1 = jx;
if (x[i__1].r != 0. || x[i__1].i != 0.) {
i__1 = jx;
temp.r = x[i__1].r, temp.i = x[i__1].i;
ix = kx;
i__1 = j + 1;
for (i = *n; i >= i__1; --i) {
i__2 = ix;
i__3 = ix;
i__4 = i + j * a_dim1;
z__2.r = temp.r * a[i__4].r - temp.i * a[i__4].i,
z__2.i = temp.r * a[i__4].i + temp.i * a[
i__4].r;
z__1.r = x[i__3].r + z__2.r, z__1.i = x[i__3].i +
z__2.i;
x[i__2].r = z__1.r, x[i__2].i = z__1.i;
ix -= *incx;
/* L70: */
}
if (nounit) {
i__1 = jx;
i__2 = jx;
i__3 = j + j * a_dim1;
z__1.r = x[i__2].r * a[i__3].r - x[i__2].i * a[
i__3].i, z__1.i = x[i__2].r * a[i__3].i +
x[i__2].i * a[i__3].r;
x[i__1].r = z__1.r, x[i__1].i = z__1.i;
}
}
jx -= *incx;
/* L80: */
}
}
}
} else {
/* Form x := A'*x or x := conjg( A' )*x. */
if (lsame_(uplo, "U", 1L, 1L)) {
if (*incx == 1) {
for (j = *n; j >= 1; --j) {
i__1 = j;
temp.r = x[i__1].r, temp.i = x[i__1].i;
if (noconj) {
if (nounit) {
i__1 = j + j * a_dim1;
z__1.r = temp.r * a[i__1].r - temp.i * a[i__1].i,
z__1.i = temp.r * a[i__1].i + temp.i * a[
i__1].r;
temp.r = z__1.r, temp.i = z__1.i;
}
for (i = j - 1; i >= 1; --i) {
i__1 = i + j * a_dim1;
i__2 = i;
z__2.r = a[i__1].r * x[i__2].r - a[i__1].i * x[
i__2].i, z__2.i = a[i__1].r * x[i__2].i +
a[i__1].i * x[i__2].r;
z__1.r = temp.r + z__2.r, z__1.i = temp.i +
z__2.i;
temp.r = z__1.r, temp.i = z__1.i;
/* L90: */
}
} else {
if (nounit) {
d_cnjg(&z__2, &a[j + j * a_dim1]);
z__1.r = temp.r * z__2.r - temp.i * z__2.i,
z__1.i = temp.r * z__2.i + temp.i *
z__2.r;
temp.r = z__1.r, temp.i = z__1.i;
}
for (i = j - 1; i >= 1; --i) {
d_cnjg(&z__3, &a[i + j * a_dim1]);
i__1 = i;
z__2.r = z__3.r * x[i__1].r - z__3.i * x[i__1].i,
z__2.i = z__3.r * x[i__1].i + z__3.i * x[
i__1].r;
z__1.r = temp.r + z__2.r, z__1.i = temp.i +
z__2.i;
temp.r = z__1.r, temp.i = z__1.i;
/* L100: */
}
}
i__1 = j;
x[i__1].r = temp.r, x[i__1].i = temp.i;
/* L110: */
}
} else {
jx = kx + (*n - 1) * *incx;
for (j = *n; j >= 1; --j) {
i__1 = jx;
temp.r = x[i__1].r, temp.i = x[i__1].i;
ix = jx;
if (noconj) {
if (nounit) {
i__1 = j + j * a_dim1;
z__1.r = temp.r * a[i__1].r - temp.i * a[i__1].i,
z__1.i = temp.r * a[i__1].i + temp.i * a[
i__1].r;
temp.r = z__1.r, temp.i = z__1.i;
}
for (i = j - 1; i >= 1; --i) {
ix -= *incx;
i__1 = i + j * a_dim1;
i__2 = ix;
z__2.r = a[i__1].r * x[i__2].r - a[i__1].i * x[
i__2].i, z__2.i = a[i__1].r * x[i__2].i +
a[i__1].i * x[i__2].r;
z__1.r = temp.r + z__2.r, z__1.i = temp.i +
z__2.i;
temp.r = z__1.r, temp.i = z__1.i;
/* L120: */
}
} else {
if (nounit) {
d_cnjg(&z__2, &a[j + j * a_dim1]);
z__1.r = temp.r * z__2.r - temp.i * z__2.i,
z__1.i = temp.r * z__2.i + temp.i *
z__2.r;
temp.r = z__1.r, temp.i = z__1.i;
}
for (i = j - 1; i >= 1; --i) {
ix -= *incx;
d_cnjg(&z__3, &a[i + j * a_dim1]);
i__1 = ix;
z__2.r = z__3.r * x[i__1].r - z__3.i * x[i__1].i,
z__2.i = z__3.r * x[i__1].i + z__3.i * x[
i__1].r;
z__1.r = temp.r + z__2.r, z__1.i = temp.i +
z__2.i;
temp.r = z__1.r, temp.i = z__1.i;
/* L130: */
}
}
i__1 = jx;
x[i__1].r = temp.r, x[i__1].i = temp.i;
jx -= *incx;
/* L140: */
}
}
} else {
if (*incx == 1) {
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
i__2 = j;
temp.r = x[i__2].r, temp.i = x[i__2].i;
if (noconj) {
if (nounit) {
i__2 = j + j * a_dim1;
z__1.r = temp.r * a[i__2].r - temp.i * a[i__2].i,
z__1.i = temp.r * a[i__2].i + temp.i * a[
i__2].r;
temp.r = z__1.r, temp.i = z__1.i;
}
i__2 = *n;
for (i = j + 1; i <= i__2; ++i) {
i__3 = i + j * a_dim1;
i__4 = i;
z__2.r = a[i__3].r * x[i__4].r - a[i__3].i * x[
i__4].i, z__2.i = a[i__3].r * x[i__4].i +
a[i__3].i * x[i__4].r;
z__1.r = temp.r + z__2.r, z__1.i = temp.i +
z__2.i;
temp.r = z__1.r, temp.i = z__1.i;
/* L150: */
}
} else {
if (nounit) {
d_cnjg(&z__2, &a[j + j * a_dim1]);
z__1.r = temp.r * z__2.r - temp.i * z__2.i,
z__1.i = temp.r * z__2.i + temp.i *
z__2.r;
temp.r = z__1.r, temp.i = z__1.i;
}
i__2 = *n;
for (i = j + 1; i <= i__2; ++i) {
d_cnjg(&z__3, &a[i + j * a_dim1]);
i__3 = i;
z__2.r = z__3.r * x[i__3].r - z__3.i * x[i__3].i,
z__2.i = z__3.r * x[i__3].i + z__3.i * x[
i__3].r;
z__1.r = temp.r + z__2.r, z__1.i = temp.i +
z__2.i;
temp.r = z__1.r, temp.i = z__1.i;
/* L160: */
}
}
i__2 = j;
x[i__2].r = temp.r, x[i__2].i = temp.i;
/* L170: */
}
} else {
jx = kx;
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
i__2 = jx;
temp.r = x[i__2].r, temp.i = x[i__2].i;
ix = jx;
if (noconj) {
if (nounit) {
i__2 = j + j * a_dim1;
z__1.r = temp.r * a[i__2].r - temp.i * a[i__2].i,
z__1.i = temp.r * a[i__2].i + temp.i * a[
i__2].r;
temp.r = z__1.r, temp.i = z__1.i;
}
i__2 = *n;
for (i = j + 1; i <= i__2; ++i) {
ix += *incx;
i__3 = i + j * a_dim1;
i__4 = ix;
z__2.r = a[i__3].r * x[i__4].r - a[i__3].i * x[
i__4].i, z__2.i = a[i__3].r * x[i__4].i +
a[i__3].i * x[i__4].r;
z__1.r = temp.r + z__2.r, z__1.i = temp.i +
z__2.i;
temp.r = z__1.r, temp.i = z__1.i;
/* L180: */
}
} else {
if (nounit) {
d_cnjg(&z__2, &a[j + j * a_dim1]);
z__1.r = temp.r * z__2.r - temp.i * z__2.i,
z__1.i = temp.r * z__2.i + temp.i *
z__2.r;
temp.r = z__1.r, temp.i = z__1.i;
}
i__2 = *n;
for (i = j + 1; i <= i__2; ++i) {
ix += *incx;
d_cnjg(&z__3, &a[i + j * a_dim1]);
i__3 = ix;
z__2.r = z__3.r * x[i__3].r - z__3.i * x[i__3].i,
z__2.i = z__3.r * x[i__3].i + z__3.i * x[
i__3].r;
z__1.r = temp.r + z__2.r, z__1.i = temp.i +
z__2.i;
temp.r = z__1.r, temp.i = z__1.i;
/* L190: */
}
}
i__2 = jx;
x[i__2].r = temp.r, x[i__2].i = temp.i;
jx += *incx;
/* L200: */
}
}
}
}
return 0;
/* End of ZTRMV . */
} /* ztrmv_ */
/* dgemm.f -- translated by f2c (version of 23 April 1993 18:34:30).
You must link the resulting object file with the libraries:
-lf2c -lm (in that order)
*/
/* Subroutine */ int dgemm_(transa, transb, m, n, k, alpha, a, lda, b, ldb,
beta, c, ldc, transa_len, transb_len)
char *transa, *transb;
integer *m, *n, *k;
doublereal *alpha, *a;
integer *lda;
doublereal *b;
integer *ldb;
doublereal *beta, *c;
integer *ldc;
ftnlen transa_len;
ftnlen transb_len;
{
/* System generated locals */
integer a_dim1, a_offset, b_dim1, b_offset, c_dim1, c_offset, i__1, i__2,
i__3;
/* Local variables */
static integer info;
static logical nota, notb;
static doublereal temp;
static integer i, j, l, ncola;
extern logical lsame_();
static integer nrowa, nrowb;
extern /* Subroutine */ int xerbla_();
/* .. Scalar Arguments .. */
/* .. Array Arguments .. */
/* .. */
/* Purpose */
/* ======= */
/* DGEMM performs one of the matrix-matrix operations */
/* C := alpha*op( A )*op( B ) + beta*C, */
/* where op( X ) is one of */
/* op( X ) = X or op( X ) = X', */
/* alpha and beta are scalars, and A, B and C are matrices, with op( A )
*/
/* an m by k matrix, op( B ) a k by n matrix and C an m by n matrix.
*/
/* Parameters */
/* ========== */
/* TRANSA - CHARACTER*1. */
/* On entry, TRANSA specifies the form of op( A ) to be used in
*/
/* the matrix multiplication as follows: */
/* TRANSA = 'N' or 'n', op( A ) = A. */
/* TRANSA = 'T' or 't', op( A ) = A'. */
/* TRANSA = 'C' or 'c', op( A ) = A'. */
/* Unchanged on exit. */
/* TRANSB - CHARACTER*1. */
/* On entry, TRANSB specifies the form of op( B ) to be used in
*/
/* the matrix multiplication as follows: */
/* TRANSB = 'N' or 'n', op( B ) = B. */
/* TRANSB = 'T' or 't', op( B ) = B'. */
/* TRANSB = 'C' or 'c', op( B ) = B'. */
/* Unchanged on exit. */
/* M - INTEGER. */
/* On entry, M specifies the number of rows of the matrix
*/
/* op( A ) and of the matrix C. M must be at least zero.
*/
/* Unchanged on exit. */
/* N - INTEGER. */
/* On entry, N specifies the number of columns of the matrix
*/
/* op( B ) and the number of columns of the matrix C. N must be
*/
/* at least zero. */
/* Unchanged on exit. */
/* K - INTEGER. */
/* On entry, K specifies the number of columns of the matrix
*/
/* op( A ) and the number of rows of the matrix op( B ). K must
*/
/* be at least zero. */
/* Unchanged on exit. */
/* ALPHA - DOUBLE PRECISION. */
/* On entry, ALPHA specifies the scalar alpha. */
/* Unchanged on exit. */
/* A - DOUBLE PRECISION array of DIMENSION ( LDA, ka ), where ka is
*/
/* k when TRANSA = 'N' or 'n', and is m otherwise. */
/* Before entry with TRANSA = 'N' or 'n', the leading m by k
*/
/* part of the array A must contain the matrix A, otherwise
*/
/* the leading k by m part of the array A must contain the
*/
/* matrix A. */
/* Unchanged on exit. */
/* LDA - INTEGER. */
/* On entry, LDA specifies the first dimension of A as declared
*/
/* in the calling (sub) program. When TRANSA = 'N' or 'n' then
*/
/* LDA must be at least max( 1, m ), otherwise LDA must be at
*/
/* least max( 1, k ). */
/* Unchanged on exit. */
/* B - DOUBLE PRECISION array of DIMENSION ( LDB, kb ), where kb is
*/
/* n when TRANSB = 'N' or 'n', and is k otherwise. */
/* Before entry with TRANSB = 'N' or 'n', the leading k by n
*/
/* part of the array B must contain the matrix B, otherwise
*/
/* the leading n by k part of the array B must contain the
*/
/* matrix B. */
/* Unchanged on exit. */
/* LDB - INTEGER. */
/* On entry, LDB specifies the first dimension of B as declared
*/
/* in the calling (sub) program. When TRANSB = 'N' or 'n' then
*/
/* LDB must be at least max( 1, k ), otherwise LDB must be at
*/
/* least max( 1, n ). */
/* Unchanged on exit. */
/* BETA - DOUBLE PRECISION. */
/* On entry, BETA specifies the scalar beta. When BETA is
*/
/* supplied as zero then C need not be set on input. */
/* Unchanged on exit. */
/* C - DOUBLE PRECISION array of DIMENSION ( LDC, n ). */
/* Before entry, the leading m by n part of the array C must
*/
/* contain the matrix C, except when beta is zero, in which
*/
/* case C need not be set on entry. */
/* On exit, the array C is overwritten by the m by n matrix
*/
/* ( alpha*op( A )*op( B ) + beta*C ). */
/* LDC - INTEGER. */
/* On entry, LDC specifies the first dimension of C as declared
*/
/* in the calling (sub) program. LDC must be at least
*/
/* max( 1, m ). */
/* Unchanged on exit. */
/* 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. */
/* .. External Functions .. */
/* .. External Subroutines .. */
/* .. Intrinsic Functions .. */
/* .. Local Scalars .. */
/* .. Parameters .. */
/* .. */
/* .. Executable Statements .. */
/* 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 */
c_dim1 = *ldc;
c_offset = c_dim1 + 1;
c -= c_offset;
b_dim1 = *ldb;
b_offset = b_dim1 + 1;
b -= b_offset;
a_dim1 = *lda;
a_offset = a_dim1 + 1;
a -= a_offset;
/* Function Body */
nota = lsame_(transa, "N", 1L, 1L);
notb = lsame_(transb, "N", 1L, 1L);
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", 1L, 1L) && ! lsame_(transa, "T", 1L,
1L)) {
info = 1;
} else if (! notb && ! lsame_(transb, "C", 1L, 1L) && ! lsame_(transb,
"T", 1L, 1L)) {
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, 6L);
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[i + j * c_dim1] = 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[i + j * c_dim1] = *beta * c[i + j * c_dim1];
/* 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[i + j * c_dim1] = 0.;
/* L50: */
}
} else if (*beta != 1.) {
i__2 = *m;
for (i = 1; i <= i__2; ++i) {
c[i + j * c_dim1] = *beta * c[i + j * c_dim1];
/* L60: */
}
}
i__2 = *k;
for (l = 1; l <= i__2; ++l) {
if (b[l + j * b_dim1] != 0.) {
temp = *alpha * b[l + j * b_dim1];
i__3 = *m;
for (i = 1; i <= i__3; ++i) {
c[i + j * c_dim1] += temp * a[i + l * a_dim1];
/* 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[l + i * a_dim1] * b[l + j * b_dim1];
/* L100: */
}
if (*beta == 0.) {
c[i + j * c_dim1] = *alpha * temp;
} else {
c[i + j * c_dim1] = *alpha * temp + *beta * c[i + j *
c_dim1];
}
/* 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[i + j * c_dim1] = 0.;
/* L130: */
}
} else if (*beta != 1.) {
i__2 = *m;
for (i = 1; i <= i__2; ++i) {
c[i + j * c_dim1] = *beta * c[i + j * c_dim1];
/* L140: */
}
}
i__2 = *k;
for (l = 1; l <= i__2; ++l) {
if (b[j + l * b_dim1] != 0.) {
temp = *alpha * b[j + l * b_dim1];
i__3 = *m;
for (i = 1; i <= i__3; ++i) {
c[i + j * c_dim1] += temp * a[i + l * a_dim1];
/* 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[l + i * a_dim1] * b[j + l * b_dim1];
/* L180: */
}
if (*beta == 0.) {
c[i + j * c_dim1] = *alpha * temp;
} else {
c[i + j * c_dim1] = *alpha * temp + *beta * c[i + j *
c_dim1];
}
/* L190: */
}
/* L200: */
}
}
}
return 0;
/* End of DGEMM . */
} /* dgemm_ */
/* dtrmv.f -- translated by f2c (version of 23 April 1993 18:34:30).
You must link the resulting object file with the libraries:
-lf2c -lm (in that order)
*/
/* Subroutine */ int dtrmv_(uplo, trans, diag, n, a, lda, x, incx, uplo_len,
trans_len, diag_len)
char *uplo, *trans, *diag;
integer *n;
doublereal *a;
integer *lda;
doublereal *x;
integer *incx;
ftnlen uplo_len;
ftnlen trans_len;
ftnlen diag_len;
{
/* System generated locals */
integer a_dim1, a_offset, i__1, i__2;
/* Local variables */
static integer info;
static doublereal temp;
static integer i, j;
extern logical lsame_();
static integer ix, jx, kx;
extern /* Subroutine */ int xerbla_();
static logical nounit;
/* .. Scalar Arguments .. */
/* .. Array Arguments .. */
/* .. */
/* Purpose */
/* ======= */
/* DTRMV performs one of the matrix-vector operations */
/* x := A*x, or x := A'*x, */
/* where x is an n element vector and A is an n by n unit, or non-unit,
*/
/* upper or lower triangular matrix. */
/* Parameters */
/* ========== */
/* UPLO - CHARACTER*1. */
/* On entry, UPLO specifies whether the matrix is an upper or */
/* lower triangular matrix as follows: */
/* UPLO = 'U' or 'u' A is an upper triangular matrix. */
/* UPLO = 'L' or 'l' A is a lower triangular matrix. */
/* Unchanged on exit. */
/* TRANS - CHARACTER*1. */
/* On entry, TRANS specifies the operation to be performed as */
/* follows: */
/* TRANS = 'N' or 'n' x := A*x. */
/* TRANS = 'T' or 't' x := A'*x. */
/* TRANS = 'C' or 'c' x := A'*x. */
/* Unchanged on exit. */
/* DIAG - CHARACTER*1. */
/* On entry, DIAG specifies whether or not A is unit */
/* triangular as follows: */
/* DIAG = 'U' or 'u' A is assumed to be unit triangular. */
/* DIAG = 'N' or 'n' A is not assumed to be unit */
/* triangular. */
/* Unchanged on exit. */
/* N - INTEGER. */
/* On entry, N specifies the order of the matrix A. */
/* N must be at least zero. */
/* Unchanged on exit. */
/* A - DOUBLE PRECISION array of DIMENSION ( LDA, n ). */
/* Before entry with UPLO = 'U' or 'u', the leading n by n */
/* upper triangular part of the array A must contain the upper
*/
/* triangular matrix and the strictly lower triangular part of
*/
/* A is not referenced. */
/* Before entry with UPLO = 'L' or 'l', the leading n by n */
/* lower triangular part of the array A must contain the lower
*/
/* triangular matrix and the strictly upper triangular part of
*/
/* A is not referenced. */
/* Note that when DIAG = 'U' or 'u', the diagonal elements of
*/
/* A are not referenced either, but are assumed to be unity. */
/* Unchanged on exit. */
/* LDA - INTEGER. */
/* On entry, LDA specifies the first dimension of A as declared
*/
/* in the calling (sub) program. LDA must be at least */
/* max( 1, n ). */
/* Unchanged on exit. */
/* X - DOUBLE PRECISION array of dimension at least */
/* ( 1 + ( n - 1 )*abs( INCX ) ). */
/* Before entry, the incremented array X must contain the n */
/* element vector x. On exit, X is overwritten with the */
/* tranformed vector x. */
/* INCX - INTEGER. */
/* On entry, INCX specifies the increment for the elements of */
/* X. INCX must not be zero. */
/* Unchanged on exit. */
/* Level 2 Blas routine. */
/* -- 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. */
/* .. Parameters .. */
/* .. Local Scalars .. */
/* .. External Functions .. */
/* .. External Subroutines .. */
/* .. Intrinsic Functions .. */
/* .. */
/* .. Executable Statements .. */
/* Test the input parameters. */
/* Parameter adjustments */
--x;
a_dim1 = *lda;
a_offset = a_dim1 + 1;
a -= a_offset;
/* Function Body */
info = 0;
if (! lsame_(uplo, "U", 1L, 1L) && ! lsame_(uplo, "L", 1L, 1L)) {
info = 1;
} else if (! lsame_(trans, "N", 1L, 1L) && ! lsame_(trans, "T", 1L, 1L) &&
! lsame_(trans, "C", 1L, 1L)) {
info = 2;
} else if (! lsame_(diag, "U", 1L, 1L) && ! lsame_(diag, "N", 1L, 1L)) {
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, 6L);
return 0;
}
/* Quick return if possible. */
if (*n == 0) {
return 0;
}
nounit = lsame_(diag, "N", 1L, 1L);
/* 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", 1L, 1L)) {
/* Form x := A*x. */
if (lsame_(uplo, "U", 1L, 1L)) {
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[i + j * a_dim1];
/* L10: */
}
if (nounit) {
x[j] *= a[j + j * a_dim1];
}
}
/* 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[i + j * a_dim1];
ix += *incx;
/* L30: */
}
if (nounit) {
x[jx] *= a[j + j * a_dim1];
}
}
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[i + j * a_dim1];
/* L50: */
}
if (nounit) {
x[j] *= a[j + j * a_dim1];
}
}
/* 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[i + j * a_dim1];
ix -= *incx;
/* L70: */
}
if (nounit) {
x[jx] *= a[j + j * a_dim1];
}
}
jx -= *incx;
/* L80: */
}
}
}
} else {
/* Form x := A'*x. */
if (lsame_(uplo, "U", 1L, 1L)) {
if (*incx == 1) {
for (j = *n; j >= 1; --j) {
temp = x[j];
if (nounit) {
temp *= a[j + j * a_dim1];
}
for (i = j - 1; i >= 1; --i) {
temp += a[i + j * a_dim1] * 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[j + j * a_dim1];
}
for (i = j - 1; i >= 1; --i) {
ix -= *incx;
temp += a[i + j * a_dim1] * 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[j + j * a_dim1];
}
i__2 = *n;
for (i = j + 1; i <= i__2; ++i) {
temp += a[i + j * a_dim1] * 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[j + j * a_dim1];
}
i__2 = *n;
for (i = j + 1; i <= i__2; ++i) {
ix += *incx;
temp += a[i + j * a_dim1] * x[ix];
/* L150: */
}
x[jx] = temp;
jx += *incx;
/* L160: */
}
}
}
}
return 0;
/* End of DTRMV . */
} /* dtrmv_ */
/* dgemv.f -- translated by f2c (version of 23 April 1993 18:34:30).
You must link the resulting object file with the libraries:
-lf2c -lm (in that order)
*/
/* Subroutine */ int dgemv_(trans, m, n, alpha, a, lda, x, incx, beta, y,
incy, trans_len)
char *trans;
integer *m, *n;
doublereal *alpha, *a;
integer *lda;
doublereal *x;
integer *incx;
doublereal *beta, *y;
integer *incy;
ftnlen trans_len;
{
/* System generated locals */
integer a_dim1, a_offset, i__1, i__2;
/* Local variables */
static integer info;
static doublereal temp;
static integer lenx, leny, i, j;
extern logical lsame_();
static integer ix, iy, jx, jy, kx, ky;
extern /* Subroutine */ int xerbla_();
/* .. Scalar Arguments .. */
/* .. Array Arguments .. */
/* .. */
/* Purpose */
/* ======= */
/* DGEMV performs one of the matrix-vector operations */
/* y := alpha*A*x + beta*y, or y := alpha*A'*x + beta*y, */
/* where alpha and beta are scalars, x and y are vectors and A is an */
/* m by n matrix. */
/* Parameters */
/* ========== */
/* TRANS - CHARACTER*1. */
/* On entry, TRANS specifies the operation to be performed as */
/* follows: */
/* TRANS = 'N' or 'n' y := alpha*A*x + beta*y. */
/* TRANS = 'T' or 't' y := alpha*A'*x + beta*y. */
/* TRANS = 'C' or 'c' y := alpha*A'*x + beta*y. */
/* Unchanged on exit. */
/* M - INTEGER. */
/* On entry, M specifies the number of rows of the matrix A. */
/* M must be at least zero. */
/* Unchanged on exit. */
/* N - INTEGER. */
/* On entry, N specifies the number of columns of the matrix A.
*/
/* N must be at least zero. */
/* Unchanged on exit. */
/* ALPHA - DOUBLE PRECISION. */
/* On entry, ALPHA specifies the scalar alpha. */
/* Unchanged on exit. */
/* A - DOUBLE PRECISION array of DIMENSION ( LDA, n ). */
/* Before entry, the leading m by n part of the array A must */
/* contain the matrix of coefficients. */
/* Unchanged on exit. */
/* LDA - INTEGER. */
/* On entry, LDA specifies the first dimension of A as declared
*/
/* in the calling (sub) program. LDA must be at least */
/* max( 1, m ). */
/* Unchanged on exit. */
/* X - DOUBLE PRECISION array of DIMENSION at least */
/* ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n' */
/* and at least */
/* ( 1 + ( m - 1 )*abs( INCX ) ) otherwise. */
/* Before entry, the incremented array X must contain the */
/* vector x. */
/* Unchanged on exit. */
/* INCX - INTEGER. */
/* On entry, INCX specifies the increment for the elements of */
/* X. INCX must not be zero. */
/* Unchanged on exit. */
/* BETA - DOUBLE PRECISION. */
/* On entry, BETA specifies the scalar beta. When BETA is */
/* supplied as zero then Y need not be set on input. */
/* Unchanged on exit. */
/* Y - DOUBLE PRECISION array of DIMENSION at least */
/* ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n' */
/* and at least */
/* ( 1 + ( n - 1 )*abs( INCY ) ) otherwise. */
/* Before entry with BETA non-zero, the incremented array Y */
/* must contain the vector y. On exit, Y is overwritten by the
*/
/* updated vector y. */
/* INCY - INTEGER. */
/* On entry, INCY specifies the increment for the elements of */
/* Y. INCY must not be zero. */
/* Unchanged on exit. */
/* Level 2 Blas routine. */
/* -- 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. */
/* .. Parameters .. */
/* .. Local Scalars .. */
/* .. External Functions .. */
/* .. External Subroutines .. */
/* .. Intrinsic Functions .. */
/* .. */
/* .. Executable Statements .. */
/* Test the input parameters. */
/* Parameter adjustments */
--y;
--x;
a_dim1 = *lda;
a_offset = a_dim1 + 1;
a -= a_offset;
/* Function Body */
info = 0;
if (! lsame_(trans, "N", 1L, 1L) && ! lsame_(trans, "T", 1L, 1L) && !
lsame_(trans, "C", 1L, 1L)) {
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, 6L);
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", 1L, 1L)) {
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", 1L, 1L)) {
/* 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[i + j * a_dim1];
/* 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[i + j * a_dim1];
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[i + j * a_dim1] * 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[i + j * a_dim1] * x[ix];
ix += *incx;
/* L110: */
}
y[jy] += *alpha * temp;
jy += *incy;
/* L120: */
}
}
}
return 0;
/* End of DGEMV . */
} /* dgemv_ */
/* zgerc.f -- translated by f2c (version of 23 April 1993 18:34:30).
You must link the resulting object file with the libraries:
-lf2c -lm (in that order)
*/
/* Subroutine */ int zgerc_(m, n, alpha, x, incx, y, incy, a, lda)
integer *m, *n;
doublecomplex *alpha, *x;
integer *incx;
doublecomplex *y;
integer *incy;
doublecomplex *a;
integer *lda;
{
/* System generated locals */
integer a_dim1, a_offset, i__1, i__2, i__3, i__4, i__5;
doublecomplex z__1, z__2;
/* Builtin functions */
void d_cnjg();
/* Local variables */
static integer info;
static doublecomplex temp;
static integer i, j, ix, jy, kx;
extern /* Subroutine */ int xerbla_();
/* .. Scalar Arguments .. */
/* .. Array Arguments .. */
/* .. */
/* Purpose */
/* ======= */
/* ZGERC performs the rank 1 operation */
/* A := alpha*x*conjg( y' ) + A, */
/* where alpha is a scalar, x is an m element vector, y is an n element
*/
/* vector and A is an m by n matrix. */
/* Parameters */
/* ========== */
/* M - INTEGER. */
/* On entry, M specifies the number of rows of the matrix A. */
/* M must be at least zero. */
/* Unchanged on exit. */
/* N - INTEGER. */
/* On entry, N specifies the number of columns of the matrix A.
*/
/* N must be at least zero. */
/* Unchanged on exit. */
/* ALPHA - COMPLEX*16 . */
/* On entry, ALPHA specifies the scalar alpha. */
/* Unchanged on exit. */
/* X - COMPLEX*16 array of dimension at least */
/* ( 1 + ( m - 1 )*abs( INCX ) ). */
/* Before entry, the incremented array X must contain the m */
/* element vector x. */
/* Unchanged on exit. */
/* INCX - INTEGER. */
/* On entry, INCX specifies the increment for the elements of */
/* X. INCX must not be zero. */
/* Unchanged on exit. */
/* Y - COMPLEX*16 array of dimension at least */
/* ( 1 + ( n - 1 )*abs( INCY ) ). */
/* Before entry, the incremented array Y must contain the n */
/* element vector y. */
/* Unchanged on exit. */
/* INCY - INTEGER. */
/* On entry, INCY specifies the increment for the elements of */
/* Y. INCY must not be zero. */
/* Unchanged on exit. */
/* A - COMPLEX*16 array of DIMENSION ( LDA, n ). */
/* Before entry, the leading m by n part of the array A must */
/* contain the matrix of coefficients. On exit, A is */
/* overwritten by the updated matrix. */
/* LDA - INTEGER. */
/* On entry, LDA specifies the first dimension of A as declared
*/
/* in the calling (sub) program. LDA must be at least */
/* max( 1, m ). */
/* Unchanged on exit. */
/* Level 2 Blas routine. */
/* -- 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. */
/* .. Parameters .. */
/* .. Local Scalars .. */
/* .. External Subroutines .. */
/* .. Intrinsic Functions .. */
/* .. */
/* .. Executable Statements .. */
/* Test the input parameters. */
/* Parameter adjustments */
a_dim1 = *lda;
a_offset = a_dim1 + 1;
a -= a_offset;
--y;
--x;
/* 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_("ZGERC ", &info, 6L);
return 0;
}
/* Quick return if possible. */
if (*m == 0 || *n == 0 || alpha->r == 0. && alpha->i == 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) {
i__2 = jy;
if (y[i__2].r != 0. || y[i__2].i != 0.) {
d_cnjg(&z__2, &y[jy]);
z__1.r = alpha->r * z__2.r - alpha->i * z__2.i, z__1.i =
alpha->r * z__2.i + alpha->i * z__2.r;
temp.r = z__1.r, temp.i = z__1.i;
i__2 = *m;
for (i = 1; i <= i__2; ++i) {
i__3 = i + j * a_dim1;
i__4 = i + j * a_dim1;
i__5 = i;
z__2.r = x[i__5].r * temp.r - x[i__5].i * temp.i, z__2.i =
x[i__5].r * temp.i + x[i__5].i * temp.r;
z__1.r = a[i__4].r + z__2.r, z__1.i = a[i__4].i + z__2.i;
a[i__3].r = z__1.r, a[i__3].i = z__1.i;
/* 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) {
i__2 = jy;
if (y[i__2].r != 0. || y[i__2].i != 0.) {
d_cnjg(&z__2, &y[jy]);
z__1.r = alpha->r * z__2.r - alpha->i * z__2.i, z__1.i =
alpha->r * z__2.i + alpha->i * z__2.r;
temp.r = z__1.r, temp.i = z__1.i;
ix = kx;
i__2 = *m;
for (i = 1; i <= i__2; ++i) {
i__3 = i + j * a_dim1;
i__4 = i + j * a_dim1;
i__5 = ix;
z__2.r = x[i__5].r * temp.r - x[i__5].i * temp.i, z__2.i =
x[i__5].r * temp.i + x[i__5].i * temp.r;
z__1.r = a[i__4].r + z__2.r, z__1.i = a[i__4].i + z__2.i;
a[i__3].r = z__1.r, a[i__3].i = z__1.i;
ix += *incx;
/* L30: */
}
}
jy += *incy;
/* L40: */
}
}
return 0;
/* End of ZGERC . */
} /* zgerc_ */
/* izamax.f -- translated by f2c (version of 23 April 1993 18:34:30).
You must link the resulting object file with the libraries:
-lf2c -lm (in that order)
*/
integer izamax_(n, zx, incx)
integer *n;
doublecomplex *zx;
integer *incx;
{
/* System generated locals */
integer ret_val, i__1;
/* Local variables */
static doublereal smax;
static integer i;
extern doublereal dcabs1_();
static integer ix;
/* finds the index of element having max. absolute value. */
/* jack dongarra, 1/15/85. */
/* modified to correct problem with negative increment, 8/21/90. */
/* Parameter adjustments */
--zx;
/* Function Body */
ret_val = 0;
if (*n < 1) {
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;
if (*incx < 0) {
ix = (-(*n) + 1) * *incx + 1;
}
smax = dcabs1_(&zx[ix]);
ix += *incx;
i__1 = *n;
for (i = 2; i <= i__1; ++i) {
if (dcabs1_(&zx[ix]) <= smax) {
goto L5;
}
ret_val = i;
smax = dcabs1_(&zx[ix]);
L5:
ix += *incx;
/* L10: */
}
return ret_val;
/* code for increment equal to 1 */
L20:
smax = dcabs1_(&zx[1]);
i__1 = *n;
for (i = 2; i <= i__1; ++i) {
if (dcabs1_(&zx[i]) <= smax) {
goto L30;
}
ret_val = i;
smax = dcabs1_(&zx[i]);
L30:
;
}
return ret_val;
} /* izamax_ */
/* lsame.f -- translated by f2c (version of 23 April 1993 18:34:30).
You must link the resulting object file with the libraries:
-lf2c -lm (in that order)
*/
logical lsame_(ca, cb, ca_len, cb_len)
char *ca, *cb;
ftnlen ca_len;
ftnlen cb_len;
{
/* System generated locals */
logical ret_val;
/* Local variables */
static integer inta, intb, zcode;
/* -- LAPACK auxiliary routine (version 1.0) -- */
/* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
/* Courant Institute, Argonne National Lab, and Rice University */
/* February 29, 1992 */
/* .. Scalar Arguments .. */
/* .. */
/* Purpose */
/* ======= */
/* LSAME returns .TRUE. if CA is the same letter as CB regardless of */
/* case. */
/* Arguments */
/* ========= */
/* CA (input) CHARACTER*1 */
/* CB (input) CHARACTER*1 */
/* CA and CB specify the single characters to be compared. */
/* .. Intrinsic Functions .. */
/* .. */
/* .. Local Scalars .. */
/* .. */
/* .. Executable Statements .. */
/* Test if the characters are equal */
ret_val = *ca == *cb;
if (ret_val) {
return ret_val;
}
/* Now test for equivalence if both characters are alphabetic. */
zcode = 'Z';
/* Use 'Z' rather than 'A' so that ASCII can be detected on Prime */
/* machines, on which ICHAR returns a value with bit 8 set. */
/* ICHAR('A') on Prime machines returns 193 which is the same as */
/* ICHAR('A') on an EBCDIC machine. */
inta = *ca;
intb = *cb;
if (zcode == 90 || zcode == 122) {
/* ASCII is assumed - ZCODE is the ASCII code of either lower o
r */
/* upper case 'Z'. */
if (inta >= 97 && inta <= 122) {
inta += -32;
}
if (intb >= 97 && intb <= 122) {
intb += -32;
}
} else if (zcode == 233 || zcode == 169) {
/* EBCDIC is assumed - ZCODE is the EBCDIC code of either lower
or */
/* upper case 'Z'. */
if (inta >= 129 && inta <= 137 || inta >= 145 && inta <= 153 || inta
>= 162 && inta <= 169) {
inta += 64;
}
if (intb >= 129 && intb <= 137 || intb >= 145 && intb <= 153 || intb
>= 162 && intb <= 169) {
intb += 64;
}
} else if (zcode == 218 || zcode == 250) {
/* ASCII is assumed, on Prime machines - ZCODE is the ASCII cod
e */
/* plus 128 of either lower or upper case 'Z'. */
if (inta >= 225 && inta <= 250) {
inta += -32;
}
if (intb >= 225 && intb <= 250) {
intb += -32;
}
}
ret_val = inta == intb;
/* RETURN */
/* End of LSAME */
return ret_val;
} /* lsame_ */
/* zaxpy.f -- translated by f2c (version of 23 April 1993 18:34:30).
You must link the resulting object file with the libraries:
-lf2c -lm (in that order)
*/
/* Subroutine */ int zaxpy_(n, za, zx, incx, zy, incy)
integer *n;
doublecomplex *za, *zx;
integer *incx;
doublecomplex *zy;
integer *incy;
{
/* System generated locals */
integer i__1, i__2, i__3, i__4;
doublecomplex z__1, z__2;
/* Local variables */
static integer i;
extern doublereal dcabs1_();
static integer ix, iy;
/* constant times a vector plus a vector. */
/* jack dongarra, 3/11/78. */
/* Parameter adjustments */
--zy;
--zx;
/* Function Body */
if (*n <= 0) {
return 0;
}
if (dcabs1_(za) == 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) {
i__2 = iy;
i__3 = iy;
i__4 = ix;
z__2.r = za->r * zx[i__4].r - za->i * zx[i__4].i, z__2.i = za->r * zx[
i__4].i + za->i * zx[i__4].r;
z__1.r = zy[i__3].r + z__2.r, z__1.i = zy[i__3].i + z__2.i;
zy[i__2].r = z__1.r, zy[i__2].i = z__1.i;
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) {
i__2 = i;
i__3 = i;
i__4 = i;
z__2.r = za->r * zx[i__4].r - za->i * zx[i__4].i, z__2.i = za->r * zx[
i__4].i + za->i * zx[i__4].r;
z__1.r = zy[i__3].r + z__2.r, z__1.i = zy[i__3].i + z__2.i;
zy[i__2].r = z__1.r, zy[i__2].i = z__1.i;
/* L30: */
}
return 0;
} /* zaxpy_ */
/* zdscal.f -- translated by f2c (version of 23 April 1993 18:34:30).
You must link the resulting object file with the libraries:
-lf2c -lm (in that order)
*/
/* Subroutine */ int zdscal_(n, da, zx, incx)
integer *n;
doublereal *da;
doublecomplex *zx;
integer *incx;
{
/* System generated locals */
integer i__1, i__2, i__3;
doublecomplex z__1, z__2;
/* Local variables */
static integer i, ix;
/* scales a vector by a constant. */
/* jack dongarra, 3/11/78. */
/* modified to correct problem with negative increment, 8/21/90. */
/* Parameter adjustments */
--zx;
/* Function Body */
if (*n <= 0) {
return 0;
}
if (*incx == 1) {
goto L20;
}
/* code for increment not equal to 1 */
ix = 1;
if (*incx < 0) {
ix = (-(*n) + 1) * *incx + 1;
}
i__1 = *n;
for (i = 1; i <= i__1; ++i) {
i__2 = ix;
z__2.r = *da, z__2.i = 0.;
i__3 = ix;
z__1.r = z__2.r * zx[i__3].r - z__2.i * zx[i__3].i, z__1.i = z__2.r *
zx[i__3].i + z__2.i * zx[i__3].r;
zx[i__2].r = z__1.r, zx[i__2].i = z__1.i;
ix += *incx;
/* L10: */
}
return 0;
/* code for increment equal to 1 */
L20:
i__1 = *n;
for (i = 1; i <= i__1; ++i) {
i__2 = i;
z__2.r = *da, z__2.i = 0.;
i__3 = i;
z__1.r = z__2.r * zx[i__3].r - z__2.i * zx[i__3].i, z__1.i = z__2.r *
zx[i__3].i + z__2.i * zx[i__3].r;
zx[i__2].r = z__1.r, zx[i__2].i = z__1.i;
/* L30: */
}
return 0;
} /* zdscal_ */
/* zdotu.f -- translated by f2c (version of 23 April 1993 18:34:30).
You must link the resulting object file with the libraries:
-lf2c -lm (in that order)
*/
/* Double Complex */ int zdotu_( ret_val, n, zx, incx, zy, incy)
doublecomplex * ret_val;
integer *n;
doublecomplex *zx;
integer *incx;
doublecomplex *zy;
integer *incy;
{
/* System generated locals */
integer i__1, i__2, i__3;
doublecomplex z__1, z__2;
/* Local variables */
static integer i;
static doublecomplex ztemp;
static integer ix, iy;
/* forms the dot product of two vectors. */
/* jack dongarra, 3/11/78. */
/* Parameter adjustments */
--zy;
--zx;
/* Function Body */
ztemp.r = 0., ztemp.i = 0.;
ret_val->r = 0., ret_val->i = 0.;
if (*n <= 0) {
return ;
}
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) {
i__2 = ix;
i__3 = iy;
z__2.r = zx[i__2].r * zy[i__3].r - zx[i__2].i * zy[i__3].i, z__2.i =
zx[i__2].r * zy[i__3].i + zx[i__2].i * zy[i__3].r;
z__1.r = ztemp.r + z__2.r, z__1.i = ztemp.i + z__2.i;
ztemp.r = z__1.r, ztemp.i = z__1.i;
ix += *incx;
iy += *incy;
/* L10: */
}
ret_val->r = ztemp.r, ret_val->i = ztemp.i;
return ;
/* code for both increments equal to 1 */
L20:
i__1 = *n;
for (i = 1; i <= i__1; ++i) {
i__2 = i;
i__3 = i;
z__2.r = zx[i__2].r * zy[i__3].r - zx[i__2].i * zy[i__3].i, z__2.i =
zx[i__2].r * zy[i__3].i + zx[i__2].i * zy[i__3].r;
z__1.r = ztemp.r + z__2.r, z__1.i = ztemp.i + z__2.i;
ztemp.r = z__1.r, ztemp.i = z__1.i;
/* L30: */
}
ret_val->r = ztemp.r, ret_val->i = ztemp.i;
return ;
} /* zdotu_ */
syntax highlighted by Code2HTML, v. 0.9.1