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